The ExoTETHyS Package: Tools for Exoplanetary Transits around Host Stars

Section 2 provides a technical description of the ExoTETHyS package and the algorithms adopted. Section 3 discusses the precision of the limb-darkening calculator for the analysis of exoplanetary transits. In particular, Section 3.1 compares various algorithms that are adopted in the other publicly available codes and their variants, Section 3.2 compares the performances of the alternative limb-darkening laws, and Section 3.3 provides a formula to estimate the potential error in the transit model based on the goodness of fit for the limb-darkening coefficients that should be compared with the noise level in the observations. Section 4 discusses the main functionality of the ExoTETHyS package, its current and future usage. Finally, Section 5 summarizes the key points discussed in this paper.

The SAIL subpackage is a generic stellar limb-darkening calculator that is not specific to a predetermined list of instruments or standard passbands. It is conceptually similar to the calculator provided by Espinoza & Jordán (2015), but with different features. A technical difference is the use of a novel fitting algorithm for obtaining the limb-darkening coefficients, specifically optimized for modeling the exoplanetary transits, instead of multiple algorithm options with unclear performances (see Sections 2.1.4 and 3.1).

2.1.1. Input and Output

2.1.2. From the Stellar Model Atmospheres to the Passband-integrated Intensities

The stellar model-atmosphere grids consist of one file for each triple of stellar parameters (T_eff, log g, [M/H]), providing the specific intensities (I_λ(μ)) in units of erg cm⁻² s⁻¹ Å⁻¹ sr⁻¹ at several positions on the sky-projected stellar disk over a given spectral range. For historical reasons, the independent variable is $\mu =\cos \theta$, where θ is the angle between the line of sight and the corresponding surface normal. The radial coordinate in the sky-projected disk is $r=\sqrt{1-{\mu }^{2}}$, where r = 1 (μ = 0) corresponds to the spherical surface radius. Table 1 reports the information about the databases available with the first release of ExoTETHyS. We refer to the relevant papers and references therein for comparisons between the models. The passband-integrated intensities are calculated as

$$\begin{eqnarray}&&{I}_{\mathrm{pass}}(\mu )\propto {\int }_{{\lambda }_{1}}^{{\lambda }_{2}}{I}_{\lambda }(\mu ){R}_{\mathrm{pass}}(\lambda )\lambda d\lambda ,\end{eqnarray} \tag{ 1 }$$

where R_pass(λ) is the spectral response of the instrument in electrons photon⁻¹, and λ₁ and λ₂ are the passband or wavelength bin limits. The passband-integrated intensities are obtained in units proportional to electrons cm⁻² s⁻¹ sr⁻¹. As the limb-darkening coefficients are not affected by the (omitted) proportionality factor in Equation (1), the final intensities are normalized such that I_pass (μ = 0) = 1.

The intensity profiles, I_λ(μ), have distinctive behaviors depending on the plane-parallel or spherical geometry adopted by the selected grid of model atmospheres. In particular, the spherical intensity profiles show a steep drop-off close to the stellar limb, which is not observed in the plane-parallel models. The explanation for the different behaviors is exhaustive in the literature (Wittkowski et al. 2004; Espinoza & Jordán 2015; Morello et al. 2017). The almost null intensities at small μ are integrated over lines of sight that intersect only the outermost atmospheric shells, which have the smallest emissivity. Here μ = 0 (r = 1) corresponds to the outermost shell of the model atmosphere, which is typically outside the stellar radius that would be observed in transit. Our algorithm calculates the photometric radius at the inflection point of the spherical intensity profile, i.e., where the gradient $| {dI}(r)/{dr}|$ is the maximum (Wittkowski et al. 2004; Espinoza & Jordán 2015). The radial coordinates are then rescaled such that r = 1 (μ = 0) at the photometric radius, and those intensities with rescaled r > 1 are rejected. No rescaling is performed for the plane-parallel models.

2.1.3. Limb-darkening Laws

2.1.4. From the Passband-integrated Intensities to the Limb-darkening Coefficients

The limb-darkening coefficients are obtained through a weighted least-squares fit of the passband-integrated intensity profile with weights proportional to the sampling interval in r, hereinafter referred to as weighted-r fit. The corresponding cost function is the weighted rms of residuals,

$$\begin{eqnarray}&&{\rm{weighted \mbox{-} }}r\,\mathrm{rms}={\left(\displaystyle \frac{{\sum }_{i=1}^{n}{w}_{i}{\left({I}_{\mathrm{pass}}({\mu }_{i})-{I}_{\mathrm{pass}}^{\mathrm{law}}({\mu }_{i}\right)}^{2}}{{\sum }_{i=1}^{n}{w}_{i}}\right)}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 9 }$$

with weights

$$\begin{eqnarray}{w}_{i}=\left\{\begin{array}{ll}(1-{r}_{1})+0.5\,({r}_{1}-{r}_{2}), & {\rm{if}}\ i=1\\ 0.5\,({r}_{i-1}-{r}_{i+1}), & {\rm{if}}\ 1\lt i\lt n\\ 0.5\,{r}_{n-1}, & {\rm{if}}\ i=n\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where the r_i are arranged in descending order, and r_n = 0. The choice of cost function is optimized for the study of exoplanet transits, as detailed in Section 3.1. The performances of the spherical model fits are further enhanced by discarding those points with r > 0.99623 (after rescaling as explained in Section 2.1.2). This cut is a generalization of that implemented in the quasi-spherical (QS) fits by Claret et al. (2012). For this reason, we rename the total fitting procedure explained here for the spherical intensity profiles as the weighted-r QS fit. Further details about the alternative fitting procedures are discussed in Section 3.1.

2.1.5. Interpolation from the Grid of Stellar Models

The TRIP subpackage is used to generate exact transit light curves by direct integration of the occulted stellar flux at given instants. It assumes a dark spherical planet transiting in front of a spherically symmetric (unspotted) star. In this simple case, the normalized flux (i.e., relative to the stellar flux) is a function of two geometric variables, as reported by Mandel & Agol (2002), and of the stellar intensity profile:

$$\begin{eqnarray}&&F(p,z,I(\mu ))=1-{\rm{\Lambda }}(p,z,I(\mu )),\end{eqnarray} \tag{ 11 }$$

where p is the planet-to-star radii ratio (p = R_p/R_*), z is the sky-projected distance between the centers of the two spheres in units of the stellar radius, and I(μ) is the stellar intensity profile. TRIP does not use an analytical approximation of the limb-darkening profile, unlike most transit light-curve calculators such as those provided by Mandel & Agol (2002), Giménez (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

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 closest transit event:

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \frac{t-{\rm{E}}.{\rm{T}}.}{P}-n,\end{eqnarray} \tag{ 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 orbits from the E.T. rounded to the nearest integer. Conventionally, Φ values are in the range of [−0.5, 0.5] or [0, 1] and Φ = 0 at mid-transit time. The projected star–planet separation is given by

$$\begin{eqnarray}z=\left\{\begin{array}{ll}{a}_{R}\sqrt{1-{\cos }^{2}(2\pi {\rm{\Phi }}){\sin }^{2}i} & {\rm{circular}}\,{\rm{orbit}}\\ {a}_{R}\displaystyle \frac{1-{e}^{2}}{1+e\cos f}\sqrt{1-{\sin }^{2}(f+\omega ){\sin }^{2}i} & \ \ \ {\rm{eccentric}}\,{\rm{orbit}}\end{array}\right.,\end{eqnarray} \tag{ 13 }$$

where a_R is the orbital semimajor axis in units of the stellar 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 from the orbital phase by solving Kepler's equation,

$$\begin{eqnarray}&&\displaystyle \frac{\pi }{2}-\omega +2\pi {\rm{\Phi }}=E-e\sin E,\end{eqnarray} \tag{ 14 }$$

then

$$\begin{eqnarray}&&f=2\arctan \left(\sqrt{\displaystyle \frac{1+e}{1-e}}\tan \displaystyle \frac{E}{2}\right).\end{eqnarray} \tag{ 15 }$$

2.2.3. Calculating the Normalized Flux

The total and occulted stellar flux are given, respectively, by the integrals

$$\begin{eqnarray}&&{F}_{* }={\int }_{0}^{1}I(r)\,2\pi r\,{dr},\end{eqnarray} \tag{ 16 }$$

and

$$\begin{eqnarray}&&{F}_{* ,\mathrm{occ}}={\int }_{0}^{1}I(r)\,2\pi r\,{f}_{p,z}(r)\,{dr},\end{eqnarray} \tag{ 17 }$$

with

$$\begin{eqnarray}&&\begin{array}{l}{f}_{p,z}(r)\ =\\ {\left.\left\{\begin{array}{ll}\displaystyle \frac{1}{\pi }\arccos \displaystyle \frac{{r}^{2}+{z}^{2}-{p}^{2}}{2{zr}} & | z-p| \lt r\lt z+p\\ 0 & r\leqslant z-p{\rm{or}}r\geqslant z+p\\ 1 & r\leqslant p-z\end{array}\right.\right|}_{0\leqslant r\leqslant 1}.\end{array}\end{eqnarray} \tag{ 18 }$$

I(r) is the specific intensity at the normalized radial coordinate $r=\sqrt{1-{\mu }^{2}}$, and f_{p, z}(r) is the fraction of circumference with radius r covered by the planet. Equations (16) and (17) rely on the assumed spherical symmetry for the star; Equation (18) also makes use of the planet sphericity. Finally, the normalized flux is given by Equation (11) with

$$\begin{eqnarray}&&{\rm{\Lambda }}(p,z,I(\mu ))=\displaystyle \frac{{F}_{* ,\mathrm{occ}}}{{F}_{* }}.\end{eqnarray} \tag{ 19 }$$

The integrals in Equations (16) and (17) are calculated numerically using the mid-point rule with a uniform partition in r. The specific intensities are evaluated at the partition radii by interpolating in μ or r from the input limb-darkening profiles. The TRIP algorithm with default settings is identical to the "tlc" described by Morello et al. (2017).

A long list of methods has been adopted in the literature for fitting the limb-darkening laws to the model intensity profiles leading to significantly different limb-darkening coefficients. The coefficients obtained with a simple least-squares fit depend on the spatial distribution of the precalculated intensities. The effect of sampling is particularly evident for the PHOENIXprofiles because of a much finer sampling near the drop-off region. For example, Figure 1 shows the case of a star similar to HD209458 in the mid-infrared, for which the simple least-squares solution presents a non-monotonic (unphysical) profile with unexpected undulations. In this paper, we compare the following fitting procedures:

• 1.  
  unweighted, i.e., simple least-squares fit;
• 2.  
  weighted-r, i.e., weighted least-squares fit with weights proportional to the sampling interval in r, as detailed in Equations (9) and (10);
• 3.  
  weighted-μ, i.e., weighted least-squares fit with weights proportional to the sampling interval in μ;
• 4.  
  interp-μ 100, i.e., least-squares fit on the intensities interpolated over 100 μ values with a uniform separation in μ, as suggested by Claret & Bloemen (2011);
• 5.  
  interp-μ 1000, i.e., least-squares fit on the intensities interpolated over 1000 μ values with a uniform separation in μ;
• 6.  
  interp-r 100, i.e., least-squares fit on the intensities interpolated over 100 r values with a uniform separation in r, as suggested by Parviainen & Aigrain (2015; with an unspecified number of interpolated values);
• 7.  
  interp-r 1000, i.e., least-squares fit on the intensities interpolated over 1000 r values with a uniform separation in r;
• 8.  
  unweighted QS, i.e., least-squares fit with a cutoff r ≤ 0.99623;
• 9.  
  weighted-r QS, i.e., analogous to weighted-r with a cutoff r ≤ 0.99623;
• 10.  
  weighted-μ QS, i.e., analogous to weighted-μ with a cutoff r ≤ 0.99623.

The cutoff is used to remove the steep drop-off characteristic of the spherical models, hence the term QS. The QS approach was first proposed by Claret et al. (2012), who applied a cutoff μ ≥ 0.1 to their library of PHOENIX models with the original μ values. In this work, we redefine the cutoff using the rescaled r, such that it corresponds to the same fraction of the photometric stellar radius for all the models (see Section 2.1.2). Our new definition with r ≤ 0.99623 is equivalent to the previous one for the majority of models, particularly for those models that may correspond to main-sequence stars. However, the libraries of PHOENIX models incorporated in the ExoTETHyS package also include models of stellar atmospheres with lower gravities than those analyzed by Claret et al. (2012), corresponding to subgiant, giant, and supergiant stars. For some of these models, the intensity drop-off occurs at μ > 0.1, so that the cutoff of μ ≥ 0.1 (not rescaled) would be ineffective.

Zoom In Zoom Out Reset image size

In order to evaluate the merits of the alternative fitting procedures to the stellar intensity profile, we generated exact synthetic transit light curves using the TRIP subpackage and compared these light curves with their best-fit solutions obtained with the various sets of claret-4 limb-darkening coefficients. Figure 2 shows the residuals obtained for a noiseless simulation of the transit of HD209458 b in the TESS passband when adopting the different sets of limb-darkening coefficients. The weighted-r QS method implemented in ExoTETHyS.SAIL gives the smallest residuals, with a peak-to-peak of 2 ppm and rms amplitude below 1 ppm. The other QS methods, weighted-μ QS and unweighted QS, lead to almost identical residuals, with a peak-to-peak of 3 ppm. Among the spherical methods, the weighted-r gives the smallest residuals with a peak-to-peak of 9 ppm and an rms amplitude of 2 ppm, followed by the interp-r 100 and interp-r 1000 with about 1.5 and 2 times larger residual amplitudes, respectively. All the other methods lead to significantly larger residuals of tens to a few hundred ppm, which are comparable with the predicted noise floor of 60 ppm for the TESS observations (Ricker et al. 2014).

Zoom In Zoom Out Reset image size

Figure 3 shows the peak-to-peak of the residuals for the same transit as a function of wavelength, based on simulated light curves with 20 nm passband widths. This spectral analysis confirms the relative ranking of the 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 longer than 1 μm, and overall below 8 ppm. The other QS methods are marginally worse than weighted-r QS at wavelengths shorter than 2 μm, but the worst case peak-to-peak of residuals is less than 13 ppm. The weighted-r method leads to a peak-to-peak of residuals in the range of 5–15 ppm, with a sawtooth-like modulation in wavelength. We noted that the small, but abrupt, jumps that occur at certain wavelengths correspond to changes of the inflection point in the stellar intensity profile as defined in Section 2.1.2. The same phenomenon occurs for all the other spherical models with larger sawtooth-like modulations. It may appear surprising that the peak-to-peak of residuals obtained with the spherical methods tends to be larger at the longer wavelengths, for which the limb-darkening effect is expected to be smaller. The cause of the poor performances of most spherical methods in the infrared is the intensity drop-off, which is typically steeper than the drop-off in the UV and visible. Such drop-off has a negligible effect in the numerically integrated transit light curves, hence the better performances of the QS fits.

Zoom In Zoom Out Reset image size

Figure 4 shows the best-fit transit parameters corresponding to the same spectral light curves, and compared with the respective input parameters corrected for the rescaled r (see Section 2.1.2). We retrieved the correct transit depth within 5 ppm, the impact parameter within 6 × 10⁻⁴, and the transit duration within 1 s at all wavelengths, when using the weighted-r or QS limb-darkening coefficients. However, slightly larger spectral trends appear in these parameters because of the wavelength-dependent stellar radius. The peak-to-peak variation in transit depth over the spectral range of 0.25–10 μm is 10 ppm. The other sets of limb-darkening coefficients introduce orders-of-magnitude larger biases in the retrieved transit parameters, also larger spectral sawtooth-like modulations in the infrared (few tens of ppm in transit depth across 1–10 μm), and severe discrepancies between the parameter values obtained in the UV/visible and those obtained in the infrared.

Zoom In Zoom Out Reset image size

Figure 5 compares the peak-to-peak of the spectral light-curve residuals when adopting the limb-darkening coefficients calculated by ExoTETHyS.SAIL for different limb-darkening laws, as well as the corresponding weighted-r QS rms of the residuals to the stellar intensity profiles. The correlation between the two goodness-of-fit measures is explored in Section 3.3. At wavelengths ≳3 μm, the precision of the power-2 and square-root limb-darkening coefficients is comparable to that of the claret-4 coefficients, resulting in light-curve residuals below 5 ppm. While the claret-4 law performs similarly well even at shorter wavelengths, the two-coefficient laws lead to larger light-curve residuals up to ∼100 ppm in the UV and visible. The quadratic law is less precise, leading to light-curve residuals above 25 ppm even at 10 μm.

Zoom In Zoom Out Reset image size

Figure 6 shows the fitted transit parameters and their expected values. Typically, the bias in transit depth is of the same order of magnitude of the light-curve residuals, but it can be both larger or smaller than their peak-to-peak amplitudes owing to parameter degeneracies. Table 2 reports the statistics of the errors in transit depth obtained with the different limb-darkening laws across given spectral ranges. The maximum bias in transit depth at 5–10 μm is within 10 ppm for any limb-darkening parameterization, which is just below the minimum photon noise floor for JWST/Mid-InfraRed Instrument (MIRI) observations (Beichman et al. 2014). At ∼1 μm, the two-coefficient laws may introduce a spectral slope of a few tens of ppm, which may have an impact in the analysis of the HST/WFC3 spectra (Tsiaras et al. 2018). At wavelengths shorter than 1 μm the two-coefficient laws are unreliable for exoplanet spectroscopy, so that the claret-4 law must be preferred. These conclusions are in agreement with previous studies based on both simulated and real data (Espinoza & Jordán 2016; Morello et al. 2017; Maxted 2018; Morello 2018).

Zoom In Zoom Out Reset image size

Figure 7 shows that, for a fixed transit geometry, the peak-to-peak of light-curve residuals is roughly proportional to the weighted-r QS rms of stellar intensity residuals. We found an approximately linear correlation between the two goodness-of-fit measures for the simulated spectral light curves and stellar intensity profiles, therefore obtaining a wavelength-independent factor. We repeated this test for analogous sets of spectral light curves with different transit parameters, then obtaining different proportionality factors. Our preliminary study suggests that

$$\begin{eqnarray}&&\begin{array}{l}{\left({\rm{peak \mbox{-} to \mbox{-} peak}}\right)}_{\mathrm{ppm}}=(k\times {10}^{6})\times {p}^{2}\\ \qquad \times \ ({\rm{weighted \mbox{-} }}r\,{\rm{QS}}\,{\rm{rms}}),\end{array}\end{eqnarray} \tag{ 20 }$$

where k is a factor of order unity (k ≳ 1). Equation (20) provides a useful tool for estimating the systematic noise in the light-curve models solely due to the limb-darkening parameterization. The systematic noise in the light-curve models should be smaller than the photon noise limit of the observation in order to avoid significant parameter biases. Note that Equation (20) does not account for uncertainties in the stellar parameters, discrepancies between real and model intensity profiles, and other contaminating signals that may increase the total systematic noise.

Zoom In Zoom Out Reset image size