# Planetary Orbit Eccentricity Trends (POET<EMAIL_ADDRESS>The Eccentricity-Metallicity
Trend for Small Planets Revealed by the LAMOST-Gaia-Kepler Sample.

Dong-Sheng An (安东升)  School of Astronomy and Space Science, Nanjing
University, Nanjing 210023, People’s Republic of China  Key Laboratory of
Modern Astronomy and Astrophysics, Ministry of Education, Nanjing 210023,
People’s Republic of China  Ji-Wei Xie (谢基伟)  School of Astronomy and Space
Science, Nanjing University, Nanjing 210023, People’s Republic of China  Key
Laboratory of Modern Astronomy and Astrophysics, Ministry of Education,
Nanjing 210023, People’s Republic of China  Yuan-Zhe Dai (戴远哲)  School of
Astronomy and Space Science, Nanjing University, Nanjing 210023, People’s
Republic of China  Key Laboratory of Modern Astronomy and Astrophysics,
Ministry of Education, Nanjing 210023, People’s Republic of China  Ji-Lin Zhou
(周济林)  School of Astronomy and Space Science, Nanjing University, Nanjing
210023, People’s Republic of China  Key Laboratory of Modern Astronomy and
Astrophysics, Ministry of Education, Nanjing 210023, People’s Republic of
China

(Received XXX; Revised YYY; Accepted ZZZ)

###### Abstract

Orbital eccentricity is one of the basic planetary properties, whose
distribution may shed light on the history of planet formation and evolution.
Here, in a series of works on Planetary Orbit Eccentricity Trends (dubbed
POET), we study the distribution of planetary eccentricities and their
dependence on stellar/planetary properties. In this paper, the first work of
the POET series, we investigate whether and how the eccentricities of small
planets depend on stellar metallicities (e.g., [Fe/H]). Previous studies on
giant planets have found a significant correlation between planetary
eccentricities and their host metallicities. Nevertheless, whether such a
correlation exists in small planets (e.g. super-Earth and sub-Neptune) remains
unclear. Here, benefiting from the large and homogeneous LAMOST-Gaia-Kepler
sample, we characterize the eccentricity distributions of 244 (286) small
planets in single (multiple) transiting systems with the transit duration
ratio method. We confirm the eccentricity-metallicity trend that
eccentricities of single small planets increase with stellar metallicities.
Interestingly, a similar trend between eccentricity and metallicity is also
found in the radial velocity (RV) sample. We also found that the mutual
inclination of multiple transiting systems increases with metallicity, which
predicts a moderate eccentricity-metallicity rising trend. Our results of the
correlation between eccentricity (inclination) and metallicity for small
planet support the core accretion model for planet formation, and they could
be footprints of self (and/or external) excitation processes during the
history of planet formation and evolution.

eccentricity, planet, stellar metallicity

††journal: AAA

##  1 Introduction

Orbital eccentricity is one of the fundamental parameters in planetary
dynamics, which provides crucial constraints on planet formation and
evolution. Based on the fact that the solar system’s planets have small
orbital inclinations (mean $\sim 3^{\circ}$) and eccentricities (mean $\sim
0.06$), Kant and Laplace in the 18th century put forward that the solar system
formed from a nebula disk, laying the foundation for the modern theory of
planet formation.

Since the discovery of 51 Pegasi b by Mayor & Queloz (1995), the radial
velocity (RV) method has been widely used to detect exoplanets and to measure
their orbital eccentricities. In contrast to the near circular orbits of solar
system planets, exoplanets detected by RV are commonly found on eccentric
orbits (mean eccentricity $\sim 0.3$), which may imply that some violent
dynamical processes, e.g., planet-planet scattering (Chatterjee et al., 2008;
Ford & Rasio, 2008; Raymond et al., 2010) may occur in the history of
exoplanet formation and evolution. Although the RV method plays an important
role in measuring exoplanet eccentricity, it suffers from some notable biases
and degeneracies which can cause considerable systematical uncertainties of
eccentricity distributions (Shen & Turner, 2008; Anglada-Escudé et al., 2010;
Zakamska et al., 2011). Furthermore, the majority of eccentricities measured
by RV are limited for Jupiter sized giant planets.

After the launch of Kepler space telescope, the transit method began
dominating the search of exoplanets. However, it is hard to measure
eccentricities of individual planets via transit observations alone, except
for under some special circumstances, such as giant planets with large
eccentricities (Dawson & Johnson, 2012), systems with significant transit
timing variations (TTVs, e.g., Hadden & Lithwick, 2014), and planet host stars
precisely characterized by asteroseismology (Van Eylen & Albrecht, 2015; Van
Eylen et al., 2019). Still, one can robustly constrain the eccentricity
distribution of a sample of transiting planets from their transit duration
ratio distribution provided the knowledge of the densities of their host stars
(Ford et al., 2008). Using the transit duration method and with the Large Sky
Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al., 2012)
spectral characterization of the transit host stars, Xie et al. (2016) found
an eccentricity dichotomy in Kepler planets, namely, Kepler singles are on
eccentric orbits with mean eccentricity $\sim$ 0.3, whereas the multiples are
on nearly circular and coplanar orbits similar to those of the solar system
planets.

Know thy star, know thy planet. Beside the numerous exoplanet discoveries in
recent years, great progresses have also been made in characterizing the host
stars of exoplanets with the help of various surveys of spectroscopy, e.g.,
Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et
al., 2017), California-Kepler Survey (CKS; Petigura et al., 2017), GALactic
Archaeology with HERMES (GALAH; Zucker et al., 2012) and LAMOST (Cui et al.,
2012) and of astrometry, e.g., Gaia (Gaia Collaboration et al., 2018).

For example, benefiting from Gaia DR2 (Lindegren et al., 2018), Berger et al.
(2020) produced the Gaia-Kepler Stellar Properties Catalog with uncertainties
of stellar radius, stellar mass, stellar density and $\log{g}$ are
$\thicksim$4%, $\thicksim$7%, $\thicksim$13% and $\thicksim$0.05 dex
respectively. In addition, LAMOST has observed more than one third of the
Kepler targets without bias towards planet host stars (De Cat et al., 2015;
Zong et al., 2018), providing a homogeneous stellar characterization (e.g.,
$\log{g}$, RV, metallicity, etc.) in the Kepler field. Recently, by combining
the data of LAMOST, Gaia and Kepler, Chen et al. (2021) has built a catalog of
kinematic properties (i.e., Galactic positions, velocities, and the relative
membership probabilities among the thin disk, thick disk, Hercules stream, and
the halo) as well as other basic stellar parameters for 35,835 Kepler stars,
i.e., the LAMOST-Gaia-Kepler (LGK) sample, paving the road towards studying
exoplanets in different Galactic environments and ages. Thanks to the advances
both in exoplanet discovery and stellar characterization, we are allowed to
carry out a series of work from here on, a project that we dub Planetary Orbit
Eccentricity Trends (POET), aiming to reveal the dependence of orbital
eccentricity on various stellar and planetary properties, and thus to provide
new insights into planet formation and evolution.

Among various stellar properties, metallicity is a crucial factor in planet
formation and evolution. On the one hand, metallicity can influence planetary
occurrence rate. The occurrence of giant planets is strongly correlated to
stellar metallicity (e.g., $\rm Fe/H$), i.e., the so called planet-metallicity
relation, which has been well established (Fischer & Valenti, 2005) and is
considered as one of key evidences to support the core accretion theory of
planet formation (Ida & Lin, 2004). Recent studies have shown that the planet-
metallicity relation may also extend to smaller planets (Wang & Fischer, 2015;
Zhu, 2019), especially to those on short period orbits (Mulders et al., 2016;
Dong et al., 2018).

On the other hand, metallicity may also influence the orbital architecture of
a planetary system, e.g., eccentricity. Both Dawson & Murray-Clay (2013) and
Buchhave et al. (2018) found that giant planets (e.g., Jovian planets) on
eccentric orbits are preferentially residing with metal-rich stars. However,
the situation is still unclear for small planets, e.g., super-Earths and/or
sub-Neptunes. Using the asteroseismology sample, Van Eylen et al. (2019)
studied the eccentricities of planets ($R_{p}<6$ $R_{\oplus}$), and found no
significant trend with stellar metallicity. Nevertheless, using the CKS
sample, Mills et al. (2019) found tentative evidence that eccentric small
planets prefer metal-rich stars.

Benefiting from the precise stellar parameters (e.g., mass, radius, and
metallicity) provided by the LAMOST-Gaia survey (Berger et al., 2020; Chen et
al., 2021), and following the method from Xie et al. (2016), here, in the
first paper of the series of POET, we revisit the dependence of small planets’
($R_{p}<4$ $R_{\oplus}$) eccentricities on stellar metallicities with the LGK
sample (Chen et al., 2021). This paper is organized as follows. In Section 2,
we describe the sample selection procedure. In Section 3, we introduce the
method of deriving eccentricity distribution. We present our results in
Section 4 and discuss the implications of the results in Section 5. Finally,
we give a summary in Section 6.

##  2 sample selection

Figure 1: Stellar properties in our sample. The left panel shows stellar
effective temperature vs. stellar radius, and the right panel shows stellar
mass vs. metallicity, i.e., [Fe/H].

We started our sample selection from the Kepler Data Release 25 (DR25;
Thompson et al., 2018; NASA Exoplanet Archive, 2020), which has 6923 host
stars with 8054 KOIs (Kepler Objects of Interest). First, we excluded false
positives in our sample, after which 4034 confirmed/candidate planets were
left around 3069 stars. To obtain a precise and homogeneous sample of
metallicity, we then cross-matched the Kepler DR25 with LAMOST Data Release 8
(DR8). Note that we also cross-matched Kepler DR25 with LAMOST DR5 in the same
way, and removed the stars whose [Fe/H] difference between DR5 and DR8 is
greater than 3 $\sigma$. After this, 1409 planets in 1049 systems were left
with a median metallicity uncertainty of $\thicksim$ 0.04 dex for the host
stars. This uncertainty of metallicity reflects only the internal uncertainty
of LAMOST measurements (see Figure S1 of Xie et al. (2016)). For the
systematic uncertainty, Xie et al. (2016) (see their Figure S2) found there is
no significant offset but a larger dispersion, increasing the median of the
total uncertainty of [Fe/H] to $\sim$0.1 dex. In the following analyses, we
will adopt a relatively large bin size of [Fe/H] ($\sim$0.15-0.6 dex) to
reduce the effect of [Fe/H] uncertainty. Subsequently, we cross-matched the
data with Berger et al. (2020) for other stellar parameters, resulting in 1343
planets in 995 systems with median uncertainties of $\thicksim$4%,
$\thicksim$7%, $\thicksim$0.05 dex in stellar radius, stellar mass, and
$\log{g}$ respectively. To exclude the influence of potential binary stars, we
adopted a cutoff of RUWE$<1.2$.111RUWE is the abbreviation of Gaia DR2’s re-
normalized unit-weight error computed by Berger et al. (2020). A large RUWE,
e.g, RUWE$>$1.2, generally indicates the target is probably a binary. In
addition, stars with GOF222GOF is the abbreviation of goodness-of-fit computed
by Berger et al. (2020), a star with GOF $\leq$ 0.99 is an outlier with
abnormal small errors. $\leq$ 0.99 should be cautious (Berger et al., 2020)
and were excluded. We then adopted the following cuts, i.e., $\log{g}$ $>$ 4,
$4700K<T_{\rm eff}<6500K$ to focus on solar type main sequence stars. After
all the above criteria on stars, there were left 899 planets in 638 systems.

We also applied the following criteria on planets to further refine the
sample. Following Borucki et al. (2011), we adopted a transit signal noise
ratio cut SNR$>$7.1 to select reliable planet candidates. Similar to Mills et
al. (2019), here we also applied a cut on the uncertainty of the radius ratio
of planet and star, i.e., relative error of
$r\equiv\frac{R_{planet}}{R_{star}}<0.3$. Following Thompson et al. (2018), we
also selected KOIs with disposition score larger than 0.9 to have a reliable
sample of planets. Furthermore, as mentioned before the dependence of
eccentricity on metallicity for gas giant planets have been relatively well
established (Dawson & Murray-Clay, 2013; Buchhave et al., 2018) but not for
small planets. Therefore, we only focus on small planets ($R_{p}<4$
$R_{\oplus}$) here. In addition, the orbit of planet with short period would
be circularized via the tide between the host and the planet. To avoid the
influence of tide, we only considered planets with orbital period $P>5$ days (
e.g., Dong et al., 2021; Van Eylen et al., 2019). Finally, we have 244 single
transiting systems with 244 small planets and 152 multiple transiting systems
with 286 small planets in our sample. The data of the sample are provided in
the Appendix B (Table 2 and Table 3).

Table 1 is a summary of the sample selection process. Some basic properties of
stars and planets in our sample are shown in Figure 1 and Figure 2,
respectively. From the right panel of the Figure 1, we can see a trend that
stellar metallicity increases with stellar mass. Thus, to study the
relationship between eccentricity and metallicity, one should remove the
potential effects of stellar mass (see Section 4.1.2).

Figure 2: Orbital period vs. radius of planets in our sample.  Table 1:
Summary of the sample selection

| Host star | Planet
---|---|---
| (Number) | (Number)
Kepler DR25 | 6923 |  | 8054
Not false positive | 3069 |  | 4034
Match with LAMOST | 1049 |  | 1409
Match with Berger et al. (2020) | 995 |  | 1343
RUWE $\leq$ 1.2 | 862 |  | 1166
GOF $>$ 0.99 | 860 |  | 1164
Main sequence $(\log{g}>4)$ | 740 |  | 1024
4700$K$ $<$ $T_{\rm eff}$ $<$ 6500K | 638 |  | 899
SNR (signal noise ratio) $>$ 7.1 | 632 |  | 891
Relative error of $r\equiv\frac{R_{p}}{R_{*}}$ $<$ 0.3 | 618 |  | 875
Single (Multiple)a | 450 (168) |  | 450 (424)
Disposition score $>$ 0.9 | 394 (166) |  | 394 (380)
$R_{p}<4$ $R_{\oplus}$ | 345 (164) |  | 345 (362)
Period $>$ 5 days | 244 (152) |  | 244 (286)

  * 1

aThe error of transit duration, stellar radius, stellar mass are needed, thus
one planet is removed for the absent of transit duration error.

##  3 method

We follow Xie et al. (2016) to derive eccentricity of our planet sample. We
briefly summarize it as the follows. And the details of the method can be
found in the supplementary of their paper (see also in Ford et al. (2008) and
Moorhead et al. (2011)).

The basic idea is to use the distribution of transit duration ratio
($TDR\equiv T/T_{0}$) to constrain the eccentricity distribution, where $T$ is
the observed transit duration and $T_{0}$ is the reference transit duration
which assumes the transit impact parameter $b=0$ and eccentricity $e=0$. For
illustrative purpose, $TDR$ is a function of $e$, $b$, and the argument of
pericenter $\omega$,

$TDR\equiv\frac{T}{T_{0}}\thicksim\frac{\sqrt{(1-b^{2})(1-e^{2})}}{1+e\sin{\omega}}.$
(1)

Since there are three unknowns ($e$, $b$ and $\omega$) in Equation 1, one
cannot solve the eccentricity individually. Nevertheless, by assuming
reasonable distributions of $b$ and $\omega$, we can constrain the
distribution of $e$ from the distribution of a sample of observed $T/T_{0}$.

In practice, we use a more precise formula to model the transit duration which
is (Kipping, 2010),

$T_{\rm mod}\ =\
\frac{P}{\pi}\frac{(1-e^{2})^{3/2}}{(1+e\sin{\omega})^{2}}\arcsin{(\frac{\sqrt{(1+r)^{2}-b^{2}}}{b\tan{i_{0}}})},$
(2)

where $P$ and $i_{0}$ are the orbital period and the inclination of the
planet, respectively. And $i_{0}$ is related to $b$ via

$b=4.2P^{2/3}(\frac{\rho_{*}}{\rho_{\odot}})^{1/3}\frac{(1-e^{2})}{1+e\sin{\omega}}\cos{i_{0}},$
(3)

where $\rho_{*}$ is the density of the host star. $T_{0}$ is calculated as

$T_{0}\ =\
13P^{\frac{1}{3}}(\frac{\rho_{*}}{\rho_{\odot}})^{-\frac{1}{3}}(1+r).$ (4)

Then, we have the modeled transit duration ratio $TDR_{\rm mod}=T_{\rm
mod}/T_{0}$, and the observed transit duration ratio $TDR_{\rm obs}=T_{\rm
obs}/T_{0}$. The distribution of the eccentricity is obtained by fitting
$TDR_{\rm obs}$ with $TDR_{\rm mod}$.

Maximum likelihood method is used to conduct the eccentricity fitting process.
The likelihood for a given model (assuming Rayleigh distribution function for
planets with a mean eccentricity $\bar{e}$) to produce an observed TDR is
(Hadden & Lithwick, 2014)

$\mathcal{L}(TDR_{\rm obs}|\bar{e},\bar{i})=\int\,P(TDR_{\rm
mod}|\bar{e},\bar{i})*e^{\frac{-(TDR-TDR_{\rm
obs})^{2}}{2\sigma_{TDR}^{2}}}dTDR.$ (5)

where $\bar{i}$ is the mean mutual inclination (relevant only for multiple
transiting systems) of the planets. And $\sigma_{\rm TDR}$ is the uncertainty
of $TDR_{\rm obs}$, calculated by
$\sigma_{TDR}=\sqrt{(\frac{\sigma_{T_{obs}}}{T_{obs}})^{2}+(\frac{\sigma_{\rho_{*}}}{3\rho_{*}})^{2}+(\frac{\sigma_{r}}{1+r})^{2}}$,
where $\sigma_{T_{obs}}$, $\sigma_{\rho_{*}}$, and $\sigma_{r}$ are the
uncertainties of $T_{obs}$, $\rho_{*}$, and $r$ respectively. The first term,
$P(TDR_{\rm mod}|\bar{e},\bar{i})$ in Equation 5 is the probability that a
modeled transit planet produces the corresponding TDR given the mean
eccentricity $\bar{e}$ and the mean inclination $\bar{i}$. And the second term
reflects the consideration of uncertainty by assuming Gaussian error. Here, we
simply assume the distribution of $TDR_{obs}$ is symmetric about $TDR$. In
Appendix A, we test the effect of asymmetric posterior distribution of
$TDR_{obs}$.

For single systems, $\bar{i}$ is not defined and there is only one fitting
parameter, $\bar{e}$. To obtain $P(TDR_{\rm mod}|\bar{e})$ for single planets,
we conduct the following simulations using the transit planets in our sample.
First, for each planet, we assign it an orbital eccentricity $e$ from a
Rayleigh distribution with a mean of $\bar{e}$, an argument of pericenter
$\omega$ and an $\cos{i_{0}}$ from uniform distributions. We repeat this step
if the assigned eccentricity is so high that let the planet hit the surface of
the star, i.e., $a(1-e)<R_{*}+R_{p}$, where $a$ is the orbital semi-major
axis. Next, we calculate the impact parameter using Equation 3. If the
absolute value of the impact parameter is too high, i.e., $|b|>1+r$ to make a
transit, then we go back to the above first step to restart the simulation.
Then, we calculate the transit duration $T_{\rm mod}$ (Equation 2) and $T_{0}$
(Equation 4). We also estimate the modeled transit signal noise ratio
($SNR_{\rm mod}$) from the observed one ($SNR_{\rm obs}$), i.e., $SNR_{\rm
mod}=\sqrt{T_{\rm mod}/T_{\rm obs}}*SNR_{\rm obs}$. We set a simple criterion
$SNR_{\rm mod}>7.1$ to ensure the modeled transit planet to be detectable
333Note, the criterion $SNR_{\rm mod}>7.1$ adopted is rather simple. In
reality, the detection efficiency function is complicated. For example, it
increases gradually above 7.1, rather than a sharp transition as one would
compute from an idealized model.. Otherwise, we go back to the above first
step to restart the simulation. We repeat above steps until each observed
planet has 300 corresponding modeled transit duration ratio ($TDR_{\rm mod}$)
in our simulated sample. Finally, we use the Gaussian Kernel Density
Estimation function to fit the distribution of all the simulated $TDR_{\rm
mod}$ to obtain the probability density function, i.e, $P(TDR_{\rm
mod}|\bar{e})$.

For multiple systems, the method to obtain $P(TDR_{\rm mod}|\bar{e},\bar{i})$
is similar to single systems, except that the mutual inclinations in the
system are correlated. We following Zhu et al. (2018),

$\cos{i_{0}}=\cos{I}\cos{i}-\sin{I}\sin{i}\cos{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}},$
(6)

where $i_{0}$ is the inclination of the planet relative to the line of sight,
$I$ is the invariable plane for the system ($\cos{I}$ is a uniform
distribution), $\phi$ is the phase angle and drawn from a uniform distribution
independently, $i$ is the inclination of the planet in the system relative to
the invariable plane and drawn from a Rayleigh distribution with a mean
inclination $\bar{i}$.

We multiply the likelihood to produce each observed transit duration ratio to
calculate the total likelihood $\mathcal{L}(TDR_{\rm obs}|\bar{e},\bar{i})$.
For single systems, we only have one parameter, we calculate the total
likelihood as a function of $\bar{e}$ , and fit it with an polynomial
function. As shown in Figure 4, the best fit $\bar{e}$ (where
$\mathcal{L}(TDR_{\rm obs}|\bar{e})$ is the maximum) and the $1\sigma$ (68.3%)
confidence of interval are shown. For multiple systems, we map the total
likelihood in the $\bar{e}$-$\bar{i}$ plane, and give the best fit of
$\bar{e}$, $\bar{i}$ and corresponding $1\sigma$ (68.3%) confidence of
intervals (as shown in Figure 10).

##  4 results

Since the distribution of the transit duration ratio can reflect the
distribution of eccentricity (Section 3), we first plot a [Fe/H]-TDR diagram
for singles (the left panel) and multiples (the right panel) respectively in
Figure 3 to have an intuitive feeling of transit duration ratio distribution
as a function of metallicity. As can be seen, the TDR distribution for singles
is apparently wider in metal-rich systems than in metal-poor systems, which
qualitatively suggests larger planetary eccentricities associate with higher
stellar metallicities. However, the TDR distribution for multiples is more
concentrated around TDR=1 (indicating low eccentricities), and there is no
significant dependence on metallicity.

Figure 3: The metallicity ([Fe/H])- transit duration ratio (TDR) diagram for
single (left panel) and multiple (right panel) transit systems.

###  4.1 Single Transit Systems

####  4.1.1 Without Parameter control: Metal Rich Stars Host High-e Planets

Figure 4: Transit duration fitting for single transiting systems. Left panels:
the cumulative distributions for modeled (gray) and observed (red) transit
duration ratios (TDR). The range of metallicity and the number of planets for
each sub-sample are also printed. Right panels: the relative likelihood of
fitting the observed TDR distribution vs. the $\bar{e}$ assumed in the modeled
TDR distributions, which gives the $1\ \sigma,\ 2\ \sigma,$ and $3\ \sigma$
confidence interval of $\bar{e}$.  Figure 5: Cumulative distributions for
$T_{\rm eff}$, $M_{*}$, P, and $R_{p}$ in different sub-samples before (top
four panels) and after (bottom four panels) parameter control, where bin1
(dotted line), bin2 (dashed line), and bin3 (solid line) represent the lowest,
the intermediate, and the highest metallicity bin respectively. P-values of
the Kolmogorov–Smirnov (KS) test between different bins are provided, where
p-value13 and p-value23 represent the p-value of bin1 VS bin3 and bin2 VS
bin3.

In order to quantify the eccentricity distributions of different
metallicities, we then divide the single sample into three sub-samples
according to [Fe/H]. Specifically, we first sort the whole single sample in
the order of [Fe/H]. We take the $\thicksim$ 20% of systems at the highest
[Fe/H] end as the metal-rich bin. And divide the rest into two bins in
approximately equal size. The reason that we artificially make the two latter
bins larger is we can apply parameter control to further analysis (Section
4.1.2).

For each of the three sub-samples, we fit the transit duration ratio
distribution to constrain the eccentricity distribution by following the
method as say in Section 3. The results are shown in Figure 4. As can be seen,
$\bar{e}=0.09_{-0.02}^{+0.02}$ in the lowest metallicity bin,
$\bar{e}=0.15_{-0.03}^{+0.04}$ in the intermediate metallicity bin, and
$\bar{e}=0.28_{-0.05}^{+0.05}$ in the highest metallicity bin. Above results
quantitatively suggest that eccentricity increases with stellar metallicity
for small single planets.

However, these results may be influenced by other parameters. For example, the
above three sub-samples may differ in other parameters, e.g., stellar
temperature, mass, planetary period, and radius (top panels of Figure 5),
which could also affect eccentricity distribution. In the next subsection we
will perform a parameter control analysis to isolate the effect of metallicity
on eccentricity.

Figure 6: Same as Figure 4, but results for after parameter control on $T_{\rm
eff}$, $M_{*}$, P, and $R_{p}$.

####  4.1.2 With Parameter control: Minimize the Effects of other Stellar and
Planetary Properties

In order to minimize the effects of other stellar and planetary properties, we
control these parameters (stellar effective temperatures $T_{\rm eff}$,
stellar masses $M_{*}$, planets’ periods $P$, and planets’ radii $R_{p}$) to
let them have similar distributions in all the bins. Specifically, for each
system in the metal-rich bin ($\thicksim$ 20% of the whole sample) we search
for the nearest two neighbors in the metal-poor and metal-intermediate bin
respectively. We calculate the Euclidean distance (D) between planet systems
as follows to find the nearest two neighbors. Specifically,
$D=((k_{1}\frac{\Delta T_{\rm eff}}{\widetilde{\Delta T_{\rm
eff}}})^{2}+(k_{2}\frac{\Delta M_{*}}{\widetilde{\Delta
M_{*}}})^{2}+(k_{3}\frac{\Delta P}{\widetilde{\Delta
P}})^{2}+(k_{4}\frac{\Delta R_{p}}{\widetilde{\Delta R_{p}}})^{2})^{1/2}$ ,
where $\Delta T_{\rm eff}$, $\Delta M_{*}$, $\Delta P$, $\Delta R_{p}$ are the
differences in stellar effective temperature, stellar mass, planetary period
and planetary radius, and $\widetilde{\Delta T_{\rm eff}}$, $\widetilde{\Delta
M_{*}}$, $\widetilde{\Delta P}$, $\widetilde{\Delta R_{p}}$ are the
corresponding typical values for scaling purpose, which are calculated from
the following procedure. For each of $N$ systems in the metal-rich bin, we
calculate the $\Delta T_{\rm eff}$ between the system and systems in other
bins, then find the smallest two $\Delta T_{\rm eff}$. The $\widetilde{\Delta
T_{\rm eff}}$ is calculated as the median of the $2\times N$ $\Delta T_{\rm
eff}$. And $\widetilde{\Delta M_{*}}$, $\widetilde{\Delta P}$,
$\widetilde{\Delta R_{p}}$ are calculated by following the same procedure.
$k_{1}$, $k_{2}$, $k_{3}$, $k_{4}$ are four weighting coefficiencies. We tried
different $k_{1}$, $k_{2}$, $k_{3}$, $k_{4}$ (range from 0.1 to 20) and
performed the Kolmogorov–Smirnov (KS) test between the metal-rich bin and the
selected neighbor systems in other bins for $T_{\rm eff}$, $M_{*}$, $P$, and
$R_{p}$ to evaluate the goodness of finding the neighbors, and adopted the
$k_{1}$, $k_{2}$, $k_{3}$, $k_{4}$ which lead to the highest p-value of KS
test.

For example, before parameter control, we have 49, and 97 systems in metal-
rich and metal-poor bin. We select the neighbors of the metal-rich bin in the
metal-poor bin by following the steps below. Step 1: we calculate the scaling
parameters $\widetilde{\Delta T_{\rm eff}}=7.400K$, $\widetilde{\Delta
M_{*}}=0.003M_{\odot}$, $\widetilde{\Delta P}=0.204$day, $\widetilde{\Delta
R_{p}}=0.017R_{\oplus}$. Step 2: for each system in the metal-rich bin, we
select its two nearest neighbors in the metal-poor bin, i.e. the two systems
with the smallest D values given a set of $k_{1}$, $k_{2}$, $k_{3}$, $k_{4}$.
Step 3: we perform KS tests in the distributions of $T_{\rm eff}$, $M_{*}$,
$P$, and $R_{p}$ between the metal-rich bin and the selected metal-poor bin,
and record the smallest p-value ($P_{KS}$) of the four KS tests. Step 4: we
repeat step 2 and step 3 for 10,000 times by adopting different sets of
$k_{1}$, $k_{2}$, $k_{3}$, $k_{4}$ and choose the one with the highest
$P_{KS}$ as the final result. After the above steps and removing some
duplicate systems, finally, we select 38 systems in the metal-poor bin as the
neighbors of systems in the metal-rich bin.

In Figure 5, we show the distributions of these parameters before (top panels)
and after (bottom panels) the parameter control process. Before the parameter
control, bins of different metallicities differ significantly in the
distributions of $T_{\rm eff}$ and $M_{*}$. Metal-rich stars tend to have
larger masses, which also can be seen in the right panel of Figure 1. After
the parameter control, the metal-poor and metal-intermediate bins both have
similar distributions in $T_{\rm eff}$, $M_{*}$, $P$, and $R_{p}$ as compared
to the metal-rich bin (with all KS test p-values $\geq$ 0.28, indicating that
the two samples are likely to be drawn from the same distribution.). To
further quantify how well the parameters have been controlled between
different metallicity bins. We compare their median values and the
corresponding 68.3$\%$ intervals. Specifically, $T_{\rm
eff}(K)=5832_{-316}^{+282}$, $5828_{-464}^{+356}$, $5762_{-398}^{+338}$ after
the parameter control for bin1, bin2 and bin3 in Figure 5. The difference of
the median, i.e., $T_{\rm eff}\sim 70k$, between different bins is much
smaller than the 68.3$\%$ interval and is even smaller than the typical
measuring error of $\sim 100k$. Similarly,
$M_{*}(M_{\odot})=1.015_{-0.150}^{+0.132}$, $1.023_{-0.145}^{+0.202}$,
$1.056_{-0.127}^{+0.201}$ after the parameter control for bin1, bin2 and bin3
in Figure 5. The difference of the median, i.e., $M_{*}\sim 0.04M_{\odot}(\sim
4\%)$, between different bins is much smaller than the 68.3$\%$ interval and
is even smaller than the typical measuring error of $\sim 7\%$. For plantary
properties, $P($days$)=16.72_{-9.65}^{+44.23}$, $13.00_{-6.51}^{+31.86}$,
$11.52_{-5.40}^{+39.94}$ after the parameter control for bin1, bin2 and bin3.
The difference of the median, i.e., $P\sim 5.2$days, between different bins is
much smaller than the 68.3$\%$ interval.
$R_{p}(R_{\oplus})=2.268_{-0.982}^{+0.816}$, $2.076_{-0.813}^{+0.744}$,
$2.105_{-0.876}^{+0.872}$ after the parameter control for bin1, bin2 and bin3.
The difference of the median, i.e., $R_{p}\sim 0.19R_{\oplus}$, between
different bins is much smaller than the 68.3$\%$ interval. So far, $T_{\rm
eff}$, $M_{*}$, $P$, and $R_{p}$ have been well controlled for different Fe/H
bins, therefore, any significant eccentricity-Fe/H trend identified after the
above parameter control should not affect by these controlled parameters.

Then for each controlled sub-sample, we preform the transit duration ratio
distribution fitting as mentioned in Section 3 to obtain mean eccentricity.
The fitting results are shown in Figure 6. The mean eccentricities of planets
from the bin of lowest to highest in metallicity (from the top to the bottom
panels) are $\bar{e}=0.07_{-0.04}^{+0.04}$, $\bar{e}=0.13_{-0.04}^{+0.06}$,
and $\bar{e}=0.28_{-0.05}^{+0.05}$ respectively. Figure 7 compare the results
before and after the parameter control. The error bars become larger than
before due to the reduction of number of planet in each sub-sample after
parameter control. As can be seen, both show a similar trend that eccentricity
increases with metallicity. This may suggest eccentricity of small planet is
not sensitive to stellar effective temperature and stellar mass (the main
controlled parameters here).

Figure 7: Mean eccentricities (median $\bar{e}$ from Section 4.1.1, 4.1.2) vs.
[Fe/H] (median [Fe/H] of each sub-sample), where error bars represent 1
$\sigma$ uncertainties of the mean eccentricity. And the gray dashed lines are
the boundaries of each sub-sample. Boundaries between each sub-sample change
less than 0.01 after parameter control, so we don’t update the boundaries
here. In addition, the result of the highest metallicity sub-sample before
parameter control has been artificially moved to the left by 0.02 index to
avoid overlap of symbols.

####  4.1.3 Effects of Binning

Figure 8: Same as Figure 4, but results for the four sub- samples after
parameter control on $T_{eff}$, $M_{*}$, P, and $R_{p}$.

In our standard method, we set the size of the last bin as $\thicksim$ 20% of
the total sample, here we test whether our results are sensitive to the bin
size . Specifically, we consider two other cases with the last bin size as
$\thicksim$ 15% and $\thicksim$ 25%, then preformed the same parameter control
procedure and transit duration ratio fitting to derive eccentricity
distribution. For the case of $\thicksim$ 15%, we find that
$\bar{e}=0.07_{-0.04}^{+0.04}$, $\bar{e}=0.10_{-0.04}^{+0.06}$, and
$\bar{e}=0.31_{-0.05}^{+0.06}$ for the metal-poor, metal-intermediate, and
metal-rich bins respectively. For the case of $\thicksim$ 25%, we find that
$\bar{e}=0.08_{-0.04}^{+0.05}$, $\bar{e}=0.17_{-0.04}^{+0.05}$, and
$\bar{e}=0.27_{-0.04}^{+0.04}$ for the metal-poor, metal-intermediate, and
metal-rich bins respectively. As can be seen, the results of different bin
sizes are consistent with each other (within 1 $\sigma$), and all show the
same trend that planetary eccentricity increases with stellar metallicity.

We also check the influence of the number of bins. We redivide all planets in
our sample into four sub-samples, then perform the same parameter control
procedure and transit duration ratio fitting to derive eccentricity
distribution. The results are shown in Figure 8:
$\bar{e}=0.07_{-0.05}^{+0.07}$, $0.09_{-0.04}^{+0.05}$,
$0.10_{-0.05}^{+0.07}$, and $0.31_{-0.05}^{+0.06}$ for metallicity from low
(top panel) to high (bottom panel) respectively. These results also show that
eccentricity increases with stellar metallicity, which is still consistent
with the result in Section 4.1.2 where we divide the sample into three sub-
samples.

Therefore, we conclude that our results are not sensitive to the choice of bin
size nor bin number.

Figure 9: Left panel: the fits for the three sub-samples, where the turquoise
dashed line represents the exponential best fit. $\Delta AIC$ is AIC
difference between the constant and exponential best fit. The light turquoise
band is 1 $\sigma$ confidence interval for the exponential best fit. Right
panel: same as the left panel, but shows the fits for the four sub-samples.

####  4.1.4 Fit the Metallicity-Eccentricity Relation

To quantitatively study the relationship between metallicity and eccentricity,
here we fit our results by considering three models: a constant model
($\rm{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bar{e}}=constant$),
a linear model
($\rm{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bar{e}}=a*[Fe/H]+b$),
and an exponential model
($\rm{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bar{e}}=a*10^{b*[Fe/H]}$)
using least square method. In order to evaluate these models, we adopt the
Akaike Information Criterion (AIC) (Akaike, 1974). Generally, a model with
smaller AIC score is statistically better. We calculate the AIC scores of best
fits and corresponding parameters for the three models. For our standard case
of three bins, $\rm AIC=14.5$, $\rm AIC=4.6$, and $\rm AIC=4.0$ for the
constant, linear, and exponential models respectively. The exponential model
is preferred with the smallest AIC.

In order to investigate the effect of the data uncertainty on the fitting
results, we performed the following re-sampling and fitting analysis. We re-
sample $\bar{e}$ according to the fitted probability distribution function of
$\bar{e}$ (black curves in the  right panels of Figure 6) and re-fit the
eccentricity-metallicity relation with the above three models. We repeat the
re-sample and re-fit procedure 1,000,000 times and record the best fit
parameters and the corresponding AIC in each time. For our standard case of
three bins, we find that compared to the constant model, the linear model is
preferred with smaller AIC in 966,932 times, and the exponential model is
preferred in 978,835 times, corresponding to confidence levels of 96.7% and
97.9%.

Since the exponential model is preferred from the AIC analysis, we adopt it as
our nominal model. The best fit parameters are obtained by using the median
$\bar{e}$ of each bins (Section 4.1.2 and Section 4.1.3) with the least square
method. The 1 $\sigma$ confidence intervals of the model parameters are taken
as $50\pm 34.1$ percentile of the 1,000,000 times re-sample fitting results.
The result for three sub-samples is

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(7)

And the result for four sub-samples is

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(8)

These results are shown in Figure 9. As can be seen they are consistent within
1 $\sigma$.

To summary this section, we find the trend that the eccentricity increases
with metallicity is robust and it is best fit with an exponential model.

###  4.2 Multiple Transit Systems

####  4.2.1 Metallicity-Eccentricity Trend

Similar to Section 4.1.1 and Section 4.1.2, we divide the multiple planets’
hosts into three sub-samples and preformed the same parameter control
procedure. For each sub-sample of the multiples, we simulate the transit
duration to constrain the mean eccentricity $\bar{e}$ and the mean inclination
$\bar{i}$ via the method as mentioned in Section 3. The results are shown in
Figure 10. As can be seen that the best fit of $\bar{e}$ are 0, 0 and 0.05 in
the lowest metallicity bin, the intermediate metallicity bin and the highest
metallicity bin respectively. If we use the median of $\bar{e}$, then
$\bar{e}=0.026_{-0.026}^{+0.045}$, $0.030_{-0.030}^{+0.031}$ and
$0.048_{-0.048}^{+0.063}$ from the bin of lowest to highest metallicity. We
also fit the metallicity-e relation with the constant, linear and exponential
model, and AIC=2.1, 4.0, 4.0 correspondingly. The constant model is preferred
with the smallest AIC. Therefore, although the median eccentricity of
multiples tend to increase slightly with metallicity, the trend is barely
significant due to the relatively large uncertainties. As for inclination, the
constrain is weak with much larger error bar compared to the eccentricity. We
will revisit the metallicity-inclination trend in Section 4.2.2.

Figure 11 shows the metallicity-eccentricity trend for singles and multiples.
While the singles shows a strong rising trend between eccentricity and
metallicity, the trend is weaker for multiples given the relatively large
uncertainties. For all the metallicity bins, the multiples have smaller
eccentricities than singles and the difference in eccentricity is larger for
higher metallicity.

Figure 10: Transit duration fitting for multiple transiting systems. Left
panels: the cumulative distributions for modeled and observed (red) transit
duration ratios. For comparison, simulations of $\bar{e}=0-0.4$ and
$\bar{i}=10^{\circ}$ are shown in gray, simulations of $\bar{e}=0$ and
$\bar{i}=0-10^{\circ}$ are shown in yellow to blue. Right panel: contours of
relative likelihood in $\bar{e}$-$\bar{i}$ plane. The best fit of $\bar{e}$,
$\bar{i}$ (also shown as black star) and corresponding confidence interval
($1\ \sigma,\ 2\ \sigma,$ and $3\ \sigma$) are also printed. Figure 11:
Similar to Figure 7, but showing mean eccentricities of singles vs. multiples
under different metallicities. The solid error bars represent $\bar{e}$
constrained via TDR (Section 4.1 and 4.2.1), while the dashed error bars
represent $\bar{e}$ calculated by assuming $\bar{e}=\bar{i}$ (see Section 5.2)
and $\bar{i}$ is constrained by the distribution of mutual transit duration
ratio, i.e., $\xi$ (Section 4.2.2). In addition, the dashed error bar of the
medium metallicity has been artificially moved to right by 0.02 dex to avoid
overlap of symbols.

####  4.2.2 Metallicity-Inclination Trend

Figure 12:  The mutual transit duration ratio ($\xi$) fitting for multiple
transiting systems. Left panels: the observed $\xi$ distribution (red) and the
best fit (turquoi). Right panels: the relative likelihood of fitting the
observed $\xi$ distribution as a function of the $\bar{i}$ and $\bar{e}$
assumed in modeling the $\xi$ distributions, which gives the $1\ \sigma$
confidence interval of $\bar{i}$ and $\bar{e}$. The metallicty increases from
the top to the bottom panels, i.e., [F/H]$<$-0.02, -0.02$\leq$[F/H]$<$0.16 and
[F/H]$\geq$0.16 respectively.

In addition to the normalized transit duration ratio (TDR, Equation 5), there
is a another metric, i.e., the mutual transit duration ratio ($\xi$), which is
more sensitive to mutual inclination of multiple transiting systems. Following
Fabrycky et al. (2014), $\xi$ is defined as

${\xi\equiv\frac{T_{in}/P_{in}^{1/3}}{T_{out}/P_{out}^{1/3}}},$ (9)

where $T$ and $P$ are the transit duration and orbital period, “in” and “out”
represents the inner and the outer planets. Similar to Fabrycky et al. (2014)
and Zhu et al. (2018), here we use $\xi$ to constrain the mean inclinations of
the multiples (the same three sub-samples as in Section 4.2.1)444After the
parameter control, some multiple systems just left one small planet, and thus
these systems were excluded in the $\xi$ fitting process.. Following Zhu et
al. (2018), the likelihood of the simulated $\xi$ produce the observed
$\xi_{j}$ under the mean inclination $\bar{i}$ and the mean eccentricity
$\bar{e}$ is defined as,

$\mathcal{L}=\prod_{j=1}^{N}\int\,P_{sim}(\ln\xi)*e^{\frac{-(\ln\xi_{j}-\ln\xi)^{2}}{2\sigma_{\ln\xi,j}^{2}}}d\ln\xi,$
(10)

where $P_{sim}(\ln\xi)$ is the probability that the model produces the
corresponding $\xi$ (in log scale) given the mean inclination $\bar{i}$ and
the mean eccentricity $\bar{e}$, $\xi_{j}$ is the observed $\xi$ for the $j$th
planet pair, and $\sigma_{\ln\xi,j}$ is the corresponding uncertainty. Similar
to modeling TDR as mentioned in section 3, the modeling of $\ln\xi$ and thus
the calculation of $P_{sim}(\ln\xi)$ also take into account the transit
geometry effect and the detection efficiency by simply adopting a signal noise
ratio SNRmod cut at 7.1. The above likelihood function implicitly assumes that
the $\ln\xi$ follows the Gaussian distribution. In order to test if a Gaussian
approximation can apply to $\ln\xi$, we randomly draw 100,000 $T_{in}$ and
$T_{out}$ from Gaussian distributions given their reported values and errors,
and take the value of $P_{in}$ and $P_{out}$ without errors (a good
approximation) to calculate 100,000 $\ln\xi$ using Equation 9. We find the
distribution of $\ln\xi$ can be well fit by a Gaussian distribution. Figure 12
and Figure 13 show the $\xi$ fitting results, where the mean mutual
inclinations are constrained to be $\bar{i}=0.6_{-0.5}^{+0.3}$,
$3.3_{-1.7}^{+1.3}$ and $5.5_{-1.1}^{+2.3}$ degrees from the bin of lowest
(the top panel) to highest (the bottom panel) metallicity. Although with large
uncertainties, the mean mutual inclination tends to increase with metallicity.

Figure 13:  Mean mutual inclinations $\bar{i}$ of multiple transiting systems
constrained by $\xi$ distribution (Section 4.2.2) as a function of [Fe/H].
Each point in the error bar represents the median of [Fe/H] and $\bar{i}$ in
each sub-sample

##  5 discussions

Figure 14: Planet eccentricity (e) as a function of the stellar metallicity
([Fe/H]) for single (left panel) and multiple (right panel) systems in the RV
sample. For the RV singles, the data prefer an exponential model and best fit
is denoted by the gray dashed lines. For comparison, we also plot the best fit
for Kepler singles (Section 4.1.4, the turquoise dashed lines). The colored
bands represent 1 $\sigma$ confidence interval of the exponential best fit.
For the RV multiples, the data prefer a flat model, $e=0.12$ (gray dashed
line). For comparison, we also plot the best fit of Kepler multiples
($\bar{e}=0.03$, the yellow dashed line).

###  5.1 Comparison to radial velocity planets

Our above results are based on the Kepler sample, here we compare them with
that of the radial velocity (RV) sample. Specifically, we download the RV
sample from the NASA Exoplanet Archive (NASA Exoplanet Archive, 2021). To
build a RV sample that is comparable to the Kepler sample, we select RV
planets with the following criteria. First, we retain planets that have both
reported eccentricities and host metallicities. Then, we set conditions:
$M\sin{i}<32M_{\oplus}(\thicksim 0.1M_{J})$ or $R_{p}<4R_{\oplus}$ to focus on
small planets. Subsequently, we adopt an eccentricity quality cut, i.e.,
relative error of eccentricity $<75\%$ and absolute error of eccentricity
$<0.1$. And we also tried different relative error and absolute error cut,
e.g., 100% and 0.2, and all the results are similar. In addition, similar to
the Kepler sample, we only consider the systems with one (more than one)
planet (planets) in $P<400$ days to build the RV single (multiple) sample. We
also exclude planets with $P<5$ days to reduce the influence of tide. Finally,
we have 18 (29) planets in the selected RV single (multiple) sample. The data
of these 18 (29) RV singles (multiples) are provided in the Appendix B (Table
4 and Table 5).

Figure 14 shows the distribution of the RV sample in the eccentricity-
metallicity diagram. Apparently for singles in the left panel of Figure 14,
there is also a trend that planetary eccentricity increases with stellar
metallicity. In order to qualify this trend, we preformed the same AIC
analysis as in the Kepler sample (Section 4.1.4). We find that $\rm
AIC=106.2$, $\rm AIC=67.0$, and $\rm AIC=66.1$ for the constant, linear, and
exponential models respectively. The exponential model is still preferred with
the smallest AIC and the best fit (as print in Figure 14) is consistent with
the result of the Kepler singles (Section 4.1.4) within 1 $\sigma$.

###  5.2 Comparison to previous studies

Van Eylen et al. (2019) studied eccentricities of small planets using the
asteroseismology sample and find no significant trend between single planetary
eccentricity and stellar metallicity. Nevertheless, we note that their sample
is small (only 30 planets with radius less than 6 $R_{\oplus}$) and the
stellar metallicities are not homogeneous but coming from different
literature.

Mills et al. (2019) also studied the eccentricity of small planets by using
the data from California-Kepler Survey (CKS). They identified 7 single planets
(6 with radius less than 4 $R_{\oplus}$) with high eccentricity and all of
them are orbiting around metal-rich stars ([Fe/H]$>$0). Therefore, they
tentatively conclude that small eccentric planets are preferentially found in
high metallicity stars. In this study, with a large and homogeneous sample
from the LAMOST-Gaia-Kepler catalog (Chen et al., 2021), we confirm the trend
that eccentricity increases with stellar metallicity for single small planets.
Encouragingly, a similar eccentricity-metallicity relation for singles is also
revealed by the RV sample (Figure 14).

Xie et al. (2016) shows that Kepler multiples and solar system objects follow
a common relation ($\bar{e}\sim(1$–$2)\times\bar{i}$) between mean
eccentricities and mutual inclinations. Given such a correlation between
eccentricity and inclination and the observed rising trend between inclination
and metallcitity for multiple systems (Fig.13), one may expect that the
eccentricity of multiples should also increase with metallicity. However, our
results show that, in multiple systems, the rising trend with metallicity is
much weaker for eccentricity than for inclination (Figure 11 vs Figure 13).
One possible reason could be the much lower precision in measuring
eccentricity than in measuring inclination for transiting systems. If we adopt
a simple correlation $\bar{e}=\bar{i}$ and use the inclination-metallicity
trend in Figure 13 to produce an expected eccentricity-metallicity trend
(dashed bars in Figure 11), we find that current eccentricity measurements
(solid bars in Figure 11) actually do not rule out such an expected trend
considering the relatively large uncertainties.

Another major finding of Xie et al. (2016) is that Kepler singles are on
eccentric orbits with $\bar{e}\sim 0.3$, whereas the multiples are on nearly
circular($\bar{e}\sim 0.04$). In this work, our measurement of $\bar{e}$ for
singles is generally lower than that of Xie et al. (2016). This is probably
because we exclude large planets (presumably have larger eccentricities than
small planets) to focus on small planets here. Nevertheless, in line with Xie
et al. (2016), our results also show that singles have larger eccentricities
than multiples (Figure 11). Furthermore, we find that the difference in
eccentricity between singles and multiples increases with metallicity. We will
further discuss the implications of this result below.

###  5.3 Implications to planet formation and evolution

The observed result that eccentricity increases with metallicity for small
planet is not unexpected, instead it has important implications to planet
formation and evolution. According to the core accretion model for planet
formation, planets form in proto-planetary disk through a bottom-up process
from dust, planetesimals, planetary embryos and finally to full planets (Ida &
Lin, 2004). Generally, a higher stellar metallicity suggests a higher disk
metallicity and thus more solid materials to form planetesimals and more
massive planets, which have stronger gravitational interactions to pump up
larger orbital eccentricities. Specifically, we expect two eccentricity
exciting mechanisms as follows.

On the one hand, eccentricities can be self excited among small planets
themselves. According to the N-body simulations by Moriarty & Ballard (2016),
the mean eccentricities of the planets increase from $\bar{e}\thicksim 0.06$
to $\thicksim 0.10$ when the total mass of the planetesimals in the disk
increases from 7 $M_{\oplus}$ to 35 $M_{\oplus}$. If the increase of the mass
is caused by the increase of the metallicity completely, we can estimate the
corresponding metallicities ($Z$) by $Z=-\lg(0.01M_{g}/M_{s})$, where $M_{s}$
and $M_{g}$ are the mass of solid and gas in the disk respectively (Greaves et
al., 2007). Here, we assume the total mass of the gas plus solid remain as a
constant and mass of the disk to be 0.01 $M_{\odot}$ (Hayashi (1981), for the
solar-like system) to estimate metallicity, e.g., $M_{s}+M_{g}=0.01M_{\odot}$.
Therefore, the mass of solid increases from 7 $M_{\oplus}$ to 35 $M_{\oplus}$
corresponds to metallicity increases from [Fe/H]$\thicksim-0.68$ to
[Fe/H]$\thicksim 0.03$. Such an increase in metallicity leads to an
eccentricity increase from $\thicksim 0.01$ to $\thicksim 0.12$ according
Equation 7. This result is comparable to the N-body simulation result
($\bar{e}$ from $\thicksim 0.06$ to $\thicksim 0.10$) except at the metal-
poorest end where Equation 7 underestimates the eccentricity somewhat. Such a
small inconsistency at the metal-poorest end is not unexpected because, in
reality, a small eccentricity of a few percent is generally common for an
essentially circular orbit. Furthermore, if taking an approximation that
mutual inclination is correlated to eccentricity ($\bar{e}\sim
1-2\times\bar{i}$, Xie et al. (2016)), then the inclination can be pumped up
to $\bar{i}\sim 0.05-0.1\sim 2.9-5.7^{\circ}$ by self-excitation, which is
comparable to our finding for Kepler multiples (the middle point of Figure
13). Based on the above simple estimation, we conclude that the self-
excitation mechanism should, more or less, have played a role in producing the
observed eccentricity (inclination)-metallicity trends.

On the other hand, eccentricities of inner small planets can be excited and
with the transiting multiplicity being reduced (producing more singles) at the
same time, by perturbations of outer giant planets (e.g. Huang et al. (2017);
Pu & Lai (2018); Poon & Nelson (2020)). Combining this with the well known
metallicity-giant planet correlation (Fischer & Valenti, 2005), one may
naturally expect that a correlation between eccentricity and metallicity for
small planets in single transiting systems. Therefore, the outer-perturbed
mechanism is complementary to the self-excitation mechanism, the combination
of which is promising to fully explain the eccentricity-metallicity trend of
singles. Unfortunately, our sample of planets are detected by Kepler with the
transit method, which strongly biases against long period planets, thus we
could not test this mechanism in this work. In the near future, the Gaia
mission would find many long period giant planets, which combined with
Kepler’s short period small planets can further explore this scenario. In
addition, RV surveys of planets found by TESS, K2 and PLATO could also help
test the predictions of this scenario.

##  6 summary

Since the discovery of 51 Pegasi b (Mayor & Queloz, 1995), the number of
exoplanets has increased dramatically. Furthermore, various surveys of
spectroscopy and astrometry provide comprehensive characterizations for the
host stars of exoplanets, allowing one to statistically study the relationship
between stars and planets. Here we start a project, Planetary Orbit
Eccentricity Trends (POET) to investigate how orbital eccentricities of
planets depend on various stellar/planetary properties. In this work, the
first paper of the POET series, we study the relationship between small
planets’ ($R_{p}<4$ $R_{\oplus}$) eccentricities and stellar metallicities
with the LAMOST-Gaia-Kepler sample (Chen et al., 2021).

We found that, in single transiting systems, eccentricities of small planets
increase with stellar metallicities (Figure 4). We excluded the influences of
$T_{\rm eff}$, $M_{*}$, $P$, and $R_{p}$ on the eccentricity-metallicity trend
by adopting a parameter control method to control these parameters (Figure 6,
7). We also explored the effects of binning and found the eccentricity-
metallicity trend is not sensitive to the size nor number of bins (Section
4.1.3 and Figure 8). Furthermore, we fitted the eccentricity-metallicity trend
and found it is best fitted with an exponential function (Section 4.1.4 and
Equation 7).

In contrast, we found that, in multiple transiting systems, the eccentricity-
metallicity rising trend is less clear. Although an inclination-metallcity
trend is seen in multiples (Figure 13) and it predicts a moderate
eccentricity-metallcity rising trend as well, such an eccentricity-metallcity
rising trend can neither be established nor be ruled out given the relatively
large uncertainty in measuring eccentricity (Figure 11).

We then compared our results with the data from RV, and found they are
consistent within 1 $\sigma$ (Figure 14). We also compared our results with
Van Eylen et al. (2019) and Mills et al. (2019) in Section 5.2, where we
emphasized the importance to have large and homogeneous stellar parameters
when studying the relation between stars and planets. Our results have shown
that the difference of the mean eccentricity between singles and multiples
increases with metallicity. Finally, we discussed the implication of the
eccentricity-metallicity trend on planet formation and evolution (Section
5.3). We identified two mechanisms (self-excitation and external excitation)
that could potentially explain the observed eccentricity-metallicity trend.
Future studies of both simulations and observations on a larger sample will
further test them.

The authors thank the referee for the constructive comments and suggestions.
The authors also thank Fei Dai, Jia-Yi Yang, Di-Chang Chen, Rui-Sheng Zhang,
Ke-Ting Shin, and Tong Bao for helpful discussions. This work is supported by
the National Key R$\&$D Program of China (No. 2019YFA0405100) and the National
Natural Science Foundation of China (NSFC; grant Nos. 12150009, 11933001, and
12273011). J.-W.X. also acknowledges the support from the National Youth
Talent Support Program and the Distinguish Youth Foundation of Jiangsu
Scientific Committee (BK20190005). This work has included data from
Guoshoujing Telescope. Guoshoujing Telescope (the Large Sky Area Multi-Object
Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project
built by the Chinese Academy of Sciences. Funding for the project has been
provided by the National Development and Reform Commission. LAMOST is operated
and managed by the National Astronomical Observatories, Chinese Academy of
Sciences.

##  Appendix A The effect of the asymmetric distribution of $TDR_{obs}$

In the main text of this paper, we have assumed the distribution of the
$TDR_{obs}$ is symmetric about $TDR$ in the likelihood function (Equation 5).
However, the percentage of stars below the zero age main sequence (ZAMS) is
small for Kepler stars, which will lead to an asymmetric distribution of
$TDR_{obs}$. Specifically, to test the effect of the asymmetry of stellar
density, we modify the uncertainty term in likelihood function as,

$\sigma_{TDR}=\sqrt{(\frac{\sigma_{T_{obs}}}{T_{obs}})^{2}+K(\frac{\sigma_{\rho_{*}}}{3\rho_{*}})^{2}+(\frac{\sigma_{r}}{1+r})^{2}}\begin{cases}K=0&\text{
$TDR<TDR_{obs}$ }\\\ K=1&\text{ $else$ }\end{cases}$ (A1)

where $K$ is a coefficient to take into account the asymmetry of stellar
density. Here, we consider an extreme case, i.e., $K=0$ for $TDR<TDR_{obs}$
and $K=1$ for $TDR\geq TDR_{obs}$.  This is equivalent to adopting a mixture
of two half-Normal distributions with different widths for each observed
$TDR_{obs}$.  Figure 15 shows the results for the 3 sub-samples (same sub-
samples as in Section 4.1.2) by using the above modified uncertainty term. As
can be seen, the mean eccentricities are $\bar{e}=0.08_{-0.04}^{+0.04}$,
$\bar{e}=0.12_{-0.04}^{+0.05}$, and $\bar{e}=0.28_{-0.04}^{+0.02}$ from the
lowest [Fe/H] bin to the highest [Fe/H] bin respectively. These results are
very close to our norminal results in Section 4.1.2, where we assumed a
symmetric $TDR_{obs}$. Therefore, we conclude that the asymmetric distribution
of $TDR_{obs}$ has little effect on our results.

##  Appendix B Data of the planet systems used in this paper

We provide the data used in this work here. Table 2 and Table 3 show a part of
the Kepler singles and Kepler multiples analysed in this work respectively.
And Table 4 and Table 5 show a part of the RV singles and RV multiples
analysed in this work respectively.

Table 2: Kepler single planet systems analysed in this work.

Kepler Name | period | $R/R_{*}$ | transit duration | SNR | …
---|---|---|---|---|---
| (day) |  | (hour) |  | …
K00049.01 | $8.313773_{-0.000042}^{+0.000042}$ | $0.02636_{-0.00081}^{+0.00124}$ | $2.964_{-0.125}^{+0.125}$ | 19.6 |
K00072.02 | $45.294223_{-0.000056}^{+0.000056}$ | $0.01981_{-0.00013}^{+0.00033}$ | $6.830_{-0.034}^{+0.034}$ | 134.8 |
K00084.01 | $9.286964_{-0.000004}^{+0.000004}$ | $0.02378_{-0.00011}^{+0.00035}$ | $3.365_{-0.015}^{+0.015}$ | 251.9 |

  * •

Note: Planet parameters, stellar parameters (except [Fe/H]) and [Fe/H] are
from Kepler DR 25 (NASA Exoplanet Archive, 2020), Berger et al. (2020) and
LAMOST DR8 respectively. Here, only a part of table is shown and the whole
table is available online.

Table 3: Kepler multiple planet systems analysed in this work.

Kepler Name | period | $R/R_{*}$ | transit duration | SNR | …
---|---|---|---|---|---
| (day) |  | (hour) |  | …
K00041.01 | $12.815904_{-0.000019}^{+0.000019}$ | $0.01407_{-0.00049}^{+0.00159}$ | $6.386_{-0.044}^{+0.044}$ | 106.7 |
K00041.02 | $6.887071_{-0.000023}^{+0.000023}$ | $0.00808_{-0.00063}^{+0.00001}$ | $4.468_{-0.076}^{+0.076}$ | 38.2 |
K00041.03 | $35.333193_{-0.000225}^{+0.000225}$ | $0.00962_{-0.00054}^{+0.00198}$ | $5.966_{-0.128}^{+0.128}$ | 28.8 |

  * •

Note: Planet parameters, stellar parameters (except [Fe/H]) and [Fe/H] are
from Kepler DR 25 (NASA Exoplanet Archive, 2020), Berger et al. (2020) and
LAMOST DR8 respectively. Here, only a part of table is shown and the whole
table is available online.

Table 4: RV single planet systems analysed in this work.

Planet Name | Radius | Mass(or Mass*$\sin i$) | e | Metallicity | Period
---|---|---|---|---|---
| (R⊕) | (M⊕) |  |  | (day)
HD 97658 b | $2.12_{-0.06}^{+0.06}$ | $8.30_{-1.10}^{+1.10}$ | $0.050_{-0.030}^{+0.040}$ | $-0.230_{-0.030}^{+0.030}$ | $9.4897116_{-0.0000008}^{+0.0000008}$
HD 95338 b | $3.89_{-0.20}^{+0.19}$ | $42.44_{-2.08}^{+2.22}$ | $0.197_{-0.024}^{+0.029}$ | $0.040_{-0.100}^{+0.100}$ | $55.0870000_{-0.0200000}^{+0.0200000}$

  * •

Note: Planet and stellar parameters are from NASA Exoplanet Archive (2021).
Here, only a part of table is shown and the whole table is available online.

Table 5: RV multiple planet systems analysed in this work.

Planet Name | Radius | Mass(or Mass*$\sin i$) | e | Metallicity | Period
---|---|---|---|---|---
| (R⊕) | (M⊕) |  |  | (day)
HD 219134 f | $1.31_{-0.02}^{+0.02}$ | $7.30_{-0.40}^{+0.40}$ | $0.148_{-0.047}^{+0.047}$ | $0.110_{-0.040}^{+0.040}$ | $22.7170000_{-0.0150000}^{+0.0150000}$
HD 219134 c | $1.51_{-0.05}^{+0.05}$ | $4.36_{-0.22}^{+0.22}$ | $0.062_{-0.039}^{+0.039}$ | $0.110_{-0.040}^{+0.040}$ | $6.7645800_{-0.0003300}^{+0.0003300}$

  * •

Note: Planet and stellar parameters are from NASA Exoplanet Archive (2021).
Here, only a part of table is shown and the whole table is available online.

Figure 15: Same as Figure 6, but results for the consideration of asymmetry
density.

## References

  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Anglada-Escudé et al. (2010) Anglada-Escudé, G., López-Morales, M., & Chambers, J. E. 2010, ApJ, 709, 168, doi: 10.1088/0004-637X/709/1/168
  * Berger et al. (2020) Berger, T. A., Huber, D., van Saders, J. L., et al. 2020, AJ, 159, 280, doi: 10.3847/1538-3881/159/6/280
  * Borucki et al. (2011) Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19, doi: 10.1088/0004-637X/736/1/19
  * Buchhave et al. (2018) Buchhave, L. A., Bitsch, B., Johansen, A., et al. 2018, ApJ, 856, 37, doi: 10.3847/1538-4357/aaafca
  * Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580, doi: 10.1086/590227
  * Chen et al. (2021) Chen, D.-C., Yang, J.-Y., Xie, J.-W., et al. 2021, AJ, 162, 100, doi: 10.3847/1538-3881/ac0f08
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197, doi: 10.1088/1674-4527/12/9/003
  * Dawson & Johnson (2012) Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122, doi: 10.1088/0004-637X/756/2/122
  * Dawson & Murray-Clay (2013) Dawson, R. I., & Murray-Clay, R. A. 2013, ApJ, 767, L24, doi: 10.1088/2041-8205/767/2/L24
  * De Cat et al. (2015) De Cat, P., Fu, J. N., Ren, A. B., et al. 2015, ApJS, 220, 19, doi: 10.1088/0067-0049/220/1/19
  * Dong et al. (2021) Dong, J., Huang, C. X., Dawson, R. I., et al. 2021, ApJS, 255, 6, doi: 10.3847/1538-4365/abf73c
  * Dong et al. (2018) Dong, S., Xie, J.-W., Zhou, J.-L., Zheng, Z., & Luo, A. 2018, Proceedings of the National Academy of Science, 115, 266, doi: 10.1073/pnas.1711406115
  * Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146, doi: 10.1088/0004-637X/790/2/146
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102, doi: 10.1086/428383
  * Ford et al. (2008) Ford, E. B., Quinn, S. N., & Veras, D. 2008, ApJ, 678, 1407, doi: 10.1086/587046
  * Ford & Rasio (2008) Ford, E. B., & Rasio, F. A. 2008, ApJ, 686, 621, doi: 10.1086/590926
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Greaves et al. (2007) Greaves, J. S., Fischer, D. A., Wyatt, M. C., Beichman, C. A., & Bryden, G. 2007, MNRAS, 378, L1, doi: 10.1111/j.1745-3933.2007.00315.x
  * Hadden & Lithwick (2014) Hadden, S., & Lithwick, Y. 2014, ApJ, 787, 80, doi: 10.1088/0004-637X/787/1/80
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Huang et al. (2017) Huang, C. X., Petrovich, C., & Deibert, E. 2017, AJ, 153, 210, doi: 10.3847/1538-3881/aa67fb
  * Ida & Lin (2004) Ida, S., & Lin, D. N. C. 2004, ApJ, 616, 567, doi: 10.1086/424830
  * Kipping (2010) Kipping, D. M. 2010, MNRAS, 407, 301, doi: 10.1111/j.1365-2966.2010.16894.x
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2, doi: 10.1051/0004-6361/201832727
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94, doi: 10.3847/1538-3881/aa784d
  * Mayor & Queloz (1995) Mayor, M., & Queloz, D. 1995, Nature, 378, 355, doi: 10.1038/378355a0
  * 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
  * Moorhead et al. (2011) Moorhead, A. V., Ford, E. B., Morehead, R. C., et al. 2011, ApJS, 197, 1, doi: 10.1088/0067-0049/197/1/1
  * Moriarty & Ballard (2016) Moriarty, J., & Ballard, S. 2016, ApJ, 832, 34, doi: 10.3847/0004-637X/832/1/34
  * Mulders et al. (2016) Mulders, G. D., Pascucci, I., Apai, D., Frasca, A., & Molenda-Żakowicz, J. 2016, AJ, 152, 187, doi: 10.3847/0004-6256/152/6/187
  * NASA Exoplanet Archive (2020) NASA Exoplanet Archive. 2020, Kepler Objects of Interest DR25, Version: 2020-Jun-4 10:37, NExScI-Caltech/IPAC, doi: 10.26133/NEA5
  * NASA Exoplanet Archive (2021) —. 2021, Planetary Systems, Version: 2021-Aug-25 19:35, NExScI-Caltech/IPAC, doi: 10.26133/NEA12
  * Petigura et al. (2017) Petigura, E. A., Howard, A. W., Marcy, G. W., et al. 2017, AJ, 154, 107, doi: 10.3847/1538-3881/aa80de
  * Poon & Nelson (2020) Poon, S. T. S., & Nelson, R. P. 2020, MNRAS, 498, 5166, doi: 10.1093/mnras/staa2755
  * Pu & Lai (2018) Pu, B., & Lai, D. 2018, MNRAS, 478, 197, doi: 10.1093/mnras/sty1098
  * Raymond et al. (2010) Raymond, S. N., Armitage, P. J., & Gorelick, N. 2010, ApJ, 711, 772, doi: 10.1088/0004-637X/711/2/772
  * Shen & Turner (2008) Shen, Y., & Turner, E. L. 2008, ApJ, 685, 553, doi: 10.1086/590548
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38, doi: 10.3847/1538-4365/aab4f9
  * 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. (2019) Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, AJ, 157, 61, doi: 10.3847/1538-3881/aaf22f
  * Wang & Fischer (2015) Wang, J., & Fischer, D. A. 2015, AJ, 149, 14, doi: 10.1088/0004-6256/149/1/14
  * 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
  * Zakamska et al. (2011) Zakamska, N. L., Pan, M., & Ford, E. B. 2011, MNRAS, 410, 1895, doi: 10.1111/j.1365-2966.2010.17570.x
  * Zhu (2019) Zhu, W. 2019, ApJ, 873, 8, doi: 10.3847/1538-4357/ab0205
  * Zhu et al. (2018) Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101, doi: 10.3847/1538-4357/aac6d5
  * Zong et al. (2018) Zong, W., Fu, J.-N., De Cat, P., et al. 2018, ApJS, 238, 30, doi: 10.3847/1538-4365/aadf81
  * Zucker et al. (2012) Zucker, D. B., de Silva, G., Freeman, K., Bland-Hawthorn, J., & Hermes Team. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 458, Galactic Archaeology: Near-Field Cosmology and the Formation of the Milky Way, ed. W. Aoki, M. Ishigaki, T. Suda, T. Tsujimoto, & N. Arimoto, 421