# LDTk: Limb Darkening Toolkit

###### Abstract

We present a Python package LDTk that automates the calculation of custom stellar limb darkening (LD) profiles and model-specific limb darkening coefficients (LDC) using the library of phoenix-generated specific intensity spectra by Husser et al. (2013). The aim of the package is to facilitate analyses requiring custom generated limb darkening profiles, such as the studies of exoplanet transits–especially transmission spectroscopy, where the transit modelling is carried out for custom narrow passbands–eclipsing binaries (EBs), interferometry, and microlensing events. First, LDTk can be used to compute custom limb darkening profiles with uncertainties propagated from the uncertainties in the stellar parameter estimates. Second, LDTk can be used to estimate the limb-darkening-model specific coefficients with uncertainties for the most common limb-darkening models. Third, LDTk can be directly integrated into the log posterior computation of any pre-existing modelling code with minimal modifications. The last approach can be used to constrain the LD model parameter space directly by the LD profile, allowing for the marginalization over the LD parameter space without the need to approximate the constraint from the LD profile using a prior.

###### keywords:

Methods: numerical–Planets and satellites: general–Gravitational lensing: micro–Binaries: eclipsing–Techniques: interferometric## 1 Introduction

Stellar limb darkening (LD) is a major source of uncertainty in the characterization of transiting exoplanets (Csizmadia et al., 2013; Espinoza & Jordan, 2015), interferometry-based stellar radius estimation (White et al., 2013; Neilson & Lester, 2013a), and analysis of microlensing events (An et al., 2002).

Limb darkening is approximated by a variety of limb darkening models, or laws (Mandel &
Agol, 2002; Giménez, 2006).
These models aim to reproduce the stellar intensity profile as a function of ^{1}^{1}1Where , is the normalized distance from the centre of the stellar disk, and is the foreshortening
angle. with a relatively small number of parameters, limb darkening coefficients.
In the context of exoplanet transit modelling, the early transiting exoplanet studies generally fixed the limb darkening
coefficients to values derived from the model and passband-specific fits to numerical stellar models by
Claret (2000, 2004). However, fixing the coefficients will introduce biases in the parameter estimates if the
models fail to reproduce the true stellar intensity profiles (Csizmadia et al., 2013; Espinoza &
Jordan, 2015), which has been shown to
be the case in several situations where inferences about the true limb darkening profile have been possible (Fields et al., 2003; Claret, 2008, 2009; Howarth, 2011, but
see also Müller et al. 2013). In addition, even with an accurate model,
fixing the coefficients would lead to underestimated uncertainties since the uncertainties in the stellar parameter
estimates are not propagated into uncertainties in the limb darkening profiles.

The problems with fixed limb darkening have led to the use of more robust approaches to transit modelling, and currently
the limb darkening is usually either constrained using (more or less informative) priors
based on the tabulated limb darkening coefficients, or using uninformative priors and marginalizing over all the limb
darkening profiles allowed by the observations (Kipping, 2013). Also, the modern limb darkening tabulations
by Claret
et al. (2013, 2014), Sing (2010), Neilson &
Lester (2013a, b), Magic
et al. (2014), and
Husser et al. (2013) are based on more advanced stellar atmosphere models than the early ones, using a spherical model
geometry instead of a plane-parallel one (Claret
et al., 2012, 2014; Neilson &
Lester, 2013a, b; Husser et al., 2013), using fully
three-dimensional models geometries (Hayek
et al., 2012), and including hydrodynamical simulations to account for the
stellar granulation (Magic
et al., 2014) .^{2}^{2}2Note that even with reliable stellar models, errors can arise from
the limb darkening model fitting. See Sect. 2.2 in Espinoza &
Jordan (2015), especially the need to renormalise the radius
from the spherical models.

Allowing the data to speak for itself–using noninformative priors and marginalizing over the limb darkening allowed by the data–yields the most reliable parameter estimates, but is not the optimal approach in some situations. For example, in transmission spectroscopy, where minute changes in the planetary radius as a function of wavelength are used to make inferences about the planet’s atmosphere, the uncertainties from unconstrained limb darkening easily drown any variation in the radius estimates. Further, the use of pretabulated coefficients is impossible due to the use of narrow non-standard passbands, and unaccounted-for small-scale variations in the limb darkening can lead to systematic trends and spurious features in the inferred transmission spectrum.

A more robust way to account for limb darkening in any modelling situation is thus to generate limb darkening profiles from the modelled stellar spectra on a case-by-case basis for each passband in a way that propagates the uncertainties in the stellar properties to uncertainties in the limb darkening profile, and allows us to set a factor telling how much we trust on the stellar atmosphere models used to create the spectra. This approach offers a compromise between tightly constrained and completely unconstrained limb darkening. The limb darkening profiles can be used to either create prior distributions for the limb darkening model coefficients by estimating the coefficient posteriors; the profiles can be used to directly constrain the limb darkening in the modelling process, fitting a limb darkening model to the profile simultaneously with the rest of the analysis; or the profiles can be used directly as is, if the modelling approach allows for arbitrary stellar limb darkening profiles.

We present a Python package LDTk that automates the calculation of custom stellar limb darkening profiles and model-specific limb darkening coefficients from the library of phoenix-generated specific intensity spectra by Husser et al. (2013). LDTk propagates the uncertainties in the stellar parameter estimates into the uncertainties in the limb darkening profile, and offers methods to create multivariate normal priors for the most common limb darkening laws. The package can also be directly integrated into the modelling code, bypassing the need to pre-calculate the priors. While the following discussion considers mainly exoplanet transit modelling, the package can be equally used in the context of eclipsing binary, microlensing, and interferometry studies, or in any work benefiting from custom-built stellar intensity profiles.

## 2 Examples

### 2.1 Calculation of model coefficients

We start with an example of limb darkening model coefficient estimation. LDTk offers methods to calculate the maximum likelihood estimates–as well as the full posterior distributions–for the coefficients of the limb darkening models listed in Table 1.

Limb darkening profiles and quadratic model coefficients for three simple passbands for a star with K, , and can be calculated as

We first import the `LDPSetCreator`

class that will create a set of limb darkening profiles for a given set of
filters. Next, we define the filters, here using simple boxcar filters that have zero transmission outside the given
minimum and maximum wavelength range. We continue by creating an instance of `LDPSetCreator`

, initialising it with
the filter set and stellar parameter estimates. `LDPSetCreator`

then creates a list of needed spectrum files,
checks them against the local cache, and downloads the uncached files from the Husser et al. (2013) FTP library. After
this, `sc.create_profiles`

is used to create the limb darkening profiles for each filter, shown in
Fig. 1, contained in an instance of `LDPSet`

class. The `LDPSet`

class is finally used to
estimate the quadratic limb darkening model coefficients `qc`

and their uncertainty estimates `qe`

,
illustrated in Fig. 2.
Figure 3 shows the quadratic model coefficients for the same star, but calculated for 19 narrow
(15 nm wide) passbands from 500 to 800 nm.

### 2.2 Integration into transit modelling

Rather than calculating the limb darkening coefficients prior to modelling, LDTk can be integrated directly into the code carrying out the log posterior evaluation. This can be done in two ways: (1) LDTk can be used to evaluate the log likelihood for the limb darkening model given model coefficients during the posterior sampling; (2) LDTk can be used to create a multivariate normal prior automatically in the initialisation before the sampling.

#### 2.2.1 Log likelihood evaluation

The log likelihood for an LD profile can be calculated as

where `cf`

is a coefficient array, and `xx`

is the model abbreviation from Table 1. The log
likelihood call can be integrated into the log posterior computation directly

where `LPFunction`

is a callable class that encapsulates the log posterior computation. The limb darkening
profiles are created during the first initialisation, and saved to a file so that we do not need to regenerate them
every time. After the initialisation, all that is required is to add a call to calculate the log likelihood to the
calculation of the log posterior. The example implements this in the `__call__`

method, where `pv`

is the
model parameter vector, and ldsl is a slice selecting the limb darkening coefficients from the parameter
vector.

#### 2.2.2 Multivariate normal prior generation

The second approach for the direct integration of LDTk is to use it to construct a multivariate normal prior in the initialisation

where `lnp_ld`

is now a function that returns the joint log density of multivariate normal priors estimated from
the profiles.

## 3 Implementation

### 3.1 Overview

The package aims to be simple to install and integrate into pre-existing Python-based modelling code, and depends only on standard Python modules, NumPy, SciPy, and PyFITS. The code with IPython-notebook-based examples is freely available from

The repository contains also examples of the package’s basic usage, ways to integrate the package into a transit light curve analysis, and notebooks detailing the implementation of the code.

### 3.2 Profile creation

The limb darkening profiles are created using the phoenix-calculated specific intensity spectra by Husser et al. (2013). The spectra are stored in an FTP server as FITS files–a single file for a single set of stellar parameters–where each file contains a model spectrum spanning the wavelength range from 50 to 2600 nm with a 0.1 nm resolution for 78 values of . Each specific intensity file is 15 MB, and the whole library has a size of hundreds of GBs. However, downloading the whole library is unnecessary, and only a small subset of spectra are required for any given star.

The profile creation starts with the calculation of the necessary spectrum files, continues with a check testing which
of the required files are already in a local cache directory, after which the missing files are downloaded from the
library FTP server. The `LDPSetCreator`

class carries out these tasks during its initialisation, including all
stellar spectra inside a -wide cube on the three stellar parameters, where are the stellar
parameter uncertainties, and is a freely set factor defaulting to three.

Next, limb darkening profiles are calculated for each filter and stellar parameter set as

(1) |

where is the detector quantum efficiency, is the filter transmission, and is the stellar spectrum. The profiles are renormalized in (see Espinoza & Jordan, 2015, Sect. 2.2), and then used to construct a () interpolation grid for each filter.

The final limb darkening profiles are calculated using Monte Carlo sampling. The `create_profiles`

method generates
stellar parameter samples from a multivariate normal distribution with the means and variances defined by the
stellar parameter estimates and their uncertainties (we assume symmetrical normal uncertainty distributions), and
calculates the limb darkening profiles for each sample. The final profile for each filter is defined by the mean of the
profile set, and the profile uncertainty by its standard deviation.

### 3.3 Limb darkening coefficient estimation

Name | Abbr. | Model |
---|---|---|

Linear | ln | |

Quadratic | qd | |

Nonlinear | nl | |

General | ge |

The `LDPSet`

class offers methods of form `coeffs_xx`

, where `xx`

is a limb darkening model abbreviation
from Table 1, to estimate the limb darkening model coefficients with their uncertainties. For example,
the quadratic model coefficients can be obtained as

where `qc`

contains the coefficient estimates and `qe`

their uncertainties for each filter. By default, the
estimates correspond to the maximum likelihood estimates, and the uncertainties to the 68% central likelihood intervals
assuming multivariate normal likelihood with zero correlation between the coefficients.
That is, we estimate the likelihood uncertainty as

(2) |

where is the likelihood density, is the maximum likelihood estimate for , its standard deviation, and the derivative is calculated numerically.

More robust uncertainty estimates can be obtained using Markov Chain Monte Carlo (MCMC) sampling

which also allows for the estimation of the covariance matrices

which can be useful when the estimates are used as priors in the modelling process.

The methods also allow the MCMC parameter estimation to be fine tuned by setting the number of iterations, the chain
thinning factor, and the number of burn-in iterations, and the MCMC chain is stored inside the `LDPSet`

instance to
allow for a more detailed convergence analysis.

### 3.4 Log likelihood evaluation

The `LDPSet`

class offers the methods `lnlike_xx`

, where `xx`

is again the model abbreviation, for the
evaluation of the model likelihoods given a set of model coefficients. The joint likelihood for a quadratic model can be
calculated as

where we evaluate the log likelihood with two pairs of coefficients, one pair for each filter in the set considered in the first example.

The likelihoods are calculated assuming that the profile values for a given are normally distributed, which leads to the usual log likelihood equation

(3) |

where is the number of datapoints, are the profile uncertainties, are the profile values, are the model values, and is a factor that can be used to increase the uncertainties based on how much we trust the stellar spectrum models.

### 3.5 Uncertainty multiplier

The uncertainty multiplier is a subjective factor that defines how strongly the limb darkening profile (or the prior created from it) constrains the final analysis (that is, how much we trust the stellar atmosphere models used to create the profiles.) The uncertainty multiplier is applied to log likelihood calculations and model coefficient uncertainty estimation as described in Eq. 3.

### 3.6 Resampling

The sampling in used by Husser et al. (2013) focuses on the detailed sampling of the stellar edge, while the inner
parts of the stellar disk are sampled sparsely. While sampling like this is a good choice for representing the stellar
intensity profile, it may not be optimal for likelihood evaluation and model fitting (at least without sample
weighting), since the fit will be dominated by the very edge. A linear sampling in , for example, may be a better
approach. LDPSet includes a method to resample the profiles to any given vector of or values,
`resample(mu=mu_array,z=None)`

, and also two utility methods `resample_linear_mu(nsamples)`

and
`resample_linear_z(nsamples)`

to resample the profiles linearly in or .

### 3.7 Filters

A filter defines a single passband, and can be either a function or a callable class that returns an array of
transmission values between zero and one given an array of wavelengths in nm. The package comes with some utility
filters, such as the `BoxcarFilter`

used in the first example, and the `TabulatedFilter`

, which allows for
the use of tabulated transmission values listed as a function of wavelength.

### 3.8 Quantum efficiency

Any function that evaluates the quantum efficiency as a function of wavelength (similar to a filter) can be used to
define the instrumental quantum efficiency. LDPSetCreator can be given a detector quantum efficiency curve in the
initialisation, and a `TabulatedQE`

class can be used the same way as `TabulatedFilter`

.

## 4 Discussion

While our ability to model stellar limb darkening is improving due to the advances in the stellar atmosphere models, we still need to develop the ways we integrate the information from these models into our analyses. Current stellar spectrum libraries allow us to calculate limb darkening profiles for custom passbands while propagating the uncertainties in the stellar parameter estimates into the limb darkening profile uncertainties, and we should take advantage of this capability. The marginalization over all the limb darkening profiles allowed by the observations yields the most robust (or, at least, conservative) parameter estimates, but the use of information from the stellar models is also justified in many situations. Transmission spectroscopy, future planet-characterization space missions such as TESS (Ricker et al., 2014) and Plato (Rauer et al., 2014), stellar interferometry, and microlensing studies will all benefit from more advanced ways to integrate the results from stellar atmosphere modelling into the analysis process.

The LDTk-calculated profiles and model coefficients cannot be directly compared with the previous tabulations by Claret et al. (2013) or others due to the differences in handling of the geometry and sampling. Especially the redefinition of the stellar edge (see Espinoza & Jordan, 2015, Sect. 2.2) and the use of uniform sampling in -space will affect the model coefficient estimates (however, these both can be changed by the user). Figure 4 shows LDTk-created profiles for a , , and star observed in , , , and filters with the corresponding quadratic models using the tabulations by Claret et al. (2013). The Claret et al. fits correspond to spherical phoenix models where the intensities corresponding to have been removed from the fit. The Claret et al. (2013) models agree with the LDTk results close to the stellar limb, but deviate significantly for a large span of (and thus for most part of the stellar disk) in the blue passbands.

## 5 Conclusions

We have presented a Python package LDTk that aims to facilitate the inclusion of information from the Husser et al. (2013) stellar atmosphere model spectrum library into astrophysical modelling problems. LDTk offers a modelling approach where the information from the stellar atmosphere models can be used to constrain the stellar limb darkening softly, the strength of the constraint being defined by the researcher’s (subjective) trust to the stellar models. The package can be used to calculate limb darkening coefficients and their uncertainties for freely defined passbands prior to the modelling, or it can be directly integrated into the modelling code. LDTk can also be used to generate limb darkening profile samples directly, bypassing the need to use simple limb darkening models if the modelling approach allows for the use of arbitrary limb darkening profiles (such as the batman transit modelling code, Kreidberg, 2015). This can be useful when the features that are not well captured by the simple limb darkening models, such as the very edge of the star, are important.

## References

- An et al. (2002) An J. H., Albrow M. D., Beaulieu J., Caldwell J. a. R., DePoy D. L., Dominik M., Gaudi B. S., Gould A., Greenhill J., Hill K., Kane S., Martin R., Menzies J., Pogge R. W., Pollard K. R., Sackett P. D., Sahu K. C., Vermaak P., Watson R., Williams A., 2002, ApJ, 572, 521
- Claret (2000) Claret A., 2000, A&A, pp 1081–1190
- Claret (2004) Claret A., 2004, A&A, 1005, 1001
- Claret (2008) Claret A., 2008, A&A, 482, 259
- Claret (2009) Claret A., 2009, A&A, 506, 1335
- Claret & Bloemen (2011) Claret A., Bloemen S., 2011, A&A, 529, A75
- Claret et al. (2014) Claret A., Dragomir D., Matthews J. M., 2014, A&A, 567, A3
- Claret et al. (2012) Claret A., Hauschildt P. H., Witte S., 2012, Astron. Astrophys., 546, A14
- Claret et al. (2013) Claret A., Hauschildt P. H., Witte S., 2013, A&A, 552, A16
- Csizmadia et al. (2013) Csizmadia S., Pasternacki T., Dreyer C., Cabrera J., Erikson A., Rauer H., 2013, A&A, 549, A9
- Espinoza & Jordan (2015) Espinoza N., Jordan A., 2015, MNRAS, 450, 1879
- Fields et al. (2003) Fields D. L., Albrow M. D., An J., Beaulieu J., Caldwell J. a. R., DePoy D. L., Dominik M., Gaudi B. S., Gould A., Greenhill J., Hill K., Jorgensen U. G., Kane S., Martin R., Menzies J., Pogge R. W., Pollard K. R., Sackett P. D., Sahu K. C., Vermaak P., Watson R., Williams A., Glicenstein J., Hauschildt P. H., 2003, ApJ, 596, 1305
- Giménez (2006) Giménez A., 2006, A&A, 450, 1231
- Hayek et al. (2012) Hayek W., Sing D. K., Pont F., Asplund M., 2012, A&A, 539, A102
- Howarth (2011) Howarth I. D., 2011, MNRAS, 418, 1165
- Husser et al. (2013) Husser T.-O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6
- Kipping (2013) Kipping D. M., 2013, MNRAS, 435, 2152
- Kreidberg (2015) Kreidberg L., 2015, ArXiv, 1507.08285
- Magic et al. (2014) Magic Z., Weiss A., Asplund M., 2014, A&A, 573, 10
- Mandel & Agol (2002) Mandel K., Agol E., 2002, ApJ, 580, L171
- Müller et al. (2013) Müller H. M., Huber K. F., Czesla S., Wolter U., Schmitt J. H. M. M., 2013, A&A, 560, A112
- Neilson & Lester (2013a) Neilson H. R., Lester J. B., 2013a, A&A, 556, A86
- Neilson & Lester (2013b) Neilson H. R., Lester J. B., 2013b, A&A, 556, A86
- Rauer et al. (2014) Rauer H., Catala C., Aerts C., Appourchaux T., Benz W., Brandeker A., Christensen-Dalsgaard J., Deleuil M., Gizon L., Goupil M.-J., Güdel M., Janot-Pacheco E., Mas-Hesse M., Pagano I., Piotto G., Pollacco D., Santos Ä., Smith A., Suárez J.-C., Szabó R., Udry S., Adibekyan V., Alibert Y., Almenara J.-M., Amaro-Seoane P., Eiff M. A.-v., Asplund M., Antonello E., Barnes S., Baudin F., Belkacem K., Bergemann M., Bihain G., Birch A. C., Bonfils X., Boisse I., Bonomo A. S., Borsa F., Brandão I. M., Brocato E., Brun S., Burleigh M., Burston R., Cabrera J., Cassisi S., Chaplin W., Charpinet S., Chiappini C., Church R. P., Csizmadia S., Cunha M., Damasso M., Davies M. B., Deeg H. J., Díaz R. F., Dreizler S., Dreyer C., Eggenberger P., Ehrenreich D., Eigmüller P., Erikson A., Farmer R., Feltzing S., de Oliveira Fialho F., Figueira P., Forveille T., Fridlund M., García R. A., Giommi P., Giuffrida G., Godolt M., da Silva J. G., Granzer T., Grenfell J. L., Grotsch-Noels A., Günther E., Haswell C. A., Hatzes A. P., Hébrard G., Hekker S., Helled R., Heng K., Jenkins J. M., Johansen A., Khodachenko M. L., Kislyakova K. G., Kley W., Kolb U., Krivova N., Kupka F., Lammer H., Lanza A. F., Lebreton Y., Magrin D., Marcos-Arenal P., Marrese P. M., Marques J. P., Martins J., Mathis S., Mathur S., Messina S., Miglio A., Montalban J., Montalto M., P. F. G. Monteiro M. J., Moradi H., Moravveji E., Mordasini C., Morel T., Mortier A., Nascimbeni V., Nelson R. P., Nielsen M. B., Noack L., Norton A. J., Ofir A., Oshagh M., Ouazzani R.-M., Pápics P., Parro V. C., Petit P., Plez B., Poretti E., Quirrenbach A., Ragazzoni R., Raimondo G., Rainer M., Reese D. R., Redmer R., Reffert S., Rojas-Ayala B., Roxburgh I. W., Salmon S., Santerne A., Schneider J., Schou J., Schuh S., Schunker H., Silva-Valio A., Silvotti R., Skillen I., Snellen I., Sohl F., Sousa S. G., Sozzetti A., Stello D., Strassmeier K. G., Švanda M., Szabó G. M., Tkachenko A., Valencia D., Van Grootel V., Vauclair S. D., Ventura P., Wagner F. W., Walton N. A., Weingrill J., Werner S. C., Wheatley P. J., Zwintz K., 2014, Exp. Astron., 38, 249
- Ricker et al. (2014) Ricker G. R., Winn J. N., Vanderspek R., Latham D. W., Bakos G. A., Bean J. L., Berta-Thompson Z. K., Brown T. M., Buchhave L., Butler N. R., Butler R. P., Chaplin W. J., Charbonneau D., Christensen-Dalsgaard J. r., Clampin M., Deming D., Doty J., De Lee N., Dressing C., Dunham E. W., Endl M., Fressin F., Ge J., Henning T., Holman M. J., Howard A. W., Ida S., Jenkins J., Jernigan G., Johnson J. A., Kaltenegger L., Kawai N., Kjeldsen H., Laughlin G., Levine A. M., Lin D., Lissauer J. J., MacQueen P., Marcy G., McCullough P. R., Morton T. D., Narita N., Paegert M., Palle E., Pepe F., Pepper J., Quirrenbach A., Rinehart S. A., Sasselov D., Sato B., Seager S., Sozzetti A., Stassun K. G., Sullivan P., Szentgyorgyi A., Torres G., Udry S., Villasenor J., 2014 Transiting Exoplanet Survey Satellite (TESS)
- Sing (2010) Sing D. K., 2010, A&A, 510, A21
- White et al. (2013) White T. R., Huber D., Maestro V., Bedding T. R., Ireland M. J., Baron F., Boyajian T. S., Che X., Monnier J. D., Pope B. J. S., Roettenbacher R. M., Stello D., Tuthill P. G., Farrington C. D., Goldfinger P. J., McAlister H. A., Schaefer G. H., Sturmann J., Sturmann L., ten Brummelaar T. A., Turner N. H., 2013, MNRAS, 433, 1262