# Different Planetary Eccentricity-Period (PEP) Distributions of Small- and
Giant-Planets

Dolev Bashi Astrophysics Group, Cavendish Laboratory, University of Cambridge,
JJ Thomson Avenue, Cambridge CB3 0HE, UK  School of Physics and Astronomy, Tel
Aviv University, Tel Aviv, 6997801, Israel  Tsevi Mazeh School of Physics and
Astronomy, Tel Aviv University, Tel Aviv, 6997801, Israel  Simchon Faigler
School of Physics and Astronomy, Tel Aviv University, Tel Aviv, 6997801,
Israel

###### Abstract

We used the database of $1040$ short-period ($1\leq P<200$ days) exoplanets
radial-velocity (RV) orbits to study the planetary eccentricity-period (PEP)
distribution. We first divided the sample into low- and high-mass exoplanet
sub-samples based on the distribution of the (minimum) planetary masses, which
displays a clear two-Gaussian distribution, separated at
$0.165M_{\mathrm{J}}$. We then selected $216$ orbits, low- and high-mass
alike, with eccentricities significantly distinct from circular orbits. The
$131$ giant-planet eccentric orbits display a clear upper envelope, which we
model quantitatively, rises monotonically from zero eccentricity and reaches
an eccentricity of 0.8 at $P\sim 100$ days. Conversely, the $85$ low-mass
planetary orbits display a flat eccentricity distribution between $0.1$ and
$0.5$, with almost no dependence on the orbital period. We show that the
striking difference between the two PEP distributions is not a result of the
detection technique used. The upper envelope of the high-mass planets, also
seen in short-period binary stars, is a clear signature of tidal
circularization, which probably took place inside the planets, while the
small-planet PEP distribution suggests that the circularization was not
effective, probably due to dynamical interactions with neighboring planets.

planets and satellites: fundamental parameters — planets and satellites:
dynamical evolution and stability — methods: statistical

††software: astropy (Astropy Collaboration et al., 2013, 2018), NumPy (Harris
et al., 2020), Matplotlib (Hunter, 2007), scikit-learn (Pedregosa et al.,
2011), SciPy (Virtanen et al., 2020), emcee (Foreman-Mackey et al., 2013)
00footnotetext: Corresponding author: Dolev Bashi<EMAIL_ADDRESS>
##  1 Introduction

The detection of thousands of exoplanets has ushered in a new era in our
understanding of planetary systems, revealing a diverse array of planetary
architectures (Mishra et al., 2023) and compositions that challenge
traditional theories of planet formation and orbital dynamics (e.g., Lissauer
et al., 2011; Zhu & Dong, 2021; Mordasini & Burn, 2024; Yang & Jin, 2024). A
central aspect in the study of these planetary systems is the role of tidal
interaction of the planet with the central star and the dynamical interaction
with neighbouring planets in shaping the orbits and eccentricities of
exoplanets (e.g., Rasio et al., 1996; Matsumura et al., 2010). Therefore, the
observed planetary eccentricity distribution attracted much attention in the
past, as it might help distinguish between the different theories of exoplanet
formation and dynamical evolution (see a study with a broader perspective by
Ke-ting et al., 2024).

For example, Kipping (2013) and Shabram et al. (2016) fitted the shape of the
eccentricity distribution of planets, either with a $\beta$-function or with a
two-component Gaussian mixture model; see also Van Eylen et al. (2019), who
considered the different eccentricity distributions of single and multiple
planets, and Dong et al. (2021), who studied the distribution of planets
discovered by TESS in its first year of operation.

The dependence of the eccentricity on stellar metallicity and multiplicity,
and in particular on the planetary mass, was considered by a few studies
(e.g., Xie et al., 2016; Mills et al., 2019; Ribas & Miralda-Escudé, 2007; He
et al., 2020; An et al., 2023; Rosenthal et al., 2024). For example, Kane et
al. (2012) used the transit durations of the Kepler data to conclude that the
eccentricity distribution is consistent with that found from radial-velocity
(RV) surveys to a high degree of confidence, and the mean eccentricity of the
Kepler candidates decreases with decreasing planet size, indicating that
smaller planets are preferentially found in low-eccentricity orbits. We will
come back to this distinction in the present work.

Quite a few works considered the eccentricity dependence on the orbital period
(hereafter, Planetary Eccentricity-Period (PEP) distribution). Early works
appeared already in 2011 (Pont et al., 2011; Husnoo et al., 2012), when the
number of planetary orbits was quite small. Later works (e.g. Giacalone et
al., 2017; Dawson & Johnson, 2018; Jackson et al., 2023) have further aimed to
confront the observed PEP with the high-eccentricity migration (HEM) scenario,
which assumes that giant planets were formed with high orbital eccentricities
at large distances from their parent star, and then migrated inward by tidal
interaction induced by the high eccentricity (Socrates et al., 2012; Bonomo et
al., 2017; Dawson & Johnson, 2018; Jackson et al., 2023).

A few studies e.g., Giacalone, Matsakos, & Königl (2017) and Dawson & Johnson
(2018, in their review paper, Figure 4 in particular); see also (Jackson et
al., 2023), focused on the upper envelope of the PEP as an important feature
of the planetary population, which is the subject of this study.

The number of derived planetary orbits grew substantially in the last few
years, especially when high-precision spectrographs dedicated to the RV study
of exoplanets went into operation (e.g., Cosentino et al., 2012; Quirrenbach
et al., 2014; Pepe et al., 2021). Therefore, it is time to revisit the PEP
distribution with the newly accumulated orbits and its dependence on the
planetary mass. We chose to use only RV orbits as the orbital eccentricities
and planetary masses are not well constrained by transit light curves alone
(but see, for example, Kipping, 2008; Ford & Rasio, 2008; Dawson & Johnson,
2012; Kipping, 2014). We fit an upper-envelope function to the PEP
distribution, using a formalism we developed recently (Bashi, Mazeh, &
Faigler, 2023) to analyze the eccentricity distribution and circularization
period cut-off of binary stars. Following previous studies, we compare the PEP
distributions of the small and giant planets, finding substantially different
shapes.

Section 2 presents our sample selection aimed to produce an unbiased subset of
low- and high-mass planets. Section 3 presents the PEP distribution of the
small and giant planets, section 4 discusses the results and section 5
summarizes this study.

##  2 The sample

For the study of the PEP distribution, we used the Exoplanet Archive
’Composite’ catalog (Akeson et al.,
2013),111https://exoplanetarchive.ipac.caltech.edu/index.html a comprehensive
and well-maintained database that aggregates data from multiple sources,
including transit surveys and RV measurements. We opted for this catalog due
to its extensive coverage and reliability, which has been cross-verified
against other databases and studies (e.g., Perryman, 2018; Bashi et al., 2018;
Adibekyan, 2019). As of the 28th of August, 2023, the catalog includes orbital
periods and eccentricity for $1987$ planets that were observed and confirmed
by RV measurements. In this study we focused on the $1040$ ’short-period’
$1\leq P<200$ days orbits.

While the Exoplanet Archive ’Composite’ Planet Data Table offers useful
preliminary insights into the population-level properties of exoplanets, such
studies require careful scrutiny to ensure statistical robustness. We
therefore have conducted a rigorous examination of the sample selection and
potential observational biases, as detailed in the following sub-sections.
This includes the difference in exoplanet detection methodologies and the
position on the color-magnitude diagram (CMD) of the parent stars.

###  2.1 Division between low- and high-mass planets

To divide the sample into small- and large-mass planets, we plot in Figure 1
the (minimum) mass histogram of the planets of all $1040$ RV orbits, eccentric
and circular. The figure clearly reveals two distinct populations of planets,
which we have modeled using a Gaussian Mixture Model (GMM) (Pedregosa et al.,
2011).222https://scikit-learn.org/stable/modules/generated/
sklearn.mixture.GaussianMixture.html One Gaussian component is centered around
$M_{\rm p}\simeq 1M_{\rm J}$, representing our high-mass (giant) planets. The
other Gaussian component is centered at $M_{\rm p}\simeq
10M_{\rm\oplus}=0.03M_{\rm J}$, representing the low-mass (Super-Earth and
Sub-Neptune) planet sample. Here, $M_{\rm J}$ and $M_{\rm\oplus}$ denote the
masses of Jupiter and Earth, respectively.

Figure 1: Minimum-mass histogram (blue bins) of all planets with RV orbits,
eccentric and circular alike. Dashed curves represent our best Gaussian-
Mixture Model of the small- and giant-planet samples. The mean and standard-
deviation of the small planets (in log $M_{\rm J}$) are $-1.53$
($=9.4M_{\rm\oplus}$) and $0.40$, and $0.00$ ($=1M_{\rm J}$) and $0.43$ for
the giant planets. The red dotted line marks the intersecting point between
the two Gaussians at $M_{\rm p}=0.165M_{\rm J}$. Overlaid, using gray bins, is
the distribution of the transiting planets in our sample.

A significant result of the GMM fitting is the identification of a division at
$M_{\rm p}=0.165M_{\rm J}$, around the Neptunian desert (e.g., Ida & Lin,
2004; Mazeh, Holczer, & Faigler, 2016), which we use to separate the two
populations into $470$ small and $570$ giant planets.

To show that our division is not dependent on the detection method, we
overlay, using gray bins, the distribution of the transiting planets in our
sample. In these cases, the planetary mass estimate derived from the RV orbit,
$M_{\mathrm{p}}\sin i$, is similar to the actual planetary mass, as $\sin i$
is close to unity. One can see that the gap and the division between the small
and giant planets is similar to the gap derived for the overall sample.

###  2.2 Significance of the orbital eccentricity

As pointed out by Lucy & Sweeney (1971) and more recently by Hara et al.
(2019) small orbital eccentricities are not always reliable; such eccentric
orbits might be indistinguishable from circular ones. Consequently, we divided
our sample into ’eccentric’ orbits, for which $e>2\Delta e$, where $\Delta e$
is the reported lower uncertainty of the eccentricity and ’circular’ orbits,
for which $e\leq 2\Delta e$, or in cases where the detection paper assumed a
circular solution. We were left with $762$ circular and $301$ eccentric
orbits.

###  2.3 Planet-host stars selection

Next, we aimed to reduce possible biases of our analysis caused by the
characteristics of the host stars, which might significantly influence the
observed orbital parameters of their planets, either by inherent physical
interactions or by some detection sensitivities (see, for example Banerjee et
al., 2024; Tang et al., 2024, whose studies are based on the characteristics
of the planet hosts). To do so, Figure 2 presents a color-magnitude-diagram
(CMD) of the planet host-star samples, using Gaia G, and BP-RP magnitudes,
separated to star hosting small and giant planets (blue points). The
background gray-scale sample is added in both panels for reference and is
based on the Gaia Catalog of Nearby Stars (GCNS Gaia Collaboration et al.,
2021), a clean catalog of objects within $100$ pc of the Sun. As expected, the
giant-planet hosts tend towards the early-type stars and evolved stars. To
match the common region on the CMD where most small- and giant-planet host
stars are situated, we selected only stars inside the red polygon333The
polygon boundaries are as follows (BP-RP, Abs. G): (0.70, 2.5), (0.92, 2.5),
(1, 4), (2, 7.8), (2, 9.3), (0.70, 4.7), (0.70, 2.5). marked in Figure 2. We
were left with a sample of $85$ small and $131$ giant planets.

Figure 2: Color-magnitude-diagram (CMD) of all planet-hosting stars with
eccentric orbits in our sample. Left: low-mass planets, Right: high-mass
planets. Background gray-scale sample is added for reference, using the Gaia
Catalog of Nearby Stars (GCNS). Red polygons in both panels mark our selected
sample of planet-host stars with similar stellar characteristics. See text for
further description.

These are the final two samples of eccentric planets we chose to analyze.

###  2.4 Discovery method

While we exclusively focus on planetary orbits derived from RV measurements,
many of these planets were initially identified through transit surveys before
subsequent RV follow-ups. The specific statistical characteristics of transit-
detected planets versus those discovered through RV surveys vary based on the
distinct survey protocols and their respective detection thresholds. The
transit detection method can exhibit a bias toward planets with shorter
orbital periods, for example.

In our case, $49$ and $36$ of the small planets were first detected by the RV
and transit methods, respectively, while in the case of the giant planets, the
numbers are $64$ and $67$, respectively. In the next section, when exploring
the PEP two distributions, we will also consider the discovery method of each
planet, showing that the different methods have no impact on our results.

We list in tables 1 and 2 the full sample of low-mass and high-mass planet
parameters of the orbital motion and host properties used in this work to
assess the PEP distributions.

Table 1: Small Planets ($M_{\rm p}<0.165M_{\rm J}$) properties used in the analysis of the PEP distributions. The table includes properties of orbital period, eccentricity and lower eccentricity uncertainties, planet mass, number of other planets in the system, discovery method and Gaia BP-RP and Abs. G host magnitudes. Sources are sorted by orbital period.  Name | Period | $e$ | $\Delta e$ | $M_{\rm p}$ | Np | Discovery Method | BP-RP | Abs. G
---|---|---|---|---|---|---|---|---
| [days] |  |  | ${\rm[}M_{\rm J}{\rm]}$ |  |  | [mag] | [mag]
GJ 2030 b | $2.3813$ | $0.239$ | $0.058$ | $0.01$ | $2$ | RV | $0.91$ | $3.17$
TOI-1288 b | $2.6998$ | $0.064$ | $0.015$ | $0.13$ | $2$ | Transit | $0.98$ | $5.15$
Kepler-411 b | $3.0052$ | $0.146$ | $0.004$ | $0.08$ | $4$ | Transit | $1.23$ | $6.30$
HD 134060 b | $3.2696$ | $0.450$ | $0.040$ | $0.03$ | $2$ | RV | $0.77$ | $4.25$
TOI-1272 b | $3.3160$ | $0.338$ | $0.062$ | $0.08$ | $2$ | Transit | $1.11$ | $5.90$
GJ 676 A d | $3.6008$ | $0.192$ | $0.079$ | $0.01$ | $4$ | RV | $1.90$ | $7.85$
HD 219828 b | $3.8348$ | $0.091$ | $0.025$ | $0.06$ | $2$ | RV | $0.80$ | $3.59$
TOI-1062 b | $4.1130$ | $0.177$ | $0.064$ | $0.03$ | $2$ | Transit | $1.00$ | $5.43$
HD 18599 b | $4.1374$ | $0.340$ | $0.090$ | $0.08$ | $1$ | Transit | $1.06$ | $5.81$
HD 12235 b | $4.3631$ | $0.376$ | $0.154$ | $0.03$ | $2$ | RV | $0.76$ | $3.22$
Table 2: Same as Table 1 but for the giant planets sample ($M_{\rm p}\geq 0.165M_{\rm J}$).  Name | Period | $e$ | $\Delta e$ | $M_{\rm p}$ | Np | Discovery Method | BP-RP | Abs. G
---|---|---|---|---|---|---|---|---
| [days] |  |  | ${\rm[}M_{\rm J}{\rm]}$ |  |  | [mag] | [mag]
WASP-12 b | $1.0914$ | $0.045$ | $0.004$ | $1.47$ | $1$ | Transit | $0.75$ | $3.42$
HAT-P-23 b | $1.2129$ | $0.106$ | $0.044$ | $2.09$ | $1$ | Transit | $0.83$ | $4.35$
HAT-P-36 b | $1.3273$ | $0.063$ | $0.023$ | $1.85$ | $1$ | Transit | $0.87$ | $4.75$
WASP-5 b | $1.6284$ | $0.038$ | $0.018$ | $1.58$ | $1$ | Transit | $0.84$ | $4.51$
HD 86081 b | $2.1378$ | $0.012$ | $0.005$ | $1.48$ | $1$ | RV | $0.76$ | $3.53$
WASP-140 b | $2.2360$ | $0.047$ | $0.004$ | $2.44$ | $1$ | Transit | $0.98$ | $5.49$
HAT-P-16 b | $2.7760$ | $0.036$ | $0.004$ | $4.19$ | $1$ | Transit | $0.72$ | $4.02$
HAT-P-20 b | $2.8753$ | $0.015$ | $0.005$ | $7.25$ | $1$ | Transit | $1.42$ | $6.72$
HAT-P-13 b | $2.9162$ | $0.013$ | $0.004$ | $0.85$ | $2$ | Transit | $0.86$ | $3.47$
HD 46375 b | $3.0236$ | $0.063$ | $0.026$ | $0.23$ | $1$ | RV | $1.00$ | $5.36$

##  3 PEP DISTRIBUTION FOR SMALL AND GIANT PLANETS

Figure 3 shows the PEP distributions of the low- and high-mass planets.

Figure 3: Separate PEP diagrams for the small and giant planets. Triangles and
squares represent planets discovered by the transit and RV methods,
respectively. Red open circles mark planets in known multi-planet systems. The
open triangles in the small-planet diagram present orbits from the study of
Feng et al. (2022), which did not publish some of their orbital elements, nor
the RV measurements. The red dashed line in the giant-planet diagram marks our
upper-envelope best-fit model, marked by a vertical dotted line, based on
planets with an orbital period shorter than $100$ days. The red area marks the
transition region of $\pm\delta$ width along the envelope (see text).Gray
dashed curve and the corresponding gray region mark the expected constant
angular momentum evolutionary track, given a final orbital period of $2\leq
P_{\mathrm{final}}<5$ days.

The difference between the two PEP diagrams is striking. The high-mass orbits
display a clear upper-envelope shape that monotonically rises from zero
eccentricity at the orbital period of $\sim 1.5$ days and reaches eccentricity
of 0.8 at a period of $\sim 100$ days (see below). On the other hand, the low-
mass planetary orbits display a flat eccentricity distribution between $\sim
0.1$ and $\sim 0.5$, with almost no dependence on the orbital period.

To see if the pronounced contrast in the PEP between small and giant planets
might be attributable to the detection method, we have segregated the data in
Figure 3 to distinguish between transit-detected planets (represented as
squares) and those discovered through RV surveys (represented as triangles).
No clear dependence on the detection method can be noticed (see below that the
numerical fitting of the RV and transit detected planets of the giant planets
yield similar upper envelopes).

The low-mass PEP distribution includes five orbits with high eccentricities,
clearly above the upper limit of the population at $\sim 0.5$. Interestingly,
all five systems have a distant sub-stellar companion, as detailed in Table 3,
while the overall fraction of systems with known companions in the small-
planet sample is $15$ out of $85$. This is based on the value marked by
’sy_pnum’ and listed in the Exoplanet Archive ’Composite’ Planet Data Table,
but note that the actual binary nature (separation, mass ratio) is not listed.
One might attribute these eccentricities to dynamical pumping induced by the
distant companions (e.g., Mazeh & Shaham, 1979; Fabrycky & Tremaine, 2007, see
discussion).

Table 3: Five inner small planets with high eccentricities and their distant
companions.

Name |  Inner low-mass planet |  Sub-stellar distant companion | Reference
---|---|---|---
| Period | $e$ | $M_{\rm p}$ | Period | $e$ | $M_{\rm p}$ |
| [days] |  | $[M_{\rm J}]$ | [days] |  | $[M_{\rm J}]$ |
GJ-3222 | 10.667 | 0.93±0.06 | 0.036±0.008 | 1429 | 0.537±0.052 | 52±5 | [1]
Kepler-1656 | 31.562 | 0.838±0.33 | 0.15 ±0.02 | 1919 | 0.527±0.051 | 0.40±0.09 | [2]
HD 144899 | 40.439 | 0.82±0.15 | 0.064±0.018 | 8113 | 0.355±0.019 | 100±10 | [1]
GJ 2056 | 69.971 | 0.72±0.10 | 0.051±0.011 | 2982 | 0.81±0.018 | 0.44±0.05 | [1]
HD 13724 | 142.228 | 0.651±0.083 | 0.16±0.023 | 15784 | 0.327±0.029 | 36±2 | [1]

Note. — References: [1] Feng et al. (2022); [2] Angelo et al. (2022).

We mark by red circle in the left panel of the low-mass planets sample of
Figure 3 the sources that are part of known multi-planet systems,
$N_{\mathrm{p}}>1$, based on the value marked by ’sy_pnum’ and listed in the
Exoplanet Archive ’Composite’ Planet Data Table. The majority of planets in
the low-mass sample ($64$ out of $85$) have at least one additional detected
planet in the same stellar system. In contrast, for the high-mass planet
sample, the fraction of ones with other planets orbiting the same host star is
much smaller ($39$ out of $131$).

To quantify the observed difference in the two PEP diagrams, we used a KS-test
to compare the eccentricity distributions of planets with orbital periods
shorter and longer than $10$ days. The p-values obtained were $0.88$ and
$9.2\times 10^{-4}$ for the low- and high-mass planetary samples respectively.
These results further confirm the observed trend: while the eccentricity
distributions of low-mass planets remain consistent across the orbital period,
those of high-mass planets exhibit significant variations depending on the
period cut. This suggests a distinct dependency of eccentricity on orbital
period for planets of higher mass.

In the following we fit an upper envelope to the giant-planet PEP diagram.

###  3.1 Modeling the upper envelope

To numerically characterize the upper envelope seen in the observed high-mass
PEP, we apply the methodology constructed by Bashi, Mazeh, & Faigler (2023),
which defined a functional form of two parameters:

$f(P)=\left\\{\begin{aligned}
&1-\left(\frac{P_{\mathrm{cut}}}{P}\right)^{\frac{1}{\tau}},&&\text{if
}P>P_{\mathrm{cut}}\ ,\\\ &0,&&\text{if }P\leq P_{\mathrm{cut}}\
,\end{aligned}\right.$ (1)

where $P_{\rm cut}$ is the cutoff period and $\tau$ is a dimensionless
parameter that determines the slope of the envelope. The function gets a zero
value at $P_{\rm cut}$ for any value of $\tau$, and asymptotically rises to
unity as $P$ gets substantially longer than $P_{\rm cut}$.

We describe the probability density function (PDF) for planets in the ($\log
P$, $e$) plane using a modified Fermi function,
$\mathcal{F}_{\mathrm{PDF}}(P,e;\,{P_{\mathrm{cut}}},\tau,\delta)$. The PDF
converges to zero above the envelope, and to a positive constant below it,
with a transition region of width $\delta$ along the envelope.

For each set of parameters of the PDF—${P}_{\rm cut}$, $\tau$ and $\delta$, we
derive the likelihood of the sample as

$\mathcal{L}=\prod_{i=1}^{N}\mathcal{F}_{\mathrm{PDF}}({P}_{i},e_{i};{P}_{\rm
cut},\tau,\delta)\ .$ (2)

where $e_{i}$ and $P_{i}$ are the eccentricity and period of the $i-$th
planet. We search the three parameters that maximize the likelihood of the
observed sample.

###  3.2 The upper Envelope of the high-mass PEP diagram

We applied our algorithm to the PEP of the high-mass planets sample, using an
MCMC routine with $50$ walkers and $10^{4}$ steps, with uninformative priors
for the three parameters. We used the Python emcee package (Foreman-Mackey et
al., 2013) to find the best parameter values and their uncertainties. Our
prior boundaries and best-fit posterior values are given in Table 4. The red
dashed line in the high-mass diagram of Figure 3 shows our best-fit upper
envelope. The red area in the figure marks the transition region of
$\pm\delta$ width along the upper envelope.

A single outlier in the figure, somewhat above our fitted envelope, is that of
CoRoT-20 b (Deleuil et al., 2012) with $e=0.59\pm 0.02$ and $P=9.2429\pm
0.0003$ days. Interestingly, in a later work, Rey et al. (2018) have found
this system hosts another sub-stellar eccentric companion, CoRoT-20 c, with
$M_{\mathrm{p}}\sin i=17\pm 1M_{\mathrm{J}}$, $e=0.60\pm 0.03$ and $P=4.59\pm
0.05$ yr, suggesting a dynamical eccentricity pumping induced by the distant
companion (see above).

In addition, there are three outliers located left to our upper envelope fit
with orbital periods shorter than $2$ days and significant eccentricities.
These planets are WASP-12 b ($e=0.0447\pm 0.0045$, Hebb et al., 2009),
HAT-P-36 b ($e=0.063\pm 0.023$, Wang et al., 2019) and HAT-P-23 b ($e=0.106\pm
0.044$, Bakos et al., 2011). The small but significant eccentricities for
planets so close to their host stars are unexpected and might be interesting
candidates for further exploration. We suggest, as discussed in Arras et al.
(2012) and more recently pointed out by Yee et al. (2020) in the case of
WASP-12 b, that these eccentricities are likely artifacts of the tidal
interactions induced by the planet on its host rather than true orbital
characteristics. In their work, Arras et al. (2012) have demonstrated how the
tidal fluid flow and the epicyclic motion of a slightly eccentric orbit
produce an RV signal at half the orbital period. These signals can
consequently mimic the effects of a non-zero eccentricity, leading to the
erroneous conclusion that the orbit is eccentric when it is, in fact,
circular.

We further note, that these three planets were detected as part of large
ground-based transit search surveys: HAT (Bakos et al., 2004) and WASP
(Pollacco et al., 2006). Given the significantly larger number of stars these
surveys can monitor compared to radial velocity (RV) surveys, we find it
reasonable that the shortest orbital period giant planets were discovered
first by transiting surveys where the transit detection probability is also
the highest given the geometric probability.

Table 4: Prior and posterior values of upper envelope model of the high-mass planet sample. Parameter | Prior | Posterior
---|---|---
$P_{\mathrm{cut}}$ | $\log\mathcal{U}(-1,1)$ | $1.66^{+0.24}_{-0.25}~{}\mathrm{day}$
$\tau$ | $\log\mathcal{U}(0,2)$ | $2.50^{+0.29}_{-0.22}$
$\delta$ | $\mathcal{U}(0.01,0.3)$ | $0.079^{+0.032}_{-0.023}$

To look for potential bias in the derived envelope due to the planetary
discovery method, we conducted separate fittings of upper envelopes for the
subsets of high-mass planets identified by the transit and RV methods. We
found the upper-envelope parameters of the two sub-samples to be consistent
with each other. The $\log P_{cut}$ values were found to be
$0.15^{+0.08}_{-0.12}$ and $0.35^{+0.17}_{-0.12}$ for the transit and RV
detected sub-samples, respectively, indicating $\sim 1.5\sigma$ difference,
probably due to the bias of transit detections towards short periods. For the
slope of the envelope, $\log\tau$, we found $0.94^{+0.19}_{-0.14}$ and
$1.21^{+0.26}_{-0.15}$, respectively, within $\sim 1\sigma$ of each other.

To further explore the impact of the uncertainties of the orbital parameters
on our upper envelope posterior estimates, we employed a bootstrap approach.
We randomly drew eccentricities for each planet in the sample, using a normal
distribution with the mean and standard deviation corresponding to that
planet’s reported mean and uncertainty values, as listed in the Exoplanet
Archive ’Composite’ Planet Data Table. In this way, we generated revised
orbital parameters for our 131 high-mass planets and then applied the
algorithm of finding the parameters of the upper envelope. By this way we got
distributions of the values of the three parameters, based on the eccentricity
uncertainties. Overall, we found that the bootstrap uncertainties for the
upper envelope posterior values were an order of magnitude smaller than the
estimated uncertainties given by the final MCMC model. This finding affirms
that our results are insensitive to the uncertainties of the planetary orbital
eccentricities.

##  4 Discussion

This study considered $1040$ planetary orbits with periods in the $1\leq
P<200$ days range with RV solutions, as appeared in the Exoplanet Archive
’Composite’ catalog database. The sample was not obtained by any systematic
search but is a compilation of results reported by many studies. Therefore, it
is subject to different observational detection thresholds and to the ability
of the different projects to derive the orbital eccentricity correctly.
Nevertheless, the trends found are highly significant and do not depend on the
completeness of the sample.

Of the 1040 orbits, only 301 have significant eccentricity, with $e>2\Delta e$
that can be used to explore the PEP distribution. After considering possible
biases in detection method and host properties, we were left with a sample of
$131$ high-mass orbits. This sample displays a clear upper envelope, which we
model using our novel formalism (Bashi, Mazeh, & Faigler, 2023). We find that
the upper envelope monotonically rises from zero eccentricity at an orbital
period of $\sim 1.7$ days and reaches an eccentricity of 0.8 at a period of
$100$ days. The $85$ low-mass planetary orbits, on the other hand, display a
flat eccentricity distribution between $\sim 0.1$ and $\sim 0.5$, with almost
no dependence on the orbital period. As seen from the diagram and supported by
our separate analysis, these results are independent of the planetary
discovery method.

###  4.1 Tidal Circularization—small versus giant planets

The giant PEP distribution has a typical upper-envelope signature of tidal
circularization, similar to the eccentricity distribution of spectroscopic
binaries, as analyzed recently by Bashi, Mazeh, & Faigler (2023), and was
interpreted as the result of circularization, either during the stellar
lifetime or at the pre-main-sequence stellar phase. The small-planet PEP
distribution, on the other hand, is very far from having the shape of a sample
that went through tidal circularization. Can we explain the substantially
different PEP distributions?

We start by comparing the circularization timescale $\tau_{\rm circ}$ of the
small- and giant-planets, using a formula derived from Goldreich & Soter
(1966) theory:

$\tau_{\rm circ}=\left(4/63\right)Q^{\prime}_{\rm p}\frac{M_{\rm
p}}{M_{*}}\left(\frac{a}{R_{\rm p}}\right)^{5}P_{\rm orb}\,,$ (3)

where $a$ is the orbital semi-major axis, $M_{*}$ is the stellar mass, and
$M_{\rm p}$, $R_{\rm p}$ and $P_{\rm orb}$ are the planet mass, radius and
orbital period, respectively. The tidal quality factor $Q^{\prime}_{\rm p}$ is
a dimensionless parameter that quantifies the efficiency of energy dissipation
due to tidal forces (Goldreich & Soter, 1966).

Consider, for example, an orbit of $4$ days, for which $a\sim 10R_{\odot}$.
For a giant planet, with $M_{p}/M_{\odot}\sim 10^{-3}$ and $R_{p}\sim
10^{-1}R_{\odot}$ and a solar-type star one gets

${\rm Giant\,Planets:}\,\,\frac{M_{\rm p}}{M_{*}}\left(\frac{a}{R_{\rm
p}}\right)^{5}\sim 10^{7}\,,$ (4)

while for small-mass planets, $M_{p}/M_{\odot}\sim 3\times\,10^{-5}$ and
$R_{p}\sim 10^{-2}R_{\odot}$ (e.g., Chabrier et al., 2009; Bashi et al., 2017;
Ulmer-Moll et al., 2019), and therefore

${\rm Small\,Planets:}\,\,\frac{M_{\rm p}}{M_{*}}\left(\frac{a}{R_{\rm
p}}\right)^{5}\sim 3\,\times 10^{10}\,.$ (5)

Previous works have argued that the typical quality factor of giant planets is
$Q_{\rm p}\sim 10^{5}$, while rocky planets of the solar systems usually have
$Q_{\rm p}\sim 10^{2}$ (Goldreich & Soter, 1966; Murray & Dermott, 1999;
Clausen & Tilgner, 2015). We adopt these values, although they represent
order-of-magnitude estimation only (e.g., Matsumura et al., 2008; Papaloizou,
2011). When we put the corresponding values in Equation 3 we get

${\rm Giant\,Planets:}\,\,\tau_{\rm circ,g}\sim 10^{11}P_{\rm orb}\sim
10^{9}{\rm yrs}\,,$ (6)

while for small-mass planets,

${\rm Small\,Planets:}\,\,\tau_{\rm circ,s}\sim 3\times\,10^{11}P_{\rm
orb}\sim 3\times\,10^{9}{\rm yrs}\,.$ (7)

Note that the ratio between the two timescales, $\tau_{\rm circ,s}/\tau_{\rm
circ,g}\sim 3$ , is independent of the orbital period. The longer
circularization timescale of the small planets indicates a weaker tidal
interaction relative to that of the giant planets by a factor of $\sim 3$, all
other binary elements being equal. The two timescales are order-of-magnitude
estimations only, because of the spread of periods and masses of the two
populations.

The small factor difference between the estimations of the two types of
planets cannot explain the striking difference between the two PEP
distributions. If the small planets went through a similar tidal
circularization, we could expect their PEP distribution to show a similar
upper-envelope shape, with a slightly longer $P_{\rm cut}$ that compensates
for the factor $3$ weakness of their tidal interaction. Apparently, the small
planets were not able to go through tidal circularization.

On the other hand, if the dissipation is dominated by friction inside the
parent star, it would provide an easy explanation of the difference we have
observed between the PEP diagrams of low-mass and high-mass planets.
Obviously, low-mass planets do not perturb the star as much as high-mass
planets.

As suggested by previous studies (e.g., Van Eylen et al., 2019), the
eccentricity distribution of the small planets might be due to the fact that
many of the small planets were formed with neighboring planets, probably at
relatively large distances from their parent star. Supporting this statement
is the fact that most of the close-in low-mass planets in our sample reside in
multi-planet systems, as suggested by Figure 3 (see, for example Obertas et
al., 2023). These planets then migrated inwards, while probably dynamically
locked with their neighbouring planets (see, for example, the discussion by
Weiss et al., 2023). These ”multis” were under strong gravitational
interaction with their neighbors, which could inflate their eccentricities. On
the other hand, eccentricities too high, with, say, $e>0,5$, in one of the
orbits, can render the configuration dynamically unstable (e.g., Tremaine,
2015; Correia et al., 2020; Goldberg & Batygin, 2022; Obertas et al., 2023;
Ghosh & Chatterjee, 2023); see also Matsumoto et al. (2015).

Note, however, that not all small planets are known to reside in planetary
multiple systems. This could be because some of the neighboring planets are
undetectable, as their mass, inclination, and/or period push them below the
observational RV detection threshold. The eccentricity distribution as a
function of the multiplicity of the planetary systems (e.g., Van Eylen &
Albrecht, 2015; Van Eylen et al., 2019; Sagear & Ballard, 2023) could be key
evidence for understanding the role of the neighboring planets in planetary
eccentricity. This is however behind the scope of this paper.

As a final remark, we note that data irregularities and sparse observational
coverage can occasionally lead to inaccurate interpretations of orbital
configurations. For instance, what appears to be a single low-mass, highly
eccentric exoplanet could be two separate planets each in a more typical,
near-circular orbit (Wittenmyer et al., 2019). This kind of misinterpretation
can occur when signals from two closely spaced or overlapping planets are not
sufficiently sampled, causing them to be mistaken for a single eccentric
orbit.

###  4.2 Tidal Circularization— giant planets versus binaries

Interestingly, the $P_{\rm cut}$ of the giant-planet distribution, $\sim 1.8$
day, is significantly shorter than the one obtained for the spectroscopic
binaries, which is about $\sim 6$ days for similar primary stars. We note that
if we apply Equation 3 to binaries with stellar companions, replacing the
planetary mass and radius by stellar mass and radius, $M_{\rm comp}$ and
$R_{\rm comp}$, we get, assuming the same $Q$,

$\tau_{\rm circ}\propto\frac{M_{\rm comp}}{R_{\rm comp}^{5}}\,.$ (8)

This leads to an estimation of the stellar circularization timescale, for
solar mass $M_{\rm comp}\sim 10^{3}M_{\rm J}$ and radius $R_{\rm comp}\sim
10R_{\rm J}$ ,

$\tau_{\rm circ,stellar}\sim 10^{-2}\tau_{\rm circ,g}\,.$ (9)

This implies that the circularization processes are more effective in
binaries, and therefore could circularize systems with longer periods, if they
were active during the same period of time, the stellar lifetime, for example.
This is consistent with what has been observed.

###  4.3 Periastron period

As the circularization is most effective during periastron passage, we might
expect that the upper envelope represents a constant periastron distance.
Planets with too short periastron passage went through fast tidal interaction.
This means that for any given period the eccentricity has a maximum value,
parameterized by $P_{\mathrm{peri}}$, such that

$e_{\rm max}=1-\left(\frac{P_{\mathrm{peri}}}{P}\right)^{2/3}\,.$ (10)

This form is the same as our upper-envelope function, which is translated to
$1/\tau\sim 2/3$. Our finding of $1/\tau\sim 0.4$, a value not far from the
one expected, given the scatter of the sample and the simplifying assumptions.

###  4.4 Exploring the constant angular momentum track

As a final test, we wish to explore the expected ’evolutionary track’ of HEM
giant planets assuming they conserve their angular momentum. Given the period
distribution of circular hot Jupiters has a strong 3-day accumulation (e.g.
Cumming et al., 1999; Udry et al., 2003; Cumming et al., 2008), we chose to
use this circularized ’final’ orbital period, $P_{\mathrm{final}}$ to rewind
the hot Jupiter’s tidal evolution history using (Dawson & Johnson, 2018;
Jackson et al., 2023)

$P_{\mathrm{final}}=P(t)\left(1-e^{2}\right)^{3/2}\,.$ (11)

We plot in the right panel of Figure 3 the expected ’evolutionary track’ of
these objects given we find most circularized hot-Jupiters in the $2\leq
P_{\mathrm{final}}<5$ days limits. Clearly, these evolutionary tracks are
different from the derived upper envelope of the sample.

###  4.5 Highly eccentric planetary orbits

Five orbits of small planets do show high eccentricities larger than $0.6$.
Interestingly, all five systems are known to have a sub-stellar companion at
an orbital period of $1000$–$16,000$ days. Those companions might be
responsible for pump-in high eccentricity into the planetary orbit during the
formation or evolution of these planets (e.g., Fabrycky & Tremaine, 2007;
Correia et al., 2012; Petrovich et al., 2014). Obviously, one needs a
systematic search for distant companions of all small planets, not only the
ones with high eccentricities, to assess the contribution of the assumed
eccentricity pumping.

In fact, there are a few known giant planets with very eccentric orbits, with
periods longer than the range considered here. One of the most well-known
systems is that of HD 80606b (also seen in our PEP diagram) with a period of
$\sim 111$ days and $e=0.927\pm 0.012$ (Naef et al., 2001). HD 20782b is
another example, with a period of $\sim 600$ days and $e=0.97\pm 0.01$
(O’Toole et al., 2009). The most recently discovered one, HIP 66074b, with a
period of $\sim 300$ days and $e=0.948\pm 0.004$, has a distant companion that
might also have a sub-stellar mass (see discussion by Sozzetti et al., 2023).
For these cases, the periastron distance, $\sim 0.05$ au, corresponds to about
the histogram peak of circular orbits of giant planets at $\sim 4$ days.

##  5 Conclusion

The present PEP probably reflects both the primordial eccentricity and the
circularization processes. The peak of the period distribution of the circular
orbits at $3$–$4$ days apparently presents the giant planets whose orbits have
already been circularized, as suggested by the the HEM model (Socrates et al.,
2012; Bonomo et al., 2017; Giacalone et al., 2017; Dawson & Johnson, 2018;
Jackson et al., 2023). To find the primordial periods and eccentricities of
these orbits, one needs to integrate back the orbits with the equations that
govern the eccentricity and period evolution. Note, however, that a backward
model might be quite complex, as it might need to include a possible orbital
misalignment of the spin axes of the parent star and the planet (e.g.,
Mardling & Lin, 2002; Vick et al., 2023), unseen small planets, and
differential stellar rotation, which could complicate the model.

Finally, it is interesting to note that the eccentricity of the exoplanet
population revealed in the last three decades surprising features (see
extensive discussion in Biazzo et al., 2022), in contrast to the Solar system,
where all planets have almost circular orbits. The dependence of the
eccentricity on the planetary mass and orbital period was studied here,
suggesting that the giant planets went through tidal circularization processes
while the small planets did not. Other features, like the dependence of the
eccentricity on the parent star metallicity (An et al., 2023) and the
planetary multiplicity (Van Eylen et al., 2016), are also extremely important.
All these features carry with them many hints on the formation and evolution
of the exoplanet population.

##  6 Acknowledgments

We wish to thank the anonymous reviewer for an enlightening report. We extend
our heartfelt gratitude to Josh N. Winn and Daniel C. Fabrycky for adding
further comments and suggestions on an early version of this work. This
research was supported by Grant No. 2016069 of the United States-Israel
Binational Science Foundation (BSF) to TM, and Grant No. I-1498-303.7/2019 of
the German-Israeli Foundation for Scientific Research and Development (GIF) to
TM. This research is partially supported by Israel Ministry of Innovation,
Science and Technology grant 0004138. D.B. acknowledges the support of the
Blavatnik family and the British council as part of the Blavatnik Cambridge
Fellowship. This research has made use of the NASA Exoplanet Archive, which is
operated by the California Institute of Technology, under contract with the
National Aeronautics and Space Administration under the Exoplanet Exploration
Program. This work has also made use of data from the European Space Agency
(ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia
Data Processing and Analysis Consortium (DPAC;
https://www.cosmos.esa.int/web/gaia/ dpac/consortium). Funding for DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

Exoplanet Archive

## References

  * Adibekyan (2019) Adibekyan, V. 2019, Geosciences, 9, 105, doi: 10.3390/geosciences9030105
  * Akeson et al. (2013) Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989, doi: 10.1086/672273
  * An et al. (2023) An, D.-S., Xie, J.-W., Dai, Y.-Z., & Zhou, J.-L. 2023, AJ, 165, 125, doi: 10.3847/1538-3881/acb533
  * Angelo et al. (2022) Angelo, I., Naoz, S., Petigura, E., et al. 2022, AJ, 163, 227, doi: 10.3847/1538-3881/ac6094
  * Arras et al. (2012) Arras, P., Burkart, J., Quataert, E., & Weinberg, N. N. 2012, MNRAS, 422, 1761, doi: 10.1111/j.1365-2966.2012.20756.x
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bakos et al. (2004) Bakos, G., Noyes, R. W., Kovács, G., et al. 2004, PASP, 116, 266, doi: 10.1086/382735
  * Bakos et al. (2011) Bakos, G. Á., Hartman, J., Torres, G., et al. 2011, ApJ, 742, 116, doi: 10.1088/0004-637X/742/2/116
  * Banerjee et al. (2024) Banerjee, B., Narang, M., Manoj, P., et al. 2024, arXiv e-prints, arXiv:2404.16499, doi: 10.48550/arXiv.2404.16499
  * Bashi et al. (2018) Bashi, D., Helled, R., & Zucker, S. 2018, Geosciences, 8, 325, doi: 10.3390/geosciences8090325
  * Bashi et al. (2017) Bashi, D., Helled, R., Zucker, S., & Mordasini, C. 2017, A&A, 604, A83, doi: 10.1051/0004-6361/201629922
  * Bashi et al. (2023) Bashi, D., Mazeh, T., & Faigler, S. 2023, MNRAS, 522, 1184, doi: 10.1093/mnras/stad999
  * Biazzo et al. (2022) Biazzo, K., Bozza, V., Mancini, L., & Sozzetti, A. 2022, in Astrophysics and Space Science Library, Vol. 466, Demographics of Exoplanetary Systems, Lecture Notes of the 3rd Advanced School on Exoplanetary Science, ed. K. Biazzo, V. Bozza, L. Mancini, & A. Sozzetti, 143–234, doi: 10.1007/978-3-030-88124-5_3
  * Bonomo et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107, doi: 10.1051/0004-6361/201629882
  * Chabrier et al. (2009) Chabrier, G., Baraffe, I., Leconte, J., Gallardo, J., & Barman, T. 2009, in American Institute of Physics Conference Series, Vol. 1094, 15th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. E. Stempels, 102–111, doi: 10.1063/1.3099078
  * Clausen & Tilgner (2015) Clausen, N., & Tilgner, A. 2015, A&A, 584, A60, doi: 10.1051/0004-6361/201526082
  * Correia et al. (2012) Correia, A. C. M., Boué, G., & Laskar, J. 2012, ApJ, 744, L23, doi: 10.1088/2041-8205/744/2/L23
  * Correia et al. (2020) Correia, A. C. M., Bourrier, V., & Delisle, J. B. 2020, A&A, 635, A37, doi: 10.1051/0004-6361/201936967
  * Cosentino et al. (2012) Cosentino, R., Lovis, C., Pepe, F., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84461V, doi: 10.1117/12.925738
  * Cumming et al. (2008) Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531, doi: 10.1086/588487
  * Cumming et al. (1999) Cumming, A., Marcy, G. W., & Butler, R. P. 1999, ApJ, 526, 890, doi: 10.1086/308020
  * Dawson & Johnson (2012) Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122, doi: 10.1088/0004-637X/756/2/122
  * Dawson & Johnson (2018) —. 2018, ARA&A, 56, 175, doi: 10.1146/annurev-astro-081817-051853
  * Deleuil et al. (2012) Deleuil, M., Bonomo, A. S., Ferraz-Mello, S., et al. 2012, A&A, 538, A145, doi: 10.1051/0004-6361/201117681
  * Dong et al. (2021) Dong, J., Huang, C. X., Zhou, G., et al. 2021, ApJ, 920, L16, doi: 10.3847/2041-8213/ac2600
  * Fabrycky & Tremaine (2007) Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298, doi: 10.1086/521702
  * Feng et al. (2022) Feng, F., Butler, R. P., Vogt, S. S., et al. 2022, ApJS, 262, 21, doi: 10.3847/1538-4365/ac7e57
  * Ford & Rasio (2008) Ford, E. B., & Rasio, F. A. 2008, ApJ, 686, 621, doi: 10.1086/590926
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, emcee: The MCMC Hammer, Astrophysics Source Code Library.  http://ascl.net/1303.002
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Smart, R. L., Sarro, L. M., et al. 2021, A&A, 649, A6, doi: 10.1051/0004-6361/202039498
  * Ghosh & Chatterjee (2023) Ghosh, T., & Chatterjee, S. 2023, arXiv e-prints, arXiv:2304.12352, doi: 10.48550/arXiv.2304.12352
  * Giacalone et al. (2017) Giacalone, S., Matsakos, T., & Königl, A. 2017, AJ, 154, 192, doi: 10.3847/1538-3881/aa8c04
  * Goldberg & Batygin (2022) Goldberg, M., & Batygin, K. 2022, AJ, 163, 201, doi: 10.3847/1538-3881/ac5961
  * Goldreich & Soter (1966) Goldreich, P., & Soter, S. 1966, Icarus, 5, 375, doi: 10.1016/0019-1035(66)90051-0
  * Hara et al. (2019) Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738, doi: 10.1093/mnras/stz1849
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * He et al. (2020) He, M. Y., Ford, E. B., Ragozzine, D., & Carrera, D. 2020, AJ, 160, 276, doi: 10.3847/1538-3881/abba18
  * Hebb et al. (2009) Hebb, L., Collier-Cameron, A., Loeillet, B., et al. 2009, ApJ, 693, 1920, doi: 10.1088/0004-637X/693/2/1920
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Husnoo et al. (2012) Husnoo, N., Pont, F., Mazeh, T., et al. 2012, MNRAS, 422, 3151, doi: 10.1111/j.1365-2966.2012.20839.x
  * Ida & Lin (2004) Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388, doi: 10.1086/381724
  * Jackson et al. (2023) Jackson, J. M., Dawson, R. I., Quarles, B., & Dong, J. 2023, AJ, 165, 82, doi: 10.3847/1538-3881/acac86
  * Kane et al. (2012) Kane, S. R., Ciardi, D. R., Gelino, D. M., & von Braun, K. 2012, MNRAS, 425, 757, doi: 10.1111/j.1365-2966.2012.21627.x
  * Ke-ting et al. (2024) Ke-ting, S., Dong-sheng, A., Ji-wei, X., & Ji-lin, Z. 2024, Chinese Astron. Astrophys., 48, 1, doi: 10.1016/j.chinastron.2024.03.006
  * Kipping (2008) Kipping, D. M. 2008, MNRAS, 389, 1383, doi: 10.1111/j.1365-2966.2008.13658.x
  * Kipping (2013) —. 2013, MNRAS, 434, L51, doi: 10.1093/mnrasl/slt075
  * Kipping (2014) —. 2014, MNRAS, 444, 2263, doi: 10.1093/mnras/stu1561
  * Lissauer et al. (2011) Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8, doi: 10.1088/0067-0049/197/1/8
  * Lucy & Sweeney (1971) Lucy, L. B., & Sweeney, M. A. 1971, AJ, 76, 544, doi: 10.1086/111159
  * Mardling & Lin (2002) Mardling, R. A., & Lin, D. N. C. 2002, ApJ, 573, 829, doi: 10.1086/340752
  * Matsumoto et al. (2015) Matsumoto, Y., Nagasawa, M., & Ida, S. 2015, ApJ, 810, 106, doi: 10.1088/0004-637X/810/2/106
  * Matsumura et al. (2010) Matsumura, S., Peale, S. J., & Rasio, F. A. 2010, ApJ, 725, 1995, doi: 10.1088/0004-637X/725/2/1995
  * Matsumura et al. (2008) Matsumura, S., Takeda, G., & Rasio, F. A. 2008, ApJ, 686, L29, doi: 10.1086/592818
  * Mazeh et al. (2016) Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75, doi: 10.1051/0004-6361/201528065
  * Mazeh & Shaham (1979) Mazeh, T., & Shaham, J. 1979, A&A, 77, 145
  * Mills et al. (2019) Mills, S. M., Howard, A. W., Petigura, E. A., et al. 2019, AJ, 157, 198, doi: 10.3847/1538-3881/ab1009
  * Mishra et al. (2023) Mishra, L., Alibert, Y., Udry, S., & Mordasini, C. 2023, A&A, 670, A68, doi: 10.1051/0004-6361/202243751
  * Mordasini & Burn (2024) Mordasini, C., & Burn, R. 2024, arXiv e-prints, arXiv:2404.15555, doi: 10.48550/arXiv.2404.15555
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics
  * Naef et al. (2001) Naef, D., Latham, D. W., Mayor, M., et al. 2001, A&A, 375, L27, doi: 10.1051/0004-6361:20010853
  * Obertas et al. (2023) Obertas, A., Tamayo, D., & Murray, N. 2023, MNRAS, doi: 10.1093/mnras/stad1921
  * O’Toole et al. (2009) O’Toole, S. J., Tinney, C. G., Jones, H. R. A., et al. 2009, MNRAS, 392, 641, doi: 10.1111/j.1365-2966.2008.14051.x
  * Papaloizou (2011) Papaloizou, J. C. B. 2011, Celestial Mechanics and Dynamical Astronomy, 111, 83, doi: 10.1007/s10569-011-9344-4
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of machine learning research, 12, 2825
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96, doi: 10.1051/0004-6361/202038306
  * Perryman (2018) Perryman, M. 2018, The Exoplanet Handbook
  * Petrovich et al. (2014) Petrovich, C., Tremaine, S., & Rafikov, R. 2014, ApJ, 786, 101, doi: 10.1088/0004-637X/786/2/101
  * Pollacco et al. (2006) Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407, doi: 10.1086/508556
  * Pont et al. (2011) Pont, F., Husnoo, N., Mazeh, T., & Fabrycky, D. 2011, MNRAS, 414, 1278, doi: 10.1111/j.1365-2966.2011.18462.x
  * Quirrenbach et al. (2014) Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91471F, doi: 10.1117/12.2056453
  * Rasio et al. (1996) Rasio, F. A., Tout, C. A., Lubow, S. H., & Livio, M. 1996, ApJ, 470, 1187, doi: 10.1086/177941
  * Rey et al. (2018) Rey, J., Bouchy, F., Stalport, M., et al. 2018, A&A, 619, A115, doi: 10.1051/0004-6361/201833180
  * Ribas & Miralda-Escudé (2007) Ribas, I., & Miralda-Escudé, J. 2007, A&A, 464, 779, doi: 10.1051/0004-6361:20065726
  * Rosenthal et al. (2024) Rosenthal, L. J., Howard, A. W., Knutson, H. A., & Fulton, B. J. 2024, ApJS, 270, 1, doi: 10.3847/1538-4365/acffc0
  * Sagear & Ballard (2023) Sagear, S., & Ballard, S. 2023, Proceedings of the National Academy of Science, 120, e2217398120, doi: 10.1073/pnas.2217398120
  * Shabram et al. (2016) Shabram, M., Demory, B.-O., Cisewski, J., Ford, E. B., & Rogers, L. 2016, ApJ, 820, 93, doi: 10.3847/0004-637X/820/2/93
  * Socrates et al. (2012) Socrates, A., Katz, B., Dong, S., & Tremaine, S. 2012, ApJ, 750, 106, doi: 10.1088/0004-637X/750/2/106
  * Sozzetti et al. (2023) Sozzetti, A., Pinamonti, M., Damasso, M., et al. 2023, A&A, 677, L15, doi: 10.1051/0004-6361/202347329
  * Tang et al. (2024) Tang, Y., Li, X., Xiao, K., et al. 2024, Universe, 10, 182, doi: 10.3390/universe10040182
  * Tremaine (2015) Tremaine, S. 2015, ApJ, 807, 157, doi: 10.1088/0004-637X/807/2/157
  * Udry et al. (2003) Udry, S., Mayor, M., & Santos, N. C. 2003, A&A, 407, 369, doi: 10.1051/0004-6361:20030843
  * Ulmer-Moll et al. (2019) Ulmer-Moll, S., Santos, N. C., Figueira, P., Brinchmann, J., & Faria, J. P. 2019, A&A, 630, A135, doi: 10.1051/0004-6361/201936049
  * Van Eylen & Albrecht (2015) Van Eylen, V., & Albrecht, S. 2015, ApJ, 808, 126, doi: 10.1088/0004-637X/808/2/126
  * Van Eylen et al. (2016) Van Eylen, V., Winn, J. N., & Albrecht, S. 2016, ApJ, 824, 15, doi: 10.3847/0004-637X/824/1/15
  * Van Eylen et al. (2019) Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, AJ, 157, 61, doi: 10.3847/1538-3881/aaf22f
  * Vick et al. (2023) Vick, M., Su, Y., & Lai, D. 2023, ApJ, 943, L13, doi: 10.3847/2041-8213/acaea6
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2019) Wang, Y.-H., Wang, S., Hinse, T. C., et al. 2019, AJ, 157, 82, doi: 10.3847/1538-3881/aaf6b6
  * Weiss et al. (2023) Weiss, L. M., Millholland, S. C., Petigura, E. A., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 863
  * Wittenmyer et al. (2019) Wittenmyer, R. A., Bergmann, C., Horner, J., Clark, J., & Kane, S. R. 2019, MNRAS, 484, 4230, doi: 10.1093/mnras/stz236
  * Xie et al. (2016) Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Science, 113, 11431, doi: 10.1073/pnas.1604692113
  * Yang & Jin (2024) Yang, J., & Jin, L. 2024, A&A, 685, A35, doi: 10.1051/0004-6361/202348347
  * Yee et al. (2020) Yee, S. W., Winn, J. N., Knutson, H. A., et al. 2020, ApJ, 888, L5, doi: 10.3847/2041-8213/ab5c16
  * Zhu & Dong (2021) Zhu, W., & Dong, S. 2021, ARA&A, 59, 291, doi: 10.1146/annurev-astro-112420-020055