1 Introduction↩︎

The canonical story of planet formation is that of a localized, isolated process [1]–[3], largely independent of galactic-scale phenomena. Within this picture, to the extent that the galactic context matters, it is through the steady enrichment of metals in the interstellar medium (ISM; [4]). Stellar metallicity traces the metallicity and mass of the protoplanetary disk [5]. This, in turn, is thought to determine planet radii and system architectures [6]–[13].

However, recent observational results challenge the notion that metallicity alone can explain variations in planet occurrence on galactic scales. The extent to which metallicity is itself deterministic, or whether it is a tracer for other processes that shape planet formation, is under active debate. Establishing the relationship between host star metallicity and planet occurrence is itself complicated. For example, metallicity and close binarity [14], [15] both affect planet outcomes, but they are related to one another (see review by [16]) and observationally entangled [17]. This has led some studies to conclude that the masses of planets are determined by factors other than the availability of solids per se [18], instead regulated by an as-yet unknown process. Stellar mass, effective temperature, and age are similar quantities in the planet formation story. They are related to planet outcomes (see e.g. [19]–[24]) but also related to one another. Disambiguating these effects poses a major challenge to understanding occurrence. For example, [25] observed that among Kepler stars, more massive stars are more metal-rich. However, this reflects a selection bias whereby more massive stars are likelier to be younger, and thus formed from more recent and enriched galactic material (see also [18]). On top of these quantities that trace galactic star formation history, planetary systems may themselves evolve over time, which can potentially masquerade as a metallicity-dependent effect. Planetary systems could potentially be self-disrupting [26], [27]. Alternatively, they may change their configuration due to stellar flybys or galatic tides: these can either directly impact outer planets and potentially propagate inward [28]–[30], or affect a binary companion to the host star, whose altered orbit might then disrupt planets [31], [32].

The Gaia mission [33] has spurred a transformation in exoplanetary studies, driven by an improved understanding of host star properties and also their location and movement within the galaxy. Recent studies have begun placing the relationship between metallicity and planet demographics in the context of the host stars’ position in the Milky Way and dynamical history [4], [34]–[38]. Within this new galactic framework for exoplanets, it is necessary to fold together our understanding of exoplanet demographics with a picture of how stars orbit the galactic potential. These orbits can change significantly over the stellar lifetime. Changes occur as stars repeatedly experience gravitational interactions with inhomogeneities in the galactic disk, such as giant molecular clouds (GMCs) [39]–[41], spiral arms [42]–[44], and the bar [45], [46], torques from misaligned stellar and gas disks [47], [48], feedback-driven fluctuations in the potential [49], as well as effects from the cosmological environment such as satellite interactions and mergers [50]–[53]. Interactions with lumps in the disk mass distribution, such as GMCs, can convert in-plane motions to vertical excursions from the plane of the disk [41], [54]–[56], and satellite bombardment can significantly increase the vertical velocity dispersion of the disk [53]. Additionally, there is growing evidence that stellar populations formed during the early stages of the disk (\({\sim}8-4\) Gyr ago in the Milky Way) were born with higher vertical velocity dispersion than in the current epoch [57]–[59]. Broadly speaking, older stellar populations are kinematically warmer than younger populations, an expectation reinforced by findings directly connecting Gaia kinematic data to stellar age (see e.g. [60]–[66]).

Situated at a key junction in this landscape is a recent result by [67] that found that planet occurrence is higher among stars with lower oscillation amplitude from the galactic midplane: it decreases by a factor of several between 100-1000 pc (dependent upon planet size). Given the known relation between stellar metallicity and planet outcomes, a natural explanation is that the negative metallicity gradient with galactic scale height is responsible. This gradient – which itself has been observed as far back as [68] and continues to be further constrained today [69]–[73] – is not sufficient to explain the decrease in planet occurrence with increasing height from the midplane of the Milky Way, at least for small planets. Using a calculated raw exoplanet occurrence among Kepler Super-Earths and Sub-Neptunes (planets with radius <4 \(R_{\oplus}\)) with a period of 1-40 days, as well as Gaia stellar galactic oscillation amplitudes, [67] fit a power law to model the slope of this relation, finding that the reduction in planet occurrence over 1 kpc is significantly greater than the expected reduction from a solely metallicity-driven occurrence-scale height trend. There is some tension between this finding and other studies, such as [74], who concluded that the planet formation rate of Neptune-sized and smaller planets increases slightly and gradually with age, in a way that is adequately explained by galactic metallicity evolution.

The authors of [67] (hereafter Z23) comment that the oscillation amplitude (hereafter denoted as \(Z_{\textrm{max}}\)) could be a proxy for another parameter like stellar age; this is suggestive of a phenomenological experiment, in which we are agnostic about the driver of this occurrence-height trend but can at least narrow down suspects by constraining the possible timescales on which such a driver may operate. Our aim in this study is to test whether and how a fiducial boost in the planet occurrence rate sometime in our Galaxy’s past could have produced such a relation among close-in small planets around Sun-like stars.

In Section 2, we describe our stellar sample (Section 2.1), an idealized, independently synthesized stellar sample (Section 2.1.3), our planet host fraction evolution models (Section 2.2), the way in which we draw our synthetic planetary system populations (Section 2.3), and our completeness calculations (Section 2.4). In Section 3, we present the detection yields for each model and highlight those that come close to matching the Z23 relation between planet occurrence and galactic scale height. In Section 4, we discuss the physical processes that might be represented by these favored models. Finally, we conclude in Section 5.

2 Methods↩︎

We aim to test whether and how a time-dependent planet occurrence rate, when folded together with typical kinematic heating of stars over Gyr, can reproduce the Z23 finding. To this end, we simulate suites of synthetic planetary system populations based on a prescription for the planet host fraction, \(f\), in the Milky Way as a function of stellar age. These models vary \(f\)over time in ways we describe in this section. With an additional prescription for the kinematic heating of stellar orbits in the galaxy, we can synthesize “present-day" planet occurrence yields that should vary with the host stars’ height above the galactic midplane. To be clear, the exact”height" quantity employed by Z23 for the host stars is not the instantaneous position above the midplane, but rather \(Z_{\textrm{max}}\), the maximum vertical oscillation amplitude above and below the galactic midplane. In this study, we use the inverse detection efficiency method (IDEM), in which the detection efficiency is pre-computed and then applied to a detected yield in order to estimate a “true" yield (see Section 2.4).

To streamline and generalize this workflow across different occurrence scenarios, we developed a planetary system population synthesis code in Python called psps¹. The code enables rapid generation of planetary systems drawn probabilistically from user-defined prescriptions for various exoplanet demographic attributes. In Section 2.1, we describe how we construct the stellar samples in this work. In Section 2.2, we define the time-dependent occurrence models used to determine which stars host planets. We then populate those systems using psps (Section 2.3) and “observe" them with geometric and instrumental completeness models that reflect Kepler’s detection sensitivity (Section 2.4).

2.1 Stellar Sample↩︎

The Z23 result was derived from a joint accounting of Kepler and K2 stars. As such, in order to reproduce their result, we must begin with a stellar sample that uniformly includes both surveys. We construct this sample in two ways. First, we populate our sample with actual Kepler and K2 stars, using their isochrone ages and computing their \(Z_{\textrm{max}}\) based on kinematic information from Gaia and a model of the Galactic potential. This approach has the benefit of matching the true host star sample exactly. However, it inherits the significant observational uncertainties associated with stellar age measurements from isochrones. For this first sample, we use the uniform catalog from [75]. We describe this sample in Section 2.1.1.

Second, we generate a synthetic Kepler and K2 stellar sample using TRILEGAL [76], a population synthesis code widely used in exoplanet demographic studies (e.g., [77]–[80]). While the TRILEGAL sample is not drawn from the Kepler and Ecliptic Plane Input Catalogs, it can be tuned to resemble them reasonably well [81], and it offers the advantage of generating arbitrarily large samples with small stellar age uncertainties. This allows us to apply time-dependent planet occurrence models with greater clarity. We describe this synthetic sample in Section 2.1.3.

2.1.1 Real Kepler and K2 sample↩︎

We begin by generating stellar samples from homogeneously derived Kepler and K2 catalogs, culling them in much the same manner as Z23 before enriching them with isochrone ages and \(Z_{\textrm{max}}\). We use the homogeneously derived joint Kepler-K2 catalog from [75] (hereafter HU25), which already precludes likely astrometric binaries via a RUWE cut at 1.4 and noisy stars via CDPP cuts (only Kepler stars with CDPP\(_{7.5hr}\) < 1000 ppm and K2 stars with CDPP\(_{8hr}\) < 1200 ppm are kept). From this catalog, we replicate the cuts made by Z23 and retain only stars that satisfy the following conditions:

1. 4000 K \(\leq\) T_eff\(\leq\) 6500 K

2. log g\(\geq\) \(log \textit{g}_{H16}\)

where H16 is the threshold from the heuristic devised by [82]: \[\label{eq:giant-cut} \frac{arctan(\frac{6300 K-\textit{T}\textsubscript{eff}}{67.172 K})}{4.671} + 3.876.\tag{1}\]

For T_eff, log g, and kinematic information, we use DR3 parameters from HU25 for both Kepler and K2 samples. For the Kepler and K2 isochrone ages, we draw from the Kepler-Gaia cross-match by [83] (hereafter B20) and the K2 portion of the Kepler-K2-TESS-Gaia cross-match by [84] (hereafter B25), respectively. From B20 and B25, we first remove stars with uninformative isochrone age posteriors (‘unReAgeFlag’ in B20 and ‘\(f_{Age}\)’ in B25) and whose best-fit isochrone ages are less than 0.5 Gyr or greater than 8 Gyr. Next, we apply the same T_effand log gcuts to B20 and B25 as we did to HU25 (and as did in Z23). We then cross-match the HU25 Kepler catalog with B20 and the HU25 K2 catalog with the K2 subset of B25. This data preparation step introduces two possible sources of bias into our final sample: 1) the B25 sample contains only planet hosts; and 2) stellar ages from B20 are derived from Gaia DR2 instead of DR3.

To address the first potential issue, in Figure 1, we verify that the age distributions of Kepler stars are essentially the same between hosts and non-hosts. Comparing the cumulative distribution functions (CDFs) of the B20 host ages and the B20 non-host (also known as “field star") ages, we find a Kolmogorov-Smirnov (KS) test statistic of 0.05 with a p-value of 7e-4, which indicates a very high probability that the two sets of ages are drawn from approximately the same population. We deem this result sufficient to move forward using B25 for the K2 ages.

To address the second potential issue, we visually inspect the component and cross-matched catalogs’ Hertzsprung-Russell (HR) and Kiel diagrams to verify that they occupy the same regions of relevant parameter space. The left panels of Figure 2 show that while B20 does not extend to the cooler stars contained in HU25, it more importantly does not leave the footprint of the HU25 sample in T_eff, log g, and luminosity space, indicating that their cross-match does not introduce a bias in stellar temperature, surface gravity, or luminosity. This is less of a concern for K2, since B25 and HU25 both use Gaia DR3, but we can still confirm in the right panels of Figure 2 that B25 and HU25 are self-consistent in their HR and Kiel diagrams. The HR diagrams show the number of stars in each catalog after the T_effand log g(and, for B20 and B25, age) cuts. These numbers, as well as the number of stars remaining at each intermediate parameter space cut, are also listed in Table 1.

In order to calculate the maximum oscillation amplitude \(Z_{\textrm{max}}\)of a star, we also require its parallax, proper motion, and radial velocity, as well as the uncertainties for those parameters. We thus make one final cut to our sample, keeping only those stars with full astrometric solutions. For K2, this removes no stars from our HU25-B25 cross-match of 141 stars; for Kepler, this takes our cross-matched HU25-B20 sample down from 50675 to 22043 (see Table 1). In total, our final sample for this study comprises 22184 stars from Kepler and K2. In the top panel of Figure 3, we show a Kiel diagram of these stars, colored by their isochrone ages, along with marginalized histograms of the T_effand log gdistributions. The median T_effin our sample is 5703 K (with a median absolute deviation, or MAD, of 215 K) and the median log gis 4.356 dex (with a MAD of 0.089 dex). Our sample differs from that of Z23 primarily in our cuts on uninformative age posteriors, which removed most main sequence stars between 4000 and 5000 K. This is not surprising due to the difficulty in obtaining isochrone age measurements for cooler main sequence stars. This difference partially informed our addition to this study of the TRILEGAL stellar sample, which does contain cooler dwarf stars (see Section 2.1.3).

2.1.2 Age-\(Z_{\textrm{max}}\)relation↩︎

We then compute the vertical oscillation amplitude for each star in our curated sample using the software package Gala [85], [86]. Gala reconstructs galactic orbits from the instantaneous Gaia kinematic information together with a model for the Milky Way galactic potential. We run Gala using a Galactic potential model of a spherical nucleus and bulge, a Miyamoto-Nagai disk, and a spherical Navarro-Frenk-White dark matter halo [87], [88]. From these initial inputs and Gaia DR3 parallax, proper motion, and radial velocity information, we integrate the orbits of each star in the sample for 250 Myr (approximately one Solar orbit around the Milky Way) with a timestep of 0.5 Myr before retrieving \(Z_{\textrm{max}}\).

For a time-dependent planet formation history to produce a change in planet occurrence over \(Z_{\textrm{max}}\), the stellar sample needs to exhibit a discernible trend between age and vertical oscillation amplitude. In the bottom panel of Figure 3, we show that using isochrone ages from our joint B20-B25-HU25 sample presents a soft positive trend consistent with previous investigations of \(Z_{\textrm{max}}\)for Kepler stars [62], [66], [89], [90].

From this seed sample of 22184 Kepler and K2 stars with ages and Galactic scale heights, we re-draw their ages and \(Z_{\textrm{max}}\), in order to ensure that we properly account for the large stellar age uncertainties, along with their T_eff, log g, stellar radii, and stellar mass. We use the Gaia DR3 isochrone-derived parameter uncertainties for each draw. We draw 30 stellar samples, around which we will synthesize 30 distinct planet populations (see Section 2.3).

2.1.3 TRILEGAL Kepler and K2 sample↩︎

We next construct a synthetic joint Kepler-K2 sample using TRILEGAL [76], a population synthesis code for Milky Way stars. For our purposes, the TRILEGAL sample allows us to consider an optimistic scenario in which stellar ages are known precisely. To first order, a simulated TRILEGAL pointing of the Kepler footprint furnishes a good match to measured ages from gyrochronology – [80] found that a two-step star formation rate (SFR) prescription, together with TRILEGAL’s built-in kinematic heating model over time, is favored over a constant SFR and no kinematic heating, at least for a sample \(<4\) Gyr. This motivates our choice of the same two-step SFR, in which the modeled thin disk produces a factor of 1.5 more stars with ages between 1–4 Gyr [76].

We compose the joint sample by synthesizing the Kepler and K2 fields separately and then pieceing them together. To generate the TRILEGAL “Kepler-like" stellar sample, we query TRILEGAL² with input pointing parameters of l=76.32 deg, b=13.5 deg, and field area of 10 \(\textrm{deg}^2\). We employ the Chabrier lognormal initial mass function (IMF), binaries toggled on with a binary fraction of 0.5 (rather than the default of 0.3, based on the more recent results of [91]) and mass ratios of 0.7 to 1, default extinction parameters, a solar position of R=8122 pc and z=20.8 pc, a thin disk population generated using a two-step SFR (with otherwise default thin-disk parameters), a thick disk population generated using default parameters, and no halo or bulge population. TRILEGAL’s prescription for kinematic heating includes an evolving scale height, with \(h_z\) = \(z_0\)(1+\(t\)/\(t_0\))^\(\alpha\), where \(z_0\)=94.7 pc, \(t_0\)=5.55 Gyr, and \(\alpha\) = 5/3).

Since TRILEGAL parameters are generated on a grid, we perturb relevant stellar parameters (specifically, age and distance) with a spread equal to the bin size, following the procedure from [80]. The bin sizes for log(age) and distance modulus are 0.02 and 0.05, respectively. This means that the average assigned age uncertainty for the TRILEGAL stellar sample is 1.0 Gyr; for comparison, the average upper and lower age uncertainties for the HU25-B20 Kepler sample are 2.8 Gyr and 1.9 Gyr, respectively, while the average upper and lower age uncertainties for the HU25-B25 K2 sample are 5.6 Gyr and 3.4 Gyr. We build our initial TRILEGAL Kepler sample by running the web tool ten times using the prescription described above, resulting in an initial sample size of 280,647 stars. This was subsequently cut to 132,951 non-binary stars, 101,635 stars with T_effbetween the maximum (6499 K) and minimum (4679 K) values in the combined Kepler-K2 sample, 67,408 stars after the Equation 1 giant removal procedure was applied, 82177 stars between 4.0 and 4.7 in log g, and 28,792 stars within the Kepler magnitude range of HU25-B20-B25 (7.40-15.07).

Unlike Kepler, the K2 mission was comprised of 19 separate campaigns (omitting Campaign 0) with different pointings. We make separate TRILEGAL queries for each of the 19 fields, two times each. Each query uses the same input parameters as the TRILEGAL Kepler sample, except for the Galactic coordinates, which vary per campaign and are taken from the Data Release Notes on MAST³. In Figure 4, we show the locations in Galactic latitude and longitude of the search fields for the Kepler and K2 TRILEGAL samples, overlaid on the Galactic coordinates of our HU25-B20-B25 sample. The queries over the 19 campaigns returned 842,457 stars.

From this initial set, we cut the sample to 496,124 non-binary stars, 131,346 stars with T_effbetween the maximum and minimum values of the HU25-B20-B25 sample, 61,736 stars after applying the giant cut with Equation 1 , and 20,317 stars with Kepler magnitude within the HU25-B20-B25 sample’s corresponding range. In Figure 5, we compare the T_eff, log g, Kepler magnitude, age, and instantaneous height (or \(Z_{\textrm{max}}\), where applicable) distributions of our TRILEGAL and HU25-B20-B25 samples. Notably, HU25-B20-B25 peaks approximately 350 K cooler than TRILEGAL, and we attribute this to the uninformative age posterior cuts described above. The TRILEGAL log gdistribution is also wider, with more puffy stars that are likely sub-giants and early post-main sequence giants, as well as more high-log gstars in the main sequence. In age (and Galactic height), HU25-B20-B25 also peaks slightly later (and higher). In total, across the Kepler and K2 field queries, our final TRILEGAL sample contains 49,109 stars.

A Kiel diagram for the TRILEGAL sample is shown in top panel of Figure 6. There are two main differences between this sample and our HU25-B20-B25 based on real Kepler and K2 stars. First, the TRILEGAL sample extends well into the K dwarf regime. Second, it omits likely sub-giants in the 5000 < T_eff< 5500 K and 4.2 < log g dex range.

Unlike the HU25-B20-B25 sample, the TRILEGAL stellar sample does not have astrometric parameters with which to calculate \(Z_{\textrm{max}}\). Instead, we are given the distance modulus, which, together with the inclination of the Kepler field, can be used to calculate each star’s instantaneous “height" from the midplane. In bottom panel of Figure 6, we show that the TRILEGAL stellar sample exhibits a positive relation between age and this height, similar to our”real" sample (bottom panel of Figure 3). By definition, one should expect the TRILEGAL instantaneous heights to be lower than the B20 \(Z_{\textrm{max}}\)values; this is borne out in our samples, although only slightly, with the TRILEGAL median height being 244 pc (MAD of 107 pc) and the HU25-B20-B25 median \(Z_{\textrm{max}}\)being 263 pc (MAD of 85 pc).

2.2 Time-dependent occurrence model↩︎

We employ two families of functions to model the planet occurrence rate history of the Milky Way. In both prescriptions, we assume that the Milky Way is 13.7 Gyr old [65], [92], and that planet formation has proceeded during that time in multiple phases. An acceptable model must also reproduce the present-day planet-host fraction of \(\sim\)​30% [93] of stars, though the path to arriving at that fraction can vary. Given that we are attempting to mimic a late-time rise in planet occurrence consistent with a higher midplane fraction, planet occurrence increases monotonically as the galaxy ages in all models.

The first family of models we explore are step functions with three parameters:

• fraction of stars hosting a planetary system before a planet formation-boosting threshold, \(f_1\)

• fraction of stars hosting a planetary system after the planet formation-boosting threshold, \(f_2\)

• lookback time at which \(f_1\)increases to \(f_2\), \(t\).

A model of {\(f_1\)=0.20, \(f_2\)=0.95, and \(t\)=1.7}, for example, means 20% of stars hosted planets up until \(\sim\)​1.7 Gyr ago (that is, in the first \(\sim\)​12 Gyr of the Milky Way’s lifetime). After that point, 95% of stars host a planetary system. In Figure 7, we show seven different step-function models, with \(t\)ranging from 1.7 to 8.2 Gyr ago and \(f_1\)and \(f_2\)chosen to yield an overall present-day planet host fraction of 30-35% [93], as well as a control (flat) model set constant at \(f\)=33%. The step function models are as follows:

1. \(f_1\)=20%, \(f_2\)=95%, \(t\)=1.7 Gyr ago

2. \(f_1\)=10%, \(f_2\)=100%, \(t\)=2.2 Gyr ago

3. \(f_1\)=1%, \(f_2\)=95%, \(t\)=2.7 Gyr ago

4. \(f_1\)=1%, \(f_2\)=75%, \(t\)=3.2 Gyr ago

5. \(f_1\)=5%, \(f_2\)=50%, \(t\)=4.2 Gyr ago

6. \(f_1\)=1%, \(f_2\)=40%, \(t\)=6.2 Gyr ago

7. \(f_1\)=1%, \(f_2\)=35%, \(t\)=8.2 Gyr ago

In the second family of models for planet occurrence with time, rather than an instantaneous jump between two constants, we employ a piecewise function. Planet occurrence begins with an initially constant planet-host fraction, \(f_1\), followed by a positively sloped linear function of \(f\)with cosmic time (or, equivalently, a negatively sloped function of planet occurrence with lookback time or stellar age), up to a constant, present-day planet-host fraction, \(f_2\). We choose the piecewise functions such that they increase until the present-day in order to construct a model archetype that is sufficiently different from the step function. As with the step function models, these piecewise models are chosen such that the total planet host fraction for stars of all ages comes to 30-35%. The piecewise models employed in this study are as follows and are illustrated in orange in Figure 7:

1. \(f_1\)=15%, \(f_2\)=100%, \(t\)=3.2 Gyr ago

2. \(f_1\)=5%, \(f_2\)=100%, \(t\)=4.2 Gyr ago

3. \(f_1\)=1%, \(f_2\)=95%, \(t\)=5.2 Gyr ago

4. \(f_1\)=5%, \(f_2\)=70%, \(t\)=6.2 Gyr ago

5. \(f_1\)=1%, \(f_2\)=65%, \(t\)=7.2 Gyr ago

6. \(f_1\)=1%, \(f_2\)=60%, \(t\)=8.2 Gyr ago,

For both model archetypes, model selection began with a wide range of \(t\)and iterated to the set presented in this study. Each \(t\)is degenerate in that there is a set of \(f_1\)and \(f_2\)that lead to a present-day \(f\)of 30-35%. As we will later see in Section 3, however, so much of the model space is ruled out due to a too-shallow occurrence-\(Z_{\textrm{max}}\)trend. Therefore, we can rule out entire \(t\)by testing the two solutions per \(t\)that generate the theoretically steepest trend – one with \(f_1\)=1% and one with \(f_2\)=100% – and finding that both are still too shallow.

2.3 Planetary System Population Synthesis↩︎

We now have in hand a stellar sample, spanning a variety of ages, with some established fraction hosting planets according to the prescription described in Section 2.2. We have already established that acceptable models must reproduce the present-day planet host fraction of \(\sim\)​30-35%. In order to maintain consistency with observed transit multiplicity as well, we adopt a heuristic for planetary system architecture. There is theoretical support for planetary systems exhibiting warmer dynamical temperatures with stellar age [26], [94], but [27] found no distinction over Gyr timescales between dynamically hot and cold systems from the same B20 sample. They concluded rather that \(18^{+15}_{-10}\)% of planet-hosting Kepler FGK dwarfs have a “system of tightly-packed inner planets" (STIPs) characterized by low angular momentum deficit [95], but with no evidence that the fraction changes with time. We designate systems as either dynamically”hot" or “cold" via their system properties such as mutual inclination, orbital eccentricity, and planet multiplicity [95]–[97]. We hold the relative hot/cold fraction to be constant, varying only the fraction of stars that host a planetary system. In other words, stellar age determines only the fraction of planet hosts in the sample and not the number of planets per system, nor their transit geometries.

We determine the fraction of dynamically “cool" systems in each stellar sample by drawing from a normal distribution, \(\sim\) \(N\)(0.18, 0.1), which approximates the corresponding posterior distribution from [27]. Planet-hosting systems designated dynamically cool are randomly assigned either 5 or 6 planets, with inclinations from the midplane drawn from \(N\)(\(\mu\)=0, \(\sigma\)=2) [27], [93], [95], [98]–[100]. Systems that are dynamically hot are randomly assigned either 1 or 2 planets, with inclinations from the midplane drawn from \(N\)(\(\mu\)=0, \(\sigma\)=8). Planet eccentricities are drawn from a Rayleigh distribution with a peak at 0.24 for singles and a peak at 0.06 for multis [101]. All planet hosts are assigned a midplane orientation randomly drawn from \(U\)(-\(\pi\)/2, \(\pi\)/2). All planets are assigned a longitude of periastron that is also randomly drawn from \(U\)(-\(\pi\)/2, \(\pi\)/2).

Planets are assigned equal probability of being either a Super-Earth (1.2 \(R_{\oplus}\) < \(R_{p}\) < 2 \(R_{\oplus}\)) or Sub-Neptune (2 \(R_{\oplus}\) < \(R_{p}\) < 4 \(R_{\oplus}\)). As in Z23, the radii are assigned following a power law: \[q(R_p) = R_p^\alpha, \label{eq:p95intact3}\tag{2}\] where \(R_{p}\) is the planet radius and \(\alpha\) is drawn once per stellar population from \(N\)(\(\mu\)=-1, \(\sigma\)=0.2) if the planet is a Super-Earth (1.2 \(R_{\oplus}\) < \(R_{p}\) < 2 \(R_{\oplus}\)) and from \(N\)(\(\mu\)=-1.5, \(\sigma\)=0.1) if the planet is a Sub-Neptune (2 \(R_{\oplus}\) < \(R_{p}\) < 4 \(R_{\oplus}\)). We then draw planet periods on either side of the radius valley as parameterized by [102], which found a slope \(\gamma\) of -0.09\(^{+0.02}_{-0.04}\), upper envelope y-intercept \(a_{upp}\) of 0.44\(^{+0.04}_{-0.03}\), and lower envelope y-intercept \(a_{low}\) of 0.29\(^{+0.04}_{-0.03}\). Planets are given a 5% chance of falling within the valley. For each of our 30 stellar population realizations, we draw \(\gamma\), \(a_{upp}\), and \(a_{low}\) from asymmetric distributions following the approximate method in [103] and [104]. Masses are ascribed using the forecaster radius-mass relation [105]. We calculate mutual Hill stability and redraw if the system is Hill unstable [106]–[108]. Finally, to place our planet sample on even footing with Z23, we keep only planets with orbital period less than 40 days, to compare against Z23.

2.4 Detection and Completeness↩︎

After generating planetary systems around the stellar sample, we “observe" the sample to determine the transit yield. We determine which of the planets geometrically transit, based on the orbital parameters drawn in Section 2.3. To do so, we use Equation 7 in [109]: \[b = \frac{a*\cos{i}}{R_*} \frac{1-e^2}{1+e*\sin{\omega}}, \label{eq:impact95parameter}\tag{3}\] where \(a\) is the semi-major axis, \(i\) is the mutual inclination of the planet’s orbit from its system’s midplane, \(e\) is the planet’s orbital eccentricity, and \(\omega\) is the longitude of periastron.

Next, we use two different prescriptions of sensitivity, to respect the difference in the sensitivity functions of Kepler and K2. For Kepler, the signal to noise ratio (SNR) is calculated following Equation 4 in [110]: \[SNR = \sqrt{\frac{t_{obs}*f_0}{P}} \frac{(R_p/R_*)^2}{CDPP_{eff}}, \label{eq:snr}\tag{4}\] where \[CDPP_{eff} = \sqrt{\frac{t_{CDPP}}{t_{dur}}}{CDPP_N}. \label{eq:cdpp95eff}\tag{5}\]

We then apply the SNR versus detection probability ramp from [111] to establish a detection likelihood (see also Section 3.3 of [27] for an identical application). Each geometrically transiting planet is assigned a detection likelihood, upon which a binary detection or non-detection is drawn.

For K2, we follow the prescription the recovery fraction, \(f(MES)\), from Equations 5 and 6 in [112]: \[MES = C\frac{depth}{CDPP_{t_{dur}}} \sqrt{N_{tr}}, \label{eq:k2-mes}\tag{6}\] and \[f(MES) = \frac{a}{1+e^{-k(MES-1l}}, \label{eq:k2-recovery-fraction}\tag{7}\] where \(CDPP_{t_{dur}}\) = \(CDPP_{6 hr}\) (obtained from Table 2 of [112]), \(N_{tr}\) is calculated using the baseline of the particular K2 campaign and the planet’s period, depth is calculated using the stellar and planet radii, \(C\) is a global correction factor of 0.9488, and a, k, and l are best-fit parameters for the logistic function describing K2’s recovery fraction. These logistic function parameters are provided in Table 3 of [112] for Campaigns 1-8 and 10-18, while for Campaigns 9 and 19, we use the summary values marginalized over AFGK dwarfs (\(a\)=0.6095, \(k\)=0.6088, \(l\)=10.8986). As with the detection likelihood in our Kepler sample, we draw a detection or non-detection for each geometrically transiting K2 planet based on its recovery fraction.

In this study, we use the inverse detection efficiency method (IDEM), in which the detection efficiency is independently computed and then applied to a detected yield in order to estimate a “true" or”adjusted" yield. The detection efficiency is computed along 9 bins each of planet period (log space from 1-40 days) and radius (linear space from 1-4 \(R_{\oplus}\)), and is a product of the correspondingly-binned detection sensitivity and geometric transit probability maps. The geometric transit probability map is constructed once, since it is just a function of the planet semi-major axis and stellar and planet radii, which we assume to vary negligibly over \(Z_{\textrm{max}}\). It is evaluated at each period and radius bin as \((R_* + R_p)/a\), where \(a\) is the planet semi-major axis, and \(R_*\) is assumed to be 1 \(R_{\odot}\). The period is converted to \(a\) using Kepler’s Third Law assuming 1 \(M_{\odot}\).

Meanwhile, a detection sensitivity map is constructed for each of the five \(Z_{\textrm{max}}\)bins (in log space from 100 to 1000 pc). We take 1000 random stars from each height bin of the Kepler sample (and all of the stars from each height bin in the K2 sample, since there are so few of such stars), injecting a planet with an impact parameter of zero at each period and radius bin. We next determine which transiting planets exceed the detection threshold by applying the Kepler and K2 sensitivity functions described above. The detection sensitivity in each period and radius bin is the number of detected planets out of 1000 (or however many K2 planets there are in that particular height bin). As fiducial examples, we show the Kepler and K2 detection sensitivity maps for the lowest and highest \(Z_{\textrm{max}}\)bins in Figure 8. The K2 sensitivity maps are dominated by shot noise due to the much smaller sample, particular in bins where both sensitivity and the planet count is low; this is a known problem in real occurrence calculations and an inherent instability of the inverse detection efficiency method [113]. We check and verify that using these versus the sensitivity maps from [114] and [112] for Kepler and K2, respectively, make no appreciable difference to the adjusted planet yields.

3 Results↩︎

3.1 Results Using Joint Kepler-K2 Sample↩︎

The uncertainty at each \(Z_{\textrm{max}}\)bin reflects the apparent change in planet yield across 30 draws from the completeness function. This uncertainty per \(Z_{\textrm{max}}\)bin depends upon the model occurrence prescription, as it ought: depending on when in time planetary systems are born, they will move between \(Z_{\textrm{max}}\)bins according to kinematic heating. In draws with smaller numbers of planets in that \(Z_{\textrm{max}}\)bin, the resulting occurrence estimate is noisier.

We use the NUTS sampler in the numpyro framework [115] to fit a simplified version of Equation 4 in Z23 to the recovered yield. The model is constructed as follows: \[y = 100 * \eta * Z_{\textrm{max}}^{\tau} * \frac{1}{C * dlnZ_{max}}, \label{eq:power-law}\tag{8}\]

with normalizing parameter \[C = \frac{1}{\tau+1}1000^{\tau+1} - \frac{1}{\tau+1}100^{\tau+1}, \label{eq:power-law-constant}\tag{9}\]

and normalizing bin-size dln\(Z_{\textrm{max}}\)= 0.0011. Our support for \(Z_{\textrm{max}}\)spans 100 to 1000 pc, while \(y\) represents the planet occurrence rate per 100 stars. The free parameters, \(\tau\) and \(\eta\), represent the slope of the planet occurrence-scale height relation and occurrence rate normalization, respectively. \(\tau\) is drawn from a uniform prior of \(U(-1, 1)\) (that is, we allow for the model to fit cases where planet occurrence increases with \(Z_{\textrm{max}}\)), while \(\eta\) is drawn from a uniform prior of \(U(0.01, 1.0)\). We run the MCMC for 1000 warm-up steps, 10000 samples, and 8 chains (making sure acceptance rates were at least 90%), sampling from a normal distribution around \(y\) with error equal to the per-height-bin standard deviation of the completeness-recovered planet occurrence yield.

In Figure 9, we show the resulting occurrence-\(Z_{\textrm{max}}\)trend for the step function models illustrated in blue in Figure 7, overlaid against the trend found by Z23 (in red). Z23 analyzed Super-Earths and Sub-Neptunes separately; in this work, we combine these into one designation of “rocky" planets from 1 to 4 \(R_{\oplus}\), eschewing our ability to comment on the age dependence of the radius valley, in favor of greater statistical power for our toy models. The physical (”true", in dark purple) and recovered (blue) planet occurrences from each model are binned in the same manner as in Z23 – five evenly log-spaced bins from 100 to 1000 pc. We show the 16th and 84th percentile envelope fits to our recovered yields with the blue envelopes in Figure 9, and we report the best-fit planet occurrence trend slope (\(\tau\)) and occurrence normalization (\(\eta\)) for each model in Table 2. For comparison, Z23 found \(\tau\) = -0.30\(\pm\)​0.06 for Super-Earths, and \(\tau\) = -0.37\(\pm\)​0.07 for Sub-Neptunes. After combining these separate observed samples ⁴ and adding their uncertainties in quadrature, we re-fit and find their \(\tau\) to be -0.28\(\pm\)​0.08 and \(\eta\) to be 0.46\(\pm\)​0.02 for small, close-in planets. The 16th and 84th percentile envelopes for this fit are also plotted in red in Figure 9.

Our control (“no change whatsoever") model, in which planet occurrence is stalled at 33% for the age of the Milky Way, results in a constant planet occurrence rate at all galactic heights and is unsurprisingly inconsistent with the Z23 results. However, we find moreover that for the majority of our step function models, the occurrence-\(Z_{\textrm{max}}\)trend is similarly insufficiently steep. In general, \(f_1\)controls the high-\(Z_{\textrm{max}}\)end of the trend and \(f_2\)controls the low-\(Z_{\textrm{max}}\)end. Since the present-day total planet host fraction is fixed, it therefore follows that the limit at \(f_1\)=1% and \(f_2\)=100% sets the maximum slope of the occurrence-\(Z_{\textrm{max}}\)trend for any given \(t\). Using this conceit, it is possible to systematically rule out \(t\)because these limiting cases fail to produce a sufficiently steep trend. The models that are not ruled out have step thresholds \(t\)=2.2-3.2 Gyr ago. The models with \(t\)=2.7 and \(t\)=3.2 Gyr are contrived to push the possible window of step function model thresholds as wide as possible, with the former reaching the theoretical maximum \(f_2\)and the latter reaching the theoretical minimum \(f_1\). In so doing, we can rule out models in which the Galactic small planet host fraction around FGK dwarfs jumped (in a relatively short time) any earlier than 2.2 Gyr ago or any later than 3.2 Gyr ago.

In general, models that fail to match Z23 do not produce the necessary overabundance of planets with maximum oscillation amplitudes close to the Galactic midplane and relative underabundance of planets that travel to higher vertical excursions from the Galactic midplane, instead producing a flatter slope across the relevant \(Z_{\textrm{max}}\)range. This issue persists among the piecewise models, which are parameterized as a flat planet host fraction \(f_1\)until some time \(t\)at which the fraction of planet hosts gradually increases to a present-day level of \(f_2\). For this class of models, we rule out those with \(t\)=6.2 Gyr ago or earlier, as well as \(t\)=3.2 Gyr ago or more recently. For most of the favored models from both archetypes, the small planet occurrence-\(Z_{\textrm{max}}\)slope, \(\tau\), is barely within one standard deviation of the Z23 slope.

3.2 Results Using TRILEGAL↩︎

Without knowing the noise properties of the TRILEGAL stellar samples, we are limited to comparing the true, un-adjusted planet occurrence to Z23. As such, the errorbars reported in Figure 11 are much smaller than those in Figures 9 and 10, driven entirely by the age and \(Z_{\textrm{max}}\)uncertainties rather than also by the sensitivity functions. Like with HU25-B20-B25, we deploy a battery of step and piecewise functions; however, since the samples differ somewhat in their stellar age distributions, we do not use exactly the same models but instead modify the \(f_1\)and \(f_2\)for each threshold. Four representative models for each archetype are shown in Figure 11. The piecewise models are anchored by a maximum \(f_2\)at \(t\)=3.2 Gyr that is unable to achieve a steep enough \(\tau\) and a minimum \(f_1\)at \(t\)=8.2 Gyr that is similarly too shallow. The step models are anchored by a minimum \(f_1\)at a relatively recent \(t\)=3.2 Gyr that is too shallow, but unlike with HU25-B20-B25, we are able here to identify models prescribing recent step increases in \(f\)that fit more closely to the Z23 result. Step function {\(f_1\)=20%, \(f_2\)=100%, \(t\)=1.7 Gyr} yielded a \(\tau\) of -0.20\(\pm\)​0.02 and \(\eta\) of 0.46\(\pm\)​0.01, and step function {\(f_1\)=5%, \(f_2\)=100%, \(t\)=2.2 Gyr} yielded a \(\tau\) of -0.24\(\pm\)​0.02 and \(\eta\) of 0.43\(\pm\)​0.00. Despite the better match, these two models are tuned to the theoretical maximum \(f_2\)=100%.

3.3 Results Using Only Kepler↩︎

Thus far, we have sought to match the Z23 joint Kepler-K2 occurrence rate, but our study is limited by a lack of K2 stars with age measurements. We now consider what the Kepler planet occurrence-\(Z_{\textrm{max}}\)relation would look like, and whether an instantaneous or gradual rise in \(f\)would be able to reproduce this trend. We prepared the stellar and planetary sample as described in Z23, but only for the Kepler field, and propagate the observed planet occurrence-\(Z_{\textrm{max}}\)trend through the same completeness pipeline in Z23. Without \(Z_{\textrm{max}}\)-disaggregated sensitivity maps, we simply bin all planets by in 9x9 radius and period space regardless of \(Z_{\textrm{max}}\)and multiply by the reliability map and divide by the sensitivity map from [114]. We divide this by the same geometric transit probability map from the joint HU25-B20-B25 sample, and then sum the adjusted planets per \(Z_{\textrm{max}}\)bin to get the Kepler planet occurrence per 100 stars as a function of \(Z_{\textrm{max}}\): {49.1\(\pm\)​3.1, 43.4\(\pm\)​1.7, 36.5\(\pm\)​1.5, 37.4\(\pm\)​2.4, 30.2\(\pm\)​3.8}. We fit this using the same power law from Z23 and get a best-fit \(\tau\) of -0.25\(\pm\)​0.05 and \(\eta\) of 0.35\(\pm\)​0.01. While the Kepler occurrence is understandably less than the Z23 result (by \(\sim\)​0.10, which suggests that that is the K2 small planet occurrence around FGK dwarfs), it is notable that the planet occurrence-\(Z_{\textrm{max}}\)trend slope is essentially the same. This suggests that the observed trend from Z23 should hold across the Kepler and K2 fields.

Using this as our new “ground truth", we attempt to model step and piecewise increases in \(f\)over time as before. We find a relatively good match in the step function \(f_1\)=20%, \(f_2\)=100%, \(t\)=2.2 Gyr} (depicted in the top panel of Figure 12), which produced a \(\tau\) of -0.20\(\pm\)​0.07 and a \(\eta\) of 0.34\(\pm\)​0.01. This model has the same \(t\)and \(f_2\)as one of the successful models for HU25-B20-B25, but it has a higher \(f_1\), likely due to differences in the age distribution between the joint and Kepler-only samples. The latter contains more older stars, requiring a boost to either \(f_1\)or \(f_2\)to maintain the correct occurrence normalization. The age histograms are shown in the bottom panel of Figure 12. This difference in age distributions may also be responsible for the slight bump in planet occurrence at the 300 pc \(Z_{\textrm{max}}\)bin. This deviation from a monotonically decreasing trend is persistent across model thresholds and archetypes.

4 Discussion↩︎

We have demonstrated in the preceding sections that an increase in the planet host fraction over the last 8 Gyr – however gradual – is insufficient to reproduce the steepness of the observed decrease in small planet occurrence with increasing maximum oscillation amplitude from the Galactic midplane. We synthesized stellar and planet populations following a battery of models in which the planet host fraction, \(f\), increased either instantaneously or over Gyr timescales, with the increase (or onset of the increase) occurring at various times in the Milky Way’s past. Of these models, we ruled out all but the following: a step increase 2.2 Gyr ago from \(f_1\)= 10% to \(f_2\)= 100%, a step increase 2.7 Gyr ago from \(f_1\)= 1% to \(f_2\)= 95%, a step increase 3.2 Gyr ago from \(f_1\)= 1% to \(f_2\)= 75%, a gradual increase starting 4.2 Gyr ago from \(f_1\)= 5% to \(f_2\)= 100%, and a gradual increase starting 5.2 Gyr ago from \(f_1\)= 1% to \(f_2\)= 95%. These models are not ruled out because their best-fit \(\tau\) and \(\eta\) are within a standard deviation from the Z23 combined and re-fit Super-Earth and Sub-Neptune \(\tau\) and \(\eta\). However, these models’ best-fit \(\tau\) lie on the edge of the acceptable \(\tau\) envelope of \(\tau\) = 0.28\(\pm\)​0.08. Moreover, for all of these models, \(f_1\)or \(f_2\)was pushed to an unrealistic minimum or maximum, respectively. In two cases – the step increase 2.7 Gyr ago and the gradual increase beginning 5.2 Gyr ago – we nearly pushed both \(f_1\)and \(f_2\)to their extrema.

It is informative that almost none of our models applied to the HU25-B20-B25 sample yield a relation that is steeper than the Z23 slope. The Z23 stellar sample did not have stellar ages, but it is possible at least to compare their \(Z_{\textrm{max}}\)distribution, which peaks between 288 and 408 pc, to that of our joint sample, which peaks at 263 pc. While we constructed our sample in much the same way as Z23, our requirement that stars in our sample must have age measurements means that our samples are not identical. Most significantly, our sample’s \(Z_{\textrm{max}}\)distribution is bunched around lower \(Z_{\textrm{max}}\), leading to proportionally larger errorbars for the last two Galactic height bins and reduced sensitivity to changes in planet occurrence over this \(Z_{\textrm{max}}\)region.

The TRILEGAL sample likewise has a lower instantaneous height peak and smaller MAD than the Kepler and K2 \(Z_{\textrm{max}}\)in Z23. However, it has much more precise stellar ages than the HU25-B20-B25 sample. We demonstrated in Figure 11 that a late-time instantaneous rise in \(f\)can match the Z23 result. While this shows that it is possible in principle for an increase in \(f\)to produce a strongly downward sloping planet occurrence with increasing \(Z_{\textrm{max}}\), the requirement among the best-fit models that either \(f_1\)= 1% or \(f_2\)= 100% suggests that \(f\)cannot single-handedly drive the overall small planet occurrence-\(Z_{\textrm{max}}\)relation.

Future work could be more successful at matching Z23 by adding complexity to our model in two ways. First, Z23 conducted their analysis on Super-Earths and Sub-Neptunes separately, while our models were agnostic about this distinction. In Z23, it was the Sub-Neptunes that drove the steepness of the planet occurrence-\(Z_{\textrm{max}}\)relation, whereas the Super-Earth occurrence was flat until 500 pc. Therefore, allowing for different \(f\)for Super-Earths and Sub-Neptunes can build in the flexibility required to match the Z23 trend. Second, while psps allows users to modulate the fraction of compact multis, in this study we have left the so-called “intact" fraction relatively constant at 18\(\pm\)​10% [27]. Allowing this to be a free parameter in addition to \(f\)could provide the means to produce greater differences in planet occurrence for stars of different ages. Since the issue with most of our models is that they are too shallow, we posit that swapping an increase in \(f\)with an increase in the intact fraction could result in models that better fit the Z23 result.

Another potential systematic cause of the diluted occurrence-\(Z_{\textrm{max}}\)trend in our models is the weak age-\(Z_{\textrm{max}}\)relation in our joint Kepler-K2 sample. This relation was measured to be 4 Gyr/kpc using asteroseismic ages [62]. ESA’s upcoming PLAnetary Transits and Oscillations of stars (PLATO) mission will measure precise solar-like oscillations for 245,000 FGK dwarf stars, enabling much more precise ages for planet hosts on a demographically relevant scale [116]–[119]. With a large, uniform sample of relatively high-precision ages, it will be possible to re-examine age-based models of Galactic-scale exoplanet demographics.

While the models we have failed to rule out are not very plausible, we would be remiss not to briefly explore their implications on planet formation in the Milky Way’s history, if only to present a framework for future work. Among all the step and piecewise models for time-varying planet occurrence that we applied to the HU25-B20-B25 and TRILEGAL joint Kepler-K2 stellar samples, all models that matched Z23 were step functions and had \(t\)no earlier than 3.2 Gyr ago. If true, this is suggestive of a single event rather than a series of events (or a more gradual process such as Galactic Chemical Evolution) that is overwhelmingly responsible for the strong correlation between planet occurrence and \(Z_{\textrm{max}}\). Such a process would have had to be relevant for not only the Kepler field but also the K2 fields. In the following paragraphs, we discuss whether certain Galactic scale events in the Milky Way’s past could have coincided with one of these late-time boosts to \(f\). These events include major and minor mergers, close passages of satellite galaxies, and global gravitational instabilities or external infall, which could drive fluctuations in the interstellar medium (ISM).

Galaxy major mergers are thought to be one of the main factors triggering star formation in galaxies and are predicted to have clear effects on their chemical evolution [120]–[123]. The timing of these dynamical events is a potential clue, but we must also consider the nature and scale of the actual processes that might affect planet occurrence. We expect mergers to affect the planet occurrence rate chemo-kinematically by adding gas of lower and heterogeneous metallicity to the planet forming budget and stirring up the interstellar medium (ISM). Mergers could also compress the ISM during passage and thus stimulate not only additional star formation but also persistent spiral structure [124]. On the other hand, close passages might contribute to ISM turbulence without modifying its chemistry. These processes could also each have competing effects for planet formation. For example, [125] showed that the injection of metal-poor gas from satellite infall could have triggered the formation of metal-poor stars as recently as \(\sim\)​3 Gyr ago, which would affect planet occurrence among younger stars.

While not all of these mechanisms directly affect planet formation, they can impact (1) the formation conditions for protoplanetary disks and (2) the kinematic signatures or even survival of mature planetary systems. To that end, we can ask, based on our timing prescriptions in this paper: does the timing line up for when we need a planet-moderating effect to occur? Sources in the literature point to a number of potential significant dynamical interactions that our Galaxy has had with external bodies over the course of its history: the Gaia-Sausage-Enceladus (GSE) dwarf galaxy merger 8-10 Gyr ago [126], [127]; the Sagittarius dwarf spheroidal (Sgr dSph) galaxy first, second, and third passages 5.7, 1.9, and 1.0 Gyr ago, respectively [123]; and the Virgo Radial Merger (VRM) 2.7 Gyr ago [128]–[130], sometimes called the Virgo Overdensity, [131]. Additionally, satellite infall of dwarf galaxies spans the Milky Way’s history, but we are primarily interested in ones large and close enough to significantly shift the planet formation or evolution conditions in the Solar neighborhood. The favorable window includes the VRM merger event [128]–[130], [132], a putative radial dwarf galaxy merger approximately 2.7 Gyr ago that produced debris in the Solar neighborhood. It is also probable that the Milky Way disk has experienced multiple gas-rich radial mergers [132]–[135] similar to the VRM that could have stimulated episodes of enhanced star formation. The Sagittarius dwarf galaxy’s second passage occurred [123] 1.9 Gyr ago and coincided with a narrow period of enhanced star formation. However, given that Sagittarius’ closest passage was at \({\gtrsim}15\) kpc from the Galactic center, it is unclear how much this would have impacted planet formation in the Solar Neighborhood.

Our step function models’ \(t\)parameters mark the timing of some of these events: \(t\)=1.7 Gyr corresponds to 200 Myr after the second passage of the Sagittarius dwarf galaxy; \(t\)=2.2 Gyr corresponds to 500 Myr after the VRM; \(t\)=6.2 Gyr corresponds to shortly before the first passage of Sgr dSph; and \(t\)=8.2 Gyr lies toward the end of the window of the Gaia-Enceladus merger. \(t\)=4.2 Gyr marks a fiducial waypoint in the Milky Way’s dynamical history. It is difficult to interpret more recent \(t\)because of our sample’s large isochrone age uncertainties, while probing \(t\)at cosmological times earlier than \(\sim\)​5.5 Gyr would require extrapolating beyond the oldest stars in our sample. Using strictly timing arguments, therefore, our analysis can rule out the imprint of the GSE merger and the first passage of Sgr dSph on planet occurrence in our Galactic neighborhood.

Mergers or close passages offer a potentially illuminating framework here in another way, which is to mediate the apparent link between planet occurrence and stellar kinematics [36]–[38]. Infalls themselves could have dynamically heated stars in the path of their trajectories, puffing up the orbits of older stars to greater scale heights. Similarly, stellar radial migration [136]–[140] could have moved stars from the metal-rich center outward, and from the metal-poor periphery inward, where there is evidence that the Sun itself is an interloper from the inner disk [141], [142]. On the other hand, older populations may very well have formed in a thick, turbulent disk [57], [143]–[145]. There is an important relationship between \(Z_{\textrm{max}}\), the subject of Z23, and vertical action, \(J_{z}\): both are dynamical quantities that capture the amplitude of a star’s oscillation above and below the Galactic midplane and generally increase with age [66], [146]–[148]. Action-angle coordinates, including \(J_{z}\), are powerful tools for tracing a star’s long-term kinematic history due to dynamical disturbances [149]. Famously, the spiral substructure of the “Antoja snail" [150] is encoded in action-angle coordinates and possibly attributable to the passage of Sagittarius referenced above [151] – in this sense, an occurrence dependent upon \(Z_{\textrm{max}}\)may be interpretable in the context of galactic disturbances.

For now, it is not clear whether information from galactic processes (e.g. mergers, satellite in-fall) can reliably propagate down to local processes such as turbulence, or whether the local ISM turbulence at planet formation sites is instead causally disconnected from Galactic chemical evolution (GCE) and kinematics. [152] suggested that ISM turbulence could prolong protoplanetary disk lifetimes among 20-70% of systems, which could affect the final accounting of mature planets. A low level of turbulence can enable local vortices at planet formation sites to efficiently trap and concentrate dust, lowering the threshold for planetesimal formation from near solar-like metallicity down to 0.08 Z_Solar for Mars-like planets [153], although higher turbulence levels may require higher metallicity thresholds for planet formation. Observational evidence was shown by [154], who found that protoplanetary disks around metal-poor Sun-like stars are longer-lived than previously thought. If the average ISM turbulence were to be more favorable among younger metal-poor stars than older ones, this could contribute to the observed trend in Z23.

We have so far focused upon some process enhancing recent star and planet formation relative to some former baseline. A seeming recent rise in recent planet occurrence could also be produced by a process progressively carving away at older planetary systems. Such a process would act by diminishing planet occurrence preferentially around older stars, perhaps in a way where the cumulative probability of destruction or ejection accrues with time. The above scenarios that we have touched upon so far invoke galactic dynamical history as a way to drive bursts in star (and potentially planet) formation, rather than directly acting on the planets after formation. Stellar flybys have been suggested as a possible direct mechanism for planet loss and the altering of system architectures [28], [29]. In concert with additional secular effects like Kozai-Lidov [155], close stellar flybys in dense cluster environments can sometimes lead to Hot Jupiters and ultra-cold Saturns, although they more often lead to planet ejections [156]. [157] show through N-body simulations that stellar flybys as far away as 1000 AU can disrupt low-order resonant chains, although smaller planets are relatively more resilient to this effect. Along with [158], they also demonstrate cases in which a delayed disruption of the architecture occurs, over 100 Myr timescales. It is also possible that planet formation was suppressed earlier on from the radiation environment during formation. [159] showed that the protoplanetary disks of thick disk stars, which formed closer to cosmic noon, experienced \(\sim\)​7 orders of magnitude more background radiation than Solar neighborhood stars, limiting their lifetimes to 0.2-0.5 Myr.

The menu of culprits listed above that could affect planet formation and survivability is long but can be divided into varying operative timescales. If we adopt the premise that a late-time apparent increase in planet formation has occurred, the balance of the processes listed above must continually result in an increasingly net gain of planets over time. Sudden jumps in planet occurrence at earlier cosmological times cannot explain the Z23 result. This means that present-day exoplanet demographics were either not set by earlier Galactic-scale events (e.g. the Gaia-Enceladus merger), or whatever imprints were left by those events have now been erased. Determining how this occurs – perhaps through some combination of deleterious effects that become ameliorated over time (e.g. elevated radiation environment in high-traffic dense birth clusters) and increasingly fertile grounds for planet formation (e.g. the average local ISM turbulence becoming more favorable among more recently-formed metal-poor stars) – would require synthesizing observational and simulation studies of star and planet formation sites encompassing the Solar neighborhood and spanning a \(\pm\)​1 kpc column on either side of the Galactic midplane.

5 Conclusion↩︎

The impact of the galactic environment on exoplanet demographics has been the subject of much recent study. [67] showed that the metallicity gradient expected from Galactic chemical evolution is insufficient to explain a decreasing short-period planet occurrence with increasing galactic height, as observed by the Kepler and K2 missions. In this work, we forward modeled two prescriptions for increasing recent planet occurrence, as well as a flat (control) model, producing observed synthetic populations of planetary systems using the software package psps. By assigning planet populations, tuned to reflect different “burst" times in planet occurrence, to a synthetic Kepler stellar sample, we show that the timing and nature of that burst produce a range of slopes in \(Z_{\textrm{max}}\)versus planet occurrence space.

Generally speaking, we find that an increase in the small planet host fraction among FGK dwarfs is insufficient to reproduce the steepness of the downward trend in planet occurrence with increasing \(Z_{\textrm{max}}\). There are some models – exclusively step increases (rather than more gradual increases in \(f\)) – at recent Galactic times that can match the Z23 result to within one standard deviation, but these models are tuned to implausibly high \(f_2\)or low \(f_1\), and they are not strong matches. Part of this may be attributable to our lack of precise stellar ages (or, indeed, ages of any quality for K2 stars), but even an idealized 1 Gyr age uncertainty from the TRILEGAL sample does not appreciably change the intrinsic planet occurrence yields of the models explored in this work. As future work, perhaps tuning the intact fraction rather than the planet host fraction will provide a better match to Z23. In this study, we also independently derived the Kepler small planet occurrence trend with \(Z_{\textrm{max}}\)and find that it is relatively easier for a step increase in \(f\)to match this slope.

We have focused in this manuscript about the timing of the effect required, rather than a particular mechanism. We center our discussion upon the ways that galactic processes could potentially drive planet formation on the right timescales, but as yet, we cannot distinguish between a process that enhances recent planet formation, versus one that mimics the same trend by progressively carving away at older planetary systems. A satisfactory theory of the role of Galactic-scale events in planet formation and survivability requires making the connection between the spatial scales of these two seemingly disparate processes. Specifically, it is important to pose one or more plausible physical mechanisms for directly affecting either the planet formation stage or the ability of mature planetary systems to keep hold of their planets. There are major conceptual challenges in connecting the spatially and temporally extended processes relevant to the Galaxy on Gyr timescales to the smaller-scale physics of planet assembly and evolution. However, the canonical picture of local metal enrichment as the only galactic-planetary connective tissue cannot fully explain planet demographics in the Gaia era. The puzzles of why thin/thick disk planet populations differ [37], [118], [159]–[161], whether and how common kinematics drive patterns in planet demographics [34]–[36], [38], [162]–[164], and whether and how planetary systems evolve over Gyr in isolation or in response to external perturbations [23], [26], [29], [31], [165], [166] merit a larger-scale contextual discussion. Continued observational and simulation studies connecting GCE and Galactic dynamics to giant molecular cloud (GMC) and ISM dynamics, and then connecting GMCs and the ISM to protoplanetary disks, will enable stronger statements about the role of the Galactic context in exoplanet demographics.

6 Acknowledgments↩︎

We wish to thank Luke Bouma, Jon Zink, Pat Tamburo, Desika Narayanan, Jamie Tayar, Adrian Price-Whelan, Carrie Filion, George Privon, Quadry Chance, Natalia Guerrero, Jason Dittmann, William Schap III, and the UF Astronomy Department Stars & Planets Journal Club for their helpful comments and suggestions. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the NASA Exoplanet Archive [167], which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This material is based upon work supported in part by the National Science Foundation GRFP under Grant No. 1842473.

KJD acknowledges support from the Heising Simons Foundation grant # 2022-3927. She also respectfully acknowledge that the University of Arizona is home to the O’odham and the Yaqui. She respects and honors the ancestral caretakers of the land, from time immemorial until now, and into the future.

We acknowledge that IPAC/Caltech resides on the traditional, ancestral and unceded territory of the Gabrielino/Tongva peoples, the original caretakers of this land. We respectfully recognize the Gabrielino/Tongva peoples who still reside in Tovaangar (the Los Angeles basin and South Channel Islands). We also acknowledge that the main campus of the University of Florida is located on the ancestral territory of the Potano and of the Seminole peoples. The Potano, of Timucua affiliation, lived here in the Alachua region from before European arrival until the destruction of their towns in the early 1700s. The Seminole, also known as the Alachua Seminole, established towns here shortly after but were forced from the land as a result of a series of wars with the United States known as the Seminole Wars. We, the authors, acknowledge our obligation to honor the past, present, and future Native residents and cultures of California and Florida.

References↩︎