The ExoTETHyS package: Tools for Exoplanetary Transits around Host Stars
Giuseppe Morello, Antonio Claret, Marine Martin-Lagarde, Christophe Cossou, Angelos Tsiaras, Pierre-Olivier Lagage
aa r X i v : . [ a s t r o - ph . E P ] F e b Draft version February 4, 2020
Typeset using L A TEX twocolumn style in AASTeX62
The ExoTETHyS package: Tools for Exoplanetary Transits around Host Stars
G. Morello,
A. Claret,
4, 5
M. Martin-Lagarde, C. Cossou, A. Tsiaras, and P.-O. Lagage AIM, CEA, CNRS, Universit´e Paris-Saclay, Universit´e Paris Diderot, Sorbonne Paris Cit´e, F-91191 Gif-sur-Yvette, France INAF–Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy Dept. of Physics & Astronomy, University College London, Gower Street, WC1E 6BT London, UK Instituto de Astrof´ısica de Andaluc´ıa, CSIC, Apartado 3004, 18080 Granada, Spain Dept. F´ısica Te´orica y del Cosmos, Universidad de Granada, Campus de Fuentenueva s/n, 10871, Granada, Spain Institut d’Astrophysique Spatiale, CNRS/Universit´e Paris-Sud, Universit´e Paris-Saclay, bˆatiment 121, Universit´e Paris-Sud, 91405Orsay Cedex, France (Published AJ)
ABSTRACTWe present here the first release of the open-source python package
ExoTETHyS (stable: https://zenodo.org/badge/latestdoi/169268509, development version:https://github.com/ucl-exoplanets/ExoTETHyS/), which aims to provide a stand-alone set oftools for modeling spectrophotometric observations of the transiting exoplanets. In particular, wedescribe: (1) a new calculator of stellar limb-darkening coefficients that outperforms the existingsoftware by one order of magnitude in terms of light curve model accuracy, i.e., down to <
10 partsper million, and (2) an exact transit light curve generator based on the entire stellar intensityprofile rather than limb-darkening coefficients. New tools will be added in later releases to modelvarious effects in exoplanetary transits and eclipsing binaries.
ExoTETHyS is a reference package forhigh-precision exoplanet atmospheric spectroscopy with the upcoming
James Webb Space Telescope and
Atmospheric Remote-sensing Infrared Exoplanet Large-survey missions.
Unified Astronomy Thesaurus concepts:
Observational astronomy (1145); Spectroscopy (1558); Hightime resolution astrophysics (740); Exoplanet atmospheres (487); Stellar atmospheres (1584); Transitphotometry (1709); Eclipsing binary stars (444); Exoplanet systems (484); Stellar astronomy (1583);Limb darkening (922)Software reviewed by the Journal of Open Source Software INTRODUCTIONMore than 3000 transiting exoplanets have been dis-covered in the last 20 years. The number of transit-ing exoplanets accounts for about three-quarters of thecurrent exoplanet census , although this large fractionis due to targeted research programs rather than beinga random sample from the exoplanet population. Thesuccess of the transit method is due to several contribut-ing factors, including its ability to characterize them ingreat detail. A transit is revealed by a decrement in fluxwhile the planet occults part of the stellar disk. Themain observables are the transit depth and durations, Corresponding author: Giuseppe [email protected] source: https://exoplanetarchive.ipac.caltech.edu leading to measurements of the exoplanet size, orbitalsemimajor axis and inclination, and stellar mean density(Seager & Mall´en-Ornelas 2003). Transit spectroscopyis now routinely used to investigate the chemistry andphysics of exoplanet atmospheres, through differencesin transit depth of ∼ Morello et al.
Name Geometry a Range T eff(K) Range log g Range [
M/H ] Range λ ( µ m) Reference ATLAS
P-P 3500-50000 0.0-5.0 –5.0-1.0 0.009-160.0 Claret (2000)
PHOENIX
PHOENIX a Geometry types: P-P=plane-parallel; S1=spherical 1D
Table 1.
Stellar model atmosphere grids available with the first release of ExoTETHyS (Akinsanmi et al. 2019; Hellard et al. 2019). Amongthe nonstellar effects, the exoplanet nightside emissioncan also play a significant role (Kipping & Tinetti 2010;Morello et al. 2019).The
ExoTETHyS package is conceived as a toolboxfor those who analyze the exoplanetary transits. Thefirst release focuses on the tools for modeling the stel-lar limb-darkening effect, the importance of which isubiquitous in transit observations, as well as in opti-cal interferometry, microlensing, and eclipsing binaryobservations. Future versions of
ExoTETHyS will in-clude useful tools for modeling other effects, as wellas for estimating their impact on specific observations,based on the astrophysical system parameters, the in-strument passband, and the noise level. Accurate mod-eling of all of the aforementioned effects proved to becrucial in the analysis of several
CoRoT and
Kepler objects (e.g., Mazeh & Faigler 2010; Barnes et al. 2011;Mazeh et al. 2012; Masuda 2015; Howarth & Morello2017; Reinhold et al. 2017; Shporer 2017; Nielsen et al.2019), because of the high-precision photometry downto the .
10 ppm level (Christiansen et al. 2012). Asimilar photometric precision is expected for someof the ongoing
Transiting Exoplanet Survey Satellite ( TESS ) observations (Ricker et al. 2014), future obser-vations with the
CHaracterising ExOPlanet Satellite ( CHEOPS ; Isaak & Benz 2019) and
PLAnetary Tran-sits and Oscillations ( PLATO ; Rauer et al. 2014), andin spectroscopic observations with the upcoming
JamesWebb Space Telescope ( JWST ; Beichman et al. 2014)and
Atmospheric Remote-sensing Infrared ExoplanetLarge-survey ( ARIEL ; Pascale et al. 2018) space mis-sions.Stellar limb-darkening is the wavelength-dependentradial decrease in specific intensity. Consequently, thetransit light curve deviates from the flat-bottomed shapethat would be observed in the case of a uniform stellardisk; the difference signal can be as large as ∼ ppmfor the transit of a hot Jupiter observed at UV or visiblewavelengths. Typically, the radial intensity distributioncomputed from specific stellar atmosphere models isparameterized by a set of limb-darkening coefficients, which are fixed in the analyses of transit light curves.Many researchers have produced multiple grids of stel-lar atmosphere models with different codes, then usedto compile precalculated tables of limb-darkening coef-ficients (e.g., Claret 2000, 2003, 2004, 2008, 2017, 2018;Sing 2010; Howarth 2011a; Claret & Bloemen 2011;Claret et al. 2012, 2013, 2014; Neilson & Lester 2013a,b;Magic et al. 2015; Reeve & Howarth 2016). The lack ofempirical validation for stellar limb-darkening preventsthe final choice of the most reliable model(s). Thepresence of unocculted stellar spots during an exoplan-etary transit may alter the effective limb-darkening co-efficients, which will be slightly different from thosecalculated for the case of unspotted stellar surface(Csizmadia et al. 2013).In some cases, significantly different parametric in-tensity profiles have been obtained from the samemodel atmosphere, depending on the sampling ofthe model intensity profile, the functional form (so-called limb-darkening law), and/or the fitting algorithmadopted (Claret 2000; Heyrovsk´y 2007; Howarth 2011b;Espinoza & Jord´an 2015). The system parameters ob-tained from the light curve fits with the alternative setsof limb-darkening coefficients can vary by more thanthe respective 1 σ error bars, typically, if the relativephotometric precision of the observations is of the orderof (or better than) 100 ppm per minute interval.In this paper, we probe an optimized fitting algorithmfor the limb-darkening coefficients that minimizes thedifference between (numerically integrated) referencelight curves and the corresponding approximated transitmodels with limb-darkening coefficients. Therefore weeliminate the degeneracy from the choice between sev-eral fitting algorithms that were leading to significantlydifferent parametric profiles for the same stellar atmo-sphere model (e.g., Espinoza & Jord´an 2015). The high-fidelity match between the stellar intensity profiles andthe transit light curve models facilitates comparativestudies of the model atmospheres, especially with theincreasing number of observations with a spectrophoto-metric precision down to ∼
10 ppm (e.g.,
CoRoT , Kepler , TESS , and
Hubble Space Telescope ( HST )/WFC3 data). ample article
Structure of the paper
Section 2 provides a technical description of the
ExoTETHyS package and the algorithms adopted. Sec-tion 3 discusses the precision of the limb-darkeningcalculator for the analysis of exoplanetary transits. Inparticular, Section 3.1 compares various algorithms thatare adopted in the other publicly available codes andtheir variants, Section 3.2 compares the performancesof the alternative limb-darkening laws, and Section 3.3provides a formula to estimate the potential error in thetransit model based on the goodness of fit for the limb-darkening coefficients that should be compared with thenoise level in the observations. Section 4 discusses themain functionality of the
ExoTETHyS package, its cur-rent and future usage. Finally, Section 5 summarizesthe key points discussed in this paper. DESCRIPTION OF THE
EXOTETHYS
PACKAGEThe first release of
ExoTETHyS includes the followingsubpackages:1. Stellar Atmosphere Intensity Limb (SAIL), whichcan calculate the limb-darkening coefficients forspecific stellar targets or over predetermined pa-rameter grids;2. Transit Ring-Integrated Profile (TRIP), which cancompute an exact transit light curve by direct in-tegration of the occulted stellar flux, without us-ing an analytical function (limb-darkening law) toapproximate the stellar intensity profile.The TRIP subpackage was conceived to model exoplane-tary transits. Following requests by users, we are addinga function to model eclipsing binaries.2.1.
The SAIL subpackage
The SAIL subpackage is a generic stellar limb-darkening calculator that is not specific to a prede-termined list of instruments or standard passbands. Itis conceptually similar to the calculator provided byEspinoza & Jord´an (2015), but with different features.A technical difference is the use of a novel fitting al-gorithm for obtaining the limb-darkening coefficients,specifically optimized for modeling the exoplanetarytransits, instead of multiple algorithm options with un-clear performances (see Sections 2.1.4 and 3.1).2.1.1.
Input and output
The SAIL subpackage requires a configuration file tospecify the desired calculation. The user can choose ei-ther “individual” or “grid” calculation type. The firstoption enables calculation of the limb-darkening coeffi-cients for a star or for a list of stars with the parameters specified by the user, while the latter will provide thelimb-darkening coefficients for a grid of precalculatedstellar model atmospheres. In both cases, the user mustselect one of the available stellar model grids, which werecomputed with different codes and settings (see Table 1and references therein). For each grid, the stellar mod-els are identified by a set of three parameters, i.e., theeffective temperature ( T eff ), the surface gravity (log g ),and the metallicity ([ M/H ]). As the limb-darkening co-efficients are mostly dependent on the effective tempera-ture, the user must provide the effective temperatures ofall the individual stars. The other parameters have de-fault values of log g =4.5 and [ M/H ] =0.0, correspond-ing to a main-sequence star with solar abundances, ifthey are not given by the user. For the grid calcula-tion type, the default option is to calculate the limb-darkening coefficients for all the stellar models in theselected database. Alternatively, the user can select asubgrid by specifying the minimum and/or maximumvalues for each stellar parameter.Another key input is the passband, i.e., the total spec-tral response of the observing instrument. For most in-struments, the spectral response is available as a tableof photon-to-electron conversion factors at given wave-lengths. The limb-darkening coefficients do not dependon the absolute values of the spectral response, so thata scaled/normalized version of the spectral response willgive identical results. The spectral responses of the mostcommon instruments for transiting exoplanets are builtinto the package. The code can accept any user-definedpassband with the same file format. It is also possi-ble to calculate the limb-darkening coefficients for mul-tiple wavelength bins within a given passband by spec-ifying the two wavelengths delimiting each bin. Thisoption is particularly useful for exoplanet spectroscopicobservations, such as those currently performed with
HST /WFC3.The last mandatory input in the configuration fileis the list of limb-darkening laws to adopt (at leastone). The code includes several built-in limb-darkeninglaws, including all of the most commonly used (see Sec-tion 2.1.3), but it can also accept user-defined laws.The “basic” outputs are python dictionaries contain-ing the best-fit limb-darkening coefficients obtained forthe required passbands, wavelength bins, and limb-darkening laws. The output dictionaries also providethe corresponding weighted rms of the fitting residualsto allow for a quick quality check (see Section 3.3). Forthe case of individual calculation type, the results ob-tained for each target are stored in separate pickle files.Optionally, the user can request a “complete” output,whic includes intermediate products such as the numeric
Morello et al. intensity profiles at various stages of the calculation (seeSections 2.1.2-2.1.5). The additional information of thecomplete output is offered, mainly, as a way to iden-tify bugs in the code and/or issues with certain stellarmodel atmospheres and wavelengths. Usually, the exo-planetary scientists will be interested to the basic outputonly.2.1.2.
From the stellar model atmospheres to thepassband-integrated intensities
The stellar model atmosphere grids consist of one filefor each triple of stellar parameters ( T eff , log g , [ M/H ]),providing the specific intensities ( I λ ( µ )) in units of ergcm − s − ˚A − sr − at several positions on the sky-projected stellar disk over a given spectral range. Forhistorical reasons, the independent variable is µ = cos θ ,where θ is the angle between the line of sight and the cor-responding surface normal. The radial coordinate in thesky-projected disk is r = p − µ , where r = 1 ( µ = 0)corresponds to the spherical surface radius. Table 1 re-ports the information about the databases available withthe first release of ExoTETHyS. We refer to the relevantpapers and references therein for comparisons betweenthe models. The passband-integrated intensities are cal-culated as I pass ( µ ) ∝ Z λ λ I λ ( µ ) R pass ( λ ) λdλ, (1)where R pass ( λ ) is the spectral response of the instrumentin electrons photon − , and λ and λ are the passbandor wavelength bin limits. The passband-integrated in-tensities are obtained in units proportional to electronscm − s − sr − . As the limb-darkening coefficients arenot affected by the (omitted) proportionality factor inEquation 1, the final intensities are normalized such that I pass ( µ = 0) = 1.The intensity profiles, I λ ( µ ), have distinctive behav-iors depending on the plane-parallel or spherical ge-ometry adopted by the selected grid of model atmo-spheres. In particular, the spherical intensity profilesshow a steep drop-off close to the stellar limb, which isnot observed in the plane-parallel models. The expla-nation for the different behaviors is exhaustive in theliterature (Wittkowski et al. 2004; Espinoza & Jord´an2015; Morello et al. 2017). The almost null intensitiesat small µ are integrated over lines of sight that inter-sect only the outermost atmospheric shells, which havethe smallest emissivity. Here µ =0 ( r =1) correspondsto the outermost shell of the model atmosphere, whichis typically outside the stellar radius that would be ob-served in transit. Our algorithm calculates the photo-metric radius at the inflection point of the spherical in-tensity profile, i.e., where the gradient | dI ( r ) /dr | is the maximum (Wittkowski et al. 2004; Espinoza & Jord´an2015). The radial coordinates are then rescaled suchthat r = 1 ( µ =0) at the photometric radius, and thoseintensities with rescaled r > Limb-darkening laws
A long list of analytical forms, so-called limb-darkening laws, has been proposed in the literature toapproximate the stellar intensity profiles. The followingoptions are built in the package:1. the linear law (Schwarzschild 1906), I λ ( µ ) = 1 − a (1 − µ ); (2)2. the quadratic law (Kopal 1950), I λ ( µ ) = 1 − u (1 − µ ) − u (1 − µ ) ; (3)3. the square-root law (Diaz-Cordoves & Gimenez1992), I λ ( µ ) = 1 − v (1 − √ µ ) − v (1 − µ ); (4)4. the power-2 law (Hestroffer 1997), I λ ( µ ) = 1 − c (1 − µ α ); (5)5. the four-coefficient law (Claret 2000), hereinafterreferred to as claret-4, I λ ( µ ) = 1 − X k =1 a n (1 − µ k/ ); (6)6. a generalized n th -degree polynomial law, I λ ( µ ) = 1 − n X k =1 b k (1 − µ k ); (7)7. a generalized claret- n law, I λ ( µ ) = 1 − n X k =1 c k (1 − µ k/ ); (8)Additionally, user-defined limb-darkening laws can beeasily implemented. We recommend using the claret-4 law to achieve a model precision of .
10 ppm in theanalysis of exoplanetary transits (see Section 3.2). Thenext release of
ExoTETHyS will include a grid of whitedwarf models, for which we have also found the claret-4 law to be significantly more accurate than the two-coefficient laws (Claret et al. 2020). ample article
From the passband-integrated intensities to thelimb-darkening coefficients
The limb-darkening coefficients are obtained througha weighted least-squares fit of the passband-integratedintensity profile with weights proportional to the sam-pling interval in r , hereinafter referred to as weighted - r fit. The corresponding cost function is the weighted rmsof residuals, weighted - r rms = P ni =1 w i ( I pass ( µ i ) − I lawpass ( µ i )) P ni =1 w i ! , (9)with weights w i = (1 − r ) + 0 . r − r ) , if i = 10 . r i − − r i +1 ) , if 1 < i < n . r n − , if i = n , (10)where the r i are arranged in descending order, and r n = 0. The choice of cost function is optimized for thestudy of exoplanet transits, as detailed in Section 3.1.The performances of the spherical model fits are furtherenhanced by discarding those points with r > . weighted - r QS fit. Further details about the alternative fittingprocedures are discussed in Section 3.1.2.1.5.
Interpolation from the grid of stellar models
The process described in Sections 2.1.2-2.1.4 enablesthe calculation of limb-darkening coefficients for thestellar-atmosphere models contained in the grid, startingfrom their precalculated specific intensities. The limb-darkening coefficients for an individual target with ageneric set of stellar parameters are obtained by sequen-tial linear interpolation through the following steps:1. identification of the neighbors in the model-grid,i.e., the vertices of the cube in parameter spacethat contains the requested model (maximum 8models);2. calculation of the limb-darkening coefficients foreach of the neighbors;3. interpolation in [
M/H ] between models with thesame T eff and log g , leading to a maximum of 4 setsof limb-darkening coefficients with the requested[ M/H ]; 4. interpolation in log g between the above calculatedsets of coefficients with the same T eff , leading to amaximum of 2 sets of limb-darkening coefficientswith the requested log g and [ M/H ];5. interpolation in T eff between the above calculatedsets of coefficients.We note that this sequential interpolation is possiblebecause of the regularity of the model grids.2.2. The TRIP subpackage
The TRIP subpackage is used to generate exact transitlight curves by direct integration of the occulted stellarflux at given instants. It assumes a dark spherical planettransiting in front of a spherically-symmetric (unspot-ted) star. In this simple case, the normalized flux (i.e.,relative to the stellar flux) is a function of two geometricvariables, as reported by Mandel & Agol (2002), and ofthe stellar intensity profile: F ( p, z, I ( µ )) = 1 − Λ( p, z, I ( µ )) , (11)where p is the planet-to-star radii ratio ( p = R p /R ∗ ), z is the sky-projected distance between the centers ofthe two spheres in units of the stellar radius, and I ( µ )is the stellar intensity profile. TRIP does not use ananalytical approximation of the limb-darkening profile,unlike most transit light curve calculators such as thoseprovided by Mandel & Agol (2002), Gim´enez (2006),Agol et al. (2019), JKTEBOP (Southworth et al. 2004),
TAP (Gazak et al. 2012),
EXOFAST (Eastman et al. 2013),
PyTransit (Parviainen 2015),
BATMAN (Kreidberg 2015),and
PYLIGHTCURVE (Tsiaras et al. 2016).2.2.1. Input and output
The TRIP subpackage requires a configuration file,where the user has to specify the name of the text filescontaining the limb-darkening profile, the phase, time,or z -series for which to calculate the normalized flux,and a list of parameter values that includes p and thoseparameters eventually needed to compute the z -series(see Section 2.2.2). The limb-darkening file consists oftwo columns with the µ or r values (first column) andthe corresponding specific intensities (second column).A list of optional parameters can be used to set the cal-culation details, i.e., the number of annuli, the interpo-lation variable, and the polynomial order for the splineinterpolation (see Section 2.2.3). It is also possible todefine simple operations on the original limb-darkeningprofile, i.e., a possible cutoff in µ or r with or without https://github.com/ucl-exoplanets/pylightcurve Morello et al. I λ ( μ ) modelunweightedweighted-μweighted-μ QS I λ ( μ ) R e s i d u a l s RMS=3.45e-01, w-μ RMS=2.72e-02 RMS=2.89e-01, w-μ RMS=2.26e-02 RMS=2.08e-01, w-μ RMS=1.96e-02 0.0 0.2 0.4 0.6 0.8 1.0μ w-μ RMS=1.92e-03, w-μ QS RMS=1.33e-05 w-μ RMS=1.69e-03, w-μ QS RMS=4.00e-04 w-μ RMS=1.03e-02, w-μ QS RMS=1.04e-02
HD209458, 7.59-7.61 μm, claret-4
Figure 1.
Example with a model intensity distribution for a star similar to HD209458 ( T eff = 6100 K, log g = 4 . µ m wavelength range, by using the PHOENIX µ from the stellar atmosphere model (black circles), unweighted (gray), weighted - r (orange), and weighted - r QS (red) model fits with claret-4 coefficients. The vertical dashed line denotes the cutoff value for the quasi-spherical fit (seeSection 3.1). Top, right panel: analogous plot vs. r . Bottom panels: residuals between the fitted and model intensity values. Thecorresponding unweighted and weighted rms amplitudes of residuals are also reported. Note that, in this case, the unweightedleast-squares fit leads to a non-monotonic radial intensity profile, which is physically unexpected. rescaling the µ or r values to the cutoff radius. The out-put is a text or pickle file containing the normalized fluxseries for the requested phase, time, or z -series.2.2.2. Computing the z -series In general, z is a function of the orbital phase (Φ),i.e., the fraction of orbital period ( P ) from the closesttransit event: Φ = t − E.T. P − n, (12)where t denotes time, E.T. is the Epoch of Transit (i.e.,a reference mid-transit time), and n is the number of or-bits from the E.T. rounded to the nearest integer. Con-ventionally, Φ values are in the range of [ − . , .
5] or[0 ,
1] and Φ = 0 at mid-transit time. The projected star–planet separation is given by z = a R q − cos (2 π Φ) sin i circular orbit a R − e e cos f q − sin ( f + ω ) sin i eccentric orbit , (13)where a R is the orbital semimajor axis in units of thestellar radius, i is the inclination, e is the eccentricity, ω is the argument of periastron, and f is the true anomaly.In the eccentric case, the true anomaly is calculated fromthe orbital phase by solving Kepler’s equation, π − ω + 2 π Φ = E − e sin E, (14)then f = 2 arctan r e − e tan E ! . (15)2.2.3. Calculating the normalized flux ample article −0.03 −0.02 −0.01 0.00 0.01 0.02 0.030.9850.9900.9951.000 referenceweighted-r QS Figure 2.
Top panel: simulated transit light curve (black)of HD209458 b as it would be observed by
TESS , and best-fit model with claret-4 limb-darkening coefficients obtained withthe weighted - r QS method (red). Bottom panels: residuals be-tween the reference light curve and the best-fit models withclaret-4 limb-darkening coefficients obtained with different limb-darkening laws and fitting methods (see Section 3.1). The peak-to-peak and rms amplitudes of the residuals are reported. −2021e−4 PEAK-TO-PEAK=2 ppm RMS=0 ppmclaret-4we ghted-r QS PEAK-TO-PEAK=48 ppm RMS=12 ppmpower2 PEAK-TO-PEAK=39 ppm RMS=9 ppmquadrat c−2021e−4 PEAK-TO-PEAK=9 ppm RMS=2 ppmwe ghted-r PEAK-TO-PEAK=47 ppm RMS=12 ppm PEAK-TO-PEAK=53 ppm RMS=13 ppm−2021e−4 PEAK-TO-PEAK=91 ppm RMS=24 ppmwe ghted-μ PEAK-TO-PEAK=75 ppm RMS=17 ppm PEAK-TO-PEAK=197 ppm RMS=57 ppm−2021e−4 PEAK-TO-PEAK=182 ppm RMS=48 ppmunwe ghted PEAK-TO-PEAK=110 ppm RMS=27 ppm PEAK-TO-PEAK=294 ppm RMS=84 ppm−202 N o r m a l z e d f l u x The total and occulted stellar flux are given, respec-tively, by the integrals F ∗ = Z I ( r ) 2 πr dr, (16)and F ∗ , occ = Z I ( r ) 2 πr f p,z ( r ) dr, (17) with f p,z ( r ) = π arccos r + z − p zr | z − p | < r < z + p r ≤ z − p or r ≥ z + p r ≤ p − z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ r ≤ . (18) Morello et al. P e a k - t o - p e a k ( pp m ) Sphericalclaret-4weighted-rweighted-μunweightedinterp-μ 100 interp-μ 1000interp-r 100interp-r 1000 0 2 4 6 8 10Wavelength (μm)0.02.55.07.510.012.515.017.520.0 Quasi-spherical claret-4weighted-r QSweighted-μ QSunweighted QSweighted-r
Best-fit spectral light-curve residuals
Figure 3.
Peak-to-peak of residuals between the reference spectral light curves for the transit of HD209458 b and the best-fitmodels with claret-4 limb-darkening coefficients obtained with different fitting methods (see Section 3.1). Left panel: resultsobtained with the spherical methods, i.e., taking into account the whole spherical intensity profiles. Right panel: resultsobtained with the quasi-spherical methods, i.e., with a cutoff of r ≤ . weighted - r method (dotted, orange line).The unweighted QS (gray line) and the weighted - µ QS (green line) overlap in the plot. Note the scale difference between thetwo panels. I ( r ) is the specific intensity at the normalized radial co-ordinate r = p − µ , and f p,z ( r ) is the fraction of cir-cumference with radius r covered by the planet. Equa-tions 16 and 17 rely on the assumed spherical symmetryfor the star; Equation 18 also makes use of the planetsphericity. Finally, the normalized flux is given by Equa-tion 11 with Λ( p, z, I ( µ )) = F ∗ , occ F ∗ . (19)The integrals in Equations 16 and 17 are calculated nu-merically by using the mid-point rule with a uniformpartition in r . The specific intensities are evaluated atthe partition radii by interpolating in µ or r from the in-put limb-darkening profiles. The TRIP algorithm withdefault settings is identical to the “tlc” described byMorello et al. (2017). PERFORMANCE OF
EXOTETHYS
Comparison between fitting algorithms for thestellar intensity profiles
A long list of methods has been adopted in the liter-ature for fitting the limb-darkening laws to the modelintensity profiles leading to significantly different limb-darkening coefficients. The coefficients obtained witha simple least-squares fit depend on the spatial distri-bution of the precalculated intensities. The effect ofsampling is particularly evident for the
PHOENIX profilesbecause of a much finer sampling near the drop-off re- gion. For example, Figure 1 shows the case of a starsimilar to HD209458 in the mid-infrared, for which thesimple least-squares solution presents a non-monotonic(unphysical) profile with unexpected undulations. Inthis paper, we compare the following fitting procedures:1. unweighted , i.e., simple least-squares fit;2. weighted - r , i.e., weighted least-squares fit withweights proportional to the sampling interval in r , as detailed in Equations 9 and 10;3. weighted- µ , i.e., weighted least-squares fit withweights proportional to the sampling interval in µ ;4. interp- µ , i.e., least-squares fit on the intensi-ties interpolated over 100 µ values with a uniformseparation in µ , as suggested by Claret & Bloemen(2011);5. interp- µ , i.e., least-squares fit on the intensi-ties interpolated over 1000 µ values with a uniformseparation in µ ;6. interp- r , i.e., least-squares fit on the in-tensities interpolated over 100 r values witha uniform separation in r , as suggested byParviainen & Aigrain (2015) (with an unspecifiednumber of interpolated values); ample article −80−60−40−20020406080100 ( p − . ) * e Sphericalclaret-4weighted-rweighted-μunweightedinterp-μ 100 interp-μ 1000interp-r 100interp-r 1000rescaled input 10.012.515.017.520.022.525.027.530.0 Quasi-spherical claret-4weighted-r QSweighted-μ QSunweighted QSweighted-rrescaled input−0.50.00.51.01.52.0 b − . μ − s Best-fit spectral transit parameters
Figure 4.
Best-fit transit parameters to the reference spectral light curves for the transit of HD209458 b assuming claret-4limb-darkening coefficients obtained with different fitting methods (see Section 3.1). The true parameter values are reported inblack. Left panels: results obtained with the spherical methods, i.e., taking into account the whole spherical intensity profiles.Right panels: Results obtained with the quasi-spherical methods, i.e., with a cutoff of r ≤ . weighted - r method(dotted, orange line). Note the scale difference between the two panels. Morello et al. interp- r , i.e., least-squares fit on the intensi-ties interpolated over 1000 r values with a uniformseparation in r ;8. unweighted QS, i.e., least-squares fit with a cutoff r ≤ . weighted - r QS, i.e., analogous to weighted - r witha cutoff r ≤ . weighted - µ QS, i.e., analogous to weighted - µ witha cutoff r ≤ . µ ≥ . PHOENIX models with the original µ values. In thiswork, we redefine the cutoff using the rescaled r , suchthat it corresponds to the same fraction of the photomet-ric stellar radius for all the models (see Section 2.1.2).Our new definition with r ≤ . PHOENIX models incor-porated in the
ExoTETHyS package also include modelsof stellar atmospheres with lower gravities than thoseanalyzed by Claret et al. (2012), corresponding to sub-giant, giant, and supergiant stars. For some of thesemodels, the intensity drop-off occurs at µ > .
1, so thatthe cutoff of µ ≥ . TESS passband when adoptingthe different sets of limb-darkening coefficients. The weighted - r QS method implemented in
ExoTETHyS .SAILgives the smallest residuals, with a peak-to-peak of2 ppm and rms amplitude below 1 ppm. The other QSmethods, weighted - µ QS and unweighted
QS, lead to al-most identical residuals, with a peak-to-peak of 3 ppm.Among the spherical methods, the weighted - r gives thesmallest residuals with a peak-to-peak of 9 ppm and rmsamplitude of 2 ppm, followed by the interp- r and interp- r with about 1.5 and 2 times larger residualamplitudes, respectively. All the other methods lead tosignificantly larger residuals of tens to a few hundredppm, which are comparable with the predicted noisefloor of 60 ppm for the TESS observations (Ricker et al.2014). Figure 3 shows the peak-to-peak of the residuals forthe same transit as a function of wavelength, based onsimulated light curves with 20 nm passband widths.This spectral analysis confirms the relative ranking ofthe fitting methods derived from the
TESS simulation.In particular, the weighted - r QS method leads to a peak-to-peak of residuals below 2 ppm at wavelengths longerthan 1 µ m, and overall below 8 ppm. The other quasi-spherical methods are marginally worse than weighted - r QS at wavelengths shorter than 2 µ m, but the worstcase peak-to-peak of residuals is less than 13 ppm. The weighted - r method leads to peak-to-peak of residualsin the range of 5-15 ppm, with a sawtooth-like mod-ulation in wavelength. We noted that the small, butabrupt, jumps that occur at certain wavelengths corre-spond to changes of the inflection point in the stellarintensity profile as defined in Section 2.1.2. The samephenomenon occurs for all the other spherical modelswith larger sawtooth-like modulations. It may appearsurprising that the peak-to-peak of residuals obtainedwith the spherical methods tends to be larger at thelonger wavelengths, for which the limb-darkening effectis expected to be smaller. The cause of the poor per-formances of most spherical methods in the infrared isthe intensity drop-off, which is typically steeper thanthe drop-off in the UV and visible. Such drop-off hasa negligible effect in the numerically integrated transitlight curves, hence the better performances of the QSfits.Figure 4 shows the best-fit transit parameters corre-sponding to the same spectral light curves, and com-pared with the respective input parameters correctedfor the rescaled r (see Section 2.1.2). We retrieved thecorrect transit depth within 5 ppm, the impact param-eter within 6 × − , and the transit duration within1 s at all wavelengths, when using the weighted - r orQS limb-darkening coefficients. However, slightly largerspectral trends appear in these parameters because ofthe wavelength-dependent stellar radius. The peak-to-peak variation in transit depth over the spectral range of0.25–10 µ m is 10 ppm. The other sets of limb-darkeningcoefficients introduce orders-of-magnitude larger biasesin the retrieved transit parameters, also larger spectralsawtooth-like modulations in the infrared (few tens ofppm in transit depth across 1-10 µ m), and severe dis-crepancies between the parameter values obtained in theUV/visible and those obtained in the infrared.3.2. Performance of the limb-darkening laws
Figure 5 compares the peak-to-peak of the spectrallight curve residuals when adopting the limb-darkeningcoefficients calculated by
ExoTETHyS .SAIL for differ- ample article Residuals W e i g h t e d - r Q S R M S P e a k - t o - p e a k ( pp m ) Figure 5.
Top, left panel: peak-to-peak of residuals between the reference spectral light curves for the transit of HD209458 band the best-fit models using the limb-darkening coefficients calculated for the different laws (see Section 2.1.3). Top, rightpanel: weighted - r QS rms of residuals to the model intensity profiles. Bottom panels: zoom-in of the panels above.Wavelength range Claret-4 Power-2 Quadratic Square-rootMaximum bias 0.25–10.0 µ m 5 165 235 174(ppm) < µ m 4 165 235 174 > µ m 5 19 27 18 > µ m 3 4 10 5Rms bias 0.25–10.0 µ m 1 20 20 23(ppm) < µ m 2 71 62 81 > µ m 1 5 11 4 > µ m 1 2 6 2Spectrum 0.25–10.0 µ m 10 177 258 341peak-to-peak < µ m 7 177 254 341(ppm) > µ m 10 27 17 25 > µ m 2 3 4 3Spectrum std 0.25–10.0 µ m 2 20 18 23(ppm) < µ m 1 45 58 64 > µ m 2 7 4 5 > µ m < < < < Table 3.
Spectral analysis of the error in transit depth when adopting different limb-darkening laws. ent limb-darkening laws, as well as the corresponding weighted - r QS rms of the residuals to the stellar intensityprofiles. The correlation between the two goodness-of-fit measures is explored in Section 3.3. At wavelengths & µ m, the precision of the power-2 and square-rootlimb-darkening coefficients is comparable to that of theclaret-4 coefficients, resulting in light curve residuals be-low 5 ppm. While the claret-4 law performs similarlywell even at shorter wavelengths, the two-coefficient lawslead to larger light curve residuals up to ∼
100 ppm inthe UV and visible. The quadratic law is less precise, leading to light curve residuals above 25 ppm even at10 µ m.Figure 6 shows the fitted transit parameters and theirexpected values. Typically, the bias in transit depth isof the same order of magnitude of the light curve resid-uals, but it can be both larger or smaller than theirpeak-to-peak amplitudes owing to parameter degenera-cies. Table 3 reports the statistics of the errors in transitdepth obtained with the different limb-darkening lawsacross given spectral ranges. The maximum bias in tran-sit depth at 5–10 µ m is within 10 ppm for any limb-2 Morello et al. −200−1000100200300 ( p − . ) * e claret-4power-2quadratic square-rootrescaled input01234 b − . μ − s Be t-fit spectral transit parameters
Figure 6.
Best-fit transit parameters to the reference spectral light curves for the transit of HD209458 b using the limb-darkening coefficients calculated for the different laws (see Section 2.1.3). The true parameter values are reported in black. ample article P e a k - t o - p e a k ( pp m ) claret-4power-2quadraticsquare-rootlinear fit Figure 7.
Weighted - r QS rms of residuals to the modelintensity profiles vs. peak-to-peak of the transit light curveresiduals for the spectral templates of HD20458 b adoptingdifferent limb-darkening laws. The black line is the globallinear fit. darkening parameterization, which is just below theminimum photon noise floor for
JWST /Mid-InfraRedInstrument (MIRI) observations (Beichman et al. 2014).At ∼ µ m, the two-coefficient laws may introduce aspectral slope of a few tens of ppm, which may havean impact in the analysis of the HST /WFC3 spec-tra (Tsiaras et al. 2018). At wavelengths shorter than1 µ m the two-coefficient laws are unreliable for exo-planet spectroscopy, so that the claret-4 law must bepreferred. These conclusions are in agreement with pre-vious studies based on both simulated and real data(Espinoza & Jord´an 2016; Morello et al. 2017; Morello2018; Maxted 2018).3.3. Predicted precision in light curves
Figure 7 shows that, for a fixed transit geometry, thepeak-to-peak of light curve residuals is roughly propor-tional to the weighted - r QS rms of stellar intensity resid-uals. We found an approximately linear correlation be-tween the two goodness-of-fit measures for the simulatedspectral light curves and stellar intensity profiles, there-fore obtaining a wavelength-independent factor. We re-peated this test for analogous sets of spectral light curveswith different transit parameters, then obtaining differ-ent proportionality factors. Our preliminary study sug-gests that(peak-to-peak) ppm = ( k × ) × p × ( weighted - r QS rms) , (20)where k is a factor of order unity ( k & USAGE OF
EXOTETHYS
Currently, the main use of the
ExoTETHyS package isto compute stellar limb-darkening coefficients throughthe SAIL subpackage. These coefficients can be adoptedto simulate the exoplanetary transit light curves, whichare largely used by the scientific consortia of the futureexoplanet missions for multiple studies. In particular,
ExoTETHyS will be linked with
ARIEL -Sim (Sarkar et al.2016) and
ExoNoodle (a generator of timeseries spectraof exoplanetary systems originally designed for
JWST observations; M. Martin-Lagarde et al., in prep.), andit has already been adopted by several members of thetwo mission consortia.It is also common practice to fix the limb-darkeningcoefficients obtained from stellar models, such as thosecalculated with
ExoTETHyS .SAIL, in the analysis of ex-oplanetary transit light curves. This approach relieson the perfect match between the model and the realstellar intensity distributions, otherwise introducing apotential bias in the derived exoplanet and orbital pa-rameters. Some authors recommended setting free limb-darkening coefficients in the light curve fits to minimizethe potential bias, but the strong parameter degenera-cies may lead to larger error bars or prevent the con-vergence of the fit (Southworth 2008; Csizmadia et al.2013). The parameter degeneracies can be mitigatedby using multiwavelength transit observations to bet-ter constrain the orbital parameters (Morello et al. 2017;Morello 2018). Here we suggest an approach to take ad-vantage of the knowledge on the stellar parameters inthe form of bayesian priors. The stellar parameters willthen be optimized in the light curve fits instead of usingfixed or fully unconstrained limb-darkening coefficients.The limb-darkening coefficients for a given set of stel-lar parameters, and a given passband or spectroscopicbin, can be interpolated from a precalculated grid. Thegrid calculation type (see Section 2.1.1) was specificallydesigned for this purpose. CONCLUSIONSWe introduced
ExoTETHyS , an open-source pythonpackage that offers accessory tools for modeling tran-siting exoplanets and eclipsing binaries. It includes aversatile stellar limb-darkening calculator with multiple4
Morello et al. choices of model atmosphere grids, parameterizations,passbands (also accepting user input), and specific user-friendly calculation settings. We demonstrated an opti-mal fitting algorithm for the limb-darkening coefficients,thus eliminating the degree of freedom associated withthe choice of fitting algorithm. The claret-4 coefficientsobtained through this algorithm ensure a precision level .
10 ppm in the relevant transit light curves at all wave-lengths. The precision achieved exceeds by one order ofmagnitude that obtained with most of the algorithmsproposed in the previous literature for stellar modelswith spherical geometry. We also proposed a simpleformula for estimating the light curve model precision,based on the goodness-of-fit for the limb-darkening co-efficients. Finally, we discussed the current and future usage of
ExoTETHyS with emphasis on exoplanet atmo-spheric spectroscopy in the era of
JWST and ARIEL.The authors would like to thank Ren´e Gastaud andDaniel Dicken for useful discussions. G. M., M. M.-L.and P.-O. L. were supported by the LabEx P2IO, theFrench ANR contract 05-BLANNT09-573739 and theEuropean Unions Horizon 2020 Research and InnovationProgramme, under Grant Agreement N ◦ Agol, E., Luger, R., & Foreman-Mackey, D. 2019, arXive-prints, arXiv:1908.03222.https://arxiv.org/abs/1908.03222Akinsanmi, B., Barros, S. C. C., Santos, N. C., et al. 2019,A&A, 621, A117, doi: 10.1051/0004-6361/201834215Ballerini, P., Micela, G., Lanza, A. F., & Pagano, I. 2012,A&A, 539, A140, doi: 10.1051/0004-6361/201117102Barnes, J. W., Linscott, E., & Shporer, A. 2011, ApJS, 197,10, doi: 10.1088/0067-0049/197/1/10Beichman, C., Benneke, B., Knutson, H., et al. 2014, PASP,126, 1134, doi: 10.1086/679566Chiavassa, A., Caldas, A., Selsis, F., et al. 2017, A&A, 597,A94, doi: 10.1051/0004-6361/201528018Christiansen, J. L., Jenkins, J. M., Caldwell, D. A., et al.2012, PASP, 124, 1279, doi: 10.1086/668847Claret, A. 2000, A&A, 363, 1081—. 2003, A&A, 401, 657, doi: 10.1051/0004-6361:20030142—. 2004, A&A, 428, 1001, doi: 10.1051/0004-6361:20041673—. 2008, A&A, 482, 259, doi: 10.1051/0004-6361:200809370—. 2017, A&A, 600, A30,doi: 10.1051/0004-6361/201629705—. 2018, A&A, 618, A20,doi: 10.1051/0004-6361/201833060Claret, A., & Bloemen, S. 2011, A&A, 529, A75,doi: 10.1051/0004-6361/201116451Claret, A., Cukanovaite, E., Burdge, K., et al. 2020, arXive-prints, arXiv:2001.07129.https://arxiv.org/abs/2001.07129Claret, A., Dragomir, D., & Matthews, J. M. 2014, A&A,567, A3, doi: 10.1051/0004-6361/201423515 Claret, A., Hauschildt, P. H., & Witte, S. 2012, A&A, 546,A14, doi: 10.1051/0004-6361/201219849—. 2013, A&A, 552, A16,doi: 10.1051/0004-6361/201220942Csizmadia, S., Pasternacki, T., Dreyer, C., et al. 2013,A&A, 549, A9, doi: 10.1051/0004-6361/201219888Diaz-Cordoves, J., & Gimenez, A. 1992, A&A, 259, 227Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83,doi: 10.1086/669497Espinoza, N., & Jord´an, A. 2015, MNRAS, 450, 1879,doi: 10.1093/mnras/stv744—. 2016, MNRAS, 457, 3573, doi: 10.1093/mnras/stw224Gazak, J. Z., Johnson, J. A., Tonry, J., et al. 2012,Advances in Astronomy, 2012, 697967,doi: 10.1155/2012/697967Gim´enez, A. 2006, A&A, 450, 1231,doi: 10.1051/0004-6361:20054445Hellard, H., Csizmadia, S., Padovan, S., et al. 2019, ApJ,878, 119, doi: 10.3847/1538-4357/ab2048Hestroffer, D. 1997, A&A, 327, 199Heyrovsk´y, D. 2007, ApJ, 656, 483, doi: 10.1086/509566Howarth, I. D. 2011a, MNRAS, 413, 1515,doi: 10.1111/j.1365-2966.2011.18122.x—. 2011b, MNRAS, 418, 1165,doi: 10.1111/j.1365-2966.2011.19568.xHowarth, I. D., & Morello, G. 2017, MNRAS, 470, 932,doi: 10.1093/mnras/stx1260Isaak, K. G., & Benz, W. 2019, Nature Astronomy, 3, 873,doi: 10.1038/s41550-019-0886-9 ample article15