A Theoretical Framework for the Mass Distribution of Gas Giant Planets Forming through the Core Accretion Paradigm

The mass distribution of fundamental objects represents an important feature of any astrophysical discipline, including studies of planets, stars, galaxies, dark matter halos, and galaxy clusters. Moreover, a complete understanding of the issue requires both the observational specification of the mass distribution and a predictive theoretical framework. The stellar initial mass function (starting with Salpeter 1955), for example, affects a large number of astronomical issues, including galactic chemical evolution, supernova rates, and feedback during star formation. With thousands of exoplanets now discovered (e.g., see Han et al. 2014; Gaudi et al. 2020), we can start to address the planetary mass function (PMF), which will ultimately inform studies of planetary system architectures, planetary habitability, and many other issues. In addition, the PMF will provide a consistency check on theories of planet formation, once they are sufficiently developed. The observational determination of this mass distribution is now coming into focus, although the database remains incomplete, and many uncertainties persist. On the other hand, our theoretical understanding of the mass function remains in its infancy. The goal of this paper is to construct simple but physically motivated models for the PMF for giant planets, i.e., companions in the mass interval M_p = 0.1–10M_J (where M_J is the mass of Jupiter).

For the range of planetary masses of interest, most objects are thought to form via the core accretion paradigm (Pollack et al. 1996), which is organized into three phases (see, e.g., Benz et al. 2014 for a review). In the context of this mechanism, the accumulation of small rocky bodies (phase I) leads to the production of planetary cores with masses of order 10M_⊕, which takes place over a timescale ${t}_{\mathrm{core}}$. Note that core formation could take place either through two-body accumulation of icy planetesimals (Hansen & Murray 2013) or through rapid accretion of small pebbles assisted by gas drag (Ormel & Klahr 2010; Lambrechts & Johansen 2012). In any case, the cores accrete gas from the surrounding disk after reaching their critical mass threshold. During the subsequent phase II, growth is relatively slow and is controlled by cooling processes (e.g., opacity effects). After the total mass reaches roughly twice the core mass (mass scale M_II) at time t_onset, gas accretion occurs rapidly, and the planets accumulate most of their mass. This phase III lasts until the planets reach their final masses in the range of M_p ≈ 0.1–10M_J, and then terminate their growth at time t_end. The final mass of a given planet is thus given by the integral expression

$$\begin{eqnarray}&&{M}_{{\rm{p}}}={M}_{\mathrm{II}}+{\int }_{{t}_{\mathrm{onset}}}^{{t}_{\mathrm{end}}}{\dot{M}}_{{\rm{p}}}(t){dt}={M}_{\mathrm{II}}+\langle {\dot{M}}_{{\rm{p}}}\rangle ({t}_{\mathrm{end}}-{t}_{\mathrm{onset}}),\end{eqnarray} \tag{ 1 }$$

where ${\dot{M}}_{{\rm{p}}}(t)$ is the accretion rate onto the planet as a function of time, and where the final equality follows from the mean value theorem. Within the context of the core accretion paradigm, the mass accretion rate onto the planet is some fraction of the total mass accretion rate ${\dot{M}}_{{\rm{d}}}$ through the disk, so that $\langle {\dot{M}}_{{\rm{p}}}\rangle$ = $\langle f{\dot{M}}_{{\rm{d}}}\rangle \equiv$ $\epsilon \langle {\dot{M}}_{{\rm{d}}}\rangle$, where the final equality defines a time-averaged efficiency factor . Since circumstellar disks have lifetimes of ∼1–10 Myr, it is natural to associate the end of planetary accretion (t_end) with the disappearance of the disk. One objective of this paper is to test the plausibility of this hypothesis.

As written, Equation (1) is exact and depends only on the four variables (M_II, $\langle {\dot{M}}_{{\rm{p}}}\rangle$, t_onset, t_end). In other words, for a particular planet-forming event, if we know the mass M_II at the onset of gas accretion, the starting and ending times (t_onset, t_end), and the mean accretion rate $\langle {\dot{M}}_{{\rm{p}}}\rangle$, the final mass of the planet is completely determined. If, in addition, we know the distributions of the variables over the collection of the planet-forming disks of interest, then the distribution of planetary masses is also determined. Unfortunately, the mass accretion history of the planet, along with the timescales t_onset and t_end, depend on a large number of physical processes that remain under study, including disk viscosity, which drives accretion, the fraction of material moving through the disk captured by the growing planet, disk and planetary magnetic fields, planetary migration, the role of circumplanetary disks, and many others. In a complete theory, where all of the relevant subprocesses are known and predictable, one could calculate t_onset, t_end, and ${\dot{M}}_{{\rm{p}}}(t)$, and their distributions, from first principles. In the absence of such a complete understanding, the more modest goal of this paper is to construct observationally motivated distributions of the relevant input variables, as encapsulated via Equation (1) and use the results to predict the PMF. We can then determine the properties of these input distributions required to produce PMFs that are consistent with observations. In addition, the framework developed herein can be readily generalized in future work.

Although no rigorous, a priori theory of the PMF has been established, a substantial amount of previous work has been carried out. One approach is to use population synthesis models, wherein planet formation is simulated using a large number of physical processes that are treated in an approximate manner (e.g., see Ida & Lin 2008; Benz et al. 2014). With the proper choice of input parameters, these models produce mass distributions for planets that are roughly consistent with the observed PMF (Thommes et al. 2008; Mordasini et al. 2009), although they tend to underproduce sub-Saturn mass planets (Suzuki et al. 2018). The approach of this paper is complementary to population synthesis models. Instead of including as many physical processes as possible, our goal is to find the minimal level of complexity necessary to reproduce the observed PMF. This present approach is semiempirical in that it relies on observational input to help specify the distributions of input parameters (see Adams & Fatuzzo 1996 for an analogous approach to the stellar initial mass function). In addition, this minimalistic approach allows for PMF expressions to be obtained semianalytically. Other related studies include a statistical model based on orbital spacing (Malhotra 2015) and the role of dynamical instabilities in shaping the PMF (Carrera et al. 2018).

The scope of this paper is limited to gaseous giant planets forming through the process of core accretion, in particular, those objects for which the gas accretion phase determines the final planetary mass. As a result, this approach is limited to planets forming in the mass range of 0.1M_J ≤ M_p ≤ 10M_J. The current observational sample, especially from the Kepler mission (Borucki et al. 2010; Batalha et al. 2011), indicates that the total planet population includes a large number of smaller bodies, with a preference for superearths with masses M_p ∼ 10M_⊕ (e.g., Zhu et al. 2018). These objects are roughly analogous to the cores of the giant planets, but they are thought to contain relatively little gas, and their mass distribution is not addressed here. These lower mass planets could represent a second population (Pascucci et al. 2018; see also Schlaufman 2015), which should be addressed in future studies. In addition, this present analysis is confined to those star/disk systems where giant planets are able to form. The probability that a given system will produce giant planets is a separate but important issue (Howard et al. 2010; M. R. Meyer et al. 2020, in preparation).

This paper is organized as follows. Section 2 discusses the observational constraints utilized in our models, including a brief overview of the observed PMF, distributions of disk lifetimes, constraints on core formation timescales, and the dependence of these quantities on the mass of the host star. Section 3 presents a framework to construct the PMF. A model using the observed distribution of disk lifetimes, the expected dependence of the mass accretion rate on the Hill radius, and an additional random variable can reproduce the nearly power-law mass function that is observed. Since disk properties, and hence planet properties, depend on the mass of the central star, Section 4 explores the effects of varying the stellar mass. The paper concludes in Section 5 with a summary of our results and discussion of their implications. For completeness, the paper includes an Appendix that considers an alternative model for the PMF.

This section outlines the observational input used to inform and test our models of the PMF. We start with current constraints on the observed PMF. Unfortunately, the observed mass distribution function for gas giant planets is not fully specified with current data. Studies carried out to date for FGK stars (starting with Cumming et al. 2008; see also Mayor et al. 2011) suggest that the mass distribution has an approximately power-law form over the mass range of interest 0.1M_J ≤ M_p ≤ 10M_J, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{dF}}{{{dM}}_{{\rm{p}}}}=\displaystyle \frac{A}{{M}_{{\rm{p}}}^{p}}\qquad \mathrm{where}\qquad p\approx 1.3\pm 0.2,\end{eqnarray} \tag{ 2 }$$

where A is a normalization constant (and more data are needed at the ends of the mass range). Although significant uncertainties remain, this paper uses the power-law mass distribution of Equation (2) as its basis for comparison to observations. Note that this form for the mass function is only applicable over the specified range of planetary masses. The distribution of planet radii displays a break near R_p = 4R_⊕ (Schlaufman 2015), corresponding to planet mass M_p = 10–15M_⊕ (Wolfgang et al. 2016), which falls just below the masses of interest here. Extending the range to smaller masses, down to M_p ∼ 1M_⊕, the planetary radius distribution has a bimodal form (e.g., Fulton et al. 2017). Another complication is that planets of a given mass can display a distribution of radii (Otegi et al. 2020). Current results from microlensing surveys (Shvartzvald et al. 2016; Suzuki et al. 2016) also find power-law distributions for the ratios of companion masses q = M₂/M_* over the range of 10⁻⁴ < q < 10⁻², with slopes p ≈ 0.9–1.5, roughly consistent with Equation (2). Keep in mind that these surveys primarily sample host stars of lower mass (which are the most common stars). For completeness, we note that the observed mass distribution for wide-orbit planets with masses M_p > 3M_J follows a similar power law (Schlaufman 2018; Wagner et al. 2019), although the Gemini Planet Imager Exoplanet Survey hints at a steeper slope p ≈ 2.3 ± 0.8 (Nielsen et al. 2019).

One should keep in mind that the observed PMF is derived from planets detected around older main-sequence stars, and this distribution could have evolved from earlier epochs immediately after planet formation. For example, planet–planet interactions can eject planets from their solar systems, and migration through the disk can lead to accretion of planets by their host stars. Since these processes are unlikely to be independent of planetary mass, their action will sculpt the PMF over time. These processing effects, and perhaps others, should thus be included in the determination of the observed PMF. In other words, we ultimately need to make the distinction between the initial PMF and the observed PMF (an analogous complication is well known in considerations of the stellar mass distribution).

The disk lifetime sets an upper limit on the time available for gas accretion onto planets. The disks associated with young stars generally have lifetimes in the range of 1–10 Myr (Meyer et al. 2007; Williams & Cieza 2011; Hartmann et al. 2016; Manzo-Martínez et al. 2020). More specifically, the distribution of disk lifetimes dF/dt is often taken to have the general form

$$\begin{eqnarray}&&\displaystyle \frac{{dF}}{{dt}}=\displaystyle \frac{1}{\tau }\exp (-t/\tau ),\end{eqnarray} \tag{ 3 }$$

which is consistent with observations of disk fractions in star-forming regions of varying ages (Haisch et al. 2001; Hernández et al. 2007; Mamajek 2009; Fedele et al. 2010). Note that different signatures are used to indicate the presence of disks, including signs of active gas accretion onto the stars, near-infrared excess emission at 1–5 μm (tracing the inner disk), infrared emission at 5–100 μm (tracing disk scales 0.3–30 au), and millimeter-wave observations of dust and gas (tracing the outer disks at 10–100 au). As written, the distribution of disk lifetimes is normalized such that

$$\begin{eqnarray}&&{\int }_{0}^{\infty }\displaystyle \frac{{dF}}{{dt}}{dt}=1.\end{eqnarray} \tag{ 4 }$$

Observations indicate that the timescale τ ≈ 4–5 Myr. Note that τ represents the time required for the population of disks to decrease by a factor of e ≈ 2.72. The corresponding "half-life," the time required for the population to decrease by a factor of 2, is only about ${\tau }_{1/2}=\tau \mathrm{ln}2\sim 3\,\mathrm{Myr}$.

Observations also indicate how disk properties vary with the mass of the host star. Disk lifetimes are observed to decrease with increasing stellar mass (Hillenbrand et al. 1998; Carpenter et al. 2006; Yasui et al. 2014; Ribas et al. 2015; Yao et al. 2018). This observational finding can be incorporated using a scaling law of the form

$$\begin{eqnarray}&&\tau =\displaystyle \frac{{\tau }_{1}}{m},\end{eqnarray} \tag{ 5 }$$

where τ₁ ∼ 5 Myr, and where we have defined a dimensionless stellar mass variable

$$\begin{eqnarray}&&m\equiv \displaystyle \frac{{M}_{* }}{1\,{M}_{\odot }}.\end{eqnarray} \tag{ 6 }$$

Mass accretion rates ${\dot{M}}_{{\rm{d}}}$ for the disks surrounding T Tauri stars depend on both time and the mass of the central star. These dependences can be modeled with a function of the form

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{d}}}={\dot{M}}_{d(0)}\displaystyle \frac{{m}^{2}}{{\left(1+t/{t}_{0}\right)}^{3/2}},\end{eqnarray} \tag{ 7 }$$

where the timescale t₀ ∼ 1 Myr and where the time dependence follows from α-disk models (e.g., Hartmann 2009). The observed m² dependence (in the numerator) follows from observational results (see the review of Hartmann et al. 2016). Although no simple explanation is currently accepted for this law, it follows if disk lifetimes decrease with stellar mass (Equation (5)), and the starting disk mass M_d ∝ M_* increases with stellar mass (see also Dullemond et al. 2006; Hartmann et al. 2006). This latter scaling law has been observed in a number of surveys (e.g., Natta et al. 2007; Andrews et al. 2013; Ansdell et al. 2016), with the slope somewhat steeper than linear in some studies (e.g., Barenfeld et al. 2016; Pascucci et al. 2016). In general, the observations find a well-defined correlation between submillimeter disk luminosity and stellar mass, whereas the conversion from the observed quantities to the inferred relationship M_d ∝ M_* requires disk modeling. In addition, the observations correspond to dust emission, so that the dust-to-gas ratio must be assumed in order to obtain estimates for the total disk mass. For completeness, note that a weaker correlation is found for actively evaporating disks in the Orion Nebula Cluster (Eisner et al. 2018).

Finally, the formation time ${t}_{\mathrm{core}}$ for the cores of giant planets, and hence the time t_onset before the start of rapid gas accretion, is also expected to depend on the mass of the central star. Here we adopt a scaling law of the form

$$\begin{eqnarray}&&{t}_{\mathrm{onset}}={t}_{{\rm{c}}0}\,{m}^{1/2},\end{eqnarray} \tag{ 8 }$$

where we expect that t_c0 ∼ 1 Myr (Pollack et al. 1996; Benz et al. 2014). This scaling law follows from theoretical considerations of solid accumulation (Safronov & Zvjagina 1969) combined with empirical scaling laws. We start by assuming that planetary cores preferentially form just outside the water ice line, where the surface density of solids is expected to increase. Using the mass–luminosity relation for pre-main-sequence stars, ³ ${L}_{* }\sim {M}_{* }^{\mu }$, where μ ≈ 5/3 − 2, and a fixed sublimation temperature T_ice, we can find the semimajor axis of the ice line from energy balance. Since ${T}_{\mathrm{ice}}^{4}\propto {L}_{* }/{a}^{2}\propto {M}_{* }^{2}/{a}^{2}$ = constant, we find that a_ice ∝ M_*. The larger radius of the formation site—for larger stars—leads to slower accumulation of the cores. The Safronov accumulation time (specifically, the doubling time in the absence of gravitational focusing) depends on the surface density Σ and orbital frequency Ω (e.g., Safronov 1977) according to t_saf ∝ Σ⁻¹Ω⁻¹. Using a surface density Σ ∝ M_*/a^1/2 (e.g., Andrews et al. 2009), along with a Keplerian rotation curve, we find the result given in Equation (8). Note that steeper surface density profiles lead to more sensitive dependence of the core formation time on stellar mass.

The particular scaling law (Equation (8)) was derived for the case of classic planetesimal accretion. For the competing picture wherein pebble accretion accounts for the formation of giant planet cores (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012; Lin et al. 2018; Rosenthal et al. 2018), the scaling exponent can be different. For example, Lambrechts & Johansen (2014) found that the pebble accretion rate depends on stellar mass according to $\dot{M}\sim {M}_{* }^{-11/36}$, which coresponds to the alternate scaling law t_onset ∼ m^11/36 ∼ m^1/3 (compare with Equation (8)). The goal of this present work is to explore the effects due to the core formation time increasing with stellar mass, but the exact scaling is not as important. For completeness, we note that pebble accretion can occur over a wider range of radial locations in the disk, compared with the classical picture. ⁴ Nonetheless, studies suggest that pebble accretion is more efficient beyond the ice line, as implicitly assumed here, due to the greater mass in pebbles and higher probability of sticking (e.g., Bitsch et al. 2015; Chambers 2016).

This section constructs a working model for the PMF based on the empirical considerations of the previous section. The models constructed here are independent of the mass of the host star. The effects of varying stellar mass on the resulting PMFs are addressed in Section 4.

3.1. Distribution of Gas Accretion Timescales

Observations provide estimates for the distribution of disk lifetimes. For purposes of determining planetary masses, however, we need the distribution of time remaining after the onset of rapid gas accretion. The starting time t_onset is determined by the core formation timescale and the slow cooling phase. It takes into account the fact that rapid gas accretion onto the planet occurs only after the body has built up a critical mass, which includes both the rocky core and enough additional gas to become self-gravitating. For a given stellar M_*, the time to build the core t_onset depends on the surface density of solids, and the cooling time depends on the gas opacity in the planetary atmosphere (e.g., see Pollack et al. 1996; Benz et al. 2014; Helled et al. 2014; Pisoet al. 2015, and references therein).

For the case of interest here, where disk lifetimes are exponentially distributed, the remaining accretion time also follows an exponential law. Consider a given value t_onset of the starting time and define a new effective time variable $\tilde{t}$ according to

$$\begin{eqnarray}&&\tilde{t}=t-{t}_{\mathrm{onset}},\end{eqnarray} \tag{ 9 }$$

which resets the zero-point of time for gas accretion. The time t is the disk lifetime, which represents the time when accretion ends, and which follows the exponential distribution of Equation (3). For any value of t_onset, the probability distribution of the corrected time variable $\tilde{t}$ is determined by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dF}}{d\tilde{t}} & = & \displaystyle \frac{{dF}}{{dt}}\displaystyle \frac{{dt}}{d\tilde{t}}=\displaystyle \frac{1}{\tau }\exp [-(\tilde{t}+{t}_{\mathrm{onset}})/\tau ]\\ & = & \displaystyle \frac{\exp [-{t}_{\mathrm{onset}}/\tau ]}{\tau }\exp [-\tilde{t}/\tau ],\end{array}\end{eqnarray} \tag{ 10 }$$

where dF/dt is the distribution for the original time variable t. The probability distributions for both t and $\tilde{t}$ thus have the same exponential form. As written, however, this expression for the probability distribution (for $\tilde{t}$) is not normalized. Its integral over all $\tilde{t}\gt 0$ is less than unity because the leading coefficient takes into account the fact that those disks with lifetimes t_disk < t_onset will not produce gas giant planets. For purposes of calculating probability distributions for the time remaining after the onset of rapid gas accretion, however, the normalized form of the distribution must be used, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{dF}}{d\tilde{t}}=\displaystyle \frac{1}{\tau }\exp [-\tilde{t}/\tau ].\end{eqnarray} \tag{ 11 }$$

After disks with short lifetimes are removed from the sample, the remaining disk population will "continue to decay" with an exponential law that has the same half-life as before (analogous to the case of radioactive decay).

If all of the disks produced planetary cores on the same timescale, then the factor $\exp [-{t}_{\mathrm{onset}}/\tau ]$ appearing in Equation (10) would determine the fraction of star/disk systems that could produce Jovian planets. In the absence of additional mechanisms that suppress planet formation, the expectation value $\langle \exp [-{t}_{\mathrm{onset}}/\tau ]\rangle$ is proportional to the planetary occurrence fraction. In practice, however, the timescale t_onset depends on stellar mass, so the corresponding planetary occurrence fraction is also mass dependent. In addition, the mass accretion rate depends on stellar mass, so the resulting masses of the planets that form also depend on M_* (see Section 4). Note that the observed occurrence rate for Jovian planets (M_p = 1–10M_J), while highly uncertain, is estimated to fall in the range of 10–30% for FGK stars (integrated over all separations; see Winn & Fabrycky 2015; M. R. Meyer et al. 2020, in preparation, and references therein). Consistency implies that the typical delay time for the onset of rapid gas accretion is bounded by t_onset ≲ τ ∼ 5 Myr. This estimate is compatible with estimates from core accretion models (Pollack et al. 1996; Benz et al. 2014), although giant planet formation is likely to face additional bottlenecks, and this issue should be explored further.

3.2. Mass Accretion Rate onto the Planet

The mass accretion rate onto a forming planet represents a fundamental variable of the problem. We can conceptually divide the process into three components. (a) On the largest scale of the circumstellar disk, an accretion flow moves mass inward at a rate ${\dot{M}}_{{\rm{d}}}$. Some fraction of this mass flow will enter the sphere of influence of the planet, where this fraction can vary from system to system. (b) Since the sphere of influence of the planet grows with planetary mass, the mass accretion rate onto the planet is expected to increase as the planet becomes more massive. ⁵ (c) Only some fraction of the material that enters the Hill sphere will actually be accreted onto the planet. This treatment includes these factors as described below.

For the latter stages of gas accretion, the rate of accretion onto the planet depends on the Hill radius R_H raised to some power, where

$$\begin{eqnarray}&&{R}_{{\rm{H}}}=a{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{3{M}_{* }}\right)}^{1/3},\end{eqnarray} \tag{ 12 }$$

where a is the semimajor axis of the forming planet and M_p is its instantaneous mass. The fiducial cross section scales as R_H ² (Zhu et al. 2011; Dürmann & Kley 2017). However, numerical simulations indicate that the mass accretion rate onto forming planets scales according to a steeper power law $\dot{M}\sim {R}_{{\rm{H}}}^{4}$. The additional factor of R_H ² arises from the shock of the inward flowing gas at the boundary defined by the Hill radius (Tanigawa & Watanabe 2002; Tanigawa & Tanaka 2016). The shock raises the density by a factor of ${\left(v/{v}_{s}\right)}^{2}$, where the incoming speed v is given by v ∼ ΩR_H and v_s is the sound speed (see also Lee 2019). As a result, the mass accretion rate must scale with the mass of the forming planet

$$\begin{eqnarray}&&\dot{M}\sim {R}_{{\rm{H}}}^{4}\sim {M}_{p}^{4/3}.\end{eqnarray} \tag{ 13 }$$

If additional processes are operational, then accretion onto the planet can be suppressed further. Possible mechanisms that come into play include suppression by gap opening (which limits the density of the incoming flow; Malik et al. 2015; Ginzburg & Chiang 2019), strong magnetic fields on the planetary surface (Batygin 2018; Cridland 2018), suppression of direct infall due to the initial angular momentum of the incoming material (Machida et al. 2008), and processes taking place within circumplanetary disks (Szulágyi 2017; Fung et al. 2019). The latter three processes are analogous to well-studied mechanisms involved in the star formation process. Strong stellar magnetic fields truncate the disk and shut down direct accretion, which occurs along field lines at an attenuated rate. The initial angular momentum of the pre-collapse gas causes most of the incoming material to fall initially onto a circumstellar disk, rather than directly onto the star. Although these processes are not well studied in the context of planet formation, they will act to suppress mass accretion and are likely to vary from source to source, thereby introducing a random variable into the problem.

Putting all of the components together, we can write the mass accretion rate onto the planet in the form

$$\begin{eqnarray}&&\dot{M}=\xi {\dot{M}}_{0}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)}^{4/3}.\end{eqnarray} \tag{ 14 }$$

In this form, the benchmark value of the mass accretion rate is that expected when the planet has a Jovian mass. In order of magnitude, for typical disk properties and accretion efficiencies, we expect ${\dot{M}}_{0}\sim 1{M}_{{\rm{J}}}$ Myr⁻¹, but a range of values are possible (Helled et al. 2014). The leading coefficient ξ is a random variable on [0,1] that takes into account both the efficiency for material to enter the Hill sphere and the efficiency with which the planet accretes the incoming material. In principle, one could keep these two efficiencies as separate variables. In this treatment, however, we consider only a single random variable ξ. Moreover, given that the mass accretion process has a well-defined scale set by the mass accretion rate through the disk, the natural starting point is to consider the random variable ξ to have a uniform distribution. ⁶

The mass accretion rate onto the planet from Equation (14) is a rapidly increasing function of planetary mass. Once the planet grows sufficiently large, the rate of accretion onto the planet can become comparable to the rate of mass accretion through the disk. At this point in time, $\dot{M}$ must saturate, and would no longer follow Equation (14). For typical parameters (ξ = 1/2 and ${\dot{M}}_{0}=1{M}_{{\rm{J}}}$ Myr⁻¹; e.g., Helled et al. 2014), and a disk accretion rate ${\dot{M}}_{{\rm{d}}}={10}^{-8}{M}_{\odot }$ yr⁻¹ (e.g., Gullbring et al. 1998; Hartmann et al. 2016), this crossover point corresponds to planetary mass M_p ≈ 10M_J. Larger disk accretion rates require larger planetary masses to reach saturation (which would depend on stellar mass). Since these mass scales are outside the range of interest, and since including this effect requires specification of additional parameters, we do not model the saturation of $\dot{M}$ in this paper. Nonetheless, if the saturation of the mass accretion rate takes a specified form, it is straightforward to incorporate this complexity. Finally, we note that including this effect leads to lower mass accretion rates for sufficiently large planets, thereby leading to somewhat lower masses.

3.3. Probability Distribution Function

Starting with the form for the mass accretion rate from Equation (14), we can integrate over time to find an expression for the planet mass

$$\begin{eqnarray}&&{M}_{{\rm{p}}}={M}_{\mathrm{II}}{\left[1-\xi \displaystyle \frac{{M}_{\mathrm{II}}^{1/3}{\dot{M}}_{0}\tilde{t}}{3{M}_{{\rm{J}}}^{4/3}}\right]}^{-3}.\end{eqnarray} \tag{ 15 }$$

The scale M_II is the mass of the planet at the end of phase II (start of phase III), when gas accretion takes place rapidly. For the sake of definiteness, we take M_II = 20M_⊕ (about twice the core mass). Expression (15) thus contains two random variables $(\xi ,\tilde{t})$. The cumulative probability that one will draw the two variables in order to obtain a planet mass less than M_p is given by the integral

$$\begin{eqnarray}&&P(m\lt {M}_{{\rm{p}}})={\int }_{0}^{{t}_{* }}\displaystyle \frac{d\tilde{t}}{\tau }\exp [-\tilde{t}/\tau ]+{\int }_{{t}_{* }}^{\infty }\displaystyle \frac{d\tilde{t}}{\tau }\exp [-\tilde{t}/\tau ]\displaystyle \frac{{t}_{* }}{\tilde{t}},\end{eqnarray} \tag{ 16 }$$

where we have defined

$$\begin{eqnarray}&&{t}_{* }\equiv \displaystyle \frac{3{M}_{{\rm{J}}}^{4/3}}{{M}_{\mathrm{II}}^{1/3}{\dot{M}}_{0}}\left[1-{\left(\displaystyle \frac{{M}_{\mathrm{II}}}{{M}_{{\rm{p}}}}\right)}^{1/3}\right].\end{eqnarray} \tag{ 17 }$$

The timescale t_*(M_p) is the time required to build a planet of mass M_p at the maximum rate (with ξ = 1). The differential probability (probability density function) is thus given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dP}}{{{dM}}_{{\rm{p}}}} & = & \displaystyle \frac{1}{\tau {\dot{M}}_{0}}{\left(\displaystyle \frac{{M}_{{\rm{J}}}}{{M}_{{\rm{p}}}}\right)}^{4/3}{\displaystyle \int }_{{t}_{* }}^{\infty }\displaystyle \frac{{dt}}{t}\exp [-t/\tau ]\\ & = & \displaystyle \frac{1}{\tau {\dot{M}}_{0}}{\left(\displaystyle \frac{{M}_{{\rm{J}}}}{{M}_{{\rm{p}}}}\right)}^{4/3}{E}_{1}({t}_{* }/\tau ),\end{array}\end{eqnarray} \tag{ 18 }$$

where E₁(z) is the exponential integral (Abramowitz & Stegun 1972) and where t_* is given by Equation (17). Note that only the product ${\dot{M}}_{0}\tau$ of the fiducial mass accretion rate and timescale τ of the disk lifetime distribution appear in the final expression. To evaluate the above expressions, it is thus useful to define a composite parameter—which is a mass scale—according to

$$\begin{eqnarray}&&{M}_{0}\equiv {\left(\displaystyle \frac{{M}_{\mathrm{II}}}{{M}_{{\rm{J}}}}\right)}^{1/3}{\dot{M}}_{0}\tau .\end{eqnarray} \tag{ 19 }$$

For typical parameter values (M_II = 20M_⊕, ${\dot{M}}_{0}\sim 2{M}_{{\rm{J}}}$ Myr⁻¹, and τ = 5 Myr), we find the scale M₀ ≈ 4M_J. Given that the parameters are uncertain, we can find the distributions of planet masses with varying values of M₀. Note that the values of the starting mass M_II and the timescale τ of the disk lifetime distribution are relatively well known. As a result, specification of the mass scale M₀ is largely determined by the baseline value of the mass accretion rate ${\dot{M}}_{0}$.

Given the above results, we can construct PMFs in two ways: first, we can sample the variables in Equation (15) using their specified distributions, calculate the masses, and construct the histogram of results. Second, we can use Equation (18). Both results are shown in Figure 1 for different choices of the fiducial mass scale in the range of M₀ = 1—10 M_J. The sampling results are shown as the colored curves and the analytic results are shown as black-dashed curves (note that the results are essentially the same). For comparison, the solid black curve shows the expected power-law mass function of the form dF/dM ∝ M^−1.3. With low accretion rates (mass scale) M₀ = 1M_J (green curve), the model produces a PMF that is steeper than the target distribution, with a deficit of high-mass planets. For larger values of M₀, the theoretical PMF approaches a power-law form over most of the relevant mass range and is in reasonable agreement with the target distribution. Nonetheless, the theoretical distributions remain slightly steeper.

Zoom In Zoom Out Reset image size

As mentioned above, the mass accretion rate from Equation (14) is an increasing function of time. If the planet mass become sufficiently large, the predicted mass accretion rate becomes larger than the rate supplied by the circumstellar disk, so that the planetary mass accretion rate must saturate and approach a constant (maximum) value. The PMF here is constructed using the increasing form of the mass accretion rate. This assumption leads to somewhat larger planets for the tail of the distribution. For completeness, we also consider the opposite limit in the Appendix, where the mass accretion rate saturates and thus becomes constant while most of the planetary mass is acquired. This Appendix also illustrates how the PMF framework constructed in the paper can be modified or generalized.

The treatment thus far does not include the mass of the host star in the analysis. This section uses the observed properties of circumstellar disks, which depend on the mass of the star, to inform the models of the PMF constructed in the previous section. Within this framework, the PMF is specified up to the characteristic mass scale M₀ defined in Equation (19). Here we use the empirical scaling laws outlined in Section 2 to define how the mass scale M₀ varies as a function of the stellar mass, and use the result to specify the PMF for varying M_*.

Consider the accretion of gas onto a forming giant planet where the rate is a constant fraction of the disk accretion rate. Including the time dependence of the disk mass accretion rate from Equation (7), the total mass accreted—and hence the mass of the planet—must be given by the integral

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{p}}} & = & {\int }_{{t}_{\mathrm{onset}}}^{{t}_{\mathrm{end}}}\epsilon {\dot{M}}_{d(0)}\displaystyle \frac{{m}^{2}{dt}}{{\left(1+t/{t}_{0}\right)}^{3/2}}\\ & = & 2\epsilon {\dot{M}}_{d(0)}{t}_{0}{m}^{2}\left[\displaystyle \frac{1}{{\left(1+{t}_{\mathrm{onset}}/{t}_{0}\right)}^{1/2}}-\displaystyle \frac{1}{{\left(1+{t}_{\mathrm{end}}/{t}_{0}\right)}^{1/2}}\right],\end{array}\end{eqnarray} \tag{ 20 }$$

where m is the dimensionless stellar mass from Equation (6). In order to define the characteristic mass scale M₀ as a function of m, we use the scaling laws from Section 2. The core formation time varies with stellar mass according to Equation (8). The typical time for accretion to end is given by the disk lifetime parameter τ, which varies with m according to Equation (5). We can then define the mass scale M₀ using the ansatz

$$\begin{eqnarray}\begin{array}{rcl}{M}_{0} & = & \left[\epsilon {\dot{M}}_{d(0)}{t}_{0}\right]{ \mathcal F }(m),\qquad \mathrm{where}\\ { \mathcal F }(m) & \equiv & 2{m}^{2}\left[\frac{1}{{\left(1+\sqrt{m}\right)}^{1/2}}-\frac{1}{{\left(1+5/m\right)}^{1/2}}\right].\end{array}\end{eqnarray} \tag{ 21 }$$

Notice that the function ${ \mathcal F }(m)\to 0$ for stellar mass m = ${\left({\tau }_{1}/{t}_{{\rm{c}}0}\right)}^{2/3}\sim 3$. Our interpretation of this finding is that relatively few planets should form in disks associated with larger stellar masses (see the discussion below).

The relative mass fraction ${ \mathcal F }(m)$ from Equation (21) is shown in Figure 2 for standard values of the input parameters (solid red curve) and corresponding cases where the parameters are allowed to vary (lighter cyan curves). All of the curves follow the same basic trend: At low stellar masses, the accretion rate through the disk ($\dot{M}\sim {m}^{2}$) becomes small, and this trend leads to small relative mass scales and hence the partial suppression of giant planet formation. In the opposite limit of large stellar mass, m ∼ 3, the relative mass becomes small due to the lack of time for accretion. The time required for core formation becomes larger, while the disk lifetimes become shorter, so that little time is left for the cores (planets) to accrete gas. The observed scaling laws thus imply that larger stars (m ≳ 3) have difficulty producing giant planets via the core accretion paradigm. In the intermediate regime, for m ∼0.5–2, the relative mass fraction varies (almost) linearly with increasing mass. This result suggests that distribution of companion mass ratios should be nearly the same for stars in this mass range. In contrast, the distribution should be steeper for stellar masses outside this range (m ≪ 0.5 and/or m ≫ 2).

Zoom In Zoom Out Reset image size

For completeness, we can also evaluate the effects of stellar mass dependence in the limit where the disk mass accretion rate does not vary appreciably during the epoch when the planet accretes most of its mass. This limit corresponds to large values of the timescale t₀. One can show that in the limit of large t₀, the quantities from Equation (21) become

$$\begin{eqnarray}&&{M}_{0}\equiv \epsilon {\dot{M}}_{d(0)}{\tau }_{1}\qquad \mathrm{and}\qquad { \mathcal F }(m)\equiv m\left[1-\displaystyle \frac{{t}_{{\rm{c}}0}}{{\tau }_{1}}{m}^{3/2}\right].\end{eqnarray} \tag{ 22 }$$

In this case, the function ${ \mathcal F }(m)\to 0$ for stellar mass $m\,={\left({\tau }_{1}/{t}_{{\rm{c}}0}\right)}^{2/3}$, which is the same as found above. ⁷

We can now construct the PMF for different values of the stellar mass. The stellar mass dependence is incorporated by assuming that the total available mass for making a giant planet is proportional to the relative mass fraction ${ \mathcal F }(m)$ defined above. For a given stellar mass, the characteristic mass scale M₀ that determines the shape of the PMF is then taken to have the form

$$\begin{eqnarray}&&{M}_{0}={ \mathcal F }(m)\times 5{M}_{{\rm{J}}}.\end{eqnarray} \tag{ 23 }$$

Note that we are assuming that the mass scale for forming giant planets is proportional to the total mass available. However, since the disk can in principle form multiple planets, the coefficient in Equation (23) is chosen to be only a fraction of the disk mass. If the mass scale M₀ is used in the PMF constructed according to Section 3, we obtain the distributions shown in Figure 3. The figure shows the mass functions for planets forming around stars with masses m = 0.25–2.5. The planetary mass distributions for host stars in the range of m = 1–2.5 display a similar shape. In contrast, stars with lower masses show a significant deficit of larger planets, as shown by the PMF for m = 0.5 (yellow curve) and especially that for m = 0.25 (red curve). These steeper distributions imply that red dwarfs are less likely to have large planets, but are still expected have some Jovian companions (see also the more detailed planet formation models of Laughlin et al. 2004).

Zoom In Zoom Out Reset image size

For even larger host stars, the available mass supply decreases due to the longer core formation times and shorter disk lifetimes (see Figure 2). Since ${ \mathcal F }(m)\to 0$ for m ≳ 3, this present framework does not provide a mass distribution for planets around higher mass stars: such planets should be rare. This implication is roughly consistent with current observational results (Reffert et al. 2015), which find that the planet occurrence rate drops significantly for stellar masses larger than M_* ∼ 2.5–3M_⊙.

Emerging observational trends (M. R. Meyer et al. 2020, in preparation) comparing low-mass brown dwarf companions to gas giant planets, as a function of stellar mass, suggest that distributions of companion mass ratios are more universal than the distributions of companion masses themselves. In addition to the PMF, we thus consider the distribution of the mass ratios q = M_p/M_* in Figure 4 (for the same host star masses as before; compare with Figure 3). The resulting distributions of mass ratio are roughly similar across the range of stellar masses for m ≥ 0.5. Once again, however, the distribution for stars with the smallest mass, m = 0.25 (red curve), is steeper and shows a deficit of large-q planets compared to the others. This deficit is due to the small fiducial mass scale, which results from the smaller available mass supply in the disks surrounding low-mass stars. Although the curves are shown down to ratios q = 3 × 10⁻⁴ for all stellar masses, for m = 0.25 this value corresponds to a planet with only M_p ∼ 3M_⊕. Planets with such small masses can accrete little gas, so that their production mechanism lies outside this framework for the formation of gas giants. Finally, we note that the similarity between the PMFs and the mass-ratio distributions (Figures 3 and 4) arises because the distributions have a nearly power-law form. In the limit where the PMF is exactly a power law, the two distributions would have the same shape.

Zoom In Zoom Out Reset image size

This finding has observational implications. Microlensing surveys measure companion mass ratios and favor low-mass stars (which are more common than larger stars). On the other hand, both RV and transit surveys have historically used our solar system as a template and favored solar-type stars. As a result, the observed distribution of companion mass ratios is predicted to be somewhat steeper for microlensing surveys than for the others. For the largest stellar mass considered here, m = 2.5 (magenta curve), the PMF is nearly the same as for solar-type stars, but even higher mass stars (m > 3) should have relatively fewer large planets (because the relative mass fraction ${ \mathcal F }\to 0$). Future observational work should determine how the distributions of planetary masses and mass ratios vary across the entire spectrum of stellar masses.

This paper constructs a series of models for the initial mass function for gas giant planets forming through the core accretion paradigm, including its dependence on the mass of the host star. This section provides a brief summary of our results (Section 5.1), followed by a discussion of their implications (Section 5.2).

5.1. Summary of Results

5.2. Discussion

We would like to thank Veenu Seri for early work that helped formulate the problem. We also thank Konstantin Batygin, Nuria Calvet, and Greg Laughlin for many useful discussions, and thank the referees for their comments on the manuscript. This work was supported by the University of Michigan, the NASA JWST NIRCam project (contract number NAS5-02105), and the NASA Exoplanets Research Program (grant number NNX16AB47G).

In this Appendix, we consider the limiting case where the mass accretion rate onto the planet approaches a constant and calculate the corresponding PMF. In this context, the planetary mass M_p is given by

$$\begin{eqnarray}&&{M}_{{\rm{p}}}={M}_{\mathrm{II}}+{\dot{M}}_{0}\xi t={M}_{\mathrm{II}}+M,\end{eqnarray} \tag{ A1 }$$

where t is the disk lifetime, ξ is a uniform-random variable as before, and the mass accretion rate ${\dot{M}}_{0}$ is constant. The second equality defines the mass M accreted after the onset of rapid accretion.

The distribution of the mass M is determined by the cumulative probability

$$\begin{eqnarray}&&P(m\lt M)={\int }_{0}^{{t}_{* }}\displaystyle \frac{{dt}}{\tau }\exp [-t/\tau ]+{\int }_{0}^{{t}_{* }}\displaystyle \frac{{dt}}{\tau }\exp [-t/\tau ]\displaystyle \frac{{t}_{* }}{t},\end{eqnarray} \tag{ A2 }$$

where the timescale t_* is defined by

$$\begin{eqnarray}&&{t}_{* }=M/{\dot{M}}_{0}.\end{eqnarray} \tag{ A3 }$$

The probability distribution function for the accreted mass M thus takes the form

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dM}}=\displaystyle \frac{1}{{\dot{M}}_{0}\tau }{E}_{1}\left(\displaystyle \frac{M}{{\dot{M}}_{0}\tau }\right)=\displaystyle \frac{1}{{M}_{0}}{E}_{1}\left(\displaystyle \frac{M}{{M}_{0}}\right),\end{eqnarray} \tag{ A4 }$$

where the second equality defines ${M}_{0}={\dot{M}}_{0}\tau$. If the starting mass M_II is constant, then the corresponding probability distribution for the planetary mass M_p has the form

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{{dM}}_{{\rm{p}}}}=\displaystyle \frac{1}{{M}_{0}}{E}_{1}\left(\displaystyle \frac{{M}_{{\rm{p}}}-{M}_{\mathrm{II}}}{{M}_{0}}\right).\end{eqnarray} \tag{ A5 }$$