# Planets Across Space and Time (PAST) IV: The Occurrence and Architecture of
Kepler Planetary Systems as a Function of Kinematic Age Revealed by the
LAMOST-Gaia-Kepler Sample

Jia-Yi Yang  School of Astronomy and Space Science, Nanjing University,
Nanjing 210023, China  Key Laboratory of Modern Astronomy and Astrophysics,
Ministry of Education, Nanjing 210023, China  Di-Chang Chen  School of
Astronomy and Space Science, Nanjing University, Nanjing 210023, China  Key
Laboratory of Modern Astronomy and Astrophysics, Ministry of Education,
Nanjing 210023, China  Ji-Wei Xie  School of Astronomy and Space Science,
Nanjing University, Nanjing 210023, China  Key Laboratory of Modern Astronomy
and Astrophysics, Ministry of Education, Nanjing 210023, China  Ji-Lin Zhou
School of Astronomy and Space Science, Nanjing University, Nanjing 210023,
China  Key Laboratory of Modern Astronomy and Astrophysics, Ministry of
Education, Nanjing 210023, China  Subo Dong  Kavli Institute for Astronomy and
Astrophysics, Peking University, Beijing 100871, People’s Republic of China
Zi Zhu  School of Astronomy and Space Science, Nanjing University, Nanjing
210023, China  Key Laboratory of Modern Astronomy and Astrophysics, Ministry
of Education, Nanjing 210023, China  Zheng Zheng  Department of Physics and
Astronomy, University of Utah, Salt Lake City, UT 84112, USA  Chao Liu
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, People’s Republic of China  Weikai Zong  Department of Astronomy,
Beijing Normal University, Beijing 100875, China  Institute for Frontiers in
Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, China
Ali Luo  National Astronomical Observatories, Chinese Academy of Sciences,
Beijing 100012, People’s Republic of China

###### Abstract

One of the fundamental questions in astronomy is how planetary systems form
and evolve. Measuring the planetary occurrence and architecture as a function
of time directly addresses this question. In the fourth paper of the Planets
Across Space and Time (PAST) series, we investigate the occurrence and
architecture of Kepler planetary systems as a function of kinematic age by
using the LAMOST-Gaia-Kepler sample. To isolate the age effect, other stellar
properties (e.g., metallicity) have been controlled. We find the following
results. (1) The fraction of stars with Kepler-like planets ($F_{\text{Kep}}$)
is about 50% for all stars; no significant trend is found between
$F_{\text{Kep}}$ and age. (2) The average planet multiplicity ($\bar{N}_{p}$)
exhibits a decreasing trend ( $\sim$2$\sigma$ significance) with age. It
decreases from $\bar{N}_{p}$ $\sim$ 3 for stars younger than 1 Gyr to
$\bar{N}_{p}$ $\sim$1.8 for stars about 8 Gyr. (3) The number of planets per
star ($\eta$=$F_{\text{Kep}}$ $\times$$\bar{N}_{p}$) also shows a decreasing
trend ($\sim$2–3$\sigma$ significance). It decreases from $\eta\sim$ 1.6–1.7
for young stars to $\eta\sim$ 1.0 for old stars. (4) The mutual orbital
inclination of the planets ($\sigma_{i,k}$) increases from $1\fdg
2^{+1.4}_{-0.5}$ to $3\fdg 5^{+8.1}_{-2.3}$ as stars aging from 0.5 to 8 Gyr
with a best fit of
$\log{\sigma_{i,k}}=0.2+0.4\times\log{\frac{\text{Age}}{\text{1Gyr}}}$.
Interestingly, the Solar System also fits such a trend. The nearly
independence of $F_{\text{Kep}}$ $\sim 50\%$ on age implies that planet
formation is robust and stable across the Galaxy history. The age dependence
of $\bar{N}_{p}$ and $\sigma_{i,k}$ demonstrates planetary architecture is
evolving, and planetary systems generally become dynamically hotter with fewer
planets as they age.

methods: statistical – planetary systems – planet–star interactions

††software: emcee (Foreman-Mackey et al., 2013), scikit-leart (Pedregosa et
al., 2012), KeplerPORTs (Burke & Catanzarite, 2017), numpy (Harris et al.,
2020), matplotlib (Hunter, 2007), astropy (Astropy Collaboration et al., 2013,
2018, 2022), evolstate (Huber, 2019), RGCCA (Girka et al., 2023)

tablenum SIXtablenum

##  1 Introduction

Thanks to various surveys from ground (e.g., Mayor et al., 2011) and space
(e.g., Borucki et al., 2010), the number of known planets has reached a
milestone (5000, NASA Exoplanet
Archive111https://exoplanetarchive.ipac.caltech.edu/, Akeson. et al., 2013).
Such a rich planetary database has enabled substantial statistical studies of
the occurrence rate and architecture of planetary systems (see the reviews by
Winn & Fabrycky, 2015; Zhu & Dong, 2021), deepening our understanding of
planet formation and evolution.

Stellar properties (e.g., mass, effective temperature, and metallicity) play
crucial roles in determining planetary occurrence rate and architecture.
Although the occurrence of giant planets (Jupiter-like gas giants) has been
found to increase with stellar mass (Johnson et al., 2010; Ghezzi et al.,
2018), the trend is opposite for small planets. For the bulk of planets
detected by the Kepler mission (so called super-Earths and sub-Neptunes, with
radii between the Earth and Neptune, hereafter dubbed as Kepler planets for
short), their occurrence rate in terms of number of planets per star ($\eta$)
has an inverse relationship with stellar temperature and mass (Howard et al.,
2012; Mulders et al., 2015; Kunimoto & Matthews, 2020). In fact, $\eta$ can be
further decomposed into two factors: the fraction of stars that have planetary
systems ($F$) and the average number of planets in a planetary system
(planetary multiplicity $\bar{N}_{p}$), and they are linked by the following
equation:

${}\eta=F\times\bar{N}_{p}.$ (1)

Further studies have shown that both $F$ and $\bar{N}_{p}$ tend to decrease as
stellar temperature and mass increase (Yang et al., 2020; He et al., 2021).

Metallicity also plays a differential role in shaping planetary systems of
giant planets and small planets. On one hand, a correlation between giant
planets and metallicity has been well established (Santos et al., 2001;
Fischer & Valenti, 2005), which provides the key evidence for the core-
accretion model of planet formation (e.g., Pollack et al., 1996; Ida & Lin,
2004). On the other hand, such a planet-metallicity correlation is generally
weaker and more complicated for small planets (Buchhave et al., 2012; Wang &
Fischer, 2015; Dong et al., 2018; Petigura et al., 2018; Zhu, 2019)

Stellar environments (e.g., stellar companions, clusters and memberships of
the Galactic thin/thick disks) also affect planetary occurrence and
architecture. There has been substantial evidence showing that planetary
occurrence is reduced and planetary architecture is modified when stellar
companions are close, with separations $\lesssim$ 100 AU (Wang et al., 2014;
Kraus et al., 2016; Moe & Kratter, 2021; Fontanive et al., 2019; Su et al.,
2021). Recently, it has been reported that both the period and radius
distributions of exoplanets exhibit dependencies on stellar clustering (Winter
et al., 2020; Kruijssen et al., 2020; Chevance et al., 2021; Longmore et al.,
2021). Dai. et al. (2021) found that stellar groups with high relative
velocities tend to have a lower occurrence rate of super-Earths and sub-
Neptunes but a higher occurrence rate of sub-Earths. The Galactic membership
and total velocity of stars are also linked with the planet occurrence rate.
It has been found that stars in the thick disk (higher total velocity)
generally have fewer planets than those in the thin disk (lower total
velocity, Bashi & Zucker, 2019, 2022; Chen et al., 2021a).

The occurrence and architecture of planets in our Galaxy could also evolve
with time. Therefore, measuring planet occurrence and architecture as a
function of time can provide crucial insights into planet formation and
evolution. For example, recent studies (e.g., Berger et al., 2020a; Sandoval
et al., 2021; David et al., 2021; Chen et al., 2022) have revealed that the
relative occurrence (ratio) of super-Earths and sub-Neptunes evolves on a time
scale of Giga years, providing crucial constraints on the formation of the
radius valley (a deficit of planets with radii of $\sim$1.7-2.1 $R_{\oplus}$,
Fulton et al., 2017). More recently, Bashi & Zucker (2022) found tentative
evidence that suggests the occurrence rate of close in super-Earths detected
by Kepler is anti-correlated with stellar age. However, such an anti-
correlation is still inconclusive, probably because they adopted the isochrone
ages which suffer from large uncertainties (56% for Kepler stars, Berger et
al., 2020b). Furthermore, they didn’t isolate the effect of age from other
stellar properties (e.g., metallicity), so it is still unclear whether the
anti-correlation is intrinsic or just a projection of other correlations.

To investigate how planet occurrence and architecture evolve with time, we
have started a series of work named Planets Across Space and Time (PAST, Chen
et al., 2021b). The first challenge of this work is to determine the age of
main-sequence stars, which make up the bulk of planet hosts. In PAST I (Chen
et al., 2021b), we revisited the kinematic method to classify Galactic
components and the Age-Velocity dispersion Relationship (AVR, Strömberg, 1946;
Wielen, 1977; Holmberg et al., 2009a), extending the viable range to 1.5 kpc
to cover the majority of the known exoplanet hosts. The deduced kinematic age
for an ensemble of stars has a typical internal uncertainty of 10%–20%. The
second challenge is to isolate the effect of stellar age, because age is
generally correlated with other properties, such as stellar mass and
metallicity. Applying the revised kinematic method of PAST I, we constructed a
catalog of kinematic and other basic properties for 35,835 Kepler stars in
PAST II (Chen et al., 2021a). Such a large and homogeneous sample enables us
to further set control samples to isolate the effect of age from other stellar
properties. In PAST III, we investigated how the radius distribution of small
planets evolves with time (Chen et al., 2022). In this work, the fourth paper
of the PAST series (PAST IV), we study the occurrence and architecture of
Kepler planets as a function of stellar age.

This paper is organized as follows: in Section 2, we present the star and
planet data used in this work. In Section 3, we describe the parameter control
method to isolate the age effect, and present the apparent occurrence of
Kepler planets. In Section 4, we adopt a forward modeling method to derive the
intrinsic occurrence rate and architecture of Kepler planets. We make the
discussions and summarize the main conclusions in Section 5 and 6.

##  2 Data sample

###  2.1 Star Sample

The LAMOST-Gaia-Kepler catalog constructed in PAST II is based on LAMOST
DR4/DR5 and GAIA DR2. Since LAMOST have updated to DR8 (Yan et al., 2022), and
Gaia has released EDR3/DR3 (Gaia Collaboration et al., 2022), therefore we
have updated the LAMOST-Gaia-Kepler catalog accordingly. We start from the
stellar properties catalogue from Berger et al. (2020b), which provides a
homogeneous calibration of effective temperature, mass, and radius for most of
the Kepler stars. The Kepler team calculated Combined Differential Photometric
Precision (CDPP, Christiansen et al., 2012; Kepler Mission, 2019; Kepler
Project, 2020) for each target, which defines the completeness of transit
searching. We restrict our sample to targets with $\sigma_{\text{CDPP}}$ value
222https://exoplanetarchive.ipac.caltech.edu/docs/Kepler_completeness_reliability.html.
Then we cross-match the sample with the recently released LAMOST DR8
333http://www.lamost.org/dr8/ low-resolution catalogue (Yan et al., 2022),
which contains more spectral observations of Kepler stars reprocessed with the
latest pipeline version. We select stars with LAMOST metallicity and radial
velocity measurements, and remove stars with [Fe/H] less than $-1.0$ due to a
lack of training set (Xiang et al., 2019). Next, we cross-match with Gaia
DR3(Gaia Collaboration et al., 2022; IRSA, 2022) catalog, which includes more
accurate measurements of positions, proper motions, and parallaxes of stars
compared to DR2. Gaia DR3 also provides the renormalised unit weight error
(RUWE, Lindegren, 2018) for identifying possible binary stars. Stars with RUWE
values greater than 1.2 are excluded from our sample (Berger et al., 2020b;
Bryson et al., 2020). We obtain 70,239 star in total, the number of our sample
after each selection step is summarized in Table 1.

Utilizing Gaia DR3, we update the kinematic method and AVR of PAST I, the
details can be seen in the Appendix A. In PAST I, the calibrations of the
kinematic method and AVR extended from the solar neighborhood to a larger
range of stars with $|Z|<1.5$ kpc, $7.5<R<10$ kpc, and distance$<$1.5 kpc,
where $Z$ and $R$ are the vertical and radial components of Galactocentric
cylindrical coordinates, respectively. Here, we adopt a similar range of
stars, but further extend the distance to 2 kpc, thanks to the improvement of
the astrometric measurements from Gaia DR2 to Gaia DR3. With the updated
kinematic method (see Appendix A), we calculate the probabilities of stars
belonging to each Galactic component, i.e, thin disk, thick disk, halo, and
Hercules stream (dub as $D$, $TD$, $H$, and $Her$). We classify stars into
different components following the commonly used method introduced by Bensby
et al. (2003, 2014), and show the results in the Toomre diagram (Figure 1).
Since AVR can only be applied to disk stars, so we limit our sample to stars
within the Galactic disk, i.e., stars with $D/Her\geqslant 2$, $TD/H\geqslant
1$, and $TD/Her\geqslant 2$.

Figure 1: Toomre diagram of stars in the updated LAMOST-Gaia-Kepler sample.
The blue, orange, grey, green, red, and grey dots present stars in the thin
disk, thick disk, in between thin and thick disk, halo, and Hercules stream,
respectively. The grey dot lines represent the total Galactic velocity at 100,
200, and 300 $\unit{km.s^{-1}}$.

In Figure 2, we show stars from the updated LAMOST-Gaia-Kepler sample in the
Hertzsprung–Russell diagram. The effective temperature and radius data are
obtained from Berger et al. (2020b), and the evolve stage is calculated by the
same method as Bryson et al. (2020) using the python package
evolstate444http://ascl.net/1905.003. We further limit our star sample to
main-sequence solar type stars, with effective temperature between 4700 to
6500 $\unit{}$. So far, the number of our star sample is 19,537.

Figure 2: Hertzsprung–Russell diagram of stars in the updated LAMOST-Gaia-
Kepler sample. The blue, orange, and red dots show stars in the main-sequence
stage, subgiant, and red giant, respectively.

###  2.2 Planet Sample

Our Kepler planet sample is based on Kepler DR25 (Thompson et al., 2018; NASA
Exoplanet Archive, 2020). We select planets/candidates within our star sample,
and exclude those flagged as ‘false positives’. We only consider planets with
periods less than 100 days, as the observed planet numbers and detection
efficiencies both drop significantly beyond this period (Burke & Catanzarite,
2017). We also exclude planetary systems with Ultra Short Period planets
(USPs, period $<$ 1 day) from our Fiducial analysis (see Section 5.1 for more
discussions on USPs). This exclusions is due to the standard Kepler pipeline
is not well conditioned to search for USPs (Sanchis-Ojeda et al., 2014), and
USPs are relatively rare (with an occurrence rate $\sim$ 1%) and may have
undergone different formation and evolution process (Dai et al., 2018). The
planet radii are derived from stellar radii (Berger et al., 2020b) and the
planet-to-star radius ratio ($R_{p}/R_{s}$, Thompson et al., 2018). Planets
with radii smaller than 0.5$R_{\oplus}$ are excluded due to their relatively
low detection efficiency (Burke & Catanzarite, 2017). Since we focus on the
occurrence rate and architecture of small planets, we exclude planet systems
with planets larger than 6$R_{\oplus}$ (see Section 5.1 for more discussions
on giant planets). The selection process of the planet sample is summarized in
Table 1. The period-radius distribution of our planet sample is shown in
Figure 3.

Figure 3: Planet sample in the period-radius diagram. The grey dots show the
whole planet sample before data selection. The blue, orange, and green dots
present planets in single, double, and three or more planet systems after we
apply all the selection criteria.

###  2.3 Kinematic Age of Planet Host and non-Host

In PAST II, we have shown that stars with different numbers of transiting
planets exhibit a drop in their kinematic age drops as a function of planet
multiplicity. Here, we also bin stars according to their transiting planet
number into four groups; in each group, stars have zero, one, two, and three
or more transiting planets, respectively. We derive the kinematic age for each
group using the updated AVR (see details in Appendix A). Since kinematic age
is calculated from the dispersion of total Galactic velocity, it can be skewed
by velocity outliers. To reduce the effect caused by outliers, we calculate
the median value and the Mean Absolute Deviation (MAD) of the total Galactic
velocity for each group. Then we remove stars with total velocities higher or
lower than $\text{Median}\pm 5\times$MAD within this group. The kinematic age
results for each group are presented in Figure 4, and we compare them to the
results of PAST II (Chen et al., 2021a). As can be seen, the kinematic ages
derived in this work are generally consistent with those of PAST II in $1\sim
2\sigma$ range. They both show a declining trend in kinematic age with
increasing planet multiplicity. Nevertheless, the ages obtained in this work
are systematically lower by about 0.2–0.5 Gyr compared to those in PAST II.
This difference is expected, because we remove outlier stars in this work,
which usually have high velocities. The removal of outliers causes a decrease
in velocity dispersion and, consequently, a lower value of kinematic age.

In total, we obtain 19,358 stars in the the star sample, and 663 planets in
467 systems. The size of our sample after each step of selection can be seen
in Table 1.

Table 1: Data selection  | Star | Planet
---|---|---
Berger et al. (2020b) | 186,301 | 3,826
With $\sigma_{\text{CDPP}}$ | 185,161 | 3,826
Match with LAMOST DR8 | 70,251 | 1,562
With RV data and $\text{[Fe/H]}\geqslant-1$ | 68,567 | 1,549
Math with Gaia DR3 | 67,922 | 1,535
RUWE$\leqslant$1.2 (Remove potential binary) | 55,332 | 1,320
$|Z|<$1.5 kpc, 7.5$<R<$10 kpc, |  |
and distance$<$2 kpc | 45,914 | 1,279
TD/H $\geqslant$1, TD/Her $\geqslant$2, D/Her $\geqslant$2 | 40,347 | 1,109
(In the thin, thick disk, or in between) |  |
Main sequence | 27,213 | 940
4700$\unit{}\leqslant$$T_{\text{eff}}$ $\leqslant 6500\unit{}$ | 19,537 | 784
Orbit Period$\leqslant$100 days | … | 720
Remove ultra short period system | … | 703
(No planet with period $<$1 day) |  |
Planet radii$\geqslant$0.5 RE | … | 698
Remove giant planet system | … | 663
(No planet with radii $>$6 RE) |  |
Median$-5\times$MAD$\leqslant V_{tot}$ |  |
and $V_{tot}\leqslant$Median$+5\times$MAD | 19,358 | 641
Figure 4: The blue dots and errorbars show the kinematic age and $\pm 1\sigma$
ranges for stars with zero, one, two, and three or more planets for this work,
and the orange dots and errorbars are values from PAST II (Chen et al.,
2021a), the x-axis is offset by 0.1 for clearance.

##  3 Apparent trend analysis from parameter control

In this section, we derive the apparent planet occurrence rate as a function
of stellar kinematic age. The apparent planet occurrence rate is defined as
the _observed_ planet multiplicity function (number of stars that have one,
two, and three or more transit planets, dubbed as $N_{1}$, $N_{2}$, and
$N_{3+}$) divided by the number of stars ($N_{star}$) in each bin.

###  3.1 Parameter Control

To isolate the effect caused by stellar age, we use the parameter control
method to reduce the influences induced by other stellar properties. In this
work, we control five properties: effective temperature, mass, metallicity,
stellar radius, and $\sigma_{\text{CDPP}}$. The former three parameters need
to be controlled because they are found to affect the intrinsic planet
occurrence rate (e.g., Buchhave et al., 2012; Yang et al., 2020; He et al.,
2021). The latter two also need to be controlled because they directly affect
the detection efficiency of transiting planets.

The basic idea of the parameter control is to let stars in different age bins
have similar distributions in the controlled stellar properties. To achieve
this goal, we apply a ‘finding star neighbors’ method, similar to the one
described in Chen et al. (2022), which involves the following steps:

  1. (1)

Grouping stars. For the whole star sample with a size of $N$, we sort the
stars according to TD/D ascendingly, which is an effective indicator of
stellar age (Chen et al., 2021b). Then we group the stars into an odd number
of bins. To implement parameter control, the middle bin contains fewer stars,
while the other bins have more stars. The farther away from the middle bin,
the more stars there are. In this study, we first consider a case of three
bins, each bin containing 40%, 20%, and 40% of the stars, to have a
qualitative view of the age trend. To further quantify the age–occurrence rate
trend, we consider a case of five bins, each containing 25%, 20%, 10%, 20%,
and 25% of the stars, respectively. Due to the limited sample size, we do not
consider cases with more bins.

  2. (2)

Choosing a standard sample. We dub the stars in the middle bin as the
‘standard sample’, and the number of stars in the central bin is denoted as
$N_{st}$.

  3. (3)

Finding the nearest neighbor stars. In each bin (except the middle one), we
select $N_{st}$ stars that are the closest neighbors of the standard sample in
the space of the controlled parameters. This is done by adopting the nearest
neighborhood method from the scikit-learn (Pedregosa et al., 2012) package.

  4. (4)

Checking parameter control result. We calculate the differences of the 25, 50,
and 75 percentiles of each controlled parameter for every two bins. If all the
differences are less than the typical errors, we consider these parameters
have been controlled. The typical errors of temperature, mass, and radius are
112, 7%, and 4%, respectively (Berger et al., 2020b). For metallicity we
choose 0.05 dex as the typical error, which is the median value of the
internal measurement uncertainties in our star sample. The Kepler team has
reported $\sigma_{\text{CDPP}}$ for different timescales. We choose
$\sigma_{\text{CDPP}}$ of 3.5 hours because it is the closest to the median
transit duration of our planet sample (3.35 hours). Since the SNR of Kepler
planet is in proportion to $(R_{p}/R_{s})^{2}/\sigma_{\text{CDPP}}$, we choose
a typical error of 10% for $\sigma_{\text{CDPP}}$, to match the uncertainty
induced by stellar radius and planet-star radius ratio.

To get an intuitive view of how well the parameters have been controlled, we
plot Figure 5 and Figure 6 in which we perform parameter control for the cases
of three age bins and five age bins. In the first row of Figure 5 and Figure
6, we show the Cumulative Distribution Function (CDF) diagrams of temperature,
mass, metallicity, radius, and $\sigma_{\text{CDPP}}$ 3.5 hours for the
observation data. Using the above method, we control all parameters and show
the CDF of controlled star samples in the bottom row. By applying the
parameter control method, we have achieved the goal to let stars in different
age bins have similar distribution in stellar temperature, mass, metallicity,
radius, and $\sigma_{\text{CDPP}}$ (Figure 5 and 6).

Figure 5: Cumulative Distribution Function (CDF) diagrams of effective
temperature, mass, metallicity, radius, and $\sigma_{\text{CDPP}}$ 3.5 hours
for the three bins method before and after parameter control. The errorbar in
the upper left corner shows the typical error of each stellar property, and
the number in the lower right corner presents the size of the standard sample.
Figure 6: Similar to Figure 5, here we show CDF diagrams for the five bins
case. The blue, orange, green, red, and purple lines represent star data
sample with TD/D in the ranges of (0,0.0427), [0.0427–0.0618),
[0.0618-0.0.0808), [0.0808–0.184], and $[0.184,+\inf$).

###  3.2 Apparent Planetary Occurrence as a Function of Age

We first consider a three bins case and calculate the kinematic age of each
bin using AVR (Appendix A). We adopt the above controlled sample and calculate
the numbers of stars that have one, two, three or more planets, i.e., the
planet multiplicity function ($N_{1}$, $N_{2}$, and $N_{3+}$). The apparent
occurrence rate of one, two, and three or more planet systems is derived by
dividing $N_{1}$, $N_{2}$, and $N_{3+}$ by the number of stars ($N_{star}$) in
each bin. In Figure 7, from the left column to the right column, we present
the apparent occurrence rate for one, two, and three or more planet systems.
As can be seen, the young stars generally have a higher apparent occurrence
than the old stars. For the original data without parameter control (first row
of Figure 7), $N_{1}/N_{star}$, $N_{2}/N_{star}$, and $N_{3+}/N_{star}$ are
$2.12^{+0.18}_{-0.17}$%, $0.61^{+0.10}_{-0.09}$%, and $0.35^{+0.08}_{-0.07}$%
for stars less than 1 Gyr, which are about 4.8$\sigma$, 4.4$\sigma$, and
4.8$\sigma$ higher than those ($1.37^{+0.15}_{-0.13}$%,
$0.26^{+0.07}_{-0.06}$%, and $0.08^{+0.05}_{-0.03}$%) for stars about 8 Gyr,
respectively. For the data after all parameter control (bottom row of Figure
7), $N_{1}/N_{star}$, $N_{2}/N_{star}$, and $N_{3+}/N_{star}$ are
$1.86^{+0.25}_{-0.22}$%, $0.46^{+0.14}_{-0.11}$%, and $0.31^{+0.12}_{-0.09}$%
for the youngest group, which are about 1.6$\sigma$, 1.4$\sigma$, and
3.3$\sigma$ higher than those ($1.50^{+0.22}_{-0.20}$%,
$0.31^{+0.12}_{-0.09}$%, and $0.05^{+0.07}_{-0.03}$%) for the oldest group,
respectively. The differences in apparent rate between young and old groups
become smaller when taking into account of parameter control. Such a change is
more prominent for low multiplicity systems ($N_{1}/N_{star}$ and
$N_{2}/N_{star}$) than for high multiplicity systems ($N_{3+}/N_{star}$). We
also calculate the Pearson correlation coefficients and $p$-values for the
correlations between age and apparent planet occurrence, which are printed in
each panel of Figure 7. The $p$-values are derived using the following steps.

  1. (1)

We calculate the Pearson correlation coefficients for the observation data as
$\rho_{obs}$.

  2. (2)

We generate simulated apparent occurrence rates for each bin assuming Poisson
error, and randomly scramble their order. Then we calculate the Pearson
correlation coefficient between age and the scrambled data ($\rho_{sim}$).

  3. (3)

We repeat Step 2 for 10,000 times and calculate the fraction of the simulated
data that produce a stronger anti-correlation, i.e., $\rho_{sim}<\rho_{obs}$.
This fraction gives the $p$-value of the Pearson correlation for the observed
data.

As we can see, the anti-correlations between age and $N_{1}/N_{star}$,
$N_{2}/N_{star}$ become weaker after parameter control, with the $p$-values
rising from 0.0978 to 0.145 and from 0.0237 to 0.181, respectively.
Nevertheless, the anti-correlation between age and $N_{3+}/N_{star}$ becomes
stronger, with the $p$-value decreasing from 0.056 to 0.0117.

Figure 7: Apparent planet occurrence rate for the three bins case with one
(left), two (middle), and three or more (right) transiting planets. The top
and bottom rows correspond to the results before and after parameter control,
as shown in Figure 5. The dots and errorbars present the median value and $\pm
1\sigma$ range, assuming Poisson error. The numbers at the top of each panel
show the corresponding planet system numbers ($N_{1}$, $N_{2}$, and $N_{3+}$).
In the corner of each panel, we also print the Pearson correlation coefficient
($\rho$) and the corresponding $p$-value.

To further trace the apparent planet occurrence trend with age, we group the
stars into five age bins, and calculate the corresponding kinematic age. In
Figure 8, we present the apparent planet occurrence rate as a function of
kinematic age for systems with one, two, and three or more planets. In line
with the results of the three bins case, we find that the apparent occurrence
rate generally has a declining trend with kinematic age. For the original data
(before parameter control), $N_{1}/N_{star}$, $N_{2}/N_{star}$, and
$N_{3+}/N_{star}$ are $2.17^{+0.23}_{-0.21}$%, $0.76^{+0.15}_{-0.13}$%, and
$0.35^{+0.11}_{-0.08}$% for the youngest stars, which are 4.6$\sigma$,
4.2$\sigma$, and 3.6$\sigma$ higher than those ($1.26^{+0.18}_{-0.16}$%,
$0.29^{+0.10}_{-0.08}$%, and $0.08^{+0.07}_{-0.04}$%) for the oldest stars,
respectively (top row of Figure 8). For the data after all parameter control,
$N_{1}/N_{star}$, $N_{2}/N_{star}$, and $N_{3+}/N_{star}$ are
$1.65^{+0.35}_{-0.29}$%, $0.57^{+0.23}_{-0.17}$%, and $0.41^{+0.20}_{-0.14}$%
for stars less than 1 Gyr, which are about 0.5$\sigma$, 1.1$\sigma$, and
2.8$\sigma$ higher than those ($1.50^{+0.33}_{-0.28}$%,
$0.36^{+0.19}_{-0.13}$%%, and $0.05^{+0.12}_{-0.04}$%) for stars about 8 Gyr,
respectively (bottom row of Figure 8). Again, being consistent with the
results in the three bins case, the differences in the apparent occurrence
rate of stars with different ages become smaller after parameter control.
Nevertheless, the differences are still significant ($\sim 3\sigma$) for
systems with high transiting multiplicities ($N_{3+}/N_{star}$, right column
of Figure 8). We also calculate the Pearson correlation coefficients and
$p$-values, as in the three bins case, and print them in each panel of Figure
8. Similar to the three bins case, the anti-correlation between age and
$N_{1}/N_{star}$ becomes weaker after parameter control, with the $p$-value
rising from 0.0075 to 0.191. Nevertheless, the anti-correlations between age
and multiple planet systems ($N_{2}/N_{star}$ and $N_{3+}/N_{star}$) become a
little stronger, with the $p$-values dropping from 0.212 to 0.0108 and from
0.0134 to 0.0011, respectively.

We have also employed Canonical Correlation Analysis (CCA) to investigate the
relationship between stellar properties and the apparent planet occurrence.
The CCA method derives similar results as shown above, indicating that stellar
age is anti-correlated with planet occurrence without the need for performing
parameter control. However, the CCA method is unable to identify the star and
planet samples required for the forward modeling method (see Section 4) in
order to derive the intrinsic planet occurrence. For more detailed
information, please refer to Appendix B.

Figure 8: Apparent planet occurrence as a function of kinematic age for
systems with one (left), two (middle), and three or more (right) transiting
planets. The top and bottom rows correspond to the results before and after
parameter control as shown in Figure 6. The dots and errorbars present the
median value and $\pm 1\sigma$ range assuming Poisson error. The numbers on
the top of each panel show the corresponding planet system numbers ($N_{1}$,
$N_{2}$, and $N_{3+}$). In the corner of each panel, we also print the Pearson
correlation coefficient ($\rho$) and the corresponding $p$-value.

##  4 INTRINSIC trend analysis from forward modeling

###  4.1 Forward Modeling Method

The above apparent occurrence rates only reflect the observed planet
population. In order to derive the intrinsic planet occurrence rates of the
underlying planet population, we adopt a forward modeling method that takes
into account of _transit_ observation bias and detection/vetting efficiencies.
The framework of the method has been described in detail in Zhu et al. (2018)
and Yang et al. (2020). In this section, we summarize the general procedure
and emphasize the modifications considered in this work.

####  4.1.1 General procedure of the modeling

We derive the observed planet multiplicity function ($N_{1}$, $N_{2}$, and
$N_{3+}$) from the star sample. Then we generate simulated planet systems, and
calculate the modeled multiplicity function ($\bar{N}_{1}$, $\bar{N}_{2}$, and
$\bar{N}_{3+}$). For our star sample with a size of $\sim 20,000$, we need to
generate about $\sim 10,000$ planet systems, assuming on average 50% of the
stars own planet systems (the true value of $F_{\text{Kep}}$ differs for each
group and is automatically adjusted in the MCMC process). The total number of
generated planet systems is about 400,000,000 when the simulation is
converged. We assume the multiplicity function follows a Poisson distribution,
and optimize the likelihood function

$\mathcal{L}=\prod\limits_{k=1}^{3+}\frac{\bar{N}_{k}^{N_{k}}\exp(-\bar{N}_{k})}{N_{k}!}$
(2)

with python package emcee (Foreman-Mackey et al., 2013) applying the Markov
Chain Monte Carlo (MCMC) method. Three free parameters are constrained in our
model, which are the fraction of stars with Kepler-like planets
($F_{\text{Kep}}$), the average planet number for stars that have Kepler-like
planets ($\bar{N}_{p}$), and the inclination slope index ($\alpha$, see
below). Details of generating the modelled multiplicity function
($\bar{N}_{1}$, $\bar{N}_{2}$, and $\bar{N}_{3+}$) can be seen in Yang et al.
(2020). We briefly summarize the general procedure as follows:

  1. (1)

Assuming the intrinsic planet occurrence. For each group of stars, we assume
that a fraction of $F_{\text{Kep}}$ percent of stars have Kepler-like planets,
and on average, each planet system has $\bar{N}_{p}$ planets. For each host
star, we generate $k$ planets following a zero-truncated Poisson distribution
of $\bar{N}_{p}$ (Fang & Margot, 2012).

  2. (2)

Assuming transit parameters and radii for planets. We generate the debiased
distributions of transit parameter ($\epsilon$, $\epsilon=R_{s}/a_{p}$,
$R_{s}$ is the stellar radius and $a_{p}$ is the semi-major axis of the
planet) and planet radii ($R_{p}$) considering three kinds of bias (Mulders et
al., 2018), namely, the transit geometry bias ($f_{\text{tra}}$), detection
efficiency bias ($f_{\text{S/N}}$), and vetting efficiency bias
($f_{\text{vet}}$). For each planet, we assign values of $\epsilon$ and
$R_{p}$ that are randomly drawn from the debiased distributions.

  3. (3)

Adjusting period ratios and radius ratios. Planets within the same system tend
to have similar period ratios (Fabrycky et al., 2014; Brakensiek & Ragozzine,
2016; Jiang et al., 2020) and radius ratios (Ciardi et al., 2013; Weiss et
al., 2018). To account for this correlation, we adjust the period ratios and
radius ratios for planets within the same system. These adjustments are based
on debiased distributions calculated by CORBITS (Brakensiek & Ragozzine, 2016)
for multiple planet systems.

  4. (4)

Checking Orbital stability. To ensure that the simulated planetary systems are
physically plausible, we assess their orbital stability. We apply the
criterion provided by Deck et al. (2013) to prevent planets within the same
system from being located too close to each other.

  5. (5)

Assigning orbital inclination to generate transits. We calculate the
inclination ($I_{p}$) of the planets with respect to the observer by

$\cos I_{p}=\cos I\cos i\ -\ \sin I\sin i\cos\phi,$ (3)

where $I$ represents the inclination of the system invariable plane, $\phi$ is
the phase angle, and $i$ is the inclination of the planet with respect to the
invariable plane. Both $I$ and $\phi$ are assumed to be isotropic. Following
Zhu et al. (2018), for a planet system with $k$ planets, the inclination
dispersion of the planets follows a power-law function,

$\sigma_{i,k}\equiv\sqrt{\left<\sin^{2}i\right>}=\sigma_{i,5}\left(\frac{k}{5}\right)^{\alpha},$
(4)

where we adopt $\sigma_{i,5}$ as a Gaussian distribution with a mean of $0\fdg
8$ and a standard deviation of $0\fdg 15$ from Zhu et al. (2018). We fit the
inclination slope index $\alpha$ as the third parameter. A planet is
considered to be a transit if its impact factor is less than 1
($|\cos{I_{p}}/\epsilon|<1$).

  6. (6)

Checking detection and vetting efficiencies. Due to detection and vetting
efficiencies, not all the transiting planets can be detected. For each
transiting planet, we generate a random number from a uniform distribution
ranging from 0 to 1. We consider this planet can be detected if the generated
random number is less than the product of the detection efficiency
($f_{\text{S/N}}$, see Appendix C) and the vetting ($f_{\text{vet}}$)
efficiency.

####  4.1.2 Emphasize the difference from Yang et al. (2020)

Comparing to our previous work, we do not consider the TTV multiplicity
function, namely, the number of systems that show a TTV signal. The TTV
function is omitted for two reasons. First, the number of stars in our sample
in this work is less than 20,000, and the number of systems that show a TTV
signal is only 31. This is much smaller compared to the 100,000 star sample
and 127 systems that show TTV signals in Yang et al. (2020). The smaller
number of the TTV multiplicity function leads to larger uncertainty. Second,
as shown in the Appendix of Yang et al. (2020), although removing the TTV
multiplicity function from the likelihood leads to less constraint on the
parameter $\alpha$, it has little impact on the results of $F_{\text{Kep}}$
and $\bar{N}_{p}$, which are the core parameters of this work.

###  4.2 Intrinsic Planetary Occurrence and Architecture as a Function of Age

We show the forward modeling results for the case of three bins in Figure 9.
From the left panels to the right panels, we present the posterior
distributions of $F_{\text{Kep}}$, $\bar{N}_{p}$, $\alpha$, and $\eta$ (which
is the product of $F_{\text{Kep}}$ and $\bar{N}_{p}$, see Equation 1). For
data without parameter control, the youngest group generally has higher
intrinsic planet occurrence rates than the oldest group. $F_{\text{Kep}}$,
$\bar{N}_{p}$, and $\eta$ are $63.2^{+11.4}_{-9.4}$%, $2.71^{+0.43}_{-0.40}$,
and $1.71^{+0.16}_{-0.16}$ for stars less than 1 Gyr, which are 1.7$\sigma$,
1.8$\sigma$, and 5.4$\sigma$ higher than those ($47.5^{+8.9}_{-8.4}$%,
$1.97^{+0.41}_{-0.32}$, and $0.93^{+0.13}_{-0.12}$) for stars about 8 Gyr,
respectively (top row of Figure 9). For the data after all parameter control,
$F_{\text{Kep}}$, $\bar{N}_{p}$, and $\eta$ are $56.9^{+13.3}_{-11.6}$%,
$2.74^{+0.64}_{-0.49}$, and $1.57^{+0.23}_{-0.22}$ for the youngest group,
which are 0.2$\sigma$, 2.1$\sigma$, and 3.0$\sigma$ higher than those
($54.3^{+12.3}_{-11.4}$%, $1.75^{+0.47}_{-0.32}$, and $0.96^{+0.19}_{-0.17}$)
for the oldest group, respectively (bottom row of Figure 9). All the three
groups have nearly the same $F_{\text{Kep}}$ when parameter control is taken
into account.

This result is expected because $F_{\text{Kep}}$ is mainly determined by the
apparent occurrence rate of single planet systems ($N_{1}/N_{star}$). The
difference in apparent occurrence rate for single planet systems between the
youngest and oldest groups drops from 4.8$\sigma$ to 1.6$\sigma$ after
parameter control (Figure 7), and consequently, the difference of
$F_{\text{Kep}}$ drops from 1.7$\sigma$ to 0.2$\sigma$.

$\bar{N}_{p}$ is largely determined by the apparent occurrence rate of high
multiplicity systems ($N_{3+}/N_{star}$). The difference in $N_{3+}/N_{star}$
between the youngest and oldest groups drops mildly from 4.8$\sigma$ to
3.3$\sigma$ after parameter control (Figure 7). Interestingly, the difference
in $\bar{N}_{p}$ increases slightly from 1.8$\sigma$ to 2.1$\sigma$. Due to
the limited number of planetary systems with three or more planets in the last
bin after parameter control (2 systems in the three bins case), the Poisson
error is relatively high ($2^{+1.8}_{-1.3}$). We may have underestimated the
declining trend of $N_{3+}/N_{star}$, and the results of the forward modeling
show that the decrease in $\bar{N}_{p}$ becomes slightly more prominent. As
$\eta$ is the product of $F_{\text{Kep}}$ and $\bar{N}_{p}$, therefore, the
decrease of the difference in $\eta$ from 5.4$\sigma$ to 3.0$\sigma$ is mainly
due to the decrease of the difference in $F_{\text{Kep}}$. We calculate the
Pearson correlation coefficients and $p$-values for the correlations between
age and the intrinsic planet occurrence ($F_{\text{Kep}}$, $\bar{N}_{p}$, and
$\eta$), and print them in each panel in Figure 9. The $p$-values are derived
using a similar method as shown in Section 3.2. For $F_{\text{Kep}}$ and
$\bar{N}_{p}$, the anti-correlations between age and them are statistically
insignificant before and after parameter control with $p$-values larger than
0.05. The anti-correlation between age and $\eta$ maintains with $p$-value
smaller than 0.05.

The parameter $\alpha$ is not well constrained in all the cases before and
after parameter control, because it is mainly constrained by TTV multiplicity
function (Zhu et al., 2018), which is ignored in this work.

Figure 9: Posterior distributions of $F_{\text{Kep}}$, $\bar{N}_{p}$,
$\alpha$, and $\eta$ for the three bins case are presented. The top and bottom
rows show the forward modeling results corresponding to samples before and
after parameter control as in Figure 5 and Figure 7. The dots and errorbars
show the 50% and $\pm 1$–$\sigma$ range. In the upper right of the panels in
the first, second, and fourth columns, we also print the Pearson correlation
coefficient ($\rho$) and the corresponding $p$-value.

To further investigate the intrinsic planet occurrence as a function of age,
we group the stars into five bins as mentioned before in Section 3.2, and
adopt the forward modeling method to derive $F_{\text{Kep}}$, $\bar{N}_{p}$,
$\alpha$, and $\eta$. The results are shown in Figure 10. For data without
parameter control (top row of Figure 10), the values of $F_{\text{Kep}}$,
$\bar{N}_{p}$, and $\eta$ are $68.0^{+12.8}_{-11.9}$%, $2.68^{+0.52}_{-0.41}$,
and $1.82^{+0.22}_{-0.20}$, respectively, for stars in the youngest group.
These values are 2.1$\sigma$, 1.2$\sigma$, and 4.8$\sigma$ higher than those
($43.7^{+10.9}_{-8.9}$%, $2.09^{+0.54}_{-0.39}$, and $0.92^{+0.17}_{-0.15}$)
for stars in the oldest group. After parameter control (bottom row of Figure
10), $F_{\text{Kep}}$ is $47.1^{+17.9}_{-16.1}$% for stars less than 1 Gyr,
which is 0.6$\sigma$ lower than stars around 8 Gyr ($56.6^{+17.4}_{-14.9}$%).
As to $\bar{N}_{p}$ and $\eta$, they are $3.69^{+1.58}_{-0.96}$ and
$1.71^{+0.35}_{-0.31}$ for stars in the first bin, which are 2.3$\sigma$ and
2.2$\sigma$ higher than the values ($1.80^{+0.67}_{-0.39}$ and
$1.04^{+0.29}_{-0.24}$) in the last bin.

The forward modeling results for the five bins are consistent with those for
the three bins (Figure 9). After parameter control, the difference in
$N_{1}/N_{star}$ between young and old stars shows a significant decrease from
4.6 $\sigma$ to 0.5 $\sigma$ (Figure 8). The difference in $F_{\text{Kep}}$
decreases from 2.1$\sigma$ to 0.6$\sigma$. The first bin has a higher
$N_{1}/N_{star}$ compared to the last bin, however, it shows a slightly lower
value of $F_{\text{Kep}}$. That is because the first bin has more intrinsic
multi-planet systems, which can also contribute to the apparent occurrence of
$N_{1}/N_{star}$. The difference in $N_{3+}/N_{star}$ between young and old
stars drops moderately from 3.6$\sigma$ to 2.8$\sigma$. Forward modeling
indicates that the difference in $\bar{N}_{p}$ increases slightly from
1.2$\sigma$ to 2.3$\sigma$. Similar to the three bins case, we have only one
planet system with three or more planets in the last bin, resulting in a high
Poisson error ($1^{+2.3}_{-0.8}$). As a consequence, we might underestimate
the declining trend of $N_{3+}/N_{star}$. As for $\eta$ (the product of
$F_{\text{Kep}}$ and $\bar{N}_{p}$) the difference drops from 4.8$\sigma$ to
2.2$\sigma$, which is mainly due to the decrease of the difference in
$F_{\text{Kep}}$.

The Pearson correlation coefficients and $p$-values for the correlations
between age and the intrinsic planet occurrence ($F_{\text{Kep}}$,
$\bar{N}_{p}$, and $\eta$) are also printed in the corner of each panel in
Figure 10. Regarding the anti-correlation between age and $F_{\text{Kep}}$, it
becomes weaker after parameter control with $p$-value rising from 0.0503 to
0.729. The anti-correlation between age and $\bar{N}_{p}$ becomes
statistically significant, with $p$-value dropping from 0.0839 to 0.0052. The
anti-correlations between age and $\eta$ maintain, with $p$-value changing
from 0.0021 to 0.0096.

Similar to the three bins case, $\alpha$ is not well constrained before and
after parameter control.

Figure 10: Similar to Figure 9, posterior distributions of $F_{\text{Kep}}$,
$\bar{N}_{p}$, $\alpha$, and $\eta$ for the five bins case are presented. The
top and bottom rows show the forward modeling results corresponding to samples
before and after parameter control, as shown in Figure 6 and Figure 8. The
dots and errorbars show the 50% and $\pm 1$–$\sigma$ range. In the upper right
of the panels in the first, second, and fourth columns, we also print the
Pearson correlation coefficient ($\rho$) and the corresponding $p$-value.

##  5 Discussions

In this paper, we first investigate the apparent planet occurrence in terms of
$N_{1}/N_{star}$, $N_{2}/N_{star}$, and $N_{3+}/N_{star}$, then from which we
derive the intrinsic planet occurrence in terms of $F_{\text{Kep}}$,
$\bar{N}_{p}$, and $\eta$ as a function of stellar age using a forward
modeling method. We have applied a parameter control method in our analyses to
remove the effects caused by other stellar properties. We find that after
parameter control, younger stars generally have higher apparent occurrence
than older stars. Specifically, the intrinsic planet occurrence in terms of
the number of planets per star ($\eta$) decreases with stellar age with a
confidence level of about $2\sim 3\sigma$. Such a declining trend is mainly
driven by the decrease in the average multiplicity ($\bar{N}_{p}$, by about
2$\sigma$), and partially by the change in the fraction of stars with planet
systems ($F_{\text{Kep}}$, by less than 1$\sigma$). In what follows, we will
compare our results to those from literature and discuss the implications of
these findings for our understanding of planet formation and evolution.

###  5.1 Giant and Ultra Short Period Planets

We exclude the giant and Ultra Short Period (USP) planets from our planet
sample in Section 2.2. We dub this planet sample as the ‘Fiducial’ sample. To
investigate the influence of giants and USPs on planet occurrence, in this
section, we re-run our simulation including the giants, the USPs, and both of
them. In Figure 11, from the top to the bottom, we show the results for the
Fiducial sample, the sample including giant planets, the sample including
USPs, and the sample including both of giant planets and USPs for the three
bins case. As we can see, after applying the parameter control method, all
results show similar trends. For $F_{\text{Kep}}$, the differences between the
youngest and oldest groups are less than 1$\sigma$. The youngest groups
generally have $\bar{N}_{p}$ $\sim 2\sigma$ higher than the oldest groups, and
$2\sim 3\sigma$ higher for $\eta$. Including giant planets and USPs adds more
planets into our sample, leading to a slightly higher value of
$F_{\text{Kep}}$. At the same time, since giant planets are more likely to be
detected in single planet systems, including giant planets causes a very small
decrease in $\bar{N}_{p}$.

Figure 11: Posterior distributions of $F_{\text{Kep}}$, $\bar{N}_{p}$,
$\alpha$, and $\eta$ for the three bins case are presented. From the top to
the bottom, each row shows the forward modeling results after parameter
control for our Fiducial planet sample, the sample including the giants, the
sample including the USPs, and the sample including both of the giants and the
USPs, respectively. The dots and errorbars show the 50% and $\pm 1\sigma$
range of the posterior distributions. In the upper right corner of the panels
in the first, second, and fourth columns, we also print the Pearson
correlation coefficient ($\rho$) along with the corresponding $p$-value.

Similar to the three bins case, the results for the five bins case are
basically unchanged after including giants and USPs. As we can see in Figure
12, the anti-correlations between age and $F_{\text{Kep}}$ are statistically
insignificant, with $p$-values higher than 0.5. For $\bar{N}_{p}$ and $\eta$,
the $p$-values are less than 0.05, showing significant anti-correlations
between them and age.

Figure 12: Similar to Figure 11, posterior distributions of $F_{\text{Kep}}$,
$\bar{N}_{p}$, $\alpha$, and $\eta$ for five bins case are shown. From the top
to the bottom rows, each row shows the forward modeling results corresponding
to the Fiducial sample, the sample including the giants, the sample including
the USPs, and the sample including both of the giants and the USPs,
respectively. The dots and errorbars show the 50% and $\pm 1\sigma$ range of
the posterior distributions. In the upper right corner of the panels in the
first, second, and fourth columns, we also print the Pearson correlation
coefficient ($\rho$) and the corresponding $p$-value.

###  5.2 Comparison with Previous Studies

McTier & Kipping (2019) have studied the dependence of planet occurrence rate
on Galactocentric velocity. After correcting for selection biases, they found
that Kepler planet hosts have a similar velocity distribution to the non-host
Kepler stars. Based on such a similarity, they inferred that the planet
occurrence rate is independent on Galactocentric velocity. Their inference is
against our results, which show that planet occurrences in terms of
$\bar{N}_{p}$ and $\eta$ are anti-correlated with kinematic age and thus with
Galactocentric velocity based on the Age Velocity dispersion Relation (AVR).
In fact, we argue that their inference may not necessarily be valid for the
following reasons.

First, since the occurrence rates of Kepler planets are generally found to be
high ($\sim$50%, Mulders et al., 2018; Yang et al., 2020; He et al., 2021), a
large fraction of the _apparent_ ‘non-host’ stars are actually hosts of
planets that are not detected by Kepler. Therefore, it is not surprising that
Kepler planet hosts have a velocity distribution similar to that of the non-
host Kepler stars as found by McTier & Kipping (2019).

Second, multiple transiting systems count multiple times when counting
planets, but only one time when counting host star. Therefore, multiples,
which play a critical role in deriving the planet occurrence (e.g.,
$\bar{N}_{p}$) have little effect on determining the velocity distribution of
host stars. In fact, as shown in Figure 8 of PAST II (Chen et al., 2021a),
Kepler planet hosts are dominated by single transiting systems, which have a
velocity distribution similar to that of non-host stars.

Our results in PAST-II have shown that multiple transiting systems have
significantly different Galactocentric velocity distributions compared to the
single transiting systems. Such differences lead to the occurrence trends with
Galactocentric velocity and thus with kinematic age seen in this work (Figure
9 and Figure 10, while still maintaining the similarity in Galactocentric
velocity distribution between the Kepler planet hosts and the non-host Kepler
stars seen in (McTier & Kipping, 2019). In one word, the similarity in
Galactocentric velocity between Kepler planet hosts and non-host Kepler stars
does not necessarily infer that the intrinsic occurrence of Kepler planets is
independent of Galactocentric velocity.

In a series of papers, Bashi & Zucker (2019, 2022) have studied the planet
occurrence rate in the Galactic context. Bashi & Zucker (2019) found that for
stars with metallicity higher than $-0.25$, the planet occurrence rate
generally decreases with the Galactocentric velocity of host stars. Bashi &
Zucker (2022) found that the planet occurrence, in terms of the fraction of
stars with planets and the number of planets per star, are higher in the
Galactic thin disk stars than in the thick disk stars. In addition, Bashi &
Zucker (2022) also showed an apparent anti-correlation between planet
occurrence and the stellar isochrone age. Generally speaking, their results
are consistent with ours. In this paper, we find that planet occurrence rate
decreases with TD/D and age. Nevertheless, we emphasize some differences in
our work as compared to theirs. First, we use the kinematic age, which has
relatively smaller internal uncertainty of $\sim$10–20% (Chen et al., 2021b)
compared to isochrone age with a typical uncertainty of up to 56% (Berger et
al., 2018). Second, we use the parameter control method to isolate the effect
of age from other stellar properties. After removing these effects, we find
that the anti-correlations are weaker between planet occurrence and age,
especially for $F_{\text{Kep}}$, though they remain significant for $\eta$
(Figure 9 and 10).

###  5.3 Implications to Planet Formation and Evolution

In this study, we have revealed observational evidence that the occurrence and
architecture (in terms of $F_{\text{Kep}}$ and $\bar{N}_{p}$) of Kepler
planetary systems evolve over time. To gain deeper insights into planet
formation and evolution, one would compare our observational results to
theoretical models. Unfortunately, we did not find any models that predict
$F_{\text{Kep}}$ or $\bar{N}_{p}$ as a function of time, which would allow for
a quantitative comparison with our findings. Nevertheless, there are still
theoretical and numerical works in the literature that allow us to make
qualitative comparisons.

Systems with more than two bodies are generally chaotic, and essentially
unstable. Planetary systems are usually formed with more than one planet, and
their architecture will be further shaped during the long-term dynamical
evolution afterwards, e.g., triggered by orbital instability. The timescale of
the orbital instability depends on many factors, such as mass, number,
eccentricity, and orbit spacing of the planets within the system. For a
planetary system born with a large number of planets and tight orbital
spacing, the orbital instability occurs quickly, which causes planet ejections
and collisions, leading to a decrease in the planet number and an increase in
orbital spacing (Zhou et al., 2007). This in turn increases the timescale of
subsequently instability, which means the system needs to evolve on a longer
time scale to trigger next instability. As the system evolves, the instability
timescale can grow to as long as billions of years. Our Solar System may have
undergone such an evolutionary process (Tsiganis et al., 2005; Liu et al.,
2022). In some models of our Solar System (e.g., Nesvorný & Morbidelli, 2012),
it is thought there were initially five or even more giant planets formed in a
tightly packed orbital configuration. The current architecture of the Solar
System was mainly shaped by an orbital instability events that ejected at
least one giant and scatted other planets into a more loosely packed
configuration. For exoplanet systems, Pu & Wu (2015) found that the orbital
spacing of Kepler multi-planet systems are clustered around the threshold of
orbital instability. Based on this observation they hypothesized that most of
the Kepler systems were formed with tighter spacing configuration, and most of
them have undergone orbit instability, leading to fewer planets left on larger
orbit spacing. Using $N$-body simulations, Izidoro et al. (2017) proposed an
evolutionary scenario for the bulk of Kepler-planet (super-Earths or mini-
Neptunes) systems. In this scenario, planets were formed in a compact resonant
chains through migration in proto-planetary gas disk in the early stage. As
the gas dissipated, the chains became dynamically unstable, which led to
planets merging, ejection, and being scattered to form a spread out
configuration.

In this work, using a forward modeling method and after applying parameter
control, we find that the planet occurrence rate in terms of planet number per
star ($\eta$) decreases by about 2–3$\sigma$ as a function of time. For three
bins case, as shown in Figure 9, $\eta$ drops from 1.57 to 0.96, and for the
five bins case in Figure 10, it decreases from 1.71 to 1.04.

The first major contribution to the $\eta$ decreasing trend comes from the
planet number in planetary system, i.e., $\bar{N}_{p}$, since $\eta$ is the
product of $\bar{N}_{p}$ and $F_{\text{Kep}}$. $\bar{N}_{p}$ shows a moderate
decline about $\sim 2\sigma$ in our fitting results. In the bottom row of
Figure 9, $\bar{N}_{p}$ drops from 2.74 for stars less than 1 Gyr to 1.75 for
stars about 8 Gyr, and in Figure 10, $\bar{N}_{p}$ drops from 3.69 for the
first age group to 1.80 for the last group. This is qualitatively consistent
with the above theories, that the dynamical evolution of planetary systems
causes the merging and ejecting of planets. Furthermore, from our results we
infer that the evolution of $\bar{N}_{p}$ can continue to several gigayears,
which implies that planetary systems keep evolving through the whole stellar
lifetime.

The second potential contribution to the $\eta$ decreasing trend comes from
the fraction of star that have planet, i.e., $F_{\text{Kep}}$. In Figure 9,
$F_{\text{Kep}}$ changes from 56.9% to 54.3%, and in Figure 10, it changes
from 47.1% to 56.6%, both are less than 1$\sigma$. Due to the limited star and
planet sample, we cannot conclude that the change in $F_{\text{Kep}}$ is
statistically significant. Future studies with larger samples of planetary
systems may help us to further constrain $F_{\text{Kep}}$ as a function of
time and unveil the planet formation rate in the history of the Milky Way,
combining the information on the star formation rate as a function of age
(Binney et al., 2000).

Not only does the number of planets in a system evolve with time, but the
orbital properties also undergo changes. Since $\bar{N}_{p}$ is related to the
orbital inclination (Equation 4), we can investigate the mutual orbital
inclination as a function of time. We calculate the posterior distributions of
$\sigma_{i,k}$ for the five bins after parameter control, and show the
distributions as well as the median value and 1$\sigma$ range of
$\sigma_{i,k}$ in Figure 13. As we can see, $\sigma_{i,k}$ gradually evolves
with time. From less than 1 Gyr to about 8 Gyr, the median value of
$\sigma_{i,k}$ grows from about $1\fdg 2$ to $3\fdg 5$, and the $1\sigma$
range expands from $0\fdg 7$–$2\fdg 6$ to $1\fdg 3$–$11\fdg 7$. To further
quantify the age-$\sigma_{i,k}$ trend, we fit $\sigma_{i,k}$ with age.
Although bearing large uncertainty (as seen from the orange shaded region),
the best fit is

$\log{\sigma_{i,k}}=0.2+0.4\times\log{\frac{\text{Age}}{\text{Gyr}}}.$ (5)

For comparison, we also plot the data points for our Solar System and Kepler
multiple transiting systems in Figure 13. They are all generally fit such an
age-$\sigma_{i,k}$ trend. This result indicates that as planetary systems get
older, they become dynamically hotter, which is consistent with the
theoretical expectation (Zhou et al., 2007). The Kepler multiples show smaller
mutual inclinations compared to our Solar System, which can be explained by
their younger ages according to the age-inclination trend shown in Figure 13.
In other words, Figure 13 may hint that planets in our Solar System were in a
flatter architecture in the early time, then gradually evolve to the current
state.

Figure 13: Inclination dispersion ($\sigma_{i,k}$) as a function of age. We
show the distribution of $\sigma_{i,k}$ for five age bins after parameter
control with blue bars. The median value and the $\pm 1\sigma$ percentile
range of $\sigma_{i,k}$ are indicated by orange dots and errorbars, with the
corresponding $\pm 1\sigma$ percentile range is displayed in the lower part of
the plot. The orange dashed line represents the best fit, and the orange
shaded region denotes the corresponding uncertainties of the fit by resampling
$\sigma_{i,k}$ 10,000 times. The $\sigma_{i,k}$ value of the Solar System
terrestrial planets is shown with the red star marker, and the $\sigma_{i,k}$
values of Kepler multiple transiting systems measured by Fabrycky et al.
(2014) and by Xie et al. (2016) are shown with the green and purple star
markers, respectively. The age of the Solar System is adopted as 4.57 Gyr
(Bouvier & Wadhwa, 2010), and the age of Kepler multiples is adopted as the
kinematic age of stars that host two or more planets in our star sample.

##  6 Summary and Conclusion

In this work, which is the fourth paper of the PAST series, we update the
LAMOST-Gaia-Kepler catalog utilizing the recently released LAMOST DR8 and Gaia
DR3. Based on this catalog, we study the occurrence rate and architecture of
Kepler-like planets as a function of stellar kinematic age. We find the
following results.

  1. (1)

Younger stars generally show higher apparent planet occurrence rates (more
than 3$\sigma$) for one, two, and three or more planets ($N_{1}/N_{star}$,
$N_{2}/N_{star}$, and $N_{3+}/N_{star}$) than older stars (top rows of Figure
7 and Figure 8).

  2. (2)

Applying a parameter control method can effectively reduce the effects caused
by other stellar properties, such as effective temperature, mass, metallicity,
radius, and $\sigma_{\text{CDPP}}$. After parameter control, the differences
in $N_{1}/N_{star}$ and $N_{2}/N_{star}$ between younger and older stars
decrease to less than 2$\sigma$, while the difference in $N_{3+}/N_{star}$
maintains a confidence level of about 3$\sigma$ (bottom rows of Figure 7 and
Figure 8).

  3. (3)

Adopting a forward modeling method can help us to investigate the intrinsic
planet occurrence in terms of the fraction of stars with planetary systems
($F_{\text{Kep}}$), the average planet multiplicity ($\bar{N}_{p}$), and the
number of planets per star ($\eta$). For stars without parameter control, we
find that the younger stars have higher $F_{\text{Kep}}$ and $\bar{N}_{p}$ by
about 2$\sigma$ than older stars. The difference in $\eta$ between younger and
older stars is more obvious, at about 5$\sigma$ level (top rows of Figure 9
and Figure 10).

  4. (4)

After parameter control, the differences in $F_{\text{Kep}}$ drops to less
than 1$\sigma$, hinting the planetary system occurrence remains at a similar
rate throughout the history of the Milky Way. The difference in $\bar{N}_{p}$
is about 2$\sigma$ between the younger and older stars. This result is
consistent with theories that planet systems keep evolving as a result of the
merging and ejecting of the planets. Younger stars have a higher $\eta$ (the
product of $F_{\text{Kep}}$ and $\bar{N}_{p}$) by about $2\sim 3\sigma$ than
older stars, which is the combining effects caused by the evolution of
$F_{\text{Kep}}$ and $\bar{N}_{p}$ (bottom rows of Figure 9 and Figure 10).

  5. (5)

The orbital properties of planet systems also evolve with time. We find that
in stars aging from less than 1 Gyr to about 8 Gyr, the mutual orbital
inclination ($\sigma_{i,k}$) between their planets increases from $1\fdg 2$ to
$3\fdg 5$, and the $\pm 1\sigma$ range of $\sigma_{i,k}$ expands from $0\fdg
7$–$2\fdg 6$ to $1\fdg 3$–$11\fdg 7$ (Figure 13), hinting that planet systems
become dynamically hotter as a function of time. Both our Solar System and
Kepler multiple transiting systems fit such a trend.

Our work qualitatively agrees with theoretical expectations that planet
occurrence decreases, and planetary systems become dynamically hotter with
age. Future dedicated theoretical and numerical modeling on the occurrence and
architecture of Kepler planets as a function of age are needed to allow us to
make quantitative comparisons to our results in this work and to place key
constraints on planet formation and evolution.

The current and upcoming missions also aid in exploring exoplanets in the
dimension of time. The TESS mission has found thousands of candidates
(Guerrero et al., 2021), covering a wide range of ages (e.g., Newton et al.,
2019; Gan et al., 2020; Weiss et al., 2021). In this paper, our studies only
rely on a portion of planets from the Kepler sample. Both the sample size and
the number of bins are still limited, leading to relatively large
uncertainties in our fitting results. In the near future, missions such as
Gaia, TESS (Ricker et al., 2015), and PLATO (Rauer et al., 2014) will detect
many more exoplanets, leading to the expansion of planet sample by one order
of magnitude or even more. With the help of more data, future studies will
further refine our measurements and test our results.

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, 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). D.-C.C. also
acknowledges the Cultivation project for LAMOST Scientific Payoff and Research
Achievement of CAMS-CAS. This work has included data from Guoshoujing
Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope,
LAMOST), which is a National Major Scientific Project built by the Chinese
Academy of Sciences. Funding for LAMOST (www.lamost.org) has been provided by
the Chinese NDRC. LAMOST is operated and managed by the National Astronomical
Observatories, CAS. This work presents results from the European Space Agency
(ESA) space mission Gaia. Gaia data are being processed by the Gaia Data
Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by
national institutions, in particular the institutions participating in the
Gaia MultiLateral Agreement (MLA).

##  Appendix A Updating kinematic Characteristics of Galactic components and
Age-velocity dispersion relation with Gaia DR3 astrometry data

In the first paper of the PAST series (PAST I, Chen et al., 2021b), we
revisited the kinematic method to classify the Galactic components and the
Age-Velocity dispersion Relation (AVR) with Gaia DR2 (Gaia Collaboration et
al., 2018a, b) and the LAMOST main-sequence turn-off and subgiant (MSTO-SG)
star sample (Xiang et al., 2017). On June 13th, 2022, Gaia Data Release 3
(DR3; Gaia Collaboration et al., 2022) was released, providing five
astrometric parameters (positions, parallaxes, and proper motions) for 1.468
billion sources. Comparing to Gaia DR2, the standard uncertainties have been
reduced for the positions, parallaxes, and proper motions, which make the
astrometric results considerably more robust and reduce the systematic errors.
Therefore, here we revisit the kinematic methods and AVR with Gaia DR3 by
adopting the same procedures shown in the Section 2 and 3 of PAST I.

###  A.1 Updated Calibration Sample

To construct the calibration sample, we first cross-match the above LAMOST
MSTO-SG catalog with the Gaia DR3 catalog using the X-match service provided
by the Centre de Donnees astronomiques de Strasbourg (CDS,
http://cdsxmatch.u-strasbg.fr). Second we carry out a angular distance cut of
1.25 arcseconds and a Gaia G-band magnitude difference cut of 2.5 mag. For
stars with multiple matches, we keep those with the smallest angular
separation.

Then we calculate the stellar kinematic properties (i.e., Galactocentric
cylindrical coordinates $(R,\theta,Z)$, and Galactic rectangular velocities
$(U_{\rm LSR},V_{\rm LSR},W_{\rm LSR})$) relative to the local standard of
rest (LSR), with the procedure detailed in Section 2.1 of PAST I. We adopt the
location of the Sun of $R_{\odot}$= 8.18 kpc (Gravity Collaboration et al.,
2019; Gaia Collaboration et al., 2021) and $Z_{\odot}$ = 25 pc (Chen et al.,
2001). The solar peculiar motions are taken as [$U_{\odot}$, $V_{\odot}$,
$W_{\odot}$] = [9.58, 10.52, 7.01] $\rm km\ s^{-1}$ (Tian et al., 2015).

After that, we apply the following filters to further clean the calibration
sample.

(1) Binary filter. We remove binary systems because their kinematics contain
additional motions (Dehnen & Binney, 1998). This is done by choosing stars
flagged as ‘Normal star’ (i.e., single stars with spectral types of AFGKM) in
the LAMOST MSTO-SG catalog (Xiang et al., 2017). We also remove potential
binaries by eliminating stars with Gaia DR3 re-normalized unit-weight error
(RUWE) $>1.4$ (Lindegren, 2018).

(2) Parallax precision filter. Following Dehnen & Binney (1998), we remove
stars with relative parallax errors larger than 10 percent as reported in the
Gaia DR3.

(3) Age precision filter. We remove stars with ages older than 14 Gyr, errors
of age exceeding 25%, or blue straggler stars ($|Z|>1.5$ kpc and ages younger
than 2 Gyr) in the LAMOST MSTO-SG catalog.

(4) Distance filter (similar to Binney et al., 1997). The majority of the
remaining stars are brighter than G mag=16, where the median parallax error is
0.0494 mas. Recalling the above 10 percent parallax precision requirement,
this translates to a distance limit $\sim 1/(0.0494/0.1)\sim 2.0$ kpc. We
therefore remove stars with distances exceeding than this limit.

After applying the above filters, we are left with 134,244 stars, which are
mainly (129,089/134,244, 96.2%) located at $7.5<R<10.0$ kpc, $|\theta|<15$
deg, and $|Z|<1.5$ kpc.

###  A.2 Revisiting the Kinematic Method To Classify the Galactic Components

With stellar kinematic and age following the same criteria in the Section
2.3.3 of PAST I, we then classify the calibration sample into different
Galactic components, i.e., thin disk ($D$), thick disk ($TD$), halo ($H$), and
Hercules stream ($Herc$). In order to calculate the characteristic kinematic
parameters for each Galactic component as a function of $(R,\ Z)$, we bin the
calibration sample as the same interval in PAST I. For $|Z|$, we set 8 bins
with boundaries at $|Z|=$ 0, 0.1, 0.2, 0.3, 0.4, 0.55, 0.75, 1.0, and 1.5 kpc.
For $R$, we set 5 bins with boundaries at $R=$ 7.5, 8.0, 8.5, 9.0, 9.5, and 10
kpc. In total, there are $5\times 8=40$ grids in the $R$–$Z$ space, and all
bins have enough ($>400$) stars.

We then revise the normalized fraction $X$ (Equations 9, 10, 11, 12 in PAST I)
and the velocity ellipsoid (i.e., $\sigma_{U}$, $\sigma_{V}$, $\sigma_{W}$,
and $V_{\rm asym}$ (Equation 3 in PAST I) of each Galactic component for each
grid in the $R$–$Z$ plane following the same procedure of Section 2.3.4 and
2.3.5 of PAST I. The calculated values of $X$, $\sigma_{U}$, $\sigma_{V}$,
$\sigma_{W}$, and $V_{\rm asym}$ are tabulated in Table 2 and visualized in
Figure 14, Figure 15, and Figure 16.

Figure 14 shows the $X$ values of various Galactic components as functions of
Galactic radius $R$ and absolute value of height, $|Z|$. As expected, $X_{\rm
D}$ ($X_{\rm TD},\ X_{\rm H}$) generally decreases (increases) with $|Z|$ in
all the $R$ bins.

With the same procedure detailed described in Section 2.3.4 of PAST I, we also
fit the velocity dispersions, $\sigma_{U}$, $\sigma_{V}$, and $\sigma_{W}$ in
the following formula according to Williams et al. (2013):

$\sigma=b_{1}+b_{2}\times\frac{R}{\rm kpc}+b_{3}\times(\frac{Z}{\rm
kpc})^{2}{\rm km\ s^{-1}}.$ (A1)

We then use the following formula to calculate $V_{\rm asym}$ according to
Robin et al. (2003); Binney & Tremaine (2008)

$V_{\rm asym}=\sigma_{U}^{2}/C_{0}.$ (A2)

The values of fitting parameters and their 1$\sigma$ uncertainties are
summarized in Table 3.

Comparing to the results obtained from Gaia DR2 and LAMOST MSTO-SG sample in
PAST I, we find that for the normalized fraction $X$, the typical (median)
relative differences are only 0.6%, -3.5%, 4.6%, and 5.2% for the thin disk,
thick disk, halo, and Hercules stream, respectively. For the velocity
ellipsoid (i.e., $\sigma_{U}$, $\sigma_{V}$, $\sigma_{W}$, and $V_{\rm asym}$)
obtained with the calibration sample using astrometry data from Gaia DR2 and
DR3, as can be seen in Figure 15 and 16, the median values and $1\sigma$
errorbars are very similar, and the best fits are nearly the same with each
other. It can also be seen from Table 4, the fitting parameters (i.e.,
$b_{1}$, $b_{2}$, $b_{3}$, and $C_{0}$) are well consistent with those of PAST
I within their $1\sigma$ errorbars. Therefore, we conclude that the $X$
factors and velocity ellipsoid obtained with the updated calibration sample
are broadly unchanged from those of PAST I.

###  A.3 Revisiting the Age-Velocity dispersion relation (AVR)

According to PAST I, we divide the calibration sample into 30 bins with
approximately equal sizes ($\sim$ 4,475 stars in each bin) according to their
ages. Then we fit the AVRs following Holmberg et al. (2009b); Aumer et al.
(2016) by using a simple power law formula, i.e.,

$\sigma=k\times\left(\frac{t}{\rm Gyr}\right)^{\beta}\,\rm km\ s^{-1},$ (A3)

where $t$ represents stellar age, $\sigma$ is the velocity dispersion, and $k$
and $\beta$ are two fitting parameters. The best fits and uncertainties
($1\sigma$ interval) of the fitting parameters $(k,\ \beta)$ are calculated
with the same procedure described in Section 3.1 of PAST I and summarized in
Table 4.

Figure 17 shows the velocity dispersion as a function of the median age of
each bin. As can be seen, the best fits (black lines) for the relationship
between age and the dispersion of velocity components ($U_{\rm LSR},V_{\rm
LSR},W_{\rm LSR}$) and the total velocity ($V_{\rm tot}$) are all
indistinguishable from those of PAST I (red lines). In Table 4, we compare the
fitting parameters $(k,\ \beta)$ of AVRs obtained from the updated calibration
sample to those from PAST I. As can be seen, the median values are nearly
identical, and the $1\sigma$ uncertainties of $k$ decrease by a factor of
$\sim 10\%$ due to the improvement on the precision of stellar parallax and
proper motion measurements. Thus, the AVRs derived with calibration sample
using astrometry data from Gaia DR3 and DR2 are nearly the same with each
other.

##  Appendix B Deriving the Age-Planet Occurrence Relationship through
Canonical Correlation Analysis

Canonical Correlation Analysis (CCA) is an effective method to discover the
correlations between different sets of variables. It was first introduced by
(HOTELLING, 1936). The basic idea is to identify linear combinations of two
sets of variables that the resulting combined variables exhibit the highest
possible correlation. KETTENRING (1971) summarized various methods for
establishing connections among multiple sets of variables. Here, we apply the
CCA method to investigate the relationship between planet and star properties.
Specifically, we use the SABSCOR method to maximize the sum of the absolute
values of the correlation among different sets.

We group the planet and star properties into three sets of variables: the
planet system property (occurrence rate), the interesting stellar property
(kinematic age), and the stellar properties we want to eliminate (effective
temperature, mass, metallicity, radius, and $\sigma_{\text{CDPP}}$),
respectively.

The planet system property we are interested in is the apparent occurrence
rate. Following Zhu (2019), we use two tracers to represent the occurrence
rate, which are defined by the following equations:

$\displaystyle\text{Tracer}(\eta)$
$\displaystyle=\frac{1}{N_{star}}\sum_{j=1}^{K}jN_{j}$ (B1)
$\displaystyle\text{Tracer}({F_{p}})$
$\displaystyle=\frac{1}{N_{star}}\left(N_{1}+\sum_{j=1}^{K}N_{j}\right).$ (B2)

Tracer($\eta$) is related to the average number of planets per star, and
Tracer($F_{p}$) is correlated with the fraction of stars processing planet
systems. In these equations, $N_{star}$ represents the number of stars in each
bin, $N_{j}\ (j=1,2...K)$ is the number of systems with $j$ planets, and $K$
is the maximum number of planets we observed in the planetary system. Given
the transit method can only detect a very small fraction of planets, in our
sample of 19,358 stars, we have observed only 641 planets in 467 planetary
systems. Therefore, to calculate these two tracers, we need to group stars
into bins. Additionally, the distribution of planet systems is not uniform,
especially for systems with three or more planets, with fewer systems found
around older stars. To reduce Poisson errors, it is necessary to limit the
number of bins. Consequently, we group the stars into 40 bins with an equal
number of stars (using 30 or 50 bins yields similar results).

The interesting stellar property, kinematic age, is determined using the Age-
Velocity dispersion Relationship (AVR, see Appendix A). This relationship
requires the calculation of velocity dispersion and, as a result, is
applicable only to a group of stars. Therefore, we once again need to divide
our star sample into bins.

For the uninteresting stellar properties, we use the median value of each bin
as the representative value. Given the limited number of bins (40), to prevent
overfitting, we select only three properties (mass, [Fe/H], and
$\sigma_{\text{CDPP}}$) instead of all five. Previous studies (e.g., Yang et
al., 2020; He et al., 2021) have shown that the mass and effective temperature
of stars have a similar influence on planet occurrence. From late type to
early type stars, the increase in mass and temperature leads to a decrease in
planet occurrence. We choose mass to represent the influence of stellar type.
Furthermore, an increase in radius and $\sigma_{\text{CDPP}}$ both result in a
reduction in detection efficiency, which lowers the probability of planet
detection. Here we choose $\sigma_{\text{CDPP}}$ to present the effect of
detection efficiency (choosing temperature or radius leads to a similar
result).

In summary, we categorize all the planet/star properties into three groups.
The first group, $\mathbf{X}_{1}$={Tracer($\eta$), Tracer($F_{p}$)},
represents planet occurrence. The second group, $\mathbf{X}_{2}$={Age},
describes the interesting stellar properties, specifically, stellar kinematic
age. The third group, $\mathbf{X}_{3}$={Mass, [Fe/H], $\sigma_{\text{CDPP}}$},
is related to uninteresting stellar properties.

We use the R package RGCCA (Girka et al., 2023), to maximize the sum of the
correlations between planet occurrence and stellar age, as well as between
planet occurrence and uninteresting star properties. This can be formulated as
solving the following optimization problem:

$\text{Maximize}\
(|\text{Cor}(\mathbf{y}_{1},\mathbf{y}_{2})|+|\text{Cor}(\mathbf{y}_{1},\mathbf{y}_{3})|).$
(B3)

where $\mathbf{y}_{1}=\mathbf{X}_{1}\mathbf{a}_{1}$,
$\mathbf{y}_{2}=\mathbf{X}_{2}\mathbf{a}_{2}$, and
$\mathbf{y}_{3}=\mathbf{X}_{3}\mathbf{a}_{3}$, and $\mathbf{a}_{1}$,
$\mathbf{a}_{2}$, and $\mathbf{a}_{3}$ are the weight vectors for each
variable set. All the variables have been standardized. The results can be
seen in the following figure.

In Figure 18, each ellipse represents a group, and each box represents a
star/planet property. We also print the weight and correlation on each line.
As we can see, both Tracer($\eta$) and Tracer($F_{p}$) have positive
contributions to the planet occurrence ($\mathbf{a}>0$), showing that as the
number of planets per star and the fraction of stars with planets increase,
the planet occurrence rises. The kinematic age is in anti-correlated with
planet occurrence (Cor=-0.802), which is consistent with our result before
parameter control (see top rows of Figure 7 and Figure 8). As stars aging from
less than 1 Gyr to about 8 Gyr, both of the fraction of stars with planetary
systems and the number of planets per star decrease. The uninteresting star
properties are correlated with planet occurrence. Among these properties, the
stellar mass shows a negative weight, which is in agreement with literature
(e.g. Yang et al., 2020; He et al., 2021). An increase in mass results in a
decrease in Kepler-like planet occurrence. Metallicity shows a positive
weight, which is generally consistent with previous studies (e.g., Zhu, 2019;
Wang & Fischer, 2015), indicating that an increase in metallicity can
stimulate the formation of planets. $\sigma_{\text{CDPP}}$ demonstrates a
negative weight because higher $\sigma_{\text{CDPP}}$ leads to lower detection
efficiency, which hinders the detection of planets.

The CCA method shows a similar result to ours before parameter control,
indicating that the increase in stellar age results in a decrease in the
planet occurrence rate. This decrease is reflected in both the fraction of
stars with planetary systems and the number of planets per star.

##  Appendix C Detection Efficiency

We show the average detection efficiencies and planet samples for both the
three and five bins cases in Figure 19 and 20. The detection efficiency
metrics are calculated by the Package KeplerPORTs (Burke & Catanzarite, 2017),
and the associated data are downloaded from NASA exoplanet archive555https://
exoplanetarchive.ipac.caltech.edu/docs/.

As we can see in the top rows of Figure 19 and 20, young stars have slightly
higher detection efficiencies than old stars. This is because first, young
stars generally have smaller stellar radii, which lead to deeper transit
depths, and second, young stars have lower noise levels
($\sigma_{\text{CDPP}}$) that increase the signal to noise ratio (top rows of
Figure 5 and Figure 6). After we apply parameter control to remove the effects
caused by stellar properties, the young and old stars have similar
distributions of stellar radii and $\sigma_{\text{CDPP}}$ (bottom rows of
Figure 5 and Figure 6). As a result, the average detection efficiencies in
each bin (red lines) are similar to the mean value of the whole sample (black
lines, bottom rows of Figure 19 and 20), showing that the influence caused by
detection efficiencies on planet occurrence is effectively removed.

Table 2: Revised characteristics at different Galactic radii ($R$) and heights ($Z$) for different Galactic components using the calibration sample updated with Gaia DR3 astrometry data. $|Z|$ | $R$ | $\sigma^{\rm D}_{\rm U}$ | $\sigma^{\rm D}_{\rm V}$ | $\sigma^{\rm D}_{\rm W}$ | $V^{\rm D}_{\rm asym}$ | $\sigma^{\rm TD}_{\rm U}$ | $\sigma^{\rm TD}_{\rm V}$ | $\sigma^{\rm TD}_{\rm W}$ | $V^{\rm TD}_{\rm asym}$ | $X_{\rm D}$ | $X_{\rm TD}$ | $X_{\rm H}$ | $X_{\rm Herc}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
(kpc) | (kpc) |  ————————— ($\rm km\ s^{-1}$) —————————
$0-0.1$ | $7.5-8.0$ | 35 | 20 | 16 | $-14$ | 61 | 37 | 31 | -41 | 0.81 | 0.13 | 0.0010 | 0.06
$8.0-8.5$ | 35 | 21 | 16 | $-15$ | 63 | 39 | 35 | $-44$ | 0.83 | 0.11 | 0.0014 | 0.06
$8.5-9.0$ | 33 | 20 | 15 | $-13$ | 69 | 37 | 34 | $-51$ | 0.85 | 0.10 | 0.0013 | 0.05
$9.0-9.5$ | 31 | 19 | 15 | $-12$ | 69 | 35 | 33 | $-52$ | 0.87 | 0.09 | 0.0012 | 0.04
$9.5-10.0$ | 28 | 19 | 15 | $-10$ | 67 | 37 | 33 | $-47$ | 0.89 | 0.08 | 0.0016 | 0.03
$0.1-0.2$ | $7.5-8.0$ | 36 | 21 | 16 | $-15$ | 63 | 37 | 36 | $-44$ | 0.78 | 0.16 | 0.0015 | 0.06
$8.0-8.5$ | 36 | 22 | 17 | $-16$ | 64 | 39 | 36 | -45 | 0.78 | 0.14 | 0.0014 | 0.08
$8.5-9.0$ | 34 | 20 | 16 | $-15$ | 70 | 38 | 36 | -53 | 0.82 | 0.12 | 0.0014 | 0.06
$9.0-9.5$ | 32 | 20 | 15 | $-13$ | 70 | 39 | 34 | $-52$ | 0.84 | 0.11 | 0.0015 | 0.05
$9.5-10.0$ | 29 | 19 | 15 | $-11$ | 69 | 37 | 34 | $-51$ | 0.86 | 0.10 | 0.0019 | 0.04
$0.2-0.3$ | $7.5-8.0$ | 37 | 23 | 17 | $-15$ | 66 | 39 | 37 | $-46$ | 0.74 | 0.19 | 0.0016 | 0.07
$8.0-8.5$ | 38 | 23 | 18 | $-16$ | 67 | 42 | 40 | $-47$ | 0.75 | 0.17 | 0.0015 | 0.08
$8.5-9.0$ | 36 | 21 | 17 | $-15$ | 71 | 41 | 37 | $-54$ | 0.79 | 0.15 | 0.0017 | 0.06
$9.0-9.5$ | 33 | 21 | 16 | $-13$ | 70 | 40 | 36 | $-52$ | 0.81 | 0.14 | 0.0016 | 0.05
$9.5-10.0$ | 30 | 20 | 15 | $-11$ | 66 | 39 | 35 | $-47$ | 0.84 | 0.12 | 0.0023 | 0.04
$0.3-0.4$ | $7.5-8.0$ | 38 | 23 | 18 | $-16$ | 67 | 41 | 38 | -48 | 0.70 | 0.23 | 0.0016 | 0.07
$8.0-8.5$ | 40 | 23 | 18 | $-18$ | 67 | 41 | 41 | $-51$ | 0.72 | 0.20 | 0.0020 | 0.08
$8.5-9.0$ | 37 | 22 | 17 | $-15$ | 70 | 41 | 38 | $-53$ | 0.75 | 0.18 | 0.0020 | 0.07
$9.0-9.5$ | 35 | 21 | 17 | $-14$ | 70 | 41 | 37 | $-52$ | 0.78 | 0.17 | 0.0020 | 0.05
$9.5-10.0$ | 34 | 20 | 16 | $-13$ | 70 | 40 | 35 | $-52$ | 0.80 | 0.15 | 0.0025 | 0.05
$0.4-0.55$ | $7.5-8.0$ | 42 | 25 | 19 | $-19$ | 70 | 42 | 40 | -53 | 0.65 | 0.28 | 0.0021 | 0.07
$8.0-8.5$ | 41 | 23 | 19 | $-18$ | 69 | 41 | 42 | $-51$ | 0.66 | 0.25 | 0.0026 | 0.09
$8.5-9.0$ | 39 | 23 | 19 | $-16$ | 70 | 42 | 38 | $-53$ | 0.70 | 0.23 | 0.0022 | 0.07
$9.0-9.5$ | 39 | 22 | 18 | $-17$ | 70 | 42 | 38 | $-53$ | 0.73 | 0.21 | 0.0024 | 0.06
$9.5-10.0$ | 35 | 21 | 17 | $-15$ | 69 | 40 | 35 | $-51$ | 0.76 | 0.19 | 0.0026 | 0.05
$0.55-0.75$ | $7.5-8.0$ | 43 | 26 | 20 | $-21$ | 70 | 45 | 41 | -52 | 0.55 | 0.36 | 0.0051 | 0.09
$8.0-8.5$ | 43 | 24 | 20 | $-21$ | 71 | 43 | 42 | $-54$ | 0.58 | 0.32 | 0.0042 | 0.09
$8.5-9.0$ | 41 | 24 | 19 | $-19$ | 71 | 43 | 40 | $-54$ | 0.62 | 0.30 | 0.0026 | 0.08
$9.0-9.5$ | 39 | 22 | 18 | $-17$ | 70 | 41 | 40 | $-53$ | 0.66 | 0.28 | 0.0026 | 0.06
$9.5-10.0$ | 38 | 21 | 17 | $-16$ | 70 | 41 | 40 | $-53$ | 0.68 | 0.25 | 0.0029 | 0.07
$0.75-1.0$ | $7.5-8.0$ | 44 | 28 | 23 | $-20$ | 70 | 43 | 44 | -53 | 0.44 | 0.46 | 0.0069 | 0.09
$8.0-8.5$ | 45 | 25 | 21 | $-22$ | 72 | 43 | 43 | $-56$ | 0.47 | 0.43 | 0.0078 | 0.09
$8.5-9.0$ | 41 | 25 | 20 | $-19$ | 72 | 45 | 41 | $-56$ | 0.50 | 0.41 | 0.0048 | 0.09
$9.0-9.5$ | 40 | 23 | 19 | $-18$ | 71 | 43 | 40 | $-54$ | 0.54 | 0.38 | 0.0034 | 0.08
$9.5-10.0$ | 39 | 22 | 19 | $-17$ | 72 | 42 | 40 | $-55$ | 0.58 | 0.35 | 0.0033 | 0.07
$1.0-1.5$ | $7.5-8.0$ | 47 | 33 | 25 | $-25$ | 71 | 45 | 46 | -54 | 0.25 | 0.63 | 0.0169 | 0.10
$8.0-8.5$ | 46 | 30 | 24 | $-22$ | 72 | 45 | 44 | $-56$ | 0.26 | 0.62 | 0.0157 | 0.10
$8.5-9.0$ | 43 | 29 | 22 | $-20$ | 72 | 45 | 43 | $-55$ | 0.32 | 0.58 | 0.0151 | 0.09
$9.0-9.5$ | 41 | 26 | 21 | $-19$ | 73 | 45 | 40 | $-57$ | 0.35 | 0.56 | 0.0104 | 0.08
$9.0-9.5$ | 40 | 25 | 20 | $-18$ | 71 | 42 | 40 | $-54$ | 0.39 | 0.54 | 0.0082 | 0.07
Table 3: Fitting parameters of the velocity dispersion as functions of ($R,\ Z$, i.e., Equation A1) and asymmetric velocity as a function of $\sigma^{2}_{U}$ (i.e., Equation A2) derived with the calibration sample using astrometry data from Gaia DR3 and DR2. |  | ——— Gaia DR3 ——— | ——— Gaia DR2 ———
---|---|---|---
|  | Thin disk | Thick disk | Thin disk | Thick disk
$\sigma_{U}$ | $b_{1}$ | $65.9^{+1.1}_{-2.9}$ | $57.7^{+5.9}_{-4.0}$ | $63.4^{+1.3}_{-3.2}$ | $58.4^{+6.7}_{-4.4}$
$b_{2}$ | $-3.5^{+0.3}_{-0.1}$ | $1.2^{+0.6}_{-0.7}$ | $-3.2^{+0.3}_{-0.2}$ | $1.2^{+0.6}_{-0.7}$
$b_{3}$ | $6.8^{+0.5}_{-1.3}$ | $3.9^{+0.4}_{-0.7}$ | $7.6^{+0.5}_{-1.3}$ | $4.1^{+0.4}_{-0.8}$
$\sigma_{V}$ | $b_{1}$ | $41.2^{+0.6}_{-2.1}$ | $46.3^{+4.4}_{-3.3}$ | $41.6^{+0.7}_{-2.3}$ | $44.9^{+5.3}_{-2.7}$
$b_{2}$ | $-2.3^{+0.3}_{-0.1}$ | $-0.8^{+0.5}_{-0.5}$ | $-2.3^{+0.3}_{-0.1}$ | $5.0^{+0.2}_{-0.7}$
$b_{3}$ | $5.8^{+0.3}_{-0.5}$ | $4.8^{+0.4}_{-1.0}$ | $5.6^{+0.2}_{-0.9}$ | $5.2^{+0.3}_{-1.0}$
$\sigma_{W}$ | $b_{1}$ | $29.5^{+1.1}_{-1.2}$ | $55.2^{+1.8}_{-2.5}$ | $27.3^{+1.3}_{-1.5}$ | $55.8^{+1.8}_{-3.1}$
$b_{2}$ | $-1.5^{+0.1}_{-0.2}$ | $-2.1^{+0.3}_{-0.5}$ | $-1.2^{+0.2}_{-0.1}$ | $-2.2^{+0.4}_{-0.2}$
$b_{3}$ | $5.0^{+0.3}_{-0.8}$ | $5.8^{+0.5}_{-0.6}$ | $5.0^{+0.2}_{-0.7}$ | $6.1^{+0.4}_{-1.0}$
$C_{0}$ | $-89.7^{+1.5}_{-1.5}$ | $-93.0^{+2.4}_{-1.6}$ | $-88.5^{+1.7}_{-1.9}$ | $-92.5^{+2.8}_{-2.1}$
Table 4: Fitting parameters of the Age-Velocity dispersion Relationship with the calibration samples using astrometry data from Gaia DR3 and DR2. |  ——— $k\ \rm(km\ s^{-1})$ ——— |  ——— $\beta$ ———
---|---|---
| value |  1 $\sigma$ interval | value |  1 $\sigma$ interval
Gaia DR3
$U$ | $23.74$ | (23.47, 24.65) | $0.34$ | $(0.32,0.35)$
$V$ | $12.87$ | (12.47, 13.37) | $0.42$ | $(0.40,0.44)$
$W$ | $8.29$ | (7.92, 8.63) | $0.56$ | $(0.54,0.58)$
$V_{\rm tot}$ | $27.74$ | (27.32, 28.67) | $0.39$ | $(0.37,0.41)$
Gaia DR2, PAST I
$U$ | $23.66$ | (23.07, 24.32) | $0.34$ | $(0.33,0.36)$
$V$ | $12.49$ | (12.05, 12.98) | $0.43$ | $(0.41,0.45)$
$W$ | $8.50$ | (8.09, 8.97) | $0.54$ | $(0.52,0.56)$
$V_{\rm tot}$ | $27.55$ | (26.84, 28.37) | $0.40$ | $(0.38,0.42)$
Figure 14: The normalisation fraction $X$ of stars for each component as
functions of Galactic radius $R$ and absolute value of height $|Z|$. The
different colours denote subsamples of stars with different Galactic radii.
Figure 15: The velocity dispersions as functions of position ($R$, $|Z|$) in
the Galaxy. For the updated calibration data conducted from the LAMOST MSTO-SG
and Gaia DR3 catalogs, the velocity dispersions are plotted as solid points,
and the line segments represent $1\sigma$ errors in two colors: blue for thin
disk and red for thick disk. The solid line in each panel denotes the result
of the best fit of Equation A1 using the coefficients in Table 3. While for
the calibration sample conducted from the LAMOST MSTO-SG and Gaia DR2 catalogs
in the PAST I, the velocity dispersions are plotted as solid points, and the
line segments represent $1\sigma$ errors in two colors: light blue for thin
disk and light red for thick disk. The dashed black line in each panel denotes
the result of the best fit of Equation A1 using the coefficients in Table 3.
Figure 16: The asymmetric velocity, $V_{\rm asym}$ as a function of
$\sigma_{U}^{2}$ for the thin disk (left panel) and thick disk (right panel).
For the updated calibration data conducted from the LAMOST MSTO-SG and Gaia
DR3 catalogs, the data is plotted as blue/red points and the blue/red line
segments represent $1\sigma$ errors. The blue/red solid lines denote the
results of the best fit using Equation A2. While for the calibration sample
conducted from the LAMOST MSTO-SG and Gaia DR2 catalogs in the PAST I, the
data and $1\sigma$ errors are plotted in light blue/red color. The light
blue/red dashed lines denote the results of the best fit using Equation A2.
Figure 17: The velocity dispersions for $U_{\rm LSR},V_{\rm LSR},W_{\rm LSR}$,
and $V_{\rm tot}$ vs. age for the selected calibration star sample. The solid
black lines denote the respective best fit of refitting AVR (Equation A3)
using the coefficients in Table 4.  Figure 18: Results of Canonical
Correlation Analysis for planet and stellar properties are provided, including
the weights assigned to each property and the correlations between properties
and groups.  Figure 19: Detection efficiencies and planet samples in the
period-radius diagram for the three bins case. From top to bottom, each row
corresponds to the control samples in Figure 5 and Figure 7. The reds lines
present the average 90%, 50%, and 10% detection efficiencies for the star
sample in each bin, the grey lines show the mean 90%, 50%, and 10% detection
efficiencies for the whole star sample. The blue, orange, and green dots show
planets in one, two, and three or more planet systems, respectively.  Figure
20: Similar to Figure 19, 90%, 50%, and 10% detection efficiencies and planet
samples for the five bins case. From top to bottom, each row corresponds to
control samples in Figure 6 and Figure 8.

## References

  * Akeson. et al. (2013) Akeson., R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989, doi: 10.1086/672273
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Aumer et al. (2016) Aumer, M., Binney, J., & Schönrich, R. 2016, MNRAS, 462, 1697, doi: 10.1093/mnras/stw1639
  * Bashi & Zucker (2019) Bashi, D., & Zucker, S. 2019, AJ, 158, 61, doi: 10.3847/1538-3881/ab27c9
  * Bashi & Zucker (2022) —. 2022, MNRAS, 510, 3449, doi: 10.1093/mnras/stab3596
  * Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527, doi: 10.1051/0004-6361:20031213
  * Bensby et al. (2014) Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71, doi: 10.1051/0004-6361/201322631
  * Berger et al. (2018) Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, ApJ, 866, 99, doi: 10.3847/1538-4357/aada83
  * Berger et al. (2020a) Berger, T. A., Huber, D., Gaidos, E., van Saders, J. L., & Weiss, L. M. 2020a, AJ, 160, 108, doi: 10.3847/1538-3881/aba18a
  * Berger et al. (2020b) Berger, T. A., Huber, D., van Saders, J. L., et al. 2020b, AJ, 159, 280, doi: 10.3847/1538-3881/159/6/280
  * Binney et al. (2000) Binney, J., Dehnen, W., & Bertelli, G. 2000, MNRAS, 318, 658, doi: 10.1046/j.1365-8711.2000.03720.x
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Binney et al. (1997) Binney, J. J., Dehnen, W., Houk, N., Murray, C. A., & Penston, M. J. 1997, in ESA Special Publication, Vol. 402, Hipparcos - Venice ’97, ed. R. M. Bonnet, E. Høg, P. L. Bernacca, L. Emiliani, A. Blaauw, C. Turon, J. Kovalevsky, L. Lindegren, H. Hassan, M. Bouffard, B. Strim, D. Heger, M. A. C. Perryman, & L. Woltjer, 473–478
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Bouvier & Wadhwa (2010) Bouvier, A., & Wadhwa, M. 2010, Nature Geoscience, 3, 637, doi: 10.1038/ngeo941
  * Brakensiek & Ragozzine (2016) Brakensiek, J., & Ragozzine, D. 2016, ApJ, 821, 47, doi: 10.3847/0004-637X/821/1/47
  * Bryson et al. (2020) Bryson, S., Coughlin, J., Batalha, N. M., et al. 2020, AJ, 159, 279, doi: 10.3847/1538-3881/ab8a30
  * Buchhave et al. (2012) Buchhave, L. A., Latham, D. W., Johansen, A., et al. 2012, Nature, 486, 375, doi: 10.1038/nature11121
  * Burke & Catanzarite (2017) Burke, C. J., & Catanzarite, J. 2017, Planet Detection Metrics: Per-Target Detection Contours for Data Release 25, Kepler Science Document KSCI-19111-002
  * Chen et al. (2001) Chen, B., Stoughton, C., Smith, J. A., et al. 2001, ApJ, 553, 184, doi: 10.1086/320647
  * Chen et al. (2021a) Chen, D.-C., Yang, J.-Y., Xie, J.-W., et al. 2021a, AJ, 162, 100, doi: 10.3847/1538-3881/ac0f08
  * Chen et al. (2021b) Chen, D.-C., Xie, J.-W., Zhou, J.-L., et al. 2021b, ApJ, 909, 115, doi: 10.3847/1538-4357/abd5be
  * Chen et al. (2022) —. 2022, AJ, 163, 249, doi: 10.3847/1538-3881/ac641f
  * Chevance et al. (2021) Chevance, M., Kruijssen, J. M. D., & Longmore, S. N. 2021, ApJ, 910, L19, doi: 10.3847/2041-8213/abee20
  * Christiansen et al. (2012) Christiansen, J. L., Jenkins, J. M., Caldwell, D. A., et al. 2012, PASP, 124, 1279, doi: 10.1086/668847
  * Ciardi et al. (2013) Ciardi, D. R., Fabrycky, D. C., Ford, E. B., et al. 2013, ApJ, 763, 41, doi: 10.1088/0004-637X/763/1/41
  * Dai et al. (2018) Dai, F., Masuda, K., & Winn, J. N. 2018, ApJ, 864, L38, doi: 10.3847/2041-8213/aadd4f
  * Dai. et al. (2021) Dai., Y.-Z., Liu, H.-G., An, D.-S., & Zhou, J.-L. 2021, AJ, 162, 46, doi: 10.3847/1538-3881/ac00ad
  * David et al. (2021) David, T. J., Contardo, G., Sandoval, A., et al. 2021, AJ, 161, 265, doi: 10.3847/1538-3881/abf439
  * Deck et al. (2013) Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129, doi: 10.1088/0004-637X/774/2/129
  * Dehnen & Binney (1998) Dehnen, W., & Binney, J. J. 1998, MNRAS, 298, 387, doi: 10.1046/j.1365-8711.1998.01600.x
  * 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
  * Fang & Margot (2012) Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92, doi: 10.1088/0004-637X/761/2/92
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102, doi: 10.1086/428383
  * Fontanive et al. (2019) Fontanive, C., Rice, K., Bonavita, M., et al. 2019, MNRAS, 485, 4967, doi: 10.1093/mnras/stz671
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018b, A&A, 616, A11, doi: 10.1051/0004-6361/201832865
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211.  https://arxiv.org/abs/2208.00211
  * Gan et al. (2020) Gan, T., Shporer, A., Livingston, J. H., et al. 2020, AJ, 159, 160, doi: 10.3847/1538-3881/ab775a
  * Ghezzi et al. (2018) Ghezzi, L., Montet, B. T., & Johnson, J. A. 2018, ApJ, 860, 109, doi: 10.3847/1538-4357/aac37c
  * Girka et al. (2023) Girka, F., Camenen, E., Peltier, C., et al. 2023, RGCCA: Regularized and Sparse Generalized Canonical Correlation Analysis for Multiblock Data.  https://CRAN.R-project.org/package=RGCCA
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A, 625, L10, doi: 10.1051/0004-6361/201935656
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39, doi: 10.3847/1538-4365/abefe1
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * He et al. (2021) He, M. Y., Ford, E. B., & Ragozzine, D. 2021, AJ, 161, 16, doi: 10.3847/1538-3881/abc68b
  * Holmberg et al. (2009a) Holmberg, J., Nordström, B., & Andersen, J. 2009a, A&A, 501, 941, doi: 10.1051/0004-6361/200811191
  * Holmberg et al. (2009b) —. 2009b, A&A, 501, 941, doi: 10.1051/0004-6361/200811191
  * HOTELLING (1936) HOTELLING, H. 1936, Biometrika, 28, 321, doi: 10.1093/biomet/28.3-4.321
  * Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15, doi: 10.1088/0067-0049/201/2/15
  * Huber (2019) Huber, D. 2019, evolstate: Assign simple evolutionary states to stars, Astrophysics Source Code Library, record ascl:1905.003.  http://ascl.net/1905.003
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ida & Lin (2004) Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388, doi: 10.1086/381724
  * IRSA (2022) IRSA. 2022, Gaia Source Catalogue DR3, IPAC, doi: 10.26131/IRSA544
  * Izidoro et al. (2017) Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750, doi: 10.1093/mnras/stx1232
  * Jiang et al. (2020) Jiang, C.-F., Xie, J.-W., & Zhou, J.-L. 2020, AJ, 160, 180, doi: 10.3847/1538-3881/abb01b
  * Johnson et al. (2010) Johnson, J. A., Aller, K. M., Howard, A. W., & Crepp, J. R. 2010, PASP, 122, 905, doi: 10.1086/655775
  * Kepler Mission (2019) Kepler Mission. 2019, Kepler Stellar Properties Table, IPAC, doi: 10.26133/NEA6
  * Kepler Project (2020) Kepler Project. 2020, Kepler Mission Completeness and Reliability Products, IPAC, doi: 10.26133/NEA14
  * KETTENRING (1971) KETTENRING, J. R. 1971, Biometrika, 58, 433, doi: 10.1093/biomet/58.3.433
  * Kraus et al. (2016) Kraus, A. L., Ireland, M. J., Huber, D., Mann, A. W., & Dupuy, T. J. 2016, AJ, 152, 8, doi: 10.3847/0004-6256/152/1/8
  * Kruijssen et al. (2020) Kruijssen, J. M. D., Longmore, S. N., & Chevance, M. 2020, ApJ, 905, L18, doi: 10.3847/2041-8213/abccc3
  * Kunimoto & Matthews (2020) Kunimoto, M., & Matthews, J. M. 2020, AJ, 159, 248, doi: 10.3847/1538-3881/ab88b0
  * Lindegren (2018) Lindegren, L. 2018.  http://www.rssd.esa.int/doc_fetch.php?id=3757412
  * Liu et al. (2022) Liu, B., Raymond, S. N., & Jacobson, S. A. 2022, Nature, 604, 643, doi: 10.1038/s41586-022-04535-1
  * Longmore et al. (2021) Longmore, S. N., Chevance, M., & Kruijssen, J. M. D. 2021, ApJ, 911, L16, doi: 10.3847/2041-8213/abeb22
  * Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv e-prints.  https://arxiv.org/abs/1109.2497
  * McTier & Kipping (2019) McTier, M. A. S., & Kipping, D. M. 2019, MNRAS, 489, 2505, doi: 10.1093/mnras/stz2088
  * Moe & Kratter (2021) Moe, M., & Kratter, K. M. 2021, MNRAS, 507, 3593, doi: 10.1093/mnras/stab2328
  * Mulders et al. (2015) Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 798, 112, doi: 10.1088/0004-637X/798/2/112
  * Mulders et al. (2018) Mulders, G. D., Pascucci, I., Apai, D., & Ciesla, F. J. 2018, AJ, 156, 24, doi: 10.3847/1538-3881/aac5ea
  * 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
  * Nesvorný & Morbidelli (2012) Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117, doi: 10.1088/0004-6256/144/4/117
  * Newton et al. (2019) Newton, E. R., Mann, A. W., Tofflemire, B. M., et al. 2019, ApJ, 880, L17, doi: 10.3847/2041-8213/ab2988
  * Pedregosa et al. (2012) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2012, arXiv e-prints, arXiv:1201.0490.  https://arxiv.org/abs/1201.0490
  * Petigura et al. (2018) Petigura, E. A., Marcy, G. W., Winn, J. N., et al. 2018, AJ, 155, 89, doi: 10.3847/1538-3881/aaa54c
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Pu & Wu (2015) Pu, B., & Wu, Y. 2015, ApJ, 807, 44, doi: 10.1088/0004-637X/807/1/44
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249, doi: 10.1007/s10686-014-9383-4
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523, doi: 10.1051/0004-6361:20031117
  * Sanchis-Ojeda et al. (2014) Sanchis-Ojeda, R., Rappaport, S., Winn, J. N., et al. 2014, ApJ, 787, 47, doi: 10.1088/0004-637X/787/1/47
  * Sandoval et al. (2021) Sandoval, A., Contardo, G., & David, T. J. 2021, ApJ, 911, 117, doi: 10.3847/1538-4357/abea9e
  * Santos et al. (2001) Santos, N. C., Israelian, G., & Mayor, M. 2001, A&A, 373, 1019, doi: 10.1051/0004-6361:20010648
  * Strömberg (1946) Strömberg, G. 1946, ApJ, 104, 12, doi: 10.1086/144830
  * Su et al. (2021) Su, X.-N., Xie, J.-W., Zhou, J.-L., & Thebault, P. 2021, AJ, 162, 272, doi: 10.3847/1538-3881/ac2ba3
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38, doi: 10.3847/1538-4365/aab4f9
  * Tian et al. (2015) Tian, H.-J., Liu, C., Carlin, J. L., et al. 2015, ApJ, 809, 145, doi: 10.1088/0004-637X/809/2/145
  * Tsiganis et al. (2005) Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459, doi: 10.1038/nature03539
  * Wang & Fischer (2015) Wang, J., & Fischer, D. A. 2015, AJ, 149, 14, doi: 10.1088/0004-6256/149/1/14
  * Wang et al. (2014) Wang, J., Xie, J.-W., Barclay, T., & Fischer, D. A. 2014, ApJ, 783, 4, doi: 10.1088/0004-637X/783/1/4
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48, doi: 10.3847/1538-3881/aa9ff6
  * Weiss et al. (2021) Weiss, L. M., Dai, F., Huber, D., et al. 2021, AJ, 161, 56, doi: 10.3847/1538-3881/abd409
  * Wielen (1977) Wielen, R. 1977, A&A, 60, 263
  * Williams et al. (2013) Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101, doi: 10.1093/mnras/stt1522
  * Winn & Fabrycky (2015) Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409, doi: 10.1146/annurev-astro-082214-122246
  * Winter et al. (2020) Winter, A. J., Kruijssen, J. M. D., Longmore, S. N., & Chevance, M. 2020, Nature, 586, 528, doi: 10.1038/s41586-020-2800-0
  * Xiang et al. (2017) Xiang, M., Liu, X., Shi, J., et al. 2017, ApJS, 232, 2, doi: 10.3847/1538-4365/aa80e4
  * Xiang et al. (2019) Xiang, M., Ting, Y.-S., Rix, H.-W., et al. 2019, ApJS, 245, 34, doi: 10.3847/1538-4365/ab5364
  * 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
  * Yan et al. (2022) Yan, H., Li, H., Wang, S., et al. 2022, The Innovation, 3, 100224, doi: 10.1016/j.xinn.2022.100224
  * Yang et al. (2020) Yang, J.-Y., Xie, J.-W., & Zhou, J.-L. 2020, AJ, 159, 164, doi: 10.3847/1538-3881/ab7373
  * Zhou et al. (2007) Zhou, J.-L., Lin, D. N. C., & Sun, Y.-S. 2007, ApJ, 666, 423, doi: 10.1086/519918
  * Zhu (2019) Zhu, W. 2019, ApJ, 873, 8, doi: 10.3847/1538-4357/ab0205
  * Zhu & Dong (2021) Zhu, W., & Dong, S. 2021, ARA&A, 59, 291, doi: 10.1146/annurev-astro-112420-020055
  * Zhu et al. (2018) Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101, doi: 10.3847/1538-4357/aac6d5