TOI-815c, we included a model where we set the atmospheric mass fraction to
zero (the no-atmosphere model). Such a model is not possible for TOI-815b, as
its radius is too large to have no atmosphere. For all models, we varied the
elemental ratios of [Fe/Si] and [Mg/Si] in the mantle and the planetary age,
and we capped the water mass fraction at 50%. As a proxy for the planet’s age,
we used the age constraints of the host star. The results are summarized in
Table 3 and Figure 11.

Figure 11: Ternary plots of the nested sampling results for the free model of TOI-815b (left) and TOI-815c (right). The axes are normalized with respect to $M_{P}-M_{\mathrm{atm}}$, and the color of the points indicates the atmospheric mass fraction in logarithmic scale. Table 3: Nested sampling results for the interior structure of the free model (top), the no-water model (middle), and the no-atmosphere model (bottom). Planet (model) | $M_{\mathrm{core}}/M_{p}$ | $M_{\mathrm{mantle}}/M_{p}$ | $M_{\mathrm{water}}/M_{p}$ | $\log(M_{\mathrm{atm}}/M_{p})$
---|---|---|---|---
free: |  |  |  |
TOI-815b | $0.29\pm 0.20$ | $0.42\pm 0.22$ | $0.28\pm 0.13$ | $-2.16\pm 0.40$
TOI-815c | $0.32\pm 0.16$ | $0.44\pm 0.26$ | $0.23\pm 0.12$ | $-7.70\pm 2.69$
no-water: |  |  |  |
TOI-815b | $0.40\pm 0.25$ | $0.58\pm 0.26$ | $0$ | $-1.64\pm 0.22$
TOI-815c | $0.30\pm 0.23$ | $0.70\pm 0.24$ | $0$ | $-2.78\pm 0.75$
no-atmosphere: |  |  |  |
TOI-815c | $0.31\pm 0.15$ | $0.43\pm 0.25$ | $0.26\pm 0.11$ | N/A

As expected, we find that a plethora of interior structures are consistent
with the observational constraints of TOI-815b and TOI-815c. These range from
water-rich structures with little to no atmosphere in the case of TOI-815c, to
structures devoid of water but with a non-negligible H-He atmosphere. However,
while the internal structure of TOI-815b is rather poorly constrained, it is
clear that the observational constraints can only be satisfied with an H-He
atmosphere. This is in alignment with Hakim et al. (2018) who calculated the
mass-radius diagram with different core and mantle compositions for planets
with Moon-like, Earth-like, or Mercury-like internal structures. Even though
they obtained a broad distribution of planetary radii for planets of the same
mass, none of the models without an atmosphere can reproduce the large radius
of TOI-815b. This further supports the conclusion that it is impossible to
construct a model of TOI-815b without a H-He atmosphere.

We emphasize that our interior structure model does not include the pollution
of the H-He envelope with heavier elements, which leads to contraction of the
atmosphere (e.g., Lozovsky et al. 2018). Furthermore, water might not only be
present in an ocean layer, but also dissolved within the mantle of the planet
(e.g., Dorn & Lichtenberg, 2021) or be present in the atmosphere (e.g., Mol
Lous et al., 2022). In addition, the assumption of distinct planetary layers
might break down at high planetary masses (Helled & Stevenson, 2017;
Bodenheimer et al., 2018). Nonetheless, in face of the observational
constraints currently available, these additions are unlikely to affect the
main conclusions of this analysis.

###  5.4 Planetary atmospheric evolution

To model the atmospheric evolution of the planets in the TOI-815 system and
constrain their primordial parameters, we present modeling conducted using the
Planetary Atmospheres and Stellar RoTation RAtes (PASTA; Bonfanti et al.,
2021b) code, which is a planetary atmospheric evolution code based on the
original code presented by Kubyshkina et al. (2019a, b). As stars can be born
with different initial rotation periods, their $L_{\rm XUV}$ values can spread
over a wide range during the first $\sim 1$ Gyr following protoplanetary disk
dispersal (e.g., Tu et al., 2015). To account for that, PASTA simultaneously
constrains the evolution of planetary atmospheres and of the stellar rotation
rate by combining a model predicting planetary atmospheric escape rates based
on hydrodynamic simulations (this has the advantage over other commonly used
analytical estimates in that it accounts for both extreme-UV-driven and core-
powered mass loss; Kubyshkina et al., 2018), a model of the evolution of the
stellar extreme-UV (XUV) flux (Bonfanti et al., 2021b), a model relating
planetary parameters and atmospheric mass (Johnstone et al., 2015b), and
stellar evolutionary tracks (Choi et al., 2016).

PASTA works under two main assumptions: (1) planet migration did not occur
after the dispersal of the protoplanetary disk; and (2) the planets hosted at
some point in the past or still host a hydrogen-dominated atmosphere. The free
parameters (i.e., subject to uniform priors) are the planetary initial
atmospheric mass fractions at the time of the dispersal of the protoplanetary
disk ($f_{\rm atm}^{\rm start}$), which we assumed occurs at an age of 5 Myr
(see for example Alexander et al., 2014; Kimura et al., 2016; Gorti et al.,
2016), and the stellar rotation period at 150 Myr, which is used as a proxy
for the stellar XUV emission. The code returns constraints on the free
parameters and on their uncertainties by implementing the atmospheric
evolution algorithm in a Bayesian framework (through the MC3 code of Cubillos
et al. 2017), using the system parameters with their uncertainties as input
priors. Details of the algorithm can be found in Bonfanti et al. (2021b).

As a proxy for the evolution of the stellar rotation period, and thus XUV
emission, Figure 12 displays the posterior distribution of the stellar
rotation period at an age of 150 Myr ($P_{\mathrm{rot,150}}$). This
distribution is then compared to that of stars member of young open clusters
that are of comparable mass, as obtained from Johnstone et al. (2015a).
Because the system is still very young ($200^{+400}_{-200}$ Myr), the
posterior distribution of $P_{\mathrm{rot,150}}$ is mostly dominated by the
measurement of the present-day rotation rate ($15.3\leavevmode\nobreak\
\pm\leavevmode\nobreak\ 1.2$ days). The distribution is, however, flatter than
the prior on $P_{\mathrm{rot}}$, because of the large relative uncertainty of
the stellar age. The result indicates that the star could not be a fast
rotator (i.e., $P_{\mathrm{rot,150}}<1\,\rm days$), as constrained by the slow
present-day rotation period.

Figure 13 shows the posterior distribution of the initial atmospheric mass
fraction for TOI-815b (left panel) and TOI-815c (right panel) in comparison to
the present-day atmospheric mass fraction. For planet b we find that PASTA
clearly prefers evolutionary scenarios in which the evolution of the
atmosphere was mostly unaffected by hydrodynamic escape. This indicates that
the planet most likely still retains the majority of its primordial
atmosphere. However, it must be noted that given the uncertainties of the
observed present-day parameters PASTA is not able to completely exclude
scenarios with significant escape.

In the case of TOI-815c, PASTA concludes that, independently of the evolution
of the stellar rotation rate, the planet has not lost a significant amount of
atmosphere via hydrodynamic escape. Given its present-day orbital position and
relatively large core mass, it would be expected for the planet to have
accreted an extensive hydrogen atmosphere if it had formed within the
protoplanetary disk (Ikoma & Hori, 2012; Lee & Chiang, 2015; Mordasini, 2020;
Venturini & Helled, 2017). However, such an envelope is not observed today and
according to PASTA the planet is unable to have lost it due to hydrodynamic
escape over the course of its evolution. One possible explanation would be
formation within a gas-poor environment (e.g., Lee et al. 2022). Another
explanation is that the planet suffered from at least one giant impact toward
the end of the disk’s life (Inamdar & Schlichting, 2015). This also aligns
with the fact that the planets of the system are not in mean motion resonance
and the expected spin-orbit misalignment. In both cases the results retrieved
by PASTA are compatible with conclusions that can be drawn by the more
simplistic approach of only evaluating the restricted Jeans parameter
$\Lambda$ (Fossati et al., 2017). For TOI-815c the high value of
$\Lambda=145^{+16}_{-16}$ would imply a stable atmosphere, while for TOI-815b
the intermediate value of $\Lambda=28.6^{+5.7}_{-5.3}$ would favor a stable
atmosphere but not exclude the possibility of a hydrodynamically escaping
atmosphere.

Figure 12: Posterior distribution (dark blue) of the stellar rotation rate of
TOI-815 after 150 Myr derived by PASTA. The solid black line represents the
distribution of the stellar rotation rate of young open cluster stars with
masses comparable to that of TOI-815 based on the collection of data provided
by Johnstone et al. (2015a). The dashed black line indicates the observed
present-day stellar rotation rate.

We note that in principle PASTA is not set up to simulate the atmospheric
evolution of planets with mean planetary densities larger than that of Earth.
This is due to the fact that the code always assumes for the rocky core of a
planet to have a mean density equal to that of the Earth. In practice, PASTA
will underestimate the planetary mass of TOI-815c as its observed density is
$1.30\pm 0.20\,\rho_{\mathrm{\oplus}}$. The planetary mass estimate for planet
c returned by PASTA is $M_{\mathrm{P,c}}=19.9\pm 1.9\,M_{\mathrm{\oplus}}$,
which implies a bare rocky-core with a mean planetary density of
$1\,\rho_{\mathrm{\oplus}}$. A larger planetary mass, however, would just
further hinder hydrodynamic escape, as escaping particles would need even more
energy to escape the gravitational well of the planet. Therefore, since
hydrodynamic escape affects only marginally the evolution of the planetary
atmosphere, the impact of the assumption on the core density has a negligible
impact on the results.

Figure 13: Posterior distribution (dark blue) for the mass of the planetary
atmosphere of TOI-815b (left) and TOI-815c (right) at the time of the
dispersal of the protoplanetary disk. The light blue line represents the
distribution of the estimated present-day atmospheric mass fraction derived by
PASTA. The orange horizontal line indicates the uninformative prior
distribution.

###  5.5 Search for transit timing variations

Although the ratio of the orbital period of TOI-815 and c
($P_{\mathrm{orb,c}}$/$P_{\mathrm{orb,b}}$ $\simeq$ 3.12) lies exterior to the
3:1 resonance, we conducted a TTV analysis to search for additional, non-
transiting planets. Using juliet, we performed the analysis by utilizing all
available photometric datasets from TESS, CHEOPS, ASTEP, and LCOGT. Instead of
fitting a single period ($P_{\mathrm{orb}}$) and time-of-transit center
($T_{0}$), juliet employs a method that seeks individual transit times. This
approach involves fitting each transit independently and determining one
transit time for each one, resulting in a more consistent and coherent
analysis. Figure 14 illustrates the results of our analysis, which compares
the observed transit times with the calculated linear ephemeris derived from
all the transits. Notably, no significant variations were detected within the
data.

Figure 14: Transit timings for each transit of TOI-815b (purple) and TOI-815c
(orange) measured from the TESS and CHEOPS light curves. The purple and orange
regions represent the 1$\sigma$ uncertainties of the inner and outer planet
ephemeris reported in Table 5.

###  5.6 Potential for atmospheric and Rossiter-McLaughlin studies

We used the metric proposed by Kempton et al. (2018) to assess the potential
of atmospheric characterization for the planets in the TOI-815 system with the
James Webb Space Telescope (JWST; Gardner et al. 2023). To calculate the
transmission spectroscopy metric for all exoplanets in the PlanetS catalog2
with a radii from 2 to 4 $\it{R_{\rm\mathrm{\oplus}}}$, we utilized the scale
factors provided in Table 1 by Kempton et al. (2018) instead of the suggested
value of 0.167 for temperate planets. Notably, the transmission spectroscopy
metric values for TOI-815b and c were found to be 97 $\pm$ 20 and 15 $\pm$ 2,
respectively. Putting these values into context, TOI-815b is located in the
top 10$\%$ of well-characterized exoplanets in this regime and therefore
stands out as an excellent candidate for atmospheric characterization via
transmission spectroscopy.

We also explored the possibility of studying the spin–orbit of both planets
via Rossiter-McLaughlin (RM) observations (Rossiter 1924; McLaughlin 1924).
According to Eq. (40) by Winn (2010), we expect the maximum amplitude of the
RM effect to be $\sim$1.13 m s-1 and 0.52 m s-1 for TOI-815b and TOI-815c,
respectively. Using the classical RM technique, the signals are small and
challenging, but the sky-projected spin–orbit angle $\lambda$ should be
detectable using newer methods such as the RM effect revolutions technique
(Bourrier et al. 2021). Moreover, since the star is seen close to pole-on and
since the planets are transiting, the two orbits are probably misaligned. This
can be seen using the following equation linking the 3D spin–orbit angle
$\psi$ to other observable quantities (e.g., Winn et al., 2007):

$\cos\psi=\sin i_{\star}\sin i\cos\lambda+\cos i_{\star}\cos i,$ (4)

yielding a small value for $\cos\psi$ for any value of $\lambda$ we measure.
The case that would disfavor the orbital misalignment picture the most is if
we observe an aligned sky-projected spin–orbit angle $\lambda\sim 0^{\circ}$.
Sampling a $\psi$ distribution from Eq. (4) with a Dirac distribution at
$0^{\circ}$ for $\lambda$ and measurement-informed Gaussians for the other
parameters, we find $\psi$ = 58${}^{+18}_{-19}$, which still represents a
misaligned orbit at a 3.2$\sigma$ level. Given the young nature of the system,
tidal realignment (e.g., Barker & Ogilvie, 2010; Attia et al., 2023) may have
not had time to occur. Hence, TOI-815 offers us a unique possibility to
capture a peculiar orbital architecture, especially given the scarcity of
misaligned orbits in multi-planet systems (e.g., Albrecht et al., 2022).

##  6 Conclusions

We have reported the detection and characterization of two sub-Neptune planets
that are transiting TOI-815. TOI-815 is a young ($200^{+400}_{-200}$ Myr) K
dwarf in a K-M binary system. The planets were initially detected by the TESS
mission and subsequently confirmed with CHEOPS observations as well as precise
RV measurements obtained with ESPRESSO, which enabled the determination of
their masses. Through the combined use of TESS, CHEOPS, LCOGT, and ASTEP
photometry and the ESPRESSO spectroscopic data, we precisely characterized
both planets: TOI-815b has a period of 11.2 days, a radius of 2.94
$\it{R_{\rm\mathrm{\oplus}}}$, and a mass of 7.6
$\it{M_{\rm\mathrm{\oplus}}}$, while TOI-815c has a period of 35 days, a
radius of 2.62 $\it{R_{\rm\mathrm{\oplus}}}$, and a mass of 23.5
$\it{M_{\rm\mathrm{\oplus}}}$. The combination of these facilities allowed us
to determine the planetary radius and mass with precision better than 2$\%$
and 20$\%$ for the inner planet and 4$\%$ and 10$\%$ for the outer planet,
respectively.

TOI-815c, with a bulk density of 7.2 $\pm$ 1.1 g cm-3, stands as the densest
well-characterized warm mini-Neptune. Our study reveals a statistically
significant difference between the bulk densities of warm sub-Neptunes and hot
sub-Neptunes, with TOI-815b and c being instrumental in establishing a 95%
confidence. Further exoplanet detections with precise mass measurements are
needed to confirm this trend. Our analysis of the TOI-815 system using
similarity metrics demonstrates that while the planets in this system have
remarkably similar radii, they exhibit significantly different masses, making
TOI-815 one of the least uniform systems in terms of planetary density.

Internal structure modeling that includes four layers (an iron core, a rocky
mantle, a water layer, and a H-He atmosphere) suggests that TOI-815b likely
has a small H-He atmosphere (on the order of a percent of the planet’s mass)
if water is included. If no water is assumed, the atmospheric mass fraction
can be an order of magnitude higher. On the contrary, the radius of TOI-815c
can be reproduced without any atmosphere. This object is interesting because
it is rather massive and yet seems to have only a tiny (or no) atmosphere.
This is a challenge for planet formation theories since at these planetary
masses strong accretion of H-He is expected (for comparison, Uranus and
Neptune have lower masses but consist of $\sim$10% of H-He). A possible
explanation would be for the planet to have lost the majority of its
primordial atmosphere. However, our atmospheric evolution model shows that
such an extensive mass loss cannot be caused by XUV-driven hydrodynamic
escape. Thus, the most likely scenario is that planet c, similar to
Kepler-107c (Bonomo et al. 2019), suffered from at least one giant impact
(which removed a substantial amount of primordial gas) during the unstable
dynamical phase of disk dispersal. This is consistent with the planets of this
system not being in resonance. Alternatively, a formation within a gas-poor
environment cannot be ruled out either.

Finally, the combination of the stellar rotation period ($P_{\rm rot}$ = 15.3
$\pm$ 1.2 days) constrained from the WASP and TESS photometry and the very low
spectroscopically measured projected rotational velocity ($v\sin i_{*}$ $<$ 1
km s-1) suggest that we are observing the star close to pole-on. Given that
the planets are on orbits with inclinations close to $\sim$90 $\degree$, we
expect a spin-orbit misalignment at a 3.2$\sigma$ level. Future RM
observations will be able to confirm this.

###### Acknowledgements.

We thank the Swiss National Science Foundation (SNSF) and the Geneva
University for their continuous support to our planet low-mass companion
search programs. This work has been carried out within the framework of the
National Centre of Competence in Research PlanetS supported by the Swiss
National Science Foundation under grants 51NF40$\\_$182901 and
51NF40$\\_$205606\. OA and VB have received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (project Spice Dune, grant agreement No 947634; grant agreement No
730890). This material reflects only the authors’ views and the Commission is
not liable for any use that may be made of the information contained therein.
The authors acknowledge the financial support of the SNSF. This publication
makes use of The Data & Analysis Center for Exoplanets (DACE), which is a
facility based at the University of Geneva (CH) dedicated to extrasolar planet
data visualization, exchange, and analysis. DACE is a platform of NCCR
$PlanetS$ and is available at https://dace.unige.ch. This paper includes data
collected by the TESS mission. Funding for the TESS mission is provided by the
NASA Explorer Program. JSJ acknowledges support by FONDECYT grant 1201371 and
partial support from the ANID Basal project FB210003. Funding for the TESS
mission is provided by NASA’s Science Mission Directorate. KAC and CNW
acknowledge support from the TESS mission via subaward s3449 from MIT. This
research has made use of the Exoplanet Follow-up Observation Program (ExoFOP;
DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute
of Technology, under contract with the National Aeronautics and Space
Administration under the Exoplanet Exploration Program. This paper made use of
data collected by the TESS mission and are publicly available from the
Mikulski Archive for Space Telescopes (MAST) operated by the Space Telescope
Science Institute (STScI). We acknowledge the use of public TESS data from
pipelines at the TESS Science Office and at the TESS Science Processing
Operations Center. 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. CHEOPS is an ESA mission in partnership with Switzerland with
important contributions to the payload and the ground segment from Austria,
Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the
United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the
support received by all the agencies, offices, universities, and industries
involved. Their flexibility and willingness to explore new approaches were
essential to the success of this mission. CHEOPS data analyzed in this article
will be made available in the CHEOPS mission archive
(https://cheops.unige.ch/archive_browser/). This work makes use of
observations from the ASTEP telescope. ASTEP benefited from the support of the
French and Italian polar agencies IPEV and PNRA in the framework of the
Concordia station program, from OCA, INSU, Idex UCAJEDI (ANR- 15-IDEX-01) and
ESA through the Science Faculty of the European Space Research and Technology
Centre (ESTEC). This work makes use of observations from the LCOGT network.
Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale
Innovations Program (MSIP). MSIP is funded by NSF. Co-funded by the European
Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those
of the author(s) only and do not necessarily reflect those of the European
Union or the European Research Council. Neither the European Union nor the
granting authority can be held responsible for them. C.B. acknowledges support
from the Swiss Space Office through the ESA PRODEX program. This work has been
carried out within the framework of the NCCR PlanetS supported by the Swiss
National Science Foundation under grants 51NF40_182901 and 51NF40_205606. ML
acknowledges support of the Swiss National Science Foundation under grant
number PCEFP2_194576. S.G.S. acknowledge support from FCT through FCT contract
nr. CEECIND/00826/2018 and POPH/FSE (EC). YAl acknowledges support from the
Swiss National Science Foundation (SNSF) under grant 200020_192038. ABr was
supported by the SNSA. DG gratefully acknowledges financial support from the
CRT foundation under Grant No. 2018.2323 “Gaseousor rocky? Unveiling the
nature of small worlds”. JV acknowledges support from the Swiss National
Science Foundation (SNSF) under grant PZ00P2_208945. MNG is the ESA CHEOPS
Project Scientist and Mission Representative, and as such also responsible for
the Guest Observers (GO) Programme. MNG does not relay proprietary information
between the GO and Guaranteed Time Observation (GTO) Programmes, and does not
decide on the definition and target selection of the GTO Programme. KGI is the
ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest
Observers Programme. She does not participate in, or contribute to, the
definition of the Guaranteed Time Programme of the CHEOPS mission through
which observations described in this paper have been taken, nor to any aspect
of target selection for the programme. We acknowledge financial support from
the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación
MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through
projects PID2019-107061GB-C61, PID2019-107061GB-C66, PID2021-125627OB-C31, and
PID2021-125627OB-C32, from the Centre of Excellence “Severo Ochoa” award to
the Instituto de Astrofísica de Canarias (CEX2019-000920-S), from the Centre
of Excellence “María de Maeztu” award to the Institut de Ciències de l’Espai
(CEX2020-001058-M), and from the Generalitat de Catalunya/CERCA programme.
S.C.C.B. acknowledges support from FCT through FCT contracts nr.
IF/01312/2014/CP1215/CT0004. XB, SC, DG, MF and JL acknowledge their role as
ESA-appointed CHEOPS science team members. LBo, GBr, VNa, IPa, GPi, RRa, GSc,
VSi, and TZi acknowledge support from CHEOPS ASI-INAF agreement n.
2019-29-HH.0. ACC acknowledges support from STFC consolidated grant numbers
ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. P.E.C. is
funded by the Austrian Science Fund (FWF) Erwin Schroedinger Fellowship,
program J4595-N. This project was supported by the CNES. The Belgian
participation to CHEOPS has been supported by the Belgian Federal Science
Policy Office (BELSPO) in the framework of the PRODEX Program, and by the
University of Liège through an ARC grant for Concerted Research Actions
financed by the Wallonia-Brussels Federation; L.D. is an F.R.S.-FNRS
Postdoctoral Researcher. This work was supported by FCT - Fundação para a
Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020
- Programa Operacional Competitividade e Internacionalização by these grants:
UIDB/04434/2020; UIDP/04434/2020; 2022.06962.PTDC. O.D.S.D. is supported in
the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds
through FCT. S.G.S acknowledges the support from FCT through Investigador FCT
contract nr. CEECIND/00826/2018 and POPH/FSE (EC). B.-O. D. acknowledges
support from the Swiss State Secretariat for Education, Research and
Innovation (SERI) under contract number MB22.00046. This project has received
funding from the Swiss National Science Foundation for project 200021_200726.
It has also been carried out within the framework of the National Centre of
Competence in Research PlanetS supported by the Swiss National Science
Foundation under grant 51NF40_205606. The authors acknowledge the financial
support of the SNSF. MF and CMP gratefully acknowledge the support of the
Swedish National Space Agency (DNR 65/19, 174/18). T.D. acknowledges support
by the McDonnell Center for the Space Sciences at Washington University in St.
Louis. M.G. is an F.R.S.-FNRS Senior Research Associate. CHe acknowledges
support from the European Union H2020-MSCA-ITN-2019 under Grant Agreement no.
860470 (CHAMELEON). SH gratefully acknowledges CNES funding through the grant
837319. K.W.F.L. was supported by Deutsche Forschungsgemeinschaft grants
RA714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of
Extrasolar Planets. This work was granted access to the HPC resources of
MesoPSL financed by the Region Ile de France and the project Equip@Meso
(reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir
supervised by the Agence Nationale pour la Recherche. PM acknowledges support
from STFC research grant number ST/M001040/1. This work was also partially
supported by a grant from the Simons Foundation (PI Queloz, grant number
327127). A. S. acknowledges support from the Swiss Space Office through the
ESA PRODEX program. TWi acknowledges support from the UKSA and the University
of Warwick. NAW acknowledges UKSA grant ST/R004838/1. V.V.G. is an F.R.S-FNRS
Research Associate. C.D. acknowledges support from the Swiss National Science
Foundation under grant TMSGI2_211313. GyMSz acknowledges the support of the
Hungarian National Research, Development and Innovation Office (NKFIH) grant
K-125015, a a PRODEX Experiment Agreement No. 4000137122, the Lendület
LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the
city of Szombathely. DR was supported by NASA under award number NNA16BD14C
for NASA Academic Mission Services. R.L. acknowledges funding from University
of La Laguna through the Margarita Salas Fellowship from the Spanish Ministry
of Universities ref. UNI/551/2021-May 26, and under the EU Next Generation
funds. R. A. is a Trottier Postdoctoral Fellow and acknowledges support from
the Trottier Family Foundation. This work was supported in part through a
grant from the Fonds de Recherche du Québec - Nature et Technologies (FRQNT).
This work was funded by the Institut Trottier de Recherche sur les Exoplanétes
(iREx). This work made use of tpfplotter by J. Lillo-Box (publicly available
in www.github.com/jlillo/tpfplotter), which also made use of the python
packages astropy, lightkurve, matplotlib and numpy.

## References

  * Abe et al. (2013) Abe, L., Gonçalves, I., Agabi, A., et al. 2013, A&A, 553, A49
  * Adibekyan et al. (2015) Adibekyan, V., Figueira, P., Santos, N. C., et al. 2015, A&A, 583, A94
  * Adibekyan et al. (2012) Adibekyan, V. Z., Sousa, S. G., Santos, N. C., et al. 2012, A&A, 545, A32
  * Aguichine et al. (2021) Aguichine, A., Mousis, O., Deleuil, M., & Marcq, E. 2021, ApJ, 914, 84
  * Akeson et al. (2013) Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989
  * Albrecht et al. (2022) Albrecht, S. H., Dawson, R. I., & Winn, J. N. 2022, PASP, 134, 082001
  * Alexander et al. (2014) Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475–496
  * Aller et al. (2020) Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128
  * Almeida-Fernandes & Rocha-Pinto (2018) Almeida-Fernandes, F. & Rocha-Pinto, H. J. 2018, MNRAS, 476, 184
  * Anthony-Twarog et al. (2018) Anthony-Twarog, B. J., Deliyannis, C. P., Harmer, D., et al. 2018, AJ, 156, 37
  * Attia et al. (2023) Attia, O., Bourrier, V., Delisle, J. B., & Eggenberger, P. 2023, A&A, 674, A120
  * Attia et al. (2021) Attia, O., Bourrier, V., Eggenberger, P., et al. 2021, A&A, 647, A40
  * Barker & Ogilvie (2010) Barker, A. J. & Ogilvie, G. I. 2010, MNRAS, 404, 1849
  * Barnes (2007) Barnes, S. A. 2007, ApJ, 669, 1167
  * Behmard et al. (2022) Behmard, A., Dai, F., & Howard, A. W. 2022, AJ, 163, 160
  * Benz et al. (2021) Benz, W., Broeg, C., Fortier, A., et al. 2021, Experimental Astronomy, 51, 109
  * Berger et al. (2020) Berger, T. A., Huber, D., Gaidos, E., van Saders, J. L., & Weiss, L. M. 2020, AJ, 160, 108
  * Blackwell & Shallis (1977) Blackwell, D. E. & Shallis, M. J. 1977, MNRAS, 180, 177
  * Bodenheimer et al. (2018) Bodenheimer, P., Stevenson, D. J., Lissauer, J. J., & D’Angelo, G. 2018, The Astrophysical Journal, 868, 138
  * Bonfanti et al. (2021a) Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021a, A&A, 646, A157
  * Bonfanti et al. (2021b) Bonfanti, A., Fossati, L., Kubyshkina, D., & Cubillos, P. E. 2021b, A&A, 656, A157
  * Bonfanti et al. (2016) Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5
  * Bonfanti et al. (2015) Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18
  * Bonomo et al. (2023) Bonomo, A. S., Dumusque, X., Massa, A., et al. 2023, A&A, 677, A33
  * Bonomo et al. (2019) Bonomo, A. S., Zeng, L., Damasso, M., et al. 2019, Nature Astronomy, 3, 416
  * Borucki & Summers (1984) Borucki, W. J. & Summers, A. L. 1984, Icarus, 58, 121
  * Bourrier et al. (2021) Bourrier, V., Lovis, C., Cretignier, M., et al. 2021, A&A, 654, A152
  * Brandeker et al. (2022) Brandeker, A., Heng, K., Lendl, M., et al. 2022, A&A, 659, L4
  * Brown et al. (2013) Brown, T. M., Baliber, N., Bianco, F. B., et al. 2013, PASP, 125, 1031
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Casagrande et al. (2011) Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138
  * Castelli & Kurucz (2003) Castelli, F. & Kurucz, R. L. 2003, in IAU Symposium, Vol. 210, Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, A20
  * Chakraborty et al. (2023) Chakraborty, H., Lendl, M., Akinsanmi, B., Petit dit de la Roche, D. J. M., & Deline, A. 2023, arXiv e-prints, arXiv:2311.16864
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Cochran et al. (2011) Cochran, W. D., Fabrycky, D. C., Torres, G., et al. 2011, ApJS, 197, 7
  * Collins (2019) Collins, K. 2019, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 140.05
  * Collins et al. (2017) Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77
  * Cubillos et al. (2017) Cubillos, P., Erkaev, N. V., Juvan, I., et al. 2017, MNRAS, 466, 1868
  * Cummings et al. (2017) Cummings, J. D., Deliyannis, C. P., Maderak, R. M., & Steinhauer, A. 2017, AJ, 153, 128
  * Delgado Mena et al. (2014) Delgado Mena, E., Israelian, G., González Hernández, J. I., et al. 2014, A&A, 562, A92
  * Deline et al. (2022) Deline, A., Hooton, M. J., Lendl, M., et al. 2022, A&A, 659, A74
  * Delrez et al. (2021) Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nature Astronomy, 5, 775
  * Dorn et al. (2015) Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83
  * Dorn & Lichtenberg (2021) Dorn, C. & Lichtenberg, T. 2021, The Astrophysical Journal Letters, 922, L4
  * Dorn et al. (2017) Dorn, C., Venturini, J., Khan, A., et al. 2017, Astronomy and Astrophysics, 597, A37
  * Doyle et al. (1984) Doyle, L. R., Wilcox, T. J., & Lorre, J. J. 1984, ApJ, 287, 307
  * Dragomir et al. (2019) Dragomir, D., Teske, J., Günther, M. N., et al. 2019, ApJ, 875, L7
  * Dransfield et al. (2022) Dransfield, G., Mékarnia, D., Triaud, A. H. M. J., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12186, Observatory Operations: Strategies, Processes, and Systems IX, ed. D. S. Adler, R. L. Seaman, & C. R. Benn, 121861F
  * Espinoza et al. (2019) Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262
  * Fausnaugh et al. (2020) Fausnaugh, M. M., Burke, C. J., Ricker, G. R., & Vanderspek, R. 2020, Research Notes of the AAS, 4, 251
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fossati et al. (2017) Fossati, L., Erkaev, N. V., Lammer, H., et al. 2017, A&A, 598, A90
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1
  * Garai et al. (2023) Garai, Z., Osborn, H. P., Gandolfi, D., et al. 2023, A&A, 674, A44
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001
  * Ginzburg et al. (2018) Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759
  * Gorti et al. (2016) Gorti, U., Liseau, R., Sándor, Z., & Clarke, C. 2016, Space Sci. Rev., 205, 125
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39
  * Guillot et al. (2015) Guillot, T., Abe, L., Agabi, A., et al. 2015, Astronomische Nachrichten, 336, 638
  * Gupta & Schlichting (2019) Gupta, A. & Schlichting, H. E. 2019, MNRAS, 487, 24
  * Hadden & Lithwick (2017) Hadden, S. & Lithwick, Y. 2017, AJ, 154, 5
  * Hakim et al. (2018) Hakim, K., Rivoldini, A., Van Hoolst, T., et al. 2018, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, 15447
  * Helled & Stevenson (2017) Helled, R. & Stevenson, D. 2017, ApJ, 840, L4
  * Henry et al. (1996) Henry, T. J., Soderblom, D. R., Donahue, R. A., & Baliunas, S. L. 1996, AJ, 111, 439
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Hoyer et al. (2020) Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24
  * Huang et al. (2020) Huang, C. X., Vanderburg, A., Pál, A., et al. 2020, Research Notes of the American Astronomical Society, 4, 206
  * Ikoma & Hori (2012) Ikoma, M. & Hori, Y. 2012, ApJ, 753, 66
  * Inamdar & Schlichting (2015) Inamdar, N. K. & Schlichting, H. E. 2015, MNRAS, 448, 1751
  * Jenkins (2002) Jenkins, J. M. 2002, ApJ, 575, 493
  * Jenkins et al. (2010) Jenkins, J. M., Chandrasekaran, H., McCauliff, S. D., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7740, Software and Cyberinfrastructure for Astronomy, ed. N. M. Radziwill & A. Bridger, 77400D
  * 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 Proc. SPIE, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, 99133E
  * Jensen (2013) Jensen, E. 2013, Tapir: A web interface for transit/eclipse observability, Astrophysics Source Code Library
  * Jin & Mordasini (2018) Jin, S. & Mordasini, C. 2018, ApJ, 853, 163
  * Johnstone et al. (2015a) Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015a, A&A, 577, A28
  * Johnstone et al. (2015b) Johnstone, C. P., Güdel, M., Stökl, A., et al. 2015b, ApJ, 815, L12
  * Kempton et al. (2018) Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401
  * Kimura et al. (2016) Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016, MNRAS, 461, 2257
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161
  * Krenn et al. (2023) Krenn, A. F., Lendl, M., Patel, J. A., et al. 2023, A&A, 672, A24
  * Kubyshkina et al. (2019a) Kubyshkina, D., Cubillos, P. E., Fossati, L., et al. 2019a, ApJ, 879, 26
  * Kubyshkina et al. (2018) Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018, A&A, 619, A151
  * Kubyshkina et al. (2019b) Kubyshkina, D., Fossati, L., Mustill, A. J., et al. 2019b, A&A, 632, A65
  * Kurucz (1993a) Kurucz, R. 1993a, SYNTHE Spectrum Synthesis Programs and Line Data. Kurucz CD-ROM No. 18. Cambridge, 18
  * Kurucz (1993b) Kurucz, R. L. 1993b, SYNTHE spectrum synthesis programs and line data (Astrophysics Source Code Library)
  * Lee & Chiang (2015) Lee, E. J. & Chiang, E. 2015, ApJ, 811, 41
  * Lee et al. (2022) Lee, E. J., Karalis, A., & Thorngren, D. P. 2022, ApJ, 941, 186
  * Leleu et al. (2021) Leleu, A., Alibert, Y., Hara, N. C., et al. 2021, Astronomy and Astrophysics, 649, A26
  * Lendl et al. (2020) Lendl, M., Bouchy, F., Gill, S., et al. 2020, MNRAS, 492, 1761
  * Li et al. (2019) Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4
  * Lopez & Fortney (2013) Lopez, E. D. & Fortney, J. J. 2013, ApJ, 776, 2
  * Lozovsky et al. (2018) Lozovsky, M., Helled, R., Dorn, C., & Venturini, J. 2018, ApJ, 866, 49
  * Lucy & Sweeney (1971) Lucy, L. B. & Sweeney, M. A. 1971, AJ, 76, 544
  * Luque et al. (2019) Luque, R., Nowak, G., Pallé, E., et al. 2019, A&A, 623, A114
  * Maciel et al. (2011) Maciel, W. J., Rodrigues, T. S., & Costa, R. D. D. 2011, A&A, 47, 401
  * Marcy et al. (2014) Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  * Masuda & Winn (2020) Masuda, K. & Winn, J. N. 2020, AJ, 159, 81
  * Maxted et al. (2022) Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77
  * 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, Proc. SPIE, 107070K
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22
  * Mékarnia et al. (2016) Mékarnia, D., Guillot, T., Rivet, J. P., et al. 2016, MNRAS, 463, 45
  * Millholland et al. (2017) Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33
  * Mishra et al. (2023) Mishra, L., Alibert, Y., Udry, S., & Mordasini, C. 2023, A&A, 670, A68
  * Mol Lous et al. (2022) Mol Lous, M., Helled, R., & Mordasini, C. 2022, Nature Astronomy, 6, 819
  * Mordasini (2020) Mordasini, C. 2020, A&A, 638, A52
  * Morris et al. (2020) Morris, R. L., Twicken, J. D., Smith, J. C., et al. 2020, Kepler Data Processing Handbook: Photometric Analysis, Kepler Science Document KSCI-19081-003
  * Naponiello et al. (2023) Naponiello, L., Mancini, L., Sozzetti, A., et al. 2023, Nature, 622, 255
  * Nordström et al. (2004) Nordström, B., Mayor, M., Andersen, J., et al. 2004, A&A, 418, 989
  * Osborn (2022) Osborn, H. P. 2022, MonoTools: Planets of uncertain periods detector and modeler, Astrophysics Source Code Library, record ascl:2204.020
  * Osborn et al. (2023) Osborn, H. P., Nowak, G., Hébrard, G., et al. 2023, MNRAS, 523, 3069
  * Otegi et al. (2020) Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43
  * Otegi et al. (2022) Otegi, J. F., Helled, R., & Bouchy, F. 2022, A&A, 658, A107
  * Owen & Wu (2017) Owen, J. E. & Wu, Y. 2017, ApJ, 847, 29
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96
  * Petigura et al. (2022) Petigura, E. A., Rogers, J. G., Isaacson, H., et al. 2022, AJ, 163, 179
  * Pollacco et al. (2006) Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rogers & Owen (2021) Rogers, J. G. & Owen, J. E. 2021, MNRAS, 503, 1526
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15
  * Santos et al. (2013) Santos, N. C., Sousa, S. G., Mortier, A., et al. 2013, A&A, 556, A150
  * Schanche et al. (2020) Schanche, N., Hébrard, G., Collier Cameron, A., et al. 2020, MNRAS, 499, 428
  * Schmider et al. (2022) Schmider, F.-X., Abe, L., Agabi, A., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12182, Ground-based and Airborne Telescopes IX, ed. H. K. Marshall, J. Spyromilio, & T. Usuda, 121822O
  * Sestito & Randich (2005) Sestito, P. & Randich, S. 2005, A&A, 442, 615
  * 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
  * Sneden (1973) Sneden, C. A. 1973, PhD thesis, THE UNIVERSITY OF TEXAS AT AUSTIN.
  * Sousa (2014) Sousa, S. G. 2014, in Determination of Atmospheric Parameters of B, 297–310
  * Sousa et al. (2021) Sousa, S. G., Adibekyan, V., Delgado-Mena, E., et al. 2021, A&A, 656, A53
  * Sousa et al. (2015) Sousa, S. G., Santos, N. C., Adibekyan, V., Delgado-Mena, E., & Israelian, G. 2015, A&A, 577, A67
  * Sousa et al. (2007) Sousa, S. G., Santos, N. C., Israelian, G., Mayor, M., & Monteiro, M. J. P. F. G. 2007, A&A, 469, 783
  * Sousa et al. (2008) Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Stassun et al. (2019) Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138
  * 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
  * Sun et al. (2019) Sun, L., Ioannidis, P., Gu, S., et al. 2019, A&A, 624, A15
  * Sun et al. (2023) Sun, Q., Deliyannis, C. P., Steinhauer, A., Anthony-Twarog, B. J., & Twarog, B. A. 2023, ApJ, 952, 71
  * Tokovinin (2018) Tokovinin, A. 2018, PASP, 130, 035002
  * Tsantaki et al. (2013) Tsantaki, M., Sousa, S. G., Adibekyan, V. Z., et al. 2013, A&A, 555, A150
  * Tu et al. (2015) Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3
  * Twicken et al. (2018) Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502
  * Twicken et al. (2010) Twicken, J. D., Clarke, B. D., Bryson, S. T., et al. 2010, in Proc. SPIE, Vol. 7740, Software and Cyberinfrastructure for Astronomy, 774023
  * Ulmer-Moll et al. (2023) Ulmer-Moll, S., Osborn, H. P., Tuson, A., et al. 2023, A&A, 674, A43
  * Vanderburg & Johnson (2014) Vanderburg, A. & Johnson, J. A. 2014, PASP, 126, 948
  * Venturini et al. (2020) Venturini, J., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, A&A, 643, L1
  * Venturini & Helled (2017) Venturini, J. & Helled, R. 2017, ApJ, 848, 95
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48
  * Wielen (1977) Wielen, R. 1977, A&A, 60, 263
  * Winn (2010) Winn, J. N. 2010, in Exoplanets, ed. S. Seager, 55–77
  * Winn et al. (2007) Winn, J. N., Holman, M. J., Henry, G. W., et al. 2007, AJ, 133, 1828
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Zechmeister & Kürster (2009) Zechmeister, M. & Kürster, M. 2009, A&A, 496, 577
  * Zeng et al. (2016) Zeng, L., Sasselov, D. D., & Jacobsen, S. B. 2016, ApJ, 819, 127
  * Ziegler et al. (2020) Ziegler, C., Tokovinin, A., Briceno, C., et al. 2020, VizieR Online Data Catalog, J/AJ/159/19

##  Appendix A TESS light curves

Figure 15: TESS SAP flux light curves with the best-fit juliet model shown
with the black line (see Sect. 4 for details on the modeling). The transits of
TOI-815b and TOI-815c are indicated with purple and orange, respectively.

##  Appendix B RVs with ESPRESSO

Table 4: ESPRESSO RV observations for TOI-815, along with activity indicator
measurements, i.e. FWHM, bisector, and logR${}^{\prime}_{\mathrm{HK}}$.

Time |  RV |  $\sigma_{\mathrm{RV}}$ |  FWHM(2) |  BIS(2) |  logR${}^{\prime}_{\mathrm{HK}}$(3) |  texp
---|---|---|---|---|---|---
[BJD-2 450 000] |  [m s-1 ] |  [m s-1 ] |  [m s-1 ] |  |  |  [s]
9308.60433 |  -2973.3950 |  0.6145 |  6474.143 |  -33.038 |  -4.7084 |  600
9310.60679 |  -2970.1704 |  0.7233 |  6470.941 |  -33.052 |  -4.7338 |  600
9313.73506(1) |  -2962.6373 |  0.5712 |  6492.966 |  -30.291 |  -4.6757 |  600
9316.52823 |  -2967.0427 |  0.7459 |  6494.718 |  -23.625 |  -4.7107 |  600
9318.53033 |  -2968.9726 |  0.5863 |  6478.690 |  -29.060 |  -4.6762 |  600
9324.54484 |  -2979.7547 |  0.8117 |  6513.799 |  -28.181 |  -4.6859 |  600
9327.50777(1) |  -2965.7154 |  0.4397 |  6482.962 |  -32.938 |  -4.6229 |  600
9329.52592(1) |  -2971.0105 |  0.9842 |  6543.347 |  -29.283 |  -4.7529 |  600
9332.57219 |  -2979.7266 |  0.6474 |  6491.382 |  -26.738 |  -4.6824 |  600
9336.49956(1) |  -2982.4246 |  1.3115 |  6515.936 |  -40.000 |  -4.9931 |  600
9339.48187 |  -2973.6498 |  0.5988 |  6468.715 |  -33.267 |  -4.7189 |  600
9341.51194(1) |  -2966.5835 |  0.9108 |  6513.090 |  -35.131 |  -4.7516 |  600
9346.54064 |  -2972.2321 |  0.5959 |  6482.420 |  -24.888 |  -4.6930 |  600
9373.49668 |  -2973.5189 |  0.5467 |  6462.905 |  -33.991 |  -4.6367 |  600
9375.51475 |  -2971.6794 |  0.6513 |  6471.909 |  -32.600 |  -4.6265 |  600
9406.50596 |  -2975.2108 |  0.7212 |  6492.234 |  -32.504 |  -4.6149 |  600
9410.48054 |  -2971.8872 |  0.6759 |  6499.208 |  -30.202 |  -4.6978 |  600
9415.50797 |  -2971.5355 |  0.8107 |  6491.252 |  -26.700 |  -4.6620 |  600
9673.54348 |  -2977.4158 |  0.5684 |  6488.390 |  -34.257 |  -4.6894 |  700
9675.51515 |  -2973.7947 |  0.9406 |  6487.917 |  -25.694 |  -4.6418 |  700
9677.65231 |  -2972.6960 |  0.8151 |  6492.405 |  -31.643 |  -4.6580 |  700
9679.55323(1) |  -2973.4334 |  0.5445 |  6497.504 |  -32.766 |  -4.6080 |  700
9683.53185 |  -2981.6090 |  0.4694 |  6497.328 |  -29.483 |  -4.6329 |  700
9685.55122 |  -2980.4065 |  0.7458 |  6496.394 |  -24.517 |  -4.6203 |  700
9687.52944 |  -2972.0070 |  0.4487 |  6495.493 |  -29.672 |  -4.6432 |  700
9695.52289 |  -2967.3200 |  0.9486 |  6512.740 |  -29.566 |  -4.7299 |  700
9697.58163 |  -2965.9727 |  0.5759 |  6511.966 |  -26.858 |  -4.6611 |  700
9784.49262 |  -2977.5647 |  0.6128 |  6525.410 |  -22.486 |  -4.6248 |  700
9904.83330 |  -2969.5919 |  0.4342 |  6487.460 |  -35.249 |  -4.6440 |  700
9905.77809 |  -2970.1154 |  0.6533 |  6491.720 |  -34.630 |  -4.6797 |  700
9906.80910 |  -2970.1976 |  0.4551 |  6492.162 |  -33.028 |  -4.6399 |  700
9907.79772 |  -2969.7455 |  1.3028 |  6534.337 |  -38.113 |  -4.8075 |  700
9910.75047 |  -2966.0686 |  0.7246 |  6503.361 |  -25.248 |  -4.6902 |  700
9917.72905 |  -2978.0604 |  0.6064 |  6481.170 |  -33.996 |  -4.6673 |  700

  * $\bullet$

Notes:

  * a

(1) These data did not pass the ESPRESSO flux quality control and therefore
were excluded from the analysis.

  * a

(2) The uncertainties of FWHM and BIS are twice the uncertainties of the RVs
(2$\sigma_{\mathrm{RV}}$).

  * a

(3) The median uncertainty of logR${}^{\prime}_{\mathrm{HK}}$(3) is 0.03.

##  Appendix C Spectral synthesis

Figure 16: Spectral synthesis around the Li doublet. The vertical lines mark
the position of the lines.

##  Appendix D Generalized Lomb-Scargle periodograms

Figure 17: GLS periodogram of the TOI-815 spectroscopic and photometric time
series. More specifically: (1) RV ESPRESSO measurements, (2) RV residuals,
after the subtraction of TOI-815c signal, (2) RV residuals, after the
subtraction of TOI-815b and c signals, (3)-(5) FWHM, bisector, and Ca II H & K
line activity indicator (logR${}^{\prime}_{\mathrm{HK}}$), (7) TESS Sector 63
photometric data, and (8)-(10) WASP photometric data (2008, 2012, and 2023).
The two black vertical lines correspond to the minimum and maximum rotation
periods (Prot,min $\sim$ 14.2 days and Prot,max $\sim$ 16.5 days). The purple
and orange vertical lines correspond to the orbital periods of planets b (Pb
$\sim$ 11.2 days) and c (Pc $\sim$ 35 days), respectively. The gray horizontal
lines indicating the 1$\%$ FAP.

##  Appendix E Priors

Table 5: Maximum values and 68$\%$ confidence intervals of the posterior
distributions from joint modeling with juliet.

Parameter |  Prior(a) |  Value(b)
---|---|---
Jump stellar parameters |  |
Stellar density, $\rho_{*}$ ($\it{\rho_{\mathrm{\odot}}}$) |  $\mathcal{N}(2396,139)$ |  $2431^{+122}_{-125}$
Jump parameters for TOI-815b |  |
Orbital period, $\it{P_{orb}}$ (days) |  $\mathcal{U}(10,12)$ |  $11.197259^{+0.000018}_{-0.000017}$
Transit epoch, $T_{0}$ (BJD) |  $\mathcal{N}(2459287.6,0.1)$ |  $2459287.6028^{+0.0009}_{-0.0009}$
Scaled planetary radius, $R_{P}$/$R_{*}$ |  $\mathcal{U}(0,1)$ |  $0.0349^{+0.0004}_{-0.0004}$
Impact parameter, b |  $\mathcal{U}(0,1)$ |  $0.279^{+0.058}_{-0.073}$
Eccentricity, $e$ |  $\mathcal{F}(0)$ |  0 (adopted, 3$\sigma$ $<$ 0.20)
Argument of periastron, $\omega_{*}$ (deg) |  $\mathcal{F}(90)$ |  -
RV semi-amplitude, K (m s-1) |  $\mathcal{U}(0,100)$ |  $2.6^{+0.5}_{-0.5}$
Jump parameters for TOI-815c |  |
Orbital period, $\it{P_{orb}}$ (days) |  $\mathcal{U}(34,36)$ |  $34.976145^{+0.000099}_{-0.000097}$
Transit epoch, $T_{0}$ (BJD) |  $\mathcal{N}(2459288.6,0.1)$ |  $2459288.6331^{+0.0023}_{-0.0024}$
Scaled planetary radius, $R_{P}$/$R_{*}$ |  $\mathcal{U}(0,1)$ |  $0.0312^{+0.0011}_{-0.0010}$
Impact parameter, b |  $\mathcal{U}(0,1)$ |  $0.839^{+0.015}_{-0.016}$
Eccentricity, $e$ |  $\mathcal{F}(0)$ |  0 (adopted, 3$\sigma$ $<$ 0.22)
Argument of periastron, $\omega_{*}$ (deg) |  $\mathcal{F}(90)$ |  -
RV semi-amplitude, K (m s-1) |  $\mathcal{U}(0,100)$ |  $5.5^{+0.5}_{-0.5}$
Instrumental photometric parameters |  |
Offset relative flux, $M_{\mathrm{TESS}_{S9}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  $-8.6^{+8.0}_{-8.2}$
Jitter, $\sigma_{w,\mathrm{TESS}_{S9}}$ (ppm) |  log$\mathcal{U}(0.01,1000)$ |  $0.5^{+6.5}_{-0.5}$
Offset relative flux, $M_{\mathrm{TESS}_{S36}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  $-72^{+411}_{-417}$
Jitter, $\sigma_{w,\mathrm{TESS}_{S36}}$ (ppm) |  log$\mathcal{U}(0.01,1000)$ |  $466^{+8}_{-8}$
Offset relative flux, $M_{\mathrm{TESS}_{S63}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  $-3594^{+3821}_{-6137}$
Jitter, $\sigma_{w,\mathrm{TESS}_{S63}}$ (ppm) |  log$\mathcal{U}(0.01,1000)$ |  $457^{+11}_{-11}$
Offset relative flux, $M_{\mathrm{CHEOPS}_{1}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  $-197^{+46}_{-44}$
Jitter, $\sigma_{w,\mathrm{CHEOPS}_{1}}$ (ppm) |  log$\mathcal{U}(0.01,1000)$ |  $299^{+23}_{-22}$
Offset relative flux, $M_{\mathrm{CHEOPS}_{2}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  $-474^{+48}_{-48}$
Jitter, $\sigma_{w,\mathrm{CHEOPS}_{2}}$ (ppm) |  log$\mathcal{U}(0.01,1000)$ |  $404^{+24}_{-23}$
Offset relative flux, $M_{\mathrm{CHEOPS}_{3}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  $-322^{+34}_{-35}$
Jitter, $\sigma_{w,\mathrm{CHEOPS}_{3}}$ (ppm) |  log$\mathcal{U}(0.01,1000)$ |  $293^{+19}_{-19}$
Offset relative flux, $M_{\mathrm{LCO}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  1089${}^{+88}_{-88}$
Jitter, $\sigma_{w,\mathrm{LCO}}$ (ppm) |  log$\mathcal{U}(0.01,5000)$ |  986${}^{+72}_{-68}$
Offset relative flux, $M_{\mathrm{ASTEP}}$ (10-6) |  $\mathcal{N}(0,0.1)$ |  761${}^{+133}_{-131}$
Jitter, $\sigma_{w,\mathrm{ASTEP}}$ (ppm) |  log$\mathcal{U}(0.01,5000)$ |  3369${}^{+89}_{-85}$
GP parameters |  |
$\sigma_{GP,\mathrm{TESS}_{S9}}$ (relative flux) |  $\mathcal{F}(0)$ |  -
$\rho_{GP,\mathrm{TESS}_{S9}}$ (days) |  $\mathcal{F}(0)$ |  -
$\sigma_{GP,\mathrm{TESS}_{S36}}$ (10-4 relative flux) |  log$\mathcal{U}$(10-6, 100) |  $12^{+2}_{-2}$
$\rho_{GP,\mathrm{TESS}_{S36}}$ (days) |  log$\mathcal{U}(0.001,100)$ |  $0.94^{+0.13}_{-0.11}$
$\sigma_{GP,\mathrm{TESS}_{S63}}$ (10-4 relative flux) |  log$\mathcal{U}$(10-6, 100) |  $57^{+44}_{22}$
$\rho_{GP,\mathrm{TESS}_{S63}}$ (days) |  log$\mathcal{U}(0.001,100)$ |  $9.11^{+4.37}_{-2.77}$
Instrumental RV parameters |  |
Jitter, $\sigma_{w,\mathrm{ESPRESSO}}$ (m s-1) |  log$\mathcal{U}(0.001,100)$ |  $1.88^{+0.29}_{-0.25}$
Systemic RV, $\mu_{\mathrm{ESPRESSO}}$ (m s-1) |  $\mathcal{U}(-3500,-1000)$ |  $-2973.61^{+0.36}_{-0.36}$
Limb darkening parameters |  |
$q_{1,\mathrm{ESPRESSO}}$ |  $\mathcal{N}(0.408,0.017)$ |  0.407${}^{+0.015}_{-0.015}$
$q_{2,\mathrm{ESPRESSO}}$ |  $\mathcal{N}(0.349,0.022)$ |  0.345${}^{+0.020}_{-0.019}$
$q_{1,\mathrm{CHEOPS}}$ |  $\mathcal{N}(0.519,0.012)$ |  0.518${}^{+0.011}_{-0.011}$
$q_{2,\mathrm{CHEOPS}}$ |  $\mathcal{N}(0.408,0.012)$ |  0.406${}^{+0.011}_{-0.011}$
$q_{1,\mathrm{LCO}}$ |  $\mathcal{N}(0.324,0.008)$ |  0.324${}^{+0.007}_{-0.007}$
$q_{2,\mathrm{LCO}}$ |  $\mathcal{N}(0.314,0.014)$ |  0.314${}^{+0.012}_{-0.012}$
$q_{1,\mathrm{ASTEP}}$ |  $\mathcal{N}(0.492,0.016)$ |  0.492${}^{+0.014}_{-0.014}$
$q_{2,\mathrm{ASTEP}}$ |  $\mathcal{N}(0.376,0.017)$ |  0.378${}^{+0.015}_{-0.016}$

  * $\bullet$

Notes: (a)For the priors, $\mathcal{N}$($\mu$, $\sigma^{2}$) indicates a
normal distribution with mean $\mu$ and variance $\sigma^{2}$,
$\mathcal{U}$(a, b) a uniform distribution between a and b,
log$\mathcal{U}$(a, b) a log-uniform distribution between a and b and
$\mathcal{F}$(a) a parameter fixed to value a. (b)The posterior estimate
indicates the median value and then error bars for the 68 $\%$ credibility
intervals.

##  Appendix F CHEOPS alias hunting

Table 6: CHEOPS photometric observations of TOI-815.

Visit | Expected planet | Start time | Duration | Eff. | RMS | Alias | Detection | File key |  Detrending(d)
---|---|---|---|---|---|---|---|---|---
|  | [UTC] | [h] |  [$\%$] | ppm | [d] |  |  |
1 | c | 2022-03-06 00:36 | 8.34 | 73 | 422 | 22.26 | no | CH_PR110048_TG018001_V0300 | ${(e)}$
2 | c | 2022-03-06 16:26 | 8.37 | 60 | 461 | 20.99 | no | CH_PR110048_TG018901_V0300 | ${(e)}$
3 | c | 2022-03-17 04:35 | 9.84 | 61 | 417 | 22.95 | no | CH_PR110048_TG018101_V0300 | ${(e)}$
4(a) | b | 2022-03-18 08:36 | 11.64 | 66 | 514 | 11.20 | yes | CH_PR100031_TG050801_V0300 |  t, x, y, y2, sin($\phi$), cos($\phi$), cos(2$\phi$)
5 | c | 2022-04-13 01:20 | 8.47 | 60 | 429 | 17.92 | no | CH_PR110048_TG025201_V0300 | ${(e)}$
6 | c | 2022-04-30 06:47 | 8.34 | 51 | 497 | 29.38 | no | CH_PR110048_TG027601_V0300 | t (clear flux trend in data)
7 | c | 2023-02-15 11:09 | 8.29 | 69 | 520 | 31.93 | no | CH_PR110048_TG036501_V0300 | ${(e)}$
8 | c | 2023-02-22 01:40 | 9.04 | 68 | 447 | 25.33 | no | CH_PR110048_TG036701_V0300 | ${(e)}$
9 | c | 2023-02-27 14:48 | 8.22 | 64 | 580 | 19.85 | no | CH_PR110048_TG036901_V0300 | ${(e)}$
10(b) | b | 2023-02-28 14:16 | 8.39 | 62 | 460 | 18.83 | yes | CH_PR110048_TG037001_V0300 |  t2, x, y, y2, bg, sin($\phi$), sin(2$\phi$), cos(2$\phi$)
11(c) | c | 2023-04-23 08:31 | 10.09 | 58 | 462 | 34.98 | yes | CH_PR110048_TG042201_V0300 |  t, t2, y, y2, bg, cos($\phi$)

  * $\bullet$

Notes:

  * a

(a) Figure 7, left.

  * b

(b) Figure 7, center.

  * c

(c) Figure 7, right.

  * d

(d) The correlation terms were determined using pycheops.

  * e

(e) Detrended with: x, y, bg, cos($\phi$), cos(2$\phi$), cos(3$\phi$),
$\Delta$T, sin($\phi$), sin(2$\phi$), sin(3$\phi$), bg2, x2, y2

Figure 18: CHEOPS unsuccessful transit observations that ruled out the most
likely period aliases of TOI-815c. The data have been detrended using the
detrending terms presented in Table 6.


# Wide twin binaries are extremely eccentric: evidence of twin binary
formation in circumbinary disks

Hsiang-Chih Hwang School of Natural Sciences, Institute for Advanced Study,
Princeton, 1 Einstein Drive, NJ 08540, USA  Kareem El-Badry Center for
Astrophysics $\mid$ Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138,
USA  Hans-Walter Rix Max-Planck-Institut für Astronomie, Königstuhl 17,
D-69117 Heidelberg, Germany  Chris Hamilton School of Natural Sciences,
Institute for Advanced Study, Princeton, 1 Einstein Drive, NJ 08540, USA
Yuan-Sen Ting Research School of Astronomy & Astrophysics, Australian National
University, Cotter Rd., Weston, ACT 2611, Australia  School of Computing,
Australian National University, Acton ACT 2601, Australia  Nadia L. Zakamska
School of Natural Sciences, Institute for Advanced Study, Princeton, 1
Einstein Drive, NJ 08540, USA  Department of Physics & Astronomy, Johns
Hopkins University, Baltimore, MD 21218, USA

###### Abstract

The Gaia mission recently revealed an excess population of equal-mass “twin”
wide binaries, with mass ratio $q\gtrsim 0.95$, extending to separations of at
least 1000 AU. The origin of this population is an enigma: twin binaries are
thought to form via correlated accretion in circumbinary disks, but the
largest observed protostellar disks have radii of $\sim 100$ AU, far smaller
than the separations of the widest twins. Here, we infer the eccentricity
distribution of wide twins from the distribution of their $v$-$r$ angles,
i.e., the angle between the components’ separation and relative velocity
vectors. We find that wide twins must be on extremely eccentric orbits. For
the excess-twin population at 400-1000 AU, we infer a near-delta function
excess of high-eccentricity system, with eccentricity $0.95\lesssim e\leq 1$.
These high eccentricities for wide twins imply pericenter distances of order
$10$ AU and suggest that their orbits were scattered via dynamical
interactions in their birth environments, consistent with a scenario in which
twins are born in circumbinary disks and subsequently widened. These results
further establish twin wide binaries as a distinct population and imply that
wide twins can be used as a probe of the dynamical history of stellar
populations.

binaries: general – stars: kinematics and dynamics – stars: formation –
protoplanetary disks

##  1 Introduction

Binary population demographics encode information about the star formation
process and subsequent dynamical processing. One striking feature of the
binary population is the existence of an excess population of equal-mass
“twin” binaries with mass ratios $0.95\lesssim q\leq 1$, where $q=M_{2}/M_{1}$
and $M_{1}$ and $M_{2}$ are the masses of the primary and secondary. This twin
excess is strongest at short periods, and was long thought to exist only at
orbital periods $P_{\rm orb}\lesssim 40$ days (e.g. Lucy & Ricco, 1979;
Tokovinin, 2000). However, studies of large samples of binaries have recently
shown that an excess population exists even at very wide separations,
extending beyond 1000 AU for solar-type stars (e.g. Soderhjelm, 2007; Moe & Di
Stefano, 2017; El-Badry et al., 2019).

The leading hypothesis for the origin of the twin excess is that twins form in
circumbinary disks. In this scenario, the lower-mass star has a higher
accretion rate from the circumbinary disk due to its wider orbit around the
common center of mass, or due to the dynamics of the accretion streams (Farris
et al., 2014; Duffell et al., 2020), thus driving the mass ratio to unity.
Torques from the disk may also shrink the binary orbit, potentially explaining
why the amplitude of the twin excess is largest at close separations
(Tokovinin & Moe, 2020). However, this picture is still under debate. First,
it is unclear whether the lower-mass star in a circumbinary disk will always
have the higher accretion rate, with different codes making qualitatively
different predictions (Bate & Bonnell, 1997; Ochi et al., 2005; Young &
Clarke, 2015). Furthermore, how the interaction between the binary and the
circumbinary disk affects a binary’s orbital evolution remains an open
question (Artymowicz et al., 1991; Artymowicz & Lubow, 1994; Pichardo et al.,
2005; Shi et al., 2012; Miranda et al., 2017; Ragusa et al., 2020; Tiede et
al., 2020; Heath & Nixon, 2020; Dittmann & Ryan, 2022). Thus, twin binaries
provide a unique observational opportunity to study the effects of accretion
from circumbinary disks.

The excess of twin binaries with separations $\gtrsim 100$ AU is particularly
puzzling (Soderhjelm, 2007). Their existence out to $>1000$ AU is now
confirmed by Gaia and is not due to selection effects (El-Badry et al., 2019).
Observed protoplanetary and circumbinary disks have typical radii of $\sim$100
AU (e.g. Andrews et al., 2018), and so wide twins cannot have formed in
circumbinary disks at their current separations. Although the excess twin
population only contains a few percent of all binaries at these separations,
the twin excess is manifest as a sharp, step function-like jump in the mass
ratio distribution above $q_{\rm twin}\approx 0.95$ (El-Badry et al., 2019).
Because $q_{\rm twin}$ does not vary between close and wide binaries, it is
natural to assume that wide twins form via the same process as close twins.
One possible formation scenario is that wide twins formed at closer
separations ($10-100$ AU) within a circumbinary disk and were subsequently
widened by dynamical interactions. However, the nature of this interaction and
widening process is poorly understood.

If twins’ formation involves dynamical widening either through gravitational
interactions with other stars, or through recoils or kicks to one of the
components, then highly eccentric orbits are expected for wide twins.
Furthermore, binary eccentricities are a direct observable from simulations of
binary-circumbinary disk interaction (e.g. Artymowicz & Lubow, 1994; Roedig et
al., 2011; Dittmann & Ryan, 2022). Therefore, the eccentricity of wide twin
binaries is a critical connection between observational and theoretical work.

In this letter, we investigate the eccentricities of wide twin binaries. The
individual eccentricities of wide binaries are challenging to measure due to
their long orbital periods ($\gtrsim 10^{3}$ yr for binaries at $>10^{2}$ AU),
but the population eccentricity distribution can be statistically constrained
by the distribution of $v$-$r$ angles, the angle between the separation vector
($r$) and the relative velocity vector ($v$) of a wide binary (Hwang et al.,
2022). This paper is structured as follows. Sec. 2 explains the sample
selection. Sec. 3 presents the main results, showing that wide twin binaries
are highly eccentric. We discuss the results and conclude in Sec. 4.

##  2 Sample selection

Figure 1: The distributions of magnitude difference $\Delta G$ among the
components of wide binaries at different separations. The vertical dashed line
marks $\Delta G=0.25$ mag, the selection adopted for _twin_ wide binaries in
this paper. The distributions are normalized to unity at $\Delta G=1.5$-2 mag.
Compared to binaries at $>10000$ AU, binaries at separations $\lesssim 3000$
AU show a significant excess of twins at $\Delta G<0.25$ mag.

The high-precision astrometry from Gaia (Gaia Collaboration et al., 2016) has
enabled large-scale wide binary search (Oh et al., 2017; El-Badry & Rix, 2018;
Tian et al., 2020; Hartman & Lépine, 2020; Hwang et al., 2021). In this work,
we use the $\sim 1$-million wide binaries within 1 kpc (El-Badry et al., 2021)
selected from Gaia early Data Release 3 (eDR3, Gaia Collaboration et al.,
2021). To avoid contamination, we require the chance-alignment probability
$<0.1$ (El-Badry et al., 2021).

We use the $v$-$r$ angle method to infer the eccentricity distribution of wide
binaries (Tokovinin, 1998; Tokovinin & Kiyaeva, 2016; Tokovinin, 2020; Hwang
et al., 2022). Qualitatively, the $v$-$r$ angle distribution of randomly-
oriented circular binaries peaks at 90∘, and the peak moves toward 0∘ and 180∘
for more eccentric orbits (e.g. Fig. 3 in Hwang et al. 2022). The observed
$v$-$r$ angles of Gaia wide binaries are computed from the angle between
projected $v$ and projected $r$ vectors, where $r$ is the vector connecting
two member stars’ sky coordinates and $v$ is the vector of two stars’ proper
motion differences. We require that all wide binaries in the sample have
angular separations $>1.5$″ to avoid Gaia’s systematics on $v$-$r$ angles in
pairs below $1.25$″ (Hwang et al., 2022). To reliably measure $v$-$r$ angles,
we further limit the sample to have parallaxes $>5$ mas (i.e. distances $<$200
pc) and proper motion differences that are $3\sigma$ from being zero. With
Gaia’s proper motion precision of $\sim 0.1$ mas yr-1, simulations in Hwang et
al. (2022) show that these selection criteria can recover 89.7% of 1000 AU
wide binaries’ proper motion differences at $>3\sigma$ at distances of 200 pc,
assuming that they are equal-solar-mass binaries with the so-called thermal
eccentricity distribution ($f(e)de=2ede$). All $v$-$r$ angle measurements and
related quantities used in this work are publicly available from Hwang et al.
(2022).

With these selection criteria, the classification based on the absolute
$G$-band magnitudes and $BP$-$RP$ colors (binary_type in El-Badry et al. 2021)
shows that 88% of the wide binaries with separations of 3000-10000 AU (hence
angular separations $>15$″) are double main-sequence (MS) binaries, 10% are
white dwarf (WD)-MS, and the rest are double-WD wide binaries. In Gaia eDR3,
$G$-band photometry has an angular resolution of $\sim 0.7$″ (Fabricius et
al., 2021). However, because $BP$ and $RP$ photometry does not have deblending
treatment, $BP$ and $RP$ are not reliable for pairs $<2$″. For this reason, we
do not use $BP$-$RP$ colors to differentiate MS from WD in this work, but we
expect the sample to be dominated ($\gtrsim 90$%) by MS-MS wide binaries.

Following El-Badry et al. (2019), we select wide twin binaries by $\Delta
G<0.25$ mag, where $\Delta G$ is the difference in Gaia’s broadband photometry
$G$ of two component stars of a wide binary. For MS-MS binaries, $\Delta
G<0.25$ mag corresponds to a mass-ratio selection of $q>0.95$, with some
slight dependence on the primary mass (El-Badry et al., 2019). Fig. 1 shows
the distributions of $\Delta G$ for different binary separations. Compared to
binaries at $>10000$ AU, wide binaries $\lesssim 3000$ AU have an excess of
twin binaries at $\Delta G<0.25$ mag. The median error of $\Delta G$ in our
sample is 0.001 mag, much smaller than the width of the excess twin component
at $\Delta G\sim 0.25$ mag.

Binaries selected by $\Delta G<0.25$ mag contain two components: one component
is from the extension of the smooth, power-law-like mass-ratio distribution at
$q\gtrsim 0.5$, and the other component is from the excess twins that only
contribute to $q>q_{twin}\sim 0.95$ (Moe & Di Stefano, 2017; El-Badry et al.,
2019). Here we use binaries at separations $>10000$ AU as the baseline where
no significant excess twins are present (El-Badry et al., 2019) to quantify
the fraction of excess twins at other separations. Specifically, binaries with
$\Delta G<0.25$ mag constitute $p_{400-1000}=7.70\pm 0.16$% of all wide
binaries at 400-1000 AU, compared to $p_{>10000}=5.89\pm 0.36$% at $>10000$
AU. Therefore, at $\Delta G<0.25$ mag and binary separations of 400-1000 AU,
the fraction of excess twins is
$\mathcal{E}_{400-1000}=(p_{400-1000}-p_{>10000})/p_{400-1000}=23.5\pm 5.1$%.
Similarly, $\mathcal{E}_{1000-3000}=4.0\pm 4.8$% for separations at 1000-3000
AU, where the excess is more significant for a certain primary mass range (El-
Badry et al., 2019). Due to the difference in the detailed sample selection
(e.g. using $BP$ and $RP$ photometry or not, different distance distribution,
and lack of corrections for incompleteness), the derived twin excess
$\mathcal{E}$ may differ from the intrinsic mass-ratio distribution models
inferred by El-Badry et al. (2019).

Figure 2: The $v$-$r$ angle distributions for wide binaries with separations
at 400-1000 AU. The solid and dashed black lines are simulated binaries with
arbitrary normalization, showing that circular orbits (solid) have a $v$-$r$
angle distribution peaking at $90^{\circ}$, and the peaks for eccentric
$e=0.9$ orbits (dashed) are closer to $0^{\circ}$ and $180^{\circ}$. The
$v$-$r$ angle distributions of wide twin binaries (blue) are strongly enhanced
at 0∘ and 180∘ compared to non-twins (orange), indicating the presence of
highly eccentric twins with $e>0.9$.  Figure 3: The power-law indices $\alpha$
of eccentricity distributions ($f(e)\propto e^{\alpha}$) as a function of
binary separations. The black symbols are the results from Hwang et al. (2022)
for all binaries, which are dominated by non-twins; the red markers show
$\alpha$ and separation bin sizes (horizontal) for the twin binaries. At
100-3000 AU, twin binaries are significantly more eccentric (larger $\alpha$)
than other wide binaries at the same separations.

Figure 4: Left: the $v$-$r$ angle distributions for wide binaries at 400-1000
AU (blue) and for simulated models (black and red). The best-fit power-law
model of $\alpha=2.14$ (dotted black) well describes the overall trend of twin
binaries (solid blue) but with some noticeable differences. The dashed red
line corresponds to 81.1% of normal wide binaries and 18.9% of
$e=0.948$-$0.992$ binaries, which agrees with the observation better than the
power-law model. Right: the modeled eccentricity distribution (the dashed red
line in the left panel) for twin binaries at separations 400-1000 AU. The
black shaded region shows the 20-80 percentiles of the model.

##  3 Wide twin binaries are eccentric

Fig. 2 shows the distribution of $v$-$r$ angles for wide binaries with binary
separations of 400-1000 AU. Wide twin binaries (blue) are selected by $\Delta
G<0.25$ mag, and wide non-twin binaries (orange) by $\Delta G>0.25$ mag. Fig.
2 also shows that the simulated $v$-$r$ angle distributions for circular
binaries (solid black) with random orientation strongly peaks at $90^{\circ}$,
and the peaks move toward $0^{\circ}$ and $180^{\circ}$ for more eccentric
orbits (e.g. $e=0.9$ for the dashed black line). The $v$-$r$ angles of wide
twin binaries are strongly enhanced at 0∘ and 180∘, indicating the presence of
highly eccentric binaries with $e>0.9$. In contrast, wide non-twin binaries in
this separation range have a flat $v$-$r$ angle distribution, suggesting an
underlying thermal eccentricity distribution (Hwang et al., 2022).

Using the Bayesian inference detailed in Hwang et al. (2022), we measure the
eccentricity distributions of wide twin binaries as a function of projected
binary separations (denoted by $s$). Specifically, we model the eccentricity
distribution $f(e)$ as a power law $f(e)=(1+\alpha)e^{\alpha}$, and obtain the
best fit $\alpha$ given the observed $v$-$r$ angle distribution.

Fig. 3 shows the best-fit power-law indices $\alpha$ for the eccentricity
distributions of twin binaries (red) as a function of binary separations. For
comparison, the black points and solid black line show the results from Hwang
et al. (2022) for all wide binaries, which are dominated by non-twin MS-MS
binaries. If one excludes twin binaries in the black points from Hwang et al.
(2022), then their $\alpha$ values are only modified by an amount smaller than
their measurement uncertainties. The red dashed horizontal error bars indicate
the bin sizes, and the vertical error bars represent the 68% credible
interval. At separations of 100-3000 AU, twin binaries have significantly
larger $\alpha$ than wide non-twin binaries with similar separations, meaning
that they are highly eccentric. There is no significant difference in $\alpha$
between twins and non-twins at $>3000$ AU, in agreement with the fact that
there is no significant twin excess at larger separations (El-Badry et al.,
2019). Interestingly, at 50-100 AU, twin binaries do not have significantly
different $\alpha$ even though the excess twin population is very significant
at these small separations (El-Badry et al., 2019), hinting that the majority
of twin binaries below 100 AU do not experience the orbit-widening processes
that alter their orbital eccentricities and they may have formed in disks at
their current separations. We remind the reader that “twin binaries” in Fig. 3
includes both “excess” twins and the smooth background population, implying
that the excess twin are even more eccentric than suggested by Fig. 3.

The left panel in Fig. 4 shows the $v$-$r$ angle distribution of wide twin
binaries at 400-1000 AU and the simulated binaries with the best-fit
$\alpha=2.14$ (dotted black line). Although the overall observed distribution
agrees with the best-fit power-law result, there are some subtle but
significant differences. Specifically, compared to the observed distribution
(blue), the power-law model (black) is lower at 90∘, 0∘, and 180∘, and is
higher at 45∘ and 135∘. These differences suggest that the eccentricity
distribution of wide twin binaries is not a perfect power law.

Alternatively, we can model the eccentricity distribution of wide twin
binaries as a sum of two populations: “normal” binaries that just happen to
have $q\approx 1$, and excess twins. We assume that a fraction $F_{\rm ecc}$
of twins have unusually eccentric orbits and the remaining $1-F_{\rm ecc}$ of
them follow the eccentricity distribution of non-twin binaries at similar
separations. Therefore, the total eccentricity distribution $f(e)$ of wide
twin binaries is

$f(e)=(1-F_{\rm ecc})f_{\rm normal}(e)+F_{\rm ecc}f_{\rm ecc}(e),$ (1)

where $f_{\rm normal}(e)$ and $f_{\rm ecc}(e)$ are the eccentricity
distributions for normal binaries and the excess eccentric binaries,
respectively. At 400-1000 AU, the non-twin ($\Delta G>0.25$ mag) wide binaries
have a flat $v$-$r$ angle distribution and its best-fit $\alpha$ is $1.05\pm
0.05$, and therefore we adopt a thermal eccentricity distribution ($\alpha=1$)
for $f_{\rm normal}(e)=2e$ (Hwang et al., 2022). For excess eccentric
binaries, we choose $f_{\rm ecc}(e)$ as a top-hat function whose value is
$1/(e_{1}-e_{0})$ between two free parameters $e_{0}$ and $e_{1}$, and zero
elsewhere. With details shown in Appendix A, we use the affine invariant
Markov chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al.,
2013) to constrain $F_{\rm ecc}$, $e_{0}$, and $e_{1}$.

The best-fit parameters are $F_{\rm ecc}=0.189^{+0.042}_{-0.033}$,
$e_{0}=0.948^{+0.021}_{-0.017}$, $e_{1}=0.992^{+0.004}_{-0.014}$. The best-fit
values are the most probable values of the marginalized posterior
distributions, and the uncertainties represent the highest posterior density
interval that includes 68 per cent of the area. Therefore, the observed
$v$-$r$ angle distribution is best fit by $F_{\rm ecc}=18.9^{+4.2}_{-3.3}$% of
highly eccentric ($e>e_{0}=0.948$) binaries among the twin binaries. The red
dashed histogram in Fig. 4 shows the simulated $v$-$r$ angle distribution of
this model, well consistent with the observed distribution. The best fit of
$f(e)$ (Eq. 1) in the right panel of Fig. 4 illustrates the presence of the
highly eccentric twin binaries.

The median error of the $v$-$r$ angles for 400-1000 AU twin binaries is 3.3∘,
meaning that the highest measurable eccentricity is
$e=\cos(3.3^{\circ})=0.9983$ (Hwang et al., 2022). Therefore, our measured
$e_{1}=0.992^{+0.004}_{-0.014}$ is approaching the precision limit. Despite
the high eccentricities close to the hyperbolic regime ($e>1$), the symmetric
observed $v$-$r$ angle distribution in Fig. 4 indicates that they are still on
stable Keplerian orbits; otherwise, disrupting binaries on hyperbolic orbits
would have $v$-$r$ angles enhanced at $0^{\circ}$ but not at 180∘.

One natural explanation is that all excess twins identified from the mass-
ratio distribution are highly eccentric. In this case, we would expect the
twin excess fraction from the mass-ratio distribution equals the fraction of
the eccentric component, i.e. $\mathcal{E}=F_{\rm ecc}$. At separations
400-1000 AU, the fraction of eccentric twins $F_{\rm ecc}=18.9^{+4.2}_{-3.3}$%
is consistent with the excess twin fraction at $\Delta G<0.25$ mag measured
from the $\Delta G$ (i.e. mass ratio) distribution,
$\mathcal{E}_{400-1000}=23.5\pm 5.1$% (Sec. 2). However, with the similar
analysis applying to twin binaries at 1000-3000 AU, the fraction of eccentric
twins is $F_{ecc,1000-3000}=21.6^{+5.9}_{-4.5}$%, which is much higher than
$\mathcal{E}_{1000-3000}=4.0\pm 4.8$%. Therefore, the connection between the
excess in the mass ratio distribution ($\mathcal{E}$) and the excess in the
eccentricity distribution ($F_{\rm ecc}$) remains not fully established, and
future work is needed to investigate if $F_{\rm ecc}$ follows the same primary
mass dependence as $\mathcal{E}$ (El-Badry et al., 2019).

Fig. 5 shows the eccentricity distribution power-law indices $\alpha$ versus
$\Delta G$. For all separation bins between 100 and 3000 AU, $\alpha$ strongly
increases at $\Delta G<0.25$ mag, the criterion used to select twin binaries.
The marginal increase in $\alpha$ at $\Delta G=0.25-0.5$ for binaries at
400-1000 AU may hint that this bin still has some contribution from the
eccentric twin binaries. At $\Delta G>0.5$ mag, $\alpha$ becomes flat for all
separation bins. Therefore, the high-eccentricity population is specifically
for twin binaries ($\Delta G<0.25$ mag), and the eccentricity distribution
does not strongly depend on mass ratios anymore at $\Delta G>0.5$ mag
($q\lesssim 0.9$).

Figure 5: The power-law indices $\alpha$ of eccentricity distributions
($f(e)\propto e^{\alpha}$) as a function of $\Delta G$. Only $\Delta G<0.25$
mag has a significant increase in $\alpha$, suggesting that the high
eccentricity is specifically associated with the twin excess, instead of a
general trend with mass ratios.

##  4 Discussion and Conclusions

Wide twin binaries are a mysterious population that has near equal-mass ratios
($q>0.95$), common among close binaries, but has large binary separations of
$\sim 1000$ AU. In this letter, using the $v$-$r$ angle method to infer the
eccentricity distribution (Hwang et al., 2022), we show that these wide twin
binaries are highly eccentric (Fig. 2, 3). Our result suggests that
$18.9^{+4.2}_{-3.3}$% of wide twin binaries (selected by $\Delta G<0.25$ mag)
at 400-1000 AU have eccentricities between $e_{0}=0.948^{+0.021}_{-0.017}$ and
$e_{1}=0.992^{+0.004}_{-0.014}$(Fig. 4). Despite the high eccentricities,
their symmetric $v$-$r$ angle distribution with respect to $90^{\circ}$
suggests that they are on stable Keplerian orbits, instead of being dissolving
binaries. This high eccentricity is specific to excess twins, as opposed to
being a smooth function of the mass ratio (Fig. 5).

Binaries on eccentric orbits will most likely be observed near their
apocenters where they spend most of their time due to the lower orbital
velocities. For apocenters $r_{apo}=1000$ AU, the inferred eccentricities of
$e=0.95$ and $0.99$ correspond to pericenters of
$r_{peri}=r_{apo}(1-e)/(1+e)=26$ AU and 5 AU, respectively. The high
eccentricities of $\lesssim 1000$ AU wide twins must originate from their
formation, since interactions with passing stars and secular torquing by the
Galactic tide mainly affect wide binaries at $\gtrsim 10^{4}$ AU, and the
latter have been shown only make the initial eccentricity distribution closer
to thermal (Hamilton, 2022). Pericenters smaller than the typical $\sim 100$
AU sizes of circumbinary disks are consistent with the scenario in which twins
are born in circumbinary disks and subsequently widen through subsequent
dynamical interactions in the birth environments.

The nature of the dynamical interactions causing the wide and eccentric orbits
remains uncertain. One approach to excite high eccentricities is an
instantaneous velocity kick. To change a circular orbit with an orbital
velocity $v_{circ}$ to an eccentric orbit with $e>0.9$, the required velocity
kick $v_{kick}$ is comparable with the original circular velocity
($v_{kick}=v_{circ}(1-\sqrt{1-e})$), where $v_{circ}=19$ km s-1 at 5 AU.
Therefore, an instantaneous kick with $v_{kick}\sim 10$ km s-1 taking place in
close binaries can produce eccentric ($e>0.9$) wide binaries. However, the
source of such strong velocity kicks is currently not known in star formation
processes. Protostellar outflows and jets can reach velocities of several
$100$ km s-1, but they are launched from the disk at the stage when the disk
mass is only $\sim 10$% of the central star (Bally, 2016), and thus are
unlikely to impose a strong kick on the star.

Chaotic three-body interactions can lead to the formation of close binaries
with highly eccentric tertiary companions (Reipurth & Mikkola, 2012). If the
close binary was unresolved, the resulting triple system would look like a
wide binary. However, radial velocity variations among wide twin binaries
suggest that close companions are not more common in wide twins than in wide
non-twins (El-Badry et al., 2019). Therefore, wide twin binaries are unlikely
to form from chaotic three-body interaction.

For hierarchical triples, the resolved outer companions would have less
eccentric orbits required by dynamical stability (Shatsky, 2001; Tokovinin &
Kiyaeva, 2016; Hwang et al., 2022). Since the presence of unresolved
companions can cause non-zero $\Delta G$, the equal-mass selection by $\Delta
G<0.25$ mag preferentially excludes systems with unresolved companions, thus
allowing more eccentric outer companions. This scenario is possible because it
is not uncommon to have unresolved companions in (predominantly non-twin) wide
pairs (El-Badry et al., 2018; Hwang et al., 2020; Fezenko et al., 2022). To
test this effect, we simulate the photometry of wide binaries using MIST
(Dotter, 2016; Choi et al., 2016) and
brutus111https://github.com/joshspeagle/brutus (Speagle et al. in prep), where
50% of them are assumed to have unresolved companions and all component stars’
masses are drawn from the Kroupa initial mass function (Kroupa, 2001). We find
that $\Delta G>0.25$ mag is only $\sim 5$% more likely to have unresolved
companions than $\Delta G<0.25$ mag, thus this potential selection effect is
not able to explain $F_{\rm ecc}=18.9^{+4.2}_{-3.3}$% of high-eccentricity
binaries at $\Delta G<0.25$ mag.

Twin binaries may form through the enhanced accretion and interaction with
circumbinary disks (Tokovinin & Moe, 2020). Then the strong interaction
between binaries and the circumbinary disk may increase binary’s eccentricity,
even though the observed $e>0.95$ is unusually high for typical disk-binary
interaction in simulations (Artymowicz et al., 1991; Cuadra et al., 2009;
Roedig et al., 2011). Alternatively, the process that widens binary orbits and
causes high eccentricities may not be specific to twins. It is possible that
such a process takes place among all close binaries, producing both eccentric
wide twins and non-twins out to $\sim 1000$ AU. Then since twins are more
common in close binaries, the high-eccentricity component is more apparent in
wide twins than wide non-twins. Future investigations are needed to establish
the connection between close binaries, wide binaries, and their mass ratios
and eccentricities.

Our results suggest that wide twin binaries have eccentricities approaching
unity. Therefore, the process making these twin binaries wide and eccentric
should also disrupt some of binaries during the star formation. These
disrupted binaries may contribute to the low-mass runaway or walkaway stars in
star-forming regions (e.g. Schoettler et al. 2020), predicting a population of
equal-mass runaway or walkaway pairs with opposite directions. Furthermore,
some of these wide twin binaries may have pericenters comparable to or smaller
than the radii of giant stars, which may lead to collision at later stellar
evolution and the formation of blue stragglers (Kaib & Raymond, 2014),
although the twin fraction is lower in more massive stars where they can
evolve to giants within the Hubble time (Moe & Di Stefano, 2017; El-Badry et
al., 2019).

## Acknowledgements

HCH appreciates the discussions with Jim Stone, Roman Rafikov, and Scott
Tremaine. HCH acknowledges the support of the Infosys Membership at the
Institute for Advanced Study. HWR acknowledges support from the GIF grant
I-95-303.5-2018. This work was supported by a grant from the Simons Foundation
(816048, CH). Y.S.T. acknowledges financial support from the Australian
Research Council through DECRA Fellowship DE220101520.

##  Appendix A Markov-Chain Monte-Carlo results

Our eccentricity model has three free parameters: $F_{\rm ecc}$, $e_{0}$, and
$e_{1}$ in Eq. 1. According to Bayes’ theorem, their posterior distributions
are

$p(\\{F_{\rm ecc},e_{0},e_{1}\\}|\\{\gamma_{obs,i}\\})\propto\\\ p(F_{\rm
ecc})p(e_{0})p(e_{1})\Pi_{i}\int
p(\gamma_{obs,i}|\gamma_{true,i})p(\gamma_{true,i}|e_{i})p(e_{i}|\\{F_{\rm
ecc},e_{0},e_{1}\\})d\gamma_{true,i}de_{i},$ (A1)

where $i$ is the index of individual binaries, $\gamma_{obs,i}$ is the
observed $v$-$r$ angle for binary $i$, $p(\gamma_{obs,i}|\gamma_{true,i})$ is
the error distribution for $v$-$r$ angles. To probe the high eccentricity
values close to 1, we use binary simulations to compute the
$p(\gamma_{true,i}|e_{i})$ for eccentricities from 0 to 1 with a step of 0.001
and $\gamma_{true,i}$ with a step of $0.5^{\circ}$. The details of these terms
can be found in Hwang et al. (2022).

We then use emcee to derive the posterior distributions. In the MCMC run, we
numerically compute the two-dimensional integral in Eq. A1 with equal spacings
of $\Delta e=0.001$ and $\Delta\gamma=0.5^{\circ}$. Uninformative flat priors
are adopted for $p(e_{0})$, $p(e_{1})$, and $p(F_{\rm ecc})$, and we further
require that $e_{1}>e_{0}$. We use the Gaussian move as the proposal function
in the MCMC to mitigate the effect of finite $\Delta e$ in the integral
calculation. The resulting posterior distributions for 400-1000 AU wide twin
binaries are shown in Fig. 1 (Foreman-Mackey, 2016).

Figure 1: The posterior distributions of the eccentricity distribution model.
The solid lines are the most probable values of the marginalized
distributions, and the dashed lines are the highest posterior density interval
that includes 64-per cent of the area.

## References

  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/AAF741
  * Artymowicz et al. (1991) Artymowicz, P., Clarke, C. J., Lubow, S. H., & Pringle, J. E. 1991, ApJ, 370, L35, doi: 10.1086/185971
  * Artymowicz & Lubow (1994) Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651, doi: 10.1086/173679
  * Bally (2016) Bally, J. 2016, ARA&A, 54, 491, doi: 10.1146/ANNUREV-ASTRO-081915-023341
  * Bate & Bonnell (1997) Bate, M. R., & Bonnell, I. A. 1997, MNRAS, 285, 33, doi: 10.1093/MNRAS/285.1.33
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637x/823/2/102
  * Cuadra et al. (2009) Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C. 2009, MNRAS, 393, 1423, doi: 10.1111/J.1365-2966.2008.14147.X
  * Dittmann & Ryan (2022) Dittmann, A. J., & Ryan, G. 2022, MNRAS, doi: 10.1093/mnras/stac935
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Duffell et al. (2020) Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2020, ApJ, 901, 25, doi: 10.3847/1538-4357/abab95
  * El-Badry & Rix (2018) El-Badry, K., & Rix, H.-W. 2018, MNRAS, 480, 4884, doi: 10.1093/mnras/sty2186
  * El-Badry et al. (2021) El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, MNRAS, 506, 2269.  https://arxiv.org/abs/2101.05282
  * El-Badry et al. (2019) El-Badry, K., Rix, H.-W., Tian, H., Duchêne, G., & Moe, M. 2019, MNRAS, 489, 5822.  https://arxiv.org/abs/1906.10128
  * El-Badry et al. (2018) El-Badry, K., Ting, Y.-S., Rix, H.-W., et al. 2018, MNRAS, 476, 528, doi: 10.1093/mnras/sty240
  * Fabricius et al. (2021) Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5, doi: 10.1051/0004-6361/202039834
  * Farris et al. (2014) Farris, B. D., Duffell, P., Macfadyen, A. I., & Haiman, Z. 2014, ApJ, 783, 134, doi: 10.1088/0004-637X/783/2/134
  * Fezenko et al. (2022) Fezenko, G. B., Hwang, H.-C., & Zakamska, N. L. 2022, MNRAS, 511, 3881, doi: 10.1093/MNRAS/STAC309
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/JOSS.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1.  https://arxiv.org/abs/2012.01533
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Hamilton (2022) Hamilton, C. 2022, ApJ, 929, L29, doi: 10.3847/2041-8213/ac6600
  * Hartman & Lépine (2020) Hartman, Z. D., & Lépine, S. 2020, ApJS, 247, 66, doi: 10.3847/1538-4365/ab79a6
  * Heath & Nixon (2020) Heath, R., & Nixon, C. 2020, A&A, 641, A64.  https://arxiv.org/abs/2007.11592
  * Hwang et al. (2020) Hwang, H.-C., Hamer, J. H., Zakamska, N. L., & Schlaufman, K. C. 2020, MNRAS, 497, 2250.  https://arxiv.org/abs/2007.03688
  * Hwang et al. (2021) Hwang, H. C., Ting, Y. S., Schlaufman, K. C., Zakamska, N. L., & Wyse, R. F. 2021, MNRAS, 501, 4329, doi: 10.1093/mnras/staa3854
  * Hwang et al. (2022) Hwang, H.-C., Ting, Y.-S., & Zakamska, N. L. 2022, MNRAS.  https://arxiv.org/abs/2111.01789
  * Kaib & Raymond (2014) Kaib, N. A., & Raymond, S. N. 2014, ApJ, 782, 60, doi: 10.1088/0004-637X/782/2/60
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Lucy & Ricco (1979) Lucy, L. B., & Ricco, E. 1979, AJ, 84, 401, doi: 10.1086/112434
  * Miranda et al. (2017) Miranda, R., Muñoz, D. J., & Lai, D. 2017, MNRAS, 466, 1170, doi: 10.1093/MNRAS/STW3189
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Ochi et al. (2005) Ochi, Y., Sugimoto, K., & Hanawa, T. 2005, ApJ, 623, 922, doi: 10.1086/428601
  * Oh et al. (2017) Oh, S., Price-Whelan, A. M., Hogg, D. W., Morton, T. D., & Spergel, D. N. 2017, AJ, 153, 257, doi: 10.3847/1538-3881/aa6ffd
  * Pichardo et al. (2005) Pichardo, B., Sparke, L. S., & Aguilar, L. A. 2005, MNRAS, 359, 521.  https://arxiv.org/abs/0501244
  * Ragusa et al. (2020) Ragusa, E., Alexander, R., Calcino, J., et al. 2020, MNRAS, 499, 3362, doi: 10.1093/MNRAS/STAA2954
  * Reipurth & Mikkola (2012) Reipurth, B., & Mikkola, S. 2012, Nature, 492, 221, doi: 10.1038/nature11662
  * Roedig et al. (2011) Roedig, C., Dotti, M., Sesana, A., Cuadra, J., & Colpi, M. 2011, MNRAS, 415, 3033, doi: 10.1111/J.1365-2966.2011.18927.X
  * Schoettler et al. (2020) Schoettler, C., De Bruijne, J., Vaher, E., & Parker, R. J. 2020, MNRAS, 495, 3104, doi: 10.1093/mnras/staa1228
  * Shatsky (2001) Shatsky, N. 2001, A&A, 380, 238, doi: 10.1051/0004-6361:20011401
  * Shi et al. (2012) Shi, J.-M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118, doi: 10.1088/0004-637X/749/2/118
  * Soderhjelm (2007) Soderhjelm, S. 2007, A&A, 463, 683, doi: 10.1051/0004-6361:20066024
  * Tian et al. (2020) Tian, H.-J., El-Badry, K., Rix, H.-W., & Gould, A. 2020, ApJS, 246, 4, doi: 10.3847/1538-4365/ab54c4
  * Tiede et al. (2020) Tiede, C., Zrake, J., MacFadyen, A., & Haiman, Z. 2020, ApJ, 900, 43, doi: 10.3847/1538-4357/aba432
  * Tokovinin (2020) Tokovinin, A. 2020, MNRAS, 496, 993.  https://arxiv.org/abs/2004.06570
  * Tokovinin & Kiyaeva (2016) Tokovinin, A., & Kiyaeva, O. 2016, MNRAS, 456, 2070, doi: 10.1093/mnras/stv2825
  * Tokovinin & Moe (2020) Tokovinin, A., & Moe, M. 2020, MNRAS, 491, 5158.  https://arxiv.org/abs/1910.01522
  * Tokovinin (1998) Tokovinin, A. A. 1998, Astronomy Letters, 24, 178
  * Tokovinin (2000) —. 2000, A&A, 360, 997
  * Young & Clarke (2015) Young, M. D., & Clarke, C. J. 2015, MNRAS, 452, 3085, doi: 10.1093/mnras/stv1512