The Mass Budgets and Spatial Scales of Exoplanet Systems and Protoplanetary Disks

We use the CORALIE–HARPS survey of Mayor et al. (2011) to calculate the solid mass distribution of radial-velocity-detected exoplanets. This survey was conducted using a volume-limited sample of stars similar in mass to the Sun. It is sensitive to planets as low in mass as ∼3 M_⊕ short orbital periods, while the largest planets are detected out to a distance of ∼10 au (Figure 1(b)). The detected planet mass distribution is shown with the purple line in Figure 3, with a median detected planet mass of ∼100 M_⊕.

Zoom In Zoom Out Reset image size

To highlight the effects of the stellar sample selection, the planet mass distribution of all detected exoplanets from the online archive exoplanet.eu that was used by Manara et al. (2018) is shown in pink in Figure 3. This sample is a heterogeneous combination of exoplanets surveyed with different sensitivities in terms of planet mass. Hot Jupiters are overrepresented in this sample, which leads to a median planet mass of 300 M_⊕ that is a factor of 3 higher than the volume-limited stellar sample of the Mayor et al. (2011) survey, as indicated by the purple arrow in Figure 3 and the horizontal lines in Figure 1.

Even within a single survey, the sensitivity to detect planets of different masses varies from star to star. The survey detection efficiency as a function of planet mass and orbital period was characterized by Mayor et al. (2011) using injection-recovery tests. For each detected exoplanet, j, we extract the survey completeness, C_j , as described in Fernandes et al. (2019). This completeness is defined as the probability that the planet is detected in the survey based on its mass and orbital period. We calculate a weight factor, ${w}_{j}=\tfrac{1}{{N}_{\star }\,{C}_{j}}$, where N_⋆ is the number of stars in the survey, to debias the mass distribution of exoplanets (Figure 1(c) and Figure 3, teal line). The median planet mass is ∼5 M_⊕, which is more than an order of magnitude lower than the median mass of the biased distribution of detections, as indicated by the teal arrow. Accounting for the combined detection and selection effects yields a median mass that is almost two orders of magnitude lower than the median mass of all known exoplanets from the online archive exoplanet.eu.

The more massive planets are gas giants, and thus a significant fraction of their measured mass is not in solids. To calculate the solid components in these exoplanets, M_c , we use the best-fit scaling relation between heavy-element abundance and planet mass from Thorngren et al. (2016), given by

$$\begin{eqnarray}{M}_{c}=\left\{\begin{array}{cc}60\,{M}_{\oplus }{\left(\displaystyle \frac{M}{{M}_{J}}\right)}^{0.6} & \mathrm{if}\,M\gt {M}_{c}\\ M & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 1 }$$

The correction for planet heavy-element mass primarily affects the more massive end of the planet mass distribution (Figure 3, teal dotted line). We ignore the intrinsic dispersion in the relation because we are primarily interested in the core mass distribution, which is less affected by scatter than estimates for individual planets.

A final correction that has to be made before exoplanets can be compared to disks is to calculate what fraction of stars have planets of a given mass. Some stars may not have planetary systems at all, or some planets may be below the survey detection limit. Assuming that the mass distribution of detected exoplanets applies to all stars can significantly bias the results, especially if planet occurrence rates are much smaller than 1. For example, giant planets are typically detected around ∼10% of stars (e.g., Cumming et al. 2008; Fernandes et al. 2019), so as long as 10% of stars have similarly massive disks, the average disk mass can be much lower without creating a conundrum.

The fraction of stars with planetary systems, here denoted by F, is related to the planet occurrence rate, f = ∑_j w_j , by

$$\begin{eqnarray}&&F=\displaystyle \frac{f}{\bar{n}},\end{eqnarray} \tag{ 2 }$$

where $\bar{n}$ is the average number of planets per planetary system. The number of planets per system has not been robustly determined from radial velocity surveys. The average detected number of planets per system is 1.5 in the sample of Mayor et al. (2011), but this number increases for the smaller planets that are intrinsically more common, so this number is likely a lower limit. Therefore, we assume a somewhat higher number of two planets per system (${\bar{n}}_{\mathrm{RV}}=2$). We note that the detection efficiency of multiplanet systems was not considered in the injection-recovery tests of Mayor et al. (2011), and thus it is difficult to determine the intrinsic planet multiplicity based on publicly available data.

The mass in solids per planetary system, M_s , is estimated from each individual planet detection as ${M}_{s}={\bar{n}}_{\mathrm{RV}}\,{M}_{c}$. The fraction of stars with at least a given amount of solids in its planetary system is then

$$\begin{eqnarray}&&F(\gt {M}_{s})=\displaystyle \frac{1}{{\bar{n}}_{\mathrm{RV}}}\sum _{j,{M}_{{sj}}\gt {M}_{s}}{M}_{{sj}}{w}_{j}.\end{eqnarray} \tag{ 3 }$$

After applying these corrections, the fraction of stars with planetary systems more massive than 10M is 32%, e.g., the radial velocity planet population represents roughly a third of solar-mass stars.

We derive the mass distribution of transiting exoplanets based on the survey completeness from Mulders et al. (2018). This completeness was calculated on a grid of orbital period and planet radius using all main-sequence stars in the Kepler survey. This stellar sample has a median mass, estimated from isochrone fitting using isoclassify (Huber et al. 2017), of 1.0 solar mass and a standard deviation of 0.3 M_⊙, and its occurrence rate is representative of that of solar-mass stars. We therefore do not make any additional mass or effective temperature cuts, as they would only reduce the sample size. The calculation uses the detected planet candidates from the Kepler DR25 planet catalog (Thompson et al. 2018), the survey detection efficiency measurements from KeplerPORTS (Burke & Catanzarite 2017), and the Gaia stellar parameters (Berger et al. 2018). For each detected planet candidate, we determine the weight factor ${w}_{j}=\tfrac{1}{{N}_{\star }\,{C}_{j}}$ by evaluating the survey completeness at the period and radius of the planet.

Next, we calculate the planet mass from the detected planet radii using the best-fit mass–radius relation from Chen & Kipping (2017):

$$\begin{eqnarray}M({M}_{\oplus })=\left\{\begin{array}{cc}{\left(\displaystyle \frac{R}{{R}_{\oplus }}\right)}^{1/0.28} & \mathrm{if}\,R\lt {R}_{t}\\ {M}_{t}{\left(\displaystyle \frac{R}{{R}_{t}}\right)}^{1/0.59} & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 4 }$$

with M_t = 2M and ${R}_{t}={({M}_{t}/{M}_{\oplus })}^{0.28}{R}_{\oplus }$. While the mass–radius relation in Chen & Kipping (2017) is probabilistic, we ignore this aspect since we are considering the wider distribution of planet properties, which is less affected by intrinsic scatter. The mass in solids per planet is again calculated using the estimates of heavy-element abundance from Thorngren et al. (2016) in Equation (1). The planet solid mass distribution is shown in Figure 4 with the solid dark-purple line. We note that the cumulative occurrence of planets adds up to a number larger than 1, showing that a correction for multiple planets per star is necessary.

Zoom In Zoom Out Reset image size

As a final step, we calculate the fraction of stars with the given amount of planetary mass in solids by correcting for the number of planets per system, n_j . The fraction of stars with planets, F, is related to the individual planet weights as

$$\begin{eqnarray}&&F=\sum _{j}\displaystyle \frac{{w}_{j}}{{n}_{j}},\end{eqnarray} \tag{ 5 }$$

where n_j is the number of planets in the system. Unfortunately, in transit surveys n_j cannot be directly measured because of geometrical biases: the probability that multiple planets transit their star is low even for low mutual inclinations of the planets' orbits. Thus, the number of detected transiting planets around a star is not a good measure of the real multiplicity because most planets in a system will not be transiting.

However, the average number of planets per system, $\bar{n}$, can be estimated using a careful treatment of the complex geometric detection biases (e.g., Lissauer et al. 2011; Tremaine & Dong 2012). The mass in solids per planetary systems is approximated as ${M}_{s}=\bar{n}\,{M}_{c}$. The previous equation can then be simplified as

$$\begin{eqnarray}&&F=\displaystyle \frac{1}{\bar{n}}\sum _{j}{w}_{j},\end{eqnarray} \tag{ 6 }$$

and the fraction of stars with at least a given amount of solids is

$$\begin{eqnarray}&&F(\gt {M}_{s})=\displaystyle \frac{1}{\bar{n}}\sum _{{M}_{{sj}}\gt {M}_{s}}{w}_{j}.\end{eqnarray} \tag{ 7 }$$

Figure 4 shows the fraction of stars with at least a given amount of solids for different values of $\bar{n}=1,2,3,6$. Forward modeling of Kepler planet populations indicates that $\bar{n}$ is typically in the range 2–6 (Mulders et al. 2018; Zhu et al. 2018; He et al. 2019), though the intrinsic distribution of planet multiplicity is uncertain and model dependent (Zhu et al. 2018; Sandford et al. 2019; Zink et al. 2019).

To identify the best value of $\bar{n}$ for the Kepler sample, we compare the different distributions to the fraction of stars with planetary systems estimated using a parametric forward model from Mulders et al. (2018) (dashed line). We find that the assumption $\bar{n}=3$ works well for the Kepler sample. This yields a cumulative fraction of Sun-like stars with detectable planets of F = 47% (Figure 6), which is also consistent with estimates from other forward models (e.g., Zhu et al. 2018; He et al. 2019; Yang et al. 2020).

Solid masses of Class II protoplanetary disks have been measured with ALMA in multiple star-forming regions (e.g., Ansdell et al. 2016; Barenfeld et al. 2016; Pascucci et al. 2016; Ansdell et al. 2017; Cieza et al. 2019). These mass measurements represent the amount of solids in millimeter-sized dust grains and are calculated from the millimeter flux density using the stellar distance and assumptions on the disk temperature, dust opacity, and optical depth. Here we use the combined Chamaeleon I (Pascucci et al. 2016) and Lupus (Ansdell et al. 2016) stellar samples because for these regions stellar masses have also been homogeneously recomputed (e.g., Pascucci et al. 2016; Alcalá et al. 2017; Manara et al. 2017). We use the disk masses from the online table of Mulders et al. (2017), calculated using an average disk temperature of T = 20 K. The cumulative distribution of disk masses is shown with the dotted line in Figure 5. From this sample, we select only disks with approximately solar-mass host stars (0.5 < M_⋆/M_⊙ < 2). This cut mainly removes low-mass stars because they are more numerous than Sun-like stars. Because of the positive correlation between protoplanetary disk mass and stellar mass, this cut increases the median disk mass by a factor of five from ∼2M to ∼10 M_⊕.

Zoom In Zoom Out Reset image size

Stars without protoplanetary disks (Class III objects) are excluded from these surveys, and a correction has to be made to estimate the fraction of stars with a given amount of solids. Here, we assume that Class III disks have a negligible amount of detectable dust compared to Class II objects, as is the case for Lupus at least (Lovell et al. 2020). The fraction of stars with protoplanetary disks decreases with time (e.g., Haisch et al. 2001) and is typically 50% at the age of Lupus and Chamaeleon I (Luhman et al. 2005). Therefore, we normalize the mass distribution of disks such that the cumulative number of stars with disks of any mass is F_disk = 50% (Figure 6). Because the disk fraction tends to decrease with time, the cumulative distribution function (CDF) of an older (younger) star-forming region would have a lower (higher) normalization.

Zoom In Zoom Out Reset image size

Figure 6 shows the distribution of solid masses derived from samples of transiting planets, radial velocity planets, and Class II protoplanetary disks.

The Kepler and radial velocity surveys overlap in the range 5 < M_s /M_⊕ < 20, where both surveys trace the same population of exoplanets in terms of mass and period (see also Figure 7). At larger planetary system solid masses, the radial velocity planets trace an additional population of colder giants outside of 1 au that is not detected with Kepler, and the distributions deviate.

Zoom In Zoom Out Reset image size

All three distributions peak near ≈10 M_⊕, consistent with previous estimates of the average amount of solids in planets and disks (e.g., Najita & Kenyon 2014; Mulders et al. 2015b; Pascucci et al. 2016). The distribution of disk masses appears wider than that of exoplanets. However, we do not see many exoplanetary systems whose masses are higher than the most massive disks as claimed by Manara et al. (2018), which we attribute to the authors not taking into account necessary bias corrections for exoplanets. While a number of planetary systems exist with solid masses above 100 M_⊕—a region where few disks are detected—these constitute a tiny fraction of the total systems in the tail of the probability distribution (see Figure 6) and do not constitute a significant mass budget issue at the population level.

We note that the fraction of stars with disks and planets never goes above 50%, so our analysis is only accounting for half the Sun-like stars. The other stars may either have exoplanets outside of the detection limits of transit and radial velocity surveys or never have formed planets at all. The observed fraction of stars with disks depends on the evolutionary phase and typically decreases with time (e.g., Haisch et al. 2001). It is therefore unlikely that protoplanetary disks in older regions like Upper Sco, with a disk fraction below 25% (Luhman & Mamajek 2012), represent the formation sites of all observed exoplanets. While planet formation is likely ongoing in the disks that survive to this age, more than half the exoplanet systems must have formed around stars that have already dispersed their disks by the age of Upper Sco, 5–10 Myr.

Younger regions, such as Taurus, with a higher disk fraction of 75% (Luhman et al. 2010), therefore include a significant fraction of disks that do not form the exoplanets that can be detected with radial velocity or transit surveys. Either some of these disks form different types of planets (such as traced by microlensing or debris), or some of them may not form exoplanets at all.

In general, there is relatively good agreement between the different solid mass distributions. This result supports the earlier conclusion from Najita & Kenyon (2014) that the mass reservoirs of solids contained in the Kepler and radial velocity exoplanets are of similar magnitude to those in protoplanetary disks.

For exoplanets, the planetary system solid mass M_s was estimated from each individual planet detection by accounting for the average number of planets per planetary system, $\bar{n}$. We continue following this approach here and estimate the spatial scale where this mass is detected from the semimajor axis of that planet, which in turn is calculated from the planet orbital period using Kepler's third law and assuming a solar-mass star.

These spatial scales are shown in Figure 7, with each symbol representing a single detected planet. The symbol area size for each planet is proportional to the weight factor, w_j , to illustrate the detection biases of these surveys. That is, large symbols correspond to planets detected despite a low survey detection efficiency, and thus those symbols represent a large number of undetected planets. Likewise, small symbols correspond to those planets that are easy to detect, and such planets are likely not underrepresented in our data.

The majority of planetary systems, those with solid masses in the range 5–20 M_⊕, are detected between 0.05 and 0.5 au. Their spatial distribution follows the known orbital period dependence of sub-Neptunes from the Kepler survey, where planet occurrence is measured to be roughly constant in the logarithm of orbital period between 10 days and 1 yr, or 0.1–1 au for a solar-mass star (e.g., Petigura et al. 2013). Inside of 0.1 au, the Kepler planet occurrence rate falls off rapidly with decreasing orbital period (e.g., Youdin 2011; Howard et al. 2012). Planetary systems with solid masses above 20 M_⊕ are primarily detected in the radial velocity survey between 1 and 10 au, consistent with estimates of Bryan et al. (2016) and Fernandes et al. (2019). The onset of this rise in giant planet occurrence with semimajor axis (e.g., Dong & Zhu 2013; Santerne et al. 2016) is also seen here in the Kepler data between 0.5 and 1 au, near the edge of the detection limit at ∼1 au.

Planetary systems with solid masses less than 5 M_⊕ are primarily detected within 0.2 au. This is, at least in part, due to the detection bias of Kepler, whose minimum detectable planet size increases with orbital period and is not sensitive to Earth-sized planets outside of ∼50 days. However, occurrence rate studies have also shown that smaller planets tend to be located closer in to the star than larger ones (Howard et al. 2012; Petigura et al. 2018), suggesting that the apparent correlation between spatial scale and planet size in Figure 7 may be a real feature of the exoplanet population and not purely be attributed to survey biases and detection limits.

Overall, the spatial scales of detected exoplanets increase with planetary system solid mass, a trend also noted by other works (e.g., Martinez et al. 2019). A quantile regression of the median planet semimajor axis at different solid system masses, taking into account the detection efficiency weight factors w, reveals that the Kepler planet population scales as a ∝ M_s ^0.5, while the radial velocity exoplanet population scales as $a\propto {M}_{s}^{1}$ (Figure 9). We note that the scaling derived here is significantly steeper than that of the minimum-mass (extra)solar nebula, which would be M_s ∝ a^0...0.5 for surface density index p = 1.5...2. This shows that this approach is indeed measuring variations in planetary system sizes among stars and not the surface density profiles of individual protoplanetary disks.

Zoom In Zoom Out Reset image size

The detection limit of the combined Kepler and radial velocity exoplanet census roughly follows a steeper diagonal line that scales as $a\propto {M}_{s}^{1.35}$. As mentioned before, the exoplanet census for small (<10 M_⊕) planets is incomplete owing to detection bias outside of 1 au (and outside of 0.2 au for Earth-sized planets). While a decrease in the giant planet distribution with semimajor axis has been detected (Fernandes et al. 2019; Nielsen et al. 2019), there is no clear indication from Kepler planet occurrence rate of a decline with semimajor axis near the detection limit at 1 au, and thus the peak of this planet population could lie at larger semimajor axes. Assuming an occurrence rate constant in the logarithm of orbital period, a similar number of planets would be present between 1 and 10 au as detected with Kepler between 0.1 and 1 au. In that case, the apparent correlation between planet size and semimajor axis in Figures 7 and 9 could be the result of these detection limits. A scenario where the peak of planet occurrence of sub-Neptunes lies near that of the RV giants can neither be excluded nor be confirmed with these data. Microlensing is sensitive to sub-Neptunes at separations between 1 and 10 au (Cassan et al. 2012), but their semimajor axis distribution has not been characterized in the context of a peaked distribution. A census of small planets at larger semimajor axes, such as with the Nancy Grace Roman Space Telescope, is needed to confirm the relation between planet size and semimajor axis.

We calculate the spatial scale of the detected solid mass reservoir based on observed outer disk radii. Note that we do not consider the surface density profiles of individual disks, but that we use disk outer radii as a proxy for the spatial scale at which most solids are located. The spatial distribution of dust in the outer regions of protoplanetary disks can be computed from the intermediate-resolution ALMA data targeting all Class II disks in Lupus and Chameleon I (Ansdell et al. 2016; Pascucci et al. 2016; Ansdell et al. 2018; Long et al. 2018).

Here, we use dust disk outer radii from the analysis of Hendler et al. (2020), who fitted radial profile models to visibility data of a large sample of Class II disks with available ALMA interferometric data and reported the radii that enclose 68% and 90% of the millimeter flux. We take the radii that enclose 68% to be a tracer of where most of the disk mass is located. The flux distribution is a proxy for the solid mass distribution, with the two related by the dust temperature and optical depth. The location that encloses 68% of the flux encloses less than 68% of the mass because temperature increases inward. Therefore, we take the R₆₈ to be a proxy of M₅₀, e.g., the point where half of the mass is detected. A full radiative transfer calculation to estimate the disk surface density is beyond the scope of this paper. However, R₆₈ and R₉₀ are tightly correlated, and their offset is only 0.1 dex (Hendler et al. 2020), and thus the choice for one over the other would not significantly affect the conclusions of this paper.

From the sample of 54 Sun-like stars (0.5 < M/M_⊙ < 2) in Lupus and Chamaeleon I used in Section 2, 36 have a disk radius measurement, of which 6 are an upper limit. We calculate the disk mass from the millimeter flux using the same disk temperature of T = 20 K as in Pascucci et al. (2016). To verify that the sample of disks with radius measurements is not biased toward bigger, brighter disks compared to the entire sample that also includes disks without radius measurements, we calculate the Kolmogorov–Smirnov (KS) distance between the two mass distributions and find a high (p_KS = 83%) probability that the two samples are drawn from the same distribution.

The disk masses and radii are shown in red in Figure 7. Upper limits are indicated with triangles. The disk mass and outer radii are positively correlated (Figure 9), with R₆₈ ∝ M_d ^0.5 and a scatter of 0.2 dex (see also Tripathi et al. 2017; Andrews et al. 2018).

The comparison between spatial scales of detected solids in protoplanetary disks and exoplanets does not significantly change the need for inward migration of pebbles or protoplanets: the spatial scales where solids in protoplanetary disks are detected are two orders of magnitude larger than those of exoplanets. Most of the solids in exoplanet systems are detected between 0.1 and 1 au, while most of the disk emission that represents a similar amount of solids is detected between 10 and 100 au. This is indicated with an arrow in Figure 8(b). If observed exoplanets indeed form from the observed mass reservoir in protoplanetary disks, then planets or their building blocks need to migrate inward over two orders of magnitude during formation. Both planet migration and a pebble accretion scenario could achieve this.

The comparison does show a hint of how this migration mechanism operates: detected exoplanets and disks both show a positive correlation between the amount of solids and the semimajor axis where they are detected (Figure 9). More massive planets are located farther from the star, and the disks that contain enough solids to form them are also larger. This suggests that the observed size scaling of the solid mass reservoir may be retained during the planet formation process, even if exoplanets or their building blocks migrate inward by a factor 100 in semimajor axis.

Within the hypothesis that planets migrate inward through type I or type II migration, it is not clear why such a relation would be preserved, as the final locations of exoplanets are determined more by the location where migrating planets are trapped and less by their starting location. Within the pebble accretion hypothesis, a similar argument can be made that the fast radial drift of pebbles could erase this scaling relation and that the locations of exoplanets are determined more by where pebbles are accreted, and not where they came from. However, since the observed Class II disk dust distributions may already be sculpted by drift of pebbles (e.g., Pinilla et al. 2020; van der Marel & Mulders 2021), this could be seen as an argument in favor of pebble accretion.