Planetary Orbit Eccentricity Trends (POET). I. The Eccentricity-Metallicity Trend for Small Planets Revealed by the LAMOST-Gaia-Kepler Sample

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.


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 ∼ 3 • ) and eccentricities (mean ∼ 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 ∼ 0.3), which may imply that some violent dynamical processes, e.g., planet-planet scattering (Chatterjee et al. 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.  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 σ. After this, 1409 planets in 1049 systems were left with a median metallicity uncertainty of ∼ 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 ∼0.1 dex. In the following analyses, we will adopt a relatively large bin size of [Fe/H] (∼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 ∼4%, ∼7%, ∼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. 1 In addition, stars with GOF 2 ≤ 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 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 ≡ R planet Rstar < 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 ⊕ ) 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).

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  and Moorhead et al. (2011)).
The basic idea is to use the distribution of transit duration ratio (T DR ≡ 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, T DR is a function of e, b, and the argument of pericenter ω, Since there are three unknowns (e, b and ω) in Equation 1, one cannot solve the eccentricity individually. Nevertheless, by assuming reasonable distributions of b and ω, 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), where P and i 0 are the orbital period and the inclination of the planet, respectively. And i 0 is related to where ρ * is the density of the host star. T 0 is calculated as Then, we have the modeled transit duration ratio T DR mod = T mod /T 0 , and the observed transit duration ratio T DR obs = T obs /T 0 . The distribution of the eccentricity is obtained by fitting T DR obs with T DR 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ē) to produce an observed TDR is (Hadden & Lithwick 2014) whereī is the mean mutual inclination (relevant only for multiple transiting systems) of the planets. And σ TDR is the uncertainty of T DR obs , calculated by σ T DR = ( where σ T obs , σ ρ * , and σ r are the uncertainties of T obs , ρ * , and r respectively. The first term, P (T DR mod |ē,ī) in Equation 5 is the probability that a modeled transit planet produces the corresponding TDR given the mean eccentricityē and the mean inclinationī. And the second term reflects the consideration of uncertainty by assuming Gaussian error. Here, we simply assume the distribution of T DR obs is symmetric about T DR. In Appendix A, we test the effect of asymmetric posterior distribution of T DR obs .
For single systems,ī is not defined and there is only one fitting parameter,ē. To obtain P (T DR mod |ē) 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ē, an argument of pericenter ω 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 mod (Equation 2) and T 0 (Equation 4). We also estimate the modeled transit signal noise ratio (SN R mod ) from the observed one (SN R obs ), i.e., SN R mod = T mod /T obs * SN R obs . We set a simple criterion SN R mod > 7.1 to ensure the modeled transit planet to be detectable 3 . 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 (T DR mod ) in our simulated sample. Finally, we use the Gaussian Kernel Density Estimation function to fit the distribution of all the simulated T DR mod to obtain the probability density function, i.e, P (T DR mod |ē).
For multiple systems, the method to obtain P (T DR mod |ē,ī) is similar to single systems, except that the mutual inclinations in the system are correlated. We following Zhu et al. (2018), 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), φ 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ī.
We multiply the likelihood to produce each observed transit duration ratio to calculate the total likelihood L(T DR obs |ē,ī). For single systems, we only have one parameter, we calculate the total likelihood as a function ofē , and fit it with an polynomial function. As shown in Figure 4, the best fitē (where L(T DR obs |ē) is the maximum) and the 1σ (68.3%) confidence of interval are shown. For multiple systems, we map the total likelihood in thē e-ī plane, and give the best fit ofē,ī and corresponding 1σ (68.3%) confidence of intervals (as shown in Figure 10).

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.

Single Transit Systems
3 Note, the criterion SN R 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.  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ē assumed in the modeled TDR distributions, which gives the 1 σ, 2 σ, and 3 σ confidence interval ofē.

Without Parameter control: Metal Rich Stars Host High-e Planets
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 ∼ 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,ē = 0.09 +0.02 −0.02 in the lowest metallicity bin,ē = 0.15 +0.04 −0.03 in the intermediate metallicity bin, andē = 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.

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 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 (∼ 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 ) 2 ) 1/2 , where ∆T eff , ∆M * , ∆P , ∆R p are the differences in stellar effective temperature, stellar mass, planetary period and planetary radius, and ∆T eff , ∆M * , ∆P , ∆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 ∆T eff between the system and systems in other bins, then find the smallest two ∆T eff . The ∆T eff is calculated as the median of the 2 × N ∆T eff . And ∆M * , ∆P , ∆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 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 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 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 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 eff , M * , P , and R p as compared to the metal-rich bin (with all KS test p-values ≥ 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 eff (K) = 5832 +282 −316 , 5828 +356 −464 , 5762 +338 −398 after the parameter control for bin1, bin2 and bin3 in Figure 5. The difference of the median, i.e., T eff ∼ 70k, between different bins is much smaller than the 68.3% interval and is even smaller than the typical measuring error of ∼ 100k. Similarly, M * (M ) = 1.015 +0.132 −0.150 , 1.023 +0.202 −0.145 , 1.056 +0.201 −0.127 after the parameter control for bin1, bin2 and bin3 in Figure 5. The difference of the median, i.e., M * ∼ 0.04M (∼ 4%), between different bins is much smaller than the 68.3% interval and is even smaller than the typical measuring error of ∼ 7%. For plantary properties, P (days) = 16.72 +44.23 −9.65 , 13.00 +31.86 −6.51 , 11.52 +39.94 −5.40 after the parameter control for bin1, bin2 and bin3. The difference of the median, i.e., P ∼ 5.2days, between different bins is much smaller than the 68.3% interval.
−0.876 after the parameter control for bin1, bin2 and bin3. The difference of the median, i.e., R p ∼ 0.19R ⊕ , between different bins is much smaller than the 68.3% interval. So far, T 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ē = 0.07 +0.04 −0.04 ,ē = 0.13 +0.06 −0.04 , andē = 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).

Effects of Binning
In our standard method, we set the size of the last bin as ∼ 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 ∼ 15% and ∼ 25%, then preformed the same parameter control procedure and transit duration ratio fitting to derive eccentricity distribution. For the case of ∼ 15%, we find thatē = 0.07 +0.04 −0.04 ,ē = 0.10 +0.06 −0.04 , andē = 0.31 +0.06 −0.05 for the metal-poor, metal-intermediate, and metal-rich bins respectively. For the case of ∼ 25%, we find thatē = 0.08 +0.05 −0.04 ,ē = 0.17 +0.05 −0.04 , andē = 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 σ), 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:ē = 0.07 +0.07 −0.05 , 0.09 +0.05 −0.04 , 0.10 +0.07 −0.05 , and 0.31 +0.06 −0.05 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.

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 (ē = constant), a linear model (ē = a * [Fe/H] + b), and an exponential model (ē = 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, AIC = 14.5, 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. AIC = 4.6, and 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ē according to the fitted probability distribution function ofē (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ē of each bins (Section 4.1.2 and Section 4.1.3) with the least square method. The 1 σ confidence intervals of the model parameters are taken as 50 ± 34.1 percentile of the 1,000,000 times re-sample fitting results. The result for three sub-samples is These results are shown in Figure 9. As can be seen they are consistent within 1 σ.
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.  ∆AIC is AIC difference between the constant and exponential best fit. The light turquoise band is 1 σ confidence interval for the exponential best fit. Right panel: same as the left panel, but shows the fits for the four sub-samples.
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ē and the mean inclinationī via the method as mentioned in Section 3. The results are shown in Figure 10. As can be seen that the best fit ofē 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ē, thenē = 0.026 +0.045 −0.026 , 0.030 +0.031 −0.030 and 0.048 +0.063 −0.048 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.

Metallicity-Inclination Trend
In addition to the normalized transit duration ratio (TDR, Equation 5), there is a another metric, i.e., the mutual transit duration ratio (ξ), which is more sensitive to mutual inclination of multiple transiting systems. Following Fabrycky et al. (2014), ξ is defined as 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 ξ to constrain the mean inclinations of the multiples The mutual transit duration ratio (ξ) fitting for multiple transiting systems. Left panels: the observed ξ distribution (red) and the best fit (turquoi). Right panels: the relative likelihood of fitting the observed ξ distribution as a function of theī andē assumed in modeling the ξ distributions, which gives the 1 σ confidence interval ofī andē. The metallicty increases from the top to the bottom panels, i.e., [ (the same three sub-samples as in Section 4.2.1) 4 . Following Zhu et al. (2018), the likelihood of the simulated ξ produce the observed ξ j under the mean inclinationī and the mean eccentricityē is defined as, where P sim (ln ξ) is the probability that the model produces the corresponding ξ (in log scale) given the mean inclination i and the mean eccentricityē, ξ j is the observed ξ for the jth planet pair, and σ ln ξ,j is the corresponding uncertainty. Similar to modeling TDR as mentioned in section 3, the modeling of ln ξ and thus the calculation of P sim (ln ξ) also take into account the transit geometry effect and the detection efficiency by simply adopting a signal noise ratio SNR mod cut at 7.1. The above likelihood function implicitly assumes that the ln ξ follows the Gaussian distribution. In order to test if a Gaussian approximation can apply to ln ξ, 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 ξ using Equation 9. We find the distribution of ln ξ can be well fit by a Gaussian distribution. Figure 12 and Figure 13 show the ξ fitting results, where the mean mutual inclinations are constrained to beī = 0.6 +0.3 −0.5 , 3.3 +1.3 −1.7 and 5.5 +2.3 −1.1 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.

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: 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 σ 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 (ē = 0.03, the yellow dashed line).
M sin i < 32M ⊕ (∼ 0.1M J ) or R p < 4R ⊕ 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 AIC = 106.2, AIC = 67.0, and 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 σ.

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 ⊕ ) 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 ⊕ ) 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 (ē ∼ (1-2) ×ī) 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ē =ī 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ē ∼ 0.3, whereas the multiples are on nearly circular(ē ∼ 0.04). In this work, our measurement ofē 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.

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ē ∼ 0.06 to ∼ 0.10 when the total mass of the planetesimals in the disk increases from 7 M ⊕ to 35 M ⊕ . 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 (Hayashi (1981), for the solar-like system) to estimate metallicity, e.g., M s + M g = 0.01M . Therefore, the mass of solid increases from 7 M ⊕ to 35 M ⊕ corresponds to metallicity increases from [Fe/H]∼ −0.68 to [Fe/H]∼ 0.03. Such an increase in metallicity leads to an eccentricity increase from ∼ 0.01 to ∼ 0.12 according Equation 7. This result is comparable to the N-body simulation result (ē from ∼ 0.06 to ∼ 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 (ē ∼ 1−2×ī, Xie et al. (2016)), then the inclination can be pumped up toī ∼ 0.05 − 0.1 ∼ 2.9 − 5.7 • 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.

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 ⊕ ) 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 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 eccentricitymetallcity 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 σ ( 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.

A. THE EFFECT OF THE ASYMMETRIC DISTRIBUTION OF T DR OBS
In the main text of this paper, we have assumed the distribution of the T DR obs is symmetric about T DR 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 T DR obs . Specifically, to test the effect of the asymmetry of stellar density, we modify the uncertainty term in likelihood function as, 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 T DR < T DR obs and K = 1 for T DR ≥ T DR obs . This is equivalent to adopting a mixture of two half-Normal distributions with different widths for each observed T DR 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ē = 0.08 +0.04 −0.04 ,ē = 0.12 +0.05 −0.04 , andē = 0.28 +0.02 −0.04 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 T DR obs . Therefore, we conclude that the asymmetric distribution of T DR obs has little effect on our results.

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.