##  5 Photometric analysis

###  5.1 Transit analysis

####  5.1.1 Derivation of the system parameters

We performed a joint fit of all the transit light curves described in Sect. 2
using the most recent version of the adaptive MCMC code presented in Gillon et
al. (2012, see also Gillon et al. 2014). This implementation makes uses of the
Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970) and
the Gibbs sampler (Geman & Geman, 1984). The data were modelled using the
quadratic limb-darkening transit model of Mandel & Agol (2002) multiplied by a
photometric baseline model, different for each light curve, aimed at
representing the photometric variations caused by other astrophysical,
instrumental, or environmental effects. Fitting the transit signals
simultaneously with the possibly underlying correlated noise (instead of pre-
detrending the data) ensures a better propagation of the uncertainties to the
derived system parameters of interest. For each light curve, we explored a
large range of baseline models consisting of low-order polynomials with
respect to, e.g., time, airmass, full width at half maximum of the point
spread function, target’s location on the detector, background, or any
combination of these parameters. The minimal baseline model is a simple
constant to account for any out-of-transit flux offset. We show in Table
LABEL:tab:baselines the baseline model selected for each light curve based on
the Bayesian Information Criterion (BIC, Schwarz 1978). This table also gives
for each light curve the two scaling factors, $\beta_{w}$ and $\beta_{r}$,
that were applied to the photometric error bars to account respectively for
over- or under-estimated white noise and the presence of residual correlated
(red) noise (see Gillon et al. 2012 for details).

The transit model parameters sampled by the MCMC for each of the two planets
were: the log of the orbital period (log $P$), the mid-transit time ($T_{0}$),
the transit depth ($\mathrm{d}F=R_{\mathrm{p}}^{2}/R_{\star}^{2}$ where
$R_{\rm{p}}$ is the radius of the planet and $R_{\star}$ the stellar radius),
the cosine of the orbital inclination (cos $i_{\rm{p}}$), and the two
parameters $\sqrt{e_{\rm{p}}}$ cos $\omega_{\rm{p}}$ and $\sqrt{e_{\rm{p}}}$
sin $\omega_{\rm{p}}$ (with $e_{\rm{p}}$ the eccentricity and
$\omega_{\rm{p}}$ the argument of periastron). We also fitted the log of the
stellar density (log $\rho_{\star}$), the log of the stellar mass (log
$M_{\star}$), and the effective temperature ($T_{\rm{eff}}$). Finally, we
fitted for each bandpass the combinations $q_{1}=(u_{1}+u_{2})^{2}$ and
$q_{2}=0.5\,u_{1}(u_{1}+u_{2})^{-1}$ of the quadratic limb-darkening
coefficients ($u_{1}$ and $u_{2}$), following the triangular sampling scheme
advocated by Kipping (2013). The coefficients of the baseline models were not
sampled by the MCMC, but determined from the residuals at each step of the
procedure using a singular value decomposition method (Press et al., 1992).

The priors used in our analysis are listed in the last column of Table 3. We
assumed normal prior distributions for $M_{\star}$ and the luminosity
$L_{\star}$ based on the values derived in Sect. 3.2 (Table 2). Normal priors
for the quadratic limb-darkening coefficients ($u_{1}$ and $u_{2}$) were
computed for each bandpass using the PyLDTk code (Parviainen & Aigrain, 2015)
and the library of PHOENIX high-resolution synthetic spectra of Husser et al.
(2013). We increased the widths of these computed priors by a conservative
factor of ten to account for model-dependent uncertainties.

We performed two different analyses, assuming either eccentric or circular
orbits for the two planets (setting in this case their respective
$\sqrt{e_{\rm{p}}}$ cos $\omega_{\rm{p}}$ and $\sqrt{e_{\rm{p}}}$ sin
$\omega_{\rm{p}}$ values to zero). For each analysis, we ran two chains of 250
000 steps (including 20% burn-in) and checked their convergence by using the
statistical test of Gelman & Rubin (1992), ensuring that the test values for
all sampled parameters were ¡1.01. The physical parameters of the system were
deduced from the above transit model parameters at each step of the MCMC. For
each planet, the semi-major axis in stellar radii ($a_{\rm{p}}/R_{\star}$) was
computed from $\rho_{\star}$ and $P$ using Kepler’s third law (Seager &
Mallén-Ornelas, 2003). $R_{\star}$ was deduced from $\rho_{\star}$ and
$M_{\star}$. $L_{\star}$ was computed from $R_{\star}$ and $T_{\rm{eff}}$. For
each planet, $R_{\rm{p}}$, $i_{\rm{p}}$, $a_{\rm{p}}$, $e_{\rm{p}}$,
$\omega_{\rm{p}}$, the stellar irradiation ($S_{\rm{p}}$), and equilibrium
temperature ($T_{\rm{eq}}$) were then easily obtained from the above
parameters. We did not find evidence for eccentric orbits for any of the two
planets based on the transit photometry, so we adopted the results of the
circular analysis as our nominal solution (this is also consistent with
simulations of the system’s tidal evolution, see Sect. 7.3).

The posterior distribution functions for the main fitted transit parameters
are displayed in Fig. 19. Table 3 gives the medians and 1-$\sigma$ credible
intervals of the posterior distribution functions obtained for the system
parameters. The corresponding best-fit transit models are shown in Fig. 9
(planet b) and 10 (planet c), together with the detrended phase-folded
photometry. Most notably, we find that the two planets have similar radii of
$1.320_{-0.027}^{+0.053}$ and $1.367_{-0.039}^{+0.055}$ $R_{\oplus}$,
respectively, and that they orbit rather close to the second-order 3:1 mean
motion resonance, but not exactly within
($P_{\mathrm{c}}/P_{\mathrm{b}}=3.098$). We also note that the stellar density
$\rho_{\star}=32.9_{-3.0}^{+1.5}\>\rho_{\odot}$ derived from our global
transit fit is in very good agreement with the value of
$31.3_{-4.7}^{+5.9}\>\rho_{\odot}$ expected based on the stellar radius from
our SED fit and the stellar mass from our evolutionary modelling (see Table
2).

####  5.1.2 Check for transit chromaticity

To assess the transit chromaticity, we repeated the same global analysis as in
Sect. 5.1.1, but allowing this time some transit depth variations between the
different bandpasses. The transit depths found for the two planets in each
observed bandpass are given in Table 4. For each planet, the transit depths
are all in agreement (at the $\sim$1.1$\sigma$ level) and do not show any
chromatic dependence.

Table 4: Transit depths returned by our global data analysis for which we allowed transit depth variations between the different bandpasses (see Sect. 5.1.2). Bandpass | Transit depth | Transit depth
---|---|---
| planet b (ppm) | planet c (ppm)
$g^{\prime}$ | $6634_{-2170}^{+2190}$ | $-$
$r^{\prime}$ | $6100_{-604}^{+628}$ | $4657_{-2180}^{+2400}$
$i^{\prime}$ | $5971_{-288}^{+296}$ | $6227_{-887}^{+870}$
$z^{\prime}$ | $6337_{-265}^{+266}$ | $7203_{-500}^{+509}$
$I+z$ | $6353_{-221}^{+227}$ | $6545_{-311}^{+323}$
blue-blocking | $6267_{-1620}^{+1650}$ | $-$
0.85–1.55 $\mu$m | $5771\pm 1250$ | $-$
TESS | $5937_{-857}^{+850}$ | $6221_{-1920}^{+1880}$
Table 5: Transit timings returned by our global data analysis for which we allowed TTVs from the linear transit ephemerides defined by the orbital period $P$ and mid-transit time $T_{0}$ given in Table 3 (see Sect. 5.1.3). Epoch | Transit timing | TTV
---|---|---
| ($\mathrm{BJD_{TDB}}-2\,450\,000$) | (min)
Planet b
-4 | $9436.90651_{-0.00034}^{+0.00033}$ | $-0.37_{-0.50}^{+0.48}$
0 | $9447.82619_{-0.00036}^{+0.00035}$ | $-0.26_{-0.51}^{+0.50}$
4 | $9458.74582_{-0.00054}^{+0.00049}$ | $-0.23_{-0.77}^{+0.70}$
11 | $9477.85672_{-0.00075}^{+0.00065}$ | $2.05_{-1.07}^{+0.93}$
14 | $9486.04532_{-0.00036}^{+0.00038}$ | $0.46_{-0.52}^{+0.54}$
15 | $9488.7757_{-0.0020}^{+0.0022}$ | $1.18_{-2.91}^{+3.19}$
23 | $9510.61356_{-0.00100}^{+0.00082}$ | $-0.81_{-1.44}^{+1.17}$
25 | $9516.07459_{-0.00056}^{+0.00054}$ | $0.96_{-0.81}^{+0.78}$
29 | $9526.99350\pm 0.00027$ | $-0.06\pm 0.39$
30 | $9529.72301_{-0.00066}^{+0.00072}$ | $-0.62_{-0.95}^{+1.03}$
37 | $9548.83271_{-0.00089}^{+0.00078}$ | $-0.07_{-1.28}^{+1.12}$
38 | $9551.56258_{-0.00046}^{+0.00047}$ | $-0.11_{-0.66}^{+0.67}$
49 | $9581.5911_{-0.0020}^{+0.0015}$ | $-0.63_{-2.92}^{+2.24}$
Planet c
-2 | $9503.79132_{-0.00063}^{+0.00067}$ | $0.40_{-0.91}^{+0.96}$
0 | $9520.70594_{-0.00059}^{+0.00057}$ | $-0.04_{-0.85}^{+0.82}$
2 | $9537.6207\pm 0.0011$ | $-0.34_{-1.53}^{+1.52}$
9 | $9596.82308_{-0.00052}^{+0.00053}$ | $-0.07_{-0.75}^{+0.76}$
Figure 11: TTVs of the two planets measured from the ground-based transit
photometry (see Sect. 5.1.3).

####  5.1.3 Search for transit timing variations

As mentioned in Sect. 5.1.1, the system is within 10$\%$ of the second-order
3:1 mean motion resonance. To explore our photometric dataset for possible
transit timing variations (TTVs), we repeated the same global analysis as in
Sect. 5.1.1, but allowing this time for each transit of the two planets a
timing offset with respect to the reference transit ephemerides defined by the
$P$ and $T_{0}$ values reported in Table 3. We limited our search for TTVs to
the ground-based transit photometry, as TESS individual transits do not have a
high enough S/N to allow precise measurements of their timings. Table 5
presents the timings and corresponding TTVs that we obtained for the 13
transits of planet b and 4 transits of planet c that were observed from the
ground. Fig. 11 shows the TTVs as a function of time. We do not detect any
significant TTVs in the current dataset, nor any trend that could suggest the
existence of additional planets in the system.

###  5.2 Searches for additional transiting planets and detection limits

####  5.2.1 TESS photometry

As mentioned in Sect. 2.1, LP 890-9 was observed by TESS in four sectors.
While the TESS Science Office issued an alert for TOI-4306.01 (LP 890-9 b), it
did not report any other transit signal that could be related to LP 890-9 c or
other planets in the system. However, due to the detection threshold
(MES=7.1$\sigma$) to issue an alert by the automatic SPOC and Quick Look (QLP)
pipelines, some lower S/N transit signals may have gone unnoticed. We thus
performed our own search for transit signals in the TESS data using our custom
SHERLOCK999The SHERLOCK (Searching for Hints of Exoplanets fRom Lightcurves Of
spaCe-based seeKers) code is fully available on GitHub:
https://github.com/franpoz/SHERLOCK pipeline (Pozuelos et al., 2020; Demory et
al., 2020).

SHERLOCK is an end-to-end pipeline that combines six different modules that
allow exploring the TESS data fast and robustly. These six modules consist of
(1) downloading and preparing the light curves from their online repositories,
(2) searching for planetary candidates, (3) performing a semi-automatic
vetting of the interesting signals, (4) computing a statistical validation,
(5) modelling the signals to refine their ephemerides, and (6) computing
observational windows from ground-based observatories to trigger a follow-up
campaign. By default, SHERLOCK works with the PDCSAP light curve and applies a
multi-detrend approach employing the bi-weight algorithm provided in the wōtan
package (Hippke et al., 2019) to optimise the transit search. This strategy
allows the user to maximise the S/N and the signal detection efficiency (SDE)
during the transit search, which is performed over the nominal PDCSAP light
curve, jointly with the new detrended light curves, employing the transit
least squares (TLS) package (Hippke & Heller, 2019). TLS is optimised for
detecting shallow periodic transits using an analytical transit model based on
stellar parameters. The transit search is carried out in a loop; once a signal
is found, it is stored and masked, and then the search keeps running until no
more signals above a user-defined S/N threshold are found in the dataset. Each
of these search-find-mask actions is called a ‘run’. Our experience points out
that results found beyond five runs are less reliable due to the accumulated
gaps in a given light curve after many mask-and-run iterations. Thus, in our
search, we allowed a maximum of five runs. In addition, SHERLOCK allows one to
apply a preliminary Savitzky-Golay filter (Savitzky & Golay, 1964), which
enhances the detection of shallow transits with the associated risk of
obtaining more false-positive detections that the user needs to vet carefully.
This filter is very useful when dealing with threshold-crossing events with
low S/N like TOI-4306.01, which was found with a MES value of 8.3 by SPOC.

We first tried to recover this candidate by running two suites of searches,
one with a pure TLS search and the other with the Savitzky-Golay filter
applied previously. In both cases, we successfully recovered the candidate
issued by TESS in the first run, with signal-to-noise ratios ranging from 3.2
to 3.9 and 6.1 to 7.1, respectively. The very low signal-to-noise ratios
reported by the pure TLS search motivated us to keep using only the Savitzky-
Golay filter strategy. Then, we performed two independent transit searches for
extra planets by considering the four sectors simultaneously. (1) We focussed
our search on orbital periods ranging from 0.5 to 40 days where a minimum of
two transits was required to claim a detection. (2) We then focussed on longer
orbital periods, ranging from 40 to 80 days, where single events could be
recovered. None of these strategies yielded positive results, with all the
signals found being attributable to systematics, noise, or variability.

Following Wells et al. (2021), several scenarios may explain the lack of extra
signals in TESS data: (1) no other planets exist in the system; (2) if they do
exist, they do not transit; or (3) they do exist, but the photometric
precision of the data is not good enough to detect them, or they have periods
longer than the ones explored in this study. If scenario (2) or (3) is true,
extra planets might be detected by means of RV follow-up as discussed in Sect.
8.2. On the other hand, we know from the ground-based photometric follow-up
that there is at least a second transiting planet in the system, LP 890-9 c.
The fact that it was undetected by both the SPOC and SHERLOCK transit searches
suggests that we are limited in this case by its detectability in the TESS
data (scenario 3).

To explore the detection limits of the TESS data, we used our MATRIX ToolKit
101010The MATRIX ToolKit (Multi-phAse Transits Recovery from Injected
eXoplanets ToolKit) code is available on GitHub:
https://github.com/PlanetHunters/tkmatrix (Dévora-Pajares & Pozuelos, 2022),
which was used successfully in a number of previous studies (see, e.g., Wells
et al., 2021; Schanche et al., 2022). MATRIX performs transit injection-and-
recovery experiments by inserting synthetic planets over the PDCSAP light
curve, combining all the sectors available. The user defines the ranges in the
$R_{\mathrm{p}}$–$P$ parameter space to be examined. In our case, we explored
the ranges of 0.5 to 3.0 $R_{\oplus}$ with steps of 0.28 $R_{\oplus}$, and
0.5–30.0 days with steps of 0.33 days. In addition, each combination of
$R_{\mathrm{p}}$–$P$ was explored using four different phases, that is,
different values of $T_{0}$. Hence, we explored a total of 3600 scenarios. For
simplicity, the synthetic planets were injected assuming their impact
parameters and eccentricities equal zero. Moreover, we detrended the light
curves using a bi-weight filter with a window size of 1 day, which was found
to be the optimal value during the SHERLOCK executions, and masked the
transits corresponding to LP 890-9 b. We considered a synthetic planet to be
recovered when its epoch was detected with 1 hour accuracy and the recovered
period was within 5 % of the injected period ($P_{\mathrm{injected}}$). In
addition, for each period found that did not match the injected period, MATRIX
checked if it could correspond to the first harmonic (2$\times
P_{\mathrm{injected}}$) or sub-harmonic ($P_{\mathrm{injected}}/2$). If this
was the case, then the signal was also considered as recovered. It is worth
noting that since we injected the synthetic signals in the PDCSAP light curve,
then these signals were not affected by the PDCSAP systematic corrections. In
this respect, the detection limits that we find here should be considered as
an optimistic scenario (see, e.g., Pozuelos et al., 2020; Eisner et al.,
2020). As demonstrated previously during the SHERLOCK executions, using a
preliminary Savitzky-Golay filter before the transit search yields higher S/N
detections. We thus used the same strategy during the injection-and-recovery
experiments. Finally, for each scenario, we performed up to three runs (as
defined above for SHERLOCK) to try to recover the injected signal. Under these
considerations, our transit search performed above via SHERLOCK, with a larger
number of runs (five) and the multi-detrend approach, can be considered more
efficient than the injection-and-recovery map obtained by MATRIX, which still
offers a reasonable estimation of the detection limits.

The results, shown in Fig. 12 (upper panel), allow us to rule out planets
whose sizes are $>$1.8 $R_{\oplus}$ and orbital periods $<$14 days, with
recovery rates ranging from 100 to $\sim$80 $\%$. The same planetary sizes but
with orbital periods between 14 and 24 days have recovery rates ranging from
$\sim$80 to 40 $\%$. Beyond an orbital period of 24 days, the sensitivity
decreases to $\sim$20–0 $\%$. Smaller planets with sizes between 1.0 and 1.8
$R_{\oplus}$ have good recovery rates, between 100 and $\sim$60 $\%$, for
short orbital periods $<$3.2 days. LP 890-9 b lies in this region, with a
recovery rate of $\sim$60 $\%$. On the other hand, when exploring longer
orbital periods $>$3.2 days, the detectability of such small planets with
sizes between 1.0 and 1.8 $R_{\oplus}$ decreases to $\sim$50–0 $\%$. In this
region lies LP 890-9 c, with a recovery rate of only $\sim$6 $\%$. This result
is consistent with the non-detection of LP 890-9 c in the TESS data during our
transit search with SHERLOCK. For the whole range of periods explored, planets
with $R_{\mathrm{p}}\leq$1.0 $R_{\oplus}$ can not be found, which hints at a
clear limitation of the TESS data to discover transiting (sub-)Earth-sized
planets in this system.

Figure 12: Results from the transit injection-and-recovery tests performed on
the TESS (upper panel) and SSO (lower panel) data to assess the detectability
of transiting planets in the LP 890-9 system. Upper panel: We explored a total
of 3600 different scenarios. Each pixel evaluated 32 scenarios, that is, 32
light curves with injected planets having different $P$, $R_{\mathrm{p}}$, and
$T_{0}$. Larger recovery rates are presented in yellow and green, while lower
recovery rates are shown in blue and darker hues. LP 890-9 b (recovery rate
$\sim$60 $\%$) and c (recovery rate $\sim$6 $\%$) are displayed as red and
blue stars, respectively (see Sect. 5.2.1 for details). Lower panel: We
injected 6000 artificial planets into the combined $I+z$-filter SSO light
curve. On average, each box in this plot contains 50 injection scenarios with
differing $P$, $R_{\rm{p}}$, $i_{\rm{p}}$, and $\phi$. Similarly to the upper
panel, we include LP 890-9 b (recovery rate of 100 $\%$) and c (recovery rate
of 96 $\%$) as red and blue stars, respectively (see Sect. 5.2.2 for details).

####  5.2.2 SPECULOOS-South photometry

After the announcement of the TOI-4306.01 planet candidate by TESS in July
2021, we initiated a blind search for other transiting planets in the system
with SSO and detected the first transit of planet c after about 2 months of
daily monitoring. We stopped our follow-up in mid-January 2022, only when we
were confident that the empirical (optimistic) habitable zone (HZ) of the star
(Kopparapu et al. 2013, 2014, see also Sect. 8.3) had been reasonably well
explored. To quantify this, we used a metric that we introduced in Sebastian
et al. (2021), the effective phase coverage. This metric gives an estimation
of how the phase of a hypothetical planet would be covered for a range of
periods. In that regard, we computed the percentage of phase covered for each
orbital period from $P=0.1$ day to a given $P=P_{\mathrm{max}}$ and took the
effective phase coverage for periods $\leq$ $P_{\mathrm{max}}$ to be the
integral over the period range. We stopped our photometric monitoring of LP
890-9 when we reached an effective coverage of at least 80% for
$P_{\mathrm{max}}=25.823$ days, which is the orbital period associated to the
outer limit of the empirical HZ (Kopparapu et al. 2013, 2014, see also Sect.
8.3). Fig. 13 shows the evolution of the phase coverage as a function of the
orbital period of a putative transiting planet based on the SSO photometry.
The periodic drops in phase coverage that are visible for all integer periods
are expected for ground-based observations from a unique site. The effective
phase coverage for the outer edge of the empirical HZ is 86%. As a comparison,
we detected the first transit of planet c when its effective phase coverage
was only 57%, and observed its second transit when the effective coverage was
78%.

Figure 13: Phase coverage as a function of the orbital period using all
photometric data obtained with SSO for LP 890-9. The effective coverage of 86%
is the integral of the black curve. The orbital periods of planets b and c are
shown as orange and red lines, respectively. The figure also shows the
empirical (optimistic) HZ in green, i.e. the region between the ‘Recent Venus’
and ‘Early Mars’ limits (Kopparapu et al. 2013, 2014, and Sect. 8.3). Figure
14: Global $I+z$-filter light curves for the SSO/Io (red) and SSO/Europa
(blue) observations. The black points correspond to 20-min bins.

We ran an automatic transit search through the SSO $I+z$ photometry to check
for potential additional transit signals that may have been missed by our
visual inspection of the data. Transit-finding tools often require a rather
‘flat’ light curve to work optimally. Therefore, we first stitched together
each night’s normalised $I+z$ differential light curve, as these light curves
are optimised for the atmospheric conditions of the night, combining light
curves from all SSO telescopes. This ‘combined’ light curve was then cleaned
for bad weather, as described in Murray et al. (2020), and frames were removed
where the atmospheric conditions would significantly affect the photometry. We
also removed all observations between 28 December 2021 and 3 January 2022
inclusive, as our autoguiding software, donuts (McCormac et al., 2013), did
not run on these nights and the corresponding light curves were significantly
more noisy. Then, to model remaining photometric variations (caused by
intrinsic stellar variability and nightly atmospheric variations) without
modelling transit features, we employed a Gaussian Process (GP) approach
(Rasmussen & Williams, 2006). We used a squared exponential kernel (Rasmussen
& Williams, 2006) implemented by the george package (Foreman-Mackey et al.,
2014; Ambikasaran et al., 2015):

$k(t_{i},t_{j})=A\exp\left(-\frac{(t_{j}-t_{i})^{2}}{l^{2}}\right)$ (1)

where $i$ and $j$ are two different data points in the light curve (taken at
times $t_{i}$ and $t_{j}$), $A$ is the amplitude of correlation, and $l$ is
the length-scale over which correlations decay. To apply the GP regression,
the light curve was binned every 10 minutes and sigma-clipped twice to lower
the significance of any infrequent, short features. In essence, binning and
sigma-clipping have the same effect; they remove short-duration structures to
prevent over-fitting with the GP. The first sigma-clip was a basic 3$\sigma$
clip of the entire light curve, to remove the most extreme outliers. We also
included a second, ‘running’ sigma-clip to account for any photometric
variability. Here, we binned the light curve into hours, longer than a typical
transit duration, to produce a running median (a series of points representing
the median of each hour bin) and a running standard deviation (defined
similarly). We then computed the median running standard deviation, $\sigma$,
and clipped the data points more than 3$\sigma$ above or below the running
median. Even though this target does not demonstrate significant photometric
variability, this method can also be used to sigma-clip rapidly-varying stars.

For a transit-finding algorithm, we used the astropy (Astropy Collaboration et
al., 2013, 2018) implementation of Box Least Squares, or BLS (Kovács et al.,
2002). We considered using Transit Least Squares, or TLS (Hippke & Heller,
2019); however, for very red stars, such as those in SPECULOOS’s target list,
the detection efficiencies of BLS and TLS converge while TLS is more demanding
in terms of computing time and power (Hippke & Heller, 2019). We limited BLS
to search for 10,000 periods between 0.5 and 30 days, with a transit duration
between 0.02–0.09 days. The highest significance period is flagged, and the
corresponding transits masked before re-running the BLS transit search on the
light curve. This automatic, iterative process is run until the S/N of the
transits is less than 5. For LP 890-9, we first successfully detected planet b
with a period of 2.730 days. Next, planet c was detected with a period of
8.457 days. All $I+z$ transits from planets b and c in Table 1 were recovered.
There were no additional, convincing detections with S/N$\geq$5.

To assess the detection efficiency of our transit-search pipeline and the
detection limits of the data themselves, we performed injection and recovery
tests on the SSO $I+z$ light curve. We generated transits for 6000 artificial
planets using PyTransit (Parviainen, 2015), a transit light curve modelling
package that implements the Mandel & Agol (2002) model. We used the limb-
darkening coefficients $u_{1}=0.384\pm 0.028$ and $u_{2}=0.303\pm 0.089$ for
the $I+z$ filter from Sect. 5.1. For each of the planets injected, their
parameters (radius $R_{\rm{p}}$, orbital period $P$, and inclination
$i_{\rm{p}}$) were drawn from the following distributions:

$\begin{split}\frac{R_{\rm{p}}}{R_{\oplus}}&\sim\mathcal{U}(0.5,3.0)\\\
\cos{i_{\rm{p}}}&\sim\mathcal{U}(\cos{i_{\rm{min}}},\cos{i_{\rm{max}}})\\\
\frac{P}{\rm{days}}&\sim\mathcal{U}(0.5,30.0)\end{split}$ (2)

where $\mathcal{U}(a,b)$ represents a uniform distribution between $a$ and
$b$. $i_{\rm{min}}$ and $i_{\rm{max}}$ are the minimum and maximum
inclinations for a transiting planet. The host mass and radius were assumed
constant, using the values in Table 2. However, the inclination limits depend
on the orbital period; therefore, when drawing the planetary parameters, we
drew individual inclinations from a range set by each period. We only
considered circular orbits (with $e$=0) as we do not know the underlying
eccentricity distribution. Close-in planets are also expected to experience
significant circularisation of their orbits (Luger et al., 2017). The time at
which the first transit was injected was also randomly drawn from a uniform
distribution, $\phi\sim\mathcal{U}$(0,1), where $\phi$ is the phase of the
period, such that the first transit was injected at $\phi P$ from the start of
observations.

We injected each of the 6000 artificial planets into the target’s differential
light curve and ran the transit-search pipeline. The planets were injected
before the GP-detrending to allow for the possibility that the GP could over-
model and remove or warp the transit signals. If we were to inject planets
after GP-detrending, this could artificially inflate the recovery results, as
it is easier to detect a transit in a ‘flat’ light curve than in a light curve
exhibiting time-dependent structure. From BLS, we only considered the most
likely transit and did not look at other peaks in the periodogram. A planet
was successfully recovered if at least one transit was detected by BLS above a
S/N threshold of 5, with an epoch detected within an hour of the injected
transit epoch. We did not include any requirements on the recovered period due
to large gaps in the data from the day-night cycle and bad weather conditions.
These gaps mean that we often recover period aliases. Additionally, when
testing periods up to 30 days, we have a high chance of only detecting a
single transit with no meaningful period information. As we only require one
‘real’ transit for recovery, it is possible that we can recover a planet even
if we miss some of its transits, or incorrectly detect noise as transits. The
lack of a period condition on planet recovery may thus result in an
overestimation of our recovery fractions. Therefore, this injection-recovery
framework is most useful for detecting individual transit signals (including
single transits) which would then be followed up with manual vetting.

For all synthetic planets with radii $R_{\rm{p}}$ = 0.5–3.0 $R_{\oplus}$ and
periods $P$ = 0.5–30 days, we recovered $\sim$72%. The lower panel of Fig. 12
shows how the recovery rates vary in the $R_{\rm{p}}$ – $P$ parameter space.
Comparing these results with the ones obtained in Sect. 5.2.1 for the TESS
data (upper panel of Fig. 12), we see that the SSO data are sensitive to
smaller planets and/or longer orbital periods. SSO has a high detection
potential for short-period ($P\,\leq\,10$ days) planets in the Earth to super-
Earth regime ($R_{\rm{p}}\,\geq\,1\,R_{\oplus}$) with recovery rates ranging
from 86 to 100%, and an average of $\sim$97%. This is consistent with the
detection of both the b and c planets in the SSO data, with respective
recovery rates of 100 and 96 %. The sensitivity is still reasonable for
smaller planets with 0.75 $R_{\oplus}\leq R_{\rm{p}}<1\,R_{\oplus}$ in the
same period range ($P$ $<$ 10 days), with an average of $\sim$57%. However, it
drops to $\sim$15% for even smaller planets with 0.5
$R_{\oplus}\leq\,R_{\rm{p}}<\,0.75\,R_{\oplus}$. Due to the long observation
span of LP 890-9 with the SSO and our decision not to include a period
criterion in our recovery, we obtain good recovery rates (above $\sim$60%) for
Earth-sized planets and super-Earths ($R_{\rm{p}}\,\geq\,1\,R_{\oplus}$) with
periods longer than 10 days, up to 30 days. For the same range of periods, the
sensitivity drops to $\sim$25–50% for planets with 0.75 $R_{\oplus}\leq
R_{\rm{p}}<1\,R_{\oplus}$ and below 15% for the smallest planets with 0.5
$R_{\oplus}\leq\,R_{\rm{p}}<\,0.75\,R_{\oplus}$. Based on these results, we
conclude that it is unlikely there are any additional transiting short-period
(super-)Earth-sized planets hidden in this system. However, it is still
possible that there could be more transiting bodies, either too small or too
far away from their host to be detected in the SSO data.

Figure 15: Lomb-Scargle periodograms for SSO/Io (top, red) and SSO/Europa
(bottom, blue). The five highest peaks are marked by black crosses in both
cases. For SSO/Europa, the most significant peak is around 50–55 days and the
next highest four are at $\sim$1 day and were determined to be aliases. The
results for SSO/Io are less clear due to large gaps and much less data. The
most significant peaks are found at 23.0, 54.0, and 4.3 days. The peak around
54 days could correspond to the one detected for SSO/Europa, but the other two
are not present in SSO/Europa’s periodogram.

###  5.3 Stellar variability

We searched for flares and rotational variability in the SSO $I+z$-filter
photometry (which has a higher precision than the TESS photometry for this
late-type target) using the techniques described in Murray et al. (2022). The
light curves used in this section were cleaned for bad weather and poor
photometric conditions, as described in Sect. 5.2.2. However, unlike in Sect.
5.2.2, we did not use the nightly normalised light curves. Instead, we used
the raw flux values from aperture photometry across all nights observed with
the same telescope as one time series, and performed our differential
photometry process on the entire time series, to obtain what we call a
‘global’ light curve (see also Murray et al. 2020). Doing so allows to
preserve the flux relationships between nights and study long-term photometric
trends, but limits our ability to optimise the photometry night-by-night
(e.g., aperture size, comparison stars). The global light curves for SSO/Io
and SSO/Europa are presented in Fig. 14.

The search for flares was done in two parts: a simple, automatic flare-
detection algorithm was run first to extract all the flare candidates (based
on the method described in Lienhard et al. 2020), followed by a manual vetting
process to confirm them. We detected no conclusive flares (flare structures
that could not be attributed to atmospheric effects) in any of the light
curves using this method.

To search for a rotation period, we applied the Lomb-Scargle periodogram
analysis (Lomb, 1976; Scargle, 1982) to the global light curves, binned every
20 minutes. The SSO light curves are not uniformly sampled as there are gaps,
not only from the day/night cycle, but also from bad weather and changes to
SSO’s observation strategy. Therefore, careful treatment of the resulting
periodograms is essential to remove aliases. We searched for periods between
the Nyquist limit (twice the bin size) and the entire observation window. We
removed all peaks in the periodogram with a false alarm probability below
3$\sigma$ (0.0027). We visually inspected the periodograms of LP 890-9, shown
in Fig. 15, and the corresponding phase-folded light curves for all possible
rotation periods (i.e., all significant peaks in this periodogram). In
addition, by comparing with periodograms of the time stamps, airmass, and
FWHM, we eliminated signals which arose from the non-uniform sampling and from
ground-based systematics. The most likely rotation period detected for the
SSO/Europa light curve (individually and combined with SSO/Io) was found to be
around 50–55 days. However, as this value is approximately a third of the
entire observation span of SSO/Europa, more observations would be required to
confirm. The Lomb-Scargle results for SSO/Io were less clear, with the most
significant peaks at 23.0, 54.0, and 4.3 days. The signal around 54 days could
correspond to the one detected for SSO/Europa, but the other two are not
present in SSO/Europa’s periodogram. These results were consistent regardless
of whether the transits of planets b and c were masked. However, when applying
the Lomb-Scargle analysis to the light curves with a longer binning of 60
minutes, we found the same result for SSO/Europa, but no significant peaks in
the periodogram for SSO/Io. Additionally, there were unexplained jumps in the
light curve on the nights of 3 December 2021 and 15 January 2022 that did not
correlate with FWHM, sky background level, target’s position on the CCD,
airmass, nor with the photometric aperture or any instrumental parameter (CCD
temperature, exposure time, …). If these jumps are due to unknown systematics,
then they will affect the reliability of our rotation analysis.

Figure 16: Subaru/IRD radial velocities obtained in the $YJ$\- (grey squares)
and $H$\- (black circles) bands. The data are plotted with both their original
error bars (thick solid lines) and the error bars enlarged by the best-fit
jitter terms (thinner transparent lines). The top panel displays the data as a
function of time. The solid red line shows the 2-planet model generated from
the posterior median, while semi-transparent red lines show 25 models randomly
drawn from the posteriors. The residuals around the median solution are shown
in the second panel. The two lower panels show the data folded on the orbital
periods of planets b and c, respectively, after removing the RV component from
the other planet.

##  6 Radial velocity analysis

We analysed the Subaru/IRD data using the juliet package (Espinoza et al.,
2019), which is built over radvel (Fulton et al., 2018) for the modelling of
radial velocities and the dynesty (Speagle, 2020) dynamic nested sampling
algorithm for estimating Bayesian posteriors and evidences. We modelled the
RVs using a sum of two Keplerians (one for each planet) and assumed circular
orbits. We imposed normal priors on the orbital period $P$ and mid-transit
time $T_{0}$ of each planet based on the results of our global transit
analysis reported in Table 3. We assumed wide uniform priors between 0 and 500
m/s for the RV semi-amplitudes $K$. We treated the $YJ$-band and $H$-band RVs
as independent datasets, thus fitting for each of them a different zero point
(systemic RV) and extra jitter term. For each dataset, this jitter is added
quadratically to the measurement errors of the data points, to account for any
underestimation of the uncertainties or any excess noise not captured by the
model. Fig. 16 shows the results of our RV fit. The data are plotted with both
their original error bars (thick solid lines) and the error bars enlarged by
the best-fit jitter terms (thinner transparent lines). These jitters are large
($33.0_{-7.4}^{+10.5}$ m/s for the $YJ$-band and $19.6_{-4.3}^{+5.7}$ m/s for
the $H$-band) and the datasets are rather limited, so that the RV semi-
amplitudes of the two planets are poorly constrained and we were only able to
derive upper limits: $K<25.1$ m/s for planet b and $K<33.0$ m/s for planet c
(2$\sigma$ upper limits). The derived upper limits on the planetary masses
(using the orbital inclinations from our transit analysis, Table 3) are
$M_{\mathrm{p}}<13.2\>M_{\oplus}$ for planet b and
$M_{\mathrm{p}}<25.3\>M_{\oplus}$ for planet c at 2$\sigma$. Both are well
within the planetary regime.

The large extra jitters returned by our 2-planet fit suggest that there is
some variability in the data that is not captured by the model. One possible
explanation is that there may be an additional planet in the system which
produces a detectable RV signal but was not detected in our photometry (either
because it is not transiting or because it is beyond the detection limits of
our photometric dataset). We thus tried adding a third planet as a third
Keplerian in the model, using a wide uniform prior between 1 and 125 days (the
time baseline of the RV data) for its orbital period. Comparing the Bayesian
evidences $Z$ of the two models, we found that the 3-planet model is only
weakly disfavoured compared to the 2-planet one ($\Delta$ ln $Z$ = $-$0.65).
The orbital period found for the third Keplerian is $64_{-9}^{+11}$ days, with
a 2$\sigma$ upper limit on its semi-amplitude of 49 m/s. Instead of adding a
third Keplerian in the model, we also tried including a linear trend that
could be created by a third planet with an orbital period longer than the time
baseline of the RV data. We found this model to be marginally disfavoured
($\Delta$ ln $Z$ = $-$3.72) compared to the 2-planet model. Since the 2-planet
model is the simplest of the three models we tested and it has the highest
Bayesian evidence, it appears to be the best model given the data at hand.
However, additional planets in the system might be revealed with more RV data.
The variability seen in the current RV dataset may also be (at least partly)
related to stellar activity. In this context, it is worth noting that the
period of $64_{-9}^{+11}$ days found for the third Keplerian in the above
3-planet fit is within the range of the possible rotation period of 50–55 days
suggested by the SSO photometry (Sect. 5.3). Finally, it is also possible that
the RV measurements are impacted by systematic effects. Indeed, large and
currently unexplained RV variabilities have been found in Subaru/IRD data
obtained for other faint and apparently quiet M dwarfs (see, e.g., Morello et
al. 2022; Mori et al. 2022).

Using the probabilistic mass–radius relationship of Chen & Kipping (2017), we
find predicted masses of $2.3_{-0.7}^{+1.7}$ and $2.5_{-0.8}^{+1.8}$
$M_{\oplus}$ for LP 890-9 b and c, respectively. Combined with the orbital
parameters of the two planets, these masses would correspond to respective
radial-velocity semi-amplitudes of $4.4_{-1.4}^{+3.2}$ and $3.3_{-1.1}^{+2.4}$
m/s, significantly smaller than the measurement errors and variability of the
current Subaru/IRD data. Additional RV measurements, possibly obtained with
other available or next-generation instruments, are required to improve the
preliminary mass constraints derived here and investigate the possible
presence of other planets in the system. We discuss some prospects for RV
follow-up in Sect. 8.2.

##  7 Dynamical analysis

###  7.1 Expected TTVs

As described in Sect. 5.1.3, we do not detect any significant TTVs in the
current photometric dataset, although the two planets are close to the second-
order 3:1 mean motion resonance (but not exactly within:
$P_{\mathrm{c}}/P_{\mathrm{b}}=3.098$). In general, the lack of observed TTVs
might be caused by either their low amplitude or the poor baseline coverage
provided by the observational data. In our case, the ground-based observations
that we used to explore TTVs were obtained over $\sim$150 days, which is much
longer than the super-period of 2.64 days (Lithwick et al., 2012) for LP 890-9
b and c. Hence, our ground-based data should show any TTVs if they were large
enough. To test the expected TTV amplitudes in the system, we used the
TTVFast2Furious package (Hadden, 2019) following the strategy presented by
Cloutier et al. (2020). We ran 103 realisations with orbital periods and times
of mid-transit sampled from the posterior distributions given in Table 3. The
planetary masses were drawn from their predictions from the probabilistic
mass–radius relationship of Chen & Kipping (2017) (see Sect. 6) and the
argument of periastron following a random distribution from 0 to 2$\pi$. We
found maximum TTV values for LP 890-9 b and c of only 1.5 and 1.0 s,
respectively. As a comparison, the typical timing uncertainties returned by
our global transit analysis are on the order of 1 minute (see Table 5). Based
on the low expected TTV amplitudes, we conclude that the LP 890-9 system is
not suitable for a determination of the planetary masses (and dynamical
studies in general) using the TTV method. This is in contrast with the
TRAPPIST-1 system, where the outer five planets are close to first-
order111111The lower the order of the mean motion resonance, the larger the
induced TTV amplitude. resonances with adjacent planets, resulting in large-
amplitude TTVs (from $\sim$5 to 100 minutes) – an important feature which made
it possible to derive precise planetary masses and eccentricities (Agol et
al., 2021).

###  7.2 Stability

From the current dataset, we concluded that the LP 890-9 system is composed of
at least two planets orbiting in nearly circular orbits and with similar
planetary sizes. To assess the global stability of the system, we used the
Mean Exponential Growth factor of Nearby Orbits (MEGNO; Cincotta & Simó, 1999,
2000; Cincotta et al., 2003), $Y(t)$ parameter. This chaos index is widely
used to explore the stability of extrasolar planetary systems (e.g., Hinse et
al., 2015; Jenkins et al., 2019; Delrez et al., 2021). In short, its time-
averaged mean value, $\langle Y(t)\rangle$, amplifies any stochastic
behaviour, which can be used to distinguish between quasi-periodic orbits, if
$\langle Y(t\rightarrow\infty)\rangle\rightarrow$ 2, or chaotic trajectories,
if $\langle Y(t\rightarrow\infty)\rangle\rightarrow\infty$. We made use of the
MEGNO implementation with the N-body integrator rebound (Rein & Liu, 2012),
which employs the Wisdom-Holman WHfast code (Rein & Tamayo, 2015). In this
context, we studied a set of scenarios, where the planetary properties of LP
890-9 b and c were drawn from the posterior distributions and derived
parameters given in Table 3. In particular, semi-major axes, orbital
inclinations, and planetary masses were sampled from uniform distributions
using their nominal values up to 3$\sigma$ of their uncertainties. In
addition, the mean anomalies and arguments of periastron were sampled randomly
from 0 to 2$\pi$. We ran two suites of simulations: first, we considered
circular orbits, and secondly, slightly eccentric orbits using a uniform
distribution of 0.10$\pm$0.05 for both planets. In each experiment, we
generated 1,000 scenarios, and they were integrated over 106 orbits of the
outer planet, with a time-step of 5$\%$ of the orbital period of the inner
planet. We found that the LP 890-9 system is highly stable in both experiments
with $\langle Y(t)\rangle$ clustering between 1.999–2.001 in 99.3$\%$ and
93.8$\%$ of cases for circular and eccentric orbits, respectively.

Having shown the stability of the LP 890-9 system on its nominal configuration
and within 3$\sigma$ of the planetary properties, we considered possible
stability limitations for the planetary masses and the orbital eccentricities,
the two main parameters that have a major impact on the orbital dynamics and
that are poorly constrained by the current data. To this end, we extended our
analyses by building stability maps following Delrez et al. (2021), that is,
by exploring the $M_{\mathrm{b}}$–$M_{\mathrm{c}}$ and
$e_{\mathrm{b}}$–$e_{\mathrm{c}}$ parameter spaces. We explored planetary
masses ranging from 1 to 20 $M_{\oplus}$ and eccentricities ranging from 0.0
to 0.6 for both planets. In these sets of simulations, we took 20 values from
each range, so each stability map has $20{\times}20$ pixels. Moreover, we set
10 different initial conditions for each realisation, randomly sampling the
argument of periastron and the mean anomaly from 0 to 2$\pi$. Hence, each
stability map consisted of 4000 scenarios. We then averaged the results of the
10 different initial conditions to obtain the averaged value of $\langle
Y\rangle$ for each pixel of the stability map. The integration time and time-
step were the same as previously used. We found that the system is fully
stable over the whole range of masses explored with $\Delta\langle
Y(t)\rangle$= 2.0 $-$ $\langle Y(t)\rangle<10^{-5}$. On the other hand, we
found that the eccentricities have to satisfy the condition
$e_{\mathrm{b}}\leq 0.39-1.11\times e_{\mathrm{c}}$. This implies that the
maximum eccentricity for planet b (when planet c has a circular orbit) is
0.39, and the maximum eccentricity for planet c (when planet b has a circular
orbit) is 0.35. For eccentricity values that do not meet this condition, the
system may become unstable. By combining the results obtained along this
subsection and assuming that there are no more planets in the system, we
conclude that the dynamical architecture of the LP 890-9 system resides in a
stable set of parameters, allowing for long-term evolution.

###  7.3 Tidal evolution

The two planets of the LP 890-9 system reside in close-in orbits and,
consequently, due to the short distance from the host star, are affected by
tides. Indeed, tidal interactions between the star and planets in closely-
packed systems have an impact on planetary orbits, altering the obliquities
($\epsilon$) and rotational periods ($P_{\mathrm{rot}}$) on short time scales,
driving the planets into a particular configuration called tidal locking,
where the tidal torques fix the rotation rate to a specific frequency, and the
planetary and stellar spins are aligned. On longer time scales, the tidal
interactions can also induce the circularisation of the orbits and reduce
their semi-major axes (see, e.g., Bolmont et al., 2014; Barnes, 2017). This
scenario is interesting since a tidally-locked planet in a circular orbit
rotates synchronously (i.e., its rotational and orbital periods are equal),
meaning that the same hemisphere of the planet is always facing the star, a
configuration that can have important implications for the planet’s global
climate.

To explore the evolution of the LP 890-9 system under tidal interactions, we
employed the constant time-lag model (CTL), where a planet is modelled as a
weakly viscous fluid deformed by the gravitational interactions (see, e.g.,
Mignard, 1979; Hut, 1981; Eggleton et al., 1998; Leconte et al., 2010). We
used the implementation of CTL in the posidonius N-body code (Blanco-Cuaresma
& Bolmont, 2017). The main factor in the CTL model that controls the tidal
dissipation of a given planet is the product of two free parameters: the
degree-2 potential Love number, $k_{2}$, and the constant time-lag,
$\Delta\tau$. The $k_{2}$ accounts for the elastic properties and the self-
gravity of the planets, with values ranging from 0 to 1.5. On the other hand,
$\Delta\tau$ refers to the lag between the line connecting the two centers of
mass (star-planet) and the direction of the tidal bulges, which can span
orders of magnitudes (Barnes, 2017). Although the masses (and thus bulk
densities) of the LP 890-9 planets are currently poorly constrained (Sect. 6),
their small sizes imply that they are most likely terrestrial. Following
Turbet et al. (2018) and Bolmont et al. (2020), who studied the tidal effects
for the TRAPPIST-1 planets, we thus used Earth’s value
$k_{2,\oplus}\Delta\tau_{\oplus}=213$ s (Neron de Surgy & Laskar, 1997) as a
reference and explored dissipation factors ranging from 0.1$\times
k_{2,\oplus}\Delta\tau_{\oplus}$ (for dry rocky planets such as Mars) to
10$\times k_{2,\oplus}\Delta\tau_{\oplus}$ (for planets dominated by fluid-
solid friction), to evaluate a range of possible tidal behaviours.

In a first step, we assessed the time needed by the LP 890-9 planets to evolve
towards the tidal locking configuration. To this end, we conducted a suite of
simulations considering different scenarios combining three different initial
planetary rotational periods $P_{\mathrm{0,rot}}$ (10, 100, and 1000 h), three
different initial obliquities $\epsilon_{0}$ (15$\degr$, 45$\degr$, and
75$\degr$) and three dissipation factors (0.1$\times$, 1$\times$, and
10$\times k_{2,\oplus}\Delta\tau_{\oplus}$). We found that the slowest
evolution towards tidal locking happened for the particular combination of
initial conditions given by $\epsilon_{0}=75\degr$ and $P_{\mathrm{0,rot}}=10$
h for dry terrestrial planets (0.1$\times k_{2,\oplus}\Delta\tau_{\oplus}$).
In such a case, we found time scales of ${\sim}1.5\times 10^{4}$ and
${\sim}1.5\times 10^{6}$ yr, for planets b and c, respectively. Hence, these
values might be considered as upper limits on the time scale required for each
planet to evolve into a tidally-locked configuration. Considering the
estimated age of the system (Sect. 3.3), the two planets are likely tidally
locked.

In a second step, we evaluated the time needed for the circularisation of the
orbits due to tidal interactions. From the global transit analysis presented
in Sect. 5.1.1, we did not find any evidence for eccentric orbits. However,
the eccentricities of the planets are poorly constrained by the transit
photometry, and high-precision RV measurements are needed for a proper
estimation. On the other hand, when studying the dynamical architecture of the
LP 890-9 system in Sect. 7.2, we found that it may tolerate a certain level of
planetary eccentricities allowing for long-term dynamical evolution. Then,
assuming initial slightly eccentric orbits of 0.01 for both planets, we
performed three simulations considering dissipation factors of 0.1$\times$,
1$\times$, and 10$\times k_{2,\oplus}\Delta\tau_{\oplus}$. We found that the
eccentricities of both planets decreased to reach mean equilibrium values of
${\sim}1.8\times 10^{-4}$ and ${\sim}1.7\times 10^{-5}$ for planets b and c,
respectively. The time scales were of the order of ${\sim}10^{9}$,
${\sim}10^{8}$ and ${\sim}10^{7}$ yr, for the different dissipation factors
considered, being longer for dry rocky planets. These small residual
eccentricities, inferior to 10-3, are the result of the competition between
tidal damping and planet-planet excitations (Turbet et al., 2018). Then,
considering that the system is likely older than 2 Gyr (Sect. 3.3), planetary
orbits may have had time to evolve towards these low eccentricities due to
tidal interactions. Thus, should future observations reveal larger
eccentricities than the equilibrium values found here, this could suggest that
other unknown planets in the system might be exciting the orbits.

Figure 17: Exoplanets found around cool dwarfs with an effective temperature
lower than 3600 K using the transit (black circles) or radial velocity (grey
triangles) techniques. Data were taken from NASA Exoplanet Archive on 20 April
2022. The planets are shown as a function of their incident stellar flux and
the effective temperature of their host star. The size of the symbols is
proportional to the planetary radius, estimated from the mass using the
probabilistic mass–radius relationship of Chen & Kipping (2017) when it is
unknown (for non-transiting RV planets). LP 890-9 b and c are shown as orange
circles. The inner (Recent Venus) and outer (Early Mars) boundaries of the
empirical (optimistic) HZ are shown as solid red and blue lines, respectively
(Kopparapu et al. 2014 and Sect. 8.3). The dashed green line shows an
alternative conservative inner edge limit from 3D models (Leconte et al.,
2013). The outer edge of the HZ in 3D models agrees with incident stellar flux
for Early Mars insulation. We note that the three limits shown do vary with
the stellar effective temperature, but this is not obvious on this plot due to
the log-scale of the $x$-axis.

##  8 Discussion

###  8.1 LP 890-9 in the context of other planetary systems around cool dwarfs

Fig. 17 shows the planets that have been found so far around cool dwarfs using
the transit and radial velocity techniques, as a function of the stellar flux
they receive and the effective temperature of their host star. Aside from LP
890-9, only five planetary systems are known around hosts with
$T_{\mathrm{eff}}<3000$ K: TRAPPIST-1 (Gillon et al., 2016, 2017), Teegarden’s
star (Zechmeister et al., 2019), GJ 1061 (Dreizler et al., 2020), LP 791-18
(Crossfield et al., 2019), and LHS 1140 (Dittmann et al., 2017; Ment et al.,
2019). With an effective temperature of $2850\pm 75$ K (Sect. 3), LP 890-9 is
the second-coolest star found to host planets after TRAPPIST-1. Interestingly,
all these systems – including LP 890-9 – are multi-planet systems, in line
with the observation that compact multi-planetary systems should be a
relatively common outcome of planet formation around low-mass stars (e.g.,
Coleman et al. 2019, Miguel et al. 2020, Mulders et al. 2021).

###  8.2 Prospects for RV follow-up

As mentioned in Sect. 6, the expected RV semi-amplitudes of the two planets
based on the mass–radius relationship of Chen & Kipping (2017) are
$4.4_{-1.4}^{+3.2}$ and $3.3_{-1.1}^{+2.4}$ m/s, respectively. Given the
faintness of the star, detecting such small signals will be challenging, even
for state-of-the-art instrumentation, and will require a substantial
investment of observing time.

Besides its late spectral type, the star is too faint to recover the expected
RV semi-amplitudes with the SPIROU infrared spectrograph at the CFHT (Donati
et al., 2020) or to be measured with the upcoming NIRPS infrared spectrograph
at the 3.6 m ESO telescope (limiting magnitude $H$=9, Wildi et al. 2017). The
expected RV semi-amplitudes are also too small to be recovered with CRIRES+ at
the VLT (RV
precision121212https://www.eso.org/sci/facilities/paranal/instruments/crires/inst.html$\sim$20
m/s, Dorn et al. 2014). Simply scaling from the existing RV constraints
derived from the current Subaru/IRD data ($K_{\mathrm{b}}<25.1$ m/s and
$K_{\mathrm{c}}<33.0$ m/s) and the expected RV semi-amplitudes given above
suggests that $\sim$31 and $\sim$100 times more data would be needed to
constrain the masses of planets b and c, respectively, with that facility.

In the optical range, currently available high-resolution spectrographs, like
ESPRESSO at the VLT (Pepe et al., 2021) and MAROON-X at Gemini-North (Seifahrt
et al., 2020), would also require substantial amounts of observing time to
detect the RV signals of the planets. Even assuming an inactive and slow
rotating star (stellar jitter $\sim$1.5 m/s which is probably optimistic, see
Sect. 6), we estimate that about 60 spectra with an average exposure time of
1800s would be required with MAROON-X to recover the RV signals of both
planets with 4–5$\sigma$. With the future ANDES spectrograph at the ELT
(Marconi et al., 2021), a 7–9$\sigma$ detection should be feasible by
investing only half of this observing time.

To summarise, while a 4–5$\sigma$ detection of the RV signals would require a
large amount of observing time with current high-resolution spectrographs,
measuring the masses of the two planets with at least 9$\sigma$ accuracy will
require waiting for next-generation instruments to become available in the
next decade. Furthermore, such an RV follow-up will also allow to investigate
the origin of the large variability seen in the Subaru/IRD data (Sect. 6) and
may reveal some additional planets in the system.

Figure 18: Expected transmission signals (see Sect. 8.4 for details) of
currently known transiting terrestrial exoplanets
($R_{\rm{p}}<1.6\>R_{\oplus}$) as a function of their incident stellar flux.
Data were taken from NASA Exoplanet Archive on 29 March 2022. We only show
planets with mild incident irradiations, lower than five times that received
by the Earth ($S_{\rm{p}}<5\>S_{\oplus}$). The symbols are colour-coded as a
function of their S/N relative to TRAPPIST-1 b, obtained by scaling the signal
amplitude with the hosts’ brightness in the $J$-band and using TRAPPIST-1 b’s
S/N as a reference. The size of the symbols is proportional to the planetary
radius. Transmission signals less than 5 ppm have been removed to enhance the
readability of the figure.

###  8.3 Potential habitability

The habitable zone (HZ) is a concept that can be used to prioritise detected
rocky exoplanets for follow-up observations. It is defined as the circular
region around one or multiple stars, where liquid water is possible on the
surface of a geologically-active rocky planet (e.g., Kasting et al. 1993;
Kopparapu et al. 2013; Kaltenegger & Haghighipour 2013; Haghighipour &
Kaltenegger 2013; Kane & Hinkel 2013; Ramirez & Kaltenegger 2016). The width
and distance of the HZ region depends in a first approximation on the incident
stellar flux and the atmospheric composition of the planet.

The limits of the conservative HZ is set by the greenhouse effect of CO2 and
H2O vapour. At the outer edge, condensation and scattering by CO2 outstrips
its greenhouse capacity. At the inner edge, the mean surface temperatures
exceed the critical point of water, triggering a runaway greenhouse
environment and rapid water loss to space on short timescales (see, e.g.,
Kasting et al. 1993; Kopparapu et al. 2014; Kaltenegger 2017). Modelling
clouds is complex especially because no data exists that can provide
comparison data sets for, e.g., fast rotating or slow rotating, very hot or
very cold Earth-like planets. However, both 1D and 3D models generally find
wider boundaries when considering cloud feedback than the original classical
HZ but models differ in their specific results (see, e.g., the review in
Kaltenegger 2017). The updated conservative HZ limits (Kopparapu et al., 2014)
have been scaled to the 3D models by Leconte et al. (2013) to account for some
cloud feedback and are used here.

The concept of the empirical (or optimistic) HZ relies on observations in our
own Solar System (Kasting et al., 1993). It is based on the incident flux on
Venus and Mars when there were no more hints of liquid water on their surface.
In our Solar System, that translates to 1.76 current Earth’s incident
irradiation ($S_{\oplus}$) for Venus around 1 billion years ago and 0.32
$S_{\oplus}$ for Mars about 3.8 billion years ago, providing empirical proxies
for the incident flux limits of habitability. The inner edge of the empirical
HZ is not well known because of the geological active history of Venus’
surface. The outer edge of the HZ in 3D models agrees with incident stellar
flux for Early Mars insulation. The empirical HZ provides an empirical dataset
for the search for habitable planets, complementary to modelling efforts.

Cool stars warm an Earth-like planet’s surface more effectively than hotter
stars (see, e.g., Kasting et al. 1993). Thus, for LP 890-9, the inner limits
of the HZ are 1.49 $S_{\oplus}$ for the empirical and 0.92 $S_{\oplus}$ for
the conservative HZ respectively (Kopparapu et al., 2013). The outer limit of
the HZ for LP 890-9 is found at a stellar flux of about 0.21 $S_{\oplus}$.
While the HZ limits can change with atmospheric composition (e.g.,
Pierrehumbert & Gaidos 2011; Ramirez & Kaltenegger 2017, 2018), here we used
the empirical HZ and the conservative HZ defined for a planet with an
N2-H2O-CO2 dominated atmosphere like Earth to identify the HZ limits for the
LP 890-9 system.

With an incident stellar flux of $4.09\pm 0.12$ $S_{\oplus}$, LP 890-9 b is
too close to the star to be within the HZ. However, LP 890-9 c receives
$0.906\pm 0.026$ $S_{\oplus}$, which places it within the limits of both the
empirical and the conservative HZ. It orbits close to the inner edge of the
conservative HZ, and in the inner third of the empirical HZ of LP 890-9 (see
Fig. 17). Due to increasing stellar irradiation towards the inner edge of the
HZ, LP 890-9 c should develop a water-rich atmosphere because more surface
water would evaporate than on Earth due to increased temperatures. Such an
increase in water vapour in the atmosphere should be detectable with JWST and
upcoming Extremely Large Telescopes (see Sect. 8.4).

While LP 890-9 c orbits within the HZ of its host star, due to its position
close to a cool M star, a potential slow rotation due to tidal forces could
maintain habitability even for higher incident fluxes than the nominal limits
of the HZ, due to an increased cloud coverage on the day-side of the planet.
This makes planets that could be synchronously locked at the inner edge of the
HZ, like LP 890-9 c, very interesting test-cases for this hypothesis (see,
e.g., Yang et al. 2013; Wolf & Toon 2014; Kopparapu et al. 2016).

###  8.4 Prospects for atmospheric characterisation

Rocky exoplanets orbiting in the HZ of nearby late-type M dwarf stars provide
unique opportunities for characterising their atmospheres as well as searching
for biosignature gases. To assess the potential of the system for atmospheric
characterisation in the context of other known transiting terrestrial planets,
we followed the same approach as for K2-315 b in Niraula et al. (2020). We
retrieved the currently known transiting exoplanets with a reported radius
below 1.6 $R_{\oplus}$ from NASA Exoplanet
Archive131313https://exoplanetarchive.ipac.caltech.edu/. We focussed here on
temperate planets with mild incident irradiations, lower than five times that
received by the Earth ($S_{\rm{p}}<5\>S_{\oplus}$), and estimated their
typical signal amplitude in transit transmission spectroscopy as:

$\displaystyle\Delta\mathrm{d}F$
$\displaystyle=\frac{2R_{\rm{p}}h_{\rm{eff}}}{R_{\star}^{2}}$ (3)

where $h_{\rm{eff}}$ is the effective atmospheric height (that is, the extent
of the atmospheric annulus probed in transmission). $h_{\rm{eff}}$ is directly
proportional to the planet’s atmospheric scale height $H$:
$h_{\rm{eff}}=N_{\rm{H}}\>H=N_{\rm{H}}\>\frac{kT}{\mu g}$, where $k$ is
Boltzmann’s constant, $T$ is the atmospheric temperature, $\mu$ is the
atmospheric mean molecular mass, and $g$ is the surface gravity. $N_{\rm{H}}$
is a number ranging typically between 6 and 10, depending on the presence of
clouds, the spectral resolution and range covered (Miller-Ricci et al., 2009;
de Wit & Seager, 2013). We assumed here $N_{\rm{H}}=7$, $\mu=20$ amu (relevant
for scenarios of secondary atmospheres), and $T$ to be equal to the
equilibrium temperature, calculated for a sphere of null Bond albedo and an
efficient heat redistribution. For planets with unknown masses, we estimated
$g$ using the probabilistic mass–radius relationship of Chen & Kipping (2017).

Fig. 18 shows the estimated transmission signals as a function of the planets’
incident stellar flux. The points are colour-coded as a function of their S/N
relative to TRAPPIST-1 b, obtained by scaling the signal amplitude with the
hosts’ brightness in the $J$-band and using TRAPPIST-1 b’s S/N as a reference.
LP 890-9 b compares well to the outer TRAPPIST-1 planets in terms of potential
for atmospheric characterisation, while LP 890-9 c stands out as the second-
most favourable habitable-zone terrestrial planet after the TRAPPIST-1
planets. The projected signals in transmission are $\sim$125 and $\sim$88 ppm
for LP 890-9 b and c, respectively, with corresponding signal-to-noise ratios
relative to TRAPPIST-1 b of 0.38 and 0.26. This means that assuming comparable
secondary atmospheres, LP 890-9 b and c would require respectively 7 and 15
times the amount of in-transit time to be characterised to the same extent as
TRAPPIST-1 b. Assessing the presence of a 10 bar CO2 atmosphere with
JWST/NIRSpec would thus require around 14 and 30 transits for planets b and c,
respectively (scaled from Lustig-Yaeger et al. 2019). Repeating the same
exercise to compare LP 890-9 c with TRAPPIST-1 d and e, which receive closer
amounts of stellar flux, we find a S/N of 0.29 relative to TRAPPIST-1 d and
0.38 relative to TRAPPIST-1 e. Characterising LP 890-9 c to the same extent as
TRAPPIST-1 d (resp. TRAPPIST-1 e) would thus require 12 (resp. 7) times the
amount of in-transit time. Based on these estimations, about 90 transits would
be needed to detect a 1 bar H2O atmosphere around LP 890-9 c with JWST/NIRSpec
(scaled again from Lustig-Yaeger et al. 2019). This rough comparison assumes a
similar atmosphere for all planets. However, since LP 890-9 c is located close
to the inner limit (runaway greenhouse) of the conservative HZ (see Sect.
8.3), it may have a thicker H2O atmosphere which would boost its atmospheric
signals.

##  9 Conclusions

We have presented the discovery and initial characterisation of the LP 890-9
system, which hosts two temperate super-Earths transiting a nearby M6 dwarf.
LP 890-9 is the second-coolest star found to host planets after TRAPPIST-1.
The seven terrestrial planets orbiting TRAPPIST-1 have garnered the interest
of a broad community across very diverse scientific disciplines: planet
formation and evolution (e.g., Ormel et al. 2017; Coleman et al. 2019; Raymond
et al. 2022), star-planet interactions (e.g., Bolmont et al. 2017; O’Malley-
James & Kaltenegger 2017; Dong et al. 2018), multi-planet dynamics (e.g.,
Grimm et al. 2018; Agol et al. 2021; Teyssandier et al. 2022), interior
modelling (e.g., Dorn et al. 2018; Barr et al. 2018; Barth et al. 2021),
atmospheric observations (e.g., de Wit et al. 2016; Bourrier et al. 2017; de
Wit et al. 2018) and modelling (e.g., Lustig-Yaeger et al. 2019; Pidhorodetska
et al. 2020; Lin et al. 2021), or climate predictions (e.g., Turbet et al.
2018, 2021; Sergeev et al. 2021), among others. The discovery of the
remarkable LP 890-9 system presented in this work offers another rare
opportunity to study temperate terrestrial planets around our smallest and
coolest neighbours.

###### Acknowledgements.

Funding for the TESS mission is provided by NASA’s Science Mission
Directorate. We acknowledge the use of public TESS data from pipelines at the
TESS Science Office and at the TESS Science Processing Operations Center. This
research has made use of the Exoplanet Follow-up Observation Program website,
which is operated by the California Institute of Technology, under contract
with the National Aeronautics and Space Administration under the Exoplanet
Exploration Program. Resources supporting this work were provided by the NASA
High-End Computing (HEC) Program through the NASA Advanced Supercomputing
(NAS) Division at Ames Research Center for the production of the SPOC data
products. This paper includes data collected by the TESS mission that are
publicly available from the Mikulski Archive for Space Telescopes (MAST). The
research leading to these results has received funding from the European
Research Council (ERC) under the FP/2007–2013 ERC grant agreement n∘ 336480,
and under the European Union’s Horizon 2020 research and innovation programme
(grants agreements n∘ 679030 & 803193/BEBOP); from an Action de Recherche
Concertée (ARC) grant, financed by the Wallonia–Brussels Federation, from the
Balzan Prize Foundation, from the BELSPO/BRAIN2.0 research program (PORTAL
project), from the Science and Technology Facilities Council (STFC; grants n∘
ST/S00193X/1, ST/00305/1, and ST/W000385/1), and from F.R.S-FNRS (Research
Project ID T010920F). This work was also partially supported by a grant from
the Simons Foundation (PI Queloz, grant number 327127), as well as by the
MERAC foundation (PI Triaud). TRAPPIST is funded by the Belgian Fund for
Scientific Research (Fond National de la Recherche Scientifique, FNRS) under
the grant PDR T.0120.21, with the participation of the Swiss National Science
Fundation (SNF). This work is partly supported by MEXT/JSPS KAKENHI Grant
Numbers JP15H02063, JP17H04574, JP18H05439, JP18H05442, JP19K14783,
JP21H00035, JP21K13975, JP21K20376, JP22000005, Grant-in-Aid for JSPS Fellows
Grant Number JP20J21872, JST CREST Grant Number JPMJCR1761, the Astrobiology
Center of National Institutes of Natural Sciences (NINS) (Grant Numbers
AB031010, AB031014), and Social welfare juridical person SHIYUKAI (chairman
MASAYUKI KAWASHIMA). This paper is based on data collected at the Subaru
Telescope, which is located atop Maunakea and operated by the National
Astronomical Observatory of Japan (NAOJ). We wish to recognise and acknowledge
the very significant cultural role and reverence that the summit of Maunakea
has always had within the indigenous Hawaiian community. This paper is based
on observations made with the MuSCAT3 instrument, developed by the
Astrobiology Center and under financial supports by JSPS KAKENHI (JP18H05439)
and JST PRESTO (JPMJPR1775), at Faulkes Telescope North on Maui, HI, operated
by the Las Cumbres Observatory. Some of the observations in the paper made use
of the High-Resolution Imaging instrument Zorro obtained under Gemini LLP
Proposal Number: GN/S-2021A-LP-105. Zorro was funded by the NASA Exoplanet
Exploration Program and built at the NASA Ames Research Center by Steve B.
Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. Zorro was mounted on
the Gemini North (and/or South) telescope of the international Gemini
Observatory, a program of NSF’s OIR Lab, which is managed by the Association
of Universities for Research in Astronomy (AURA) under a cooperative agreement
with the National Science Foundation on behalf of the Gemini partnership: the
National Science Foundation (United States), National Research Council
(Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio
de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência,
Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space
Science Institute (Republic of Korea). We acknowledge funding from the
European Research Council under the ERC Grant Agreement n. 337591-ExTrA. This
work has been carried out within the framework of the NCCR PlanetS supported
by the Swiss National Science Foundation. This work is based upon observations
carried out at the Observatorio Astronómico Nacional at the Sierra de San
Pedro Mártir (OAN-SPM), Baja California, México. We warmly thank the entire
technical staff of the Observatorio Astronómico Nacional at San Pedro Mártir
for their unfailing support to SAINT-EX operations. Research at Lick
Observatory is partially supported by a generous gift from Google. L.D. is an
F.R.S.-FNRS Postdoctoral Researcher. M.G. and E.J. are F.R.S.-FNRS Senior
Research Associates. V.V.G. is an F.R.S.-FNRS Research Associate. B.V.R.
thanks the Heising-Simons Foundation for support. Y.G.M.C. acknowledges
support from UNAM-PAPIIT IG-101321. B.-O.D. acknowledges support from the
Swiss National Science Foundation (PP00P2-163967 and PP00P2-190080). M.N.G.
acknowledges support from the European Space Agency (ESA) as an ESA Research
Fellow. A.H.M.J.T acknowledges funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement n∘ 803193/BEBOP), from the MERAC foundation, and
from the Science and Technology Facilities Council (STFC; grants n∘
ST/S00193X/1, ST/00305/1, and ST/W000385/1). E.D. acknowledges support from
the innovation and research Horizon 2020 program in the context of the Marie
Sklodowska-Curie subvention 945298. V.K. acknowledges support from NSF award
AST2009343. This publication benefits from the support of the French Community
of Belgium in the context of the FRIA Doctoral Grant awarded to M.T.

## References

  * Agol et al. (2021) Agol, E., Dorn, C., Grimm, S. L., et al. 2021, The Planetary Science Journal, 2, 1
  * Aller et al. (2020) Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128
  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bailey (2012) Bailey, S. 2012, PASP, 124, 1015
  * Bakos et al. (2006) Bakos, G. Á., Pál, A., Latham, D. W., Noyes, R. W., & Stefanik, R. P. 2006, ApJ, 641, L57
  * Barclay et al. (2018) Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2
  * Barnes (2017) Barnes, R. 2017, Celestial Mechanics and Dynamical Astronomy, 129, 509
  * Barr et al. (2018) Barr, A. C., Dobos, V., & Kiss, L. L. 2018, A&A, 613, A37
  * Barth et al. (2021) Barth, P., Carone, L., Barnes, R., et al. 2021, Astrobiology, 21, 1325
  * Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527
  * Blanco-Cuaresma & Bolmont (2017) Blanco-Cuaresma, S. & Bolmont, E. 2017, in IAU Symposium, Vol. 325, Astroinformatics, ed. M. Brescia, S. G. Djorgovski, E. D. Feigelson, G. Longo, & S. Cavuoti, 341–344
  * Bochanski et al. (2007) Bochanski, J. J., West, A. A., Hawley, S. L., & Covey, K. R. 2007, AJ, 133, 531
  * Bolmont et al. (2020) Bolmont, E., Demory, B. O., Blanco-Cuaresma, S., et al. 2020, A&A, 635, A117
  * Bolmont et al. (2014) Bolmont, E., Raymond, S. N., von Paris, P., et al. 2014, The Astrophysical Journal, 793, 3
  * Bolmont et al. (2017) Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728
  * Bonfils et al. (2015) Bonfils, X., Almenara, J. M., Jocou, L., et al. 2015, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9605, Techniques and Instrumentation for Detection of Exoplanets VII, ed. S. Shaklan, 96051L
  * Bourrier et al. (2017) Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121
  * Buder et al. (2021) Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150
  * Burdanov et al. (2018) Burdanov, A., Delrez, L., Gillon, M., & Jehin, E. 2018, SPECULOOS Exoplanet Search and Its Prototype on TRAPPIST, ed. H. J. Deeg & J. A. Belmonte, 130
  * Burgasser & Mamajek (2017) Burgasser, A. J. & Mamajek, E. E. 2017, ApJ, 845, 110
  * Burgasser & Splat Development Team (2017) Burgasser, A. J. & Splat Development Team. 2017, in Astronomical Society of India Conference Series, Vol. 14, Astronomical Society of India Conference Series, 7–12
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560
  * Chen & Kipping (2017) Chen, J. & Kipping, D. 2017, ApJ, 834, 17
  * Ciardi et al. (2015) Ciardi, D. R., Beichman, C. A., Horch, E. P., & Howell, S. B. 2015, ApJ, 805, 16
  * Cincotta & Simó (1999) Cincotta, P. & Simó, C. 1999, Celestial Mechanics and Dynamical Astronomy, 73, 195
  * Cincotta et al. (2003) Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151
  * Cincotta & Simó (2000) Cincotta, P. M. & Simó, C. 2000, Astronomy and Astrophysics Supplement Series, 147, 205
  * Cloutier et al. (2020) Cloutier, R., Eastman, J. D., Rodriguez, J. E., et al. 2020, AJ, 160, 3
  * Cointepas et al. (2021) Cointepas, M., Almenara, J. M., Bonfils, X., et al. 2021, A&A, 650, A145
  * Coleman et al. (2019) Coleman, G. A. L., Leleu, A., Alibert, Y., & Benz, W. 2019, A&A, 631, A7
  * Crossfield et al. (2019) Crossfield, I. J. M., Waalkes, W., Newton, E. R., et al. 2019, ApJ, 883, L16
  * Cruz & Reid (2002) Cruz, K. L. & Reid, I. N. 2002, AJ, 123, 2828
  * Cushing et al. (2004) Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * de Wit & Seager (2013) de Wit, J. & Seager, S. 2013, Science, 342, 1473
  * de Wit et al. (2016) de Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Nature, 537, 69
  * de Wit et al. (2018) de Wit, J., Wakeford, H. R., Lewis, N. K., et al. 2018, Nature Astronomy, 2, 214
  * Delrez et al. (2021) Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nature Astronomy, 5, 775
  * Delrez et al. (2018) Delrez, L., Gillon, M., Queloz, D., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10700, Ground-based and Airborne Telescopes VII, ed. H. K. Marshall & J. Spyromilio, 107001I
  * Demory et al. (2020) Demory, B. O., Pozuelos, F. J., Gómez Maqueo Chew, Y., et al. 2020, A&A, 642, A49
  * Dhital et al. (2012) Dhital, S., West, A. A., Stassun, K. G., et al. 2012, AJ, 143, 67
  * Dittmann et al. (2017) Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333
  * Donati et al. (2020) Donati, J. F., Kouach, D., Moutou, C., et al. 2020, MNRAS, 498, 5684
  * Dong et al. (2018) Dong, C., Jin, M., Lingam, M., et al. 2018, Proceedings of the National Academy of Science, 115, 260
  * Dorn et al. (2018) Dorn, C., Mosegaard, K., Grimm, S. L., & Alibert, Y. 2018, ApJ, 865, 20
  * Dorn et al. (2014) Dorn, R. J., Anglada-Escude, G., Baade, D., et al. 2014, The Messenger, 156, 7
  * Douglas et al. (2014) Douglas, S. T., Agüeros, M. A., Covey, K. R., et al. 2014, ApJ, 795, 161
  * Dreizler et al. (2020) Dreizler, S., Jeffers, S. V., Rodríguez, E., et al. 2020, MNRAS, 493, 536
  * Dévora-Pajares & Pozuelos (2022) Dévora-Pajares, M. & Pozuelos, F. J. 2022, MATRIX: Multi-phAse Transits Recovery from Injected eXoplanets
  * Eggleton et al. (1998) Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, The Astrophysical Journal, 499, 853
  * Eisner et al. (2020) Eisner, N. L., Barragán, O., Aigrain, S., et al. 2020, MNRAS, 494, 750
  * Espinoza et al. (2019) Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262
  * Fernandes et al. (2019) Fernandes, C. S., Van Grootel, V., Salmon, S. J. A. J., et al. 2019, ApJ, 879, 94
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64
  * Fukui et al. (2011) Fukui, A., Narita, N., Tristram, P. J., et al. 2011, PASJ, 63, 287
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 650, C3
  * Gan et al. (2022) Gan, T., Soubkiou, A., Wang, S. X., et al. 2022, MNRAS, 514, 4120
  * Garcia et al. (2022) Garcia, L. J., Timmermans, M., Pozuelos, F. J., et al. 2022, MNRAS, 509, 4817
  * Gelman & Rubin (1992) Gelman, A. & Rubin, D. B. 1992, Statistical Science, 7, 457
  * Geman & Geman (1984) Geman, S. & Geman, D. 1984, IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6, 721
  * Giacalone & Dressing (2020) Giacalone, S. & Dressing, C. D. 2020, triceratops: Candidate exoplanet rating tool
  * Giacalone et al. (2021) Giacalone, S., Dressing, C. D., Jensen, E. L. N., et al. 2021, AJ, 161, 24
  * Gillon (2018) Gillon, M. 2018, Nature Astronomy, 2, 344
  * Gillon et al. (2014) Gillon, M., Demory, B. O., Madhusudhan, N., et al. 2014, A&A, 563, A21
  * Gillon et al. (2013) Gillon, M., Jehin, E., Fumel, A., Magain, P., & Queloz, D. 2013, in European Physical Journal Web of Conferences, Vol. 47, European Physical Journal Web of Conferences, 03001
  * Gillon et al. (2016) Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221
  * Gillon et al. (2011) Gillon, M., Jehin, E., Magain, P., et al. 2011, in European Physical Journal Web of Conferences, Vol. 11, European Physical Journal Web of Conferences, 06002
  * Gillon et al. (2020) Gillon, M., Meadows, V., Agol, E., et al. 2020, arXiv e-prints, arXiv:2002.04798
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456
  * Gillon et al. (2012) Gillon, M., Triaud, A. H. M. J., Fortney, J. J., et al. 2012, A&A, 542, A4
  * Grimm et al. (2018) Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39
  * Günther et al. (2019) Günther, M. N., Pozuelos, F. J., Dittmann, J. A., et al. 2019, Nature Astronomy, 3, 1099
  * Hadden (2019) Hadden, S. 2019, shadden/TTV2Fast2Furious: First release of TTV2Fast2Furious
  * Haghighipour & Kaltenegger (2013) Haghighipour, N. & Kaltenegger, L. 2013, The Astrophysical Journal, 777, 166
  * Hamuy et al. (1994) Hamuy, M., Suntzeff, N. B., Heathcote, S. R., et al. 1994, PASP, 106, 566
  * Hamuy et al. (1992) Hamuy, M., Walker, A. R., Suntzeff, N. B., et al. 1992, PASP, 104, 533
  * Hastings (1970) Hastings, W. K. 1970, Biometrika, 57, 97
  * Hinse et al. (2015) Hinse, T. C., Haghighipour, N., Kostov, V. B., & Goździewski, K. 2015, ApJ, 799, 88
  * Hippke et al. (2019) Hippke, M., David, T. J., Mulders, G. D., & Heller, R. 2019, AJ, 158, 143
  * Hippke & Heller (2019) Hippke, M. & Heller, R. 2019, A&A, 623, A39
  * Hirano et al. (2020) Hirano, T., Kuzuhara, M., Kotani, T., et al. 2020, PASJ, 72, 93
  * Howell et al. (2016) Howell, S. B., Everett, M. E., Horch, E. P., et al. 2016, ApJ, 829, L2
  * Howell et al. (2011) Howell, S. B., Everett, M. E., Sherry, W., Horch, E., & Ciardi, D. R. 2011, AJ, 142, 19
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Hut (1981) Hut, P. 1981, Astronomy and Astrophysics, 99, 126
  * Irwin et al. (2004) Irwin, M. J., Lewis, J., Hodgkin, S., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5493, Optimizing Scientific Return for Astronomy through Information Technologies, ed. P. J. Quinn & A. Bridger, 411–422
  * Jehin et al. (2018) Jehin, E., Gillon, M., Queloz, D., et al. 2018, The Messenger, 174, 2
  * Jehin et al. (2011) Jehin, E., Gillon, M., Queloz, D., et al. 2011, The Messenger, 145, 2
  * Jenkins (2002) Jenkins, J. M. 2002, ApJ, 575, 493
  * Jenkins et al. (2010) Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al. 2010, ApJ, 713, L87
  * Jenkins et al. (2020) Jenkins, J. M., Tenenbaum, P., Seader, S., et al. 2020, Kepler Data Processing Handbook: Transiting Planet Search, Kepler Science Document KSCI-19081-003
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, ed. G. Chiozzi & J. C. Guzman, 99133E
  * Jenkins et al. (2019) Jenkins, J. S., Pozuelos, F. J., Tuomi, M., et al. 2019, MNRAS, 490, 5585
  * Kaltenegger (2017) Kaltenegger, L. 2017, Annual Review of Astronomy and Astrophysics, 55, 433
  * Kaltenegger & Haghighipour (2013) Kaltenegger, L. & Haghighipour, N. 2013, The Astrophysical Journal, 777, 165
  * Kaltenegger & Traub (2009) Kaltenegger, L. & Traub, W. A. 2009, ApJ, 698, 519
  * Kane & Hinkel (2013) Kane, S. R. & Hinkel, N. R. 2013, The Astrophysical Journal, 762, 7
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152
  * Kirkpatrick et al. (2010) Kirkpatrick, J. D., Looper, D. L., Burgasser, A. J., et al. 2010, ApJS, 190, 100
  * Kopparapu et al. (2016) Kopparapu, Ravi Kumar, R., Wolf, E. T., Haqq-Misra, J., et al. 2016, The Astrophysical Journal, 819, 84
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131
  * Kopparapu et al. (2014) Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, The Astrophysical Journal, 787, L29
  * Kostov et al. (2019) Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32
  * Kotani et al. (2018) Kotani, T., Tamura, M., Nishikawa, J., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 1070211
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
  * Leconte et al. (2010) Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 2010, Astronomy and Astrophysics, 516, A64
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Nature, 504, 268
  * Lépine et al. (2013) Lépine, S., Hilton, E. J., Mann, A. W., et al. 2013, AJ, 145, 102
  * Lépine et al. (2003) Lépine, S., Rich, R. M., & Shara, M. M. 2003, AJ, 125, 1598
  * Lépine et al. (2007) Lépine, S., Rich, R. M., & Shara, M. M. 2007, ApJ, 669, 1235
  * Lester et al. (2021) Lester, K. V., Matson, R. A., Howell, S. B., et al. 2021, AJ, 162, 75
  * Li et al. (2019) Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506
  * Lienhard et al. (2020) Lienhard, F., Queloz, D., Gillon, M., et al. 2020, MNRAS, 497, 3790
  * Lin et al. (2021) Lin, Z., MacDonald, R. J., Kaltenegger, L., & Wilson, D. J. 2021, MNRAS, 505, 3562
  * Lissauer et al. (2012) Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, ApJ, 750, 112
  * Lithwick et al. (2012) Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447
  * Luger et al. (2017) Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nature Astronomy, 1, 0129
  * Lustig-Yaeger et al. (2019) Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, AJ, 158, 27
  * Luyten (1979) Luyten, W. J. 1979, New Luyten catalogue of stars with proper motions larger than two tenths of an arcsecond; and first supplement; NLTT. (Minneapolis (1979)); Label 12 = short description; Label 13 = documentation by Warren; Label 14 = catalogue
  * Mandel & Agol (2002) Mandel, K. & Agol, E. 2002, ApJ, 580, L171
  * Mann et al. (2013) Mann, A. W., Brewer, J. M., Gaidos, E., Lépine, S., & Hilton, E. J. 2013, AJ, 145, 52
  * Mann et al. (2014) Mann, A. W., Deacon, N. R., Gaidos, E., et al. 2014, AJ, 147, 160
  * Marconi et al. (2021) Marconi, A., Abreu, M., Adibekyan, V., et al. 2021, The Messenger, 182, 27
  * McCormac et al. (2013) McCormac, J., Pollacco, D., Skillen, I., et al. 2013, PASP, 125, 548
  * McCully et al. (2018) McCully, C., Volgenau, N. H., Harbeck, D.-R., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10707, Software and Cyberinfrastructure for Astronomy V, ed. J. C. Guzman & J. Ibsen, 107070K
  * Ment et al. (2019) Ment, K., Dittmann, J. A., Astudillo-Defru, N., et al. 2019, AJ, 157, 32
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, J. Chem. Phys., 21, 1087
  * Mignard (1979) Mignard, F. 1979, The Moon and the Planets, 20, 301
  * Miguel et al. (2020) Miguel, Y., Cridland, A., Ormel, C. W., Fortney, J. J., & Ida, S. 2020, MNRAS, 491, 1998
  * Miller-Ricci et al. (2009) Miller-Ricci, E., Seager, S., & Sasselov, D. 2009, ApJ, 690, 1056
  * Minkowski & Abell (1963) Minkowski, R. L. & Abell, G. O. 1963, The National Geographic Society-Palomar Observatory Sky Survey, ed. K. A. Strand, 481
  * Morello et al. (2022) Morello, G., Parviainen, H., Murgas, F., et al. 2022, arXiv e-prints, arXiv:2201.13274
  * Mori et al. (2022) Mori, M., Livingston, J. H., Leon, J. d., et al. 2022, AJ, 163, 298
  * Muirhead et al. (2018) Muirhead, P. S., Dressing, C. D., Mann, A. W., et al. 2018, AJ, 155, 180
  * Mulders et al. (2021) Mulders, G. D., Drażkowska, J., van der Marel, N., Ciesla, F. J., & Pascucci, I. 2021, ApJ, 920, L1
  * Murray et al. (2020) Murray, C. A., Delrez, L., Pedersen, P. P., et al. 2020, MNRAS, 495, 2446
  * Murray et al. (2022) Murray, C. A., Queloz, D., Gillon, M., et al. 2022, MNRAS, 513, 2615
  * Narita et al. (2020) Narita, N., Fukui, A., Yamamuro, T., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114475K
  * Neron de Surgy & Laskar (1997) Neron de Surgy, O. & Laskar, J. 1997, A&A, 318, 975
  * Newton et al. (2017) Newton, E. R., Irwin, J., Charbonneau, D., et al. 2017, ApJ, 834, 85
  * Niraula et al. (2020) Niraula, P., Wit, J. d., Rackham, B. V., et al. 2020, AJ, 160, 172
  * O’Malley-James & Kaltenegger (2017) O’Malley-James, J. T. & Kaltenegger, L. 2017, MNRAS, 469, L26
  * Ormel et al. (2017) Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1
  * Parviainen (2015) Parviainen, H. 2015, MNRAS, 450, 3233
  * Parviainen & Aigrain (2015) Parviainen, H. & Aigrain, S. 2015, MNRAS, 453, 3821
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96
  * Pidhorodetska et al. (2020) Pidhorodetska, D., Fauchez, T. J., Villanueva, G. L., Domagal-Goldman, S. D., & Kopparapu, R. K. 2020, ApJ, 898, L33
  * Pierrehumbert & Gaidos (2011) Pierrehumbert, R. & Gaidos, E. 2011, The Astrophysical Journal, 734, L13
  * Pozuelos et al. (2020) Pozuelos, F. J., Suárez, J. C., de Elía, G. C., et al. 2020, A&A, 641, A23
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing
  * Ramirez & Kaltenegger (2016) Ramirez, R. M. & Kaltenegger, L. 2016, The Astrophysical Journal, 823, 6
  * Ramirez & Kaltenegger (2017) Ramirez, R. M. & Kaltenegger, L. 2017, The Astrophysical Journal, 837, L4
  * Ramirez & Kaltenegger (2018) Ramirez, R. M. & Kaltenegger, L. 2018, The Astrophysical Journal, 858, 72
  * Rasmussen & Williams (2006) Rasmussen, C. E. & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning
  * Raymond et al. (2022) Raymond, S. N., Izidoro, A., Bolmont, E., et al. 2022, Nature Astronomy, 6, 80
  * Rayner et al. (2003) Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362
  * Reid et al. (1991) Reid, I. N., Brewer, C., Brucato, R. J., et al. 1991, PASP, 103, 661
  * Rein & Liu (2012) Rein, H. & Liu, S.-F. 2012, 537, A128
  * Rein & Tamayo (2015) Rein, H. & Tamayo, D. 2015, 452, 376
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Riddick et al. (2007) Riddick, F. C., Roche, P. F., & Lucas, P. W. 2007, MNRAS, 381, 1067
  * Rojas-Ayala et al. (2012) Rojas-Ayala, B., Covey, K. R., Muirhead, P. S., & Lloyd, J. P. 2012, ApJ, 748, 93
  * Sabin et al. (2018) Sabin, L., Gómez Maqueo Chew, Y., Demory, B.-O., Petrucci, R., & Saint-Ex Consortium. 2018, in 20th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 59
  * Savitzky & Golay (1964) Savitzky, A. & Golay, M. J. E. 1964, Analytical Chemistry, 36, 1627
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Schanche et al. (2022) Schanche, N., Pozuelos, F. J., Günther, M. N., et al. 2022, A&A, 657, A45
  * Schwarz (1978) Schwarz, G. 1978, Annals of Statistics, 6, 461
  * Scott et al. (2021) Scott, N. J., Howell, S. B., Gnilka, C. L., et al. 2021, Frontiers in Astronomy and Space Sciences, 8, 138
  * Seager & Mallén-Ornelas (2003) Seager, S. & Mallén-Ornelas, G. 2003, ApJ, 585, 1038
  * Sebastian et al. (2021) Sebastian, D., Gillon, M., Ducrot, E., et al. 2021, A&A, 645, A100
  * Seifahrt et al. (2020) Seifahrt, A., Bean, J. L., Stürmer, J., et al. 2020, in Ground-based and Airborne Instrumentation for Astronomy VIII, ed. C. J. Evans, J. J. Bryant, & K. Motohara, Vol. 11447, International Society for Optics and Photonics (SPIE), 305 – 325
  * Sergeev et al. (2021) Sergeev, D. E., Fauchez, T. J., Turbet, M., et al. 2021, arXiv e-prints, arXiv:2109.11459
  * Sharma et al. (2018) Sharma, S., Stello, D., Buder, S., et al. 2018, MNRAS, 473, 2004
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136
  * Stassun et al. (2018a) Stassun, K. G., Corsaro, E., Pepper, J. A., & Gaudi, B. S. 2018a, AJ, 155, 22
  * Stassun et al. (2018b) Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018b, AJ, 156, 102
  * Stassun & Torres (2016) Stassun, K. G. & Torres, G. 2016, AJ, 152, 180
  * Stassun & Torres (2021) Stassun, K. G. & Torres, G. 2021, ApJ, 907, L33
  * Stumpe et al. (2014) Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985
  * Sullivan et al. (2015) Sullivan, P. W., Winn, J. N., Berta-Thompson, Z. K., et al. 2015, ApJ, 809, 77
  * Tamura et al. (2012) Tamura, M., Suto, H., Nishikawa, J., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84461T
  * Teyssandier et al. (2022) Teyssandier, J., Libert, A. S., & Agol, E. 2022, A&A, 658, A170
  * Turbet et al. (2018) Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86
  * Turbet et al. (2021) Turbet, M., Fauchez, T. J., Sergeev, D. E., et al. 2021, arXiv e-prints, arXiv:2109.11457
  * Twicken et al. (2018) Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502
  * Van Grootel et al. (2018) Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30
  * Vanderspek et al. (2019) Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, ApJ, 871, L24
  * Wells et al. (2021) Wells, R. D., Rackham, B. V., Schanche, N., et al. 2021, A&A, 653, A97
  * Wildi et al. (2017) Wildi, F., Blind, N., Reshetov, V., et al. 2017, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10400, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1040018
  * Wolf & Toon (2014) Wolf, E. T. & Toon, O. B. 2014, Geophysical Research Letters, 41, 167
  * Yang et al. (2013) Yang, J., Cowan, N. B., & Abbot, D. S. 2013, The Astrophysical Journal, 771, L45
  * Zechmeister et al. (2019) Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49
  * Zhang et al. (2019) Zhang, S., Luo, A. L., Comte, G., et al. 2019, ApJS, 240, 31

##  Appendix A Radial velocities

$\mathrm{BJD_{TDB}}-2\,450\,000$ | RV (km/s) |  $\sigma_{\mathrm{RV}}$ (km/s)
---|---|---
$YJ$-band |  |
9467.104595 | 8.975819 | 0.01608
9467.125868 | 8.926950 | 0.01971
9488.104150 | 8.864151 | 0.00970
9508.011182 | 8.920489 | 0.01211
9509.985105 | 8.945039 | 0.01015
9510.006380 | 8.924036 | 0.01078
9515.018797 | 8.906155 | 0.01468
9515.040059 | 8.958196 | 0.01124
9515.061315 | 8.898919 | 0.01441
9530.936633 | 8.911857 | 0.01190
9530.957919 | 8.834334 | 0.01502
$H$-band |  |
9467.104595 | 3.504429 | 0.01054
9467.125868 | 3.520200 | 0.01299
9488.104150 | 3.452051 | 0.00735
9508.011182 | 3.450579 | 0.00872
9509.985105 | 3.524839 | 0.00767
9510.006380 | 3.500336 | 0.00786
9515.018797 | 3.458665 | 0.01036
9515.040059 | 3.508536 | 0.00820
9515.061315 | 3.492029 | 0.00974
9530.936633 | 3.489667 | 0.00914
9530.957919 | 3.464614 | 0.01166
9588.854504 | 3.497201 | 0.00797
9588.875855 | 3.508872 | 0.00809
9588.897185 | 3.495941 | 0.00792
Table 6: Subaru/IRD radial velocity measurements. We recommend treating the
$YJ$-band and $H$-band RVs as independent datasets, as done in the analysis
presented in Sect. 6.

##  Appendix B Photometric baseline models and error scaling factors

Planet | Facility | Bandpass |  $T_{0}$ ($\mathrm{BJD_{TDB}}$ | Baseline model | Residual RMS | $\beta_{w}$ | $\beta_{r}$
---|---|---|---|---|---|---|---
|  |  |  $-$ 2 450 000) |  | (exp. time) |  |
b | TESS Sector 4 | TESS | 8413.1930 | offset | 0.0197 (120 s) | 0.85 | 1.25
| TESS Sector 4 | TESS | 8415.9229 | offset | 0.0218 (120 s) | 0.94 | 1.46
| TESS Sector 4 | TESS | 8426.8425 | offset | 0.0231 (120 s) | 1.00 | 1.00
| TESS Sector 4 | TESS | 8429.5724 | offset | 0.0217 (120 s) | 0.95 | 1.23
| TESS Sector 4 | TESS | 8432.3023 | offset | 0.0219 (120 s) | 0.95 | 1.00
| TESS Sector 4 | TESS | 8435.0322 | offset | 0.0200 (120 s) | 0.87 | 1.30
| TESS Sector 5 | TESS | 8440.4920 | offset | 0.0160 (120 s) | 0.90 | 1.67
| TESS Sector 5 | TESS | 8443.2219 | offset | 0.0145 (120 s) | 0.82 | 1.44
| TESS Sector 5 | TESS | 8445.9518 | offset | 0.0151 (120 s) | 0.86 | 1.06
| TESS Sector 5 | TESS | 8448.6817 | offset | 0.0158 (120 s) | 0.89 | 1.30
| TESS Sector 5 | TESS | 8454.1415 | offset | 0.0170 (120 s) | 0.96 | 1.12
| TESS Sector 5 | TESS | 8456.8714 | offset | 0.0156 (120 s) | 0.89 | 1.11
| TESS Sector 5 | TESS | 8459.6013 | offset | 0.0162 (120 s) | 0.92 | 1.29
| TESS Sector 5 | TESS | 8462.3312 | p($t^{1}$) | 0.0160 (120 s) | 0.91 | 1.10
| TESS Sector 31 | TESS | 9144.8071 | offset | 0.0194 (120 s) | 1.03 | 1.18
| TESS Sector 31 | TESS | 9147.5370 | offset | 0.0181 (120 s) | 0.97 | 1.29
| TESS Sector 31 | TESS | 9150.2669 | offset | 0.0178 (120 s) | 0.95 | 1.17
| TESS Sector 31 | TESS | 9152.9968 | offset | 0.0157 (120 s) | 0.84 | 1.26
| TESS Sector 31 | TESS | 9155.7267 | p($t^{1}$) | 0.0180 (120 s) | 0.93 | 1.62
| TESS Sector 31 | TESS | 9161.1865 | offset | 0.0158 (120 s) | 0.85 | 1.59
| TESS Sector 31 | TESS | 9163.9164 | offset | 0.0167 (120 s) | 0.90 | 1.78
| TESS Sector 31 | TESS | 9166.6463 | offset | 0.0189 (120 s) | 1.01 | 1.12
| TESS Sector 31 | TESS | 9169.3762 | offset | 0.0170 (120 s) | 0.91 | 1.16
| TESS Sector 32 | TESS | 9174.8360 | offset | 0.0217 (120 s) | 1.10 | 1.00
| TESS Sector 32 | TESS | 9177.5659 | offset | 0.0178 (120 s) | 0.90 | 1.00
| TESS Sector 32 | TESS | 9180.2958 | offset | 0.0183 (120 s) | 0.93 | 1.15
| TESS Sector 32 | TESS | 9183.0257 | offset | 0.0179 (120 s) | 0.91 | 1.12
| TESS Sector 32 | TESS | 9188.4855 | offset | 0.0184 (120 s) | 0.93 | 1.16
| TESS Sector 32 | TESS | 9191.2154 | offset | 0.0170 (120 s) | 0.86 | 1.56
| TESS Sector 32 | TESS | 9193.9453 | offset | 0.0180 (120 s) | 0.92 | 1.42
| TESS Sector 32 | TESS | 9196.6752 | offset | 0.0168 (120 s) | 0.86 | 1.13
| TESS Sector 32 | TESS | 9199.4051 | offset | 0.0173 (120 s) | 0.88 | 1.43
| SSO/Europa | $I+z$ | 9436.9067 | p($y^{1}$) | 0.00238 (35 s) | 0.77 | 1.10
| TRAPPIST-South | blue-blocking | 9436.9067 | p($a^{1}$) | 0.00329 (120 s) | 0.62 | 1.83
| SSO/Callisto | $i^{\prime}$ | 9447.8263 | p($a^{1}$) | 0.00236 (92 s) | 0.83 | 1.15
| SSO/Io | $I+z$ | 9447.8263 | p($a^{1}$) | 0.00236 (39 s) | 0.64 | 1.29
| SSO/Europa | $i^{\prime}$ | 9458.7460 | p($t^{1}$) | 0.00255 (92 s) | 0.58 | 1.78
| SSO/Io | $z^{\prime}$ | 9458.7460 | p($a^{1}$, $y^{1}$) | 0.00252 (49 s) | 0.76 | 1.10
| SSO/Ganymede | $r^{\prime}$ | 9477.8553 | p($a^{1}$) | 0.00616 (180 s) | 1.12 | 1.11
| SSO/Io | $r^{\prime}$ | 9477.8553 | p($t^{1}$) | 0.00584 (180 s) | 0.94 | 1.23
| SSO/Europa | $r^{\prime}$ | 9477.8553 | p($a^{1}$) | 0.00630 (180 s) | 1.10 | 1.00
| ExTrA1 | 0.85$-$1.55 $\mu$m | 9477.8553 | p($t^{4}$) | 0.00540 (60 s) | 1.02 | 1.43
| ExTrA2 | 0.85$-$1.55 $\mu$m | 9477.8553 | p($t^{4}$) | 0.00593 (60 s) | 0.95 | 1.00
| MuSCAT3 | $g^{\prime}$ | 9486.0450 | p($t^{1}$) | 0.0139 (120 s) | 0.89 | 1.22
| MuSCAT3 | $r^{\prime}$ | 9486.0450 | offset | 0.00430 (120 s) | 0.90 | 1.15
| MuSCAT3 | $i^{\prime}$ | 9486.0450 | p($y^{2}$) | 0.00214 (110 s) | 1.14 | 1.50
| MuSCAT3 | $z^{\prime}$ | 9486.0450 | p($a^{1}$, $y^{1}$) | 0.00152 (100 s) | 1.14 | 1.23
| SSO/Europa | $I+z$ | 9488.7749 | offset | 0.00256 (39 s) | 0.56 | 3.46
| SSO/Europa | $I+z$ | 9510.6141 | p($f^{1}$) | 0.00246 (39 s) | 0.79 | 1.42
| MuSCAT3 | $g^{\prime}$ | 9516.0739 | offset | 0.0158 (120 s) | 0.86 | 1.00
| MuSCAT3 | $r^{\prime}$ | 9516.0739 | p($a^{1}$) | 0.00558 (120 s) | 0.98 | 1.35
| MuSCAT3 | $i^{\prime}$ | 9516.0739 | offset | 0.00338 (55 s) | 1.09 | 1.35
| MuSCAT3 | $z^{\prime}$ | 9516.0739 | p($f^{1}$) | 0.00303 (29 s) | 0.92 | 1.55
| MuSCAT3 | $g^{\prime}$ | 9526.9935 | p($t^{1}$) | 0.0120 (120 s) | 1.10 | 2.00
| MuSCAT3 | $r^{\prime}$ | 9526.9935 | p($t^{1}$) | 0.00422 (120 s) | 1.08 | 1.72
| MuSCAT3 | $i^{\prime}$ | 9526.9935 | p($t^{1}$) | 0.00246 (55 s) | 0.96 | 1.08
| MuSCAT3 | $z^{\prime}$ | 9526.9935 | offset | 0.00259 (29 s) | 0.97 | 1.22
| SSO/Europa | $z^{\prime}$ | 9529.7235 | p($f^{1}$) | 0.00259 (49 s) | 0.77 | 1.63
| SSO/Europa | $I+z$ | 9548.8328 | p($t^{1}$) | 0.00283 (39 s) | 0.47 | 1.91
| SAINT-EX | $I+z$ | 9548.8328 | p($t^{2}$,$y^{1}$) | 0.00214 (128 s) | 1.66 | 1.14
| SSO/Europa | $I+z$ | 9551.5627 | p($t^{2}$) | 0.00253 (39 s) | 0.50 | 1.36
| SSO/Europa | $I+z$ | 9581.5916 | p($t^{2}$,$y^{1}$) | 0.00294 (39 s) | 0.66 | 3.01
c | TESS Sector 4 | TESS | 8412.7784 | offset | 0.0213 (120 s) | 0.93 | 1.08
| TESS Sector 4 | TESS | 8429.6934 | offset | 0.0217 (120 s) | 0.94 | 1.43
| TESS Sector 5 | TESS | 8438.1508 | offset | 0.0186 (120 s) | 0.94 | 1.89
| TESS Sector 5 | TESS | 8446.6083 | offset | 0.0160 (120 s) | 0.91 | 1.00
| TESS Sector 5 | TESS | 8455.0658 | offset | 0.0161 (120 s) | 0.92 | 1.17
| TESS Sector 5 | TESS | 8463.5232 | offset | 0.0161 (120 s) | 0.90 | 1.57
| TESS Sector 31 | TESS | 9148.5776 | offset | 0.0164 (120 s) | 0.87 | 1.43
| TESS Sector 31 | TESS | 9165.4926 | offset | 0.0168 (120 s) | 0.90 | 1.29
| TESS Sector 32 | TESS | 9182.4075 | offset | 0.0187 (120 s) | 0.95 | 1.22
| TESS Sector 32 | TESS | 9190.8650 | offset | 0.0180 (120 s) | 0.92 | 1.45
| TESS Sector 32 | TESS | 9199.3224 | offset | 0.0173 (120 s) | 0.88 | 1.43
| SSO/Io | $I+z$ | 9503.7910 | p($a^{1}$, $f^{1}$) | 0.00231 (39 s) | 0.83 | 1.47
| SSO/Europa | $I+z$ | 9520.7060 | p($a^{1}$) | 0.00232 (39 s) | 0.74 | 1.13
| SSO/Europa | $I+z$ | 9537.6209 | p($a^{1}$) | 0.00287 (39 s) | 0.43 | 1.70
| MuSCAT3 | $r^{\prime}$ | 9596.8231 | p($a^{1}$, $x^{1}$) | 0.00682 (120 s) | 1.03 | 1.30
| MuSCAT3 | $i^{\prime}$ | 9596.8231 | p($t^{1}$) | 0.00311 (55 s) | 1.06 | 1.56
| MuSCAT3 | $z^{\prime}$ | 9596.8231 | p($t^{1}$) | 0.00268 (29 s) | 0.96 | 1.38
Table 7: Photometric baseline models, residual RMS, and error scaling factors
$\beta_{w}$ and $\beta_{r}$ for each transit light curve used in our global
transit analysis (see Sect. 5.1.1). For the baseline function, p($\alpha^{n}$)
denotes a $n$-order polynomial function of the parameter $\alpha$, with
$\alpha$ that can be $t=$ time, $a=$ airmass, $x$ and $y=$ $x$\- and $y$\-
position of the target on the detector, or $f=$ full width at half maximum of
the point spread function.

##  Appendix C Transit fit posterior distributions

Figure 19: MCMC posterior distributions for the main fitted transit parameters
(see Sect. 5.1.1): the log of the orbital period (log $P$), the mid-transit
time ($T_{0}$), the transit depth ($\mathrm{d}F$), and the cosine of the
orbital inclination (cos $i_{\rm{p}}$) for each of the two planets; and the
log of the stellar density (log $\rho_{\star}$). This plot was made using the
corner python package (Foreman-Mackey 2016).