## Modelling the Transit Light Curve

The Sun (Credit: NASA)

We looked at a simplistic understanding of planetary transits and transit light curves (here). In the interests of accuracy, it would be worth investigating a more formal understanding of the transit light curve. Specifically, this post will be aimed at actually modelling transit light curve while taking into account some more complicated dynamics such as limb darkening.

The actual shape of the transit light curve can be mathematically described by a function of the sky-projected centre-to-centre distance between the star and the planet, z, where z = 0 is understood to be the mid-transit of a centrally-crossing transit and the impact parameter, b is understood to be the minimum value of z throughout the transit, and the transit starts at about z ~ 1. More accurately, first and fourth contact start at $z = R_* + R_p$, second and third contacts start at $z = R_* - R_p$. Furthermore, let’s define p as the ratio of the radii of the planet and the star, $p = R_p/R_*$.

Let us now define a function of the total fractional area of the star that is not obscured by the planet. Normally you might be inclined to feel that $(R_p/R_*)^2$ suffices, however it describes only the area ratio between second and third contacts, and is thus inadequate to describe the entire transit. We therefore will construct the amount of obscured stellar area as a piecewise function of p and z whose intervals describe different parts of the transit. It will be given as the ratio between the total area and the unobscured area subtracted away from the whole total area. Specifically, $A(p,z) = 1 - \lambda(p,z)$.

$\displaystyle \lambda(p,z) = \begin{cases} 0, & 1 + p < z\\ \frac{1}{\pi} \left[ p^2\kappa_0+\kappa_1 - \sqrt{\frac{4z^2-(1+z^2-p^2)^2}{4}} \right], & \lvert 1 - p\lvert < z \le 1 + p\\ p^2, & z \le 1 - p\\ 1, & z \le p - 1 \end{cases}$

Where $\kappa_1 = cos^{-1}[(1 - p^2 + z^2)/2z]$ and $\kappa_0 = cos^{-1}[(p^2 + z^2 - 1)/2pz]$. There’s a lot going on here but it isn’t too complicated to understand. The first case, $1 + p < z$, describes when the planet is completely off the stellar disc, before first contact hence there is no (zero) obscured area. The second case is the most complicated and describes the time between first and second (third and fourth) contacts. The third case is what we covered earlier, where the amount of obscured area is simply the square of the ratio of their radii, and therefore simply $p^2$. The final case is only satisfied in the (unlikely) condition that the planet-to-star radius ratio > 1, and the entire area of the stellar disc is obscured.

Calculating z as a function of time may be done with knowledge of the impact parameter and use of Pythagorean theorem.

For a transiting planet with an impact parameter of zero and a radius ratio of $R_p/R_* = 0.1$ (producing a 1% transit depth), the above model produces the following light curve plot.

Simulated Light Curve

Real light curves are actually curved, however, whereas this model is clearly much more “rigid” in shape. The reason for this is an effect known as limb darkening, where the outer layers of the star are not able to isometrically scatter light from underneath, causing less light from the limb of the star to reach the observer than from the centre of the stellar disc. For this reason, stars generally appear dimmer at their limbs and brighter in their centres. The magnitude of this difference will depend on a number of factors intrinsic to each star.

Proper handling of limb darkening will be required to obtain high-fidelity fits to high-accuracy data. A particularly popular and effective model for limb darkening was presented by Claret (2000). We can consider a four-parameter limb-darkening law describing the relative brightness of a point on the star expressed as a function of the angle between the observer, the centre of the star, and a line between the stellar centre and the given point on the star.

$\displaystyle I(r) = 1 - \sum^4_{n=1} c_n(1 - \mu^{n/2})$

Where $\mu$ is given as $\cos \theta = \sqrt{1 - r^2}$, and each value of c is called a “limb darkening coefficient.” In the case of the centre of the stellar disc, $I(0) = 1$. Each star will have a unique set of limb darkening coefficients, and they will change depending on what wavelength the transit is being observed in (due to how the properties of the stellar photosphere’s light scattering varies with wavelength). One could choose to have as many limb darkening parameters (and thus coefficients) as desired, but a high number is rarely required to adequately fit transit light curve data.

Limb darkening must be taken into account to get accurate measurements of radius ratios of transiting objects, and thus planetary radii. Because limb darkening will concentrate the brightness of the star into the centre of the disk, a planet that would cause a 1% transit depth in a uniformly bright star would cause a deeper transit depth in a central transit around limb-darkened star, and a shallower transit depth in a grazing transit around a limb-darkened star. Proper care in modelling the limb-darkening behaviour of the star is therefore required to ensure good measurement of the planetary radius.

We see therefore that the brightness difference is not a pure function of the ratio of the stellar and planetary radii. $\Delta F = (R_p/R_*)^2$ is only an approximation, and only correct for the unrealistic case of a uniformly bright stellar disc.

So what can we do to get a more comprehensive representation of a light curve? Let’s modify the equation for $\Delta F$ we’re familiar with and consider brightnesses instead of just purely area ratios. Let’s take our limb-darkening equation as the radius of the star and determine the area, since it describes the brightness of the star from centre to limb. With this equation as the stellar “radius,” we use some simple calculus to expand it to the “brightness area” of the circle that is the stellar disc. This calculates the total brightness of the star while taking into consideration the changing brightness behaviour from centre to limb. The actual transit depth of a planet at any given distance to the centre of the star (represented as z again) is simply the ratio of the total stellar brightness in the area of the star under the planet to the total brightness of the stellar disc. Therefore, we draft a new function based off our original understanding of $\Delta F$, while taking into effect the above model for how much of the star’s area is covered at any given time in the transit, as well as bringing in our limb-darkening model to describe the brightness distribution of the stellar disc:

$\displaystyle F(p,z) = \left[ \int^1_0 dr 2rI(r) \right]^{-1} \int^1_0 dr I(r)\frac{d[A(p/r,z/r)r^2]}{dr}$

And it makes quite a difference, too. By incorporating limb darkening coefficients, we’re able to much more accurately simulate the transit of a planet across a limb-darkened star. The diagram below shows the same planet (with p = 0.1) transiting stars of different limb-darkening parameters. The shape and modification to the transit depth is readily apparent.

Effects of Limb-Darkening coefficients

Detecting planets around other stars is hard. One of the easiest ways is to detect the apparent dimming of a star as a planet crosses between the observer and the star. The planet will block some of the photons from reaching the telescope, resulting in an apparent dimming of the star. This event is known as a “transit,” and the planets that do this are called “transiting planets.” The amount of light that is blocked is easy to calculate when one considers the problem from the perspective of simple 2-dimensional geometry. The amount of light blocked will be directly proportional to the amount of surface area of the star covered by the planet. For a star of radius $R_{star}$, and a planet of radius $R_{pl}$, the amount of the change in flux, $\Delta F$, is simply ratio of the area of the two bodies:

$\displaystyle \Delta F = \left( \frac{R_*}{R_pl} \right) ^2$

Clearly, then, larger planets will produce larger ΔF than smaller planets, which directly translates to a more detectable planet. But the difficulty does not scale linearly with the radius of the planet, as one can tell from the exponent in the above equation. Consider two planets, b and c, with radii $R_b$ and $4_c$, such that $R_b$ = 2$R_c$, essentially planet b is twice the radius of planet c. With the above equation it is seen that $\Delta F_b = (2^2)\Delta F_c = 4\Delta F_c$ – One-half the planet produces one-fourth the dip in brightness. Detecting Earth-sized planets can be understood to be far more difficult than detecting Jupiter-sized planets. The typical Jupiter-sized planet will be about 10 times the radius of the typical Earth-sized planet. This translates to 100 times the $\Delta F$. Can you detect a transiting Jupiter-sized planet? You need a hundred-fold increase in sensitivity to detect Earths.

Not all planets will transit. Assuming a random planetary inclination, the geometric probability that a planet will transit its star may be expressed (somewhat over-simplistically*) as

$\displaystyle P_{transit} = \frac{a}{2R_*}$

Where a is the semi-major axis of the planet’s orbit. We can see, therefore, that planets which orbit closer to their stars are more likely to transit than those further away. For Mercury, this works out to a ~1.19% transit probability. For Earth, this is an even more disappointing ~0.5%. This also means that if we assume that other stars are randomly distributed across the sky (which is not completely unreasonable out to distances where the structure of the galaxy does not become apparent), then we can say that ~0.5% of stars will have the right perspective to view Earth as a transiting planet. Similarly, we might also say that if all solar-radius stars have Earth-radius planets, then 0.5% of them are detectable through this method.

An important feature in determining how long transits last, aside from their orbital period, is the “impact parameter,” b, which is a measure of how far the transit chord is from the centre of the stellar disc, as measured in units of the stellar radius. A transit with b = 0 will be a perfectly edge-on orbit with the transiting planet passing straight through the centre of the stellar disk, with higher values of b being less dead-on transits.

Values of b > 1 imply a non-transiting planet, with the maximum attainable value for b being achieved for a face-on orbit (i = 0° or 180°) at $a/R_*$. There is a very small range of values of b > 1 for which transits still occur, depending on the radius of the planet. While the centre of the planetary disc may not intrude upon the stellar disc for these values of b, planets are not points, and so have a radius of their own. Transits with values of b close to 1 are called “grazing” transits, because the planet just grazes the stellar disk.

Mathematically, the impact parameter may be calculated by $b \equiv a \cos{i}$.

These transits may be plotted as brightness as a function of time, leading to a “transit light curve.” A transit event will have four events called “contacts.” The “first contact” is when the planetary disc first reaches the stellar disc. The “second contact” occurs when the entire planetary disc has moved onto the stellar disc. “Third contact” occurs when the planetary disc has reached the other edge of the stellar disc on its way out of the transit. Finally, “fourth contact” occurs when the planetary disc has moved completely off the stellar disc. The time from first to second contact is characterised by a significant drop flux, followed by a comparatively constant flux from second to third contacts. From third contact to fourth contact, the flux jumps up to the pre-transit value. The time between first and fourth contacts is the total transit duration, tt, while the time between second and third contacts is the “full transit duration,” tf, denoting the amount of time the planet is fully transiting the star.

For the image below, two transits are shown, one with a high impact parameter and one with an impact parameter of b = 0. Notice how differing the value of b changes both the duration and shape of the transit light curve. For both, all four contacts are labelled as vertical lines.

How does this look for the case of a real planet transiting a real star? Below is the transit light curve of the planet HD 189733 b.

The transit of HD 189733 b

Notice that the light curve between second and third contact is curved. This is because of stellar limb darkening, where the light coming from the limb of the star is darkened compared to the light from the centre of the stellar disc.

* Taking the radius of the planet, the eccentricity of its orbit, and the longitude of periastron into account, the geometric probability that a planet with a randomly oriented orbit will transit is expressed as

$\displaystyle P_{transit} = \left( \frac{a}{2(R_*+R_{pl})} \right) \left( \frac{1-e \cos{\pi/2 - \omega}}{1-e^2} \right)$

This dependence on the longitude of perihelion can be understood from the consideration of eccentric orbits. In reality, what is the dominant driver as to the probability that a planet will transit is more of the planet’s distance from the star during the transit window, as opposed to the planet’s distance from the star as measured by the semi-major axis. Below shows two planets with identical orbits, except for the latter having a higher value of ω. Note that because of this, the latter planet’s distance from the star during the transit window, shown in light-grey, is much further away than the first planet, and so its transit probability is considerably lower.