California-Kepler Survey. IX. Revisiting the Minimum-mass Extrasolar Nebula with Precise Stellar Parameters

Our study is mostly based on the CKS (Petigura et al. 2017) given its uniformity and precision. We started with the whole CKS catalog, which consists of the 960 bright (Kepler magnitude <14.2) planet hosts, the multiplanet host stars (484), hosts of ultra-short-period planets (USP: P_orb < 1 day; 127), and 109 host stars of various other types of planets. We incorporated the new stellar/planetary properties from Fulton & Petigura (2018) that included Gaia parallax constraints (Gaia Collaboration et al. 2018) and dilution factors from neighboring stars. We imposed several filters to identify a more uniform subsample:

• 1.  
  We restricted our attention to the brightest host stars (Kepler Magnitude <14.2).
• 2.  
  We removed any system that is designated as a false positive by the Kepler team.
• 3.  
  Because planets around evolved stars are harder to detect and usually have larger radius uncertainty, we excluded any evolved star with the ad hoc relation suggested by Fulton et al. (2017):

  $$\begin{eqnarray}&&{R}_{\star }/{R}_{\odot }\gt {10}^{0.00025({T}_{\mathrm{eff}}/{\rm{K}}-5500)+0.20}.\end{eqnarray} \tag{ 1 }$$
• 4.  
  We only included host stars with T_eff between 4700 K and 6500 K to avoid the known limitations of synthetic stellar atmospheric models (Petigura et al. 2017),
• 5.  
  We focused on the sub-Neptune planets ($1{R}_{\oplus }\lt {R}_{p}\,\lt 4{R}_{\oplus }$) in this work. As most of the Kepler planets do not have mass constraints, we have to rely on the measured radii and a mass–radius relationship to infer the planetary masses. Previous works (e.g., Wolfgang & Lopez 2015) showed that planets in the sub-Neptune radius range likely has a rocky core and a H/He envelope of $\lesssim 10$% in mass. In other words, their mass is dominated by the mass of the core. Therefore, the bulk mass of the sub-Neptune planets roughly corresponds to the amount of solid disk materials in them. In contrast, the mass of giant planets is dominated by H/He gas. There is a large uncertainty in the core mass of giant planets; even the core mass of Jupiter is uncertain (Wahl et al. 2017). The cut at $4{R}_{\oplus }$ also excludes the emerging class of "super-puff" planets, where mass–radius relationships fail because the transit radii are likely inflated by high-altitude dusts or aerosols (Gao & Zhang 2020; Libby-Roberts et al. 2020; Wang & Dai 2019).
• 6.  
  We excluded the USP planets because most of them were discovered by an independent study (Sanchis-Ojeda et al. 2014) with a Fourier-based transit detection method. On the other hand, most Kepler planets were identified by the Kepler team with a wavelet-based matched filter (Jenkins et al. 2002). The USPs likely have different detection probability from the other Kepler planets. Moreover, the larger mutual inclination and orbital period spacing of USP planets both suggest orbital migration (Dai et al. 2018), thus making it difficult to incorporate the USPs in the fundamentally in situ MMEN framework.
• 7.  
  Planets with grazing transits (impact parameter b > 0.9) were excluded. A grazing orbit often produces V-shaped transits. V-shaped transits suffer from a higher false-positive probability. In addition, transit parameters tend to be degenerate with each other, leading to a larger uncertainty in the planetary radius.

After the various cuts described above, we are left with a sample size of 1041 planets around 712 stars. Figure 1 shows the distribution of the effective temperature T_eff, metallicity [Fe/H] of the host stars, as well as the measured radii R_p, inferred masses M_p, and orbital periods ${P}_{\mathrm{orb}}$ of the planets.

Zoom In Zoom Out Reset image size

The CKS stellar properties are most precise for stars with $4700\,{\rm{K}}\lt {T}_{{\rm{eff}}}\lt 6500\,{\rm{K}}$, where synthetic spectral models perform best (Petigura et al. 2017). We therefore used the Kepler Object of Interest (KOI) planets (Mathur et al. 2017) for a sample of broader dynamical range in spectral types, M_⋆, and [Fe/H]. We note that the stellar properties in the KOI sample are less precise than the CKS sample and are derived from a hodgepodge of methods including asteroseismology, photometric bands, spectroscopy, etc. Therefore, the conclusions of this paper are based on the analysis of the CKS sample; the KOI sample merely served as a consistency check. We applied the same sample selections as in the previous section except that we removed the cut on T_eff. This resulted in a total of 1197 planets around 851 stars for KOI sample.

For both the CKS and the KOI sample, we had to rely on mass–radius relationships to estimate the masses of the planets from the transit radii. For more definite mass constraints, we used planets with transit-timing variations (TTV) or radial-velocity (RV) mass constraints. We used the TTV sample of 145 planets in 55 systems from the uniform analysis of Hadden & Lithwick (2017). We adopted the CKS stellar parameters for these systems. For the RV sample, we queried the Exoplanet Archive⁹ for the stellar and planetary parameters of planets with a ${M}_{{\rm{p}}}$ sin$i\lt 30{M}_{\oplus }$ and at least 3Σ confidence. The query resulted in 135 planets from 51 systems.

Assuming that close-in sub-Neptune planets formed in situ, we can use each planet to sample the local solid surface density of its disk. For the CKS and KOI sample, we start by converting the observed transit radii to planetary masses using various mass–radius relationships reported in the literature.

(1) Lissauer et al. (2011) proposed a simple power law based on the six solar system planets in the mass range between Mars and Saturn:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{p}}}{{M}_{\oplus }}={\left(\displaystyle \frac{{R}_{{p}}}{{R}_{\oplus }}\right)}^{2.06}.\end{eqnarray} \tag{ 2 }$$

(2) Weiss & Marcy (2014) presented a sample of 65 sub-Neptunes with orbital period <100 days. The masses of these planets were derived from both RV follow-up and TTV. Weiss & Marcy (2014) proposed a separate mass–radius relationship for super-Earths (R_p < 1.5R_⊕) and mini-Neptunes (1.5R_⊕ < R_p < 4R_⊕):

$$\begin{eqnarray}&&{R}_{{p}}\lt 1.5{R}_{\oplus }:{\rho }_{{\rm{p}}}=2.43+3.39\left(\displaystyle \frac{{R}_{{p}}}{{R}_{\oplus }}\right)\,{\rm{g}}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&1.5{R}_{\oplus }\lt {R}_{{p}}\lt 4{R}_{\oplus }:\displaystyle \frac{{M}_{{p}}}{{M}_{\oplus }}=2.69{\left(\displaystyle \frac{{R}_{{p}}}{{R}_{\oplus }}\right)}^{0.93}.\end{eqnarray} \tag{ 4 }$$

(3) Wolfgang et al. (2016) calibrated a mass–radius relationship using ∼80 sub-Neptune planets (<4R_⊕) from RV mass measurements. The analysis was done in a hierarchical Bayesian framework:

$$\begin{eqnarray}&&{M}_{{p}}/{M}_{\oplus }=2.7{\left({R}_{{p}}/{R}_{\oplus }\right)}^{1.3}.\end{eqnarray} \tag{ 5 }$$

They also reported a 1.9 ${M}_{\oplus }$ Gaussian dispersion as the intrinsic scatter in the mass distribution.

(4) The mass–radius relationship proposed by Chen & Kipping (2017) is designed to be applicable to asteroid-sized objects all the way to stars. Calibrated with more than 300 solar system and extrasolar objects, the mass–radius relationship spans more than eight orders of magnitude in mass. Specifically, in the sub-Neptune regime, Chen & Kipping (2017) drew a distinction between planets smaller or larger than ${2.0}_{-0.6}^{+0.7}{M}_{\oplus }$:

$$\begin{eqnarray}&&{M}_{{p}}\lt 2{M}_{\oplus }:{R}_{{p}}\sim {M}_{{p}}^{0.28},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{M}_{{p}}\gt 2{M}_{\oplus }:{R}_{{p}}\sim {M}_{{p}}^{0.59}.\end{eqnarray} \tag{ 7 }$$

We used the Forecaster¹⁰ code provided by Chen & Kipping (2017) to sample the masses of the CKS planets from their measured radii.

(5) More recently, using an updated list of both RV mass measurements and TTV, Mills & Mazeh (2017) proposed the following power-law relation:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{p}}}{{M}_{\oplus }}={10}^{0.45}{\left(\displaystyle \frac{{R}_{{p}}}{{R}_{\oplus }}\right)}^{0.99}.\end{eqnarray} \tag{ 8 }$$

We employed the mass–radius relationship from Wolfgang et al. (2016) for the purpose of generating figures, but our quantitative results account for the systematic uncertainties between these mass–radius relationships. For the TTV and RV samples, we skipped this step and used the reported mass constraints directly.

Given that sub-Neptunes usually have $\lesssim 10$% of their mass in H/He envelopes (Wolfgang & Lopez 2015), the bulk mass of these planets is also approximately the mass of solid materials they contain. We remind the readers that the MMEN analysis implicitly assumes that the planets themselves formed and remained near their current-day orbits, and they only accreted solid materials from their neighborhood: the feeding zone. We then calculated a solid surface density Σ assuming different widths of local feeding zones:

(1) The first prescription assumes the planet accreted solids from a feeding zone whose width is proportional to the planet's Hill radius:

$$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{{M}_{{p}}}{2\pi a{\rm{\Delta }}a},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&a={\left(\displaystyle \frac{{{GM}}_{\star }{P}_{\mathrm{orb}}^{2}}{4{\pi }^{2}}\right)}^{1/3},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}a={{kR}}_{\mathrm{Hill}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{Hill}}\equiv a{\left(\displaystyle \frac{{M}_{{p}}}{3{M}_{\star }}\right)}^{1/3},\end{eqnarray} \tag{ 12 }$$

where Σ is the solid surface density, a is the semimajor axis of the planet, and k is a constant. We chose k = 10 because Weiss et al. (2018b) found empirically that typical orbital spacing between Kepler multiplanet systems is 10 mutual Hill radii or larger.

(2) Schlichting (2014) argued that planets can have a much larger effective feeding zone if we account for giant impact collisions. This entails swapping Equation (11) with

$$\begin{eqnarray}&&{\rm{\Delta }}a={2}^{3/2}a{\left(\displaystyle \frac{a}{{R}_{p}}\displaystyle \frac{{M}_{p}}{{M}_{\star }}\right)}^{1/2}.\end{eqnarray} \tag{ 13 }$$

(3) We also used the simple prescription of Chiang & Laughlin (2013) for easier comparison with their results:

$$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{{M}_{{p}}}{2\pi {a}^{2}}.\end{eqnarray} \tag{ 14 }$$

In this prescription, each planet is spread out into an annulus of width equal to its semimajor axis.

(4) Finally, we tried the prescription of Raymond & Cossou (2014). Here, neighboring planets' feeding zones are separated by the geometric means of their semimajor axes. For the innermost and outermost planets, the feeding zone extends to ${a}_{\mathrm{in}}/1.5$ and $1.5{a}_{\mathrm{out}}$. We found that this prescription is more susceptible to missing (nontransiting) planets. With a missing planet, the neighboring planets have much larger feeding zones and in turn give rise to anomalously low Σ. We will return to discuss this point further in Section 4.6.

We found that these prescriptions do not change the qualitative conclusions of this work, it generally only affects the overall normalization of Σ. In the subsequent part of the paper, we will stick to the first prescription for its simplicity and its robustness against missing planets. In the MMSN work (Hayashi 1981), the solid surface density Σ was augmented by a gaseous component such that the resultant protoplanetary disk has solar abundance. However, recent ALMA observations (Ansdell et al. 2016) showed that the gas-to-dust ratio in a protoplanetary disk, particularly in the inner disk, deviates significantly from the 100–200 assumed for the interstellar medium (ISM). Therefore, we focused on the solid surface density Σ without augmenting it with a gaseous component.

The detection of transiting planets suffers from two major biases: transit probability and detection bias. If unaccounted for, these two effects result in a bias toward shorter-period, larger planets. The transit probability is a geometric effect. In order for an observer to see a transit, the inclination of a planet's orbit must be close to 90°. In the limit of small planets (${R}_{{p}}\ll {R}_{\star }$) and low eccentricity, the transit probability is given by

$$\begin{eqnarray}&&{p}_{\mathrm{tra}}=0.9\displaystyle \frac{{R}_{\star }}{a}.\end{eqnarray} \tag{ 15 }$$

Note that we have included a factor of 0.9 because we only considered planets with impact parameter b < 0.9.

The detection bias quantifies how complete a particular transit search pipeline is in detecting bona fide planets. For Kepler, given the complexity of the entire pipeline, the completeness of the pipeline was determined empirically with injection-recovery tests. Previous works (Petigura et al. 2013; Christiansen et al. 2015) have done extensive injection-recovery tests to characterize the behavior of the Kepler pipeline. Here, we adopt Fulton et al.'s (2017) formulae for the detection probability as a function of signal-to-noise ratio p_det (see their Equations (2) and (3)).

To debias the CKS, KOI, and TTV samples, we simply combined the transit probability and the pipeline detection probability. In our subsequent analyses, each planet in our sample was given a statistical weight inversely proportional to the combined probability:

$$\begin{eqnarray}&&w=\displaystyle \frac{1}{{p}_{\mathrm{tra}}{p}_{\det }}.\end{eqnarray} \tag{ 16 }$$

For the RV sample, we caution the readers that we could not find a simple way to correct for detection bias given the inhomogeneity of RV surveys. We chose to leave it uncorrected and assigned statistical weight only based on measurement uncertainties. In Section 4.4, we will see that the Σ inferred from the RV sample is very similar to but slightly denser than that of the transiting planets. This may relate to a detection or at least a publication bias toward heavier planets in RV surveys.

We plot the solid surface density Σ reconstructed from the CKS planets as a function of orbital distance a in Figure 2. Here, we employed the mass–radius relationship from Wolfgang et al. (2016) and assumed feeding zone widths given by prescription (1) in Section 3.2. The solid surface density follows a tight relation with the semimajor axis a similar to the power-law profile previously reported for MMSN (Hayashi 1981) and MMEN (Chiang & Laughlin 2013). This tight correlation with a partially stems from the construction of MMEN itself. If all planets had equal masses, the solid surface density would be ${\rm{\Sigma }}\propto {a}^{-2}$ with the (1) and (3) prescription of feeding zone width in Section 3.2 or ${\rm{\Sigma }}\propto {a}^{-2.5}$ with the (2) prescription.

Zoom In Zoom Out Reset image size

We color-code the planets with the host star mass M_⋆ and metallicity [Fe/H]. Visually, both M_⋆ and [Fe/H] are correlated with the residual variation in Σ. This trend is more obvious if we divide the sample into quartiles of M_⋆ and [Fe/H] in Figure 3. The metal-rich systems and the more massive systems tend to have a larger Σ. To quantify the above statement, we fit the correlations between Σ, a, M_⋆, and [Fe/H] with a simple power-law model in the form

$$\begin{eqnarray}&&{\rm{\Sigma }}={{\rm{\Sigma }}}_{0}\,{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{l}\,{10}^{m[\mathrm{Fe}/{\rm{H}}]}\,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{n},\end{eqnarray} \tag{ 17 }$$

where Σ is the solid surface density, ${{\rm{\Sigma }}}_{0}$ is a normalization constant, a is the orbital distance, [Fe/H] is the host star metallicity proxied by the iron content, and M_⋆ is the mass of the host star. The equation is easier to work with in logarithmic space:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({\rm{\Sigma }}) & = & \mathrm{log}({{\rm{\Sigma }}}_{0})+l\,\mathrm{log}(a/\mathrm{au})\\ & & +\,m\,[\mathrm{Fe}/{\rm{H}}]+n\,\mathrm{log}({M}_{\star }/{M}_{\odot });\end{array}\end{eqnarray} \tag{ 18 }$$

the set of free parameters in this model is log(Σ₀), l, m, and n. We also included a Gaussian dispersion Δlog(Σ₀) around log(Σ₀) to account for any intrinsic scatter in the solid surface density in addition to the measurement uncertainties.

Zoom In Zoom Out Reset image size

We note that when the independent variable has measurement uncertainty, the ordinary least-square estimator (OLS) tends to underestimate any correlation between the dependent and independent variables (e.g., Tremaine et al. 2002; Kelly 2007). In our case, the independent variables are a = $a({M}_{\star },{P}_{\mathrm{orb}})$, M_⋆, and [Fe/H], and they all have measurement uncertainties associated with them. It is necessary to use unbiased estimators rather than OLS. We first experimented with the orthogonal distance regression (ODR) implemented in scipy.odr, which is described in detail by Boggs et al. (1988). Kelly (2007) derived another Bayesian estimator that accounts for the measurement uncertainty in the independent variables if the independent variables can be described as a mixture of Gaussian distributions. This estimator has been successfully employed in a wide range of astronomical contexts. We coupled the likelihood from the estimators to MultiNest (Feroz et al. 2009) for both model selection and sampling parameter posterior distribution. We imposed uniform priors [−10, 10] on the parameters log(Σ₀), Δlog(Σ₀), l, m, and n. We ran MultiNest with default settings. Both estimators gave consistent results. We chose to report the results from ODR for its simplicity and because Kelly's (2007) assumption of Gaussian-mixture-independent variables does not strictly apply in this case. The results are presented in Tables 1 to 3 and discussed in the next section.

Previous works made contradictory claims on whether M_⋆ or [Fe/H] leave a stronger imprint on the observed distribution of sub-Neptune planets (see Owen & Murray-Clay 2018; Wu 2019). Using the location of the radius gap as a function of host star properties, Wu (2019) argued for a linear scaling between planetary mass and stellar mass, and no correlation between planet size and host star metallicity. However, Owen & Murray-Clay (2018) seem to disagree. They found that stellar metallicity is correlated to planetary radius (mass). More importantly, they showed for planets with <2.5 day orbits (where one expects complete photoevaporation) and planets with >25 day orbits (where one expects limited photoevaporation), it is still the case that higher-metallicity stars host larger planets. They suggest that some effects, other than the metallicity-dependent efficiency of photoevaporation, must be at work in deciding the size of sub-Neptune planets.

Here we offer an alternative and also more quantitative perspective using the MMEN framework. One caveat before going to the results is that the stellar mass M_⋆ and metallicity [Fe/H] in the CKS sample are strongly correlated with each other (Figure 4). Previous works have attributed this covariance to Galactic chemical evolution, i.e., more massive stars are younger and hence they likely formed from the more metal-enriched Galaxy (Fulton & Petigura 2018). However, several works directly looked for a possible correlation between age and metallicity in the young, thin-disk stars and did not find a compelling trend (Bensby et al. 2014; Silva Aguirre et al. 2018). While investigating this covariance further, we found that it is also related to sample selection and model degeneracy. In Figure 5, we plotted the effective temperature T_eff and mass M_⋆ of the CKS and KOI stars while color-coding the points with metallicity [Fe/H]. At the same stellar mass, a lower metallicity translates to lower opacity in the stellar atmosphere. The star would appear smaller but hotter. The variation of [Fe/H] by 1 dex results in a change of T_eff by hundreds of kelvins. The CKS sample was constructed with a hard boundary on the effective temperature 4700 K$\lt {T}_{\mathrm{eff}}\lt$6500 K; therefore, the low-metallicity stars at the higher mass end (>1.2${M}_{\star };$ [Fe/H] < 0) were excluded by the survey. An analogous effect occurs at the lower mass end. These effects contribute to the M_⋆–[Fe/H] covariance in the CKS sample.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We tried to disentangle the influence of M_⋆ and [Fe/H] by modeling them separately in Equation (18) (see also Figure 6) and used the Bayesian evidences log(Z) computed by MultiNest to gauge which of them has a stronger effect. The parameter posterior distribution and the Bayesian evidence as inferred from MultiNest are summarized in Tables 1 to 3. We can see that M_⋆ generally has a much stronger explanatory power in the variation of MMEN Σ compared to [Fe/H]. This is manifested by the steeper correlation slope m < n. Furthermore, adding the M_⋆ dependence led to several orders of magnitude improvement in the Bayesian evidence Z compared to the [Fe/H] dependence. For a summary of the results, we calculated the average MMEN profile with the CKS multitransiting planets, thanks to its uniformity and dynamical quietness (we will discuss why we excluded single-transiting systems in Section 4.2):

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Sigma }} & = & {50}_{-20}^{+33}\,{\rm{g}}\,{\mathrm{cm}}^{-2}{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{-1.75\pm 0.07}\\ & & \times {10}^{0.22\pm 0.05[\mathrm{Fe}/{\rm{H}}]}\,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1.04\pm 0.22}.\end{array}\end{eqnarray} \tag{ 19 }$$

We note that the systematic variation between different mass–radius relationships outweighs the internal uncertainties from sampling the parameter posterior distribution (see Tables 1 to 3). To capture this systematic uncertainty, we report here the MMEN profile using the mean and the standard deviation of the model parameters inferred from different mass–radius relationships (Section 3.1).

Zoom In Zoom Out Reset image size

This strong dependence on M_⋆ and the weaker dependence on [Fe/H] are seen in the KOI and RV samples, too (Tables 1 and 3). Moreover, we restrict our attention to Sun-like stars (0.9${M}_{\odot }\lt {M}_{\star }\lt 1.2{M}_{\odot };$ −0.2 < [Fe/H] < 0.2), where there is limited covariance between M_⋆ and [Fe/H] (Figure 5), and repeated the analysis. We arrived at a result similar to Equation (19), although the smaller dynamical range weakens the statistical significance. As another check of the applicability of Equation (19) across different stellar types, we repeated the MMEN analysis for the TRAPPIST-1 system: one of the most characterized, low stellar-mass planetary systems (0.09${M}_{\odot }$; Gillon et al. 2017). We arrived at a solid disk surface density ${\rm{\Sigma }}\approx 6\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ (see Figure 7) at 1 au that is about 10% of that in Equation (19), confirming the linear dependence on M_⋆.

Zoom In Zoom Out Reset image size

These observations gave us more confidence that M_⋆ plays a more important role in the MMEN framework. The ${\rm{\Sigma }}\propto {M}_{\star }^{1.04\pm 0.22}$ dependence fits the expectation of a simple in situ formation scenario. Millimeter observations of the dust continuum by Andrews et al. (2013) revealed that the dust mass of a protoplanetary disk scales more or less linearly with host star mass: ${M}_{\mathrm{dust}}\propto {M}_{\star }$. It is perhaps not surprising that this linear relation continues to the <1 au innermost disk where Kepler-like planets reside: ${\rm{\Sigma }}\propto$ ${M}_{\mathrm{dust}}\propto {M}_{\star }$.

Invoking the more sophisticated α-disk model and assuming steady accretion in the inner disks, the total surface density (${{\rm{\Sigma }}}_{\mathrm{tot}}$) is linked to mass accretion rates (${\dot{M}}_{\mathrm{acc}}$) by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}\propto \alpha {c}_{{\rm{s}}}^{2}{{\rm{\Sigma }}}_{\mathrm{tot}}{{\rm{\Omega }}}^{-1},\end{eqnarray} \tag{ 20 }$$

where α is the viscous parameter, c_s is the sound speed, and Ω is the Keplerian angular velocity. The sound speed c_s is related to the gas temperature (${T}_{\mathrm{gas}}^{0.5}$) in the inner disk regions with

$$\begin{eqnarray}&&{T}_{\mathrm{gas}}^{4}={T}_{\mathrm{acc}}^{4}+{T}_{\mathrm{irr}}^{4},\end{eqnarray} \tag{ 21 }$$

where ${T}_{\mathrm{acc}}$ is due to heating from viscous dissipation and ${T}_{\mathrm{irr}}$ is due to the stellar irradiation. If viscous dissipation is the dominant source of heating, ${T}_{\mathrm{gas}}^{4}\propto {T}_{\mathrm{acc}}^{4}\propto {M}_{\star }{\dot{M}}_{\mathrm{acc}}$ at a given disk radius (Hartmann & Bae 2018) and ${\rm{\Omega }}\propto {M}_{\star }^{0.5}$. Furthermore, from the observations made by Vorobyov & Basu (2009) and Alcalá et al. (2017), ${\dot{M}}_{\mathrm{acc}}\propto {M}_{\star }^{1.3\pm 0.3}$ for $0.2{M}_{\odot }\leqslant {M}_{\star }\lt 3{M}_{\odot }$. Thus, ${{\rm{\Sigma }}}_{\mathrm{tot}}\,\propto {M}_{\star }^{1.2\pm 0.3}$. Conversely if stellar irradiation dominates disk heating, ${T}_{\mathrm{irr}}^{4}\propto {R}_{\star }{L}_{\star }$ at a given disk radius (Dullemond et al. 2007), where R_⋆ and L_⋆ are the stellar radius and luminosity, respectively, of the central young stars. For 1–3 Myr-old stars, pre-main-evolutionary models give ${L}_{\star }\propto {M}_{\star }^{1.5}$ and ${R}_{\star }\propto {M}_{\star }^{0.4}$ for 0.6 M_⊙ $\leqslant \,{M}_{\star }\,\leqslant$1.4 M_⊙ (Siess et al. 2000), which leads to ${{\rm{\Sigma }}}_{\mathrm{tot}}\propto {M}_{\star }^{1.3\pm 0.3}$. Now, if we marginalize over [Fe/H], the MMEN solid surface density inferred from CKS multiplanet systems has a profile quite close to the predictions of α-disk models:

$$\begin{eqnarray}&&{\rm{\Sigma }}={50}_{-19}^{+32}\,{\rm{g}}\,{\mathrm{cm}}^{-2}{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{-1.76\pm 0.07}\,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1.39\pm 0.29}.\end{eqnarray} \tag{ 22 }$$

If the in situ planet formation is only limited by the availability of solid materials in the disk, one might also expect ${\rm{\Sigma }}\propto {10}^{[\mathrm{Fe}/{\rm{H}}]}$, i.e., assuming the disk solid mass scales with the host star metallicity. Alternatively, if in situ planet formation is limited by the coagulation rate of planetesimals, one might expect a steeper dependence, e.g., $\propto {10}^{2[\mathrm{Fe}/{\rm{H}}]}$, often seen in collisional problems. This latter scenario is often invoked to explain the strong correlation between giant-planet occurrence and host star metallicity (Fischer & Valenti 2005). However, we only found a sublinear relation ${\rm{\Sigma }}\propto {10}^{0.22\pm 0.05[\mathrm{Fe}/{\rm{H}}]}$. This relation confirms the picture that the formation of sub-Neptune planets occurs readily in lower metallicity ([Fe/H] ≈ −0.4) systems while the formation of giant planets strongly favors metal-rich environments (Petigura et al. 2018). As we mentioned in Section 2, our analysis excluded gas giants because we could not estimate the solid mass accurately with only the planetary radius. By excluding the giant planets (>4 R_⊕), we have excluded the most massive cores from this analysis and thus weakened the correlation between MMEN Σ and [Fe/H]. On a related note, Zhu (2019) found that the fraction of stars with sub-Neptune planets increased steadily from low- to high-[Fe/H] systems, whereas the average number of sub-Neptune planets per star increases with [Fe/H] initially but plateaus at [Fe/H] ∼ 0.1, possibly caused by the emergence of giant planets in higher-metallicity systems. Could the emergence of a giant planet dynamically disrupt its sub-Neptune companions? We will discuss this point in more detail in the next Section.

Several observations indicate that the orbital architectures of Kepler single-transiting and multitransiting systems are distinct. Xie et al. (2016), van Eylen et al. (2019) and Mills et al. (2019) arrived at consistent conclusions that single-transiting systems tend to have more eccentric orbits, while multitransiting planets favor circular orbits. Fang & Margot (2012) and Zhu et al. (2018) showed that systems with fewer planets have substantially larger mutual inclination dispersion. Moriarty & Ballard (2016) attributed this architectural difference to the different disk properties of single versus multiplanet host stars. However, later work by Weiss et al. (2018a) showed that singles and multis are very similar in terms of both planetary and stellar properties, betraying a common origin. Another possible explanation for the single–multi difference is the dynamical interaction between the sub-Neptune planets themselves or with giant planets in the same system. Zhu & Wu (2018) and Bryan et al. (2019) showed that the occurrence of Kepler-like sub-Neptune planets and cold Jupiters (>1 au) are correlated: as much as 30% of Kepler-like systems have a cold Jupiter whereas a cold Jupiter almost certainly has inner Kepler-like planetary companions. A follow-up study by Masuda et al. (2020) revealed that the mutual inclination between the cold Jupiter and inner planetary system is drastically larger if the inner planetary system only has one transiting planet. The proposed explanation is that an inclined cold Jupiter dynamically disrupts or scatters the inner planets, thereby reducing the number of planets and/or inducing large mutual inclinations and eccentricities. This dynamically hot history significantly modifies the initial orbital architecture of a system, which is manifested today as the lower multiplicity, higher mutual inclination, and higher orbital eccentricities of the single-transiting systems.

Knowing this potential architectural difference, we split the CKS sample into single-transiting and multitransiting systems and fitted their MMENs separately. The results are summarized in Table 2. The slopes of correlation between Σ, [Fe/H], and M_⋆ (m and n) are weaker in the single-transiting systems (see Figure 8). In addition, looking at the Bayesian evidence Z as a statistical measure of model improvement, including [Fe/H] and M_⋆ dependence in MMENs for single-transiting systems resulted in very limited model improvement (Bayesian evidence Z improved by no more than two orders of magnitude, or a 2σ significance). On the other hand, it is more than 15 orders of magnitude for the multitransiting systems (or >5σ significance). Our interpretation is that the dynamically hot history of single-transiting systems has at least partially eliminated the memory of formation out of the disk. Imagine that both the single-transiting and multitransiting systems form a gaseous disk. After the disk dissipates, dynamical interaction between the planets or with a distant cold Jupiter starts to incite a dynamical upheaval in some of the systems. The dynamically hotter systems have a reduced number of planets on more eccentric, mutually inclined, and scrambled orbits. These systems are naturally observed to be single-transiting and they have also lost the memory of the initial orbital architecture from in situ formation. Therefore, the singles poorly fit an MMEN framework. On the other hand, the dynamically colder systems are more likely observed as multis and largely preserved the initial orbital architecture. With the MMEN framework, we have rediscovered the different dynamical history that underpins the observed architectural differences between the single-transiting and multitransiting systems (Xie et al. 2016; Mills & Mazeh 2017; Weiss et al. 2018a; Zhu et al. 2018; van Eylen et al. 2019).

Zoom In Zoom Out Reset image size

A prominent feature displayed by sub-Neptune planets is that the occurrence rate as a function of orbital distance follows a broken power law: the occurrence rate increases steadily with orbital distance until about 0.1 au (orbital period of ∼10 days) after which the occurrence rate plateaus out to 1 au across FGKM hosts (e.g., Dressing & Charbonneau 2015; Petigura et al. 2018). In Figure 9, we reproduced this result using the CKS sample.

Zoom In Zoom Out Reset image size

The coagulation rate of planetesimals cannot be blamed for this decline of occurrence rate because the coagulation rate increases with orbital velocities and hence are higher in the inner disk (Chiang & Laughlin 2013). One possible explanation for the occurrence rate decline <0.1 au is that the amount of planet-forming solid materials may be significantly reduced at the inner disk due to dust sublimation, rapid radial drift of dusts, or some other processes. We call this possible scenario "solid disk thinning." If this were true, the inner disk should host fewer and smaller planets, and one should observe a corresponding drop in the MMEN solid surface density Σ within 0.1 au. We note that if the occurrence rate declined by about 1.5 orders of magnitude, a corresponding change in MMEN Σ would have been visually obvious in Figure 2, but was not seen.

We also performed a more quantitative test. We compared two models: (1) a single power law of Σ as a function of orbital distance ${\rm{\Sigma }}\propto {a}^{l};$ and (2) a broken power law where the power-law index l is allowed to switch to a different value below a critical distance a_crit. Using the same procedure described in Section 3, we performed model selection using the Bayesian evidence Δlog(Z) computed by MultiNest. The broken power law, with the introduction of two additional parameters, was not favored by the data (Figure 10). This statement is true for any of the mass–radius relationships in Section 3.1. The a_crit parameter is also poorly constrained rather than being pinned down to ≈0.1 au as in the occurrence rate model.

Zoom In Zoom Out Reset image size

We therefore disfavor the "solid disk thinning" scenario and argue for the alternative explanation, "solid disk truncation" (Lee & Chiang 2017). In this case, the protoplanetary disk maintains a smooth solid density profile until it is truncated by, e.g., the magnetosphere of the host star. In a truncated disk, planets can form in situ relatively undisturbed beyond the truncation radius. On the other hand, planet formation is totally quenched within the truncation radius. Because truncation occurs at the stellar corotation radius, the rotation period of the host star sets the inner edge of planet formation. For an ensemble of stars, one therefore observes a decline of occurrence rate toward the host star. On the other hand, the MMEN disk profile, thanks to closer-in planets around faster-rotating stars, remains smooth.

In short, the "solid disk thinning" scenarios should manifest as a significant drop of Σ within 0.1 au, i.e., the upper-left corner of Figure 10, while "solid disk truncation" should give rise to fewer planets in the upper-left corner that lie around a continuous trend from outer disk. This truncation scenario was discussed in detail by Lee & Chiang (2017). Their magnetosphere truncation scenario can very well reproduce the observed broken power-law occurrence rate profile with the observed rotation periods of young stars in clusters and the equilibrium tidal decay of planets after disk dissipation. We note that the disk truncation creates a pressure bump in the inner disk and the resultant inverted pressure gradient prevents the rapid loss of dusty materials due to radial drift. This maintains the solid density in the inner disk and also allows more time for the formation of close-in sub-Neptune planets.

The main conclusions of this paper are drawn from analyzing the uniform and precise CKS sample. However, as most CKS planets do not have measured masses, we had to rely on mass–radius relationships. We repeated our analyses on systems with direct planetary mass measurements: 145 TTV planets in 55 systems from Hadden & Lithwick (2017) and 135 RV planets in 51 systems from Exoplanet Archive (see Figure 7).

One word of caution before going into the results of the RV sample is that we do not have a simple way to quantify detection bias. The RV planets were derived from a host of RV surveys with different target selection criteria, observation strategies, and instrument characteristics. Rather than embarking on a tour de force of quantifying the detection bias ourselves, we simply left the detection bias of the RV planets uncorrected in this work. Nonetheless, the inferred MMEN profile from RV planets is very similar to that of the CKS multiplanets:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Sigma }} & = & {100}_{-48}^{+90}\,{\rm{g}}\,{\mathrm{cm}}^{-2}{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{-1.77\pm 0.11}\\ & & \times {10}^{0.07\pm 0.16[\mathrm{Fe}/{\rm{H}}]}\,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{0.87\pm 0.23},\end{array}\end{eqnarray} \tag{ 23 }$$

i.e., a stronger almost linear correlation with M_⋆ but a weaker dependence on [Fe/H]. We also note that the inferred solid surface density is somewhat higher in the RV sample compared to the CKS sample. The mean solid disk density is about 2σ higher in logarithmic space: log(${{\rm{\Sigma }}}_{0}/{\rm{g}}\,{\mathrm{cm}}^{-3}$) = ${10}^{1.70\pm 0.06}$ versus ${10}^{1.99\pm 0.10}$ (Table 3). One may be tempted to attribute this difference to a detection or publication bias toward heavier planets in RV surveys (Mills & Mazeh 2017). We join previous authors (Burt et al. 2018; Montet 2018) in urging the community to publish upper limits of mass constraints for statistical unbiased analyses.

For the TTV sample, we did not robustly detect the M_⋆ or [Fe/H] dependence of MMEN. The model favored by Bayesian evidence comparison only contains a dependence on orbital distance: ${\rm{\Sigma }}={72}_{-21}^{+30}\,{\rm{g}}\,{{\rm{cm}}}^{-2}{\left(a/{\rm{au}}\right)}^{-1.69\pm 0.09}$. The simple explanation is that the much smaller sample size (145 planets) and limited dynamical range in the TTV sample was not able to reveal the more subtle M_⋆ and [Fe/H] dependence. However, this is probably not the whole story. The RV sample contains a similar number of planets yet Equation (23) agreed well with the CKS sample Equation (19). We performed a bootstrap test by randomly selecting a sample of 145 planets from the CKS sample and repeating the MMEN analysis. We found that the size of the 145 TTV planets should have been more than sufficient to recover the M_⋆ or [Fe/H] dependence of MMEN as in Equation (19). Our bootstrap analysis gives a Bayesian evidence improvement of ∼7 orders of magnitudes or 4σ confidence. We thus provide an interesting but more speculative explanation that involves a different formation pathway of TTV planets. TTV signals are preferentially measured for planets near mean-motion resonances. It is proposed that these near-resonant systems formed farther out in the protoplanetary disk followed by convergent disk migration that locked the planets into resonance (Lee & Chiang 2015; Mills et al. 2016). The colder outer disk facilitates the accretion of gaseous envelopes and partially explains the observed inflated radii of planets near resonance (Millholland 2019). With such a formation pathway, TTV planets must have undergone some amount of disk migration; they are misplaced in the MMEN frame, which is an inherently in situ formation idea. The TTV planets simply do not encode the disk properties at their current-day orbits; it is hence not surprising that the Σ inferred from TTV planets do not show M_⋆ and [Fe/H] dependence.

Another interesting observation is that the occurrence rate of close-in sub-Neptune planets increases steadily toward lower mass stars (Dressing & Charbonneau 2015; Mulders et al. 2015a). In Figure 9, we reproduced this result using the CKS sample. More intriguingly, Mulders et al. (2015b) showed that the total solid mass locked in planets also increases steadily toward lower mass host stars. For FGKM stars, the amount of solid mass in planets are respectively $3.6\pm 0.1{M}_{\oplus }$, $5.0\pm 0.1{M}_{\oplus }$, $5.4\pm 0.2{M}_{\oplus }$, and $7.3\pm 0.7{M}_{\oplus }$. Here, Mulders et al. (2015b) directly summed the mass of detected Kepler planet multiplied by a correction factor equal to the inverse of the detection probability and transit probability. At this point, the readers may be confused: how does this result compare to our MMEN analysis? Briefly, our MMEN analysis spreads the mass of each planet into its feeding zone to probe the local disk surface density. We then fitted the profile of surface density with a power-law assumption. Mulders et al. (2015b) strictly only counted the amount of solid mass that is locked in existing planets, whereas the MMEN framework uses planets as anchor points to sample the solid disk profile in an in situ formation setting.

With that in mind, we can now calculate the efficiency of planet formation. We integrated the solid disk profile (Equation (19)) from 0.02 au to 0.7 au (the boundaries are chosen to match those of Mulders et al. 2015b). We found that for FGKM stars, the average total solid mass in the inner disk is $27{M}_{\oplus }$, $23{M}_{\oplus }$, $19{M}_{\oplus }$, and $11{M}_{\oplus }$, subject to a scatter of about 0.2 dex. These translate to an in situ formation efficiency of 3.6/27 ≈ 13%, 5.0/23 ≈ 22%, 5.4/19 ≈ 28%, and 7.3/11 ≈ 67% for FGKM stars. Alternatively, one can also think of this ratio as a fraction of surviving planets after dynamical evolution. In both cases, this fraction increases monotonically toward lower mass stars. We caution that the result for the M stars is mostly based on the RV sample (Equation (23)); the more uniform and more precise CKS sample has no M stars due to difficulties in modeling M-star atmospheres. This situation will be improved with the CKS-Cool project (E. Petigura et al. 2020, in preparation) that extends to lower mass stars using more empirical spectral methods. One caveat of this result is that the width of the feeding zone may somehow depend on the spectral type of the host star. If true, this would bias the inferred MMEN surface density. Moreover, we have also assumed the same innermost radius of the disk and the same radial dependence (l) for different stellar types. The first assumption on innermost radius can be partially justified by the fact that planets on extremely short orbital periods ($\approx 0.02$ au) have been observed across FGKM stars (Sanchis-Ojeda et al. 2014). We repeated our MMEN analysis in Section 3 for stars of different spectral types separately. We found that radial dependence (l) is similar across stellar types, as we assumed; it shows a standard deviation of 0.06 consistent with the estimated uncertainty in Equation (19).

Yang et al. (2020) analyzed in conjunction the occurrence and architecture (multiplicity and mutual inclination) of Kepler planets following the methods of Zhu et al. (2016). They found that for late-K, early-M stars (<4000 K), the fraction of stars hosting sub-Neptune planets is ∼55% with a typical intrinsic multiplicity of three planets, whereas the numbers drops steadily for early-type stars (∼15% F stars host Kepler-like systems with a multiplicity of ∼2). This is additional evidence for the higher efficiency of forming/retaining planets for lower mass stars. What could be the underlying reason? We know that M dwarfs are much less likely to host a giant planet (Cumming et al. 2008; Clanton & Gaudi 2014). For giant planets within <2000 day period, the suppression factor is about 3–10 times lower than for Sun-like stars. Masuda et al. (2020) showed that the presence of a misaligned giant planet dynamically disturbs the inner planetary systems, possibly reducing the number of planets and inducing higher mutual inclinations. The lack of giant planets around lower mass stars may be a boon for the close-in sub-Neptune planets, which are much likely to survive around M dwarfs.

A recent work by Moe & Kratter (2019) showed that stellar companions with a < 1 au completely suppresses the formation of S-type planets around the primary, while those at about 10 au suppresses planet occurrence to 15% that of single stars. Moe & Kratter (2019) further argue that the higher occurrence rate of sub-Neptune planets around lower mass stars can be partially accounted for by the corresponding lower binary fraction of their hosts.

Another explanation for the higher efficiency of planet formation around lower mass stars is their longer disk lifetimes. By examining the existing disk fraction across clusters of different ages, Kennedy & Kenyon (2009) suggested that lower mass stars have longer lasting disks. Silverberg et al. (2020) reported several "Peter Pan" disks around >20 Myr M dwarfs. Lee (2019) showed that the disk lifetime is a major factor in the assembly of planetesimals into planetary cores. The longer disk lifetime of lower mass stars might be instrumental for converting a larger fraction of the disk solids into sub-Neptune planets.

In a separate line of argument, there is an entire literature on the possibility of planet engulfment after planet formation (see the review by Ramírez et al. 2019). Comoving binaries are often assumed to be coeval and chemically similar; however, detailed abundance studies revealed that sometimes one of the stars can be preferentially enhanced in refractory materials, e.g., ∼0.2 dex (Oh et al. 2018). A possible explanation is the ingestion of $\sim 15{M}_{\oplus }$ worth of Earth-like materials in the star's convective layer. Planet formation and evolution can be an intrinsically lossy process wherein a significant fraction of disk solid materials do not end up in planets. We note that so far the detected possible engulfment signatures are all around G- and F-type stars (Oh et al. 2018; Nagar et al. 2019; Ramírez et al. 2019). However, the existing engulfment surveys are far from being uniform across stellar types; it is hence premature to claim whether planet engulfment happens more frequently around early-type stars.

Raymond & Cossou (2014) constructed the MMEN using Kepler multitransiting systems with three to six transiting planets. They used the geometric means of the planets' semimajor axes as the boundaries of the feeding zones (see Section 3.2). This prescription is very susceptible to nontransiting planets. Take the solar system, for example; it is very unlikely, even with infinite observation time, that any external observer could see three planets transit his/her line of sight (Wells et al. 2018). The MMEN this observer derives with the prescription of Raymond & Cossou (2014) would be drastically different from the MMSN one derives knowing the complete planetary inventory of the solar system. Applying their prescription, Raymond & Cossou (2014) found that the MMEN profile for individual Kepler systems can range between Σ $\propto \,{a}^{-3.2}$ to ${a}^{0.5}$. To correct for the effect of missing planets, Raymond & Cossou (2014) used a six-planet model with a 5° dispersion in mutual inclination to mimic Kepler observations. Their model failed to obtain a satisfactory fit to the observed Kepler systems. They hence cast doubt on a universal MMEN and argued against in situ formation. As we mentioned in Section 3.2, using the geometric mean of the planets' semimajor axes as the feeding zone is particularly susceptible to missing planets. A missing planet produces a large gap in orbital separation, which in turn produces anomalously low surface densities for its neighbors. Moreover, a 5° mutual inclination is simply too restrictive; recent works showed that the mutual inclination distribution depends on both multiplicity and orbital distance (Dai et al. 2018; Zhu et al. 2018).

Schlichting (2014) also called the MMEN construction into question. They argued that the total surface density (gas and dust, assuming a typical gas-to-dust ratio of 200 in the ISM) of the MMEN constructed from Kepler planets is close to the gravitational instability threshold (Toomre 1964). Schlichting (2014) then suggested that the majority of Kepler sub-Neptunes likely formed at >2 au before disk migration brought them in. However, the ALMA survey of Lupus by Ansdell et al. (2016) revealed that the enhancement of dusty materials in the inner disk is a ubiquitous phenomenon, with the majority of systems showing gas-to-dust ratios closer to 10 than 100. Moreover, Schlichting (2014) assumed a fixed temperature of 1000K throughout the whole disk. Here we use a more realistic temperature profile of $T\approx 150K{(a/1\mathrm{au})}^{-0.6}$ (Kenyon & Hartmann 1987; Chiang & Goldreich 1997). After plugging in the Toomre stability criterion $Q\equiv {c}_{s}{\rm{\Omega }}/\pi G{{\rm{\Sigma }}}_{\mathrm{tot}}$, the critical total surface density is about Σ_tot = 1.2 × 10⁵ g cm⁻² ${(a/1\mathrm{au})}^{-1.8}$. This is safely above the MMEN surface density (Equations (19) and (23)) with a reasonable choice of gas-to-dust ratio.

The analyses we presented in this work favored the MMEN idea particularly because the independent CKS and RV samples gave consistent results (Equation (19) and (23)). Furthermore, the concerns with MMEN pointed out by Schlichting (2014) and Raymond & Cossou (2014) seem to be largely resolved in light of the new observational results (Ansdell et al. 2016; Dai et al. 2018; Zhu et al. 2018).

How should we understand the solar system in the context of the MMEN? First of all, the MMSN Σ ≈ 7.1 ${(a/\mathrm{au})}^{-1.5}$ g cm${}^{2}$ (Hayashi 1981) is about an order of magnitude less dense than the MMEN (Equation (19)). If we include the four terrestrial planets of our solar system in the MMEN analysis, they appear roughly one order of magnitude lower than the cloud of points from Kepler (see Figure 7). Is our solar system somehow unusual? Although the Kepler spacecraft achieved exquisite photometric precision, it is probably unable to discover the analogs of our terrestrial planets. This is demonstrated both in injection-recovery tests (Christiansen et al. 2015) and in the large uncertainty of ${\eta }_{\oplus }$ calculation (e.g., Hsu et al. 2018). The detection of our terrestrial planets is even more challenging in existing RV surveys. In short, the analogs of solar system terrestrial planets are not well represented in the CKS or RV sample. The solar system may be on the less massive end of the MMEN distribution that is yet to be explored observationally. However, our solar system is not much of an outlier provided that only ∼30% of Sun-like stars host Kepler-like close-in super-Earths/sub-Neptunes (1–4R_⊕; Zhu et al. 2018). The hosts of these Kepler-like systems are probably at the higher end of the solid surface density distribution.

It is also worth noting that the formation of our terrestrial planets might have been heavily influenced by the migration of Jupiter and Saturn in the Grand Tack model (e.g., Walsh et al. 2011). Such a migration history might not have happened for many Kepler systems. In addition, if pebble accretion (e.g., Lambrechts & Johansen 2012) was important, the in situ formation narrative adopted by MMEN/MMSN analysis has to be revised in a future work.