Development and Application of Tools to Characterize Transiting Astrophysical Systems

Since the discovery of the first exoplanets (planets outside our Solar System) more than 20 years ago, there has been an increasing need for photometric and spectroscopic models to characterize these systems. While imaging has been used extensively for Solar System bodies and extended objects like galaxies, the small angular extent of typical planetary systems makes it difficult or impossible to resolve them. Spatially integrated observations like measuring the total brightness or spectrum, however, can be conducted at a resonable cost. This thesis focuses on photometric models in the context of transiting systems, which exhibit a number of phenomena that can be exploited for characterization.

First, we showcase the popular methods of transiting exoplanet discovery and characterization by ground based observations on the hot Jupiter HAT-P-27b. We demonstrate how transits allow us to constrain planetary mass, radius, and orbital inclination, which would not be possible based only on, for example, radial velocity measurements.

Next, we perform reflection spectroscopy on HAT-P-1b, another hot Jupiter, using the binary companion of the host star as a reference to remove systematic errors from the signal. Here the transiting nature of the system allows us to look for the very faint light reflected by the planet.

We also apply the idea of planetary transits to investigate the feasibility of transit observations in astrophysical systems of very different scale: stars in galactic nuclei potentially transiting the accretion disk of the supermassive black hole in the galactic center.

Finally, we focus on mapping spots on the stellar surface using transits. This method has been used for a decade, and helped constrain stellar rotation or orbital geometry in a number systems. We study starspots on HAT-P-11 not only to learn more about stellar rotation, but also to investigate the size and contrast of the spots themselves.

See DAC.pdf

Page intentionally left blank.

Development and Application of Tools

to Characterize Transiting Astrophysical Systems

[1.5em] A dissertation presented

[1.5em] by

[1.5em] Bence Béky

[1.5em] to

[1.5em] The Department of Astronomy

in partial fulfillment of the requirements

for the degree of

Doctor of Philosophy

in the subject of

Astronomy and Astrophysics

[1.5em] Harvard University

Cambridge, Massachusetts

[1.5em] 2014 April

Copyright © 2014 Bence Béky

All rights reserved.

List of Figures


I shall start by giving due credit to people who helped me on the journey before I started my career as an astronomer. I thank my parents, my sisters, and other members of my extended but close family; fellow students and educators during my formal and informal education, especially high school; and friends. These people play a key role in my life, and are responsible for that when I enrolled to Harvard four years ago, I felt I was ready.

I would like to express gratitude to a large number of astronomers who accepted my younger self with no formal training in their discipline as one of them, and helped me along my apprenticeship. First of all, I thank Gáspár Bakos, who recruited me as a research assistant and introduced me to exoplanets and people in this field. I thank Matt Holman for taking over the role of my advisor after Gáspár left for another university and was not successful in convincing me to leave the CfA behind.

I am also grateful to a large number of people who volunteered their time and advice without being labeled advisors: Dave Kipping and Bob Noyes, with whom I had regular weekly meetings throughout the last two years, Bence Kocsis, who gave me the opportunity to lead the research paper based on his idea, and many, many other people, whom I had the opportunity to consult with regularily or occasionally. This includes collaborators and coauthors, mentors, lecturers, and people I met at SSP coffee, exoplanet pizza lunch, and other scientific gatherings. I am especially indebted to the Mount Hopkins Observatory staff for their infinite patience with me as I made my debut as an observer.

When one thinks about whom they learned the most from, it is easy to think only about senior people and forget about peers. However, my confusion was often so great that I would have been ashamed to seek help of anyone other than the folks who still remember vividly their own struggles with astronomy: fellow graduate students. They did not leave a single question of mine unanswered, and every interaction left me humbled by how much more there was for me to learn.

I would like to acknowledge the hard work of Peg Herlihy, Robb Scholten, and Donna Adams, without whom I would not have had an office, an advisor, stipend, healthcare, I could not have taken a single class, I could not have taught, I could not have loaned a single book from the library, I could not have entered the building after 17:30, and I would have been deported from the US, just to name a few things.

Being an astronomer is about more than just astronomy, and being a graduate student is about more than just being in graduate school. I am grateful for all the people who kept reminding me of this in the past years, including friends that I met through CouchSurfing, through the local Hungarian community, and through other stochastic processes; fellow graduate students who showed me that they also have a non-astronomer side; but most of all, Katherine Beaty, who has earned a special place in my heart and on my ring finger. And I thank Owen Gingerich and Charles Alcock for making it possible that she put that ring on my finger in the dome of the Great Refractor.

So long, and thanks for all the fish

Douglas Adams (1984)

To my family

Chapter 1 Introduction

1.1 History of Solar System planets

History always has a few tricks up its frayed sleeve. It’s been around a long time.

Terry Pratchett (1987)

The first planets known to humanity are those visible to the naked eye: Mercury, Venus, Mars, Jupiter, and Saturn. The earliest written observations of these objects come from the Babylonian civilization, dating back to second millenium BC (Sachs, 1974). The word “planet” originates from ancient Greek πλανήτης (planētēs), meaning wanderer, as they move across the sky with respect to stars. Mercury, Venus, Mars, Jupiter, and Saturn were the five planets in Greek astronomy after Pythagoras or Parmenides correctly identified the evening star and morning star as the same object, Venus (Aphrodite in Greek) (Burnet, 2007). Paradoxially, Earth, the planet closest to humankind, with the largest apparent size, has only been proposed to be a planet in the early 1500s, by the Polish astronomer Mikołaj Kopernik (Latinized as Nicolaus Copernicus). His framework, known as the heliocentric system, only became widely accepted two centuries later, with the advent of telescopic observations.

The rest of the Solar System planets were discovered using what were considered large telescopes at the time: Uranus in 1781 by William Herschel, and Neptune by Urbain Le Verrier and Johann Gottfried Galle in 1846.

Pluto was discovered in 1930 by Clyde Tombaugh, although it had been captured on photographic plates fifteen times before since as early as 1909. It was reclassified as a dwarf planet in 2006, leaving us with eight Solar System planets. However, it has not been the first object demoted from being a planet: the Sun and the Moon were thought of as planets until the acceptance of the heliocentric system, although the concept of planet had necessarily been less refined at the time. After the discovery of Uranus, a number of other celestial bodies were discovered and labelled as planets, most notably Ceres, Pallas, Juno, and Vesta, between 1801 and 1807. At the time of the discovery of Neptune, it was considered the 13th planet, until around 1850, when the others were reclassified as asteroids. Ceres was reclassified again in 2006, as a dwarf planet, together with Pluto (Forbes, 1971; Hughes & Marsden, 2007).

The quest for discovering exoplanets (or extrasolar planets, that is, planets outside our Solar System) was held back until very recently by their distance: Centauri Bb, the closest known potential exoplanet, is ten thousand times as far from the Earth as Neptune, the furthest planet within the Solar System. It is not possible to look up the night sky and spot such a remote planet with the naked eye, and direct observations even with very powerful telescopes are rare.

However, not all Solar System planets were discovered by scouring around the sky: the presence of Neptune was deducted by the irregular motion of Uranus. In 1843, John Couch Adams started to examine observations of Uranus to predict the position of Neptune, which was discovered a few years later based on the independent calculations of Urbain Le Verrier. This teaches us an important lesson: looking into your telescope pointed at a random field on the sky is not the only way to look for planets. Sometimes, observing another object like a planet or a star will provide indirect evidence of a new planet, even without being able to directly see it.

Another important lesson from history is that sometimes what was thought to be a planet turns out to be something else. In case of the Solar System, only the meaning of the word “planet” evolved, without actually changing the physical nature of the Sun, the Moon, Ceres, Pallas, Juno, Vesta, Pluto, and many others. For indirectly discovered exoplanets, however, reanalysis of data or gathering more observations might reveal that what was thought to be a planet is too massive to be one, does not exist at all, or is indeed a planet but with a different orbital period.

1.2 Exoplanet detection methods

L’essentiel est invisible pour les yeux.

Antoine de Saint Exupéry (1943)

Exoplanets can be detected with a number of methods. We briefly review the most important ones in this section. The takeaway message is that no method is superior: they have different biases, limitations, and costs.

Biases are important, because we need to account for them when inferring the ensemble distribution of planet parameters from the distribution of known planets. They are also important when choosing a method, designing an instrument, and setting up an observing strategy for the specific goal of finding planets of certain properties, for example, ones that might potentially harbor life.

1.2.1 Direct imaging

Just like Uranus was discovered as a small bright object in the field of the telescope, it is possible to resolve an exoplanet, that is, separate the light from two objects: the planet and its host star. For a telescope with an optical element of diameter , at wavelength , diffraction of light limits the angular resolution to

A quick order-of-magnitude calculation tells us that for a nearby star at , observed by a 10 meter telescope at , this translates to a projected star–planet distance of . At this distance, the starlight reflected by the planet is too faint compared to the star to be detected with current technology, therefore direct imaging is limited to the thermal radiation originating from the planets themselves, typically observed in the infrared. This is essentially different from the situation in our Solar System, where every planet on the sky was first detected using light that originated from the Sun and reflected off the planet. Direct imaging observations favor young, massive planets, because they have more residual heat, and a larger surface to radiate. These biases are reflected by the fact that until recently, all directly imaged exoplanets were younger than (Kuzuhara et al., 2013), and with the exception of Fomalhaut b, all have a inferred planetary mass at least four times that of Jupiter. So far, eight planets have been discovered by direct imaging.\@footnotemark

1.2.2 Timing variations

While direct imaging is the natural extension of observations of Solar System planets, the first exoplanets were discovered using different methods. One of the first discoveries is based on timing variation of radio signals. The millisecond pulsar PSR 1527 + 12 is a rapidly rotating object that emits radio waves synchronously to its rotation at a highly regular rate, every . As a planet orbits the pulsar, the pulsar does not remain stationary: in fact, both of them orbit their common barycenter. During this orbit, both objects periodically approach and recede as seen by the observer. Since the radio emission rate of the pulsar is very stable, this radial displacement can be detected in the timing of the radio signals, because it takes light a little longer to reach the observer when the pulsar is furher away. Two planets were detected around PSR 1527 + 12 via this effect, even though they change the pulsation period by only (Wolszczan & Frail, 1992).

Timing effects are also used to discover non-transiting planets via the transit timing variations (TTV) of transiting planets in the same system. For example, the transit times of Kepler-19b deviate from strictly periodic by up to approximately five minutes, revealing the presence of Kepler-19c (Ballard et al., 2011). Dynamical arguments give a lower limit of 0.1 for the mass of Kepler-19c. Another example of analyzing TTVs to constrain planetary masses is the Kepler-88 system (Nesvorný et al., 2013). This method can sometimes also be applied to targets that are too faint for radial velocity measurements.

Timing variations are limited to systems exhibiting very regularly periodic behaviour. In the above cases, we are limited to planets transiting a millisecond pulsar, and planets that strongly interact with already discovered transiting ones. To date, less than two dozen planets have been discovered via some kind of timing phenomenon.\@footnotemark

1.2.3 Radial velocity

Just like with millisecond pulsar timing, planets orbiting stars can be detected by their gravitational pull on the star, if we can measure the radial velocity of the host star precisely enough. This can be done via spectroscopic observations: the observed wavelength of absorption lines in the stellar spectrum are influenced by the radial velocity of the star via the Doppler effect (Marcy & Butler, 1992; Seager, 2011). Targeted surveys of main sequence stars are fruitful methods of discovering new planets, albeit very expensive in terms of observations, because typically dozens of observations with medium to large telescopes are required for discovering such planets. So far, hundreds of planets have been discovered via radial velocity observations.\@footnotemark

Radial velocity measurements are important not only for discovering planets, but also for the mass determination of planets discovered by other methods, like transits. However, when other methods are not available for characterizing the planet, radial velocity observations by themselves are not able to resolve the degeneracy between planetary mass and orbital inclination , as they only constrain the product .

The biases associated with radial velocity measurements are the following: first, since the primary measured quantity is radial velocity, massive planets and those with short periods are favored as they result in larger stellar radial velocity amplitude. Second, as rapidly rotating stars have broad spectral features that hinder precise radial velocity determination, rapidly rotating stars like massive main sequence A stars and very young stars are usually excluded from observing campaigns. Third, spectroscopic reference lines from ionine cells and ThAr lamps are densest at optical wavelengths, therefore RV surveys were historically biased against cooler stars like M dwarfs, because their emission mostly falls in the infrared.

1.2.4 Transits

The orbital orientation of a planet might be such that it gets in between the star and the observer once every orbit. In this case, the planet blocks part of the star’s light, which we observe as a periodical dip in the stellar lightcurve (brightness as a function of time). This method is fairly unique in that the depth of the transit reveals the planetary radius relative to the stellar radius. In addition, the transit lightcurve by itself allows for measuring the density of the host star (Seager & Mallén-Ornelas, 2003; Seager, 2011).

Compared to other methods, transits are relatively inexpensive to observe: a transiting survey can monitor tens of thousands of stars in a single field simultaneously for brightness variations. It is therefore not suprising that most known transiting exoplanets have been discovered by transits, more than a thousand by now.\@footnotemark

However, special configurations of background stars or hierachical triple systems might mimic the behaviour a transiting system, thus usually other observations are required to confirm the planetary nature of the system. Therefore it is not uncommon for transiting exoplanet surveys to focus their efforts on targets that are known to be feasible targets for radial velocity measurements, that is, relatively bright stars that are not fast rotators. This is also beneficial because transits tell us about the planetary size, and radial velocity measurements confine the mass, providing us with the unique opportunity to determine planetary density and thus potentially learn about the composition of the planet.

Ground-based transiting exoplanet surveys will be discussed in detail in Section 1.3.

1.2.5 Gravitational microlensing

Gravitational microlensing happens when two stars are almost perfectly aligned as seen by the observer. In this case, the gravitational field of the foreground star acts a magnifying lens, enhancing the light from the background star. Such events usually last for days of weeks, because of the motion of the background star, the foreground star, and the Solar System with the observer in it, in space. If the lensing star has a planet, it may be revealed as a second spike in the microlensing light curve if the planet is in a favorable location during the event. To date, 18 planets have been discovered\@footnotemark via this phenomenon: 10 by the MOA project (Street et al., 2013), and 8 by OGLE (Kains et al., 2013).

The unique charactestic of microlensing is that these are rare events, usually only happening once to a system. Therefore this method cannot constrain orbital period directly, only planetary mass and projected separation.

1.2.6 Orbital brightness modulation

A planet might influence the total brightness of the system even if it does not transit the star. For example, as the planet goes through phases during the orbit, the amount of reflected starlight changes periodically. This has been measured for a number of known planets, but the first planets discovered through this effect are Kepler-70b and Kepler-70c (Charpinet et al., 2011). This method is biased toward large, close-in planets with high albedos, because they reflect more starlight.

In addition, radial motion of stars result in tiny brightness variations due to relativistic beaming, allowing for detection of planets. The first such discovery was Kepler-76b by Faigler et al. (2013). In the lightcurve analysis, the authors also account for starlight reflected by the planet, and variations due to the ellipsoidal shape of the star caused by tidal forces.

1.2.7 Polarimetry

Starlight is unpolarized, whereas reflection from the planetary atmosphere is partially polarized. This effect has been observed, for example, for the exoplanet HD 189733b (Berdyugina et al., 2008), but no previously unknown planets have been discovered yet with this method.

1.2.8 Astrometry

Just like a planet induces periodic motion of its host star in the radial direction (unless the orbit is in the sky plane), it does so in the sky plane too. Therefore astrometry, that is, careful observation of the position of the star, could potentially reveal exoplanets. This effect has not been observed for any planetary host yet.

However, astrometry sometimes plays an important role in ruling out transiting binary systems mimicking transit behaviour. For the transit depth to be reasonable for a planet, the light of a binary system must be blended together with light from a field star. This system might not be resolved by the instrument, but astrometric measurements could reveal that the light is coming more from the field star during transit than at other times. Note, however, that it is not physical motion of a star that is detected in this case, but an artifact due to the structure of the blended source. Not detecting such “astrometric” variations decreases the likelihood of a false positive (Torres et al., 2011).

1.3 Planet factories

The most successful transiting exoplanet discovery project by the number of planets is the Kepler Space Satellite, with hundreds of confirmed and validated detections to date (Borucki et al., 2010, 2011; Batalha et al., 2013; Burke et al., 2014; Rowe et al., 2014). While there are other satellites looking for exoplanets, the next most fruitful efforts are ground-based surveys: WASP (Pollacco et al., 2006) with 64 planets, and HATNet (Bakos et al., 2004) and HATSouth (Bakos et al., 2013) with a total of 43 planets discovered to date.111, retrieved 2014-03-25 In this section, we highlight a few aspects of these ground-based surveys.

1.3.1 Instrument design

The HATNet, HATSouth, and WASP surveys use commercially available telephoto lenses or telescopes measuring 11 cm, 18 cm, and 20 cm in diameter, respectively. While part of the reason for the relatively small size compared to other professional telescopes is to keep the design simple and costs low, the planetary discovery rate does not seem to suffer from the relatively small apertures. The reason is that transiting system candidates are typically confirmed using radial velocity (RV) follow-up measurements, which would become costly for stars fainter than visual magnitude, and for brighter stars, these apertures already collect enough light for successful identification. Brighter stars require shorter exposures not only for RV measurements, but also for spectroscopic observations of transits, which would allow us to learn about the composition of the planetary atmosphere. With such observations, longer exposures can be overly complicated: with the Hubble Space Telescope, for example, each exposure is limited to approximately 40 minutes, before the spacecraft hides behind the Earth on its orbit. Also, if the required exposure time exceeds the length of the transit, observations taken at multiple transits need to be combined.

Given that the scientific goal of HATNet, HATSouth, and WASP is to find a large number of transiting exoplanets that are suitable for various follow-up observations, and that the entire sky has not been surveyed yet for transiting exoplanets, they were designed to be wide angle, shallow surveys. They employ fast optics, with focal ratios between and , resulting in a field between 4 and 8 on the side, with tens of thousands of stars down to fourteenth magnitude for each telescope.

In a ground-based survey, the observed brightness of targets is influenced by systematic effects due to the atmosphere. The amount of light taken away by the atmosphere (extinction) depends on the amount of Earth’s atmosphere the light from the star has to traverse to reach the telescope (airmass). Extinction is function of humidity, but also zenith angle, so it changes throughout the night as the target star rises and sets, and as weather conditions change. To separate this effect from the astrophysical brightness variations, a number of stable stars are selected in the same field, and their brightness is compared. To complicate things, the wavelength-dependent scattering in air results in different extinction for stars of different color, a phenomenon known as differential extinction. This necessitates a large number of stars in the field among which appropriate ones can be chosen as photometric reference, also justifying wide fields for photometric surveys. The HATNet and HATSouth surveys use the Trend Filtering Algorithm (Kovács et al., 2005) to remove systematics using lightcurves of reference stars, and the External Parameter Decorrelation (Bakos et al., 2007b) to account for effects based on where the telescope was pointing and where the star fell on the detector.

1.3.2 Observing strategy

After deciding on the hardware: shallow, wide field versus deep, narrow field, the next major decision is in the observing strategy: observe one field for many months, or a new field every week. At the two extremes, a very short campaign could not cover two consecutive transits that are indispensible for determining the orbital period, and a very long campaign would have a small marginal discovery rate, because all short period planets are already discovered, while long period transiting planets are a lot less frequent due to the lower geometric probability of a transit.

It is important to note that typical transits last for a few hours, therefore a single-longitude survey has a good chance of missing any particular one. In light of this, it is not surprising that the continuous sky coverage (weather permitted) provided by the three HATSouth sites spread out in longitude results in an expected planet recovery fraction exceeding three times that of a single site (Bakos et al., 2013).

1.3.3 Follow-up observations

Once a planet host candidate is identified from the small aperture photometric observations, typically three kinds of follow-up observations are performed: reconnaissance spectroscopy to determine the spectral type of the star, photometry to refine transit parameters, and RV measurements to constrain the mass and orbital eccentricity of the planet.

Reconnaissance spectroscopy is typically done on a one meter class telescope. A small number of exposures are usually enough to tell a main sequence star from a giant. This is very important, because the object that causes a 1% deep transit in the lightcurve of a giant star is too large to be a planet. Many eclipsing binary systems can also be ruled out based on their double-lined spectra.

At this step, the spectral type of the star is also identified. Usually, A stars are not followed up, because they rotate fast, therefore their spectral lines broaden because of the Doppler shift, making them unsuitable to detect RV variations due to low-mass planets. However, evolved A stars eventually spin down, making them suitable targets. These “retired” A stars are among the most massive known planetary hosts (Johnson et al., 2007, 2008a; Bowler et al., 2010; Johnson et al., 2010a, b, 2011a, 2011b). Recently, however, it has been suggested that the evolved stars targeted in these surveys are not as massive as originally thought (Lloyd, 2011; Schlaufman & Winn, 2013).

Follow-up photometry can also be performed with a one meter class telescope. Surveying stars looking for transits would be too expensive with such telescopes, but once there is a transit prediction, high precision photometry is valuable in measuring the shape and depth of the transit, sometimes in multiple bands, to constrain orbital parameters. In addition, the shape of the transit might provide further confirmation to the planetary nature of the system, as those exhibit transits with mostly flat bottoms, as opposed to the V-shaped transits of eclipsing binaries.

Once a planetary candidate passes these two tests with modest size telescopes, the last step is to perform radial velocity measurements. Depending on the brightness of the target and the instruments available to the research group, this might be done with telescopes with an aperture diameter starting from a few meters. The amplitude of the RV variations is determined by the planet to star mass ratio, which together with the planet size inferred from transit depth allows for constraining the planetary density. In addition, RV data points help determine the orbital eccentricity and the argument of periastron. However, measurements are typically taken at quadrature, that is, at the predicted peak radial velocity values, because such data points will yield the best constaints on the mass ratio. Unfortunately, such observations are the least useful in determining the orbital eccentricity.

1.4 Planetary system characterization

Nature shows us only the tail of the lion. But there is no doubt in my mind that the lion belongs with it even if he cannot reveal himself to the eye all at once because of his huge dimension.

Albert Einstein (1914)

In this section, we highlight a few very specific aspects of characterizing transiting planetary systems.

1.4.1 Lightcurve models

A central theme in exoplanetary explorations is the inability to resolve the planetary system (except for planets at large orbital separation). Usually, all information we have about the planet and the host star has to be inferred from observations without spatial resolution. Measuring the lightcurve, that is, brightness as a function of time, plays an important role in the field of exoplanets, partially because brightness measurements are less costly than spectroscopic ones.

It is important to note that even though the observer cannot resolve the transited object, a good understanding of its spatial structure is crutial for generating an accurate lightcurve model.

Since the first observation of planetary transits (HD 209458, Henry et al., 2000; Charbonneau et al., 2000), there has been a need for fast and accurate modeling of transit lightcurves. Speed is an important issue, because parameters of the system (radius ratio, impact parameter, etc.) can be constrained by generating a large number of model lightcurves and comparing them to the observations. A popular and robust method for exploring parameter uncertainties and correlations is Monte Carlo Markov Chain (MCMC) calculations, which require an even larger number of model evaluation.

To model a transit lightcurve, one has to calculate not only what fraction of the planet is in front of the star at any given moment, but also account for the fact that the stellar radiation is not uniform across the stellar disk. Most photons we see come from a certain distance under the stellar surface, and at the center of the star, this distance corresponds to a larger depth than closer to the limb, where we see the surface at an angle. Since the star is hotter the deeper we go under the surface, and hotter material emits more, the center of the disk will be seen to be brighter, and the limb darker. This phenomenon is called limb darkening.

A pioneering paper for lightcurve models is written by Mandel & Agol (2002). They deal with the nonlinear limb darkening in the form

and its special case the quadratic limb darkening, when (Claret, 2000). Here is the intensity of the stellar disk as seen by the observer as a function of , the distance from the center, in stellar radius units, so that .

The main idea in their paper is essentially a change of variables: the stellar disk can be thought of a continuous family of concentric disks with radii ranging from zero to one, each with uniform intensity distribution, superimposed. For the smaller disks, the planet-star separation and the planetary radius are proportionally blown up when expressed in units of disk radii. This treatment allows the authors to evaluate the integral analytically, creating the fast lightcurve model that has became known as the Mandel–Agol model.

Building on this model, lightcurve calculations have been developed for other situations or applications. Pál (2008) presents analytical formulae for the partial derivatives of flux, which speeds up certain fitting algorithms by almost an order of magnitude. Pál (2012) investigates elaborate configurations of multiple transiting bodies. Abubekerov & Gostev (2013) derive efficient analytic and numerical formula for a large range of limb darkening laws.

Models that require numerical integration in one dimension are more computationally expensive than analytic ones, and models that require numerical integration in two dimensions are much, much more computationally expensive. In Chapter 6, we present a one-dimensional numerical integral for transits of spotted stars, a problem for which only two-dimensional numerical integral models existed before.

Chapter 4 applies some methods and understanding from the domain of planetary transits to a different scenario: stars transiting accretion disks around supermassive black holes in galactic centers. An important difference from the usual treatment of planetary transits is that in this case, the source has a very complicated structure, so most computation time is spent on generating the image as the observer would see the accretion disk unocculted.

1.4.2 Obliquity

In our Solar System, planetary orbits are aligned with the equator of the Sun within 7. The angle between the orbit of an exoplanet and the equator of its host star is commonly called obliquity, and from the Solar System, one might guess that it is typically low.

The most popular method for measuring obliquity is the Rossiter–McLaughlin effect (Rossiter, 1924; McLaughlin, 1924). It exploits the redshift of the receding half of the star and the blueshift of the approaching half, to infer obliquity from time-resolved spectroscopic measurements during transit. It follows from the nature of this phenomenon that it requires a large telescope for taking exposures shorter than the transit, that its sensitivity depends on the rotation rate and the inclination of the host star and the transit impact parameter, and that it in fact measures projected obliquity (that is, projected onto the sky plane).

The first observation of the Rossiter–McLaughlin effect on a transiting exoplanet is described by Queloz et al. (2000), who find that HD 209458 has a low projected obliquity with respect to the orbit of HD 209458b, as expected by extrapolation of the Solar System observations. However, later observations unveiled a number of planets on polar and retrograde orbits, which has important implications on planetary formation and migration theories. Winn et al. (2010a) review projected obliquity measurements, find a correlation with host star effective temperature, and explore possible explanations.

1.4.3 Eccentricity

Measurement results of most physical parameters, like planetary radius or semi-major axis, are often given as a Gaussian posterior distribution. Even if the underlying distribution is log-normal, a normal distribution usually provides an adequate approximation as long as the relative uncertainty is small. In such cases, the result can be parametrized with a single error bar, like . In some cases, the distribution is skewed, and characterization like is customary.

However, orbital eccentricity is usually more complicated: for bound orbits, , so the distribution is truncated. Also, highly eccentric orbits are rare, so the distribution is usually skewed. Last, but not at least, close systems often circularize due to tidal interaction, yielding in a bimodal eccentricity distribution: for a positive fraction of objects, and for the rest of them.

Since it is an important question to decide whether a planetary orbit is circular, sometimes two models are investigated: one where the orbit is restricted to be circular, and another with two extra parameters: eccentricity and argument of periastron (observer-star-periastron angle) . The two models can be compared by statistical methods, for example, using the Bayesian Information Criterion (BIC, see e.g., Schwarz, 1978), or calculating the Bayesian evidence (see e.g., Feroz et al., 2009).

However, occasionally only the eccentric case is analyzed, and the resulting eccentricity is quoted with a single error bar. Lucy & Sweeney (1971) investigate this case and provide a method to estimate the posterior likelihood of a circular orbit.

When evaluating the eccentric model, there is still a need to quantify the uncertainty of eccentricity. Instead of citing the maximum likelihood value and standard deviation for itself, a better proxy is to give the results in terms of the orbital parameters and (Ford, 2005). These describe the case seamlessly, without either parameters being truncated at this limit. Such characterization has been used for planetary orbits, for example, by Bakos et al. (2009) and Bakos et al. (2010).

Finally, it is important to note that when performing Monte Carlo Markov Chain (MCMC) calculations, this method in fact provides a biased prior for the eccentricity: in space, each eccentricity value has a measure proportional to . Using and instead of and to describe and result in a flat prior in , which is preferable to one proportional to . In the context of exoplanet orbits, this characterization has been used, for example, by Anderson et al. (2011a); Gillon et al. (2011).

1.5 Atmospheres

With insufficient data it is easy to go wrong.

Carl Sagan (1980)

The first order approximation of transits assumes that the planet is a completely dark sphere. This treatment is sufficient for measuring planetary radius, impact parameter, and other characteristics based on the lightcurve of the star.

In reality, however, planets are not completely dark. They emit thermal radiation, and extinguish part of the starlight traveling through their atmosphere. We already encountered the first phenomenon when discussing direct imaging of young planets in the infrared. Just like those planets exhibit thermal radiation due to residual heat since their formation, planets on close orbits exhibit thermal radiation due to their elevated equilibrium temperature caused by the large stellar flux.

When discussing equilibrium temperature, it is usually assumed that heating from incoming stellar flux get efficiently redistributed throughout the entire planetary atmosphere, unless measurements indicate otherwise. In this case, the effective temperature of the planet can be calculated from that of the star according to the energy balance:

Here is the Stefan–Boltzmann constant, is the orbital radius (we assume a circular orbit in this derivation), and is the Bond albedo that tells us what fraction of the incoming electromagnetic radiation power is reflected by the planet. The left hand side of the equation is the total electromagnetic power emitted by the star times a geometric factor that tells us what fraction of that radiation falls on the planet times a factor that tells us what fraction is absorbed. The right hand side is the total power emitted by the planet.

However, many close planets are presumably tidally locked, and in this case, it is possible that not all incoming power get redistributed over the entire planetary surface. Therefore the dayside temperature, that is, of the side facing the star, is more elevated than the nightside temperature. See, for example, Perez-Becker & Showman (2013) for a detailed theoretical model of heat redistribution by the atmosphere, and Cowan & Agol (2011) for heat redistribution coefficient measurements based on photometric observations.

If the planet has a constant effective temperature on its entire surface, the only way to distinguish its thermal radiation from that of the star is observing the occultation of the planet by the star, also called a secondary eclipse. On the other hand, if there is a difference between the dayside and nightside temperature, that shows up as an orbital modulation of the lightcurve.

Depending on the wavelength of the observations, the thermal emission of the planet can be negligible, and the secondary eclipse depth can be dominated by the reflected light. Such observations allow direct measurement of the geometric albedo, which in turns reveals information about the atmosphere. For example, a cloud-free atmosphere is expected to have low albedo due to pressure-broadened absorption lines of neutral sodium and potassium (Sudarsky et al., 2000). On the other hand, high altitude clouds reflect off most incoming starlight, resulting in a high albedo value.

When performing multiband or spectroscopic observations of secondary eclipses, emission lines can also be detected (e.g., Knutson et al., 2008).

Another exciting method for probing planetary atmospheres is transmission spectroscopy: detecting the extinction of starlight as it travels through the atmosphere, as first theoretically described by Seager & Sasselov (2000). Since most scattering phenomena depend smoothly on wavelength, the extinction spectrum will in fact reveal absorption features, allowing us to learn about atmospheric composition. One way to think about this phenomenon is to consider that observed at a strong absorption feature, the atmosphere blocks light, therefore the planet effectively seems to be larger, resulting in a deeper transit. On the other hand, a flat transmission spectrum, like the one observed for GJ1214b (Kreidberg et al., 2014), might indicate a grey opacity source, suggesting the presence of clouds.

1.6 Stellar activity

The sun comes up just about as often as it goes down, in the long run, but this doesn’t make its motion random.

Donald E. Knuth (1969)

An idealized transit lightcurve model might assume a spherical, homogeneously radiating star. An idealized radial velocity model might assume narrow lines emitted uniformly from the stellar surface. In reality, limb darkening results in an axisymmetric intensity profile of the stellar disk, and spectral lines are broadened by temperature and by rotation. Stellar activity, which includes flares, granulation, and starspots, further complicates the lightcurve or spectrum of the star. In the rest of this section, we focus on the effect of starspots.

Stellar activity is linked to magnetic processes. In fully convective stars or stars with a convective envelope, the magnetic field interacts with material flow close to the surface. Stronger magnetic fields can effectively block convection in areas, resulting in dark features on the surface known as spots. The best studied starspots are those on the Sun, but observations show that spots can differ greatly in their size, lifetime, latitude, and other characteristics, even on a Sun-like star like EK Draconis (see, e.g., Strassmeier & Rice, 1998). These observations are based on separating spots on the approaching (blueshifted) side of the star from those on the receding (redshifted) side, a technique called Doppler imaging (Deutsch, 1958; Vogt et al., 1987; Rice et al., 1989).

Starspots rotate with the stellar surface, causing variability on the timescale of stellar rotation. While this is a nuisance, for example, when modeling planetary transits, one can actually use this variation to model spots (see, e.g., Kipping, 2012; Bonomo & Lanza, 2012), or to measure the stellar rotational period (see, e.g., McQuillan et al., 2013).

Spots on transiting planet hosts can actually be eclipsed by the planet, resulting in an anomalous rebrightening of the transit lightcurve. This effect has been first observed by Silva (2003), and has been extensively studied since. This phenomenon is the central point of Chapters 56.

Chapter 2 Discovery of the transiting exoplanet HAT-P-27b

This thesis chapter has been originally published as .


We report the discovery of HAT-P-27b, an exoplanet transiting the moderately bright G8 dwarf star GSC 0333-00351 (). The orbital period is d, the reference epoch of transit is (BJD), and the transit duration is d. The host star with its effective temperature

Chapter 3 High precision relative photometry with HST

This thesis chapter has been originally published as .


We present HST STIS observations of two occultations of the transiting exoplanet HAT-P-1b. By measuring the planet to star flux ratio near opposition, we constrain the geometric albedo of the planet, which is strongly linked to its atmospheric temperature gradient. An advantage of HAT-P-1 as a target is its binary companion ADS 16402 A, which provides an excellent photometric reference, simplifying the usual steps in removing instrumental artifacts from HST time-series photometry. We find that without this reference star, we would need to detrend the lightcurve with the time of the exposures as well as the first three powers of HST orbital phase, and this would introduce a strong bias in the results for the albedo. However, with this reference star, we only need to detrend the data with the time of the exposures to achieve the same per-point scatter, therefore we can avoid most of the bias associated with detrending. Our final result is a upper limit of 0.64 for the geometric albedo of HAT-P-1b between 577 and 947 nm.

3.1 Introduction

The effective temperature of Jovian extrasolar planets in close orbits is strongly influenced by irradiation from their host stars. It is the Bond albedo, defined as the reflected fraction of incident electromagnetic power, that determines how much of this irradiation contributes to the thermal balance of the planet. Atmospheric temperature, in turn, determines – for a given composition – the presence and position of absorbers, clouds and other structures determining the albedo. Once the albedo is inferred from observations, the goal is to find a self-consistent atmospheric model and temperature.

Even though it is the Bond albedo that can be used directly in the calculation of effective temperature, the geometric albedo is more accessible observationally. It is defined as the ratio of reflected flux at opposition (zero phase angle) to the reflected flux by a hypothetical flat, fully reflecting, diffusely scattering (Lambertian) surface of the same cross-section at the same position. The geometric albedo of Lambertian surfaces is at most one: for example, a fully reflecting Lambertian sphere has a geometric albedo of (Hanel et al., 1992). However, some surfaces reflect light preferentially in the direction where it came from (a phenomenon known as opposition surge), and thus can exhibit geometric albedos exceeding one (e.g. Hapke, 1990).

A simple way to measure albedo is to observe an occultation (secondary eclipse) of a transiting exoplanet and directly compare the brightness of the star only (while the planet is occulted) to the total brightness of the star and the planet near opposition (shortly before or after the occultation). Examples of other phenomena that can be used to constrain albedo but are not discussed in this paper are phase variations (e.g. Harrington et al., 2006), Doppler-shift of reflected starlight (e.g. Liu et al., 2007), and polarization of reflected starlight (e.g. Hough et al., 2006).

To mention specific examples, the exoplanet HD 209458b has mass , radius , orbital period days, and zero-albedo equilibrium temperature K (Torres et al., 2008). When calculating , perfect heat redistribution on the planetary surface is assumed. Rowe et al. (2008, 2009) observe the occultation of HD 209458b in the 400–700 nm bandpass with the Microvariability and Oscillations of Stars (MOST) satellite. They find a geometric albedo of , and infer a upper limit of 0.12 on the Bond albedo, indicating the absence of reflective clouds. Based on atmospheric models, this constrains the atmospheric temperature to between 1400 K and 1650 K. Normally, a cloud-free atmosphere exhibits low albedo due to the strong pressure-broadened absorption lines of neutral sodium and potassium (Sudarsky et al., 2000). However, Spitzer Space Telescope observations of HD 209458b at 3.6, 4.5, 5.8, and 8.0 m indicate water emission features, suggesting temperature inversion in the higher atmosphere, which hints that an unknown absorber is present at low pressure (Knutson et al., 2008; Burrows et al., 2007b).

Another well-studied example is HD 189733 (, , days, K, Torres et al., 2008). Grillmair et al. (2008) and Charbonneau et al. (2008) observe its occultations with the Spitzer Space Telescope. They find strong water absorption features, indicating the lack of temperature inversion in the atmosphere. Sing et al. (2011) perform transmission spectroscopy on this planet, and interpret the results as indicating high altitude haze. This would cause the planet to exhibit high geometric albedo in the visible. Rowe et al. (2009) report on MOST observations of this planet’s occultation, but cannot constrain the albedo due to the high activity of the host star. Note that HD 209458b and HD 189733b only differ by a few hundred kelvins in terms of , yet seem to exhibit very different atmospheres.

The subject of this work is HAT-P-1b, a transiting exoplanet with mass , radius , orbital period days (Johnson et al., 2008b), and zero-albedo equilibrium temperature K (Torres et al., 2008). This last value is between those of HD 209458b and HD 189733b. Indeed, Todorov et al. (2010) observe two occultations of HAT-P-1b with Spitzer and infer a modest temperature inversion in the atmosphere from the occultation depths at 3.6, 4.5, 5.8, and 8.0 m. In the light of these observations, constraining the geometric albedo of this planet in the visible and near infrared would be useful for better understanding its atmospheric structure, and for refining atmospheric models.

In this paper, we report on observations of two occultations of HAT-P-1b. Section 3.2 presents the details of the observations, flux extraction, and detrending. We describe how we calculate the geometric albedo in Section 3.3, with special attention to handling each uncertainty source, and comparing relative photometry results to those without using the reference star. We discuss our findings in Section 3.4.

3.2 Observations and data analysis

The exoplanetary host star HAT-P-1 is a member of the wide binary system ADS 16402, with projected separation at a distance of 139 pc. HAT-P-1, also known as ADS 16402 B, with is only 0.4 magnitude fainter than its binary companion ADS 16402 A of the same G0V spectral type. This allows for relative (differential) photometry to mitigate the effect of systematic errors.

3.2.1 Observations and data preprocessing

A proposal was accepted as GO 11069 to observe two occultations (secondary eclipses) of HAT-P-1b with the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) High Resolution Channel (HRC). However, ACS failed in 2007 January, before the observations would have been carried out, and HRC remains inoperational to date. Instead, another HST instrument, the Space Telescope Imaging Spectrograph (STIS) carried out the observations, as program GO 11617. Two occultations of HAT-P-1b were observed during two visits, including two orbits before and one after each occultation. Because of the brightness of the targets, spectroscopy was required to allow for reasonably long exposures. To capture both stars without very tight constraints on spacecraft orientation, we did not use a slit. Table 3.1 summerizes the details of the observations. With these settings, the largest electron count in each exposure was about half the well size, well below saturation. However, it is interesting to note that longer exposures would not have posed a problem in terms of linearity either: when STIS pixels get saturated, the excees charge bleeds into surrounding pixels with virtually no loss, and summing these pixel counts still results in a linear response (Gilliland et al., 1999).

visit 1 visit 2
HJD at beginning of first exposure 2 455 544.8269 2 455 888.6602
HJD at end of last exposure 2 455 545.1181 2 455 888.9497
number of orbits 5 5
number of spectra for each orbit 19+23+23+23+23 19+23+23+23+23
number of spectra in total 111 111
grating G750L
slit slitless
exposure time 100 s
cadence 128 s
subarray size 380x1024 pixels
gain 4
Table 3.1: HST/STIS program GO 11617 observation parameters

We identify hot pixels and exclude them from our apertures, and identify cosmic rays and substitute them by the average of the previous and next frame values. Then we perform rectangular aperture photometry on the two stars, and define an aperture in the entire length of the detector in the dispersion direction for sky background estimation. For each exposure, we subtract the background aperture photon count from the stellar aperture photon counts, scaled by the number of pixels. Figure 3.1 shows a typical exposure from each visit, with the apertures used. The blue end of the stellar apertures is 564 nm for the first visit and 557 nm for the second. The difference is caused by different orientations of the telescope around its optical axis, resulting in the STIS detector edge cutting the HAT-P-1 spectrum at different positions for the two visits. These wavelength values are calculated after identifying the H and Na I D lines in the stellar spectra. We extract the same spectral range for the two stars within a visit to fight wavelength-dependent systematics. However, we allow for different blue end cuts between visits, otherwise we would lose too many photons in the second visit due to the more restrictive wavelength limit of the first one. We do not expect the geometric albedo to vary significantly due to this small change in blue end wavelength cut.

Figure 3.1: The first spectrum used in the final analysis (the second exposure of the second orbit) of the first and second visits (left and right panels, respectively). The upper star is ADS 16402 A, the lower one is the planetary host ADS 16402 B. The rectangular apertures around the spectra are also shown. The bottom rectangle is the background aperture. Note the different cross-dispersion distance and dispersion direction shift between the two stars’ spectra for the two visits due to spacecraft orientation. The blue end of the stellar apertures is determined by where the detector edge cuts the specrum of ADS 16402 B.

3.2.2 Detrending

The next step is to detrend the data, that is, to mitigate instrumental effects by subtracting multiples of vectors describing circumstances of the observations. We try detrending with time (to remove the overall linear trend within a visit), HST orbital phase and its powers (to remove orbitwise periodic variations), CCD housing temperature (CCD chip temperature is not available with the current Side-2 electronics), a focus model provided by HST Observatory Support, fine pointing data available from telemetry, and fine pointing data based on a fit for the position of the spectra on the CCD. To detrend, we simultaneously fit for free parameters of the lightcurve model (reference flux and planet-to-star flux ratio) as well as detrending vector coefficients using a linear algebraic least square method. We actually use the magnitude of HAT-P-1 or the magnitude difference of the two stars, that is, we assume that instrumental effects are multiplicative. Since both the planet-to-star flux ratio and detrending corrections are very small, this is equivalent to assuming additive effects and fitting for the flux or flux ratio. Each detrending vector is mean subtracted so that they do not change the average stellar magnitude. To quantify the effect of detrending and avoid overfitting, we minimize the Bayesian Information Criterion (BIC), which is the sum of and a term penalizing extra model parameters (Schwarz, 1978).

When analyzing STIS data to perform photometry on exoplanetary host stars, Knutson et al. (2007) find it justified to detrend with a cubic polynomial of HST orbital phase, whereas Brown et al. (2001) and Sing et al. (2011) use fourth order polynomials. Indeed, if we only consider the lightcurve of HAT-P-1, we find the lowest BIC when detrending with time at mid-exposure and the first three powers of HST orbital phase. However, if we divide the lightcurve of the planetary host by that of the reference star ADS 16402 A, we find the lowest BIC when detrending only with mid-exposure time. This shows that relative photometry is less sensitive to systematics, and can mitigate HST orbital effects enough that detrending with orbital phase is not justified. The Bayesian Information Criterion does not justify detrending with the temperature, focus, or jitter vectors in either case.

Figure 3.2 presents the raw lightcurves (panels a–d), without detrending, which indeed show strong orbitwise periodic variations. Panels (e, f) show the lightcurve of HAT-P-1 detrended with time at mid-exposure and the first three powers of HST phase, demonstrating that this procedure corrects for the overall linear trend and most of the orbitwise periodic variation. On the other hand, the raw relative lightcurves shown on panels (g, h) do not exhibit such large variations, and we only need to remove the linear trend (panels i, j). These panels all have logaritmic vertical axes with the same scaling, so that relative scatter is directly comparable. Note that since we perform a simultaneous fit of the occultation lightcurve and detrending vectors, the resulting tells us how close the data are to a model accounting for both the occultation and systematics, without the danger of misinterpreting the occultation as scatter.

The strength of these observations is the presence of the reference star, which already proves to be advantageous. In order to quantify how much it improves the albedo limits, we perform a full analysis both without and with this reference data, independently tuning all extraction parameters.

3.2.3 Aperture parameters and data omission

If we only extract the flux of HAT-P-1, we find that we obtain the least scatter after detrending (with mid-exposure time and first three powers of HST orbital phase) if the stellar apertures are 14 and 16 pixels wide in cross-dispersion direction for the two visits, respectively, and the stellar spectra are cut off at 788 nm. Redward of this wavelength there is extra scatter due to fringing, which is difficult to combat in slitless mode. If we use the flux ratio of the two stars, we get the least scatter after detrending (this time with mid-exposure time only) if the stellar apertures are 16 and 18 pixels wide in cross-dispersion direction for the two visits respectively, and the stellar spectra are cut off further in the near infrared at 947 nm. Using the reference star thus allows us to extract photons from a larger aperture. The optimal background aperture is 27 pixels wide for the first visit and 22 for the second in both cases.

We estimate sky background using an aperture placed as far from the stellar spectra as possible, to avoid contamination by starlight. Varying the background aperture width by a few tens of pixels introduces scatter of 0.02 in the best fit albedo as long as the aperture is not too narrow and not too close to the stars either. Given the other uncertainty sources discussed in Section 3.3.2, this means that the results are practically insensitive to the exact choice of the background aperture. However, if we place the background aperture between the two spectra on the detector, or expand it on the side to get within 50 pixels of HAT-P-1, we find a larger average error, of 0.05, in the geometric albedo due to stray starlight.

Finally, we investigate whether it is justified to omit data points. For example, Sing et al. (2011) and Knutson et al. (2007) both omit the first orbit of each five orbit visit, also the first exposure of each subsequent orbit, because these data points exhibit larger scatter. The scatter of the first orbit might be attributed to the thermal settling of the spacecraft after its new pointing. We therefore calculate the scatter per data point for all twelve possible combinations of omitting the first orbit or not, omitting the first or first two exposures of each orbit or not, and omitting the last exposure of each orbit or not. Comparing these results, we find that both for the analysis of HAT-P-1b only and for relative photometry using the reference star, it justified to omit the first orbit of each visit and the first exposure of each subsequent orbit, but not more, consistently with Sing et al. (2011) and Knutson et al. (2007). We therefore keep 88 data points per visit. These data points are represented with filled circles on Figure 3.2, whereas omitted data points are represented with empty ones.

Figure 3.2: Panels (a, b): background-subtracted photon count per exposure for reference star ADS 16402 A, in million photons. Panels (c, d): same for planetary host star HAT-P-1. Panels (e, f): photon count of HAT-P-1 detrended with time at mid-exposure and first three powers of HST orbital phase, in million photons. The model lightcurve of a fully reflecting planet is overplotted with a solid line. Panels (g, h): relative flux, that is, background-subtracted photon count of HAT-P-1 divided by background-subtracted photon count of reference star. Panels (i, j): relative flux detrended with time at mid-exposure. The model lightcurve of a fully reflecting planet is overplotted with a solid line. Panels (k, l): the fraction of the planetary surface that is illuminated and unobscured, with the values for each exposure overlaid on the continuous curve. The bottom horizontal axes show mid-exposure time in HJD, the top horizontal axes show , the planetary orbital phase since midoccultation. Panels (a–j) have logarithmic vertical axes with the same scaling, so that relative scatter is directly comparable. Panels (a), (c), (e), (g), (i) and (k) display the first visit, (b), (d), (f), (h), (j), and (l) the second. Filled circles represent data points included in the analysis, empty circles the omitted ones.

Note that the data analysis parameter space has many dimensions: stellar and background aperture geometry, data omission, and choice of detrending parameters. The choices presented above were found after multiple iterations, and they yield the minimum residual sum of squares in each dimension, but strictly speaking, they might not represent the global minimum.

3.3 Upper limit on geometric albedo

3.3.1 Lightcurve model

To model the occultation lightcurve, we adopt the transit ephemeris, the planetary radius and orbital semi-major axis in stellar radius units, and transit impact parameter values and uncertainties reported by Johnson et al. (2008b). We assume that the orbit is circular, consistently with theoretical expectations, radial velocity measurements (Bakos et al., 2007a), and occultation timing (Todorov et al., 2010). We account for the light travel time of 55 s across the planetary orbit, and neglect the thermal radiation of the planet. We calculate the projected area of the unobscured part of the dayside of the planet by assuming that both the stellar and planetary disks are circular, and the terminator line on the planet is an arc of an ellipse. Let denote this area relative to the total planetary disk. Then the observed flux at any given time is , where is the stellar flux, and is the flux of the planet in opposition. Note, however, that we can never observe the theoretical maximum flux from the planet, because opposition happens during occultation (therefore is always smaller than one). We derive the following model for :


Here is the stellar radius, is the planetary radius, is the planetary orbital radius, is the planetary orbital phase, is the orbital phase at midoccultation, is the impact parameter, is the projected distance of the center of the planetary and stellar disks, and are auxiliary functions for the case of ingress and egress. The bottom panels of Figure 3.2 show as a continuous function of time, with the value in the middle of each exposure overlaid.

This model gives the exact result as long as the boundary of the stellar disk and the terminator line on the planet do not intersect in projection, and is extended in a monotonic and continuous manner at the end of ingress and beginning of egress when they do. This introduces a small error which only affects a few data points. Note, however, that there are other sources of error: when determining the terminator line, we assume that the planet is irradiated by parallel rays from the direction of the center of the star, instead of properly calculating irradiation in the belt from where the star is partially seen on the horizon of the planet. Also, we do not account for that the substellar point is closer to the star, therefore it is exposed to more irradiation. These errors are in the order of , that is, they bias our geometric albedo estimate by the order of 1%.

We assume Lambertian reflectance of the incident light. Note that the maximum angle of stellar irradiance on the planet to the line of sight during our science exposures used in the final analysis is . In case of reflection from materials like regolith, this would be very different from reflection in opposition, but we expect the reflected flux for a gas giant like HAT-P-1b to depend only weakly on the incident angle.

Using this model lightcurve, we then fit simultaneously for the planet to star flux ratio at opposition and the coefficients of the detrending vectors. Then the geometric albedo can be calculated using the expression (e.g. Rowe et al., 2008)


To illustrate the magnitude of the effect with respect to the scatter of the data, Figure 3.2 panels (e, f) and (i, j) feature a plot of a model lightcurve, assuming that the planet is a fully reflecting Lambertian sphere with a geometric albedo of .

3.3.2 Uncertainty sources

There are two sources of errors in the inferred geometric albedo: observational errors and contribution of parametric uncertainties of the system. We estimate the total observational uncertainty by bootstrapping. Since instrument parameters like temperature might not have a normal distrubtion or might not influence the data in a linear fashion, we cannot assume that the error distribution is normal. Therefore we apply the bias corrected accelerated bootstrap method (Efron, 1987), a generalized bootstrap algorithm that partially corrects for effects due to non-normal error distributions.

The total observational uncertainty is due to photon noise and other uncertainties due to instrument thermal instability, pointing jitter, and stellar activity in both stars. (Quantization noise due to the gain set to 4 can be neglected.) We estimate the uncertainty due to photon noise by independently redrawing every data point (HAT-P-1, reference star, and background for each exposure) from a Poisson distribution with a parameter given by the original photon count, and recalculating the geometric albedo. Finally, we get the estimate of the uncertainty due to other sources by subtracting the photon noise estimate from the total observational uncertainty estimate in quadrature.

When quantifying the uncertainty contributions of the planetary system parameters, we divide them into two groups: ephemeris (mid-transit time and period), and geometry (planetary radius and orbital semi-major axis relative to the stellar radius, and impact parameter). This division is important for two reasons: first, the two visits took place 1181 and 1525 days after the reference mid-transit ephemeris of Johnson et al. (2008b), respectively, therefore we expect that the second visit will suffer more from the uncertainty in mid-occultation time. Second, it is easier to refine the ephemeris by photometric follow-up transit observations than to refine the geometric parameters, so we want to assess how much this would improve the results. Note that geometric parameters not only influence the occultation lightcurve shape that we use for fitting, but also factor in when converting the planet-to-star flux ratio to geometric albedo. We estimate the uncertainties in the geometric albedo due to uncertainties in ephemeris, due to those in geometric parameters, and their total contribution due to all parametric uncertainties by redrawing the respective parameters from independent normal distributions defined by their best fit values and uncertainties, and recalculating the geometric albedo in each case. We confirm that and add up in quadrature to . We use 1 000 000 bootstrap iterations to estimate the observational uncertainty , and 1 000 000 random drawings to estimate each of , , , and .

Finally, we add the observational and parametric uncertainty estimates in quadrature to estimate the total uncertainty , and add this to the best fit geometric albedo to calculate the upper limits. The breakdown of uncertainty sources is represented in a tree structure in Tables 3.2 and 3.3: the upper limit is the sum of its children nodes best fit and total uncertainty, and uncertainties are quadrature sums of their children nodes.

We find that most lower albedo limits of the final analysis are negative, therefore we only present upper limits for the geometric albedo. The and upper limits are determined so that the probability of the geometric albedo being smaller than this is and , respectively. As a comparison, the one-sided and upper limits of a normal distribution occur at and above the mean, respectively, and we expect our corresponding uncertainties to have a similar ratio.

Table 3.2 presents the best fit values, uncertainties due to various phenomena, and upper limits for the geometric albedo, based on the lightcurve of HAT-P-1 only. We performed the calculations for the two visits separately, and also for a joint model, where a single geometric albedo value was fit for the two visits, but we allowed for different coefficients of the detrending vectors for the two visits. The upper limits highlight the weakness of detrending: the upper limit based on the first visit data is meaningless (greater than one), and even worse, the limit based on the second visit data only is unphysical (negative). The reason for this is that the detrending vectors are not orthogonal to the occultation signal, that is, systematic effects have a component that mimics the occultation lightcurve. Therefore even though detrending is justified by the Bayesian Information Criterion, it introduces a bias in the geometric albedo values.

On the other hand, if we feed the flux ratio of the planetary host star HAT-P-1 and the reference star ADS 16402 A into the occultation lightcurve model, much less detrending is justified, thus we expect less bias in the result. Indeed, Table 3.3 shows that the results from the two visits are much closer to each other, and all upper limits are positive. Even though the uncertainties in the two cases are very similar, as we can tell by comparing values in Tables 3.2 and 3.3, in the second case, this is achieved by using only one detrending vector instead of four. This shows the enormous advantage of the reference star: to diminish the need for detrending, therefore arrive at similar uncertainties with much less bias. We adopt the upper limits of the joint fit as our final result.

visit 1 visit 2 joint fit
upper limit 1.44 1.94 0.03
– best fit
0.19 0.69 0.57
0.19 0.68 0.57
0.15 0.52 0.36
0.12 0.44 0.44
0.03 0.11 0.05
0.00 0.00 0.04
0.03 0.11 0.03
Table 3.2: Best fit values, uncertainties due to different sources, and upper limits for the geometric albedo, based on HAT-P-1 lightcurve only, for comparison
visit 1 visit 2 joint fit
upper limit 0.08 0.59 0.55 1.14
– best fit
0.20 0.71 0.23 0.83 0.15 0.54
0.19 0.71 0.23 0.83 0.15 0.53
0.15 0.57 0.16 0.56 0.11 0.39
0.12 0.42 0.17 0.61 0.10 0.36
0.03 0.10 0.01 0.07 0.01 0.07
0.02 0.08 0.00 0.05 0.01 0.05
0.01 0.05 0.01 0.06 0.01 0.05
Table 3.3: Best fit values, uncertainties due to different sources, and upper limits for the geometric albedo, based on the lightcurves of HAT-P-1 and reference star ADS 16402 A, our final results

By comparing the different contributions to the uncertainty of the geometric albedo, we see that the observational uncertainties are much larger than the parametric ones. This means that performing further photometric observations of occultations could significantly improve the albedo upper limit, whereas using additional transit observations to refine the ephemeris and geometric parameters and use them to reanalyze these data would not. As for the observational uncertainties, residual systematic uncertainty and photon noise contribution are within a factor of two, which tells us that further improving our data analysis could push the upper albedo limits down only by a small amount. Surprisingly, we do not always find to be larger for the second visit than for the first, but we confirm that and uncertainties of each kind have a ratio close to what is expected for a normal distribution.

3.4 Discussion and summary

Based on HST STIS observations, we established as the , and as the upper limit for the geometric albedo of HAT-P-1b in the 557–947 nm band. Unfortunately, this limit is not tight enough to determine whether there is temperature inversion in the atmosphere. This question is relevant because HAT-P-1b has an equilibrium temperature between that of HD 209458b (thought to exhibit temperature inversion) and HD 189733b (thought not to). In addition, a better constrained albedo would provide information about the actual atmospheric temperature of the planet, as well as indicate the presence or absence of reflective clouds or high-altitude haze.

Our data analysis demonstrates that the reference star ADS 16402 A helps us greatly to reduce systematic effects. Even though detrending with powers of HST orbital phase would equally reduce scatter in the signal (just like demonstrated by Brown et al., 2001; Knutson et al., 2007; Sing et al., 2011), we find that it introduces a bias in the geometric albedo estimate. We attribute this effect to the fact that these detrending vectors are not orthogonal to the occultation signal. Most of this bias can be avoided by performing relative photometry, so that much less detrending is necessary to mitigate systematic effects. This is possible because ADS 16402 A has a similar brightness and same spectral type as HAT-P-1, and their angular distance is small enough so that they fit in the STIS field of view, but large enough so that their PSFs do not overlap.

We found that the uncertainties of the system parameters have a negligible effect on the geometric albedo uncertainty. The dominant uncertainty sources are photon noise and other noise effects (thermal instability, pointing jitter, other systematics, and astrophysical noise of the two stars), the contributions of which scale inversely with the square root of the number of observations. This also means that additional observations would improve the upper limit, without being too limited by how precisely we know the geometry and ephemeris of the planetary system. For example, if the geometric albedo of HAT-P-1b was 0.1, then approximately seven times more observations would be required to arrive at 0.4 as a upper limit, enough to infer the absence of an omnipresent reflective cloud layer.

It is interesting to note that these observations were originally proposed for the grism instrument ACS/HRC, which has a total throughput of 0.15–0.25 in this wavelength range, as opposed to 0.04–0.08 for STIS with the G750L grating. Thus ACS/HRC observations would have resulted in roughly three times more photons. Assuming the same and values, this approximately translates to a of instead of for the limit, and instead of for .

The binary companion star can help data analysis of further HAT-P-1 observations (e.g., Wakeford et al., submitted). Also, similar methods could be used for other planetary hosts in binary systems. A suitable example is XO-2, of magnitude , with the companion XO-2 S of at 31” separation. This companion star has been used as a reference for transmission spectroscopy both from the ground with GTC (Sing et al., 2012), and with HST NICMOS (Crouzet et al., 2012). XO-2 is a potential target for relative photometry during occultation with HST STIS in slitless mode: the day orbital period of XO-2b (Burke et al., 2007) would mean a larger planet-to-flux ratio than in case of HAT-P-1b for the same geometric albedo, and its zero-albedo equilibrium temperature K (Torres et al., 2008), similar to that of HAT-P-1b, would make such a measurement interesting in terms of atmospheric models.


M.J.H. and J.N.W. gratefully acknowledge support from NASA Origins grant NNX09AB33G. G.Á.B. acknowledges support from NSF grant AST-1108686 and NASA grant NNX12AH91H.

Chapter 4 Modeling prospective AGN transits

This thesis chapter has been originally published as .


Supermassive black holes (SMBH) are typically surrounded by a dense stellar population in galactic nuclei. Stars crossing the line of site in active galactic nuclei (AGN) produce a characteristic transit lightcurve, just like extrasolar planets do when they transit their host star. We examine the possibility of finding such AGN transits in deep optical, UV, and X-ray surveys. We calculate transit lightcurves using the Novikov–Thorne thin accretion disk model, including general relativistic effects. Based on the expected properties of stellar cusps, we find that around solar mass SMBHs, transits of red giants are most common for stars on close orbits with transit durations of a few weeks and orbital periods of a few years. We find that detecting AGN transits requires repeated observations of thousands of low mass AGNs to 1% photometric accuracy in optical, or in UV bands or soft X-ray. It may be possible to identify stellar transits in the Pan-STARRS and LSST optical and the eROSITA X-ray surveys. Such observations could be used to constrain black hole mass, spin, inclination and accretion rate. Transit rates and durations could give valuable information on the circumnuclear stellar clusters as well. Transit lightcurves could be used to image accretion disks with unprecedented resolution, allowing to resolve the SMBH silhouette in distant AGNs.

4.1 Introduction

Photometric transits have been successfully used to identify and characterize transiting extrasolar planets since the discovery of the first one, HD 209458b (Charbonneau et al., 2000; Henry et al., 2000). Transit shape is a telltale of planetary, orbital, and stellar parameters. Moreover, the apparent time-dependent redshift of the star due to the planet covering its receding or approaching side during transit can reveal the projected angle between the planetary orbital axis and the stellar rotational axis (Rossiter–McLaughlin effect, Rossiter, 1924; McLaughlin, 1924). These methods show that transits are very powerful in probing planetary systems.

Similarly, active galactic nuclei (AGN) accretion disks can be probed by observing occultations in X-ray due to broad line region clouds: optically thick clouds orbiting the supermassive black hole (SMBH) in the region from where broad emission lines of the AGN are thought to originate. These occultations have a large covering factor of to (see e.g. McKernan & Yaqoob, 1998; Turner et al., 2008; Bianchi et al., 2009; Maiolino et al., 2010; Risaliti et al., 2007, 2009a, 2009b, 2011b). Risaliti et al. (2011a) pointed out that analogously to the Rossiter–McLaughlin effect, temporally resolved spectroscopic observations of an eclipse could be used to probe the apparent temperature structure of the accretion disk and the origin of the iron K emission line, allowing to constrain the black hole spin and inclination.

In this paper, we examine the possibility of stars on close orbits transiting their host AGN. There are multiple coincidences that make it possible to detect such events: the radius of large stars is comparable to the characteristic size of an accretion disk around a SMBH, resulting in transits deep enough to detect. The orbital period is a few years in the innermost regions where the stellar number density is highest, short enough for repeated observations. Finally, the typical transit duration for these innermost stars is an hour to a few weeks (depending on the stellar population and the observing band), which is longer than the typical cadence required to observe AGNs, but still accessible on human timescales.

In order to produce a detectable signature in the lightcurve, the transiting object has to be a main sequence O or B star, a Wolf–Rayet (WR) star, an AGB star with a surrounding dust cloud (OH/IR star), a young supermassive star, or an evolved giant. Late type main sequence stars are too small to cause a transit detectable from the ground in an AGN with black hole mass , unless they are “bloated” by irradiation of the AGN (Lohfink et al., 2012).

High photometric and temporal resolution observations of AGN stellar transits have the potential to map the accretion disk structure with an unprecedented accuracy. Stars are optically thick in all electromagnetic bands, and unlike broad line region clouds (e.g. Maiolino et al., 2010), they have a simple spherical geometry, making it easier to interpret lightcurves and to reconstruct the image of the accretion disk along the transit chord (projected stellar trajectory). Furthermore, such transits offer unique observations of individual stars in distant galaxies. Detections of multiple transits could reveal valuable information on the number density of stellar cusp central regions.

In this paper, we calculate transit depths, durations, periods, rates, and instantaneous probabilities based on physical models of nuclear stellar clusters. We also present simulated transit lightcurves and transit depth maps in multiple frequency bands. The shape of the lightcurve depends on the observing wavelength, the mass and spin of the SMBH, the accretion rate, the inclination of the accretion disk, the projected position of the transit chord, and the orbital velocity of the transiting object. Therefore observations with sufficient photometric accuracy and time resolution allow us to constrain these parameters, and to test the employed accretion and general relativity models.

In Section 4.2, we review observations of large stars closely orbiting Sgr A*, the SMBH at the center of the Galaxy (§4.2.1), and summarize theoretical models of semi-major axis distribution and mass segregation (§4.2.2). We set up models of the stellar population and radial distribution (§4.2.3), and determine the minimum (§4.2.4) and maximum (§4.2.5) orbital radii for transits. We present accretion disk thermal emission models and ray-tracing simulations in Section 4.3, showing transit spectra (§4.3.1), transit depth maps (§4.3.2) and lightcurves (§4.3.3). We calculate transit probabilities in Section 4.4. In Section 4.5, we determine the local density of AGNs of interest (§4.5.1), and study the feasibility of transit detections with optical (§4.5.2, §4.5.3) and X-ray instruments (§4.5.4). The most important simplifying assumptions, caveats, and implications are discussed in Section 4.6. Finally, we conclude our findings in Section 4.7.

4.2 Nuclear stellar clusters

In this section we review observations and models of stellar distribution in galactic nuclei.

Let denote the mass of the SMBH, and define and the gravitational radius . Then


where . For non-spinning black holes, the Schwarzschild radius is , and the radius of the innermost stable circular orbit is . For maximally spinning black holes, , and for prograde orbits in the equatorial plane.

4.2.1 Stellar populations

Many galaxies host a dense stellar cusp in their nucleus. Seth et al. (2008a) find that galaxies with a massive nuclear cluster are more likely to be active. Stars captured and transported inwards by the accretion disk may serve to fuel the AGN (Miralda-Escudé & Kollmeier, 2005).

In the Galaxy, Eisenhauer et al. (2005) report the orbital parameters of six B type main sequence stars orbiting the central SMBH on highly eccentric orbits with semi-major axes less than . These stars may be remnants of stellar binaries tidally disrupted by the SMBH, as first proposed by Hills (1988) (see also Yu & Tremaine, 2003). Candidates for the binary counterparts ejected with high velocity have been identified by Brown et al. (2009).

Bartko et al. (2009) find more than a hundred O and WR stars further from the Galactic Center, within projected distance. About one half of them is orbiting the central SMBH in a disc-like structure (the so-called clockwise disc), while the other half is on apparently randomly oriented orbits. These stars could have formed in a massive self-gravitating gaseous disk (e.g. Levin, 2007; Hobbs & Nayakshin, 2009), or formed further away and captured in close orbits by the Hills (1988) mechanism or by a cluster of stellar mass black holes (Alexander & Livio, 2004). Young stars could also be delivered to this region by an infalling globular cluster (Tremaine et al., 1975; Milosavljević & Merritt, 2001).

Most main sequence stars in the vicinity of a SMBH eventually evolve into red giants or supergiants, our candidates for transiting the AGN. Many such giants have been observed within 1 pc of the Galactic Center (e.g. Genzel et al., 2010; Paumard et al., 2006; Bartko et al., 2009, 2010). The fraction of stars in the giant phase within a stellar population depends strongly on its initial mass function (IMF) and formation history.

A more exotic possibility is the formation of supermassive stars in the accretion disk as suggested by Goodman & Tan (2004). Such a star would form outside 1000 , have a radius of approximately in case of solar metallicity, and would radiate at its Eddington limit, being luminous enough to have a detectable optical photometric signature not only when it transits the AGN but also when it is occulted behind it.

Observations of the Galactic nucleus show that the innermost cluster of young stars (so-called S-stars) is isotropically distributed (Genzel et al., 2010), as predicted by theoretical arguments. Even if stars are formed on an orbit coplanar with the accretion disk, Rauch & Tremaine (1996) provide a relaxation mechanism that could rapidly randomize the orbital orientations. The possible presence of intermediate mass black holes may help catalyze this process (Gualandris & Merritt, 2009). And even if stars cluster in disks, this coherent behaviour averages out when observing multiple galactic cores as long as the stellar disks and the accretion disks have independent orientations. Therefore we assume an isotropic distribution of stellar orbits in the cluster for the purpose of our probability estimates.

4.2.2 Density profile

We assume circular orbits for simplicity, and denote orbital radius by . A star on an eccentric orbit with the same semi-major axis would produce a transit of the same depth, with a transit probability and transit length within a factor of order unity.

Bahcall & Wolf (1976) showed that the equilibrium spatial number density distribution of a stellar cluster around SMBH is proportional to with if the stars in the cluster have the same mass. Analytical and numerical investigations of multiple mass populations show that the distribution for each mass bin is still likely to follow a power law. The value of is predicted to be smaller (down to ) for lower mass stars (Bahcall & Wolf, 1977; Freitag et al., 2006; Hopman & Alexander, 2006). For massive stars representing a small mass fraction in the stellar cluster, can be between 2 and (Alexander & Hopman, 2009), or as large as 3 (Keshet et al., 2009). Observations of the Galaxy by Bartko et al. (2010) show that WR/O stars from a distance of 30 mpc out to 600 mpc form an interesting anisotropic structure called the clockwise disk in the Galacit nucleus. These stars are distributed with a density exponent , while the B star population from 30 mpc to 1 pc exhibits . Note, however, that the age of main sequence O stars and WR stars is less than the two-body relaxation time at (O’Leary et al., 2009), therefore they are not expected to represent the equilibrium state. The observed mass distribution of solar mass stars in the Galactic nucleus is fit by a broken power law with and inside and outside of , respectively (Schödel et al., 2007).

To estimate the total number of stars, we first define the radius of influence as the radius of the sphere centered on the SMBH containing a total mass of in stars and stellar remnants. In case of a singular isothermal sphere (), this equals to , where is the velocity dispersion of stars further than (Merritt, 2004). In order to get an estimate of the number of stars, we set , independently of .

Using the relation now allows us to express the radius of influence as a function the supermassive black hole mass only, in the form


Here and depend on the coefficients of the relation. In particular, , where is the slope of , as defined by e.g. Tremaine et al. (2002). For example, the best fit of Tremaine et al. (2002) (with ) results in and , and the best fit of Ferrarese & Ford (2005) (with ) results in and . For our numerical results, we adopt the best fit values of Gültekin et al. (2009) (with ): and (but use 1 pc to normalize in our parametric expressions).

Now let us consider the stellar population in the vicinity of the SMBH. We assume that this population contains a species of stars with mass , having a radius large enough to produce detectable transits. We also assume that the number density of these stars follows the power law . Let denote the mass fraction of these large stars within the sphere of influence relative to the total mass of all stars and stellar remnants.

Typically the inner cutoff radius for the stellar distribution is much smaller than the radius of influence (see § 4.2.4 for numerical justification). Under this assumption, the spatial number density of the stars in question as a function of orbital radius is


Note that this argument has two weaknesses: first, the relation has a relatively large scatter. For a given SMBH mass, the intrinsic scatter of the bulge velocity dispersion is dex (Gültekin et al., 2009). Second, we applied results for the isothermal sphere to power law distributions with different exponents. This limits the accuracy of the transit rate estimates presented in Section 4.4.

4.2.3 Adopted models

(mpc) (mpc) (mpc) (mpc) (mpc) (mpc) (yr) (yr)
O stars 11 51.8 30 2.5 0.006 0.16 0.16 11 320 11 0.57 340
red giants 110 518 1.5 1.75 0.16 0.54 0.54 22 000 320 320 3.7 53 000
O stars 11 5.18 30 2.5 0.013 0.40 0.40 11 1 075 11 0.75 110
red giants 110 51.8 1.5 1.75 0.35 1.3 1.3 22 000 1 075 1 075 4.2 100 000
O stars 11 0.518 30 2.5 0.028 1.0 1.0 11 3 600 11 1.0 34
red giants 110 5.18 1.5 1.75 0.76 3.0 3.0 22 000 3 600 3 600 4.9 210 000
Table 4.1: Orbital radius and period limits for AGN transits for the adopted stellar population models ( for each model)

We consider two simple models for the transiting stellar populations around AGNs, summarized in Table 4.1.

First, we assume a young stellar population, motivated by the observations of a young population with an extremely top-heavy initial mass function of mean stellar mass following a density profile with in the central parsec of the Galaxy (Bartko et al., 2010). Thus in our model, we assume that a fraction of the stars are O type, with stellar mass , radius , and . This exponent is consistent with observations and theoretical predictions reviewed in Section 4.2.2. We assume that these O stars constitute a mass fraction of the population. This is consistent with the estimated total mass of WR/O stars in the Galactic Center if we consider that these stars are confined in the center part of the sphere of influence. However, note that this mass fraction depends sensitively on recent star formation rate and initial mass function of stars in the neighborhood of the SMBH.

For the second model, we consider an evolved, relaxed population of stars, and optimistically assume that a mass fraction of them are giants (or main sequence stars otherwise enlarged, like “bloated” or surrounded by a dust cloud), which are large enough to produce detectable AGN transits. For these giants, we assume , , and a Bahcall–Wolf equilibrium value of (see Section 4.2.2).

Note that depending on the star formation history, O stars and red giants might coexist in the cusp, in which case their contributions to transits add up.

4.2.4 Minimum orbital radius

The inner orbital radius cutoff of the stellar distribution, , is set by gravitational wave inspiral, and tidal and collisional disruption. While gravitational wave inspiral is the limiting factor for compact objects (Peters, 1964), tidal or collisional disruption sets a tighter constraint for stars.

The tidal disruption radius is


where ranges from to depending on the compressibility and polytropic mass distribution index of the star (Diener et al., 1995). We adopt the moderate value . The tidal disruption radius is given here normalized for a massive main sequence star.

Collisional disruption might be responsible for the depletion of luminous giant stars in the inner 80 mpc of the Galactic Center (Alexander, 1999). Following Hopman et al. (2007), one can write that the rate at which collisions decrease stellar density is


where is the collisional cross-section, and on average, collisions are required to disrupt a star (Freitag et al., 2006). We define the radius limit for collisional disruption as the radius where the stellar density -folds in time if the initial value is given by Equation 4.3:


We set the timescale to be : this is within an order of magnitude of both the main sequence lifetime of early type stars, and the lifetime of the giant phase for less massive stars. Note that here we only consider collisions within the large species, not with other, much smaller stars, which are less likely.

The minimum orbital radius is


Table 4.1 lists the minimum and maximum orbital radii for the two stellar species and three different SMBH masses. We find that in every case, collisions set a tighter constraint than tidal disruption for the potentially transiting stars.

Table 4.1 also lists the Keplerian orbital periods for stars at the minimum and maximum radii, which can be calculated as


We find that the closest main sequence stars have orbital periods of approximately one year, making repeated transit observations feasible. However, collisions set a larger radius limit for giants, resulting in longer orbits. This inner radius limit depends on the number fraction of giants: smaller implies less frequent collisions and thus allows closer orbits. However, a smaller value for would also mean smaller transit probabilities for a given AGN (see Section 4.4 below).

4.2.5 Maximum orbital radius

For the maximum orbital radius of stars capable of producing a transit, we have to consider two factors: gravitational microlensing due to the star, and the validity range of the presumed number density power law.

Gravitational microlensing caused by the transiting star can bend the light rays of the AGN which may significantly distort the transit lightcurve (Paczynski, 1986). This can happen if the transiting object is farther from the AGN than the radius at which the Einstein radius equals the apparent angular radius of the transiting object:


Therefore we restrict our transit probability calculations to the contribution of stars within orbital radius .

Note that this is a different configuration than a galaxy microlensing a distant quasar, which can also be used to probe the spatial structure of accretion disks around AGNs (e.g. for the case of Q2237+0305, see Kochanek, 2004, and references therein).

The power law distribution discussed above only applies to the stellar populations within the radius of influence from the SMBH, where its gravitational field dominates. In this paper, we do not investigate the distribution of stars outside the sphere of influence, but conservatively ignore their contribution to transit probabilities.

The maximum orbital radius is thus the smaller of the microlensing radius and the radius of influence:


Table 4.1 shows that typically microlensing is the limiting factor among main sequence O stars, whereas the radius of influence limits transits of the red giants. The reason for this is that the factor that determines the microlensing limiting radius accoding to Equation (4.9) is very different for these two species of stars.

4.3 Transit observables: spectra and lightcurves

Next we derive the AGN spectra and the transit observables.

The AGN luminosity is bounded by the Eddington limit (Shapiro & Teukolsky, 1983). We assume an Eddington ratio of as our fiducial value, following (Kollmeier et al., 2006; Shankar et al., 2012). Then the AGN luminosity is


4.3.1 AGN spectra

We adopt the general relativistic, radiatively efficient, steady-state thin accretion disk model of Novikov & Thorne (1973). This model describes the flux of thermal radiation from the disk as a function of radius in Equation (5.6.14b) (see Page & Thorne, 1974, for the explicit form of ). We assume no radiation from within . In addition to the thermal component, AGN spectra typically exhibit emission lines, excess infrared radiation from dust reprocessing UV emission, and a hard X-ray component usually assigned to inverse Compton scattering in a hot corona (Haardt & Maraschi, 1993). We do not account for these phenomena, but choose our observing bands so that their effect is minimal: an observation window around 200 eV is low enough so that the thermal component dominates over coronal emission, but it is higher than helium Lyman absorption and detector lower energy limits. It is important to note that little is known about the geometry of the corona, and simultaneously observing a stellar transit in hard X-ray might provide feedback to the models.

We follow the accretion disk photosphere model described by Milosavljević & Phinney (2005) and Tanaka & Menou (2010): we assume that the dominant absorption mechanism is the bound-free process, with an opacity that depends on the frequency and temperature. We assume that the temperature and thus the absorption opacity are constant down to an optical depth of one (the “bottom”) in the photosphere. We can calculate the total flux from the absorption and electron scattering opacities and the photosphere bottom temperature, using Equations (A13–A15) and (A17) of Tanaka & Menou (2010), but using the relativistic angular frequency given by Novikov & Thorne (1973) instead of Keplerian velocity. However, as the flux is known and the temperature is sought for, we have to use a simple iterative process to solve this implicit equation for the temperature. We find that usually electron scattering dominates the total opacity, but in the hottest parts of the accretion disk, absorption takes over. Note that the functional form of the specific flux differs from a blackbody spectrum, because the absorption opacity depends on frequency. For simplicity, we assume that the emerging radiation is isotropic in the frame comoving with the accretion disk.

Given the specific intensity of the accretion disk as a function of radius and frequency, the observed spectrum is determined by Doppler shift, gravitational redshift, and gravitational lensing. To account for these effects, we apply the transfer function method as described by Cunningham (1975) and implemented by Speith et al. (1995)111available at We calculate the radiative efficiency as a function of spin as described by Shapiro et al. (1983) to convert luminosity to mass accretion rate, which is the input parameter of this code. We fix the inclination angle (the angle between the spin axis and the line of sight) at , so that . The calculated spectrum of the accretion disk is shown in Figure 4.1 for a non-spinning black hole (with dimensionless spin ) in the top panel, and for a highly spinning black hole with a prograde disk with a conservative value in the bottom panel. The frequency is in the source rest-frame accounting for gravitational redshift, but a possible cosmological redshift for distant AGNs is not considered. The specific luminosity value displayed here is what an isotropic source would have to have in order to produce the same flux as the AGN does at this specific inclination.

Transit depth is defined as the blocked flux to out-of-transit flux ratio in a given band. Therefore the transit depth is between zero and one: zero if the transiting object does not cover any part of the accretion disk; one if the object completely blocks radiation (in which case it is called an occultation or eclipse). The transit depth varies with frequency and the location of the transiting object in projection. To be able to efficiently calculate transit depths at different positions, we implement a linear approximation to the radius–gravitational redshift grid generated by the above code to determine these values for a light ray parametrized by its projected position far from the AGN. Then we employ high order numerical approximation222code available from to integrate over the stellar disk in the projection plane. This gives the blocked specific flux, which we then divide by the total specific flux to obtain the narrow-band transit depth.

Figure 4.1: Spectrum and transit depth of an accretion disk around a black hole. Solid red curve: equivalent isotropic luminosity per logarithmic bins of frequency. Dashed blue curve: maximum possible narrow-band transit depth caused by a star with as a function of source frequency. Top (bottom) panel presents the case of a Schwarzschild (Kerr) BH with spin ().

In addition to AGN spectra, Figure 4.1 also depicts the maximum possible narrow-band transit depth caused by an early type main sequence star with as a function of frequency. The maximum depth of a general transit can be smaller if the star does not transit the most luminous part of the projected accretion disk. At high frequencies, the most luminous region is more compact, therefore the transit is deeper. The transit depth becomes constant at infrared frequencies less than the peak frequency of a black body spectrum with temperature of the outer edge of the disk (we set in the simulations), because here the Rayleigh–Jeans approximation applies to every part of the disk. Larger stars cause deeper transits (see below).

4.3.2 Transit maps

Figure 4.2: The instantaneous narrow-band transit depth as a function of projected position of the transiting object at three different frequencies optical (top panel), EUV (middle panel), and soft X-rays (bottom panel). A transit lightcurve corresponds to the values along the stellar trajectory in the image plane shown. A non-spinning black hole is at the origin, the observation angle is relative to the accretion disk, and transit depths are shown for a main sequence O-type star with . The top half of the disk image is more severely distorted by gravitational lensing since it is farther from the observer than the black hole. The left side of the disk rotates towards the observer, thus appearing more luminous due to beaming, which results in a deeper transit. The transit is deepest if the projected position of the star crosses the most luminous regions closest to the SMBH. The black hole silhouette and the curved spacetime geometry becomes visible in the X-ray transit map.
Figure 4.3: Same as Figure 4.2, but for a giant star with . Here, the star is so large that the transit is much deeper, and the image of the accretion disk is smoothed out (e.g. the left-right asymmetry of beaming is not apparent).
Figure 4.4: Soft X-ray transit depth map for an AGN with a Kerr-BH spin for , same as the bottom panel of Figure 4.2 but for a spinning SMBH. The details of the accretion disk are smoothed since the ISCO radius in this case is smaller than the stellar radius. The accretion disk image is further modified because of the Kerr geometry.

As stated in Section 4.1, AGN transit observations can be used to map distant AGNs along the transit chord with unprecedented resolution. We now elaborate on this point. Figures 4.2 and 4.3 show the transit depth as a function the projected position of the transiting object in front of the accretion disk around a non-spinning BH for an early type main sequence star () and a giant star (), respectively. The three panels in both figures show the transit depth maps for different observing frequency: optical band (top panel, , ), extreme ultraviolet (EUV, middle panel, Å, ), and soft X-ray (bottom panel, , , ). The spatial variations in the transit depth maps imply time-variations of the observed AGN brightness as a transiting object moves across the image. A transit lightcurve corresponds to the values along the projected stellar trajectory in Figures 4.2 and 4.3. Conversely, observations of the transit lightcurve can be used to reveal the geometry of the accretion disk along a line in this image.

The maximum possible resolution of such a reconstructed image is set by the size of the transiting object and the spatial variations of the disk brightness. Since the disk temperature increases inwards, the emission at higher frequencies is confined to the innermost regions, implying that transit measurements at higher frequencies can give a higher resolution image (see different panels in Figures 4.2 and 4.3). Transits of smaller objects also yield a higher resolution image. However, transits of smaller objects are less deep and hence more difficult to detect.

The left-right asymmetry in the figure is due to different Doppler shifts for regions of the disk moving towards or away from the observer. The spacetime geometry leaves an imprint on the image, the top half of the disk image is distorted by gravitational lensing close to the BH. Remarkably, the black hole silhouette (i.e., the lack of emission within the ISCO) becomes directly visible in the X-ray image of a transiting O star (see Figure 4.2 bottom panel).

The transit depth map is also sensitive to the spacetime geometry both directly through gravitational lensing and indirectly through the change in the ISCO radius. Figure 4.4 shows the soft X-ray transit map of a Novikov–Thorne accretion disk around a Kerr BH with spin (c.f. bottom panel of Figure 4.2). In this case, the ISCO is smaller than the radius of the transiting object and the BH silhouette does not appear in the image.

4.3.3 Transit lightcurves

To get a handle on the characteristic transit duration, let us consider the timescale for the center of a star to cross a disk of radius centered on the most luminous part of the accretion disk image for a given frequency. We choose the value of based on the transit depth map, depending on the desired transit depth. Typically .

The transit duration is


The transit duration is of the order of hours for massive main sequence stars closest to the AGN, and weeks for the closest giants further out in the circumnuclear star cluster (see Table 4.2 below).

Figure 4.5: Lightcurves of a transit: transit depth on logarithmic scale as a function of the horizontal position of the transiting star in front of the AGN in units of (bottom axis), or time in case (top axis). Stellar radius is (top and middle panels) and (bottom panel). The SMBH spin is (top and bottom panels) and 0.9 (middle panel). The BH silhouette is larger than the star in the top panel only.

Figure 4.5 shows multiband transit lightcurves of an accretion disk due to stars of radius (top and middle panel) and (bottom panel). The horizontal axis shows the horizontal position of the transiting star in the image plane of Figures 4.24.4. Physical distance in the projection plane is displayed in units on the bottom axis, and the top axis shows time in hours for a star with minimum orbital radius as given in Table 4.1. Note that further stars will cause a similar lightcurve, only scaled in time. The dimensionless black hole spin is zero on the top and bottom panels and in the middle panel. The projected stellar trajectory is horizontal with , which approximately passes through the hottest part of the accretion disk. The black hole is behind the origin. The curves show the transit depth for the same three frequencies as the transit depth maps in Figures 4.24.3. Note that these graphs are different from the usual planetary transit lightcurves showing flux, but plotting transit depth is more convenient when using a logarithmic scale. A larger value corresponds to a deeper transit, that is, a larger photometric signature. Figure 4.5 is consistent with the expectations on the frequency dependence described above: at higher frequencies, the AGN is more compact, implying shorter and deeper transits. The sharp features near the center of the transit on the top and middle panels are due to the inner edge of the disk at the ISCO, and the left-right asymmetry is mostly due to Doppler shift of the approaching and receding parts of the disk. However, a transiting giant filters out this spatial feature due to its size, as seen on the bottom panel. The black hole silhouette is clearly visible in the top panel where the ISCO radius is larger than the transiting object. Comparing the top two panels, we infer that if the mass of the SMBH is known, observing a transit light curve allows us to put an upper limit on the spin.

We conclude that for an accretion disk around a SMBH, photometric accuracy is required in the optical, accuracy in the extreme ultraviolet, and accuracy in soft X-ray to detect a transit due to a 11 star. In order to detect a transit due to a 110 giant, accuracy is sufficient in the optical, and accuracy in extreme ultraviolet to X-ray (see Section 4.5.1 for a comparison with intrinsic variability).

For comparison, we ran simulations for two different SMBH masses: and . Table 4.2 shows the maximum possible transit depth in each case for different sizes of transiting stars. The transit depth for a fixed stellar size decreases with the increasing spatial extent of the accretion disk around black holes with increasing mass, making the transit feature more prominent for smaller BH masses. However, AGNs with less massive SMBHs are also intrinsically fainter, and AGNs also have a smaller local number density which makes low mass AGN transit observations more challenging (see Section 4.5 for a discussion).

4.4 Transit rates

In this section, we estimate the expected rates of stellar transits in AGNs. Recall that denotes the radius of the circular area that the center of the star has to transit in the projection plane for a given transit depth. Let denote the probability of a single star on a circular orbit with radius transiting this circular area at a given instance, and let denote the probability that the orbit is aligned so that the star transits this circular area at some point during its orbit. By geometric arguments, these probabilities are


For a single AGN, the expected value of the number of transits at a given instance, , the expected value of the transit rate (number of transits observed per time for asymptotically long observations), , and the expected value of the number of stars on orbits such that they ever transit, , can be calculated by integrating for the entire stellar population:


where is the Keplerian orbital period.

To interpret these probability indicators, we have to understand the relationship between the expected value of the transit rate and the expected value of the number of stars that ever transit a given AGN. If we monitor a single target for which , then the actual transit rate is zero with probability and with a small probability . (The probability of multiple ever transiting stars in a single AGN is negligible in this case.) However, when monitoring a large number of identical targets, these effects average out: the observed total transit rate is and there are stars causing transits among all targets in total.

We substitute Equations (4.24.3) into Equations (4.154.17), and use and (justified by Table 4.1) to obtain


Here , , and denote the value of for the corresponding population density exponents. Stars on close orbits transit more frequently, and they dominate and . However, is dominated by stars on close or wide orbits for and , respectively. This is determined by the exponent of in the integrand of Equations (4.154.17).

model transit transit
(nm) depthaaThe maximum transit depth is determined by model parameters: SMBH mass and spin, circumnuclear stellar population, and observing wavelength. depthbbThis transit depth is chosen to be less than the maximum transit depth. Subsequent columns show values dependent on this parameter. The same transit depth is chosen for corresponding and cases to demonstrate how little probability estimates depend on spin. (yr) (h) (h)
0 O star 477 0.038 180 280 0.009 89 750
145 0.076 235 210 0.012 120 980
6.2 45 1 100 0.002 22 190
0.9 O star 477 0.049 200 250 0.010 99 830
145 0.12 240 210 0.012 120 1000
6.2 40