Analytic solutions to the maximum and average exoplanet transit depth for common stellar limb darkening laws

The depth of an exoplanetary transit in the light curve of a distant star is commonly approximated as the squared planet-to-star radius ratio, (R_p/R_s)^2. Stellar limb darkening, however, results in significantly deeper transits. Here we derive analytical solutions to the overshoot of the mid-transit depth caused by stellar limb darkening compared to the (R_p/R_s)^2 estimate for arbitrary transit impact parameters. In turn, this allows us to compute the true planet-to-star radius ratio from the transit depth for a given parameterization of a limb darkening law and for a known transit impact parameter. We calculate the maximum emerging specific stellar intensity covered by the planet in transit and derive analytic solutions for the transit depth overshoot. Solutions are presented for the linear, quadratic, square-root, logarithmic, and non-linear stellar limb darkening with arbitrary transit impact parameters. We also derive formulae to calculate the average intensity along the transit chord, which allows us to estimate the actual transit depth (and therefore R_p/R_s) from the mean in-transit flux. The transit depth overshoot of exoplanets compared to the (R_p/R_s)^2 estimate increases from about 15% for A main-sequence stars to roughly 20% for sun-like stars and some 30% for K and M stars. The error in our analytical solutions for R_p/R_s from the small planet approximation is orders of magnitude smaller than the uncertainties arising from typical noise in real light curves and from the uncertain limb darkening. Our equations can be used to predict with high accuracy the expected transit depth of extrasolar planets. The actual planet radius can be calculated from the measured transit depth or from the mean in-transit flux if the stellar limb darkening can be properly parameterized and if the transit impact parameter is known. Light curve fitting is not required.


Introduction
The planetary radius (R p ) is one of the key properties that are currently being derived by several exoplanet hunting surveys. The most successful method to determine the radius of an exoplanet is the transit method, which measures the slight decrease of the brightness of a star as the planet traverses the stellar disk as seen from Earth (Struve 1952).
In a simple picture, the star's appearance can be modeled as a circle with uniform brightness. Then the ratio of the planetary radius and stellar radius (R s ) can be estimated from the constant 1 transit depth (δ) via (R p /R s ) = √ δ. In reality, however, stars show center-to-limb brightness variations that affect the estimated planet-to-star radius ratio (Csizmadia et al. 2013).
The common variable used throughout the literature to describe limb darkening is µ = cos(γ), where γ is the angle between the emerging specific intensity of the star and the line of sight of the observer. Early attempts to model stellar limb darkening used equations that are linear in µ (Milne 1921) but 1 Even with the neglect of stellar limb darkening, however, the transit would still have an ingress and egress, both of which lead to a gradual transition between the in-and out-of-transit apparent stellar brightness (a mathematical description is given in Appendix A of Heller 2014). these became insufficient as computer models could be used to model stellar atmospheres. The quadratic (Manduca et al. 1977), square root (Diaz-Cordoves & Gimenez 1992), and logarithmic limb darkening laws (Klinglesmith & Sobieski 1970) provided better agreement with the observations. The advent of exoplanet transit observations (Charbonneau et al. 2000) and new highaccuracy space-based transit photometry (Brown et al. 2001), however, required even better precision. Claret (2000) presented a four parameter nonlinear limb darkening law that has now widely been adapted in the exoplanet community in combination with proper limb darkening parameters computed with stellar atmosphere models (Claret & Bloemen 2011). Nevertheless, even though the nonlinear law might be somewhat more precise in certain regimes of the parameter space and for certain stars (Espinoza & Jordán 2016), the quadratic law is commonly used because Mandel & Agol (2002) provided an analytic solution to the resulting transit light curve for the quadratic law.
Although it is well known that the transit depth is not equal to (R p /R s ) 2 , no framework exists that conclusively derives the actual deviations from that estimate for stars with limb darkening. Here we derive the correction factors to the (R p /R s ) ≈ √ δ approximation for each of the above-mentioned stellar limb darkening laws.

Methods
We present two methods to derive the actual planet-to-star radius ratio for arbitrary parameterizations of the common limb darkening laws. The first method is based on a measurement of the transit depth (Sect. 2.1) and the second method requires a measurement of the arithmetic mean of the in-transit flux (Sect. 2.2).

Transit depth overshoot
In Fig. 1 we show an example transit light curve using the python implementation 2 of the Mandel & Agol (2002) analytic model for quadratic stellar limb darkening by Ian Crossfield. We arbitrarily chose a = 0.4 and b = 0.4, irrespective of the stellar spectral type that this choice could imply, a transit impact parameter of zero, and a circular orbit. The solid line shows the numerical computation and the horizontal dashed line shows the (R p /R s ) 2 estimate. The vertical arrow at the bottom of the transit denotes what we refer to as the overshoot (o) between the actual transit depth and the (R p /R s ) 2 approximation. There are two ways to calculate o and we can show that they are equivalent to very high precision.

Transit depth overshoot in the light curve
The first way to calculate the overshoot uses the light curve. We consider a light curve with an out-of-transit flux normalized to 1, then the transit depth is simply the difference between 1 and the minimum in-transit flux, δ = 1 − f min . We then define the transit depth overshoot as determined from the light curve as the difference between the transit depth and the (R p /R s ) 2 estimate in units of (R p /R s ) 2 , This is the definition of the overshoot shown in Fig. 1 and it illustrates our problem. From an observational perspective, δ can be readily measured, but (R p /R s ) 2 is what we actually want to determine. Hence, if we knew the transit depth overshoot for a given host star of a transiting planet, then we could correct the simple (R p /R s ) ≈ √ δ approximation and derive the true planetto-star radius ratio as Fortunately, there is a second way to calculate the overshoot.

Overshoot of the stellar limb darkening profile
This second approach relates to the star's emerging specific intensity, I, which we simply refer to as intensity from now on. The only simplification that we require is to assume that the intensity profile covered by the planet in transit is constant, which is equivalent to the assumption that the planet is infinitesmal. As shown in Sect. 3.3, the resulting error from this small planet approximation is extremely small. We consider a star that is transited by a planet (Fig. 2). The minimum in-transit distance between the planet and the center of the stellar disk is referred to as the transit impact parameter 0 ≤ p ≤ 1. The radial coordinate, measured from the disk center in units of R s , is r and we orient the abscissa, or x coordinate, of our reference system to be parallel to the transit chord of the planet and to have its origin in the center of the stellar disk, so that r 2 = p 2 + x 2 . Introducing x 1 and x 2 as the distances traversed by the planet during the first and second halves of its transit, respectively, we find 1 = p 2 + x 2 1 = p 2 + x 2 2 . The key idea of this study is that the overshoot occurs at the minimum radial distance of the planet to the stellar disk center, that is, when x = 0 and r = p. Without loss of generality, we set the intensity at the disk center I 0 = 1 to simplify our notation. We only consider radially symmetric intensity profiles I(r) in this paper and so we refer to the intensity covered by the planet at mid-transit as I(r = p) = I p .
Irrespective of the actual limb darkening law, we can now define the expected transit depth overshoot as the difference between the intensity at mid-transit and the average intensity across the entire stellar disk area, I A , as Replacing in Eq.
(2) the overshoot in the light curve o LC , which is unknown a priori, with the overshoot of the intensity profile from Eq. (3) gives (4)

Mean in-transit flux
Alternatively, we can compute the planet-to-star radius ratio from the average in-transit flux f in and from the average intensity sampled by the transit chord I x . From the mere definition of these two quantities it is clear that Plugging this expression for δ into Eq. (4) we obtain which allows us to estimate the actual planet-to-star radius ratio based on the measured mean in-transit flux of a limb-darkened transit light curve and based on the disk average, radial average, and mid-transit intensity of the limb darkening profile.

Limb darkening laws
Now we set out to derive analytical expressions for the average intensity across the stellar disk I A and the average intensity along the planetary in-transit path parallel to the x coordinate I x .

Linear limb darkening
We start by rewriting the linear limb darkening law as in Claret (2000), originally derived by Milne (1921, Eq. (45) therein) but now as a function of the angle (γ) between the line of sight and the emerging intensity, where u is the limb darkening parameter. With r = sin(γ) as the distance from the star disk center normalized by the stellar radius (see Fig. 2) and using cos ( arcsin(r) ) = √ 1 − r 2 we have We then calculate the total intensity across the area elements (dϕ dr r) of the stellar disk, where 0 ≤ ϕ ≤ 2π is the azimuth angle, as The disk-averaged intensity then is where the factor π in Eq. (11) is the area of the apparent stellar disk with a normalized radius of r = 1. Using Eq. (11) as I A and I lin (r = p) from Eq. (8) as I p , we can estimate the planet-tostar radius ratio from the observed transit depth via Eq. (4) and based on the linear limb darkening law. Alternatively, we can calculate the transit depth overshoot from Eq. (3). Moving on to the mean intensity across the transit chord, we use the relationship x 1 = 1 − p 2 = −x 2 to parameterize the lower and upper boundaries of integration. Instead of integrating the intensity from ingress (at x 2 = − 1 − p 2 ) to egress (at x 1 = + 1 − p 2 ), we take advantage of the fact that the average intensity transited by the planet during the first half is the same as during the second half of the transit, which again is equal to the average intensity transited during the entire transit. Hence, which can be used as I x , for example, in Eq. (6).

Quadratic limb darkening
Rewriting the quadratic limb darkening law as a function of the angle and the normalized radius to the disk center gives where a and b are the two limb darkening parameters. The total disk-integrated intensity I tot,quad can then be calculated as per the right-hand side of Eq. (9) and the disk-averaged intensity then turns out as I tot,quad /π, giving A137, page 3 of 7 A&A 623, A137 (2019) An equivalent expression to Eq. (16) has been presented by Csizmadia et al. (2013, Appendix B therein) to calculate the transit depth as a function of µ. Hippke & Heller (2019) have used this formula for their Transit Least Squares transit detection algorithm (TLS). The transit chord average of the intensity in the quadratic limb darkening law can be obtained in analogy to the previous section using the right-hand side of Eq. (12). This results in

Square root limb darkening
The square root limb darkening law as a function γ or r appears as respectively, where c and d are the corresponding limb darkening parameters. We derive the total disk-integrated intensity I tot,sqr as per Eq. (9) and then divide I tot,sqr by the normalized disk area π to obtain the disk-averaged intensity as The transit chord average of the intensity in the square root limb darkening then follows from the right-hand side of Eq. (12),

Logarithmic limb darkening
The logarithmic limb darkening law can be written as a function of γ or r in the following ways, where e and f are the limb darkening parameters. The total diskintegrated intensity I tot,log then follows via Eq. (9) and the diskaveraged intensity is equal to I tot,log /π, resulting in We note that in contrast to all other limb darkening laws treated in this paper, Eq. (26) can be larger than 1, which means that the star can be brighter on average than in the disk center. In the transit light curve this would manifest itself by a mid-transit depth that is not the deepest point in the transit trough.
The average intensity along the planetary transit chord in the case of logarithmic limb darkening law is computed from the right-hand side of Eq. (12) as

Nonlinear limb darkening
The nonlinear limb darkening law as a function γ and r can be written as − a 3 1 − (cos(γ)) 3/2 − a 4 1 − (cos(γ)) 2 , with a 1 to a 4 as the limb darkening parameters. The total diskintegrated intensity for the nonlinear limb darkening law I tot,nl is then obtained from the right-hand side of Eq. (9), and dividing I tot,nl by the normalized apparent disk area π gives Finally, radial averaging of the specific intensity in the case of the nonlinear limb darkening law is obtained according to Eq. (12) as where Γ 9 4 = 1.13300 . . .. This derivation requires solving an integral that does not have a closed form solution, similar to the case of the square root limb darkening law. We nevertheless find the elementary antiderivative shown above by applying a substitution, see the analogous case for the square root limb darkening law in Appendix A.

Comparing transit depths between limb darkening laws
For the linear limb darkening law, we can derive the single limb darkening parameter u o , which produces the same transit depth overshoot as any of the other limb darkening laws. The linear limb darkening law and the quadratic, square, and logarithmic limb darkening laws produce equal overshoots (o LD = o LD ) if A137, page 4 of 7 Of course, the relations between the other limb darkening parameters can be derived analogously by equating the expressions for I quad A and I sqr A , I quad A and I log A , and I quad A and I log A , respectively, and solving for whichever limb darkening parameters are of interest.

Transit depth overshoot
We apply our framework to the numerically computed light curve in Fig. 1 as an example. First, we measure δ = 100 − 99.9875% = 0.0125% = 125 parts per million. As this light curve has been computed using the quadratic limb darkening law and using a = 0.4 and b = 0.4, we use the same parameterization of Eq. (16) and calculate I quad A = 0.8. Moreover, with p = 0 in Fig. 1, we have I p = 1 and therefore Eq. (4) gives the correct value of (R p /R s ) = 0.01 used for this simulation. The corresponding overshoot is calculated from Eq. (3) as o LD = 0.25. For an analogous application of this procedure to an observed light curve the limb darkening parameters and the transit impact parameter need to be known. The limb darkening parameters can be read from pre-computed tables if the stellar spectral type is known and if the transit impact parameter can be estimated. The latter can be done analytically based on measurements of the transit depth, orbital period, the time between first and fourth contact, and the time between second and third contact (Seager & Mallén-Ornelas 2003). Figure 3 illustrates the transit depth overshoot of the quadratic limb darkening law for any combination of a and b and using p = 0. We use the tables of Claret & Bloemen (2011) to obtain the limb darkening coefficients in the Kepler band for several main-sequence stars with metallicity [Fe/H] = 0, surface gravity log(g) = 4.5, and zero microturbulent velocity. The spectral types of these stars are indicated as A5, F5, sun, K5, and M0 in Fig. 3, respectively. As a result, we find that the quadratic limb darkening law predicts an overshoot of o LD = 0.218 times (R p /R s ) 2 for a small planet around a sun-like star. This exact solution can be derived by plugging a and b in Eq. (16) and the resulting I quad A into Eq. (3), where I p = 1 for a transit impact parameter of zero. Table 1 lists the overshoots for the other template stars as well, together with the effective temperatures and limb darkening coefficients.

Mean in-transit flux
For the numerically computed example light curve in Fig. 1, where the simulated planet-to-star ratio is 0.01 and the limb darkening parameters are a = 0.4 and b = 0.4 for I quad A and I quad x in Eqs. (16) and (17), respectively, we measure f in = 0.999891. Equation (6) then predicts (R p /R s ) = 0.009996, the source of the small error in the planetary radius of 4 × 10 −6 R s ≈ 2.8 km being in the ingress and egress, as explained in the next section.

Small planet approximation
The error in the planetary radius (ε p ) from the small planet approximation can be calculated by comparing the planet-to-star radius ratio (R p /R s ) obs derived from the measured transit depth in a simulated light curve as per Eq. (4) with the actual planetto-star radius ratio that went into the simulation, (R p /R s ) sim . The simulation of a Jupiter transit with p = 0 (hence I p = 1) across For this plot, a transit impact parameter p = 0 has been chosen. The positions of several main-sequence stars as per Claret & Bloemen (2011) are indicated with filled circles in the a-b plane (see Table 1).  (4) then predicts (R p /R s ) obs = 0.099336 and the resulting error in the planetary radius of a Jupiter-like planet around a sun-like star is ε p = 5 × 10 −5 R ≈ 35 km. For or an Earth-like planet, this same exercise yields an error of ε p = 4 × 10 −8 R ≈ 28 m.
The main source of error in (R p /R s ) obs as predicted from transit depth is in the variation of the intensity over the small area of the stellar disk that is temporarily occulted by the planet. The larger the planet, the larger the variation, and the larger the error.
For both of these test cases, we also calculate the mean in-transit flux f in to predict (R p /R s ) obs via Eq. (6). For the Jupiter-sized test planet, we measure f in = 0.9894685 and with I A = 0.8211 and I x = 0.887745 I obtain (R p /R s ) obs = 0.0986959. The resulting deviation from the injected planetto-star radius ratio is ε p = 6.9 × 10 −4 R ≈ 480 km. The same exercise for the Earth-like planet yields an error of ε δ = 4.2 × 10 −6 R ≈ 2.9 km.
The main source of error in R p /R s as predicted from the mean flux is in the ingress and egress of the planet, which is not taken into account in our calculations of the radially averaged emerging specific intensity of the star. The larger the planet, the longer the times of ingress and egress, and the larger the error in the resulting (R p /R s ) obs estimate from the chord-averaged intensity.
These tests show that the effect from the small planet approximation is orders of magnitude smaller than any uncertainties arising from the total noise budget in real light curves. Even high-accuracy space-based observations of the brightest and most photometrically quiet stars (Gilliland et al. 2011) have a noise floor of several times 10 parts per million, which translates into an error in the derived radius ratio of >10 −4 for R p /R s = 0.1 and an error >10 −3 for R p /R s = 0.01 (the relative error increases for smaller planets). Our small planet approximation produces errors in R p /R s that are orders of magnitude smaller than the error coming from the uncertain limb darkening coefficients (Müller et al. 2013); this is a crucial issue for any fitting of observed light curves with numerical methods.

Conclusions
In this paper we present analytical expressions for the overshoot of the transit depth for small planets (R p R s ) with arbitrary transit impact parameter for the linear, quadratic, squareroot, logarithmic, and nonlinear stellar limb darkening laws. Equation (3) can be used to calculate the overshoot of the stellar emerging intensity with respect to the disk-averaged intensity, which translates into an overshoot of the observed transit depth. Equation (4) gives the actual planet-to-star radius ratio based on the measured transit transit depth (δ) and based on the limb darkening law and parameterization. For sun-like stars, exoplanet transits can be ∼ 20 % deeper than the (R p /R s ) 2 estimate.
We also derive analytical expressions to calculate the maximum transit depth δ from the mean in-transit flux f in for any of the limb darkening laws mentioned above. This approach is based on calculating the average stellar intensity along the transit chord of the planet across the stellar disk I x , which relates to f in and δ as per Eq. (5). Combined with our results for the overshoot of the transit depth, we derive Eq. (6) which gives an estimate for the actual planet-to-star radius ratio based only on the measured in-transit mean flux and on the parameterization of the respective limb darkening law, which requires knowledge of the transit impact parameter.
These expressions can be used to calculate the expected maximum transit depth for a given planet-star system with any parameterization of the common stellar limb darkening laws. This approach is not meant to replace the numerical fitting of light curves, the latter of which can constrain limb darkening and transit impact parameter. Our approach is certainly less accurate than the use of elliptical integrals (Abubekerov & Gostev 2013) or than the full modeling of the stellar intensity as a sum of spherical harmonics (Luger et al. 2019), let alone the numerical simulation of transits with model stellar atmospheres (Neilson et al. 2018). But our approach requires no numerical modeling and it could be more intuitive to understand.