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

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 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 Kep) is about 50% for all stars; no significant trend is found between F Kep and age. (2) The average planet multiplicity ( N¯p ) exhibits a decreasing trend (∼2σ significance) with age. It decreases from N¯p ∼ 3 for stars younger than 1 Gyr to N¯p ∼ 1.8 for stars of about 8 Gyr. (3) The number of planets per star (η = F Kep× N¯p ) also shows a decreasing trend (∼2σ–3σ significance). It decreases from η ∼ 1.6–1.7 for young stars to η ∼ 1.0 for old stars. (4) The mutual orbital inclination of the planets (σ i,k ) increases from 1.°2−0.5+1.4 to 3.°5−2.3+8.1 as the stars age from 0.5 to 8 Gyr with a best fit of logσi,k=0.2+0.4×logAge1Gyr . Interestingly, the solar system also fits such a trend. The fact that F Kep remains relatively constant at approximately ∼ 50% across different ages suggests the robustness of planet formation throughout the history of the Galaxy. The age dependence of N¯p and σ i,k demonstrates that the planetary architecture is evolving, and planetary systems generally become dynamically hotter with fewer planets as they age.


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 Archive; 8 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 the 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 "Kepler planets"), their occurrence rate in terms of number of planets per star (η) has an inverse relationship with stellar temperature and mass (Howard et al. 2012;Mulders et al. 2015;Kunimoto & Matthews 2020).In fact, η 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 N p ¯), and they are linked by the following equation: studies have shown that both F and 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 the 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) and provides the key evidence for the coreaccretion 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 members 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 of 100 au (Wang et al. 2014;Kraus et al. 2016;Fontanive et al. 2019;Moe & Kratter 2021;Su et al. 2021).Recently, it has been reported that both the period and radius distributions of exoplanets exhibit dependencies on stellar clustering (Kruijssen et al. 2020;Winter 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 the 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. 2021b).
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;David et al. 2021;Sandoval et al. 2021;Chen et al. 2022) have revealed that the relative occurrence (ratio) of super-Earths and sub-Neptunes evolves on a timescale of gigayears, providing crucial constraints on the formation of the radius valley (a deficit of planets with radii of ∼1.7-2.1 R ⊕ ; Fulton et al. 2017).More recently, Bashi & Zucker (2022) found tentative evidence that suggests that the occurrence rate of close-in super-Earths detected by Kepler is anticorrelated with stellar age.However, such an anticorrelation is still inconclusive, probably because they adopted isochrone ages that suffer from large uncertainties (56% for Kepler stars; Berger et al. 2020b).Furthermore, they did not isolate the effect of age from other stellar properties (e.g., metallicity), so it is still unclear whether the anticorrelation 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 papers called Planets Across Space and Time (PAST; Chen et al. 2021a).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. 2021a), we revisited the kinematic method to classify Galactic components and the age-velocity dispersion relation (AVR; Holmberg et al. 2009;Strömberg 1946;Wielen 1977), 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. 2021b).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 offer discussions and summarize the main conclusions in Sections 5 and 6.

Star Sample
The LAMOST-Gaia-Kepler catalog constructed in PAST II is based on LAMOST DR4/DR5 and Gaia DR2.Since LAMOST has updated to DR8 (Yan et al. 2022), and Gaia has released EDR3/DR3 (Gaia Collaboration et al. 2023), we have updated the LAMOST-Gaia-Kepler catalog accordingly.We start from the stellar properties catalog 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 the 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 a σ CDPP value.9Then we cross-match the sample with the recently released LAMOST DR810 lowresolution catalog (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 sets (Xiang et al. 2019).Next, we cross-match with the Gaia DR3 (Gaia Collaboration et al. 2023; IRSA 2022) catalog, which includes more accurate measurements of positions, proper motions, and parallaxes of stars compared to DR2.Gaia DR3 also provides the renormalized 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 stars in total; the number in 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 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 kpc < 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 (dubbed D, TD, H, and Her).We classify stars into different components following the commonly used method introduced by Bensby et al. (2003Bensby et al. ( , 2014) ) and show the results in the Toomre diagram (Figure 1).Since AVR can only be applied to disk stars, we limit our sample to stars within the Galactic disk, i.e., stars with D/Her 2, TD/H 1, and TD/Her 2.
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 evolved stage is calculated by the same method as Bryson et al. (2020) using the python package evolstate. 11We further limit our star sample to mainsequence solar-type stars with effective temperatures between 4700 and 6500 K.So far, the number of stars in our sample is 19,537.

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 of 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 ultrashort-period planets (USPs; period < 1 day) from our fiducial analysis (see Section 5.1 for more discussion of USPs).This exclusion is because 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 of ∼1%) and may have undergone different formation and evolution processes (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 ⊕ 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 ⊕ (see Section 5.1 for more discussion of 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.

Kinematic Age of Planet Host and Nonhost
In PAST II, we showed that stars with different numbers of transiting planets exhibit a drop in their kinematic age as a function of planet multiplicity.Here we also bin stars into four groups according to their number of transiting planets; 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 Median ± 5 × 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. 2021b).As can be seen, the kinematic ages derived in this work are generally consistent with those of   PAST II in 1σ-2σ 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.

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 N 1 , N 2 , and N 3+ ) divided by the number of stars (N star ) in each bin.

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 σ 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 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. 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. 2021a).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.Due to the limited sample size, we do not consider cases with more bins.2. Choosing a standard sample.We dub the stars in the middle bin the "standard sample," and the number of stars in the central bin is denoted as N st .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. Checking parameter control result.We calculate the differences of the 25th, 50th, and 75th percentiles of each controlled parameter for every two bins.If all of the differences are less than the typical errors, we consider these parameters to have been controlled.The typical errors of temperature, mass, and radius are 112 K, 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 σ CDPP for different timescales.We choose a σ CDPP of 3.5 hr because it is the closest to the median transit duration of our planet sample (3.35 hr).Since the signal-to-noise ratio (S/N) of Kepler planets is in proportion to R R p s 2 CDPP ( ) s , we choose a typical error of 10% for σ CDPP to match the uncertainty induced by the stellar radius and planet-star radius ratio.
To get an intuitive view of how well the parameters have been controlled, we plot Figures 5 and 6, in which we perform parameter control for the cases of three and five age bins.In the top row of Figures 5 and 6, we show the cumulative distribution function (CDF) diagrams of temperature, mass,  metallicity, radius, and σ CDPP of 3.5 hr for the observation data.Using the above method, we control all parameters and show the CDF of the 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 distributions in stellar temperature, mass, metallicity, radius, and σ CDPP (Figures 5 and 6).

Apparent Planetary Occurrence as a Function of Age
We first consider a three-bin case and calculate the kinematic age of each bin using AVR (Appendix A).We adopt the above controlled sample and calculate the number of stars that have one, two, and three or more planets, i.e., the planet multiplicity function (N 1 , N 2 , and N 3+ ).The apparent occurrence rate of systems with one, two, and three or more planets is derived by dividing N 1 , N 2 , and N 3+ by the number of stars (N star ) in each bin.In Figure 7, in the columns from left to right, we present the apparent occurrence rate for systems with one, two, and three or more planets.As can be seen, the young stars generally have a higher apparent occurrence than the old stars.For the original data without parameter control (top row of Figure 7), N 1 /N star , N 2 /N star , and N 3+ /N star are 2.12 0.17 0.18 -+ %, 0.61 0.09 0.10 -+ %, and 0.35 0.07 0.08 -+ % for stars of less than 1 Gyr, which are about 4.8σ, 4.4σ, and 4.8σ higher than those (1.37 0.13 0.15 -+ %, 0.26 0.06 0.07 -+ %, and 0.08 0.03 0.05 -+ %) for stars of about 8 Gyr, respectively.For the data after all parameter controls (bottom row of Figure 7), N 1 /N star , N 2 /N star , and N 3+ /N star are 1.86 0.22 0.25 -+ %, 0.46 0.11 0.14 -+ %, and 0.31 0.09 0.12 -+ % for the youngest group, which are about 1.6σ, 1.4σ, and 3.3σ higher than those (1.50 0.20 0.22 -+ %, 0.31 0.09 0.12 -+ %, and 0.05 0.03 0.07 -+ %) for the oldest group, respectively.The differences in the apparent rate between young and old groups become smaller when taking parameter control into account.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 given in each panel of Figure 7.The p-values are derived using the following steps.
1. We calculate the Pearson correlation coefficients for the observation data as ρ obs .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 the age and the scrambled data ( sim r ). 3. We repeat step 2 10,000 times and calculate the fraction of the simulated data that produce a stronger anticorrelation, i.e., sim obs r r < . This fraction gives the p-value of the Pearson correlation for the observed data.
As we can see, the anticorrelations between age and N 1 /N star , N 2 /N star become weaker after parameter control, with the pvalues rising from 0.0978 to 0.145 and from 0.0237 to 0.181, respectively.Nevertheless, the anticorrelation between age and N 3+ /N star becomes stronger, with the p-value decreasing from 0.056 to 0.0117.
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-bin 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 %) for the oldest stars, respectively (top row of Figure 8).For the data after all parameter controls, N 1 /N star , N 2 /N star , and N 3+ /N star are 1.65 0.29 0.35 -+ %, 0.57 0.17 0.23 -+ %, and 0.41 0.14 0.20 -+ % for stars of less than 1 Gyr, which are about 0.5σ, 1.1σ, and 2.8σ higher than those (1.50 0.28 0.33 -+ %, 0.36 0.13 0.19 -+ %%, and 0.05 0.04 0.12 -+ %) for stars of about 8 Gyr, respectively (bottom row of Figure 8).Again, being consistent with the results in the three-bin case, the differences in the apparent occurrence rate of stars with different ages become smaller after parameter control.Nevertheless, the differences are still significant (∼3σ) 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-bin case, and list them in each panel of Figure 8. Similar to the three-bin case, the anticorrelation 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 anticorrelations 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 results similar to those shown above, indicating that stellar age is anticorrelated 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, refer to Appendix B.

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 the 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.¯, N 2 ¯, and N 3 ¯+).For our star sample with a size of ∼20,000, we need to generate about ∼10,000 planet systems, assuming that, on average, 50% of the stars own planet systems (the true value of F Kep differs for each group and is automatically adjusted in the Markov Chain Monte Carlo, MCMC, process).The total number of generated planet systems is about 400,000,000 when the simulation is converged.We assume that the multiplicity function follows a Poisson distribution and optimize the likelihood function with python package emcee (Foreman-Mackey et al. 2013) applying the MCMC method.Three free parameters are constrained in our model: the fraction of stars with Keplerlike planets (F Kep ), the average number of planets for stars that have Kepler-like planets (N p ¯), and the inclination slope index (α; see below).The details of generating the modeled multiplicity function (N 1 ¯, N 2 ¯, and N 3 ¯+) can be seen in Yang et al. (2020).We briefly summarize the general procedure as follows.
1. Assuming the intrinsic planet occurrence.For each group of stars, we assume that F Kep percent of stars have Kepler-like planets, and, on average, each planet system has N p ¯planets.For each host star, we generate k planets following a zero-truncated Poisson distribution of N p (Fang & Margot 2012).2. Assuming transit parameters and radii for planets.We generate the debiased distributions of transit parameters (ò, ò = R s /a p , where R s is the stellar radius, and a p is the semimajor axis of the planet) and planet radii (R p ) considering three kinds of bias (Mulders et al. 2018), namely, the transit geometry bias ( f tra ), detection efficiency bias ( f S/N ), and vetting efficiency bias ( f vet ).
For each planet, we assign values of ò and R p that are randomly drawn from the debiased distributions.3. Adjusting period and radius ratios.Planets within the same system tend to have similar period (Fabrycky et al. 2014;Brakensiek & Ragozzine 2016;Jiang et al. 2020) and radius (Ciardi et al. 2013;Weiss et al. 2018) ratios.
To account for this correlation, we adjust the period 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. 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. Assigning orbital inclination to generate transits.We calculate the inclination (I p ) of the planets with respect to the observer by where I represents the inclination of the system invariable plane, f is the phase angle, and i is the inclination of the planet with respect to the invariable plane.Both I and f 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, where we adopt σ i,5 as a Gaussian distribution with a mean of 0°. 8 and a standard deviation of 0°.15 from Zhu et al. (2018).We fit the inclination slope index α as the third parameter.A planet is considered to be a transit if its impact factor is less than 1 ( ). 6. Checking detection and vetting efficiencies.Due to detection and vetting efficiencies, not all of the transiting planets can be detected.For each transiting planet, we generate a random number from a uniform distribution ranging from zero to 1.We consider this planet to be detected if the generated random number is less than the product of the detection efficiency ( f S/N ; see Appendix C) and the vetting ( f vet ) efficiency.Compared to our previous work, we do not consider the transit timing variation (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 than the 100,000 star sample and 127 systems that show TTV signals in Yang et al. (2020).The smaller number for the TTV multiplicity function leads to a 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 α, it has little impact on the results of F Kep and N p ¯, which are the core parameters of this work.

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 left to right, we present the posterior distributions of F Kep , N p ¯, α, and η (which is the product of F Kep and N p ¯; see Equation (1)).For data without parameter control, the youngest group generally has higher intrinsic planet occurrence rates than the oldest group.for stars of less than 1 Gyr, which are 1.7σ, 1.8σ, and 5.4σ higher than those (47.5 8.4 8.9 -+ %, 1.97 0.32 0.41 -+ , and 0.93 0.12 0.13 -+ ) for stars of about 8 Gyr, respectively (top row of Figure 9).For data after all parameter controls, F Kep , N p ¯, and η are 56.9 11.6 ) for the oldest group, respectively (bottom row of Figure 9).All three groups have nearly the same F Kep when parameter control is taken into account.
This result is expected because F 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σ to 1.6σ after parameter control (Figure 7), and, consequently, the difference of F Kep drops from 1.7σ to 0.2σ.
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σ to 3.3σ after parameter control (Figure 7).Interestingly, the difference in N p ¯increases slightly from 1.8σ to 2.1σ.Due to the limited number of planetary systems with three or more planets in the last bin after parameter control (two systems in the three-bin case), the Poisson error is relatively high (2 1.3 1.8 -+ ).We may have underestimated the declining trend of N 3+ /N star , and the results of the forward modeling show that the decrease in N p ¯becomes slightly more prominent.As η is the product of F Kep and N p ¯, therefore, the decrease of the difference in η from 5.4σ to 3.0σ is mainly due to the decrease of the difference in F Kep .We calculate the Pearson correlation coefficients and p-values for the correlations between the age and the intrinsic planet occurrence (F Kep , N p ¯, and η) and list them in each panel in Figure 9.The p-values are derived using a similar method as shown in Section 3.2.The anticorrelations between age and both F Kep and N p ¯are statistically insignificant before and after parameter control with p-values larger than 0.05.The anticorrelation between age and η is maintained with p-values smaller than 0.05.
The parameter α is not well constrained in all cases before and after parameter control because it is mainly constrained by the TTV multiplicity function (Zhu et al. 2018), which is ignored in this work.¯, α, and η 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 Figures 5 and  7.The dots and error bars show the 50% and ±1σ range.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) and corresponding p-value.
To further investigate the intrinsic planet occurrence as a function of age, we group the stars into five bins as mentioned in Section 3.2 and adopt the forward-modeling method to derive F Kep , N p ¯, α, and η.The results are shown in Figure 10.For data without parameter control (top row of Figure 10), the values of F Kep , N p ¯, and η are 68.0 11.9 12.8 -+ %, 2.68 0.41 0.52 -+ , and 1.82 0.20 0.22 -+ , respectively, for stars in the youngest group.These values are 2.1σ, 1.2σ, and 4.8σ higher than those (43.7 8.9 10.9 -+ %, 2.09 0.39 0.54 -+ , and 0.92 0.15 0.17 -+ ) for stars in the oldest group.After parameter control (bottom row of Figure 10), F Kep is 47.1 16.1 17.9 -+ % for stars of less than 1 Gyr, which is 0.6σ lower than for stars of around 8 Gyr (56.6 14.9 for stars in the first bin, which is 2.3σ and 2.2σ higher than the values (1.80 0.39 0.67 -+ and 1.04 0.24 0.29 -+ ) 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σ to 0.5σ (Figure 8).The difference in F Kep decreases from 2.1σ to 0.6σ.The first bin has a higher N 1 /N star compared to the last bin; however, it shows a slightly lower value of F Kep .That is because the first bin has more intrinsic multiplanet 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σ to 2.8σ.Forward modeling indicates that the difference in N p ¯increases slightly from 1.2σ to 2.3σ.Similar to the three-bin case, we have only one planet system with three or more planets in the last bin, resulting in a high Poisson error (1 0.8

-+ ).
As a consequence, we might underestimate the declining trend of N 3+ /N star .As for η (the product of F Kep and N p ¯), the difference drops from 4.8σ to 2.2σ, which is mainly due to the decrease of the difference in F Kep .
The Pearson correlation coefficients and p-values for the correlations between age and intrinsic planet occurrence (F Kep , N p ¯, and η) are also given in the corner of each panel in Figure 10.Regarding the anticorrelation between age and F Kep , it becomes weaker after parameter control, with p-values rising from 0.0503 to 0.729.The anticorrelation between age and N p ¯becomes statistically significant, with p-values dropping from 0.0839 to 0.0052.The anticorrelations between age and η are maintained, with p-values changing from 0.0021 to 0.0096.Similar to the three-bin case, α is not well constrained before and after parameter control.

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 , from which we then derive the intrinsic planet occurrence in terms of F Kep , N p ¯, and η as a function of stellar age using a forwardmodeling 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 a higher apparent occurrence than older stars.Specifically, the intrinsic planet occurrence in terms of the number of planets per star (η) decreases with stellar age with a confidence level of about 2σ-3σ.Such a declining trend is mainly driven by the decrease in the average multiplicity (N p ¯, by about 2σ) and partially by the change in the fraction of stars with planet systems (F Kep , by less than 1σ).In what follows, we will compare our results to those from the literature and discuss the implications of these findings for our understanding of planet formation and evolution.

Giant Planets and USPs
We exclude the giant planets and USPs from our planet sample in Section 2.2.We dub this planet sample the "fiducial" sample.To investigate the influence of giants and USPs on planet occurrence, in this section, we rerun our simulation including the giants, the USPs, and both of them.In Figure 11, from top to bottom, we show the results for the fiducial sample, the sample including giant planets, the sample including USPs, and the sample including both giant planets and USPs for the three-bin case.As we can see, after applying the parameter control method, all results show similar trends.For F Kep , the differences between the youngest and oldest groups are less than 1σ.The youngest groups generally have N p ¯∼ 2σ higher than the oldest groups, and 2σ-3σ higher for η.Including giant planets and USPs adds more planets into our sample, leading to a slightly higher value of F 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 N p ¯. Similar to the three-bin case, the results for the five-bin case are basically unchanged after including giants and USPs.As we can see in Figure 12, the anticorrelations between age and F Kep are statistically insignificant, with p-values higher than 0.5.For N p ¯and η, the p-values are less than 0.05, showing significant anticorrelations between them and age.

Comparison with Previous Studies
McTier & Kipping (2019) 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 nonhost Kepler stars.Based on such a similarity, they inferred that the planet occurrence rate is independent of galactocentric velocity.Their inference is against our results, which show that planet occurrences in terms of N p ¯and η are anticorrelated with kinematic age and thus with galactocentric velocity ¯, α, and η for the five-bin case are presented.The top and bottom rows show the forwardmodeling results corresponding to samples before and after parameter control, as shown in Figures 6 and 8.The dots and error bars show the 50% and ±1σ range.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) and corresponding p-value.
based on 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 (∼50%; Mulders et al. 2018;Yang et al. 2020;He et al. 2021), a large fraction of the apparent "nonhost" 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 nonhost 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., 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. 2021b), Kepler planet hosts are dominated by single-transiting systems, which have a velocity distribution similar to that of nonhost 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 kinematic age seen in this work (Figures 9 and 10) while still maintaining the similarity in galactocentric velocity distribution between the Kepler planet hosts and the nonhost Kepler stars seen in McTier & Kipping (2019).In short, the similarity in galactocentric velocity between Kepler planet hosts and nonhost 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) 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 the 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, is higher in the Galactic thin disk stars than in the thick disk stars.In addition, Bashi & Zucker (2022) also showed an apparent anticorrelation between planet occurrence and the stellar isochrone age.Generally speaking, their results are consistent with ours.In this paper, we find that the 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 a relatively smaller internal uncertainty of ∼10%-20% (Chen et al. 2021a) compared to the 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 anticorrelations are weaker between planet occurrence and age, especially for F Kep , though they remain significant for η (Figures 9 and 10).¯, α, and η for the three-bin case are presented.From top to bottom, the rows show 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 the giants and the USPs.The dots and error bars show the 50% and ±1σ range of the posterior distributions.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) along with the corresponding p-value.¯, α, and η for the five-bin case are shown.From top to bottom, the rows show the forward-modeling results corresponding to the fiducial sample, the sample including the giants, the sample including the USPs, and the sample including both the giants and the USPs.The dots and error bars show the 50% and ±1σ range of the posterior distributions.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) and corresponding p-value.

Implications for Planet Formation and Evolution
In this study, we have revealed observational evidence that the occurrence and architecture (in terms of F Kep and 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 Kep or 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 afterward, 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 number of planets and an increase in orbital spacing (Zhou et al. 2007).This in turn increases the timescale of subsequent instability, which means that the system needs to evolve on a longer timescale to trigger the 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 that 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 event 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 multiplanet systems is 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 compact resonant chains through migration in a protoplanetary gas disk in the early stage.As the gas dissipated, the chains became dynamically unstable, which led to planets merging, ejecting, 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 number of planets per star (η) decreases by about 2σ-3σ as a function of time.For the three-bin case, as shown in Figure 9, η drops from 1.57 to 0.96, and for the fivebin case in Figure 10, it decreases from 1.71 to 1.04.
The first major contribution to the η-decreasing trend comes from the number of planets in the planetary system, i.e., N p ¯, since η is the product of N p ¯and F Kep .N p ¯shows a moderate decline of about ∼2σ in our fitting results.In the bottom row of Figure 9, N p ¯drops from 2.74 for stars of less than 1 Gyr to 1.75 for stars of about 8 Gyr, and in Figure 10, 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 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 η-decreasing trend comes from the fraction of stars that have planets, i.e., F Kep .In Figure 9, F 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σ.Due to the limited star and planet sample, we cannot conclude that the change in F Kep is statistically significant.Future studies with larger samples of planetary systems may help us to further constrain F 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 N p īs 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 σ i,k for the five bins after parameter control and show the distributions, as well as the median value and 1σ range of σ i,k , in Figure 13.As we can see, σ i,k gradually evolves with time.From less than 1 Gyr to about 8 Gyr, the median value of σ i,k grows from about 1°.2 to 3°. 5, and the 1σ range expands from 0°.7-2°.6 to 1°.3-11°.7.To further quantify the age-σ i,k trend, we fit σ i,k with age.Although bearing a large uncertainty (as seen from the orange shaded region), the best fit is shown with green and purple stars, respectively.The age of the solar system is adopted as 4.57 Gyr (Bouvier & Wadhwa 2010), and the age of the Kepler multiples is adopted as the kinematic age of the stars that host two or more planets in our star sample.
For comparison, we also plot the data points for our solar system and Kepler multiple-transiting systems in Figure 13.They all generally fit such an age-σ 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 the planets in our solar system were in a flatter architecture in early times, then gradually evolved to their current state.

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. Younger stars generally show higher apparent planet occurrence rates (more than 3σ) 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 Figures 7 and 8). 2. Applying a parameter control method can effectively reduce the effects caused by other stellar properties, such as effective temperature, mass, metallicity, radius, and σ 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σ, while the difference in N 3+ /N star maintains a confidence level of about 3σ (bottom rows of Figures 7 and 8). 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 Kep ), average planet multiplicity (N p ¯), and number of planets per star (η).For stars without parameter control, we find that the younger stars have higher F Kep and N p ¯than older stars by about 2σ.The difference in η between younger and older stars is more obvious, at about the 5σ level (top rows of Figures 9 and 10).4.After parameter control, the difference in F Kep drops to less than 1σ, hinting that the planetary system occurrence remains at a similar rate throughout the history of the Milky Way.The difference in N p ¯is about 2σ 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 η (the product of F Kep and N p ¯) by about 2σ-3σ than older stars, which is the combined effect caused by the evolution of F Kep and N p ¯(bottom rows of Figures 9 and 10). 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 (σ i,k ) between their planets increases from 1°.2 to 3°. 5, and the ±1σ range of σ i,k expands from 0°.7-2°.6 to 1°.3-11°.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 of 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 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 the planet sample by 1 order of magnitude or even more.With the help of more data, future studies will further refine our measurements and test our results.
providing five astrometric parameters (positions, parallaxes, and proper motions) for 1.468 billion sources.Compared 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 Sections 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 (http://cdsxmatch.u-strasbg.fr).Second, we carry out an angular distance cut of 1 25 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, θ, Z and Galactic rectangular velocities U LSR , V LSR , W 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 e = 8.18 kpc (Gravity Collaboration et al. 2019;Gaia Collaboration et al. 2021) and Z e = 25 pc (Chen et al. 2001).The solar peculiar motions are taken as [U e , V e , W e ] = [9.58,10.52, 7.01] 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 RUWE >1.4 (Lindegren 2018).(2) Parallax precision filter.Following Dehnen & Binney (1998), we remove stars with relative parallax errors larger than 10% 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% parallax precision requirement, this translates to a distance limit of ∼1/ (0.0494/0.1) ∼ 2.0 kpc.We therefore remove stars with distances exceeding this limit.

A.2. Revisiting the Kinematic Method to Classify the Galactic Components
With stellar kinematic age following the same criteria in 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 (Her).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 eight 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 five bins with boundaries at R = 7.5, 8.0, 8.5, 9.0, 9.5, and 10 kpc.In total, there are 5 × 8 = 40 grids in the R-Z space, and all bins have enough (>400) stars.
We then revise the normalized fraction X (Equations 9-12 in PAST I) and the velocity ellipsoid (i.e., σ U , σ V , σ W , and V asym ; Equation (3) in PAST I) of each Galactic component for each grid in the R-Z plane following the same procedure as in Sections 2.3.4 and 2.3.5 of PAST I.The calculated values of X, σ U , σ V , σ W , and V asym are tabulated in Table 2 and visualized in Figures 15-17.
Figure 15 shows the X values of various Galactic components as functions of Galactic radius R and absolute value of height |Z|.As expected, X D (X TD , X H ) generally decreases (increases) with |Z| in all R bins.
With the same procedure as described in Section 2.3.4 of PAST I, we also fit the velocity dispersions, σ U , σ V , and σ W , in the following formula according to Williams et al. (2013): We then use the following formula to calculate V asym according to Robin et al. (2003) and Binney & Tremaine (2008): The values of the fitting parameters and their 1σ uncertainties are summarized in Table 3.
Compared to the results obtained from Gaia DR2 and the 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., σ U , σ V , σ W , and V asym ) obtained with the calibration sample using astrometry data from Gaia DR2 and DR3, as can be seen in Figures 16 and 17, the median values and 1σ error bars are very similar, and the best fits are nearly the same.It can also be seen from Table 4 that 4.
Figure 18 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 the velocity components (U LSR , V LSR , W LSR ) and total velocity (V tot ) are indistinguishable from those of PAST I (red lines).In Table 4, we compare the fitting parameters (k, β) 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σ uncertainties of k decrease by a factor of ∼10% due to the improvement in the precision of the stellar parallax and  Figure 16.The velocity dispersions as functions of position (R, |Z|) in the Galaxy.For the updated calibration data constructed from the LAMOST MSTO-SG and Gaia DR3 catalogs, the velocity dispersions are plotted as filled points, and the line segments represent 1σ 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.For the calibration sample constructed from the LAMOST MSTO-SG and Gaia DR2 catalogs in PAST I, the velocity dispersions are plotted as filled points, and the line segments represent 1σ errors in two colors: light blue for the thin disk and light red for the 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 17.The asymmetric velocity, V asym , as a function of U 2 s for the thin (left panel) and thick (right panel) disks.For the updated calibration data constructed from the LAMOST MSTO-SG and Gaia DR3 catalogs, the data are plotted as blue/red points, and the blue/red line segments represent 1σ errors.The blue/red solid lines denote the results of the best fit using Equation (A2).For the calibration sample constructed from the LAMOST MSTO-SG and Gaia DR2 catalogs in PAST I, the data and 1σ errors are plotted in light blue/red.The light blue/red dashed lines denote the results of the best fit using Equation (A2).
Figure 18.The velocity dispersions for U LSR , V LSR , W LSR , and V 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. proper-motion measurements.Thus, the AVRs derived with the calibration sample using astrometry data from Gaia DR3 and DR2 are nearly the identical to each other.
planets we observed in the planetary system.Given that 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 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 σ CDPP ) instead of all five.
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-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 both radius and σ CDPP results in a reduction in detection efficiency, which lowers the probability of planet detection.Here we choose σ CDPP to present the effect of detection efficiency (choosing temperature or radius leads to a similar result).
In summary, we categorize all of the planet/star properties into three groups.The first group, X 1 = {Tracer(η), Tracer(F p )}, represents planet occurrence.The second group, X 2 = {Age}, describes the interesting stellar properties, specifically, stellar kinematic age.The third group, X 3 = {Mass, [Fe/H], σ 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: where y 1 = X 1 a 1 , y 2 = X 2 a 2 , and y 3 = X 3 a 3 , and a 1 , a 2 , and a 3 are the weight vectors for each variable set.All of the variables have been standardized.The results can be seen in the following figure.
In Figure 14, each ellipse represents a group, and each box represents a star/planet property.We also give the weight and correlation on each line.As we can see, both Tracer(η) and Tracer(F p ) have positive contributions to the planet occurrence (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 anticorrelated with planet occurrence (Cor = −0.802),which is consistent with our result before parameter control (see top rows of Figures 7 and 8).As stars age from less than 1 Gyr to about 8 Gyr, both 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 the 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., Wang & Fischer 2015;Zhu 2019), indicating that an increase in metallicity can stimulate the formation of planets.σ CDPP demonstrates a negative weight because higher σ 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 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-bin cases in Figures 19 and  20.The detection efficiency metrics are calculated by the package KeplerPORTs (Burke & Catanzarite 2017), and the associated data are downloaded from the NASA Exoplanet Archive. 12s we can see in the top rows of Figures 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 (σ CDPP ) that increase the S/N (top rows of Figures 5 and 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 σ CDPP (bottom rows of Figures 5 and 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 Figures 19 and 20), showing that the influence caused by detection efficiencies on planet occurrence is effectively removed.

Figure 1 .
Figure 1.Toomre diagram of stars in the updated LAMOST-Gaia-Kepler sample.The blue, orange, gray, green, and red dots present stars in the thin disk, thick disk, in between thin and thick disk, halo, and Hercules stream, respectively.The gray dashed lines represent the total Galactic velocity at 100, 200, and 300 km s −1 .

Figure 2 .
Figure 2. Hertzsprung-Russell diagram of stars in the updated LAMOST-Gaia-Kepler sample.The blue, orange, and red dots show stars in the mainsequence stage, subgiants, and red giants, respectively.

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

Figure 4 .
Figure 4.The blue dots and error bars show the kinematic age and ±1σ ranges for stars with zero, one, two, and three or more planets for this work, and the orange dots and error bars are values from PAST II (Chen et al. 2021b).The xaxis is offset by 0.1 for clearance.

Figure 5 .
Figure5.CDF diagrams of effective temperature, mass, metallicity, radius, and σ CDPP of 3.5 hr for the three-bin method before and after parameter control.The error bar 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 7 .
Figure7.Apparent planet occurrence rate for the three-bin 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 Figure5.The dots and error bars represent the median value and ±1σ 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 give the Pearson correlation coefficient (ρ) and the corresponding p-value.
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 (N 1

Figure 8 .
Figure8.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 Figure6.The dots and error bars represent the median value and ±1σ 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 give the Pearson correlation coefficient (ρ) and corresponding p-value.

F
Kep , N p ¯, and η are 63

Figure 9 .
Figure 9. Posterior distributions of F Kep , N p¯, α, and η 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 Figures5 and 7.The dots and error bars show the 50% and ±1σ range.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) and corresponding p-value.

Figure 10 .
Figure 10.Similar to Figure 9, posterior distributions of F Kep , N p¯, α, and η for the five-bin case are presented.The top and bottom rows show the forwardmodeling results corresponding to samples before and after parameter control, as shown in Figures6 and 8.The dots and error bars show the 50% and ±1σ range.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) and corresponding p-value.

Figure 11 .
Figure 11.Posterior distributions of F Kep , N p¯, α, and η for the three-bin case are presented.From top to bottom, the rows show 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 the giants and the USPs.The dots and error bars show the 50% and ±1σ range of the posterior distributions.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) along with the corresponding p-value.

Figure 12 .
Figure 12.Similar to Figure 11, the posterior distributions of F Kep , N p¯, α, and η for the five-bin case are shown.From top to bottom, the rows show the forward-modeling results corresponding to the fiducial sample, the sample including the giants, the sample including the USPs, and the sample including both the giants and the USPs.The dots and error bars show the 50% and ±1σ range of the posterior distributions.In the upper right corner of the panels in the first, second, and fourth columns, we also give the Pearson correlation coefficient (ρ) and corresponding p-value.

13 .
Inclination dispersion (σ i,k ) as a function of age.We show the distribution of σ i,k for five age bins after parameter control with blue bars.The median value and the ±1σ percentile range of σ i,k are indicated by orange dots and error bars, with the corresponding ±1σ percentile range 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 σ i,k 10,000 times.The σ i,k value of the solar system terrestrial planets is shown with a red star, and the σ i,k values of Kepler multipletransiting systems measured by Fabrycky et al. (2014) and Xie et al. (2016) are 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σ error bars.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 AVRAccording to PAST I, we divide the calibration sample into 30 bins with approximately equal sizes (∼4475 stars in each bin) according to their ages.Then we fit the AVRs followingHolmberg et al. (2009) andAumer et al. (2016) by using a simple power-law formula, i.e., t represents the stellar age, σ is the velocity dispersion, and k and β are two fitting parameters.The best fits and uncertainties (1σ interval) of the fitting parameters (k, β) are calculated with the same procedure as described in Section 3.1 of PAST I and summarized in Table

Figure 14 .
Figure14.The results of CCA for planet and stellar properties are provided, including the weights assigned to each property and the correlations between properties and groups.

Figure 15 .
Figure 15.The normalization fraction X of stars for each component as functions of Galactic radius R and absolute value of height |Z|.The different colors denote subsamples of stars with different Galactic radii.

Figure 19 .
Figure19.Detection efficiencies and planet samples in the period-radius diagram for the three-bin case.From top to bottom, each row corresponds to the control samples in Figures5 and 7.The red lines present the average 90%, 50%, and 10% detection efficiencies for the star sample in each bin, and the gray lines show the mean 90%, 50%, and 10% detection efficiencies for the whole star sample.The blue, orange, and green dots show planets in systems of one, two, and three or more planets, respectively.

Figure 20 .
Figure 20.Similar to Figure 19, 90%, 50%, and 10% detection efficiencies and planet samples for the five-bin case.From top to bottom, each row corresponds to control samples in Figures 6 and 8.
Figure 20.Similar to Figure 19, 90%, 50%, and 10% detection efficiencies and planet samples for the five-bin case.From top to bottom, each row corresponds to control samples in Figures 6 and 8.

Table 2
Revised Characteristics at Different Galactic Radii (R) and Heights (Z) for Different Galactic Components Using the Calibration Sample Updated with Gaia DR3