Limb darkening and exoplanets: testing stellar model atmospheres and identifying biases in transit parameters

Limb-darkening is fundamental in determining transit lightcurve shapes, and is typically modeled by a variety of laws that parametrize the intensity profile of the star that is being transited. Confronted with a transit lightcurve, some authors fix the parameters of these laws, the so-called limb-darkening coefficients (LDCs), while others prefer to let them float in the lightcurve fitting procedure. Which of these is the best strategy, however, is still unclear, as well as how and by how much each of these can bias the retrieved transit parameters. In this work we attempt to clarify those points by first re-calculating these LDCs, comparing them to measured values from Kepler transit lightcurves using an algorithm that takes into account uncertainties in both the geometry of the transit and the parameters of the stellar host. We show there are significant departures from predicted model values, suggesting that our understanding of limb-darkening still needs to improve. Then, we show through simulations that if one uses the quadratic limb-darkening law to parametrize limb-darkening, fixing and fitting the LDCs can lead to significant biases -up to $\sim 3\%$ and $\sim 1\%$ in $R_p/R_*$, respectively-, which are important for several confirmed and candidate exoplanets. We conclude that, in this case, the best approach is to let the LDCs be free in the fitting procedure. Strategies to avoid biases in data from present and future missions involving high precision measurements of transit parameters are described.


INTRODUCTION
It has been known since the first observations of our Sun's surface that its observed intensity decreases towards the limb. This effect, termed limb darkening, crucially affects the shape of the transit signature a planet imprints in the observed stellar flux when passing in front of its host star. In practice, limb-darkening is parametrized by laws which depend on µ = cos(θ), where θ is the angle between the line of sight and the normal to a given point of the stellar surface. Some of the most widely used limb darkening laws in exoplanet transit lightcurve fitting are given by cn(1 − µ n/2 ) (non-linear law).
⋆ E-mail:nespino@astro.puc.cl † E-mail: ajordan@astro.puc.cl The non-linear law has as a special case the square-root law proposed by Díaz-Cordovéz & Giménez (1992), which is obtained by setting c3 = c4 = 0, and the variant, threeparameter law, introduced by Sing et al. (2009), with c1 = 0. The laws listed above do not include all the parametrizations available for limb-darkening. Klinglesmith & Sobieski (1970), for example, introduced the logarithmic law for early type stars, which is given by while Claret (2003) introduced an exponential law given by but these laws are less used for transit fitting in practice, mainly due to their complex forms which are harder to deal with computationally and are not implemented in the most widely used transit modelling codes to date (Mandel &  Confronted with measurements of an exoplanetary transit, observers usually model the effects of limb darkening by either fitting the limb-darkening coefficients (LDCs) of a given law or fixing some (or all) of them using tabulated values (see e.g., Claret & Bloemen 2011;Sing 2010, for some recent tables using the Kepler bandpass). However, exactly which strategy is optimal, and in which situations, is still unclear. Previous works (e.g., Csizmadia et al. 2013;Müller et al. 2013) have discussed this issue focusing on the uncertainty introduced on the parameters retrieved from transit observations but, to our knowledge, no study has yet addressed the potential bias introduced by them. Such study is called for as such biases could be limiting (1) the instruments currently obtaining high precision photometry like Kepler, (2) exoplanet population studies which by definition are based on averaging out random errors but not systematic ones (e.g., Schlaufman 2015) and/or (3) techniques that require high-precision measurements like transmission spectroscopy. Studying the biases introduced by the treatment of limb darkening and whether or not they are significant is the main aim of this work.
The sources of potential biases introduced by assuming a given limb-darkening law can be separated in four. One issue is the differences in the tabulated values of the limb-darkening coefficients, which are large even when the same model stellar atmospheres are used. For example, Csizmadia et al. (2013) has shown that the differences between different approaches at tabulating limb-darkening coefficients for the quadratic law with the ATLAS9 (Kurucz 1979) stellar atmosphere models can be as high as 20%, which if incorporated in the modelling can lead to significant increases in the uncertainties of the retrieved transit parameters. The second issue is the fact that, as shown by Howarth (2011), limb-darkening coefficients obtained directly from the intensity profiles cannot be compared directly to coefficients observed in transit photometry due to the fact that the optimization procedures are different in each setting. This in turn implies that using LDCs obtained from an intensity profile of the star as inputs in transit photometry should lead to biases in the retrieved transit parameters. The third issue is related to model complexity: we usually model limb-darkening, which in models is best represented by the non-linear law, with low-order laws such as the quadratic law. Finally, the fourth issue, and perhaps the most complicated to tackle, is the fact that we still do not know with certainty if available models do a good job at reproducing real intensity profiles of stars.
The paper is organized as follows. In §2 we detail our methodology for obtaining LDCs from ATLAS9 (Kurucz 1979) and PHOENIX (Husser et al. 2013) stellar model atmospheres, and compare our results with the literature in order to try understand the discrepancy between previously published tables. In §3 we compare our model LDCs to Kepler estimates for a sample of stars, taking in consideration the effects mentioned by Howarth (2011), and compare the performance of the ATLAS and PHOENIX models at predicting the observed limb-darkening effect. In §4 we explore the biases introduced on different transit parameters by fixing or fitting the LDCs, §5 presents a discussion of our results and §6 the conclusions of our work.

FITTING LIMB-DARKENING MODELS
The fitting of the limb-darkening laws to intensity profiles obtained from stellar model atmospheres is, in principle, a relatively straightforward procedure. Given a normalized response function for a given telescope/detector S λ (λ) (e.g., S λ (λ)dλ = 1, although the normalization doesn't really matter in practice for the calculation of limb-darkening coefficients), one must integrate the specific intensity I λ (λ, µ) at each angle µi = cos(θi) multiplied by the response function, i.e., where λ1 and λ2 define the wavelength limits of the band. This gives the observed (by the instrument) intensity, which can then be fitted by any of the laws cited earlier after normalizing by I(1). Of course, for CCDs, the recorded quantity are photons and therefore we have to divide the integrand in equation (1) by hc/λ, where h is Planck's constant and c is the speed of light.
In this work we make use of two widely used model intensity libraries to calculate LDCs: the ATLAS9 model atmospheres, available from Robert L. Kurucz's webpage 1 and the 1D PHOENIX model atmospheres (Husser et al. 2013). These models differ both in the geometry used to solve the stellar atmosphere (with ATLAS models using a plane-parallel approximation and the PHOENIX models using a spherically symmetric atmosphere), and on the actual physics used on each of them (see, e.g., Plez 2011, for a short up-to-date discussion on this matter), so a difference between the LDCs computed from both models is expected. We will derive and compare LDCs from these stellar atmosphere models using the Kepler high-resolution response function 2 in order to compare our results with those in the literature.

Fitting limb-darkening laws with the ATLAS models
We first note that the ATLAS intensities are given per unit frequency and not per unit wavelength, so a a c/λ 2 term has to be added in Equation (1), and so the integral in this equation now reads where the limits of integration, λ1 = 348 nm and λ2 = 970 nm are, in our case, the first and last wavelengths present on the Kepler response function. We recall, however, that in order to fit the laws cited in the introduction what we want is the intensity normalized by I(1), i.e., I(µi) I(1) = λ 2 λ 1 Iν(λ, µi)S λ (λ)/λ dλ λ 2 λ 1 Iν(λ, 1)S λ (λ)/λ dλ .
(2) In the upper panel, solid lines correspond to fits obtained following the method of S10, while the dashed lines correspond to fits obtained by following the method of CB11 (note these overlap in the leftmost panels). The lower panel shows the (percentual) difference between these fits. The vertical dashed black line marks µ = 0.05 (see text).
After using numerical integration and interpolation to perform these integrals, the resulting normalized intensity profiles were fitted by a least-squares procedure using two different approaches in order to compare our results with previous works. The first approach was to fit all the angles for the non-linear law, while fitting only intensities at µ 0.05 for the rest of the laws, a procedure followed by Sing (2010), hereafter referred to as S10. The second approach was to fit to 100-µ points obtained by interpolating the 17 angles given by the ATLAS models 3 in a linear grid from µ = 0.01 to µ = 1.0 in 0.01 steps using a cubic spline, a procedure followed by Claret & Bloemen (2011), hereafter refered as CB11 4 .
We note that, because the coefficients are linear with respect to the parameters in all of the laws considered, the least-squares solutions are unique and therefore have analytic solutions given a set of intensities and angles. The exact solution, outlined in Appendix A, was used to perform the least-square fits. Figure 1 shows the intensity profiles for a G5V star of solar metallicity (i.e., with an effective temperature of T eff = 5500, log-gravity of log g = 4.5, [M/H] = 0.0 and microturbulent velocity of v turb = 2 km/s) where, for illustration, fits to the most popular limb-darkening laws used in the literature have been plotted following the methods of S10 (solid lines) and CB11 (dashed lines). The lower panel shows the (percentual) difference between those two methods, illustrating that they actually differ by a small (but, as we will show, significant) amount, with the differences being larger for low order laws such as the linear (with a median difference of ∼ 1%, reaching a ∼ 2.5% difference near the limb) and smaller for high-order laws like the nonlinear (with a median difference of ∼ 0.1%, and reaching a ∼ 1% difference near the extreme limb). As one would expect, the difference between the methods of S10 and CB11 is more important for the low order laws.

Comparing LDCs to previous results
The large panels in Figure 2 show the limb-darkening coefficients u1 and u2 for the quadratic law, for a range of effective temperatures, T eff , calculated by us using the AT-LAS models for stars with solar metallicity ([M/H] = 0.0), log g = 4.5 and v turb = 2 km/s using the S10 and CB11 procedures (blue and black points), along with the same calculations made by S10 5 (red points), and CB11 (green points). The small panels just below each large one depict the absolute difference between our values and the ones obtained by those previous studies.
Overall, it can be seen that the coefficients we calculate follow similar trends as the ones presented in previous studies, except for the region between T eff ≈ 7500 − 8500 (gray bands in Figure 2) where both our results and the ones of CB11 show a smooth decrease in the u1 coefficient and a smooth increase in the u2 coefficient, while the S10 results show a sharp increase and decrease in u1 and u2 respectively. We found that this is due to S10 using old versions This work, S10 method This work, CB11 method S10 Figure 2. Limb-darkening coefficients for the quadratic law obtained through the methods described in this work (blue and black), along with previous results by Sing (2010, red) and Claret & Bloemen (2011, green). Below each graph the absolute difference between our results and each previous study is shown, while the difference between those studies is depicted by gray points and between the S10 and CB11 methods (but with our procedures) in black points. Note the large differences in the range T eff ≈ 7500 − 8500 (region denoted by the gray bands) with the Sing (2010) results; this is due to Sing (2010) using old versions of the ATLAS stellar model atmospheres.
of the ATLAS stellar atmospheres 6 . Using those old models (available in Robert L. Kurucz's webpage) we recover the overall shape of the coefficients presented by S10 but cannot, however, eliminate the differences between S10 and the LDCs we calculate using the S10 methods. Despite the similarity of the overall trend, there are significant differences between the LDCs we calculate and previous studies, even though we use their methods as described in their works and the same stellar atmospheres. The differences are of the same order of magnitude (∼ 10%) for both coefficients of the quadratic law, with larger differences in the u2 coefficient for cooler stars (T eff 5000 K), and in the u1 coefficient for hotter stars (T eff 30000 K). Around solar temperatures (5000 K T eff 6000 K) the differences are smaller, ∼ 0.1 − 1% between our results and CB11, and ∼ 1 − 10% between our results and S10. Although we show here only the differences with the quadratic LDCs, there are very significant deviations between the coefficients of other limb-darkening laws too. In particular, we note that the limb-darkening coefficients of the non-linear law obtained by the different methods vary widely between works and, thus, these must be used with caution. We rule out that the differences arise from different interpolation and/or integration methods by performing the same calculations in different ways and in different programming environments.
Given the results above, it is clear that the differences with previous studies do not arise from different numerical approaches but, rather, they arise either from the actual method used to obtain the LDCs (i.e., from the very definition of the integrals used to estimate the integrated intensities) or from differences in the model atmospheres used to obtain those coefficients (e.g., older versions of the ATLAS model atmospheres). Unveiling the actual source of these discrepancies is currently not feasible as the actual codes used to obtain these tables are not publicly available. We believe that in the era of high precision measurements making all details available for scrutiny in order to try understand discrepancies is fundamental. Following this logic, we provide all the codes needed to reproduce the results in this paper through GitHub 7 . Our codes can be used to obtain LDCs for arbitrary response functions.

Comparison of the methods used to fit the limb-darkening laws
We now discuss which of the two fitting methods, i.e., that of S10 or that of CB11, is the most appropriate for obtaining LDCs using stellar model atmospheres, as they clearly show significant differences in Figure 2 (black points on the small panels).
To assess which of the two methods is better, we define below a criterion that relies on an accurate analytic description of the intensity profile given by the ATLAS models. To this end, we follow Howarth (2011) and use the non-linear 6 See the note on http://kurucz.harvard.edu/grids.html; the new models are the ones ending in *new.pck (e.g., im01k2new.pck). The old ones (e.g., im01k2.pck19) are also available on the webpage (Kurucz, private comm.). 7 http://www.github.com/nespinoza/limb-darkening Fit to 17 angles (S10 method) Fit to interpolated 100 angles (CB11 method) Figure 3. Differences between LDCs for the quadratic law obtained following the different methods discussed in this work (S10 in red, CB11 in green) and the "real" underlying quadratic LDCs.
limb-darkening law as an accurate descriptor of the intensity profile for the same stellar parameters used in Figure 2. We note that this was done not to try to exactly reproduce the model intensity profiles but, rather, to emulate an intensity profile that is similar to the one that would be actually observed in a real stellar atmosphere. From this profile, we then obtain quadratic LDCs by sampling N = 17 points at the same µ values as the ones given by the ATLAS models and, with this, fit the intensity profiles using the methods of CB11 and a method similar to that of S10 (for this experiment, we choose to fit all the µ values and not only the values for which µ 0.05, which is the original method of S10, in order to have a fair comparison between the coefficients obtained by these two methods).
In order to assess how good the retrieved coefficients are, we obtain the "real" underlying quadratic LDCs by defining them as the coefficients one would recover when sampling a very large number of points N from the (real) profile and then fitting those points with the quadratic limbdarkening law. Under this definition, when N → ∞, a set of "limiting coefficients" should be recovered, which are the "real", underlying quadratic LDCs that best-fit the underlying profile. We thus take the limit as N → ∞ in our least-squares procedure, imposing that the data are sampled from an underlying profile described by the non-linear law (see Appendix B). The resulting limiting coefficients for the quadratic law are given in terms of the coefficients of the underlying non-linear law by (3) Figure 3 shows the results of this experiment. The fitting method of CB11 is clearly better at obtaining the "real" quadratic LDCs as we have defined them. The reason is that a cubic interpolation of the profile does a good job at retrieving the underlying profile (the non-linear law in the case of our experiment) and, thus, sampling points from this interpolation is very similar to sampling more points from the profile, and thus a result closer to the "real" coefficients follows. It is also interesting to see that both methods fail at retrieving the "real" u1 coefficients for very hot stars. This is due to the fact that limb darkening for hotter stars is sharper than for cooler stars, changing abruptly towards the limb.
Because the 17 angles given by the ATLAS models are more densely sampled close to the limb, this gives more weight to that region of the profile and thus the fit is dominated by it. Sampling uniformly across the profile, which is what the method of CB11 does, alleviates this problem and gives an overall better fit to the whole profile. As a final note on this topic, we would like to note that the limiting coefficients we have introduced are the best by construction because to obtain them we use a sampling scheme that weights the whole profile uniformly. To obtain the limiting coefficients one has to make use of the non-linear LDCs, which are obtained by fitting the model intensity profiles using one of the discussed methods. However, as shown in the lower panels of Figure 1, in the case of the non-linear law the fitting method assumed is not very relevant, producing median offsets on the order of only 0.1%. Because of this, we believe that the "best" quadratic LDCs that one can choose are the "limiting coefficients", whose input non-linear limb-darkening coefficients might be obtained from either of the S10 and CB11 methods, or some other similar sampling scheme. For practical applications, we recommend using the CB11 method to fit the non-linear law to a given intensity profile and then use the resulting coefficients to obtain the limiting coefficients according to equation (3). The codes we provide are able to carry out these calculations.

Fitting limb-darkening laws to the PHOENIX models
After performing the analysis for the ATLAS stellar atmospheres, we now describe how we obtain LDCs with the PHOENIX models. The method we used is very similar to the one described for the ATLAS models. except for the fact that now we have 78 angles and, unlike ATLAS models, the intensities are given per unit wavelength. There is an additional, very important difference: unlike previous studies, we are careful in taking in consideration that the geometry of the PHOENIX models is different from that of the ATLAS models and, thus, the original intensity profiles obtained from these two model atmospheres are not directly comparable.
Spherically symmetric models like the PHOENIX ones extend the atmosphere by a small but important fraction (≃ 0.4% for older versions of PHOENIX, Aufdenberg, Ludwig & Kervella 2004) from what one usually calls the "stellar radius", which is where the Rosseland mean optical depth, τR, satisfies τR ∼ 1. Plane-parallel model atmospheres like ATLAS, on the other hand, by definition have a very thin extended atmosphere outside the stellar radius and, therefore, have τR ∼ 1 effectively at (or very close to) the limb (µ = 0). If we recall that the normalized stellar radius is given by r = 1 − µ 2 , this implies that for plane-parallel model atmospheres this outer "Rosseland radius" is effectively at r = 1 (i.e., at µ = 0), while for spherically symmetric models this radius is at r < 1 (i.e., at µ > 0). This fact is evident when plotting the profiles of plane-parallel and spherically symmetric model atmospheres, where the PHOENIX models have a sudden decrease in intensity around µ ∼ 0.05 that is not present in the ATLAS profiles. This, of course, is not a fundamental difference between the models but, rather, a difference arising from the fact that these two models are defining in a different way what the "outer stellar radius" or the "limb" is.
One way of accounting for the different treatment of the stellar limb in spherical models is to search for the point rmax at which τR ∼ 1 and re-define that point as having µ = 0 before fitting any law to the intensity profile. This procedure allows to have meaningful and comparable profiles between plane-parallel and spherically symmetric models. According to empirical results with the PHOENIX models, the point rmax in the intensity profile can be found by searching the value at which the derivative of the intensity profile with respect to the radial intensity profile, |dI/dr|, is maximum (see, e.g., Wittkowski, Aufdenberg & Kervella 2004). Once found, it suffices then to re-normalize the radial profile by dividing by rmax, thus re-defining then the point at which r = 1 (µ = 0). Figure 4 shows the results of applying such corrections to the PHOENIX profiles for a G5V type star, where we have plotted the original PHOENIX profile, the "shifted" version of it and the ATLAS intensity profile for comparison, which also illustrates the sudden decrease in intensity already mentioned for the PHOENIX models near µ = 0.05 (r ≈ 0.9987). The lower panel of Figure 4 shows, for illustration, the derivative of I(r)/I(0), where the maximum is indicated by the dashed line. It is evident from the upper panel of Figure 4 that the (shifted) PHOENIX and ATLAS profiles agree for this star at the limb (which was expected because now both models are modelling the same portions of the stellar disk). It is also important to see the significant changes between the original PHOENIX profile and the corrected ones: it is clear that the original PHOENIX models were not modelling the same portions of the stellar disk as the plane-parallel models and, thus, any LDCs derived from it couldn't be directly compared to the AT-LAS models. We note that this correction was not made by Claret, Hauschildt & Witte (2012, hereafter referred as CHW, when deriving their "quasi-spherical" LDCs which were obtained in order to be compared to ATLAS models, nor by Neilson & Lester (2013a,b), who also make direct comparisons between the spherical version of the ATLAS models and their plane-parallel counterparts. Figure 5 shows fits to the shifted PHOENIX intensity  profiles to the most popular limb-darkening laws used in the literature for the same model shown in Figure 4, following the methods of S10 (solid lines) and CB11 (dashed lines). Overall, it can be seen that in the case of the PHOENIX models a larger deviation is observed between the methods used to fit the profile than for the ATLAS models (compare the lower panels between this and Figure 1). In terms of following the profile, the CB11 method does a better job than the S10 method at following the whole profile in the case of the non-linear law. This was expected in light of our discussion in Section 2.1.2. As for the quadratic and linear law fits, both methods seem to do an adequate job at describing the profile given the low flexibility of the models. Figure 6 shows the analogue of Figure 2 for the quadratic LDCs obtained using the PHOENIX model atmospheres and the different methods discussed in this work, where a solar metallicity, log g = 4.5 and v turb = 2 km/s has been used. For comparison, the coefficients of CHW have been plotted in orange in the big panels 8 . We can see in the Figure that the overall trend for stars hotter than ∼ 4000 K is similar for all cases except in the region between ∼ 7500 − 9000 The panels show the same information as Figure 1, but for the "shifted" PHOENIX models.

Comparing LDCs with previous results
K, where the differences are on the order of ∼ 50% and for stars cooler than ∼ 4000 K, where the differences are larger than 100%. The overall differences in the coefficients are large in comparison to the ones seen for the ATLAS models, with median differences on the order of ∼ 20% for both coefficients with previous works. This is mainly due to our approach using a new parametrization of the stellar disk for obtaining limb-darkening coefficients from spherical model atmospheres such as PHOENIX.

MEASURING THE LIMB-DARKENING EFFECT FROM TRANSIT LIGHTCURVES
A quantitative determination of the limb-darkening effect from transit photometry requires careful thought and has been the subject of some scrutiny in recent years. For example, Csizmadia et al. (2013) has shown that the presence of unnoculted spots can lead to significant changes in the LDCs, due to the fact that a spotted stars' intensity distribution is more complex than the simple laws given in the introduction. Furthermore, even in the case of unspotted stars the interpretation of LDCs is subtle. Howarth (2011) showed that the geometry of the transit biases the LDCs obtained from photometry. In particular, he shows that high-impact parameter transits highly bias the observed LDCs because they sample chords of the stellar surface which are closer to the limb and, thus, sample a very different part of the intensity profile as the one sampled when fitting a whole model intensity profile. This is an unavoidable problem and, therefore, LDCs obtained from high-impact parameter transits are not directly comparable to the ones obtained by model intensities.
Howarth (2011) also stresses that, because the optimization processes are different when fitting an intensity profile directly from model stellar atmospheres than when fitting transit lightcurves, the coefficients obtained by those two procedures are not directly comparable even if the tran-sit is central. This is a very important and often overlooked fact: it implies that limb-darkening coefficients obtained from transit lightcurves such as the ones obtained by Müller et al. (2013) should not be compared directly to LDCs obtained from intensity profiles derived from model stellar atmospheres. That this is relevant can be verified in a straightforward fashion with a simple simulation study: (i) Select a good representation of a model intensity profile for a given set of stellar parameters derived from stellar model atmospheres, such as the non-linear limb-darkening law (i.e., select a set of coefficients c1, c2, c3, c4).
(ii) Generate a (noiseless) synthetic transit lightcurve by feeding the chosen representation of the model intensity profile of the star and using any set of geometric parameters for the transit.
(iii) Fit this synthetic transit lightcurve with the same code as the one used to generate it, but now using a quadratic limb-darkening law parametrization, fixing all the geometric parameters of the transit to its input values (i.e., letting just the limb-darkening coefficients to float).
The result of this experiment is always the same and, at first, counter-intuitive: the quadratic LDCs obtained from the input model intensity profiles (obtained through, e.g., the limiting coefficients obtained directly from the nonlinear coefficients using equation 3), denoted in what follows as (u1, u2), are always different from the ones obtained from the fit to the synthetic lightcurve, which we denote as (u * 1 , u * 2 ). Figure 7 shows the results for this simple experiment, where we used a simulated transit lightcurve of a planet with a period of P = 3 days, semi-major axis to stellar radius ratio a/R * = 0.1, planet-to-star radius ratio Rp/R * = 0.1 on a circular (e = 0, ω = 0), edge-on (inclination i = π/2) orbit (a typical "hot-Jupiter"). The non-linear law LDCs (c1, c2, c3 and c4) were obtained for stars with log g = 4.5, [M/H] = 0.0 and v turb = 2 km/s for different temperatures of interest for exoplanet studies using the ATLAS models via the method of CB11 and the Kepler re- This work, S10 method This work, CB11 method Figure 6. Limb-darkening coefficients for the quadratic law obtained through the methods described in this work (blue and black), along with previous results by Claret, Hauschildt & Wittle (2012 in orange (here we plot the coefficients obtained with the "quasi-spherical models", which only fit values of µ 0.1 using the original PHOENIX intensity profiles); below each graph the absolute difference between our results and this previous study is shown.  (2011), where the quadratic LDCs obtained by fitting the intensity profiles from model stellar atmospheres (u 1 , u 2 ) for stars with log g = 4.5, [M/H] = 0.0 and v turb = 2 km/s using the ATLAS models for different temperatures are plotted against these same coefficients (u * 1 , u * 2 ) obtained from synthetic transit light-curves generated using these same intensity profiles. The gray line follows the temperatures of the simulated systems in increasing order for better visualization. The dashed lines depict the u i = u * i lines, which the points should follow if the coefficients from model intensity profiles and transit lightcurves were directly comparable. sponse function, and from these we derived the limiting coefficients (u1, u2). The synthetic transit lightcurves were generated using the formalism of Mandel & Agol (2002), and the coefficients (u * 1 , u * 2 ) obtained via non-linear least-squares using the Levenberg-Marquardt algorithm. 1000 points between phases −0.05 and 0.05 were generated in each simulated lightcurve.
From the experiment we can see that, qualitatively, the effect of observing the LDCs from photometry is to overestimate the "underlying" u1 coefficient, while the effect for the u2 coefficient is to underestimate it. The exact mapping (u1, u2) → (u * 1 , u * 2 ), however, is non-trivial and, thus, the comparison between LDCs obtained from transit photometry, (u * 1 , u * 2 ), and from intensity profiles obtained from model stellar atmospheres, (u1, u2), is also non-trivial. Furthermore, this mapping is in theory also dependent on the transit parameters, because the optimization process used to fit a transit lightcurve is dependent on them, although we find this dependence to be usually negligible. There is thus no simple rule to compare LDCs obtained from transit photometry with the ones obtained by model intensity profiles, and the comparison has to be done on a case-by-case basis.

Comparing photometrically obtained limb-darkening coefficients to model values
In order to compare the LDCs obtained from photometric measurements to values obtained from model intensity profiles, Howarth (2011) introduced the SPAM (Syntheticphotometry/atmosphere-model) algorithm, which can be summarised in three steps: (i) Fit a transit lightcurve and obtain the best-fit geometric parameters θp = {a/R * , Rp/R * , e, ω, i} and best-fit quadratic LDCs Obtain the stellar parameters that best represent the stellar host (e.g., via spectro- (ii) Generate a synthetic transit lightcurve using the bestfit geometrical parameters θp, and using an accurate representation of the model intensity profile of a star with stellar parameters θ * , which we choose to be the non-linear law representation.
(iii) Fit the synthetic transit lightcurve to obtain estimates of the geometric parameters and LDCs, but now using a quadratic limb-darkening law to parametrize the limbdarkening effect. With this, obtain the coefficients u * = {u * 1 , u * 2 } which can now be directly compared to the best-fit coefficients obtained from photometry u f = {u f 1 , u f 2 }. There are three very interesting things to note about this algorithm. First of all, note that this algorithm actually tests the performance of the non-linear law, i.e., of one of the best representations of the model intensity profiles that we have available, because those coefficients are the ones that are used as inputs for the algorithm which are then compared to observations. This is very interesting, as it gives us a simple tool to asses how well our model atmospheres do at predicting real, observed, intensity profiles through transit lightcurves. The second thing to note is that, because the non-linear law is the law being used as input, the actual method for obtaining the LDCs (i.e., flux-conservation method, least-squares, etc., see Howarth 2011, for a discussion on the differences on the obtained coefficients by those methods for the case of the quadratic law) is irrelevant in practice, because all those methods give the same coefficients in the case of the non-linear law, as shown by Claret (2000). Finally, it is very important to note that this algorithm deals mainly with the bias introduced by modelling a complex intensity profile (e.g., one following a law close to the non-linear) with a simple, two-parameter law such as the quadratic one in the transit fitting procedure. This can be verified by noting that small changes in the input non-linear LDCs (e.g., obtained from stars with similar temperatures, metallicities or gravities) give rise to larger changes between the modelled (ui) and photometrically observed (u * i ) LDCs than changes in the transit parameters. This means that although it is true that this algorithm deals with the bias associated with transit geometry, differences in this mapping between systems are currently mainly dominated by stellar parameters rather than by the actual transit geometry of a given system.
Although a big step forward in how to properly compare limb-darkening coefficients, the SPAM algorithm has a few shortcomings: (1) the geometric parameters are fitted in step (iii), which is inconvenient as one would want to fix them, because the objective of the algorithm is to perform the mapping (u1, u2) → (u * 1 , u * 2 ) given a geometric setting; and (2) the algorithm does not account for uncertainties on the geometric parameters of the transit lightcurve, θp, and/or for uncertainties on the stellar parameters, θ * . The latter is a problem because, as stated above, changes in the geometry of the transit lead to changes in the limb-darkening coefficients obtained from the SPAM algorithm, while errors on the stellar parameters can lead to significantly different coefficients c1, c2, c3, c4 used in step (ii) to generate the synthetic lightcurves. Ideally one does not obtain only best-fit transit parameters for the geometric parameters, but also their posterior probability distribution functions (PDFs) given the data via, e.g., a Markov Chain Monte Carlo (MCMC) algorithm, which is information that we want to use. In most cases an MCMC approach for obtaining stellar parameters is not practical and therefore posterior distributions for these parameters are not usually available. One can still fit an adhoc distribution to a given set of stellar parameters and their respective errors (e.g., a Gaussian for the case of a parameter with symmetric error bars), and use that as an approximation to the posterior distribution of the stellar parameters. We propose here a modified version of the SPAM algorithm that takes into account uncertainty information. We term this algorithm Monte-Carlo SPAM algorithm (MC-SPAM), and it consists of the following three steps: (i) Fit a transit lightcurve and obtain the posterior PDFs given the data for the geometric parameters, i.e., p( θp|data), with the corresponding posterior for the LDCs, p( u f |data). Obtain the posterior distribution of the stellar parameters (e.g., via spectroscopic observations), p( θ * |data2) 9 .
(iii) Repeat step (ii) to obtain a Monte-Carlo sample of the model/photometric LDCs u * , which can now be directly compared to the observed LDCs u f .
Because in most cases no MCMC chains are available for the stellar parameters, in order to sample directly from p( θ * |data2), we assume independence among the stellar parameters and sample the parameters directly from their respective posterior marginal distributions, i.e., For the modelling of the posterior marginal distributions, we assume the quoted estimates and errors come from an analysis done on a χ 2 surface, which is essentially an analysis of the likelihood L(χ 2 ) and, thus, can be approximated by an analysis of the posterior distribution of the parameters, p( θ * |data2) ∼ L(χ 2 ), if one assumes flat priors for them. In practice, we model each marginal distribution by a Gaussian in the case of symmetric error bars, where the mean is set to the best-fit value of the parameter and its variance to the square of the error. For asymmetrical error bars, we assume a skew-normal distribution (Azzalini 1985), which is the natural choice for a "normal-like" distribution with lack of symmetry. Details on the method we use to fit the parameters of this distribution given an estimate of a parameter and a set of asymmetrical error bars, as well as how to sample from the resulting distribution are given in Appendix C. We note that although log g, v turb and T eff are positive quantities and our treatment can in principle allow negative values, in practice, the parameter uncertainties do not allow such values to be sampled and even if they did, our algorithm is capable of detecting and discarding those values as invalid.
Our MC-SPAM algorithm is available at GitHub 10 . We now use it to compare our estimates of the observed limbdarkening coefficients to a set of Kepler planets.

Comparing model to observed LDCs using
Kepler data Photometrically derived LDCs from fits to Kepler transit lightcurves were obtained from various sources in the literature in order to retrieve LDCs for stars with different parameters that host transiting planets with different geometries. Because the posterior distributions for the parameters of those systems are not published, we choose to use the same parametrization used for the stellar parameters in order to sample values given their published quantiles (i.e., in case of symmetrical errorbars we assume the posterior, marginalized distribution of each parameter is a Gaussian, while for the case of asymmetrical errobars we assume the posterior is best described by a skew-normal distribution). We note that the uncertainties associated with the transit parameters had negligible influence on the retrieved estimates of the limb-darkening coefficients obtained with our MC-SPAM algorithm, which were mainly dominated by the stellar parameter uncertainties. We took data from the high signal-to-noise, low-impact parameter sample of Müller et al. (2013) where we additionally removed the objects that had impact parameters larger than b = 0.5, which were Kepler-43b ( In order to be as conservative as possible, we also removed 10 http://www.github.com/nespinoza/mc-spam all the objects which have not been confirmed as planets to date, including Kepler-71b, for which no spectroscopic confirmation has been published so far that can rule out the possibility of it being a low-mass star, as Howell et al. (2010) can only constrain its mass to be less than 0.1 solar masses. We also remove Kepler-3b (HAT-P-11b) from our sample because its host star has been shown to have a significant amount of activity (Fraine et al. 2014) and thus the LDCs might be significantly biased ). In addition, we add Kepler-93b to our sample (Ballard et al. 2014), the recently validated Jupiter-sized planets KOI-206b and KOI-680b (Almenara et al. 2015) and, at the expense of larger errors on the estimated limb-darkening coefficients but in order to expand the effective temperatures sampled in this work, we also add the newly confirmed planets with low impact-parameters Kepler-186f, Kepler-296f, Kepler-296e, Kepler-436b, Kepler-439b, Kepler-440b, Kepler-441b, Kepler-442b and Kepler-443b (Torres et al. 2015). Table 1 summarizes the best-fit transit parameters for each of these planets, while Table 2 summarizes the host-star parameters. Table 3 presents both the model LDCs, (u1, u2), and the u * = (u * 1 , u * 2 ) coefficients, obtained using 1000 samples from our MC-SPAM algorithm for each target. We used both the ATLAS and PHOENIX models with the methods discussed in Section 2. A microturbulent velocity of 2 km/s was assumed for stars for which no measurement of this parameter was published. Figure 8 shows the published limb-darkening coefficients, (u f 1 , u f 2 ), as white points with errorbars. For easier visualization, the planets have been divided into two groups: the ones with low precision (upper panels) and high precision (lower panels) limb-darkening measurements; note that these also separate the cooler (upper panels, T eff = 3572 − 5431 K) and the hotter (T eff = 5520 − 7650 K) stars in the sample. At the sides of each datapoint, the MC-SPAM results, (u * 1 , u * 2 ) for each system using the ATLAS models (left, blue points) and PHOENIX models (right, red points) are shown. The arrows show the changes from the median of the model quadratic LDCs sampled by the MC-SPAM algorithm, (u1, u2) (i.e., the values obtained directly from the intensity profiles; in this case, they are the "limiting coefficients" obtained using the non-linear law), to the median of the resulting MC-SPAM values, (u * 1 , u * 2 ). Figure 9 shows the differences between the measured LDCs and the MC-SPAM values, with colored bands indicating the 68% band of the mean of those differences for the case of the ATLAS (blue bands) and PHOENIX (red bands) models. Kepler-296e was omitted from the calculation of the distribution of the mean of the differences (see below).
In general, just as we observed in our simulation shown in Figure 7, the shifts between the (u1, u2) and (u * 1 , u * 2 ) coefficients produced by the MC-SPAM are larger for the u2 coefficients than for the u1 coefficients, and the effect is to correct the model underestimation of the u1 coefficients and the model overestimation of u2; in other words, the effect of the mapping u1 → u * 1 is to increase the value of u1, while the effect of the mapping u2 → u * 2 goes in the opposite direction. As observed in Figures 8 and 9, this mapping is in general very effective at better describing the observations, especially if the ATLAS model intensity profiles are used as inputs. It is interesting to note that the errors given by MC-SPAM are always dominated by the errors on the stel-    Figure 8. High (bottom panels) and low (upper panels) precision quadratic LDCs derived from transit photometry for several exoplanets (white datapoints) and MC-SPAM model LDCs (u * 1 , u * 2 ) using ATLAS (blue, to the left of each datapoint) and PHOENIX (red, to the right of each datapoint) models. The blue and red arrows next to the MC-SPAM results represent the mapping u i → u * i , i.e., from the original model LDCs obtained from fits to the intensity profiles (in practice obtained using the non-linear coefficients through equation 3) and our MC-SPAM estimates. The temperature of the host star of each system is indicated above each of the planet names, inside the figures. Note the change in scale between the upper and lower panels.    Figure 9. Differences between the observed LDCs, (u f 1 , u f 2 ), shown in Figure 8 as white datapoints and the MC-SPAM estimates, (u * 1 , u * 2 ), using the ATLAS (blue) and PHOENIX (red) models, which are the blue and red points next to the white datapoints in Figure 8. The arrows represent the effect of the MC-SPAM algorithm in the observed difference between (u f 1 , u f 2 ) and the original model LDCs obtained from fits to the intensity profiles. The blue bands indicate the 68% bands around the mean of the differences using the ATLAS results (with the median of this indicated by the blue solid line), while the red bands show the same for the PHOENIX models (with the red solid line indicating the median of this distribution). These medians and the associated 68% values are also indicated next to each band (note that in the upper panels both bands overlap). The dashed black line marks zero, for reference. lar parameters; this fact is especially evident for Kepler-77, whose uncertainty in the microturbulent velocity of ±0.3 km/s leads to larger errorbars on the results using the AT-LAS models. This highlights the importance of estimating this parameter directly from spectroscopic measurements.
From the sample of low precision LDCs, which is also the cooler star sample (upper panels of Figures 8 and 9), one can see that, without taking into account Kepler-296e, the agreement with the data is very good for ATLAS and PHOENIX model atmospheres for the u1 coefficients, which shows a mean difference of −0.046 +0.133 −0.143 if one uses the AT-LAS models and −0.049 +0.138 −0.140 if one uses the PHOENIX models, both of which are consistent with a zero mean difference. There is a barely significant offset of the u2 coefficients (2.8σ for the ATLAS models, 2.9σ for the PHOENIX models), with the model values apparently systematically overestimating those coefficients by a mean value of ∼ 0.3. For the sample of high precision limb-darkening coefficients, which is also the hotter star sample (lower panels in those figures), the agreement of the u1 coefficients is very good for the ATLAS models, with a mean difference of −0.008 +0.012 −0.012 which is consistent with zero, but poor for the PHOENIX models, whose mean difference is 0.083 +0.008 −0.008 , inconsistent with zero by more than 10σ. Note that this bias is very evident and more prominent for stars hotter than 6000 K, where the model values overestimate the LDCs by ∼ 0.1. For the u2 coefficients, both model atmospheres show slightly significant biases, with the PHOENIX models doing a better job at predicting the observed LDCs with a mean difference of −0.024 +0.012 −0.011 , which is 2σ away from zero, and with the AT-LAS models showing a mean difference of 0.042 +0.012 −0.011 , which is 3.8σ away from zero.
As a final note, there is one particular object worth discussing in detail, Kepler-296e, which has LDCs which deviate more than 2σ from those of systems with host stars of similar stellar parameters, including Kepler-296f, which according to Torres et al. (2015) orbits the same host star in an orbit almost two times farther away from it. As we can see in Figures 8 and 9, the coefficient changes induced in Kepler-296e due to the geometry of the transit are not expected to be very different from those of Kepler-296f, so the geometry of the system cannot explain the differences on the observed LDCs. Activity could, in principle, produce significant biases on the LDCs through unnoculted spots ), but it seems unlikely that activity affected only one of the observed transits. Because Kepler-296 is known to be a tight binary (Lissauer et al. 2014), one might be tempted to think that Kepler-296e maybe did not orbit the same star as Kepler-296f as claimed by Torres et al. (2015), but its companion. However, both of the stars in the system have actually very similar spectral types and, thus, very similar LDCs, which implies that the observed LDCs are actually very different compared to both stars in the system. An analysis of other alternative hypotheses is out of the scope of this work, but we note that these are the kind of analyses that can be performed from measuring, comparing and interpreting LDCs from transit lightcurves using our MC-SPAM algorithm.

THE EFFECT OF USING FIXED LDCS IN TRANSIT FITTING
In the past section, we showed that, as first noted by Howarth (2011), LDCs extracted from fits to intensity profiles of stellar model atmospheres, (u1, u2), are not directly comparable to the LDCs obtained from transit photometry, (u f 1 , u f 2 ). This is due to the fact that the two optimization procedures are significantly different from one another and, thus, a geometry dependent mapping using synthetic lightcurves has to be carried out in order to obtain the coefficients (u * 1 , u * 2 ) that can then be compared to the observed LDCs. This implies that even if one could measure with excellent precision the intensity profile of a given star, obtain its LDCs with that profile, and then measure those coefficients from transit photometry, also with excellent precision, there is an expected bias between the two sets of coefficients. This in turn means that if one fixes the LDCs obtained from the intensity profile in the transit fitting procedure, then one is using potentially biased coefficients that can then lead to biased transit parameters 11 . A strategy to avoid this bias could be to let the LDCs as free parameters in the fit. However, we also expect a bias on the transit parameters in this case if, as it is usually done, the intensity profile is modelled with a quadratic law which we we have seen is not able to accurately describe the full intensity profiles.
In addition to the above mentioned problems, there is the issue related to our imperfect knowledge of the underlying, "true", intensity profile. As we saw in §3, this is currently an issue as our models are not able to reproduce the observed LDCs with sufficient accuracy. On top of this, according to our results in §2, there are differences even between our own modelling of those profiles both between different model atmospheres and between the different methods used to derive the LDCs from them.
In order to explore these sources of bias, in this section we perform simulations to study the possible biases introduced by our limb-darkening assumptions on the retrieved transit parameters, using transit lightcurves generated with the formalism of Mandel & Agol (2002).

A simulation study
In order explore the effect on the retrieved transit parameters of fixing or having the the LDCs as free parameters, we simulate transit lightcurves with unit period, circular orbits and an assumed intensity distribution for the host star of the transiting planet. The choice of units such that P = 1 is just for convenience of sampling directly in phase and has no consequences for what follows. The geometric parameters of the transit we can retrieve from our simulated light curves are the planet-to-star radius ratio, p = Rp/R * , the semi-major axis to stellar radius ratio, aR = a/R * , and the inclination of the orbit, i. The simulations were performed as follows. First, based on the data from all transiting planets discovered to date, we choose to generate synthetic transit lightcurves for planets with all the combinations of parameters {aR, p} with values aR = {3. 27, 3.92, 4.87, 6.45, 9.52, 18.18, 200} and p = {0.01, 0.06, 0.11, 0.16, 0.21}. In order to explore the effect of different impact parameters, b = cos(i)aR, we also varied b from 0 to 0.9 in steps of 0.1. This defines 350 different orbital configurations for our simulations. For each orbital configuration, we simulated 100 noiseless, uniformly sampled lightcurves with 1000 in-transit points and 400 out-of-transit points each, whose initial times were randomly perturbed. We first assume perfect knowledge of the underlying intensity profile by generating the transits using the non-linear law with the coefficients {c1, c2, c3, c4} obtained in Section 2 for models with log g = 4.5, solar metallicity, v turb = 2 km/s and effective temperatures between 3500 K and 9000 K. We use the ATLAS models, the fitting method of CB11 and the Kepler bandpass. Once we generate a simulated light curve, we retrieve its transit geometric parameters using constrained non-linear least squares with the Levenberg-Marquardt algorithm using the lmfit package 12 . We perform the fit in two ways: (1) fixing the limb-darkening coefficients to the ones given by the quadratic law for the given star (using the limiting coefficients defined in Section 2); and (2) leaving the LDCs as free parameters in the fit. In the latter case, we fit for the parameters q1 = (u1 + u2) 2 and q2 = u1/2(u1 + u2), where the parameters q1 and q2 are constrained to be in (0, 1), in order to obtain physically sound solutions in our non-linear least-square fits (Kipping 2013). For each combination of the geometrical parameters, the median of the fitted parameters obtained from the 100 generated lightcurves, along with the corresponding errors, are reported. Figure 10 shows the results of our simulations for the simplest case of a central transit, where the percentual bias induced on p and aR is shown as a function of the effective temperature of the host star used to extract the limb-darkening coefficients; the upper panel shows the bias induced when fixing the LDCs to model values and the lower panel shows the bias induced when letting those be free parameters in the optimization procedure (the retrieved inclinations in this case were not reported as they all arrived at the input values). The sizes and colors of the points represent the input values of p and aR, respectively, with smaller, bluer points representing the lower values of both p and aR (p = 0.01, aR = 3.27), and big, red points representing the highest input values of those parameters (p = 0.21, aR = 200). Note that the error bars are plotted, but are smaller than the smallest points in the plot. Note also that some points overlap.

The case of central transits
We infer from our simulations that the biases are small for central transits (maximum of ∼ 0.2% on p in the case of fitting for the coefficients, ∼ −0.05% when fixing them), although significant for several exoplanets: from a query done to the NASA Exoplanet Archive 13 , 326 Kepler Objects of Interest (KOIs) have quoted uncertainties on p lower than 0.2% and for which the effects of this bias are important. Of those systems 57 are confirmed exoplanets and 113 show quoted uncertainties lower than 0.05% (22 of which are confirmed exoplanets). It is interesting to note that the bias seems to be more important for deeper transits (bigger points in the Figure), and the effect is to retrieve deeper transits and larger distances to the host star when fitting for the coefficients, while the opposite effect is introduced when fixing the coefficients, estimating slightly shallower transits and smaller distances to the host star than the real ones. Furthermore, the bias is clearly larger when fitting for the coefficients than when fixing them, which suggests that, if the underlying limb-darkening model is accurate at a higher level than the error made by fitting the coefficients, then the best strategy is fixing the LDCs to their model values.
We note that the shape that the biases have as a function of effective temperature are very similar between the different strategies employed to fit the transit lightcurves. For the Kepler bandpass, the difference between the ui (LDCs obtained from model intensity profiles) and the u * i (same coefficients recovered from transit fitting) observed in our experiment in the introduction of §3 (Figure 7) follow a similar shape with temperature as that observed for the biases, that is, starting at the cooler temperatures the offset is larger, decreasing until it reaches a minimum offset at T eff = 4750 K, then gradually increasing again for hotter host stars, slightly decreasing around T eff = 8500 K. This means that the actual mapping ui → u * i is less severe for temperatures around T eff = 4750 K, which is one of the reasons why we observe a smaller bias in the retrieved transit parameters around those temperatures in our experiments. Additionally, the performance of the quadratic law fit is also optimal around T eff = 4750 for this particular bandpass, becoming slightly worse for cooler and hotter stars, a fact already discussed in §2.1.2 (Figure 3). The two points above serve to explain the observed shape of the biases with temperature and in particular the observed minimum of the biases at T eff = 4750 K. Figures 11 and 12 show the results of our simulations for the cases of low (b = 0.3) and high (b = 0.8) impact parameters. We can see that the value of b strongly affects the observed bias in the retrieved transit parameters, showing biases as large as ∼ 1% for p, and ∼ 2% for aR and i. These biases are larger than for the case of central transits and, therefore, more important. A query to the NASA Exoplanet Archive returns 933 KOIs with uncertainties better than 1% on p, 164 of which are confirmed exoplanets. Furthermore, the impact parameter not only modifies the order of magnitude of the observed bias, but also modifies the trends observed in the case of central transits. For example, in this case, the effect is more important for shallow transits than for deep transits at most temperatures.

The case of low and high impact parameter transits
One very interesting fact about our simulations is that although for low impact parameters the bias seems to be larger for aR and i when fitting for the LDCs, for high impact parameters this effect is reversed, showing higher biases in those parameters when fixing the coefficients. The same effect can be seen for the bias in p in the case of small planets around low temperature stellar hosts. This is in agreement with the results of Howarth (2011), who showed that  Figure 10. Biases in the recovered planet-to-star radius ratio p = Rp/R * and the semi-major axis-to-stellar radius ratio a R = a/R * as a function of temperature obtained from the simulations described in the text. The upper panels show results when fixing the LDCs, while the lower panels show the results when letting them to float in the transit light curve fit. The size of the points denote the input value of p (p = 0.01, small, p = 0.21 big points), while the color of the points represent the input value of a R (a R = 3.27, blue points, a R = 200 red points).
the difference between LDCs obtained from model atmospheres and the ones obtained from transit photometry increases as one increases the impact parameter of the transit. This implies that the fixed limb-darkening coefficient strategy should worsen as one increases the impact parameter, which is what we observe in our simulations. Therefore, for high impact parameter transits one should fit for the LDCs if one is interested in decreasing the bias, which contrasts with the suggestion of Müller et al. (2013) of fixing them for high impact parameter transits based on the observed increased uncertainty in the retrieved transit parameters when fitting for the coefficients. In our view, the best strategy strongly depends on the quality of the photometry and the observed geometry: if the retrieved uncertainties when fixing the coefficients are smaller than the biases shown here, then the best strategy should be to let the coefficients be free parameters.

The effect of an unknown stellar intensity profile
In order to explore deviations from the known intensity profile, we make use of previously published limb-darkening co-efficients in order to simulate systematic offsets from the limb-darkening profile being modelled. To this end, we used the same transit lightcurves generated in the past subsections, generated with our non-linear limb-darkening coefficients, but fitted them fixing the limb-darkening coefficients to the values published by CB11. In this way, we can explore the impact that different methods for obtaining the LDCs can have on the retrieved transit parameters which, we note, are only lower limits on the actual offsets between the modelled and real profiles, as both using our limb-darkening coefficients or the ones published by CB11 in our analysis in §3.2 lead to offsets between observed and modelled LDCs of the same order. Figures 13 and 14 show the results of our experiments.
As can be seen from our results, the biases follow a similar shape to the difference of the observed quadratic LDCs between our work and that of CB11 in Figure 2; that is, as expected, the biases are larger where the calculated LDCs show the largest differences. Also, the effect is again more pronounced for smaller values of p, with maximum biases on the order of 3% for planets around stars with T eff = 4000 K with moderate impact parameters (b = 0.3), temperature  Figure 11, but with an impact parameter of 0.8. Note that the scale is different than that of Figure 10, but the same as that of Figure 11. and geometry at which we also observe biases on aR on the order of ∼ 10%, as well as biases on the order of 5% on i. A query to the NASA Exoplanet Archive returns 2221 KOIs with uncertainties better than 3% on p, 491 of which are confirmed exoplanets. This is a clear illustration of the importance of our limb-darkening assumptions.

Lessons learned from the comparison of model-to-observed LDCs
In Section §3, we compared observed LDCs to model LDCs using our new MC-SPAM algorithm, which not only takes into account the geometry of the transit, but also the uncertainties in the transit and host star parameters. From these results, it is was apparent that overall the MC-SPAM results using the ATLAS models do a better job at predicting Host star T eff (K) Figure 13. Same as Figure 10, but with an impact parameter of 0.3 and using LDCs obtained from a different intensity profile (obtained from the work of CB11) to the underlying one (generated using non-linear LDCs from this work). Note that the scale is different from that of previous figures. Host star T eff (K) Figure 14. Same as Figure 13, but with an impact parameter of 0.8. Note that the scale between that figure and this one is also different.
the observed u1 LDCs for hotter stars than the PHOENIX models, with the tables slightly turned in the case of the u2 coefficients. Both models seem to do a good job at predicting the u1 coefficients of cool stars, while slightly overestimating their u2 coefficients. These results suggest that the non-linear law in general seems unable to predict the limb-darkening effect with the accuracy necessary given the quality of available data and, because this is the best representation that we have for the ATLAS and PHOENIX model atmospheres, this suggests that those models are still not reproducing the observed limb-darkening effect with sufficient accuracy. As such, these model results have to be taken with care and one must be careful in trusting too much the LDCs derived directly from model atmospheres. Our view is that one should always be aware of the possible biases introduced by fixing those coefficients in practical applications.

Limb-darkening and biases in transit parameters: implications for exoplanetary science
The biases uncovered by our experiments in §4 are, in general, non-negligible. In particular, the large biases shown for the cases of non-central transits both when assuming a perfectly known intensity profile ( §4.3) and when assuming deviations from it ( §4.4) have several implications for exoplanetary science. First, the biases shown for p have a direct impact on population studies that rely on accurate measurements of the planetary radius and their associated uncertainties which, as recently shown by Schlaufman (2015), impact the derived conclusions of these studies. It also has an impact on present and future studies that aim to constrain the interior composition of exoplanets from precise measurements of the planetary radius and mass. For example, Dorn et al. (2015) recently showed that a precision better than ∼ 2 % on the planetary radius is one of the key ingredients in order to be able to constrain the interior structure of rocky exoplanets. In addition, the biases for aR have an impact on derived quantities such as the calculated incident fluxes on exoplanets. This parameter, combined with the inclination i of the orbit, defines in turn the impact parameter, which according to our results can vary significantly depending on the different limb-darkening assumptions.

Can the biases be corrected and/or avoided?
As shown in in §4, the biases introduced by simple parametrizations of the limb-darkening effect such as the quadratic law are significant in comparison to the published parameters uncertainties for several confirmed and candidate Kepler exoplanets. A natural question to ask is: how can one correct for those biases and how can one avoid them?
The correction of biases in the transit parameters retrieved from a given fit using the quadratic law would necessarily have to involve simulations such as the ones shown here around the parameters of interest. However, an even simpler method would be to fit the original lightcurve with higher order limb-darkening laws like the three-parameter or non-linear laws. The problem with this approach is that, although it seems obvious to switch to higher order laws once the order of magnitude of the bias surpasses the uncertainties in a given parameter, there is always an intrinsic bias on the retrieved transit parameters related to our potentially imperfect knowledge of the exact shape of the, underlying intensity profiles in stars which, as shown in §4.4, can lead to large biases if the coefficients are fixed in the transit fitting procedure. As shown in Section 3, it appears that better modelling is indeed necessary given that the observed LDCs are not well reproduced by the models. This will be an important endeavour to minimize biases in the parameters estimated from transit light curves, and will be important to fully exploit data from both present high-precision missions like Kepler, and future high-precision missions like the Transiting Exoplanet Survey Satellite (TESS) or the James Webb Space Telescope (JWST).
Although a better modelling of stellar atmospheres is an important task, it is also a long-term one. Given the large amount of data already obtained by missions such as Kepler and the advent of future missions that will focus on retrieving high-precision transit parameters, it is important to discuss short-term strategies that can help overcome the biases due to our limb-darkening assumptions shown in this work. The most natural approach would be to fit the transit lightcurves with high order laws letting their parameters to float. However, this will most likely be a bad idea for low and medium signal to noise lightcurves due to the degeneracy between the parameters fitted in this procedure. One strategy to overcome this and the problems stated above could be to select a set of high-precision transit lightcurves with precisely known stellar host parameters that have enough signal-to-noise as to be able to be fitted with high order laws such as the three-parameter or non-linear law. This would in turn allow one to obtain fitted LDCs which could then be used to predict or constrain those coefficients for other stars. Although in theory those fitted LDCs are dependent on the transit geometry, this dependance as discussed in Section 3 is rather weak in comparison with the typical uncertainties of the host star parameters and, thus, this strategy seems to be a promising one for achieving high accuracy and precision measurements of transit parameters.

CONCLUSIONS
In this paper, we calculated LDCs using the Kepler bandpass for both ATLAS and PHOENIX stellar atmosphere models, and showed that we cannot reproduce previously published values. We could not fully resolve the difference due to the fact the procedures used in the literature are not openly available, but we believe the differences between the different sources of limb darkening coefficients available in the literature are due to the use of different input model atmospheres. We make our codes available for users to reproduce our results and calculate limb darkening coefficients for any bandpass in a flexible way.
We showed that different methods used in the fitting stage of the intensity profiles lead to different LDCs. We define a set of quadratic law limiting coefficients that are the best description of the "real" underlying intensity distribution in the sense that they are obtained by sampling uniformly the whole profile. Among previously published methods using the ATLAS models, the best is that of Claret & Bloemen (2011) which by generating additional points to the ones directly available from the models via spline interpolation provides a good approximation to a denser sampling of the underlying profile. We also point out an important correction that needs to be applied when dealing with the spherically symmetric PHOENIX models so that the stellar radius definition is consistent with the planeparallel ATLAS models. We provide updated limb darkening coefficients using PHOENIX models with this correction applied which are now directly comparable to those derived using ATLAS models.
In order to map model LDCs to those determined by transit photometry we introduce an algorithm called MC-SPAM (Monte-Carlo Synthetic-Photometry/Atmosphere-Model), which builds upon the SPAM algorithm proposed by Howarth (2011). The algorithm takes into consideration not only the fact that LDCs obtained from intensity profiles of model stellar atmospheres are not directly comparable to the coefficients estimated from transit photometry, but also the fact that all the stellar parameters and geometrical parameters of the transit have measurement errors. We use MC-SPAM to compare limb darkening coefficients predicted by models for a sample of systems which have had their limb darkening coefficients determined from Kepler transit photometry. We show that the ATLAS and PHOENIX models are not able to fully reproduce the observed limb-darkening effect for systems with a wide range of stellar parameters. Finally, we showed that, when using the quadratic law, fixing and letting the LDCs be free in the transit fitting procedures induce biases on the retrieved transit parameters that are significant for several candidate and confirmed Kepler planets, even if one assumes a perfectly known stellar intensity profile.
Given our results, we conclude that if one is confronted with a transit lightcurve and one decides to use the quadratic limb-darkening law to model the limb-darkening effect, the best strategy in order to minimize the bias introduced in the transit parameters is to let them float in the transit fitting procedure. Although a natural strategy to follow once the biases introduced by a given limb-darkening law are important for a given system geometry and precision would be to switch to higher order laws, one has to be careful if the strategy is to fix the limb-darkening coefficients. This is because the strategy would also lead to an underlying bias due to our still incomplete knowledge of the intensity distribution of real stars; our understanding of stellar model atmospheres is not good enough to avoid biases by switching to using fixed model coefficients from higher order laws. As a short-term solution to this problem, we propose to fit a set of high signal-to-noise lightcurves with precisely known stellar host parameters using a high order law parametrization such as the three-parameter or non-linear laws in order to obtain the LDCs directly from data, which could then be used to predict or at least constrain the values of those coefficients for other exoplanetary systems.
The results shown in this work imply that the current achievable precision is smaller than the current achievable accuracy if one considers only our limb-darkening assumptions. Because of this, we call the attention to observers to be careful about their limb-darkening assumptions, and suggest always leaving room for flexibility on the LDCs in the transit fitting procedures even when fitting high order laws in order to shield against biases in the retrieved transit parameters. If precision is more important than accuracy (e.g., in applications such as transmission spectroscopy), we suggest performing simulations similar to the ones shown in this work in order to quantify the order of magnitude of any biases if the choice is to leave the limb-darkening coefficients fixed.

ACKNOWLEDGMENTS
N.E. is supported by CONICYT-PCHA/Doctorado Nacional. N.E. & A.J. acknowledge support from the Ministry for the Economy, Development, and Tourisms Programa Iniciativa Científica Milenio through grant IC 120009, awarded to the Millennium Institute of Astrophysics (MAS). A.J. acknowledges support from FONDECYT project 1130857 and from BASAL CATA PFB-06. We would like to thank the anonymous referee for helpful comments, questions and suggestions that helped to improve this work. We would also like to thank Pierre Kervella & Antoine Mérand who pointed out to us the fact that PHOENIX models have a limb definition inconsistent with that of plane-parallel ATLAS models, Benjamin Rackham for his assistance on writing initial versions of the fitting procedures for the PHOENIX models, Antonio Claret, David Sing & Tom Evans for useful discussions regarding the procedures they used for obtaining LDCs, and Ashley Villar, Jontahan Fraine, Amaury Triaud and Rafael Brahm for useful discussions regarding the results of this paper.
Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate.

APPENDIX A: LEAST-SQUARES FITS TO LIMB-DARKENING LAWS
Fits of stellar model atmosphere intensities with the laws mentioned in the introduction require a least-squares procedure to be followed which, in general, minimizes the quantity where yi are the datapoints to be fitted, f (xi) the model for those datapoints, wi the weight given to each squared residual (yi − f (xi)) 2 and N the number of datapoints to be fitted. In our case we set wi = 1 and the objective is to minimize eq. (A1) using the various models presented in the introduction of this work. If no constraint is provided, then the problem is easily solved for all the laws presented because the coefficients are linear in the parameters. Our objective then, is to minimize eq. (A1) with yi = I(µi)/I(1), and xi = µi (the angles at which each of those integrations are done), for the different models f (µi) presented in the introduction, obtaining the optimal coefficients in a least-squares sense. The models have the form where θn are the parameters of the laws (e.g., θn = cn for the non-linear law) and gn(µi) are the functions which make up each law, e.g., in the case of the non-linear law (which contains the linear, square-root and three-parameter laws depending on which coefficient θn one sets to zero), gn(µi) = (1 − µ n/2 i ). In order to minimize eq. (A1), we calculate the partial derivatives of χ 2 with respect to the different coefficients θn and set them to zero. The calculation is easily found to give 2(f (µi) − I(µi)/I(1)) ∂f (µi) ∂θ k = 0, with ∂f (µi) ∂θ k = −g k (µi), which, after rearranging terms, gives the system of k linear equations for the n coefficients θn n θnα n,k = β k , k = 1, 2, ..., n, with, gn(µi)g k (µi), g k (µi)(1 − I(µi)/I(1)), which are trivial to solve. Note also that if we write the linear system as A θ = b, with A n,k = α n,k , θ = {θ1, θ2, ..., θn} T and b = {β1, β2, ..., βn} T , in this case the matrix A is symmetric, so the system is not only linear but very fast to compute. As an example of how to use the above result, for the linear law the only function is g1(µi) = 1 − µi so, in this case, there is only one equation for the parameter θ1 = a with parameters (1 − µi)(1 − I(µi)/I(1)), which gives, (A2)

APPENDIX B: LIMITING CASES FOR KNOWN TARGET LIMB-DARKENING LAWS
If we try to fit an arbitrary law when knowing that we are sampling from a law of the form I(µi)/I(1) (e.g., the non-linear law), then "limiting coefficients" corresponding to N → ∞ can be obtained. To perform this calculation, we sample N uniform µi points by defining µi = (i−1)/(N −1), with i = 1, 2, ..., N (note that this samples µi angles from µi = 1 for i = N to µi = 0 for i = 1). We introduce this sampling sums that give the parameters (e.g., eq. (A2) for the parameter of the linear law) and then take the limit as N → ∞. with known coefficients cn. We now fit the profile with a linear law and determine the limiting coefficient a that follows as N → ∞. We first find the numerator and the denominator in equation (A2). The denominator is easily found to be and we need to take the limit as N → ∞ The main challenge for taking this limit is to obtain a closed form expression for the limit of the ratio between the sums of non-integer powers of µi and h(N ), terms which appear in the ratios S1/h(N ) and S3/h(N ). To evaluate this, we first obtain an expression for the sums of non-integer powers of µi. In order to do so, we note that, for (integer and non-integer) exponent k = ±1, where in the last step we have used the Euler-Maclaurin summation formula. This implies that whose limit as N → ∞ is easily found as all terms to the right of the first term in this expression vanish in that limit, i.e., B2 Limiting coefficients u1 and u2 for the quadratic law when sampling from the non-linear law Following the results in Appendix A, in the general case of all the two-parameter limb-darkening laws there are two equations for the parameters θ1 and θ2, with θ1 = β2α2,1 − β1α2,2 α1,2α2,1 − α1,1α2,2 , (B1) θ2 = β1α1,2 − β2α1,1 α1,2α2,1 − α1,1α2,2 , Specializing to the case of the quadratic law (i.e. θ1 = u1 and θ2 = u2), we obtain (1 − µi) n+k , (1 − µi) k (1 − I(µi)/I(1)), for n = 1, 2 and k = 1, 2. We note that in this case α n,k = α k,n , with (1 − µi) 3 = N 2 4(N − 1) , α2,2 = with the Sn given in the past sub-section and, The denominator of the expressions for u1 and u2, i.e., the denominator of equations (B1) and (B2), is given by α = α 2 1,2 − α1,1α2,2 = (3N 2 − 3N + 2)(N + 1)(2 − N )N 2 720(N − 1) 4 .
We now note that the limits we want to obtain can be written as Given an estimate of a parameter θ in the formθ σ 1 −σ 2 , where in general σ1 = σ2, we want to sample points from the posterior distribution of θ, given the data, only knowing that the distribution is asymmetric. One choice for performing such sampling is to assume that the distribution of θ is a skewnormal distribution with parameters µ, σ and α, which is given by where, pα(y) = 2φ(y)Φ(αy), and, φ(y) = exp(−y 2 /2)/ √ 2π, Φ(αy) = αy −∞ φ(t)dt.
In this distribution, µ is the mean, σ 2 its variance and α, also called the shape parameter, defines its skewness. Note that when α = 0, we recover a normal distribution of mean µ and variance σ 2 .

C1 Fitting a skew-normal distribution to observed estimates
A simple way to obtain the parameters {µ, σ, α} of the distribution for each parameter is to assume thatθ − σ2,θ andθ + σ1 define the 0.16, 0.5 and 0.84 quantiles of the parameter distribution (i.e., the quoted value of the parameter defines the median and the errors define the 68% credibility bands around it). We can now easily formulate the problem as a non-linear least-squares problem where the independent, observed, variables are the observed values at the given quantiles (i.e., x = {θ − σ2,θ,θ + σ1}) and the dependent variables are the quantiles (i.e., y = {0.16, 0.5, 0.84}). Then, we minimize r = || y − m( x, µ, σ, α)||, where the i-th element of m is given by mi = x i −∞ p(θ|µ, σ, α)dθ, which can be solved with any non-linear least-squares algorithm such as Levenberg-Marquardt.

C2 Sampling from a skew-normal with known parameters
Once the parameters for a given skew-normal distribution are known, we want to sample values from it in a simple and efficient fashion. The following algorithm has been published by A. Azzalini in his personal webpage 14 , but we quote it here for completeness. First, one has to compute the parameter δ = α/ √ 1 + α 2 . With this, one now samples the random variables u0, v, both of which are independent and have standard normal distributions, and generate the random variable u1 = δu0 + √ 1 − δ 2 v which has correlation δ with u0. Then, sample the random variable z which equals u1 if u0 > 0 and −u1 otherwise. With this, z has a skew-normal distribution with zero mean, σ = 1 and shape parameter α. Finally, the random variable SN = µ + σz has a skew-normal distribution with mean µ, variance σ 2 and shape parameter α. This paper has been typeset from a T E X/ L A T E X file prepared by the author.