# Empirical limb-darkening coefficients & transit parameters of known
exoplanets from TESS

Jayshil A. Patel Department of Astronomy, Stockholm University, AlbaNova
University Center, 10691 Stockholm, Sweden  Néstor Espinoza Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Department
of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

###### Abstract

Although the main goal of the Transiting Exoplanet Survey Satellite (TESS) is
to search for new transiting exoplanets, its data can also be used to study in
further detail already known systems. The TESS bandpass is particularly
interesting to study the limb-darkening effect of the stellar host which is
imprinted in transit lightcurves, as the widely used phoenix and atlas stellar
models predict different limb-darkening profiles. Here we study this effect by
fitting the transit lightcurves of 176 known exoplanetary systems observed by
TESS, which allows us to extract empirical limb-darkening coefficients (LDCs)
for the widely used quadratic law, but also updated transit parameters
(including ephemerides refinements) as a byproduct. Comparing our empirically
obtained LDCs with theoretical predictions, we find significant offsets when
using tabulated TESS LDCs. Specifically, the $u_{2}$ coefficients obtained
using phoenix models show the largest discrepancies depending on the method
used to derive them, with offsets that can reach up to $\Delta u_{2}\approx
0.2$ on average. Most of those average offsets disappear, however, if one uses
the SPAM algorithm introduced by Howarth (2011) to calculate the LDCs instead.
Our results suggest, however, that for stars cooler than about 5000 K, no
methodology is good enough to explain the limb-darkening effect: we observe a
sharp deviation between measured and predicted LDCs on both quadratic LDCs of
order $\Delta u_{1},\Delta u_{2}\approx 0.2$ for those cool stars. We
recommend caution when assuming limb-darkening coefficients as perfectly known
thus, in particular for these cooler stars when analyzing TESS transit
lightcurves.

editorials, notices — miscellaneous — catalogs — surveys

††journal: AJ

##  1 Introduction

The phenomenon of limb-darkening on stars — the observed intensity decrease
towards their limb — has for long been recognized as a natural consequence of
temperature gradients in stellar atmospheres (Schwarzschild, 1906). Although a
complex effect in nature, parametric formulations of it have been shown to be
very practical to use when the effect needs to be modelled. These so-called
“limb-darkening laws” have been mainly motivated from an observational
perspective and have been formulated with certain use-cases in mind (see, for
example Claret, 2000; Diaz-Cordoves & Gimenez, 1992; Sing et al., 2009;
Klinglesmith & Sobieski, 1970; Claret & Hauschildt, 2003). In exoplanetary
science, in particular, the most popular laws to date are the linear law
(Schwarzschild, 1906), the quadratic law (Kopal, 1950) and the non-linear law
(Claret, 2000). Their functional forms are,

$\displaystyle\frac{I(\mu)}{I(1)}$ $\displaystyle=$ $\displaystyle\
1-a(1-\mu);\ \textnormal{(the linear law),}$
$\displaystyle\frac{I(\mu)}{I(1)}$ $\displaystyle=$ $\displaystyle\
1-u_{1}(1-\mu)-u_{2}{(1-\mu)}^{2};\ \textnormal{(the quadratic law),}$
$\displaystyle\frac{I(\mu)}{I(1)}$ $\displaystyle=$ $\displaystyle\
1-\sum_{n=1}^{4}c_{n}(1-\mu^{n/2});\ \textnormal{(the non-linear law)}.$

Here $a$, $(u_{1},u_{2})$ and $(c_{1},c_{2},c_{3},c_{4})$ are the so-called
limb darkening coefficients (LDCs) of these respective laws and $\mu$
represents the cosine of the angle $\theta$ \- the angle between the line of
sight and the normal to the surface. The reason for the popularity of these
laws is most likely related with the fact that these give rise to very fast
lightcurve computations (Mandel & Agol, 2002), which in turn makes them the
most practical to use. Among them, the quadratic law is by far the most widely
used in the exoplanet literature mainly due to the fact that it gives rise to
fast and efficient computation of lightcurves in some of the most popular
lightcurve computing methods/algorithms (e.g., batman; Kreidberg, 2015). The
simplicity of this law also makes it an attractive one to use for researchers:
it is defined by only two parameters which can even be fitted along with the
rest of the other transit parameters if stellar model atmospheres are
suspected to not be appropriate for a given use-case. While it might not
always be the most accurate law to use (see, e.g., Southworth, 2008; Hayek et
al., 2012; Espinoza & Jordán, 2016), its usage is nonetheless widespread
across transiting exoplanet studies because of the aforementioned reasons.

The details of the limb-darkening effect become important when one wants to
fit an exoplanet transit lightcurve because its shape is strongly modulated by
the effect (Mandel & Agol, 2002; Seager & Mallén-Ornelas, 2003). Incorrect
assumptions about the effect, thus, have the potential to lead to biases on
the retrieved transit parameters. Indeed, simulations have studied in detail
the impact of either fitting or fixing the LDCs (see, e.g., Csizmadia et al.,
2013; Espinoza & Jordán, 2015; Neilson et al., 2017; Morello et al., 2017) and
even the impact the selection of a given limb-darkening law has on the
retrieval of the transit parameters (Espinoza & Jordán, 2016). Most of these
studies conclude that the most conservative assumption is to fit for the
coefficients in the transit lightcurve fitting procedure, taking special care
on the selection of the limb-darkening law to use for very precise (better
than about 200 ppm) transit lightcurves. This suggestion, however, is
impractical in cases on which the lightcurve is poorly sampled and/or
dominated by systematic noise such as in, e.g., Hubble Space Telescope (HST)
spectro-photometric measurements (see, e.g., Kreidberg, 2018). It is also not
clear from the literature how far off empirical LDCs are from stellar model
atmospheres predictions, which of all the stellar models available in the
literature are the best to use, and which exact method/table one should use if
one needed a set of LDCs to either fix or use as priors in the analysis (see,
e.g., Hayek et al., 2012; Espinoza & Jordán, 2015; Sandford & Kipping, 2017).
Studies such as that of, e.g., Hayek et al. (2012), Müller et al. (2013a),
Espinoza & Jordán (2015) and Maxted (2018), have approached the problem by
comparing LDCs empirically determined from transit lightcurves to that of
predictions from stellar model atmospheres. This intercompraison between
empirical LDCs and those obtained from actual data are critical not only to
inform stellar modellers, but also the community on best practices and uses of
stellar model atmospheres for exoplanet transit lightcurves. Such a comparison
has the power to aid us in the search for better models, and empirical LDCs
could provide good priors especially useful when confronted with low precision
lightcurves.

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al., 2014) mission
does not only provide exquisite, almost all-sky photometry with which to look
for new transiting exoplanets, but also provides a unique dataset of already
known exoplanets which is ideal to perform tests on stellar model atmospheres
through stellar limb-darkening. Given that most known transiting exoplanetary
systems have been thoroughly characterized with ground-based instruments,
precise spectroscopic parameters are available for their host stars and thus
limb-darkening predictions can be extracted from stellar model atmospheres.
The precise photometry being collected by the TESS mission, in turn, should
allow us to tightly constrain the limb-darkening effect through the LDCs, thus
providing a rich dataset to compare against theoretical predictions.
Interestingly (and as it will be shown in this work), LDCs from the most
widely used stellar model atmospheres to model limb-darkening in exoplanet
transit lightcurves, namely, the atlas (Castelli & Kurucz, 2003) and phoenix
(Husser et al., 2013) model atmospheres, have very different values in the
TESS bandpass depending on which methods one uses to obtain them, which
provides a perfect opportunity to test these two stellar model atmospheres
(and the methods used to obtain LDCs) against empirical data. In this work, we
use recently released TESS transit lightcurve data in order to extract limb
darkening coefficients for a subset of the known systems observed by the
mission in order to perform these tests.

We have organized this work as follows. In Section 2, we describe how we
select our targets and download the data for them, along with a description of
the modelling procedure used to analyze the lightcurves and to obtain
theoretical predictions for our LDCs. We present our findings in Section 3,
which are followed by a Discussion in Section 4 and Conclusions and Future
Work in Section 5.

##  2 Data, sample selection and modelling

###  2.1 Selection of targets

Because the main objective of the present work was to characterize the stellar
hosts, we chose TESS targets that had already known and confirmed transiting
exoplanets, and that had precise follow-up observations that allowed for
precise stellar (e.g., stellar effective temperatures, gravity, metallicity)
and planetary (e.g., via radial-velocity variations) characterization using
high-resolution spectroscopy. In addition, only exoplanets observed in
2-minute cadence by TESS were selected in order to minimize lightcurve
distortions due to the binning of the 30-minute cadence lightcurves (Kipping,
2010).

Figure 1: Range and distribution of Effective Temperature ($T_{eff}$) and
Surface Gravity ($\log g$) of target exoplanet host stars.

According to the NASA Exoplanet
Archive111https://exoplanetarchive.ipac.caltech.edu/cgi-bin/TblView/nph-
tblView?app=ExoTbls&config=planets, when we started our analysis, more than
4000 exoplanets were known. Among them, we chose transiting exoplanets which
were going to be or already were observed by TESS. We filtered these planets
using the Web TESS Viewing tool (WTV)222https://heasarc.gsfc.nasa.gov/cgi-
bin/tess/webtess/wtv.py. At the time of our analysis, the TESS mission had
observed targets up to sector 32, so we selected targets only up to this
sector. There were a total of 1745 of them. Among these systems, we removed
from our sample:

  * •

Systems that did not have follow-up high-resolution spectroscopic observations
that allowed us to confirm their orbits in an independent way to the transits.

  * •

Multi-planetary systems, which would involve extra complications in the
analyses (TTVs, joint fitting of some properties, checks for planet-planet
transits).

  * •

Systems that had known physically associated stellar companions, as they could
significantly dilute the exoplanet transit lightcurves.

When we analysed these filtered systems, we found that some of them had very
low signal-to-noise ratio transit lightcurves. We decided to measure the
signal-to-noise ratio of the transits through their transit depths, and
decided that systems for which the signal-to-noise ratio of the transit depths
were less than 5 would also be removed from our study. After this, a total of
176 targets remained in our sample. Figure 1 shows the range and the
distribution of effective temperature ($T_{eff}$) and surface gravity
($\log(g)$ \- in cgs units) of the exoplanet host stars of our target systems.
Table 3 shows the values of various stellar properties of each of our host
target stars.

###  2.2 Data and Modelling

We used juliet (Espinoza et al., 2019a) to download and fit the 2-minute
cadence lightcurve data directly from the Mikulski Archive for Space
Telescopes
(MAST333https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html) portal
in an automated fashion. Times, fluxes and errors were extracted from the PDC-
corrected lightcurve products.

The transit lightcurves were modelled within juliet using the batman
(Kreidberg, 2015) model. In addition to it, each lightcurve fit also included
a Gaussian Process (GP) in order to marginalize either astrophysical or
instrumental systematic trends present in the TESS lightcurves. In particular,
for each lightcurve, we fitted a GP using two kernels: an Exponential-Matèrn
kernel and a Quasi-periodic kernel. The first one is a result of multiplying
an exponential kernel, which typically samples smooth functions, and a Matèrn
3/2 kernel, which allows to model abrupt changes in the lightcurves, making
the resultant multiplied kernel a very flexible one to model instrumental
systematics or non-periodic phenomena. The Quasi-periodic kernel we opted to
use, on the other hand, was the one introduced in Foreman-Mackey et al.
(2017), and of the form,

$k(x_{l},x_{m})=\frac{B}{2+C}e^{-\tau/L}\left[\cos\left(\frac{2\pi\tau}{P_{\textnormal{rot}}}\right)+(1+C)\right],$
(1)

where $B$, $C$, $L$ and $P_{\textnormal{rot}}$ are the hyperparameters of the
kernel. This is a very useful kernel when it comes to modelling quasi-periodic
phenomena such as the one observed in, e.g., starspot rotational modulations.
These two kernels are implemented within juliet using celerite (Foreman-Mackey
et al., 2017).

In practice, our fits were done in a two-step fashion for each target. First,
the out-of-transit data for each sector was analyzed separately, and fitted
with a GP, plus a jitter term and a flux normalization factor. We then
compared the bayesian evidence on each sector for each of the kernels outlined
above, and the one that was preferred for most of the sectors was defined as
the kernel to use for that target. Then, the posteriors on each of these
parameters for each sector were used as priors for a multi-sector joint fit of
the in-transit data, to which we added the transit model. This two-step
process for fitting these transit lightcurves had two main advantages. The
first advantage is that this provides a much more efficient sampling of the
parameter space given we use nested sampling algorithms to explore it (which
are much better behaved in low dimensional spaces). The second advantage,
which is a byproduct of the first one, is that it allows for a much faster
convergence of the algorithms. We follow the recommendation in Espinoza et al.
(2019a), and use importance nested sampling with the MultiNest algorithm
(Feroz et al., 2009) via the PyMultinest library (Buchner et al., 2014) to
perform fits with up to about 20 free parameters, and use dynamic nested
sampling via the dynesty library (Speagle, 2020) to explore higher dimensional
parameter spaces. In practice, this means that we use MultiNest for multi-
sector fits involving less than 4 sectors, and dynesty for multi-sector fits
involving 4 or more sectors. Our automated routines to perform this two-step
fitting procedure is available at Github444https://github.com/nespinoza/tess-
limb-darkening/.

Table 1: Summary of priors used in our analysis. $\mathcal{U}(a,b)$ and $\mathcal{J}(a,b)$ show the uniform and log-uniform distribution between $a$ and $b$, respectively; $\mathcal{N}(\mu,\sigma^{2})$ gives the normal distribution with mean $\mu$ and the standard deviation $\sigma$. The “$lit$” subscript to symbols means that the value of that parameter is taken from the literature. A single value indicates that the parameter was fixed. Parameter | Unit | Prior
---|---|---
Planetary Parameters
$P$ | days | $\mathcal{N}(\mu_{lit},\sigma^{2}_{lit})$
$T_{0}$ | BJD | $\mathcal{N}(\mu_{lit},0.1)$
$R_{p}/R_{*}$ |  | $\mathcal{U}(0.,1.)$
$b$ |  | $\mathcal{U}(0.,1.)$
$a/R_{*}$ |  | $\mathcal{J}(1.,100.)$
eccentricity |  | $ecc_{lit}$
$\omega$ | deg | $\omega_{lit}$
Instrumental Parameters
$q_{1}$ (transformed LDC) |  | $\mathcal{U}(0.,1.)$
$q_{2}$ (transformed LDC) |  | $\mathcal{U}(0.,1.)$
$m_{flux}$ | ppm | $\mathcal{N}(0.,10^{5})$
$m_{dilution}$ |  | 1
$\sigma_{w}$ (jitter) | ppm | $\mathcal{J}(0.1,10000)$
Exponential-Matèrn kernel
Amplitude of the GP | ppm | $\mathcal{J}(10^{-5},10000)$
Timescale (Exp part) | days | $\mathcal{J}(0.001,100)$
Timescale (Matèrn part) | days | $\mathcal{J}(0.001,100)$
Quasi-Periodic kernel
B |  | $\mathcal{J}(10^{-5},10^{3})$
C |  | $\mathcal{J}(10^{-5},10^{4})$
L | days | $\mathcal{J}(10^{-3},10^{3})$
$P_{rot}$ | days | $\mathcal{J}(1.,100.)$

The priors used in our fits were, in general, wide enough to explore a vast
part of the parameter space. The only exception was the period of the orbits,
whose priors were normal distributions with means and standard deviations
taken from previous works in the literature — our posteriors on this
parameter, then, can be seen effectively as an update over those previous
estimates. The time-of-transit center also had a more or less constrained
prior; we also defined a normal distribution as a prior on this parameter
based on previous works, but the standard deviation of this distribution was
set to 0.1 days (i.e., about 2.4 hours). To parametrize the limb-darkening we
decided to use the widely used quadratic law via the uninformative sampling
scheme proposed by Kipping (2013). In this setup, instead of fitting for the
limb-darkening coefficients directly, we fit for the transformed parameters
$q_{1}$ and $q_{2}$, each of which we define to have a uniform prior
distribution between 0 and 1. We defined wide priors between 0 and 1 for the
planet-to-star radius ratio and the impact parameter of the orbits. For
$a/R_{*}$, the scaled semi-major axis, we decided to use a wide log-uniform
prior between 1 and 100. In order to account for possible out-of-transit
offsets, we also fitted for a mean out-of-transit flux in our transit fits,
$m_{\textnormal{flux}}$, which normalizes the lightcurve via
$1/(1+m_{\textnormal{flux}})$ (see Espinoza et al., 2019a, for details); we
defined a normal prior for this parameter with zero-mean and a standard
deviation of $10^{5}$ ppm. A dilution factor was set to 1 for all of our fits
(i.e., assuming no dilution by nearby sources), as most of our targets are
bright and the PDC algorithm is supposed to correct for the dilution of nearby
objects anyways. Finally, for our GP kernels, we also defined wide priors. For
the Exponential-Matèrn kernel, we used a log-uniform prior on the amplitude of
the GP from $10^{-5}$ to 10,000 ppm, and time-scales for both the exponential
and Matern parts of the kernel between 0.001 and 100 days. For the quasi-
periodic kernel, the B, C and L parameters had log-uniform priors between
$10^{-5}$ to $10^{3}$, $10^{-5}$ to $10^{4}$, and $10^{-3}$ to $10^{3}$. The
rotation period of the kernel had a log-uniform prior between 1 and 100 days.
All lightcurve models also included a white-gaussian zero-mean noise
component, whose standard deviation was also fitted in our procedure. The
prior for that “jitter” parameter was defined to be log-uniform between 0.1
and 10,000 ppm. The eccentricity and argument of periastron on our fits was
fixed to literature values.

We also note here that although our target list consist of already known
exoplanets observed by TESS up to Sector 32, we downloaded all available data
for these targets, which may include, for some of the systems, data from the
most recent data release of Sector 34.

###  2.3 Theoretical calculation of Limb Darkening Coefficients

One of the key parts of our work involves computing theoretical LDCs in order
to compare them with our empirical estimates obtained from TESS lightcurve
fits. In this work, we used two main methods to compute these. The first, and
the most popular method used by the community, is to use various tables/codes
published in the literature to directly retrieve those coefficients. These
coefficients, in turn, are obtained by fitting the intensity profiles of model
stellar atmospheres, weighted by the instrumental bandpasses (see, e.g.,
Claret et al., 2012; Espinoza & Jordán, 2015). The second method we use is the
one suggested by Howarth (2011): the Synthetic-Photometry/Atmosphere-Modelling
(or, SPAM) technique. This is a technique explicitly designed to work with
exoplanetary transit lightcurves, and involves generating synthetic transit
lightcurves with limb-darkening coefficients estimated from the same
tables/codes as the ones described above, which are then fitted with a limb-
darkening law of choice to retrieve the LDCs. We detail how we explicitly
obtain the LDCs we use in this work with these two methods below.

####  2.3.1 LDCs from Tables/Codes

In order to perform theoretical predictions for LDCs obtained from our
lightcurve fits, we made use of limb-darkening tables and/or algorithms
present in the literature that use both atlas and phoenix model stellar
atmospheres as inputs. To put all those predictions on equal footing and to
maximize the amount of information extracted from those tables, instead of
extracting tabulated quadratic limb-darkening coefficients we followed the
approach of Espinoza & Jordán (2015), in which “limiting” limb darkening
coefficients are derived using the non-linear law. This approach guards
against the fact that different models and methods to fit for the intensity
profiles of model stellar atmospheres give rise to slightly different results
depending on the number of points sampled from those profiles. In essence, the
method simply performs $\chi^{2}$ \- minimization between the non-linear law
and the quadratic law, to find a relationship that gives the quadratic limb-
darkening coefficients as a function of the non-linear limb-darkening
coefficients (see Espinoza & Jordán, 2015, for a derivation), which is given
by:

$\begin{split}u_{1}&=\frac{12}{35}c_{1}+c_{2}+\frac{164}{105}c_{3}+2c_{4}\\\
u_{2}&=\frac{10}{21}c_{1}-\frac{34}{63}c_{3}-c_{4}.\end{split}$ (2)

Here $c_{1}$, $c_{2}$, $c_{3}$ and $c_{4}$ are the LDCs of the non-linear law,
and $u_{1}$ and $u_{2}$ are the resultant “limiting” LDCs for the quadratic
law. We use these latter ones to compare our retrieved limb-darkening
coefficients from our transit fits.

We retrieve theoretical LDC calculations on the TESS bandpass for our targets
using two different sources/methods. The first set is obtained using the limb-
darkening555https://github.com/nespinoza/limb-darkening code outlined in
Espinoza & Jordán (2015). The second are the tables published by Claret
(2017), which are arguably very popular among researchers that analyze TESS
data.

There are several details and/or assumptions one has to be careful with when
dealing with those sources of LDCs, which we here pay attention to when
comparing them against our empirically determined LDCs using TESS lightcurves.
The first is that, as put forward by Wittkowski et al. (2004), atlas and
phoenix model stellar atmospheres disagree about where the stellar atmosphere
ends, which in turn implies that different results between both are
practically guaranteed if not accounted for. Different authors have approached
this problem in different ways. On the one hand, Espinoza & Jordán (2015)
follow Wittkowski et al. (2004) and look for an inflection point in the
intensity profile of spherically symmetric models (i.e., where the derivative
of the intensity as a function of $r=\sqrt{1-\mu^{2}}$ reaches its maximum),
redefining that point to be at $r=1$ (i.e., $\mu=0$). Claret (2017), on the
other hand, suggests that a “better” approach — where “better” is defined as
attaining a lower residual sum of squares with respect to the fitted intensity
profile (see their Figure 1) — is to simply not use values where $\mu<0.1$, a
method they refer to as the “quasi-spherical” method. Those two methods are
fundamentally different approaches, and it is furthermore unclear which one of
those is better (or if they can be even distinguished with actual data). In
order to attempt at testing their efficacy at predicting empirically
determined LDCs, here we extract LDCs from the tables of Claret (2017) using
both methods. Following that same study, in this work we refer to the former
as the r-method, and to the latter as the q-method.

The second important detail has to do with the fact that Espinoza & Jordán
(2015) used the “vanilla” PHOENIX - COND models, available from the phoenix
website666ftp://phoenix.astro.physik.uni-goettingen.de/SpecIntFITS/. Claret
(2017), however, suggested that their work uses slightly different versions of
those models (the ones from Claret et al., 2012), and that furthermore,
discrepancies between the works of Claret et al. (2012) and Espinoza & Jordán
(2015) for stars with effective temperatures lower than 3000 K could be
explained due to Claret et al. (2012) using PHOENIX-DRIFT models. While the
latter cannot be tested with the data obtained in this work (as none of our
stars has effective temperatures smaller than 3000 K), the impact of the
slightly different PHOENIX-COND models used by these authors can be tested, by
paying attention to the prediction error between the tables of Claret (2017)
using the q-method and the predictions made using the limb-darkening code of
Espinoza & Jordán (2015), which for phoenix use the r-method by default. Here,
thus, we also take the opportunity to test which of these versions of the
phoenix models is actually preferred to minimize mismatches with actual
empirical LDCs obtained from TESS lightcurves.

####  2.3.2 Synthetic-Photometry/Atmosphere-Model LDCs

Obtaining LDCs from fits to model intensity profiles is an inherently
different optimization problem to that of obtaining LDCs from fits to
transiting exoplanet lightcurves. They have different optimization functions
as they are optimizing different supposed observables. As shown by Howarth
(2011), this fact may give rise to completely different LDCs, depending on the
adopted fitting method, even if the same coefficients are used to generate
both the intensity profile and the transiting exoplanet lightcurve. In order
to take this into account when comparing theoretical to empirically obtained
LDCs through transit lightcurve fitting, Howarth (2011) proposed an
alternative way of computing theoretical LDCs: the Synthetic-
Photometry/Atmosphere-Model (SPAM).

The first step of the SPAM algorithm is to calculate a “synthetic” transit
lightcurve with an accurate representation of the stellar intensity profile
(in the case of the present work, described by the non-linear law) assuming
full knowledge of planetary and stellar properties (although this can be
relaxed; see Espinoza & Jordán, 2015). This synthetic transit lightcurve is
then fitted by a lightcurve model with fixed planetary properties and any
limb-darkening law of interest with free LDCs. The retrieved LDCs through this
procedure are the SPAM LDCs. Here, we follow this latter procedure by fitting
the synthetic lightcurves with a quadratic limb-darkening law. We use the
ExoCTK777https://github.com/ExoCTK/exoctk
see also, https://github.com/Jayshil/ld-spam/blob/main/p2.py for our
implementation (Bourque et al., 2021) package to compute the SPAM LDCs using
the non-linear LDCs derived from the various stellar model atmospheres
discussed in Section 2.3.1.

##  3 Results

Having fitted the TESS transit lightcurves of the 176 transiting exoplanets in
our sample, we are now in a position to compare the empirically obtained LDCs
with the theoretical ones described in the previous section. In what follows,
we present these results in two parts. In the first part, we compare the
retrieved planetary parameters obtained from our juliet fits to the TESS data
with their corresponding literature values in order to validate our modelling
procedure. Then, we compare the retrieved LDCs from those transit fits with
the ones calculated theoretically from model stellar atmospheres obtained as
discussed in Section 2.3. The codes which we used to produce all the results
in this section are publicly available on
Github888https://github.com/nespinoza/tess-limb-darkening,
https://github.com/Jayshil/ld-project-updated, and,
https://github.com/Jayshil/my_thesis.

###  3.1 Comparison with literature values

Our juliet fits allow us to estimate various planetary parameters along with
the limb darkening coefficients. The former, in turn, provide us with an
excellent dataset which is useful not only to refine these planetary
parameters, but also to validate our procedures by comparing them to the ones
already estimated in the literature. While such a validation in principle
assumes these literature values are unbiased, we believe the risk of this not
being the case is somewhat mitigated by the wide range of instruments,
analyses and assumptions made by different authors that analyzed these systems
in the past. To perform this comparison, we chose three planetary parameters
to compare in this work, which are the ones for which we obtain the most
precise constraints: the scaled semi-major axis ($a/R_{*}$), the planet-to-
star radius ratio ($R_{p}/R_{*}$) and the time of transit center ($t_{c}$).

Figure 2:  (Left) Phase-folded lightcurve for the systematics-removed WASP-62b
system with data taken by TESS combining the observations of all sectors (top
left plot; blue datapoints), and the residuals from the best-fit model (bottom
left plot). The white circles with black errorbars are binned datapoints in
2-minute intervals in this phased lightcurve. (Right) Power spectral density
(PSD) of the residuals of the best fit, which doesn’t show any significant
signal on the frequencies relevant to the transit event. Dashed lines from top
to bottom show where 99.9%, 99% and 90% of the PSD of simulated gaussian noise
with the same properties at those of the errorbars on our data falls; the
observed PSD of the residuals falling below these lines suggests it is
consistent with that of white gaussian noise). Figure 3: Corner plot of
various retrieved parameters for WASP-62b system

Before presenting the main results of our study for all the planets in our
sample, we first present a representative fit for the WASP-62b system as
observed by TESS (including all sectors on which this exoplanet was observed).
The summary of this fit is presented in Figure 2. As can be seen, the phase-
folded transit lightcurve has exquisite precision, and doesn’t give rise to
any noticeable leftover signal in the residuals. Indeed, a power spectral
density analysis reveals no significant signals are present in the residual
time-series, suggesting our Exponential Matern kernel did take care of any
correlated noise structure in the data. In terms of achieving the photon-noise
level of these observations, for some sectors a significant jitter term of a
couple of hundreds of ppm is found. This is the case for Sector 1, 4, 7, 8 and
Sector 27 through 34\. This probably points to the fact that the PDC algorithm
has some extra uncorrected systematics in those sectors, but the amplitude
seems to be overall small and, judging from our power spectral density
analysis, it is unimportant in terms of correlated noise in the time-scales of
interest for our analyses ($\sim$ hours). The posterior distribution of the
transit lightcurve parameters for this target are presented in Figure 3, where
we have transformed the $q_{1}$ and $q_{2}$ limb-darkening parameters back to
the $u_{1}$ and $u_{2}$ plane using the transformations in Kipping (2013).
These posterior parameters, in turn, are in very good agreement with previous
values in the literature. For instance, the planet-to-star radius ratio
($R_{p}/R_{*}$) from our fits is $0.1111\pm 0.0001$, which is within
$3-\sigma$ from the value found by Maciejewski (2020) of $0.1105\pm 0.0003$. A
similar agreement is found for the rest of the transit parameters: our
estimated scaled semi-major axis ($a/R_{*}$) is $9.72\pm 0.03$, which is
within 1-sigma from the value found by Stassun et al. (2017a) of $9.55\pm
0.41$. Finally, the difference between the transit center predicted by that
same work on the time-span of the TESS observations and the one estimated by
our transit fits are also in very good agreement agreement — they show a
difference of $28.757$ seconds, which is consistent with zero within 1-sigma.
The predicted ephemerides from the work of Bonomo et al. (2017) has
considerably deteriorated since its publication, showing a bias of around 12
minutes for this target. Our TESS updated ephemerides, however, significantly
improve its precision to only 3 seconds.

Figure 4: Comparison of the scaled semi-major axis ($a/R_{*}$) calculated
using juliet with its corresponding literature values. The observed values and
the values from the literature are on the y-axis and x-axis respectively in
the upper panel, while the lower panel shows the residuals between those two
values. Figure 5: Comparison of the planet-to-star radius ratio
($R_{p}/R_{*}$) calculated using juliet with its corresponding literature
values. The observed values and the values from the literature are on the
y-axis and x-axis respectively in the upper panel, while the lower panel shows
the residuals between those two values. Figure 6: Comparison of the time of
transit center ($t_{c}$) calculated using juliet with its corresponding
literature values. The observed values and the values from the literature are
on the y-axis and x-axis respectively in the upper panel, while the lower
panel shows the residuals between those two values.

The same comparison on the retrieved transit parameters done above for
WASP-62b is presented for all the targets in our sample in Figures 4, 5 and 6.
Literature values for the planetary parameters were obtained from the NASA
Exoplanet Archive, as queried on February 23, 2021. The comparison is also
presented in Table 6. Our results show in general very good agreement with
literature values. There are, however, a set of 25 exoplanetary systems in our
sample for which we find that one or more than one of the planetary parameters
compared in Figures 4, 5 and 6 (i.e., the time-of-transit center, $a/R_{*}$
and/or $R_{p}/R_{*}$) are at least $3-\sigma$ away from what is published in
the literature. A detailed, case-by-case analysis is presented in Appendix A
for those systems, but we found that for 19 of them we are confident our
retrieved values are updates over previous published values for these
parameters, and that for 5 of them the most likely explanation is either
stellar activity variations producing slight transit-depth variations, real
planetary variability in the transit depths and/or mismatches on the dilution
factors applied by the TESS PDC pipeline. The only target in which we consider
the methods here presented failed to retrieve the correct properties is
LTT9779 b, for which we sample a discrepant solution to the one described in
Jenkins et al. (2020) with our methodology. This solution was actually briefly
discussed in Jenkins et al. (2020), but was discarded as the implied stellar
density does not agree with the spectroscopic one observed in that work.
Overall, we consider having only 1 outlier out of a sample of 176 targets is
in fact a very good result, which in turns gives us confidence in our
lightcurve analysis procedures.

Before moving on with our results, we would like to highlight the curious case
of the WASP-140 b exoplanetary system, which has a nearby but not physically
associated bright star. There are many systems in our sample that have such
companions, and our hypothesis on such targets was that, while they would
dilute the transit signal, the PDC lightcurves should have accounted for these
dilutions in the final reported fluxes. However, it seems the correction in
the case of WASP-140 b was either not appropriate or not sufficient to account
for it based on a comparison of our original lightcurve fits with those
published in the literature by Hellier et al. (2017). We thus decided to take
care of this dilution via the so-called dilution factor in the modelling
process (see Espinoza et al., 2019a). Instead of fixing this dilution factor
in the modelling procedure to 1, as done for the rest of the targets, we fit
it together with the rest of the transit parameters. Our results, however,
still were inconsistent with those of Hellier et al. (2017) even if accounting
directly for this dilution. Our fixed dilution fit gives rise to a planet-to-
star radius ratio which is over 60% larger than the one reported by Hellier et
al. (2017) — and to a much larger radius if the dilution is left as a free
parameter. In both cases, the impact parameter we retrieve is $b>1$ at
5-sigma, which although a couple of sigmas away from the value reported in
that study ($b=0.93^{+0.07}_{-0.03}$) is still consistent with it in our fixed
dilution case. We believe that a full joint analysis of the entire photometric
and radial-velocity datasets for this system is needed in order to solve this
discrepancy, but we leave such an analysis for future work. Here, we simply
discard this system from our analysis.

###  3.2 Analysis of targets in multiple sectors

Figure 7: Transit lightcurves (top; black datapoints original data, white
points binned data, blue curves best-fit lightcurves) and transit depths
(bottom; blue points with errorbars) on each sector for one of our best
constrained targets in this work: WASP-62b. As can be seen in the bottom plot,
the retrieved transit depths are consistent between sectors, and they converge
to a value ($12341^{+32}_{-33}$ ppm; grey band representing this 68%
credibility band) which is in agreement with the literature value ($12210\pm
66$ ppm). Figure 8: Limb-darkening coefficients extracted for our best-
precision targets (the ones spanning data in multiple sectors), compared
against their theoretical predictions (blue for atlas models, red for
phoenix). As can be seen, a clear offset is observed for the $u_{2}$ limb-
darkening coefficient. The bands show the uncertainties in the mean offset
with the red and blue colors representing the phoenix and atlas models. The
method to compute theoretical LDCs is the one from Espinoza & Jordán (2015).

An additional good consistency check for the results presented in this work is
to perform a thorough analysis of exoplanetary systems that had data spanning
multiple sectors. This allows us to test our methods in virtually independent
datasets of the same systems which not only allows us to compare the retrieved
transit parameters among different sectors, but also extract the most precise
LDCs in our sample by combining the datasets at hand. For our analysis we
decided to use all the data available for those systems up to Sector 34, which
implies some of them had observations in sectors more recent than the ones we
selected to define our target sample. In our case, these datasets were the
ones for WASP-62 (20 sectors), TOI-481 (9 sectors), WASP-119 (12 sectors),
LHS1815 (20 sectors), TOI-157 (12 sectors), WASP-79 (4 sectors) and HATS-46 (4
sectors).

In Figure 7, we present our lightcurve fit results for WASP-62b in different
sectors, alongside with the corresponding transit depths. The sector-by-sector
fits were performed with the Exponential-Matèrn GP kernel. As can be seen, the
transit depths are mostly consistent between sectors and they actually
converge to the same value that is found in the literature. The only
significantly discrepant value is that of Sector 31, which has a transit depth
of $12687\pm 97$ ppm — $\sim 3-$sigma away from the combined transit depth.
Interestingly, our analyses show that the discrepancy in this case is mostly
driven by the selection of the GP kernel: a quasi-periodic kernel fit on this
particular sector (which gives a much better bayesian evidence) gives back a
transit depth of $12376\pm 129$ ppm; consistent with the combined depth. While
this would hint that we should perhaps allow different GP kernels to be fit on
different sectors, we found that for the combined multi-sector analysis we
perform in this work this extra complexity is not important, as we obtain the
same results either way. Not considering Sector 31, we ran a chi-square test
comparing the observed depths and errors from all sectors to that of the mean
depth across all of them. This gave back a p-value of 0.34 — with which we
fail to reject the null hypothesis that the data is consistent with gaussian
noise. This is good evidence that, indeed, the transit depths are constant
accross sectors, and consistent with the value found in the literature. As can
be seen in Table 6, the same applies to all the other planets observed in
multiple sectors. This gives us confidence that our analysis is also well
behaved between sectors.

The corresponding limb-darkening coefficients extracted from the transit
lightcurves for those targets are presented in Figure 8, where we plot the
observed minus the predicted theoretical limb-darkening coefficients, which we
obtained from Espinoza & Jordán (2015), for our targets both using phoenix and
atlas models. As can be seen, the limb-darkening coefficients are largely
consistent for the $u_{1}$ coefficient, for both atlas and phoenix models
though a small offset is still present. However, for the $u_{2}$ coefficients,
both models are evidently a poor match to the observations, producing a
comparatively large offset. As we will see below, such discrepancy is not
exclusive of the multiple-sector data: this happens with the vast majority of
the systems we analyzed.

###  3.3 Comparison between theoretical and empirical LDCs from TESS

Figure 9: Mean offset between empirically determined LDCs (through TESS transiting exoplanet lightcurves) and the theoretically derived LDCs when using a quadratic law. Different colors of the bars represent the different methodologies to compute theoretical LDCs, from LDCs Tables/Codes (Section 2.3.1, orange) and SPAM LDCs (described in Section 2.3.2, blue). The x-axis shows the various origins of these theoretical LDCs: the work of Espinoza & Jordán (2015, EJ15) and the work of Claret (2017, C17). The additional sufix for this latter work (q and r) represent two different sets of assumptions used by Claret (2017) to obtain these LDCs (see text for details). The top panel shows the offset for the $u_{2}$ coefficient while the lower panel shows the offset for the $u_{1}$ coefficient of the quadratic law. Table 2: Mean offset present in the Tabular/Code and SPAM LDCs using phoenix and atlas model stellar atmospheres with different calculation methods. Method | $\Delta u_{1}$ | $\Delta u_{2}$
---|---|---
Tabular/Code LDCs
atlas (Espinoza & Jordán, 2015) | $-0.096\pm 0.018$ | $0.145\pm 0.022$
phoenix (Espinoza & Jordán, 2015) | $-0.058\pm 0.018$ | $0.107\pm 0.022$
atlas (Claret, 2017) | $-0.101\pm 0.018$ | $0.128\pm 0.021$
phoenix \- q-method (Claret, 2017) | $-0.132\pm 0.018$ | $0.243\pm 0.022$
phoenix \- r-method (Claret, 2017) | $-0.072\pm 0.018$ | $0.107\pm 0.022$
SPAM LDCs
atlas (Espinoza & Jordán, 2015) | $-0.075\pm 0.018$ | $0.101\pm 0.022$
phoenix (Espinoza & Jordán, 2015) | $-0.019\pm 0.018$ | $0.044\pm 0.021$
atlas (Claret, 2017) | $-0.084\pm 0.017$ | $0.093\pm 0.022$
phoenix \- q-method (Claret, 2017) | $-0.034\pm 0.018$ | $0.071\pm 0.022$
phoenix \- r-method (Claret, 2017) | $-0.020\pm 0.018$ | $0.042\pm 0.022$

In the previous sections we have performed a detailed comparison between the
retrieved planetary parameters with their corresponding literature values,
finding very good agreement between the two. We consider these results as a
validation that our fits to TESS data are indeed adequately constraining the
transit lightcurve shapes.

Having validated our results with literature data on the physical planetary
parameters (relative to those of the star), we now switch our focus in order
to compare the retrieved LDCs with their corresponding theoretical values,
obtained as described in Section 2.3. The empirically determined LDCs for each
of the targets analyzed in this work as well as their theoretical predictions
are presented in Table 4 and Table 5.

We present a system-by-system comparison between the empirical and the
theoretical LDCs using the various here outlined procedures in Figures 13
through 18: from Figure 13 through 15 we provide the comparison against
theoretical LDCs obtained directly from previously published tables/codes
(i.e., following the method described in Section 2.3.1), and from Figures 16
through 18 we present the results following the SPAM methodology for obtaining
the LDCs (i.e., following the method described in Section 2.3.2). To further
summarize those results, in Figure 9 we present the mean offset between the
TESS retrieved LDCs and the theoretically predicted ones for both sets of
LDCs; these mean offsets are in turn also presented in Table 2.

As can be observed from Figure 9, the general behaviour of the offsets in the
LDCs remain the same independent of which method one decides to use to compute
the theoretical LDCs. The $u_{1}$ coefficients are consistently under-
predicted by those theoretical calculations, while the $u_{2}$ coefficients
are over-predicted. This latter coefficient, however, is the one that shows
the largest (absolute) offset accross methods. Perhaps one of the most
interesting features of these results, however, is the significantly lower
offset the SPAM LDCs show when compared against the empirically obtained LDCs
— similarly to what was observed by Kepler in the analyses of Howarth (2011)
and Espinoza & Jordán (2015). As in those works, this suggests that this is,
indeed, on average the correct way of extracting theoretical LDCs for usage in
transit lightcurve modelling. LDCs extracted without this SPAM algorithm all
show significant offsets ($>3\sigma$) with respect to empirically determined
LDCs in at least one coefficient, with the worst performing method being the
q-method of Claret (2017) which shows a mean offset on the $u_{2}$ coefficient
of $\Delta u_{2}=0.243\pm 0.022$ — a very large offset when one considers the
space of all possible coefficients for $u_{2}$ spans the range from -1 to 1.
When applying the SPAM methodology, however, most of these offsets become
consistent with zero at about 2-3 $\sigma$ levels for both $u_{1}$ and
$u_{2}$. The exceptions are the LDCs calculated using the q-method of Claret
(2017, $\Delta u_{2}=0.071\pm 0.022$) and the ATLAS LDCs, which show
significant offsets on both coefficients for both, the calculations made using
the Claret (2017) tables or the limb-darkening library of Espinoza & Jordán
(2015).

##  4 Discussion

The results presented in the previous section and summarized in Figure 9 seem
to suggest two main take-home messages. First, when using quadratic LDCs
extracted directly from limb-darkening tables, large offsets as large as
$\Delta u_{1}\approx-0.1$ and $\Delta u_{2}\approx 0.2$ can be expected, on
average, with empirically determined LDCs in the TESS bandpass, with the
offsets being worse for the PHOENIX models than for the ATLAS models (top set
of values in Table 2). If one uses the SPAM methodology to obtain these LDCs,
however, these offsets reduce, on average, to virtually zero for LDCs obtained
using the PHOENIX stellar models and the r-method with both the tables
published by Claret (2017) and the limb-darkening library of Espinoza & Jordán
(2015) (see bottom set of values in Table 2). This, in turn, suggests the SPAM
methodology is the one that should be adopted by default when using
theoretical LDCs to perform transit lightcurve modelling in the TESS bandpass,
in particular using that technique together with PHOENIX stellar model
atmospheres.

These take-home messages, however, do not tell the entire story of the
observed offsets between theoretically and empirically determined LDCs. Given
different stars have different limb-darkening profiles, we must be careful
with interpreting these results especially considering the overabundance of F
and G-type stars in our sample. With this motivation in mind, we discuss the
variation of these offsets with host star temperature in the next sub-section.

###  4.1 Variation in offset with effective temperature of the host star

For Tabular/Code LDCs

Figure 10: Temperature variation of offsets in $u_{1}$ and $u_{2}$, when,
theoretically, LDCs are calculated using code provided by Espinoza & Jordán
(2015). The dashed-dotted lines show the best fitted model to the residuals.

For SPAM LDCs

Figure 11: Same as Figure 10, but now using SPAM LDCs.

As described in the previous section, there are suggestions of mean offsets
when comparing empirical to theoretically determined LDCS. The works of Müller
et al. (2013b) and Espinoza & Jordán (2015), in turn, when performing a
similar analysis on Kepler lightcurves, observed a possible dependence of
those offsets with stellar temperature. Motivated by this, we here explore
that dependence using our TESS LDCs in what follows.

To check for any such possible dependencies, we plot the residuals between the
theoretical and the TESS retrieved LDCs as a function of the effective
temperature for both sets of theoretical LDCs (those derived from limb-
darkening Tables/Codes and the SPAM LDCs) in Figures 10 and 11 for the
representative case of the Espinoza & Jordán (2015) LDCs (which are fairly
similar to the ones tabulated by Claret (2017) using the r-method999For
completeness, we show the same plots for the LDCs obtained using the tables in
Claret (2017) in Figures 19 through Figure 22) in the Appendix.) To find any
possible correlation between these offsets and the effective temperature of
the host stars, we fit those residuals with polynomials ranging from zeroth
(i.e., a constant model) to second order. We used the Bayesian Information
Criterion (BIC; Schwarz, 1978) to determine the best fitting model among those
three models, and present the best of them (i.e., the one with the minimum
BIC) in Figures 10 and 11 with dashed-dotted lines.

In almost all of those residual plots we find a quadratic function seems to
fit the residuals best, with the vertex of the parabola touching (or being
close to) zero at about 6000 K. The residuals, in turn, seem to increase their
absolute values as they go to cooler and hotter stars. For instance, for the
PHOENIX models, these quadratic functions predict offsets of order $\approx
0.2$ for both SPAM and non-SPAM coefficients. Indeed, the SPAM LDCs seem to
significantly flatten those parabolas, but indeed, the LDCs especially for
cooler stars don’t seem to be consistent with zero. This suggest that
independent of the methods, LDCs for stars cooler than about 5000 K are not
predicted very precisely by our modelling efforts, and thus care must be taken
when using theoretical LDCs for stars at those temperatures.

While interpreting these correlation between the effective temperature and the
offsets, one needs to keep in mind that the distribution of our sample is not
uniform in the effective temperature, as shown in Figure 1 – these
distribution shows the deficiency of targets at low and very high
temperatures. Still, it is interesting to see how the residuals between
theoretical and observed LDCs seem to be better at temperatures close to that
of the Sun, which makes intuitive sense given most stellar model atmospheres
have been extensively tested against the Sun.

###  4.2 Limitations of our study

An important limitation of our study to have in mind is the fact that our work
requires, by construction, a layered approach to obtaining limb-darkening
coefficients from stellar model atmospheres. In other words, the comparison is
not direct to the models, but involves a series of assumptions and methods
that lead to the final coefficients we compare against the data. As such, it
is somewhat complicated to make direct statements about the actual validity of
the stellar model atmospheres tested in our work. For instance, while we
observe that the PHOENIX coefficients obtained through the r-method behave
better when compared against our retrieved empirical LDCs, we cannot confirm
if this is because some bias cancels out between the methodology and the
actual stellar model atmospheres, or because the stellar models are inherently
reasonably well behaved in the TESS bandpass, and the r-method is
intrinsically better. We can, however, make the practical suggestion that,
when confronted with the data, and if willing to use tabulated LDCs, one
should prefer this method if using the PHOENIX models over the q-method
proposed by Claret (2017).

##  5 Conclusion and Future Work

In this paper, our main goal was to study the limb darkening of exoplanet host
stars using TESS transiting exoplanet lightcurves. To achieve this goal, we
used precise exoplanet transit lightcurves obtained by TESS for 176 known
exoplanetary systems, for which we computed limb-darkening coefficients — the
primary products needed to meet our work’s goal. A secondary but not less
important product of our work are updated transit parameters for these
exoplanetary systems. A subset of those ($R_{p}/R_{*}$, $a/R_{*}$ and the
time-of-transit centers) were compared with their corresponding values from
the literature with which we found very good agreement. This not only provides
a valuable validation of our transit lightcurve fits, but also a rich dataset
for the community to use to, e.g., use as prior in precise follow-up
observations of these systems, or to even plan observations of these systems
given the improved precision of our ephemerides for most of our studied
systems.

We used these results to test how well the limb darkening effect is modelled
by current methodologies and stellar model atmospheres. We compared the
retrieved LDCs with both ATLAS and PHOENIX model predictions calculated by two
authors, Claret (2017) and Espinoza & Jordán (2015) for one of the most widely
use limb-darkening laws: the quadratic law. We found significant offsets on
the limb-darkening coefficients $u_{1}$ and $u_{2}$ of the quadratic law
between empirically and theoretically determined LDCs, which can be as large
as $\Delta u_{1}\approx 0.1$ and $\Delta u_{2}\approx 0.2$, depending on the
methodology and the stellar model atmospheres used for the calculations. Our
main take-home message is that the most accurate LDCs for the TESS bandpass
are the ones that use the SPAM algorithm proposed in Howarth (2011), and which
use the PHOENIX model atmospheres using the technique introduced in Wittkowski
et al. (2004) and used by Espinoza & Jordán (2015) to calculate limb-darkening
coefficients for these stellar model atmospheres — the r-method in the
terminology used in the published tables of Claret (2017).

Our analysis, however, suggests that even with the “right” methodology, large
offsets between theoretical and empirically determined LDCs could be observed
for stars cooler than about 5000 K. Below these temperatures, there seem to be
offsets as large as $\Delta u_{i}\sim 0.2$ for the LDCs of the quadratic law.
We believe an important future avenue would be to extend our work to a larger
sample of exoplanetary systems orbiting cooler stars in order to confirm this
trend. The newest of TESS discoveries are indeed paving the way to form such a
sample, which we believe would be fundamental not only to understand the
physical properties and architectures of these exoplanetary systems, but also
to understand the limitations of current models and methodologies at
predicting the imprint of the limb-darkening effect on real transit
lightcurves.

## Acknowledgements

We would like to thank an anonymous referee for their thorough feedback which
significantly improved the presentation of our results. JAP wants to thank Dr.
Néstor Espinoza and Mr. Ashok Patel for providing the financial support to
visit the Max-Planck-Institut für Astronomie (MPIA), Heidelberg in Germany to
perform the present work. A significant part of the present work constitute
the final year project work of JAP’s Master’s degree. JAP wants to thank his
parent institute – Sardar Vallabhbhai National Institute of Technology,
Surat-7, Gujarat, India – and his academic supervisor Prof. K. N. Pathak to
give permission to work at MPIA during this period. NE would like to thank M.
Günther and T. Daylan for useful discussions regarding uniform analyses of
TESS known exoplanets.

This research has made use of the NASA Exoplanet Archive, which is operated by
the California Institute of Technology, under contract with the National
Aeronautics and Space Administration under the Exoplanet Exploration Program.

This research made use of the open source Python package exoctk, the Exoplanet
Characterization Toolkit (Bourque et al., 2021).

##  Appendix A Discrepant planetary systems with literature values

As noted in Section 3.1, our lightcurve analysis showed that for 25
exoplanetary systems, one or more than one of the planetary parameters are at
least $3-\sigma$ away from what is published in the literature. These systems
are LTT9779 b, XO-3 b, WASP-161 b, WASP-121 b, HAT-P-2 b, HAT-P-69 b, KELT-11
b, CoRoT-18 b, WASP-4 b, K2-237 b, WASP-131 b, WASP-17 b, KELT-20 b, WASP-31
b, WASP-46 b, WASP-19 b, WASP-92 b, XO-6 b, TOI-157 b, WASP-22 b, K2-260 b,
WASP-7 b, WASP-95 b, CoRoT-32 b and HATS-18 b. Here, we perform a more in-
depth discussion on what might be producing these offsets on the planetary
parameters.

For 11 of those “discrepant” systems (HAT-P-69 b, KELT-11 b, WASP-4 b,
WASP-131 b, KELT-20 b, WASP-31 b, WASP-46 b, WASP-22 b, K2-260 b, WASP-95 b
and CoRoT-32 b) the offset is in the predicted time-of-transit center versus
our computed time-of-transit center, which could either point to possible
transit-timing variations (TTVs) in these systems (e.g., WASP-4b is known to
have TTVs, which we recover here; see Bouma et al., 2019) or simply act as
updates with respect to previous ephemerides which might be outdated. Similar
offsets have been found for HAT-P-69 b and WASP-95 b in Shan et al. (2021);
however, the timing offsets for KELT-11 b, WASP-131 b, KELT-20 b, WASP-31 b,
WASP-46 b, WASP-22 b, K2-260 b and CoRoT-32 b reported here have not been
reported elsewhere to our knowledge.

The remaining 14 exoplanetary systems (LTT9779 b, XO-3 b, WASP-161 b, WASP-121
b, HAT-P-2 b, CoRoT-18 b, K2-237 b, WASP-17 b, WASP-19 b, WASP-92 b, XO-6 b,
TOI-157 b, WASP-7 b and HATS-18 b) all show planet-to-star radius ratio
($R_{p}/R_{*}$) and/or semi-major axis-to-stellar radius ratio ($a/R_{*}$)
discrepancies between 3 to 5 sigma from the literature values; some also show
time-of-transit center offsets together with those.

From the planetary systems that show discrepancies on $a/R_{*}$ (5
exoplanetary systems: XO-3 b, HAT-P-2 b, K2-237 b, XO-6 b and TOI-157 b),
there are different explanations for the offsets:

  * •

For XO-3 b, while the reported value in this work of
$a/R_{*}=7.09^{+0.24}_{-0.23}$ is inconsistent with the value in the discovery
paper of (Johns-Krull et al., 2008, $a/R_{*}=4.95^{+0.18}_{-0.18}$), our value
is consistent with that of Wong et al. (2014) of
$a/R_{*}=7.052^{+0.076}_{-0.097}$. The value in Wong et al. (2014) has much
better precision as the orbit is constrained through the entire phase-curve of
the exoplanet.

  * •

For HAT-P-2 b, we obtain $a/R_{*}=9.04^{+0.19}_{-0.18}$. This value is, in
fact, consistent with various values in the literature (Pál et al., 2010;
Southworth, 2010; Stassun et al., 2017b), including the discovery paper
($a/R_{*}=9.770^{+1.100}_{-0.020}$; Bakos et al., 2007). However, it is
inconsistent with the value found by Loeillet et al. (2008) of
$a/R_{*}=10.28^{+0.12}_{-0.19}$. That value, however, was obtained by fitting
a radial-velocity only dataset for HAT-P-2 b which includes a measurement of
the Rossiter-McLaughlin (RM) effect. Pál et al. (2010) used the same dataset
but not considering the radial-velocities obtained in-transit, and considering
a wide array of ground-based photometric follow-up transits to obtain their
value of $a/R_{*}=8.99^{+0.39}_{-0.41}$ which is fully consistent with our
solution — this suggests the Loeillet et al. (2008) might be biased given its
constrain mainly comes from the RM effect. We consider our parameter here an
update to the parameters of this system given these previous attempts.

  * •

K2-237 b was observed by Kepler/K2 back in 2016, and independant analyses were
performed by the teams of Soto et al. (2018, $a/R_{*}=5.50^{+0.15}_{-0.11}$)
and Smith et al. (2019, $a/R_{*}=5.503^{+0.015}_{-0.207}$), both of which give
consistent $a/R_{*}$ values with each other, but inconsistent values at
3-sigma with the ones reported in our present work
($a/R_{*}=6.17^{+0.13}_{-0.19}$; see Table 6), especially when compared to the
work of Smith et al. (2019). Interestingly, our value of $a/R_{*}$ is
perfectly matched with the recent full re-analysis of the system made by
Ikwut-Ukwa et al. (2020, $a/R_{*}=6.07^{+0.14}_{-0.18}$). We believe that our
offset with the work of Smith et al. (2019) is the product of a simple typo in
their upper errorbars ($0.015$), which are one order of magnitude better than
what can be reasonable achieved with the K2 data-quality; it is interesting,
however, that the two analyses made on the (long-cadence) K2 photometry are
systematically lower than the TESS short-cadence analyses presented in this
work and that of Ikwut-Ukwa et al. (2020).

  * •

For XO-6 b, we also believe our value corresponds to an update with respect to
previous constraints on this parameter. First, the value we obtain in our work
($a/R_{*}=8.17^{+0.07}_{-0.07}$) is consistent with an independent analysis
made on the same TESS data by Ridden-Harper et al. (2020, $a/R_{*}=8.383\pm
0.074$). These values are however inconsistent with the one reported in
Crouzet et al. (2017) who from the photometry alone conclude on
$a/R_{*}=9.20\pm 0.19$, but when combining that data with Doppler tomographic
results settle on $a/R_{*}=9.08\pm 0.17$. Interestingly, the constrain on this
parameter using only the Doppler tomography in that work is
$a/R_{*}=8.30^{+1.2}_{-0.8}$, which is fully consistent with the value that
both our work and that of Ridden-Harper et al. (2020) retrieve using the TESS
photometry. It is likely, thus, that the photometric analysis in Crouzet et
al. (2017) is somewhat biasing their result towards a larger $a/R_{*}$, but it
is unclear which part of their analysis could give rise to such a large bias.

  * •

For TOI-157 b, our value for $a/R_{*}$ of $6.30\pm 0.09$ is inconsistent at
3-sigma with that reported in Nielsen et al. (2020) of $5.79^{+0.07}_{-0.07}$.
One of the differences between our analyses is that Nielsen et al. (2020) use,
along some ground-based photometric follow-up, 8-sectors-worth of 30-min
cadence data and 4-sectors-worth of 2-min cadence data, which arguably mostly
define the planetary properties. In our analysis, however, we use 12-sectors-
worth of 2-minute cadence data. It is interesting to note, in addition, that
the differences on this parameter between the Nielsen et al. (2020) estimate
and our work go in the same direction as those observed and discussed for
K2-237 b before: our analysis obtains a larger value for $a/R_{*}$. The only
observational similarity between the TESS observations of TOI-157 b analyzed
by Nielsen et al. (2020) and the K2 observations of K2-237 b analyzed by the
teams of Soto et al. (2018) and Smith et al. (2019) is that both datasets rely
heavily on long-cadence data. This is suspiciously consistent with what is
expected by morphological lightcurve distortions due to finite integration
time (Kipping, 2010): not properly accounting/resampling the lightcurves in
long-cadence observations would give rise to smaller values of $a/R_{*}$ for a
fixed period. We suggest, thus, that perhaps the value obtained in Nielsen et
al. (2020) is due to the fact that lightcurve resampling following procedures
similar to those outlined in Kipping (2010) were either not applied or not
performed with sufficient precision to properly account for the effect.
Similarly, perhaps this latter option is the case as well for the K2
observations of K2-237 b analyzed by both Soto et al. (2018) and Smith et al.
(2019).

Finally, for the 9 systems that show $R_{p}/R_{*}$ discrepancies with
literature values (LTT9779 b, WASP-161 b, WASP-121 b, CoRoT-18 b, WASP-17 b,
WASP-19 b, WASP-92 b, WASP-7 b and HATS-18 b), different explanations exist:

  * •

For LTT9779 b, we retrieve $R_{p}/R_{*}=0.0337^{+0.0011}_{-0.0009}$, which is
inconsistent with the value found by Jenkins et al. (2020) of
$R_{p}/R_{*}=0.0455^{+0.0022}_{-0.0017}$. In fact, for this system, our
retrieved value for $a/R_{*}=7.62^{+0.30}_{-0.36}$ is also significantly
discrepant with that of Jenkins et al. (2020, $a/R_{*}=3.88\pm 0.09$). Our
retrieved parameter values for this system, however, are most likely wrong and
should not be used — we only present them here for completeness and
transparency of our process. While the transit fits from both our solutions
and the ones reported in Jenkins et al. (2020) both fit the data equally well,
the combination of our retrieved $a/R_{*}$ with the planetary period yield a
stellar density of $13342\pm 1744$ kg/m3 which is completely inconsistent with
that obtained through spectroscopy by Jenkins et al. (2020) of $1810\pm 130$
kg/m3. This two-solution problem was in fact briefly studied by Jenkins et al.
(2020) as well in their Methods section.

  * •

For WASP-161 b, we retrieve $R_{p}/R_{*}=0.0751^{+0.0009}_{-0.0008}$, which is
much more precise but discrepant at 3-sigma with the one reported by Barkaoui
et al. (2019, $R_{p}/R_{*}=0.0671^{+0.0017}_{-0.0017}$). While this target
does have a nearby companion about 16” to the SE, the PDC algorithm takes this
dilution into account on the photometry and as such it is unlikely this is the
cause of the discrepancy. The work of Barkaoui et al. (2019) has only one full
ground-based transit (with the rest of the follow-up photometry being partial
transits) and, as such, we believe our value for $R_{p}/R_{*}$ is effectively
an update on this parameter. In addition, we also find a discrepant time-of-
transit center, which has also been reported by Shan et al. (2021).

  * •

For WASP-121 b, we retrieve $R_{p}/R_{*}=0.1217^{+0.0003}_{-0.0003}$, which is
discrepant at about 3-sigma with the value of
$R_{p}/R_{*}=0.1245^{+0.0005}_{-0.0005}$ reported in Delrez et al. (2016).
Interestingly, the Delrez et al. (2016) value is consistent with analysis of
Sector 7 TESS data from other teams (see, e.g., Yang et al., 2020; Bourrier et
al., 2020; Daylan et al., 2021). In our work, however, we use additional data
from Sectors 33 and 34. If we run our analyses on Sector 7 TESS data only, our
resulting value of $R_{p}/R_{*}=0.12394^{+0.00047}_{-0.00043}$ is consistent
with both, the value of Delrez et al. (2016) and the rest of the TESS analyses
in the literature. However, individual analyses on Sectors 33 and 34 reveal
that the transit depth on Sector 7 is significantly larger than the ones
observed in those: 500 ppm larger. This might, indeed, be true variability in
the transit depth caused either by the star or the planetary atmosphere itself
— this has already been suggested by ground-based observations by Wilson et
al. (2021). We note that the transit depth reported in this work is consistent
with the average transit depth in the HST/WFC3 transit spectrum presented by
Evans et al. (2018). We also observe a 3-sigma time-of-transit offset from
that of Delrez et al. (2016), which might also hint to possible long-term
TTVs.

  * •

For CoRoT-18 b, we obtain $R_{p}/R_{*}=0.1579^{+0.0042}_{-0.0046}$, which is
significantly discrepant with the $0.1341^{+0.0019}_{-0.0019}$ value reported
in Hébrard et al. (2011). This latter reported value, in turn, agrees well
with follow-up ground-based photometry of Southworth (2012). While the
discrepant transit depth in the TESS photometry could be due to variability
either in the star or the planet itself, given the field containing this
target is so crowded, it is also possible that the discrepancy in the planet-
to-star radius ratio is due to an over-correction of the dilution from nearby
sources, similar to the case if WASP-140 b in Section 3.1.

  * •

For WASP-17 b, we retrieve $R_{p}/R_{*}=0.1218^{+0.0016}_{-0.0014}$, which is
inconsistent at 3-sigma with the value reported in Anderson et al. (2010,
$R_{p}/R_{*}=0.1293^{+0.0008}_{-0.0008}$). We believe, however, that our value
is an update on this parameter given the agreement of our value with recent
precise HST and Spitzer transit spectroscopy reported by Saba et al. (2021).

  * •

For WASP-19 b, we obtain $R_{p}/R_{*}=0.1519^{+0.0018}_{-0.0020}$, which
agrees with some results in the literature, but not with others. In
particular, this value agrees with the originally reported value by Hebb et
al. (2010) and the recent TESS analysis performed by Wong et al. (2020), but
does not with the values reported by Tregloan-Reed et al. (2013), Mancini et
al. (2013) and Espinoza et al. (2019b) by more than 3-sigma. This is, in turn,
not surprising, given the high levels of stellar activity the star is known to
have, which has been observed to directly contaminate the observed transit
depths both due to occulted and unnoculted spots (see, e.g., Espinoza et al.,
2019b).

  * •

For WASP-92 b, we retrieve $R_{p}/R_{*}=0.1068^{+0.0015}_{-0.0022}$, which is
inconsistent at 3-sigma with the value reported in Gajdoš et al. (2019) of
$R_{p}/R_{*}=0.0963^{+0.0017}_{-0.0017}$. While the there seems to be an
extended object (the galaxy LEDA 2387340; TIC 10000706526) about 14.8” to the
NW of WASP-92, this object is indeed in the TIC catalog, and correctly
identified as an extended object in it, which means it is correctly being
introduced in the PDC dilution correction. It is, thus, unlikely this is the
source of the discrepancy. We believe our value is an update to the one
reported in Gajdoš et al. (2019), however, because that work does not have
full coverage of the transit event with all their ground-based photometric
follow-up. The TESS dataset we use here, however, is composed of 3 sectors
full of data which comprise a couple of tens of transits in total.

  * •

For WASP-7 b, we obtain $R_{p}/R_{*}=0.0790^{+0.0007}_{-0.0006}$, which is
inconsistent at more than 5-sigma from the
$R_{p}/R_{*}=0.0956^{+0.0016}_{-0.0016}$ value reported by Southworth et al.
(2011). Interestingly, our value is consistent with the discovery paper
$R_{p}/R_{*}$ reported by Hellier et al. (2009). The larger $R_{p}/R_{*}$
reported by Southworth et al. (2011) comes from only one (high-precision)
transit, whereas the value reported in this work comes from four high-
precision transit events from Sector 27. We would, thus, be tempted to think
of our value as the “correct” one. However, it is important to note that the
third transit observed by TESS shows a possible spot-crossing event feature
(see Figure 12), suggesting that perhaps the varying transit depths accross
different studies is indeed real, and produced by occulted and unnoculted
stellar heterogeneities such as the ones observed for WASP-19 b.

  * •

Finally, for HATS-18 b, we obtain $R_{p}/R_{*}=0.1456^{+0.0021}_{-0.0024}$,
which is inconsistent with the $R_{p}/R_{*}=0.1320^{+0.0004}_{-0.0004}$ value
reported in Penev et al. (2016) by more than 5-sigma. This target is quite
special because while only observed by TESS in Sector 10, its short, 20-hour
period implies 26 full transits were used for the analysis presented in this
work. All the rest of the parameters are consistent with those of Penev et al.
(2016). This target does not have any significant, nearby contaminant — so it
is isn’t likely the offset in $R_{p}/R_{*}$ is due to miscalculated dilution.
The target is, however, active — Penev et al. (2016) measure rotational
variability with a period of about 10 days. Our best explanation for this
discrepancy, thus, is that the different transit depths could arise due to
different levels of stellar activity in the TESS observations as compared to
those observed by Penev et al. (2016) — i.e., varying transit depths due to
both occulted and unnoculted spots, similar to the WASP-19 b case. There is
evidence, in fact, of varying levels of stellar activity between epochs for
HATS-18. While Penev et al. (2016) measure a peak-to-peak amplitude of this
variability of order 20 mmag (20,000 ppm) in the HAT-South photometry, we see
peak-to-peak amplitudes of this variability in our TESS photometry at least
1/3 of that — with a peak-to-peak amplitude of about 6,000 ppm. HATS-18 b
transiting brighter regions of the star during the TESS observations as
compared to the previous transits observed by Penev et al. (2016) could
explain, for instance, the larger transit depths observed in our work.

Figure 12: Third transit of WASP-7b observed by TESS in Sector 27. An event,
suggestive of a spot crossing event, is seen at about 24 minutes before mid-
transit.

From the above analyses, thus, it seems that 24 out of 25 discrepant systems
are really mostly updates to existing reported planetary parameters in the
literature for the exoplanets under consideration in this work. Our only real
outlier is LTT9799 b. We consider having one confirmed outlier out of a sample
of 176 targets is perfectly consistent with random chance and, thus, we
consider this a very good check that our procedures are giving consistent
(and/or updated) results to previous studies.

##  Appendix B Tables

Table 3: Various stellar properties of our targets. The data was retrieved from the NASA Exoplanet Achieve on February 23, 2021. Only first 20 rows are shown here; the full table is available on the website in a machine readable ASCII form. Star Name | $T_{eff}$ | [M/H] | $\log{g}$ | $V_{turb}$
---|---|---|---|---
| (K) |  | (cgs) | (km/s)
WASP-61 | 6250.0 | -0.10 | 4.26 | …
WASP-130 | 5625.0 | 0.26 | 4.49 | …
HATS-13 | 5523.0 | 0.05 | 4.52 | …
WASP-156 | 4910.0 | 0.24 | 4.60 | …
NGTS-4 | 5143.0 | -0.28 | 4.50 | …
WASP-190 | 6400.0 | -0.02 | 4.17 | …
WASP-44 | 5420.0 | -0.003 | 4.49 | …
TOI-892 | 6261.0 | 0.24 | 4.26 | …
WASP-82 | 6480.0 | 0.12 | 3.96 | …
TOI-540 | 3216.0 | 0.00 | 4.44 | …
WASP-22 | 6000.0 | -0.05 | 4.50 | …
TOI-905 | 5570.0 | 0.14 | 4.50 | …
HAT-P-42 | 5743.0 | 0.27 | 4.14 | …
WASP-62 | 6230.0 | 0.04 | 4.45 | …
HATS-3 | 6351.0 | -0.16 | 4.22 | …
LTT9779 | 5443.0 | 0.27 | 4.35 | …
K2-260 | 6367.0 | -0.14 | 4.15 | …
TrES-3 | 5650.0 | -0.20 | 4.57 | …
Qatar-10 | 6124.0 | 0.016 | 4.30 | …
WASP-144 | 5200.0 | 0.18 | 4.53 | …

Table 4: Limb darkening coefficients calculated using different methods for
quadratic law; here EJ15 means the LDCs calculated using the method by
Espinoza & Jordán (2015) and C17 implies the tabulated values from Claret
(2017), with their r-method (C17r) and q-method (C17q) when using phoenix
models
∗ dash indicates that the LDCs for the given target could not be calculated
because the effective temperature of the host star is out of the range of
Claret (2017) ATLAS tables.
Only first 20 rows are shown here; the full table is available on the website in a machine readable ASCII form. Star name | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$ | $u_{1}^{*}$ | $u_{2}^{*}$ | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$
---|---|---|---|---|---|---|---|---|---|---|---|---
| (Empirical) | (Empirical) | (EJ15 | (EJ15 | (EJ15 | (EJ15 | (C17 | (C17 | (C17q | (C17q | (C17r | (C17r
|  |  |  atlas) |  atlas) |  phoenix) |  phoenix) |  atlas) |  atlas) |  phoenix) |  phoenix) |  phoenix) |  phoenix)
WASP-61 | $0.26^{+0.12}_{-0.13}$ | $0.19^{+0.24}_{-0.22}$ | 0.23 | 0.31 | 0.30 | 0.25 | 0.23 | 0.31 | 0.24 | 0.38 | 0.31 | 0.25
WASP-130 | $0.25^{+0.18}_{-0.16}$ | $0.43^{+0.27}_{-0.35}$ | 0.36 | 0.26 | 0.36 | 0.23 | 0.33 | 0.27 | 0.30 | 0.34 | 0.36 | 0.23
HATS-13 | $0.31^{+0.17}_{-0.17}$ | $0.29^{+0.34}_{-0.33}$ | 0.34 | 0.26 | 0.37 | 0.23 | 0.34 | 0.26 | 0.31 | 0.33 | 0.37 | 0.22
WASP-156 | $0.52^{+0.14}_{-0.15}$ | $0.08^{+0.28}_{-0.25}$ | 0.43 | 0.21 | 0.44 | 0.21 | 0.45 | 0.20 | 0.39 | 0.30 | 0.43 | 0.21
NGTS-4 | $0.31^{+0.35}_{-0.22}$ | $0.10^{+0.31}_{-0.24}$ | 0.36 | 0.25 | 0.40 | 0.21 | 0.38 | 0.23 | 0.34 | 0.33 | 0.40 | 0.22
WASP-190 | $0.42^{+0.20}_{-0.20}$ | $0.30^{+0.30}_{-0.32}$ | 0.21 | 0.32 | 0.31 | 0.24 | 0.22 | 0.32 | 0.22 | 0.39 | 0.30 | 0.25
WASP-44 | $0.44^{+0.21}_{-0.22}$ | $0.30^{+0.30}_{-0.37}$ | 0.34 | 0.26 | 0.38 | 0.23 | 0.36 | 0.25 | 0.32 | 0.33 | 0.38 | 0.22
TOI-892 | $0.21^{+0.12}_{-0.12}$ | $-0.00^{+0.17}_{-0.11}$ | 0.24 | 0.32 | 0.30 | 0.25 | 0.24 | 0.32 | 0.24 | 0.38 | 0.31 | 0.25
WASP-82 | $0.35^{+0.06}_{-0.07}$ | $0.07^{+0.12}_{-0.11}$ | 0.21 | 0.33 | 0.30 | 0.24 | 0.21 | 0.33 | 0.20 | 0.44 | 0.30 | 0.24
TOI-540 | $0.96^{+0.26}_{-0.24}$ | $-0.11^{+0.23}_{-0.22}$ | 0.24 | 0.39 | 0.14 | 0.52 | — | — | 0.02 | 0.72 | 0.12 | 0.53
WASP-22 | $0.19^{+0.11}_{-0.11}$ | $0.30^{+0.21}_{-0.21}$ | 0.26 | 0.30 | 0.32 | 0.25 | 0.26 | 0.30 | 0.26 | 0.36 | 0.32 | 0.25
TOI-905 | $0.59^{+0.38}_{-0.37}$ | $0.19^{+0.41}_{-0.44}$ | 0.34 | 0.26 | 0.36 | 0.23 | 0.34 | 0.27 | 0.31 | 0.33 | 0.37 | 0.22
HAT-P-42 | $0.25^{+0.17}_{-0.15}$ | $0.39^{+0.29}_{-0.31}$ | 0.31 | 0.29 | 0.35 | 0.24 | 0.31 | 0.28 | 0.27 | 0.39 | 0.34 | 0.24
WASP-62 | $0.30^{+0.01}_{-0.01}$ | $0.13^{+0.03}_{-0.03}$ | 0.24 | 0.31 | 0.30 | 0.25 | 0.24 | 0.31 | 0.24 | 0.37 | 0.31 | 0.25
HATS-3 | $0.22^{+0.15}_{-0.13}$ | $0.20^{+0.24}_{-0.24}$ | 0.23 | 0.31 | 0.31 | 0.24 | 0.22 | 0.31 | 0.23 | 0.39 | 0.31 | 0.24
LTT9779 | $1.35^{+0.10}_{-0.11}$ | $-0.41^{+0.12}_{-0.13}$ | 0.36 | 0.26 | 0.38 | 0.23 | 0.36 | 0.25 | 0.31 | 0.35 | 0.37 | 0.23
K2-260 | $0.29^{+0.15}_{-0.16}$ | $0.14^{+0.31}_{-0.23}$ | 0.23 | 0.31 | 0.31 | 0.24 | 0.22 | 0.31 | 0.23 | 0.39 | 0.31 | 0.24
TrES-3 | $0.45^{+0.41}_{-0.30}$ | $0.11^{+0.35}_{-0.33}$ | 0.29 | 0.29 | 0.36 | 0.23 | 0.31 | 0.28 | 0.30 | 0.33 | 0.36 | 0.23
Qatar-10 | $0.15^{+0.10}_{-0.08}$ | $0.55^{+0.20}_{-0.22}$ | 0.27 | 0.30 | 0.31 | 0.25 | 0.25 | 0.31 | 0.25 | 0.38 | 0.32 | 0.25
WASP-144 | $0.26^{+0.23}_{-0.17}$ | $0.29^{+0.34}_{-0.32}$ | 0.39 | 0.23 | 0.40 | 0.22 | 0.40 | 0.23 | 0.34 | 0.32 | 0.40 | 0.22

Table 5: SPAM limb darkening coefficients calculated using different methods
for quadratic law; here EJ15 means the LDCs calculated using the method by
Espinoza & Jordán (2015) and C17 implies the tabulated values from Claret
(2017), with their r-method (C17r) and q-method (C17q) when using phoenix
models.
Only first 20 rows are shown here; the full table is available on the website in a machine readable ASCII form. Star name | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$ | $u_{1}$ | $u_{2}$
---|---|---|---|---|---|---|---|---|---|---
| (EJ15 | (EJ15 | (EJ15 | (EJ15 | (C17 | (C17 | (C17q | (C17q | (C17r | (C17r
|  atlas) |  atlas) |  phoenix) |  phoenix) |  atlas) |  atlas) |  phoenix) |  phoenix) |  phoenix) |  phoenix)
WASP-61 | 0.30 | 0.20 | 0.35 | 0.17 | 0.30 | 0.20 | 0.37 | 0.16 | 0.37 | 0.15
WASP-130 | 0.40 | 0.19 | 0.42 | 0.14 | 0.37 | 0.20 | 0.40 | 0.16 | 0.41 | 0.14
HATS-13 | 0.38 | 0.18 | 0.42 | 0.14 | 0.38 | 0.19 | 0.42 | 0.15 | 0.42 | 0.14
WASP-156 | 0.48 | 0.13 | 0.50 | 0.11 | 0.50 | 0.11 | 0.49 | 0.12 | 0.49 | 0.11
NGTS-4 | 0.31 | 0.35 | 0.32 | 0.35 | 0.32 | 0.35 | 0.32 | 0.36 | 0.32 | 0.35
WASP-190 | 0.21 | 0.29 | 0.24 | 0.31 | 0.21 | 0.29 | 0.25 | 0.32 | 0.24 | 0.31
WASP-44 | 0.38 | 0.20 | 0.43 | 0.14 | 0.38 | 0.20 | 0.42 | 0.16 | 0.42 | 0.14
TOI-892 | 0.23 | 0.31 | 0.24 | 0.32 | 0.23 | 0.31 | 0.25 | 0.32 | 0.37 | 0.15
WASP-82 | 0.21 | 0.32 | 0.37 | 0.14 | 0.21 | 0.32 | 0.36 | 0.17 | 0.37 | 0.14
TOI-540 | 0.25 | 0.36 | 0.22 | 0.35 | -0.04 | 0.20 | 0.21 | 0.34 | 0.21 | 0.34
WASP-22 | 0.31 | 0.21 | 0.37 | 0.16 | 0.32 | 0.21 | 0.37 | 0.17 | 0.37 | 0.16
TOI-905 | 0.37 | 0.21 | 0.41 | 0.14 | 0.36 | 0.21 | 0.39 | 0.18 | 0.41 | 0.15
HAT-P-42 | 0.26 | 0.32 | 0.41 | 0.14 | 0.26 | 0.32 | 0.38 | 0.18 | 0.41 | 0.14
WASP-62 | 0.29 | 0.22 | 0.36 | 0.16 | 0.29 | 0.21 | 0.35 | 0.18 | 0.36 | 0.16
HATS-3 | 0.22 | 0.29 | 0.38 | 0.14 | 0.21 | 0.29 | 0.35 | 0.18 | 0.37 | 0.15
LTT9779 | 0.28 | 0.34 | 0.28 | 0.34 | 0.28 | 0.34 | 0.28 | 0.34 | 0.28 | 0.33
K2-260 | 0.29 | 0.21 | 0.37 | 0.14 | 0.28 | 0.21 | 0.36 | 0.16 | 0.36 | 0.15
TrES-3 | 0.31 | 0.24 | 0.41 | 0.15 | 0.33 | 0.23 | 0.38 | 0.19 | 0.40 | 0.16
Qatar-10 | 0.31 | 0.22 | 0.37 | 0.16 | 0.30 | 0.22 | 0.37 | 0.17 | 0.37 | 0.15
WASP-144 | 0.43 | 0.16 | 0.46 | 0.13 | 0.44 | 0.16 | 0.44 | 0.15 | 0.45 | 0.13

Table 6: Retrieved planetary parameters along with their literature values.
The data of literature planetary parameters was extracted from the NASA
Exoplanet Archive on February 23, 2021.
Only first 20 rows are shown here; the full table is available on the website in a machine readable ASCII form. Planet name | $R_{p}/R_{*}$ | $a/R_{*}$ |  $t_{c}$ \- 2458000 | $R_{p}/R_{*}$ | $a/R_{*}$ |  $t_{c}$ \- 2458000
---|---|---|---|---|---|---
| (This work) | (This work) | (This work) | (Literature) | (Literature) | (Literature)
WASP-61b | $0.0934^{+0.0008}_{-0.0009}$ | $8.09^{+0.11}_{-0.24}$ | $1198.73545^{+0.00033}_{-0.00034}$ | $0.0937^{+0.0031}_{-0.0031}$ | $8.15^{+0.24}_{-0.24}$ | $1198.738\pm 0.003$
WASP-130b | $0.0955^{+0.0025}_{-0.0021}$ | $24.15^{+1.42}_{-1.62}$ | $607.58558^{+0.00042}_{-0.00041}$ | $0.0955^{+0.0044}_{-0.0044}$ | $22.67^{+0.77}_{-0.77}$ | $607.586\pm 0.001$
HATS-13b | $0.1427^{+0.0033}_{-0.0028}$ | $9.76^{+0.30}_{-0.51}$ | $1083.00912^{+0.00053}_{-0.00053}$ | $0.1402^{+0.0016}_{-0.0016}$ | $9.82^{+0.18}_{-0.18}$ | $1083.006\pm 0.002$
WASP-156b | $0.0672^{+0.0009}_{-0.0008}$ | $12.86^{+0.23}_{-0.47}$ | $1169.86001^{+0.00030}_{-0.00028}$ | $0.0685^{+0.0012}_{-0.0012}$ | $12.80^{+0.30}_{-0.30}$ | $1169.861\pm 0.004$
NGTS-4b | $0.0341^{+0.0017}_{-0.0018}$ | $5.69^{+0.45}_{-1.01}$ | $1227.52707^{+0.00111}_{-0.00098}$ | $0.0350^{+0.0030}_{-0.0030}$ | $4.79^{+1.21}_{-1.21}$ | $1227.529\pm 0.010$
WASP-190b | $0.0781^{+0.0029}_{-0.0020}$ | $9.05^{+0.68}_{-0.88}$ | $1098.12518^{+0.00064}_{-0.00065}$ | $0.0741^{+0.0075}_{-0.0075}$ | $8.95^{+0.57}_{-0.57}$ | $1098.122\pm 0.001$
WASP-44b | $0.1174^{+0.0040}_{-0.0036}$ | $8.47^{+0.67}_{-0.70}$ | $393.84978^{+0.00035}_{-0.00036}$ | $0.1248^{+0.0023}_{-0.0023}$ | $8.20^{+0.43}_{-0.43}$ | $393.8498\pm 0.0004$
TOI-892b | $0.0800^{+0.0009}_{-0.0009}$ | $15.77^{+0.27}_{-0.64}$ | $1208.92591^{+0.00064}_{-0.00062}$ | $0.0790^{+0.0010}_{-0.0010}$ | $14.20^{+0.80}_{-0.80}$ | $1208.922\pm 0.005$
WASP-82b | $0.0773^{+0.0004}_{-0.0004}$ | $4.42^{+0.04}_{-0.07}$ | $449.78954^{+0.00017}_{-0.00016}$ | $0.0788^{+0.0037}_{-0.0037}$ | $4.43^{+0.15}_{-0.15}$ | $449.794\pm 0.002$
TOI-540b | $0.0462^{+0.0028}_{-0.0051}$ | $13.31^{+4.68}_{-2.43}$ | $1199.92424^{+0.00030}_{-0.00029}$ | $0.0436^{+0.0012}_{-0.0012}$ | $13.90^{+0.72}_{-0.72}$ | $1199.925\pm 0.001$
WASP-22b | $0.1001^{+0.0011}_{-0.0011}$ | $8.28^{+0.22}_{-0.30}$ | $422.50014^{+0.00023}_{-0.00023}$ | $0.0978^{+0.0012}_{-0.0012}$ | $8.94^{+0.85}_{-0.85}$ | $422.468\pm 0.008$
TOI-905b | $0.1296^{+0.0027}_{-0.0024}$ | $11.47^{+0.42}_{-0.40}$ | $639.56977^{+0.00021}_{-0.00021}$ | $0.1312^{+0.0079}_{-0.0079}$ | $10.94^{+0.50}_{-0.50}$ | $639.5695\pm 0.0002$
HAT-P-42b | $0.0811^{+0.0013}_{-0.0012}$ | $9.75^{+0.23}_{-0.48}$ | $1252.87532^{+0.00052}_{-0.00052}$ | $0.0860^{+0.0033}_{-0.0033}$ | $8.08^{+0.82}_{-0.82}$ | $1252.901\pm 0.023$
WASP-62b | $0.1111^{+0.0001}_{-0.0001}$ | $9.72^{+0.03}_{-0.03}$ | $1252.58570^{+0.00004}_{-0.00004}$ | $0.1105^{+0.0003}_{-0.0003}$ | $9.55^{+0.41}_{-0.41}$ | $1252.594\pm 0.002$
HATS-3b | $0.0973^{+0.0011}_{-0.0012}$ | $8.01^{+0.34}_{-0.34}$ | $337.89608^{+0.00031}_{-0.00031}$ | $0.1011^{+0.0006}_{-0.0006}$ | $7.42^{+0.12}_{-0.12}$ | $337.896\pm 0.003$
LTT9779b | $0.0337^{+0.0011}_{-0.0009}$ | $7.62^{+0.30}_{-0.36}$ | $1112.22021^{+0.00017}_{-0.00017}$ | $0.0455^{+0.0022}_{-0.0022}$ | $3.88^{+0.09}_{-0.09}$ | $1112.208\pm 0.009$
K2-260b | $0.0947^{+0.0016}_{-0.0017}$ | $5.11^{+0.11}_{-0.20}$ | $1186.61999^{+0.00056}_{-0.00057}$ | $0.0973^{+0.0003}_{-0.0003}$ | $5.29^{+0.03}_{-0.03}$ | $1186.604\pm 0.001$
TrES-3b | $0.1706^{+0.0058}_{-0.0039}$ | $5.82^{+0.12}_{-0.13}$ | $1034.47459^{+0.00009}_{-0.00009}$ | $0.1660^{+0.0024}_{-0.0024}$ | $5.95^{+0.05}_{-0.05}$ | $1034.472\pm 0.001$
Qatar-10b | $0.1255^{+0.0015}_{-0.0015}$ | $5.05^{+0.08}_{-0.12}$ | $1034.37308^{+0.00021}_{-0.00020}$ | $0.1265^{+0.0010}_{-0.0010}$ | $4.90^{+0.12}_{-0.12}$ | $1034.371\pm 0.005$
WASP-144b | $0.1128^{+0.0029}_{-0.0041}$ | $7.32^{+0.73}_{-0.60}$ | $1071.05883^{+0.00038}_{-0.00037}$ | $0.1079^{+0.0013}_{-0.0013}$ | $8.39^{+0.23}_{-0.23}$ | $1071.060\pm 0.001$

##  Appendix C Figures

Figure 13: Comparison of the retrieved LDCs $u_{1}$ and $u_{2}$ from the TESS
data with the one calculated from the model stellar atmospheres using the code
provided by Espinoza & Jordán (2015). In upper panel, y-axis represents the
theoretical LDCs while x-axis shows the empirical values of LDCs. The black
dashed line is the line of equality between the x- and y- axis. On the other
hand, the lower panel shows the distribution of residuals between the
theoretical and empirical LDCs.

Figure 14: Comparison of the retrieved $u_{1}$ from the TESS data with the one
calculated from the model stellar atmospheres using the tabulated values from
Claret (2017). In upper panel, y-axis represents the theoretical LDCs while
x-axis shows the empirical values of LDCs. The black dashed line is the line
of equality between the x- and y- axis. On the other hand, the lower panel
shows the distribution of residuals between the theoretical and empirical
LDCs.

Figure 15: Comparison of the retrieved $u_{1}$ and $u_{2}$ coefficients from
the TESS data with the one calculated from the model stellar atmospheres using
tabulated values from Claret (2017) using r-method and q-method.

Figure 16: Comparison of the retrieved LDCs $u_{1}$ and $u_{2}$ from the TESS
data with the SPAM LDCs. The later made use of non-linear LDCs from the code
provided by Espinoza & Jordán (2015). In upper panel, y-axis represents the
SPAM LDCs while x-axis shows the empirical values of LDCs. The black dashed
line is the line of equality between the x- and y- axis. On the other hand,
the lower panel shows the distribution of residuals between the SPAM and
empirical LDCs.

Figure 17: Comparison of the retrieved LDCs $u_{1}$ and $u_{2}$ from the TESS
data with the SPAM LDCs. The later made use of non-linear LDCs from the code
provided by Claret (2017). In upper panel, y-axis represents the SPAM LDCs
while x-axis shows the empirical values of LDCs. The black dashed line is the
line of equality between the x- and y- axis. On the other hand, the lower
panel shows the distribution of residuals between the SPAM and empirical LDCs.

Figure 18: Comparison of the retrieved $u_{1}$ and $u_{2}$ coefficients from
the TESS data with the SPAM LDCs which made use of tabulated values of non-
linear LDCs from Claret (2017) using r-method and q-method.

For Tabular/Code LDCs

Figure 19: Temperature variation of offsets in $u_{1}$ and $u_{2}$, when,
theoretically, LDCs are calculated using LDCs provided by Claret (2017). The
dashed-dotted lines show the best fitted model to the residuals.

For Tabular/Code LDCs

Figure 20: Temperature variation of offsets in $u_{1}$ and $u_{2}$, when,
theoretically, LDCs are calculated using tables provided by Claret (2017),
using q-method and r-method. The dashed-dotted lines show the best fitted
model to the residuals

For SPAM LDCs

Figure 21: Same as Figure 19, but now using SPAM LDcs.

For SPAM LDCs

Figure 22: Same as Figure 20, but now using SPAM LDCs.

## References

  * Anderson et al. (2010) Anderson, D. R., Hellier, C., Gillon, M., et al. 2010, ApJ, 709, 159, doi: 10.1088/0004-637X/709/1/159
  * Bakos et al. (2007) Bakos, G. Á., Kovács, G., Torres, G., et al. 2007, ApJ, 670, 826, doi: 10.1086/521866
  * Barkaoui et al. (2019) Barkaoui, K., Burdanov, A., Hellier, C., et al. 2019, AJ, 157, 43, doi: 10.3847/1538-3881/aaf422
  * Bonomo et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107, doi: 10.1051/0004-6361/201629882
  * Bouma et al. (2019) Bouma, L. G., Winn, J. N., Baxter, C., et al. 2019, AJ, 157, 217, doi: 10.3847/1538-3881/ab189f
  * Bourque et al. (2021) Bourque, M., Espinoza, N., Filippazzo, J., et al. 2021, The Exoplanet Characterization Toolkit (ExoCTK), 1.0.0, Zenodo, doi: 10.5281/zenodo.4556063
  * Bourrier et al. (2020) Bourrier, V., Ehrenreich, D., Lendl, M., et al. 2020, A&A, 635, A205, doi: 10.1051/0004-6361/201936640
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Castelli & Kurucz (2003) Castelli, F., & Kurucz, R. L. 2003, in IAU Symposium, Vol. 210, Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, A20.  https://arxiv.org/abs/astro-ph/0405087
  * Claret (2000) Claret, A. 2000, A&A, 363, 1081
  * Claret (2017) —. 2017, A&A, 600, A30, doi: 10.1051/0004-6361/201629705
  * Claret & Hauschildt (2003) Claret, A., & Hauschildt, P. H. 2003, A&A, 412, 241, doi: 10.1051/0004-6361:20031405
  * Claret et al. (2012) Claret, A., Hauschildt, P. H., & Witte, S. 2012, A&A, 546, A14, doi: 10.1051/0004-6361/201219849
  * Crouzet et al. (2017) Crouzet, N., McCullough, P. R., Long, D., et al. 2017, AJ, 153, 94, doi: 10.3847/1538-3881/153/3/94
  * Csizmadia et al. (2013) Csizmadia, S., Pasternacki, T., Dreyer, C., et al. 2013, A&A, 549, A9, doi: 10.1051/0004-6361/201219888
  * Daylan et al. (2021) Daylan, T., Günther, M. N., Mikal-Evans, T., et al. 2021, AJ, 161, 131, doi: 10.3847/1538-3881/abd8d2
  * Delrez et al. (2016) Delrez, L., Santerne, A., Almenara, J. M., et al. 2016, MNRAS, 458, 4025, doi: 10.1093/mnras/stw522
  * Diaz-Cordoves & Gimenez (1992) Diaz-Cordoves, J., & Gimenez, A. 1992, A&A, 259, 227
  * Espinoza & Jordán (2015) Espinoza, N., & Jordán, A. 2015, MNRAS, 450, 1879, doi: 10.1093/mnras/stv744
  * Espinoza & Jordán (2016) —. 2016, MNRAS, 457, 3573, doi: 10.1093/mnras/stw224
  * Espinoza et al. (2019a) Espinoza, N., Kossakowski, D., & Brahm, R. 2019a, MNRAS, 490, 2262, doi: 10.1093/mnras/stz2688
  * Espinoza et al. (2019b) Espinoza, N., Rackham, B. V., Jordán, A., et al. 2019b, MNRAS, 482, 2065, doi: 10.1093/mnras/sty2691
  * Evans et al. (2018) Evans, T. M., Sing, D. K., Goyal, J. M., et al. 2018, AJ, 156, 283, doi: 10.3847/1538-3881/aaebff
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Gajdoš et al. (2019) Gajdoš, P., Vaňko, M., Jakubík, M., et al. 2019, MNRAS, 485, 3580, doi: 10.1093/mnras/stz676
  * Hayek et al. (2012) Hayek, W., Sing, D., Pont, F., & Asplund, M. 2012, A&A, 539, A102, doi: 10.1051/0004-6361/201117868
  * Hebb et al. (2010) Hebb, L., Collier-Cameron, A., Triaud, A. H. M. J., et al. 2010, ApJ, 708, 224, doi: 10.1088/0004-637X/708/1/224
  * Hébrard et al. (2011) Hébrard, G., Evans, T. M., Alonso, R., et al. 2011, A&A, 533, A130, doi: 10.1051/0004-6361/201117192
  * Hellier et al. (2009) Hellier, C., Anderson, D. R., Gillon, M., et al. 2009, ApJ, 690, L89, doi: 10.1088/0004-637X/690/1/L89
  * Hellier et al. (2017) Hellier, C., Anderson, D. R., Cameron, A. C., et al. 2017, MNRAS, 465, 3693, doi: 10.1093/mnras/stw3005
  * Howarth (2011) Howarth, I. D. 2011, MNRAS, 418, 1165, doi: 10.1111/j.1365-2966.2011.19568.x
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Ikwut-Ukwa et al. (2020) Ikwut-Ukwa, M., Rodriguez, J. E., Bieryla, A., et al. 2020, AJ, 160, 209, doi: 10.3847/1538-3881/aba964
  * Jenkins et al. (2020) Jenkins, J. S., Díaz, M. R., Kurtovic, N. T., et al. 2020, Nature Astronomy, 4, 1148, doi: 10.1038/s41550-020-1142-z
  * Johns-Krull et al. (2008) Johns-Krull, C. M., McCullough, P. R., Burke, C. J., et al. 2008, ApJ, 677, 657, doi: 10.1086/528950
  * Kipping (2010) Kipping, D. M. 2010, MNRAS, 408, 1758, doi: 10.1111/j.1365-2966.2010.17242.x
  * Kipping (2013) —. 2013, MNRAS, 435, 2152, doi: 10.1093/mnras/stt1435
  * Klinglesmith & Sobieski (1970) Klinglesmith, D. A., & Sobieski, S. 1970, AJ, 75, 175, doi: 10.1086/110960
  * Kopal (1950) Kopal, Z. 1950, Harvard College Observatory Circular, 454, 1
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161, doi: 10.1086/683602
  * Kreidberg (2018) —. 2018, Exoplanet Atmosphere Measurements from Transmission Spectroscopy and Other Planet Star Combined Light Observations (Springer International Publishing), 100, doi: 10.1007/978-3-319-55333-7_100
  * Loeillet et al. (2008) Loeillet, B., Shporer, A., Bouchy, F., et al. 2008, A&A, 481, 529, doi: 10.1051/0004-6361:20078167
  * Maciejewski (2020) Maciejewski, G. 2020, Acta Astron., 70, 181, doi: 10.32023/0001-5237/70.3.2
  * Mancini et al. (2013) Mancini, L., Ciceri, S., Chen, G., et al. 2013, MNRAS, 436, 2, doi: 10.1093/mnras/stt1394
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171, doi: 10.1086/345520
  * Maxted (2018) Maxted, P. F. L. 2018, A&A, 616, A39, doi: 10.1051/0004-6361/201832944
  * Morello et al. (2017) Morello, G., Tsiaras, A., Howarth, I. D., & Homeier, D. 2017, AJ, 154, 111, doi: 10.3847/1538-3881/aa8405
  * Müller et al. (2013a) Müller, H. M., Huber, K. F., Czesla, S., Wolter, U., & Schmitt, J. H. M. M. 2013a, A&A, 560, A112, doi: 10.1051/0004-6361/201322079
  * Müller et al. (2013b) —. 2013b, A&A, 560, A112, doi: 10.1051/0004-6361/201322079
  * Neilson et al. (2017) Neilson, H. R., McNeil, J. T., Ignace, R., & Lester, J. B. 2017, ApJ, 845, 65, doi: 10.3847/1538-4357/aa7edf
  * Nielsen et al. (2020) Nielsen, L. D., Brahm, R., Bouchy, F., et al. 2020, A&A, 639, A76, doi: 10.1051/0004-6361/202037941
  * Pál et al. (2010) Pál, A., Bakos, G. Á., Torres, G., et al. 2010, MNRAS, 401, 2665, doi: 10.1111/j.1365-2966.2009.15849.x
  * Penev et al. (2016) Penev, K., Hartman, J. D., Bakos, G. Á., et al. 2016, AJ, 152, 127, doi: 10.3847/0004-6256/152/5/127
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 1 , doi: 10.1117/1.JATIS.1.1.014003
  * Ridden-Harper et al. (2020) Ridden-Harper, A., Turner, J. D., & Jayawardhana, R. 2020, AJ, 160, 249, doi: 10.3847/1538-3881/abba1e
  * Saba et al. (2021) Saba, A., Tsiaras, A., Morvan, M., et al. 2021, arXiv e-prints, arXiv:2108.13721.  https://arxiv.org/abs/2108.13721
  * Sandford & Kipping (2017) Sandford, E., & Kipping, D. 2017, AJ, 154, 228, doi: 10.3847/1538-3881/aa94bf
  * Schwarz (1978) Schwarz, G. 1978, Ann. Statist., 6, 461, doi: 10.1214/aos/1176344136
  * Schwarzschild (1906) Schwarzschild, K. 1906, Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen. Math.-phys. Klasse, 195, 41
  * Seager & Mallén-Ornelas (2003) Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038, doi: 10.1086/346105
  * Shan et al. (2021) Shan, S.-S., Yang, F., Lu, Y.-J., et al. 2021, arXiv e-prints, arXiv:2111.06678.  https://arxiv.org/abs/2111.06678
  * Sing et al. (2009) Sing, D. K., Désert, J. M., Lecavelier Des Etangs, A., et al. 2009, A&A, 505, 891, doi: 10.1051/0004-6361/200912776
  * Smith et al. (2019) Smith, A. M. S., Csizmadia, S., Gandolfi, D., et al. 2019, Acta Astron., 69, 135, doi: 10.32023/0001-5237/69.2.3
  * Soto et al. (2018) Soto, M. G., Díaz, M. R., Jenkins, J. S., et al. 2018, MNRAS, 478, 5356, doi: 10.1093/mnras/sty1334
  * Southworth (2008) Southworth, J. 2008, MNRAS, 386, 1644, doi: 10.1111/j.1365-2966.2008.13145.x
  * Southworth (2010) —. 2010, MNRAS, 408, 1689, doi: 10.1111/j.1365-2966.2010.17231.x
  * Southworth (2012) —. 2012, MNRAS, 426, 1291, doi: 10.1111/j.1365-2966.2012.21756.x
  * Southworth et al. (2011) Southworth, J., Dominik, M., Jørgensen, U. G., et al. 2011, A&A, 527, A8, doi: 10.1051/0004-6361/201016183
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Stassun et al. (2017a) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017a, AJ, 153, 136, doi: 10.3847/1538-3881/aa5df3
  * Stassun et al. (2017b) —. 2017b, AJ, 153, 136, doi: 10.3847/1538-3881/aa5df3
  * Tregloan-Reed et al. (2013) Tregloan-Reed, J., Southworth, J., & Tappert, C. 2013, MNRAS, 428, 3671, doi: 10.1093/mnras/sts306
  * Wilson et al. (2021) Wilson, J., Gibson, N. P., Lothringer, J. D., et al. 2021, MNRAS, 503, 4787, doi: 10.1093/mnras/stab797
  * Wittkowski et al. (2004) Wittkowski, M., Aufdenberg, J. P., & Kervella, P. 2004, A&A, 413, 711, doi: 10.1051/0004-6361:20034149
  * Wong et al. (2014) Wong, I., Knutson, H. A., Cowan, N. B., et al. 2014, ApJ, 794, 134, doi: 10.1088/0004-637X/794/2/134
  * Wong et al. (2020) Wong, I., Benneke, B., Shporer, A., et al. 2020, AJ, 159, 104, doi: 10.3847/1538-3881/ab6d6e
  * Yang et al. (2020) Yang, F., Chary, R.-R., & Liu, J.-F. 2020, arXiv e-prints, arXiv:2012.08744.  https://arxiv.org/abs/2012.08744