Planet Mass Function around M Stars at 1–10 au: A Plethora of Sub-Earth Mass Objects

The abundance of small planets (≲1 M_⊕) remains a terra incognita in the field of exoplanetary science. Given their low masses and radii, they are incredibly difficult to find using techniques that rely on a planet's influence on its host star. Consequently, planets akin to the terrestrial planets in the solar system have only been found around other stars on extremely short orbits. Their distribution is largely unknown even at Mercury-like distances (0.3 au, ∼100 days) and it is completely unknown for distances ≳1 au (e.g., D. C. Hsu et al. 2019, 2020; W. Zhu & S. Dong 2021; A. Dattilo et al. 2023). Answering the question of whether our solar system is common or not therefore remains largely out of reach.

Fortunately, a dramatic change in our knowledge of small planets at large orbital distances is imminent. The upcoming Nancy Grace Roman Space Telescope will open new frontiers in the characterization of exoplanet populations at distances of 1–10 au. Using the microlensing technique, which detects planets that gravitationally lens background stars and cause a brief jump in their brightness (B. Paczynski 1986; A. Gould & A. Loeb 1992), it will enable us to detect objects as small as Ganymede (M. T. Penny et al. 2019; W. Zhu & S. Dong 2021)! Hitherto, the lowest-mass planets that microlensing surveys from the ground have managed to detect are super-Earths (e.g., W. Zang et al. 2023). The impact of pushing the detection limit to planets that are 2 orders of magnitude less massive will likely be similar to that of the Kepler Space Telescope more than a decade ago.

The number of small planets that we expect to detect with the Nancy Grace Roman Space Telescope is heavily dependent on how abundant such planets are, i.e., the planet mass function. Previous planet yield estimates have used constant or power-law planet occurrence rate models for low-mass planets (M. T. Penny et al. 2019) that either extrapolate from ground-based microlensing detections (A. Cassan et al. 2012) or are inferred from Kepler planets assuming a fixed mass–radius relation (J. J. Lissauer et al. 2011). The aim of our study is to provide a physically motivated planet mass function at 1–10 au for planets less massive than super-Earths. In particular, we quantify the planet mass function expected at birth as a result of pebble accretion (Section 2). Using pebble accretion efficiency and the measured disk mass function, we calculate the fraction of stars that host low-mass planets in Section 3.

Dynamical evolution over long timescales can alter the natal planet mass function of both bound and unbound ("free-floating") planets. We perform a suite of N-body simulations for low-mass planets to characterize the effects of dynamical evolution on the planet mass function for a range of initial conditions. The details and results of the simulations are presented in Section 4. Finally, in Section 5, we summarize our results and discuss the implications of our work for observations by the upcoming Nancy Grace Roman Space Telescope.

We hypothesize that planetary seeds emerge in protoplanetary disks as a result of some solid clumping (e.g., streaming instability; A. Johansen & M. Lambrechts 2017; R. Li & A. N. Youdin 2021) by the beginning of the Class I stage. These seeds then grow by accreting the remaining millimeter-to-centimeter-sized dust (pebbles) that is drifting past them. The final mass to which these seeds grow depends on the amount of dust present in the disk and the efficiency with which this dust is accreted.

The efficiency with which the drifting pebbles are accreted is given by the ratio of the pebble accretion rate and the rate at which pebbles drift

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{\dot{M}}_{\mathrm{peb}}}{{\dot{M}}_{\mathrm{drift}}}.\end{eqnarray} \tag{ 1 }$$

Detailed calculation of is outlined in Y. Chachan & E. J. Lee (2023), which we adopt. Here, we only highlight the key results. The drift rate of pebbles ${\dot{M}}_{\mathrm{drift}}=2\pi {v}_{{\rm{r}}}{{\rm{\Sigma }}}_{{\rm{d}},\ \mathrm{St}}$, where v_r is the radial velocity of the pebbles and Σ_{d, St} is the surface density of pebbles with Stokes number St. The radial velocity of pebbles is the sum of coupled radial motion due to the gas and the radial drift of pebbles due to slightly sub-Keplerian velocity of the pressure-supported gas:

$$\begin{eqnarray}&&{v}_{{\rm{r}}}=-\displaystyle \frac{3}{2}\displaystyle \frac{\nu }{r}\displaystyle \frac{1}{1+{\mathrm{St}}^{2}}-2\eta {v}_{{\rm{K}}}\displaystyle \frac{\mathrm{St}}{1+{\mathrm{St}}^{2}},\end{eqnarray} \tag{ 2 }$$

where r is the radial distance from the star, ν = α_t c_s H_g is the kinematic viscosity of the gas, α_t is the Shakura–Sunyaev parameter, c_s and H_g = c_s/Ω_K are, respectively, the local sound speed and scale height of the gas, Ω_K is the Keplerian frequency, v_K = Ω_K r is the Keplerian velocity, and the negative sign indicates inward motion toward the central star. The quantity $\eta =-0.5\gamma {({c}_{{\rm{s}}}/{v}_{{\rm{K}}})}^{2}$ quantifies the deviation of gas's motion from Keplerian and $\gamma =d\mathrm{ln}{P}_{{\rm{g}}}/d\mathrm{ln}r$ is the logarithmic pressure gradient.

The pebble accretion rate ${\dot{M}}_{\mathrm{peb}}$ is given by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{peb}}=2{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{St}}{R}_{\mathrm{acc}}{v}_{\mathrm{acc}}\times \min (1,{R}_{\mathrm{acc}}/{H}_{{\rm{d}}}),\end{eqnarray} \tag{ 3 }$$

where pebbles traveling at a relative speed of v_acc are accreted onto the planet if they approach closer than R_acc, and ${H}_{{\rm{d}}}\equiv {H}_{{\rm{g}}}{[{\alpha }_{{\rm{t}}}/({\alpha }_{{\rm{t}}}+\mathrm{St})]}^{1/2}$ is the scale height of the pebble disk. We then have

$$\begin{eqnarray}\begin{array}{l}{\dot{M}}_{\mathrm{peb}}\sim \left\{\begin{array}{cc}{\left(192\right)}^{1/3}{\mathrm{St}}^{2/3}{q}^{2/3}{{\rm{\Omega }}}_{{\rm{K}}}{{\rm{\Sigma }}}_{{\rm{d}},\ \mathrm{St}}{r}^{2} & (2{\rm{D}},\ \mathrm{hw})\\ 4{\eta }^{1/2}{\mathrm{St}}^{1/2}{q}^{1/2}{{\rm{\Omega }}}_{{\rm{K}}}{{\rm{\Sigma }}}_{{\rm{d}},\ \mathrm{St}}{r}^{2} & (2{\rm{D}},\ \mathrm{sh})\\ \displaystyle \frac{8{{\rm{\Sigma }}}_{{\rm{d}},\ \mathrm{St}}{{GM}}_{{\rm{p}}}\mathrm{St}}{{c}_{s}}{\left(\displaystyle \frac{{\alpha }_{{\rm{t}}}+\mathrm{St}}{{\alpha }_{{\rm{t}}}}\right)}^{1/2} & (3{\rm{D}})\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

where G is the gravitational constant, M_p is the protoplanet mass, and q = M_p/M_⋆ is the protoplanet mass ratio. Here, (2D, hw) case refers to R_acc > H_d , v_acc = η v_k (headwind); (2D, sh) to R_acc > H_d , v_acc = 3Ω_k R_acc/2 (shear); and (3D) to R_acc < H_d.

Putting the expressions for ${\dot{M}}_{\mathrm{peb}}$ and ${\dot{M}}_{\mathrm{drift}}$ together and for pebbles with St ≲ 1 such that (1 + St²) ∼ 1, we obtain the following expressions for accretion efficiency in different regimes:

$$\begin{eqnarray}&&\begin{array}{l}\epsilon =\\ \left\{\begin{array}{cc}\displaystyle \frac{2}{3\pi }{\left(\displaystyle \frac{q\mathrm{St}}{\eta }\right)}^{1/2}| \gamma | {\left({\alpha }_{t}+\displaystyle \frac{2}{3}| \gamma | \mathrm{St}\right)}^{-1} & (2{\rm{D}},\ \mathrm{hw})\\ \displaystyle \frac{{\left(192\right)}^{1/3}}{6\pi }{\left(q\mathrm{St}\right)}^{2/3}\left(\displaystyle \frac{| \gamma | }{\eta }\right){\left({\alpha }_{t}+\displaystyle \frac{2}{3}| \gamma | \mathrm{St}\right)}^{-1} & (2{\rm{D}},\ \mathrm{sh})\\ \displaystyle \frac{4}{3\pi }q\displaystyle \frac{r}{{H}_{{\rm{d}}}}\displaystyle \frac{| \gamma | \mathrm{St}}{\eta }{\left({\alpha }_{t}+\displaystyle \frac{2}{3}| \gamma | \mathrm{St}\right)}^{-1}. & (3{\rm{D}})\end{array}\right.\end{array}\end{eqnarray} \tag{ 5 }$$

This expression is only valid for planets on circular and coplanar orbits—we neglect the effects of nonzero eccentricity or inclination on (e.g., see C. W. Ormel & B. Liu 2018). Given that there is a range of pebble sizes present in protoplanetary disks, we calculate the mass-averaged accretion efficiency $\bar{\epsilon }$ assuming dust grains are in the Epstein drag regime (St ∝ a ⇒ m ∝ St³, where a and m are the grain size and mass, respectively):

$$\begin{eqnarray}&&\bar{\epsilon }({M}_{{\rm{p}}})=\displaystyle \frac{\int \epsilon ({M}_{{\rm{p}}},\mathrm{St})\,n(\mathrm{St})\,m(\mathrm{St})\,d\mathrm{St}}{\int n(\mathrm{St})\,m(\mathrm{St})\,d\mathrm{St}},\end{eqnarray} \tag{ 6 }$$

with the grain size distribution given by $n(\mathrm{St})={dn}/d\mathrm{St}$. For most of our parameter space, fragmentation limits the maximum grain size to Stokes number St${}_{\mathrm{frag}}={v}_{{\rm{f}}}^{2}/3{\alpha }_{{\rm{t}}}{c}_{{\rm{s}}}^{2}$, where v_f is the fragmentation velocity of dust grains (F. Brauer et al. 2008). The size distribution n(St)m(St) ∝ St^−0.5 for St < 2α_t/π (turbulent regime) and n(St)m(St) ∝ St^−0.75 for St > 2α_t/π (settling regime) (T. Birnstiel et al. 2011). Radial drift is the limiting process only at large distances (≳5 au) when St_frag is large (α_t = 10⁻⁴, v_f = 10 m s⁻¹). For radial drift dominated regime, the maximum Stokes number St${}_{\mathrm{drift}}=\delta {v}_{{\rm{K}}}^{2}/| \gamma | {c}_{{\rm{s}}}^{2}$, where δ = 0.01 is the local dust-to-gas mass ratio. The size distribution is more top heavy than in the fragmentation dominated regime with n(St)m(St) ∝ St^0.5 (this corresponds to n(a) ∝ a^−2.5; T. Birnstiel 2024). We integrate over St in the range [10⁻⁶, (min (St_frag, St_drift)] to obtain $\bar{\epsilon }$, where the lower limit is chosen so as to capture grains that are well coupled to the gas and that move radially due to gas's radial motion.

Equipped with mass-averaged accretion efficiency $\bar{\epsilon }$, we can calculate the amount of dust mass needed to grow a planet from an initial mass M₀ to a given mass M_p:

$$\begin{eqnarray}&&{M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}={\int }_{{M}_{0}}^{{M}_{{\rm{p}}}}\displaystyle \frac{1}{\bar{\epsilon }}\,{{dM}}_{{\rm{p}}}.\end{eqnarray} \tag{ 7 }$$

This differs from our calculations in Y. Chachan & E. J. Lee (2023) where we calculated the dust mass needed to grow a seed to pebble isolation mass M_iso (B. Bitsch et al. 2018):

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{iso}} & = & 25{M}_{\oplus }\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right){\left(\displaystyle \frac{{H}_{{\rm{g}}}/r}{0.05}\right)}^{3}\left[1-\displaystyle \frac{\gamma +2.5}{6}\right]\\ & & \times \left[0.34{\left(\displaystyle \frac{-3}{{\mathrm{log}}_{10}{\alpha }_{{\rm{t}}}}\right)}^{4}+0.66\right].\end{array}\end{eqnarray} \tag{ 8 }$$

Here, we focus on planet mass M_p ≤ M_iso. Of particular note and importance in this work is the efficiency's dependence on planet mass. For low-mass planets that are the focus of this work, pebble accretion primarily happens in the 3D regime. In the 3D accretion regime, ∝ M_p and, as a result, ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}\propto \mathrm{log}({M}_{{\rm{p}}}/{M}_{0})$ (see left panel of Figure 1), leading to the difference in ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ at regular logarithmic interval of M_p being the same (see how the span of required disk mass between mass bin edges are linearly uniform in the right panel of Figure 1). ⁶ In the 2D headwind regime, $\epsilon \propto {M}_{{\rm{p}}}^{2/3}$ and thus ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}\propto ({M}_{{\rm{p}}}^{1/3}-{M}_{0}^{1/3})$, which leads to a similar qualitative result. This dependence of ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ on M_p and M₀ combined with the disk mass distribution sets the planet mass function at birth. For our fiducial model, we adopt a Shakura–Sunyaev viscosity parameter α_t = 10⁻⁴ and a fragmentation velocity of dust grains v_f = 1 m s⁻¹, but we vary these parameters to quantify their effect on the planet mass function in Section 3.

Zoom In Zoom Out Reset image size

Equation (5) indicates that the temperature structure of the disk plays a central role in setting the accretion efficiency (through H_g and η) and the maximum pebble size (St${}_{\mathrm{frag}/\mathrm{drift}}\propto 1/{c}_{{\rm{s}}}^{2}$). We adopt the analytical prescription for temperature structure from S. Ida et al. (2016) based on results of A. Oka et al. (2011) and P. Garaud & D. N. C. Lin (2007):

$$\begin{eqnarray}&&{T}_{\mathrm{disk}}=\min (2000\,{\rm{K}},\max ({T}_{\mathrm{vis}},{T}_{\mathrm{irr}})),\end{eqnarray} \tag{ 9 }$$

where the temperature T_vis in the viscously heated region is given by (equivalent to assuming an optically thick region with a constant opacity of 0.087 cm² g⁻¹)

$$\begin{eqnarray}&&\begin{array}{l}{T}_{\mathrm{vis}}\simeq 200\,{\rm{K}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{3/10}{\left(\displaystyle \frac{{\dot{M}}_{\star }}{{10}^{-8\,}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{2/5}\\ \,\times {\left(\displaystyle \frac{{\alpha }_{{\rm{t}}}}{{10}^{-3}}\right)}^{-1/5}{\left(\displaystyle \frac{r}{1\,\mathrm{au}}\right)}^{-9/10},\end{array}\end{eqnarray} \tag{ 10 }$$

and the temperature T_irr in the irradiation-dominated region is

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}\simeq 150\,{\rm{K}}{\left(\displaystyle \frac{{L}_{\star }}{{L}_{\odot }}\right)}^{2/7}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-1/7}{\left(\displaystyle \frac{r}{1\,\mathrm{au}}\right)}^{-3/7}.\end{eqnarray} \tag{ 11 }$$

In this study, we are interested in stars to which microlensing surveys are most sensitive so we set M_⋆ = 0.5 M_⊙. Adopting the ${\dot{M}}_{\star }\propto {M}_{\star }^{2}$ dependence that is supported by observations and the median ${\dot{M}}_{\star }={10}^{-8}\,{M}_{\odot }$ yr⁻¹ for a Sun-like star (L. Hartmann et al. 2016; C. F. Manara et al. 2023), we set ${\dot{M}}_{\star }=2.5\times {10}^{-9}\,{M}_{\odot }$ yr⁻¹ for our 0.5 M_⊙ star. The stellar luminosity's dependence on stellar mass varies from $d\mathrm{ln}{L}_{\star }/d\mathrm{ln}{M}_{\star }\simeq 1.4-1.9$ during the pre-main-sequence phase and we adopt a value of 1.5 along with L_⋆ = L_⊙ for a Sun-like star (J. Choi et al. 2016; A. Dotter 2016). With these choices and our fiducial α_t = 10⁻⁴, the disk transitions from T_vis to T_irr at 1.5 au.

We calculate the fraction of stars with masses in the range of 0.4–0.6 M_⊙ that have planets in a given mass bin at distances of 1–10 au following the procedure laid out in Y. Chachan & E. J. Lee (2023). Using the pebble accretion efficiency $\bar{\epsilon }$, we estimate the dust mass ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ needed to grow planets from an initial mass M₀ to a final mass in the range $[{M}_{{\rm{p}},\min }$, ${M}_{{\rm{p}},\max }$) where ${M}_{{\rm{p}},\min }$ and ${M}_{{\rm{p}},\max }$ delineate the limits of a given mass bin we use to construct our model planet mass function. Subsequently, we use the cumulative distribution function (CDF) of observed disk dust masses for stars in the mass range 0.4–0.6 M_⊙ to obtain the fraction of stars that have enough dust mass to grow planet seeds to the given planet mass range $[{M}_{{\rm{p}},\min }$, ${M}_{{\rm{p}},\max }$). Since masses for stars that host Class I disks have not been measured, we use stellar and disk masses for Class II stars (C. F. Manara et al. 2023) to build the CDF and rescale the disk masses by a factor of 3 to account for the larger disk masses observed in Class I stars. This factor is based on a comparison of the disk masses in the Class 0/I sample of C. F. Manara et al. (2023) and the Class II sample in J. J. Tobin et al. (2020; see Y. Chachan & E. J. Lee 2023 for more details) and is in agreement with previous studies that compared these two populations (e.g., Ł. Tychoniec et al. 2020).

Each of our planet mass bins spans 0.75 dex in log planet mass from 10⁻² M_⊕ to the pebble isolation mass M_iso evaluated at 3 au, for a total of four bins. We limit ourselves to a minimum mass of 10⁻² M_⊕ as this is the smallest mass to which Roman is expected to be sensitive (M. T. Penny et al. 2019; W. Zhu & S. Dong 2021). Since the minimum mass of interest is 10⁻² M_⊕, we choose a slightly lower fiducial seed mass log (M₀/M_⊕) = −2.06 and vary it down to −3. This range of seed mass is motivated by the results of numerical simulations of streaming instability (J. B. Simon et al. 2016; U. Schäfer et al. 2017; C. P. Abod et al. 2019; R. Li et al. 2019) combined with subsequent growth of the born planetesimals (H. Jang et al. 2022; S. Lorek & A. Johansen 2022). The planetesimal mass distribution at birth is top heavy and the largest planetesimal that is produced tends to dominate subsequent pebble accretion. Our adopted range captures the range of seed masses for a 0.5 M_⊙ star in B. Liu et al. (2020) and S. Lorek & A. Johansen (2022). We broadly classify the planet mass bins as follows: Ganymede-like planets with M_p ∈ [10⁻², 10^−1.25) M_⊕, Mars-like planets with M_p ∈ [10^−1.25, 10^−0.5) M_⊕, Earth-like planets with M_p ∈ [10^−0.5, 10^0.25), and super-Earths with M_p ∈ [10^0.25 M_⊕, M_iso). The last bin's width can therefore be smaller than 0.75 dex if M_iso < 10 M_⊕, which is always the case for our disk parameters at 3 au.

In Figure 1 (left panel), we show ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ as a function of planet mass M_p for different orbital distances. We verify that ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ rises linearly with log M_p for almost the entire range of M_p, as we would expect from Equation (7) in the 3D accretion regime. At higher planet masses (M_p ≳ 1 M_⊕) and larger orbital distances, the gradient of ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ with M_p increases slightly due to a gradual shift from 3D to 2D accretion regime. In addition, the gradient of ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ with M_p increases at larger orbital distances due to lower accretion efficiency further away from the star, requiring more dust mass to create a target M_p.

As illustrated in the right panel of Figure 1, the disk mass function is bottom heavy for stars in the mass range 0.4–0.6 M_⊙. In addition, because the required ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ within a given mass bin is linearly uniform (see the discussion below Equation (8)), our choice of logarithmic spacing of the mass bin to construct the initial planet mass function leads to the widest span in required disk mass to populate the lowest mass bin (see vertical lines in the right panel). Visually, this span of required disk mass in the lowest M_p bin appears to narrow considerably when we choose lower M₀, only because the required disk mass to grow $\mathrm{log}{M}_{0}=-2.5$ into 10⁻² M_⊕ is much higher than starting from $\mathrm{log}{M}_{0}=-2.06$ while the absolute value of difference between ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ evaluated at the bin edges remains the same.

Figure 2 shows that the fraction of stars that produce planets in a given mass bin rises as planet mass declines, largely due to the bottom-heavy disk mass function. Around a large fraction of stars, we expect planetary seeds should only grow to Ganymede and Mars mass by pebble accretion. As demonstrated in the left panel, adopting lower seed mass flattens the initial planet mass function as the fraction of stars that can nucleate Ganymede-like objects decreases substantially. Such a reduction arises because substantially more disk mass is required to grow lighter seeds to moon-mass objects, narrowing the span of ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ that corresponds to the lowest M_p bin, probing a smaller dynamic range in the disk mass CDF f(>M_disk).

Zoom In Zoom Out Reset image size

The effect of the turbulence strength α_t and fragmentation velocity v_f on the initial planet mass function is shown in Figure 2, right panel. In the 3D accretion regime, increasing α_t and decreasing v_f decreases the accretion efficiency and increases the dust mass ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ required to grow a planet (Y. Chachan & E. J. Lee 2023), leading to the analogous effect of decreasing M₀ and therefore flattening the initial planet mass function. The difference in ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ for two neighboring planet mass values also increases, although its effect on the initial planet mass function is muted as ${M}_{\mathrm{dust}\to {M}_{{\rm{p}}}}$ is already large enough to sample the tail of the disk mass CDF.

4.1. Simulation Setup

In the previous section, we have quantified the fraction of stars that harbor a planet of a given mass at birth (i.e., initial coagulation). We now consider further dynamical evolution and its effect on the final observable planet mass function. Multiple protoplanets are likely to emerge in any given protoplanetary disk, and over a gigayear timescale, the protoplanets are expected to undergo a series of orbital instabilities that can lead to mergers, shifting the population of low-mass planets to higher masses. Close encounters between planets can also lead to scattering and ejection. Whether a given encounter between a pair of planets may lead to merger or ejection can be quantified with the ratio of the escape velocity at the surface of the protoplanet versus the circular orbital velocity of the planet (C. Petrovich et al. 2014):

$$\begin{eqnarray}\begin{array}{rcl}{\theta }^{2} & = & {\left(\displaystyle \frac{{v}_{\mathrm{esc}}}{{v}_{{\rm{K}}}}\right)}^{2}\\ & = & \left(\displaystyle \frac{2\,{{GM}}_{p}}{{R}_{p}}\right)\left(\displaystyle \frac{a}{{{GM}}_{\star }}\right)\end{array}\end{eqnarray} \tag{ 12 }$$

where R_p is the radius of the planet set by the mass–radius relationship of L. Zeng et al. (2016) assuming 50:50 water-rock composition. The relation quoted by L. Zeng et al. (2016) does not extend below 0.0625 M_⊕, so for these planets we calculate planet radii assuming a constant density 1.89 g cm⁻³ (close to the observed bulk density of Ganymede), equal to the density of the 0.0625 M_⊕ planet. As demonstrated in Figure 3, ejection is more likely (θ² > 1) for scattering between higher-mass protocores and at larger orbital distances, as expected. Between 1 and 10 au, we would generally expect dynamical interactions between cores lighter than ∼0.1–1 M_⊕ to lead to mergers rather than ejections.

Zoom In Zoom Out Reset image size

We use REBOUND (H. Rein & S. F. Liu 2012) to simulate a statistical ensemble of low-mass planetary systems to quantify the extent to which dynamical evolution alters the planet mass function. As our starting point, we use our fiducial parameters (α_t = 10⁻⁴, log (M₀/M_⊕) = −2.06, v_f = 1 m s⁻¹) to set the initial conditions for the N-body simulations. We envision a scenario in which planets have some initial eccentricity that is damped by disk gas for 1 Myr, our chosen disk dissipation time, after which the gas damping is removed and the planetary system is evolved to 100 Myr. At this point, we use SPOCK (D. Tamayo et al. 2020a) to evaluate the stability probability of the planetary systems over 10⁹ orbits (1.4 Gyr). For some simulation sets that have a high probability of long-term instability, we extend their run time to 1 Gyr to verify SPOCK's output and to evaluate the system's long-term behavior by N-body integration to optimize our use of computational resources, as 1 Gyr dynamical simulations can take a while.

For each set of simulations, we specify four properties: the masses, separations, and initial eccentricities of the planets and the depletion factor of the disk gas that damps the planets' eccentricities for the first 1 Myr. We simulate 500 realizations for a given configuration in which the mean longitudes and arguments of pericenter are randomly drawn from a uniform distribution U[0, 2π) and all the planets are coplanar. For the initial eccentricities of planets (e_initial), the two bracketing values of 10⁻⁸ (physically representing zero; technically, a nonzero value is required for numerical stability) and the disk aspect ratio h_g = H_g/r are chosen.

For each mass bin, we take the median value to which we set all initial planet masses (M_initial) in a given system to be equal (note: different from the seed mass M₀, which characterizes planet mass prior to pebble accretion). Exceptionally, for the lowest mass bin, we consider a nonuniform mass in a given system. First, we sample 500 different disk masses uniformly in log M_disk from the observed CDF within range that produces planets in the lowest mass bin at 3 au (see Figure 1). For each of the 500 disk masses, we compute M_initial at chosen orbital distances following the location-dependent pebble accretion efficiency (Equation (5)).

We place the innermost planet at 1 au and subsequent planets are emplaced following a chosen initial planet separation Δ until the we reach 10 au. In our fiducial setup, we choose Δ = 3 H_g. This choice is motivated by a few considerations: (i) it ensures that the sum of planet masses is lower than the disk masses that are hosting such systems; (ii) the typical length scale of perturbations in gas disks due to disk–planet interaction is on the order of the local disk scale height (P. Goldreich & S. Tremaine 1980; D. N. C. Lin & J. Papaloizou 1986); (iii) given that we simulate systems with initial eccentricities equal to the disk aspect ratio (H_g/r), a separation of 3 H_g likely places the systems in a configuration of marginal stability rather than dynamical isolation or assured mergers. Other values of Δ are also explored. Systems that are more tightly packed and composed of more massive planets are expected to undergo more vigorous orbital instabilities on a shorter timescale (see, e.g., J.-L. Zhou et al. 2007; B. Pu & Y. Wu 2015, and references therein). For a given mass bin, if there are no mergers for a separation of 3 H_g, we run simulations for a smaller separation of 2 H_g. For the highest mass bin, we only run simulations for a planet separation of 4.5 H_g because we already observe vigorous planet merging and growth in the lower mass bin for a separation of 3 H_g.

For our fiducial disk parameters ⁷ (α_t = 10⁻⁴, ${\dot{M}}_{\star }=2.5\,\times {10}^{-9}\,{M}_{\odot }$ yr⁻¹) assuming steady-state accretion (${{\rm{\Sigma }}}_{{\rm{g}}}={\dot{M}}_{\star }/3\pi \nu$), we find

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}\approx 5300\left(\displaystyle \frac{1}{d}\right){\left(\displaystyle \frac{r}{1\,\mathrm{au}}\right)}^{-15/14}\,{\rm{g}}\,{\mathrm{cm}}^{-2},\end{eqnarray} \tag{ 13 }$$

evaluated for the passively heated regions of the disk (r ≳ 1.5 au). The resulting Σ_g is comparable to the minimum-mass extrasolar nebula: Σ_g ∼ 10⁴ g cm⁻² at 1 au (E. Chiang & G. Laughlin 2013). The factor d ∈ [1, 10, 100] quantifies the depletion factor of the disk gas, where higher d implies more significant depletion.

Eccentricity damping in the first 1 Myr is computed using the "modify_orbits" routine in REBOUNDx (D. Tamayo et al. 2020b), whereby planet eccentricities are exponentially damped on an e-folding timescale of J. Kominami & S. Ida (2002), R. I. Dawson et al. (2016):

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{{\rm{d}}} & = & \left(\displaystyle \frac{{M}_{\star }}{{M}_{{\rm{p}}}}\right)\left(\displaystyle \frac{{M}_{\star }}{{{\rm{\Sigma }}}_{{\rm{g}}}{r}^{2}}\right){\left(\displaystyle \frac{{c}_{{\rm{s}}}}{{v}_{{\rm{K}}}}\right)}^{4}{{\rm{\Omega }}}_{{\rm{K}}}^{-1}\\ & \approx & 30\,\mathrm{yr}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\oplus }}\right)}^{-1}{\left(\displaystyle \frac{r}{1\,\mathrm{au}}\right)}^{12/7}\,d,\end{array}\end{eqnarray} \tag{ 14 }$$

where the second relation has been simplified for the stellar and disk properties in the irradiated region adopted in this study (≳1.5 au). We use this simplified expression for all planets in our simulation domain (1–10 au). Our approach underestimates τ_d for planets in the viscously heated region (<1.5 au; four such planets for Δ = 3 H_g) by a factor of ∼0.55 at most. This order unity difference is negligible compared to the range of d we consider, and we check a posteriori that the initial effect of the error in τ_d on the planet eccentricities is erased by subsequent dynamical evolution. For initial eccentricities of ∼0, we find that varying the value of the depletion factor d has no notable effect on the planet eccentricities after 1 Myr. As a result, we limit ourselves to d = 1 for simulations with initial eccentricities of 0.

Planets may migrate during the disk gas phase to the inner disk and leave the region of interest (1–10 au) for microlensing observations. We quantify the importance of type-I migration for different planet masses and disk gas densities by using the prescription for migration timescale t_mig from P. Cresswell & R. P. Nelson (2008):

$$\begin{eqnarray}&&{t}_{\mathrm{mig}}=\displaystyle \frac{2}{2.7+1.1p}\left(\displaystyle \frac{{M}_{\star }}{{M}_{{\rm{p}}}}\right)\left(\displaystyle \frac{{M}_{\star }}{{{\rm{\Sigma }}}_{{\rm{g}}}{r}^{2}}\right){\left(\displaystyle \frac{{H}_{{\rm{g}}}}{r}\right)}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{-1},\end{eqnarray} \tag{ 15 }$$

and integrating $\dot{a}=-a/{t}_{\mathrm{mig}}$ over 1 Myr, where $p\,=d\,\mathrm{ln}\,{{\rm{\Sigma }}}_{{\rm{g}}}/d\,\mathrm{ln}\,r$. The migration rate is proportional to the gas surface density: a higher value of the gas depletion factor d therefore implies a slower migration rate. For nearly inviscid disks (our fiducial α_t = 10⁻⁴), feedback from piled-up gas can halt migration ⁸ (R. R. Rafikov 2002; J. Fung & E. J. Lee 2018) when planet mass reaches:

$$\begin{eqnarray}&&{M}_{\mathrm{fb}}=4\,{M}_{\oplus }{\left(\displaystyle \frac{{H}_{{\rm{g}}}/r}{0.035}\right)}^{3}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}{r}^{2}/{M}_{\star }}{{10}^{-3}}\right)}^{5/13}.\end{eqnarray} \tag{ 16 }$$

For a given planet mass bin, we limit ourselves to the gas depletion factor d that would still retain planets in the 1–10 au region in our setup. Planets in the lowest mass bins ("Ganymede-like") do not migrate significantly for d ∈ [10, 100], and for d = 1 they stay beyond 1 au as long as their initial position ≳2 au. For the next mass bin ("Mars-like"), planets do not migrate significantly for d ≥ 10. For higher-mass planets ("Earth-like" and super-Earths), the combined effect of larger t_mig and lower M_fb when d = 100 keeps the planets from migrating significantly. Table 1 lists all the parameters and simulations that are part of this study. The properties of the simulations make up their assigned name in the following order: planet mass bin, separation in terms of disk scale height, initial eccentricity, and gas depletion factor.

Table 1. Dynamical Simulations Setup

Simulation a	Minitial	einitial b	Δ	d
	(M⊕)		(Hg)
G_3_0_1	10−1.625	0	3	1
G_3_hg_1	10−1.625	hg	3	1
G_3_hg_10	10−1.625	hg	3	10 c
G_3_hg_100	10−1.625	hg	3	100 c
GVM_3_0_1	10−2−10−1.25 (at 3 au)	0	3	1
GVM_3_hg_1	10−2−10−1.25 (at 3 au)	hg	3	1
GVM_3_hg_10	10−2−10−1.25 (at 3 au)	hg	3	10
GVM_3_hg_100	10−2−10−1.25 (at 3 au)	hg	3	100
G_2_0_1	10−1.625	0	2	1
G_2_hg_1	10−1.625	hg	2	1
G_2_hg_10	10−1.625	hg	2	10
G_2_hg_100	10−1.625	hg	2	100
M_3_0_100	10−0.875	0	3	100
M_3_hg_10	10−0.875	hg	3	10
M_3_hg_100	10−0.875	hg	3	100
E_3_0_100	10−0.125	0	3	100
E_3_hg_100	10−0.125	hg	3	100
SE_4.5_0_100	100.53	0	4.5	100
SE_4.5_hg_100	100.53	hg	4.5	100

Notes.

^a Simulation names indicate the planet mass bin (G: Ganymede-like, GVM: Ganymede-like with variable planet mass, M: Mars-like, E: Earth-like, SE: super-Earth), separation in terms of scale height, initial eccentricity, and gas depletion factor. ^b When 0, it is set to 10⁻⁸ for numerical stability. Variable h_g is the disk aspect ratio. ^c Simulations run to 1 Gyr.

Download table as:  ASCIITypeset image

4.2. Simulation Results

4.2.1. Ganymede-like Planets

For our fiducial set of simulations with equal-mass planets separated by 3 H_g, we find that there are no mergers in all our simulations within 100 Myr. Stability analysis with SPOCK after 100 Myr of integration indicates the probability of collisions is very low, except for a tail of simulations with initial eccentricity e_initial = h_g and d = 100 (G_3_hg_100, Figure 4, left panel). We extend the integration time of these simulations to 1 Gyr and find that only 14.8% of our 500 simulations have one or more pairwise mergers of the initial planets (Figure 5). We also integrate the simulation set G_3_hg_10 to 1 Gyr to verify SPOCK's results and find that there are no mergers in this case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even in the case of the small fraction of planets undergoing mergers, they are all pairwise, increasing the mass of the resulting planet by a factor of 2, and so these planets, after merging, do not jump to the next mass bin. At a separation of 3 H_g, Ganymede-mass planets therefore do not grow significantly by mergers regardless of the chosen value of e_initial and d.

We also test if having variable planet masses rather than equal-mass planets separated by 3 H_g changes the outcome of dynamical evolution (see Section 4.1 for details on setting the planet mass). After running the simulations for 100 Myr, we find that there are no mergers in any of the simulations. Stability analysis with SPOCK again indicates that collision probability in these systems over 10⁹ orbits is very low (Figure 4, middle panel). Having variable rather than equal-mass planets does not change the outcome of dynamical evolution, likely because all protocores, even if their masses were varied, are small (≲0.2 M_⊕). Ganymede-like planets stay as such over gigayear timescales if they are separated by ≥3 H_g. In addition, since these planets are small enough to undergo a negligible amount of migration, we expect M dwarfs to be teeming with Ganymede-like objects at 1–10 au, at least for initial Δ ≥ 3 H_g.

Finally, to determine the fate of these planets for tighter orbital separations, we run simulations for equal-mass planets spaced apart by 2 H_g (Figure 6). After 100 Myr of integration, we find that all of the simulations in the dynamically hottest setup (G_2_hg_100) have had at least one merger. For e_initial = h_g and d = 10, 35.2% of our 500 simulations show at least one merger. For d = 1, gas damping is still effective enough to lower planet eccentricities and prevent mergers in the vast majority of our simulations (no mergers for G_2_0_1, only one merger in G_2_hg_1). Results from SPOCK indicate that most systems have a collision probability higher than 0.5 (Figure 4, right panel). The mergers in simulation set G_2_hg_10 are all pairwise, and so the mass of the resulting planet is just twice the initial mass. However, for G_2_hg_100, 169 planets (1.68% of all planets at 100 Myr in 500 simulations) in 103 different systems (20.6% of 500 simulations) end up with masses ≥3 × M_initial, i.e., they end up crossing the threshold to the next mass bin and would be classified as Mars-like planets (Figure 6). The fraction of systems in which we expect the member planets to jump to the next mass bin by mergers over 100 Myr is still too small to significantly alter the shape of the planet mass function, especially for initially large seed mass where the initial fraction of stars harboring Ganymede-like planets is significantly higher than that of the Mars-like planets. For small seed mass M₀ where the initial planet mass function is more flat, the dynamical sculpting would likely be more pronounced and even cause a more peaked final mass function. Fully simulating the orbital dynamics over 1 Gyr and longer would be required to determine quantitatively the changes to the planet mass function.

Zoom In Zoom Out Reset image size

4.2.2. Mars-like Planets

Since Ganymede-like planets already show vigorous merging when they are separated by Δ = 2 H_g and d = 100, we expect Mars-like planets to undergo dynamical instabilities and merge for this separation, too. We therefore set Δ = 3 H_g and larger values of d = 10 and 100 to limit ourselves to the case where migration has a negligible impact on the final location of Mars-like planets. After 100 Myr of evolution, we find that none of the systems in M_3_0_100 and M_3_hg_10 undergo any mergers. For M_3_hg_100, one simulation has two pairwise mergers and another simulation has one pairwise merger, with no mergers seen for any of the other simulations.

Estimates of the collision probability from SPOCK (Figure 7) for e_initial = h_g indicate that 90% of M_3_hg_10 simulations and 86% of M_3_hg_100 simulations have a collision probability <0.5. Counterintuitively, only 58% of M_3_0_100 simulations have a collision probability <0.5 despite the dynamically colder initial conditions (e_initial = 0). We find that the eccentricities of planets in the M_3_0_100 simulation suite rapidly rise after gas disk dissipation to match the eccentricities of planets in the M_3_hg_10 simulations, at times even exceeding them for the closer-in planets. Clearly, Mars-like planets dynamically stir each other up for separations Δ = 3 H_g.

Zoom In Zoom Out Reset image size

To investigate why the M_3_0_100 simulations have higher collision probabilities, we examine which "features" of the system are driving the SPOCK predictions. SPOCK calculates a number of metrics (e.g., closeness to mean-motion resonances, ratio of planet eccentricities to orbit crossing eccentricities, chaos indicators) for a given planetary system by considering nearest planets in sets of three and by taking the maximum value of the metrics (i.e., the lowest stability probability) among all such sets. We find that the high collision probabilities for the M_3_0_100 simulations are primarily driven by the chaos indicator MEGNO and MEGNOstd (P. M. Cincotta et al. 2003), computed based on the 10⁴ yr of orbital integration. MEGNO measures chaos on short timescales by calculating the rate at which nearby orbits diverge. We find that, for systems with e_initial = 0, the rise in eccentricity is steeper with time for more of their member planets compared to simulations with e_initial = h_g (see also Figure 5 in R. I. Dawson et al. 2016). This likely drives the MEGNO chaos indicator to higher values.

Systems flagged as chaotic by MEGNO can nonetheless take a long time to go unstable depending on the time required for planet eccentricities to diffuse to orbit crossing values. To test if these systems are chaotic but not necessarily catastrophically unstable, we extend the run time of our M_3_0_100 simulations beyond 100 Myr. For 172 simulations with collision probability <0.5 and 226 simulations with collision probability >0.5 with run times in the range 250–345 Myr, no mergers are observed, so at least on a few 100 Myr timescale these systems are stable. Whether they remain so over 1 Gyr requires longer simulations. Similar to Ganymede-like objects, more than pairwise mergers are required for dynamical evolution to cause a significant change to the planet mass function between Mars-like and Earth-like mass bins.

4.2.3. Earth-like Planets

Both Ganymede-like and Mars-like planets do not undergo large-scale mergers for separations Δ ≥ 3 H_g. We therefore choose a value of Δ = 3 H_g for Earth-like planets to determine if dynamical evolution at such planet separations significantly affects their population. Figures 8 and 9 show the collision probability calculated from SPOCK and the planet multiplicity, mass, and semimajor axis distribution after 100 Myr. We exclude simulations (five from E_3_hg_100, one from E_3_0_100) in which a planet gets scattered to semimajor axis <0.1 au as we do not evolve close-in planets that would prohibitively extend the simulation timescale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 8 and 9 show that the state of these systems after 100 Myr is weakly dependent on the initial eccentricities of the planets, suggesting a high degree of orbital instabilities and series of mergers in both cases. The e_initial = h_g simulations have slightly more systems at higher multiplicity, implying that these systems are marginally dynamically colder, in agreement with the collision probabilities CDF, as well. We see that the final planet multiplicity peaks at 4–7 planets per system compared to the initial multiplicity of 18, with the most common planet mass being twice the initial mass. For E_3_0_100 simulations, 1598 of the 2935 planets (54.4%) remaining at 100 Myr in 498 out of 499 systems have masses ≥3 M_initial, i.e., they would be qualified as super-Earths, jumping to the next mass bin. Similarly, for E_3_hg_100 simulations, we find that 1545 of the 3031 planets (50.1%) remaining at 100 Myr in 490 out of 495 systems have masses ≥3 M_initial. Further evolution will likely increase the fraction of Earth-like planets that would jump to the next mass bin (39.2% of E_3_0_100 and 36.0% E_3_hg_100 simulations have collision probability >0.5). Earth-like planets separated by Δ = 3 H_g are therefore likely to be heavily dynamically sculpted, potentially leading to a gap in the final planet mass function at ∼Earth mass. The mass at which such a gap appears, if observed, can inform the initial orbital separations. Observation of a significant population of Earth-like planets would imply wider initial separations. The gap in the mass function at Earth-like planet masses could partially be filled in by mergers of Mars-like planets if they are born separated by Δ ≲ 2 H_g.

The bottom panel of Figure 9 shows the mass and orbital semimajor axis of all bound planets at 100 Myr (note: it does not show which planets are in the same system). Low-mass planets are scattered out far beyond 10 au on highly eccentric orbits, at times to such large distances that the planets would appear to be free floating although physically bound. We also find truly unbound planets with eccentricities >1. In the E_3_0_100 simulations, 149 out of 2935 (5.1%) total planets at 100 Myr are unbound, of which 86 retain their initial mass. The corresponding number for E_3_hg_100 simulations is 121 out of 3031 (4%) planets, with 64 planets that retain their initial mass. The observed dynamical interactions of Earth-like planets that lead to ejections and large-scale scattering is in line with expectations from Figure 3.

4.2.4. Super-Earths

Since Earth-like planets merge vigorously for separation Δ = 3 H_g, we simulate super-Earths with a larger separation of Δ = 4.5 H_g. This separation lies in between the bracketing values obtained by scaling Δ for super-Earths from Earth-like planets by using the relations for dynamical separation from S. Hadden & Y. Lithwick (2018; ${\rm{\Delta }}\propto {M}_{{\rm{p}}}^{1/4}$) and separation in terms of Hill radius (${\rm{\Delta }}\propto {M}_{{\rm{p}}}^{1/3}$; B. Gladman 1993; A. C. Petit et al. 2018). In the following analysis, we again exclude simulations (26 from SE_3_hg_100, 17 from SE_3_0_100) in which a planet gets scattered to semimajor axis <0.1 au. Figures 10 and 11 show the collision probability calculated from SPOCK and the planet multiplicity, mass, and semimajor axis distribution at 100 Myr. We see that the memory of initial eccentricity is even more erased as compared to the Earth-like planet case, suggesting vigorous dynamical sculpting by planet–planet interactions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This vigorous sculpting is evidenced by the sharper peak of planet multiplicity at ∼4 (top panel of Figure 11) compared to Earth-like planets at Δ = 3 H_g. Like Earth-like planet case, most mergers are pairwise leading to the modal planet mass being twice the initial mass after 100 Myr. For SE_3_0_100 simulations, we find that 972 of the 2168 planets (44.8%) remaining at 100 Myr in 471 out of 483 systems have masses ≥3 M_initial, large enough to jump to the next mass bin. Similarly, for SE_3_hg_100 simulations, 980 of the 2120 planets (46.2%) remaining at 100 Myr in 467 out of 474 systems have masses ≥3 M_initial. Mergers of super-Earths can therefore easily produce planets that are akin to Neptune. Further evolution will likely increase the fraction of super-Earths that grow beyond M_iso and evacuate the initial mass bin (44.9% of systems for sets of simulations have collision probability >0.5). Combined with most systems with Earth-like planets undergoing mergers, we would expect an overall shift of the initial mass function for Earth mass and up toward higher masses, unless the initial Δ > 4.5 H_g.

The mass and orbital semimajor axis of all bound planets at 100 Myr are shown in the bottom panel of Figure 11 (note: it does not show which planets are in the same system). Similar to the case for Earth-like planets, super-Earths that retain their initial mass are scattered out far beyond 10 au on highly eccentric orbits and would ostensibly be considered free floating. The fraction of unbound planets is higher for super-Earths, as one might expect from Figure 3. For the SE_3_0_100 simulations, 244 out of 2168 (11.3%) total planets at 100 Myr are unbound, of which 165 have a mass = M_initial. For the SE_3_hg_100 simulations, 224 out of 2120 (10.6%) planets are unbound, of which 157 planets have a mass of M_initial.

The Nancy Grace Roman Space Telescope will enable us to find Ganymede-mass to Earth-mass objects at 1–10 au for the first time. The sensitivity of the telescope to the lowest planet masses is, however, low, and therefore the number of Ganymede and Mars-like objects that we expect to find is highly sensitive to the underlying planet mass function.

Combining mass growth by pebble accretion with the observed disk mass function, we find that:

• 1.  
  The initial planet mass function is generally bottom heavy and increasingly so for larger initial seed mass approaching that of Ganymede.
• 2.  
  Consequent orbital evolution for a few 100 Myr leads to minimal change in the mass function for M_p ≤ 10^−0.5 M_⊕ and vigorous instabilities and mergers that lead to a potential gap at ∼1 M_⊕, if the planets were initially spaced at Δ = 3 H_g. The location of this gap is expected to shift to lower mass for a tighter initial separation and to higher mass for a wider initial separation.
• 3.  
  For planets with M_p ≳ M_⊕, we expect a small fraction of planets (∼4%–5% Earth-like and ∼11% for super-Earths among the total simulated planets) being scattered to wide orbits (≳100 au) or becoming unbound (with the unbound population dominating over the wide-orbit population by 2 orders of magnitude at least at 100 Myr), contributing to the population of free-floating planets.

Our findings suggest that the shape of the observed planet mass function can be used to infer the initial orbital architecture of the planetary systems. A more steeply bottom-heavy mass function would imply a higher initial seed mass, and if we find a signature of a gap in the mass function the location of the gap can be leveraged to infer the initial planet–planet spacing.

Ejection from bound orbits by dynamical scattering is easier for systems with higher-mass objects (see Figure 3). Given our setup of the problem where all planets within a given system start with the same mass, we would expect the resulting mass function of free-floating planets to be either top heavy or flat beyond ≳1 M_⊕ (with a sharp decline below), with the exact mass being dependent on the initial planet separations. Measurements of the relative abundance of bound and free-floating planets and the planet mass function of free-floating planets by Roman will therefore serve as a probe of the primordial dynamical separation of massive (≳1 M_⊕) planets.

Our work generally expects a bottom-heavy bound planet mass function and top-heavy free-floating planet mass function, which is potentially at odds with the current ground-based microlensing surveys (see T. Sumi et al. 2023, their Figure 6). Under our calculations, starting with a significantly lighter initial seed mass followed by dynamical evolution can reduce the population of Ganymede-like objects compared to that of Mars-like objects, driving the bound planet mass function to be more top heavy. Tighter initial orbital spacing (Δ < 3 H_g) could further evacuate the population of Mars-like objects, adding to the population of their more massive counterparts. Creating a bottom-heavy free-floating planet mass function may be possible if the systems with Ganymede-like and Mars-like objects also harbor more massive planets (≳1 M_⊕) and/or are orbited by a binary companion that can dynamically eject these small objects (G. A. L. Coleman & W. DeRocco 2024). Given that the observed disk mass fraction is bottom heavy, most disks that can nucleate objects of ∼Mars mass and below do not have enough material to also create Earth-like and more massive planets. It is possible that such low-mass disks do not even create the initial seed mass such that they should not be considered as planet-forming disks (see, e.g., R. Li & A. N. Youdin 2021 for a condition on the local dust-to-gas ratio to trigger clumping). In this case, the large population of free-floating low-mass objects would have originated from more massive disks as part of a family with more massive planets. Whether the free-floating planet mass function is truly bottom heavy or not remains to be verified with space-based missions.

We thank the anonymous referee for a report that improved the manuscript. This research was enabled in part by support provided by Calcul Québec and the Digital Research Alliance of Canada (alliance.can.ca). Y.C. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through the CITA National Fellowship and the Trottier Space Institute through the TSI Fellowship. Y.C. is grateful for discussions with Eugene Chiang, Sam Hadden, and Yanqin Wu at the 2023 CITA Planet Day Meeting that inspired this work. E.J.L. gratefully acknowledges support by NSERC, by FRQNT, by the Trottier Space Institute, and by the William Dawson Scholarship from McGill University.

Software: astropy (Astropy Collaboration et al. 2013, 2018), REBOUND (H. Rein & S. F. Liu 2012), REBOUNDx (D. Tamayo et al. 2020b), SPOCK (D. Tamayo et al. 2020a), Matplotlib (J. D. Hunter 2007), Numpy (C. R. Harris et al. 2020).