Primordial Radius Gap and Potentially Broad Core Mass Distributions of Super-Earths and Sub-Neptunes

We begin with the ansatz that the underlying sub-Neptune/super-Earth core mass distribution follows a power-law distribution:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\mathrm{log}{M}_{\mathrm{core}}}\propto {M}_{\mathrm{core}}^{1-\xi },\end{eqnarray} \tag{ 1 }$$

where M_core is the mass of the core and we choose ξ ∈ [0.7, 1.0, 1.3]; ξ = 0.7 is the best-fit power-law slope to the distribution of peak posterior masses of sub-Neptunes from the radial-velocity follow-up by Marcy et al. (2014). Given that radial-velocity measurements are biased against low-mass objects, we do not choose ξ lower than 0.7. We note that ξ = 0.7 is top-heavy, ξ = 1.0 is neutral, and ξ = 1.3 is bottom-heavy. We also experimented with exponential distribution in linear and logarithm of M_core and found them to provide a poor match to the data. The minimum and the maximum core masses are set to 0.01 M_⊕ and 20 M_⊕, respectively. Figure 1 demonstrates the difference between the assumed power-law mass distributions with respect to the best-fit lognormal distributions presented in literature.

Zoom In Zoom Out Reset image size

For each core, its initial envelope mass fraction is calculated using the analytic scaling relationship derived by Lee & Chiang (2015) appropriate for gas accretion by cooling (equivalent to Phase II of the core accretion theory; Pollack et al. 1996; see also Ginzburg et al. 2016). We modify the expressions for the weak dependence on the nebular density (Lee & Chiang 2016) and for the expected decrease in the bound radius due to three-dimensional hydrodynamic effects (Lambrechts & Lega 2017; Fung et al. 2019). Shrinking the outer bound radius decreases the rate of accretion in a linear fashion (Lee et al. 2014; see also Ali-Dib et al. (2020) for understanding this effect in terms of entropy delivery). We verified that the expressions we provide here match the numerical calculations.

First, cores need to be sufficiently massive to accrete gas. Their gravitational sphere of influence must encompass the cores themselves so that the speed of gas is below the cores' surface escape speed. We calculate the envelope mass only for cores that satisfy

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{core}} & \leqslant & {R}_{\mathrm{out}}\equiv {f}_{R}\,\min ({R}_{\mathrm{Hill}},{R}_{\mathrm{Bondi}})\\ {M}_{\mathrm{core}} & \geqslant & 0.02\,{M}_{\oplus }{\left(\displaystyle \frac{{T}_{\mathrm{disk}}}{1000K}\right)}^{4/3},\end{array}\end{eqnarray} \tag{ 2 }$$

where ${R}_{\mathrm{core}}\equiv {R}_{\oplus }{({M}_{\mathrm{core}}/{M}_{\oplus })}^{1/4}$ (Valencia et al. 2006), R_out is the outer radius of the bound envelope, f_R < 1 is a numerical factor that takes into account the effect of three-dimensional advective flows, R_Hill is the Hill radius, R_Bondi is the Bondi radius, T_disk = 1000 K f_T(a/0.1 au)^−3/7 is the disk temperature, a is the orbital distance, and f_T is a numerical coefficient and a free parameter. We note that for these small cores, R_Bondi < R_Hill inside 1 au. We have implicitly assumed a passively heated disk (i.e., heating dominated by stellar irradiation). In Sections 3.1 and 4.1, we will explore the effect of adopting temperature profiles relevant for active disks (i.e., heating dominated by viscous accretion).

For dusty accretion, the envelope mass fraction

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{M}_{\mathrm{env}}}{{M}_{\mathrm{core}}} & = & 0.06\,{f}_{R}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{5{M}_{\oplus }}\right)}^{1.8}{\left(\displaystyle \frac{t}{1\mathrm{Myr}}\right)}^{0.4}\\ & & \times {\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{2000{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{0.12}{\left(\displaystyle \frac{0.02}{Z}\right)}^{0.4}{\left(\displaystyle \frac{\mu }{2.37}\right)}^{3.4}\end{array}\end{eqnarray} \tag{ 3 }$$

where M_env is the mass of the gaseous envelope, t is the accretion time, Σ_gas = 1.3 × 10⁵ g cm⁻² f_dep(a/0.2 au)^−1.6 is the local disk gas surface density (Chiang & Laughlin 2013), f_dep is the disk gas depletion factor, Z is the envelope metallicity, and μ is the envelope mean molecular weight. Similarly, for dust-free accretion,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{M}_{\mathrm{env}}}{{M}_{\mathrm{core}}} & = & 0.25\,{f}_{R}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{5{M}_{\oplus }}\right)}^{1.1}{\left(\displaystyle \frac{t}{1\mathrm{kyr}}\right)}^{0.4}{\left(\displaystyle \frac{200{\rm{K}}}{{T}_{\mathrm{disk}}}\right)}^{1.5}\\ & \times & {\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{4\times {10}^{5}{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{0.12}{\left(\displaystyle \frac{0.02}{Z}\right)}^{0.4}{\left(\displaystyle \frac{\mu }{2.37}\right)}^{2.2}.\end{array}\end{eqnarray} \tag{ 4 }$$

We express Equation (4) with the disk temperature T_disk. More precisely, the relevant temperature is that at the envelope radiative–convective boundary. The outer layers of dust-free envelopes are nearly isothermal so adopting T_disk obtains the same answer. We note that Equations (3) and (4) have been adjusted for ${R}_{\mathrm{core}}\propto {M}_{\mathrm{core}}^{1/4}$ compared to Lee & Chiang (2015), who used ${R}_{\mathrm{core}}\propto {M}_{\mathrm{core}}^{1/3}$; this adjustment makes no significant difference in our conclusions.

Throughout this paper, Z = 0.02 (solar metallicity), μ = 2.37, and t is drawn from a logarithmically uniform distribution that range 0.01 and 1 Myr, consistent with the late-time formation scenario (Lee & Chiang 2016). The log-uniform distribution is chosen to account for core formation by a series of collisional mergers whose doubling timescales lengthen with time (e.g., Dawson et al. 2015). Motivated by Figure 11 of Fung et al. (2019), we explore f_R = 0.1 and 0.2. We choose f_dep = 0.001 throughout, prompted by the required level of gas depletion to reproduce the observed peaks in period ratios just outside of first-order mean motion resonances (Choksi & Chiang 2020).

For a given core mass, the maximum possible envelope mass that can be accreted is given by a fully isothermal profile (e.g., Lee & Chiang 2015). No cores are allowed to accrete more than this maximally cooled isothermal mass:

$$\begin{eqnarray}&&{M}_{\mathrm{iso}}=4\pi {\rho }_{\mathrm{disk}}{\int }_{{R}_{\mathrm{core}}}^{{R}_{\mathrm{out}}}{r}^{2}\,\mathrm{Exp}\left[\displaystyle \frac{{{GM}}_{\mathrm{core}}}{{c}_{s,\mathrm{disk}}^{2}}\left(\displaystyle \frac{1}{r}-\displaystyle \frac{1}{{R}_{\mathrm{out}}}\right)\right]{dr},\end{eqnarray} \tag{ 5 }$$

where ρ_disk ≡ Σ_gasΩ/c_s,disk is the local nebular volumetric density, Ω is the Keplerian orbital frequency, c_s,disk = kT_disk/μm_H is the local disk sound speed, k is the Boltzmann constant, and m_H is the mass of the hydrogen atom. The nebular mean molecular weight μ is assumed to be the same as that of the envelope.

While the masses of sub-Neptunes are dominated by the cores, their radii are largely determined by their envelope mass fraction (Lopez & Fortney 2014). We follow closely the procedure devised by Owen & Wu (2017) in converting envelope mass fractions to radii. Only the essential elements are shown here.

First, we assume that after the disk gas is completely dissipated and planets are laid bare to stellar insolation, their outer layers become isothermal and volumetrically thin (∼6 scale heights above the radiative–convective boundary; Lopez & Fortney 2014). From the density profile given by the inner adiabat

$$\begin{eqnarray}&&\rho (r)\simeq {\rho }_{\mathrm{rcb}}\left[{{\rm{\nabla }}}_{\mathrm{ad}}\displaystyle \frac{{{GM}}_{\mathrm{core}}}{{c}_{s}^{2}}\left(\displaystyle \frac{1}{r}-\displaystyle \frac{1}{{R}_{\mathrm{rcb}}}\right)\right],\end{eqnarray} \tag{ 6 }$$

the total envelope mass

$$\begin{eqnarray}&&{M}_{\mathrm{env}}\simeq 4\pi {\rho }_{\mathrm{rcb}}{R}_{\mathrm{rcb}}^{3}{\left({{\rm{\nabla }}}_{\mathrm{ad}}\displaystyle \frac{{{GM}}_{\mathrm{core}}}{{c}_{s}^{2}{R}_{\mathrm{rcb}}}\right)}^{1/(\gamma -1)}{I}_{2},\end{eqnarray} \tag{ 7 }$$

where ρ_rcb is the density at the radiative–convective boundary (rcb), ∇_ad ≡ (γ − 1)/γ is the adiabatic gradient, γ is the adiabatic index of the interior, G is the gravitational constant, c_s ≡ kT_eq/μm_H is the sound speed evaluated at the location of the planet, ${T}_{\mathrm{eq}}\equiv {T}_{\mathrm{eff},\odot }{({R}_{\odot }/a)}^{0.5}$ is the equilibrium temperature of the planet, R_rcb is the radius at the radiative–convective boundary, and I₂ is the structure integral that follows the form

$$\begin{eqnarray}&&{I}_{n}\equiv {\int }_{{R}_{\mathrm{core}}/{R}_{\mathrm{rcb}}}^{1}{x}^{n}{\left({x}^{-1}-1\right)}^{1/(\gamma -1)}{dx}.\end{eqnarray} \tag{ 8 }$$

To eliminate ρ_rcb, we use temperature gradient at the rcb so that

$$\begin{eqnarray}&&{\rho }_{\mathrm{rcb}}=\displaystyle \frac{64\pi {\sigma }_{\mathrm{sb}}\mu {{\rm{m}}}_{H}}{3k\kappa }{{\rm{\nabla }}}_{\mathrm{ad}}\displaystyle \frac{{{GM}}_{\mathrm{core}}{T}_{\mathrm{eq}}^{3}}{L},\end{eqnarray} \tag{ 9 }$$

where σ_sb is the Stefan–Boltzmann constant, $\kappa \equiv {10}^{C}{\rho }_{\mathrm{rcb}}^{\alpha }{(k/\mu {{\rm{m}}}_{H})}^{\alpha }{T}_{\mathrm{eq}}^{\alpha +\beta }$ is the opacity at the rcb, and L is the cooling luminosity, which can be written as

$$\begin{eqnarray}&&L\simeq \displaystyle \frac{{{GM}}_{\mathrm{core}}{M}_{\mathrm{env}}}{{\tau }_{\mathrm{KH}}{R}_{\mathrm{rcb}}}\displaystyle \frac{{I}_{1}}{{I}_{2}},\end{eqnarray} \tag{ 10 }$$

where τ_KH is the Kelvin–Helmholtz cooling time of the envelope, and I₁ again follows the structure integral given by Equation (8). We vary τ_KH ∈ (100, 300) Myr. The former choice provides a good agreement between the analytically derived planetary radius presented here with the numerically computed thermally evolving sub-Neptunes by Lopez & Fortney (2014) at 100 Myr. The latter choice agrees well with the numerical solutions of Lopez & Fortney (2014) at 1 Gyr. Substituting Equation (10) into Equation (9),

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{rcb}}^{1+\alpha } & = & \displaystyle \frac{64\pi {\sigma }_{\mathrm{sb}}\mu {{\rm{m}}}_{H}}{3k}{{\rm{\nabla }}}_{\mathrm{ad}}{10}^{-C}{\left(\displaystyle \frac{\mu {{\rm{m}}}_{H}}{k}\right)}^{\alpha }{T}_{\mathrm{eq}}^{3-\alpha -\beta }\displaystyle \frac{{I}_{2}}{{I}_{1}}\\ & & \times \displaystyle \frac{{\tau }_{\mathrm{KH}}}{{M}_{\mathrm{env}}}\left(\displaystyle \frac{{R}_{\mathrm{rcb}}}{{R}_{\mathrm{core}}}\right){R}_{\mathrm{core}}.\end{array}\end{eqnarray} \tag{ 11 }$$

By rearranging Equation (7), we find another equation for ρ_rcb:

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{rcb}} & = & \displaystyle \frac{{M}_{\mathrm{env}}}{4\pi }{\left(\displaystyle \frac{{R}_{\mathrm{rcb}}}{{R}_{\mathrm{core}}}\right)}^{-3+1/(\gamma -1)}{R}_{\mathrm{core}}^{-3+1/(\gamma -1)}\\ & & \times {\left({{\rm{\nabla }}}_{\mathrm{ad}}\displaystyle \frac{{{GM}}_{\mathrm{core}}}{{c}_{s}^{2}}\right)}^{1/(1-\gamma )}{I}_{2}^{-1}.\end{array}\end{eqnarray} \tag{ 12 }$$

We numerically solve for R_rcb/R_core that obtains ρ_rcb satisfying both Equations (11) and (12), using the root_scalar function from the SciPy optimize package. Throughout the paper, we adopt γ = 7/5,⁵ , C = −7.32, α = 0.68, and β = 0.45 (Rogers & Seager 2010).⁶ To save computation time, we set R_rcb/R_core = 1 for any M_env/M_core that gives R_rcb/R_core < 1.05, motivated by the ∼5% error in Kepler transit depth measurements (e.g., Fulton & Petigura 2018). This limit can be found easily by taking the limit of ${R}_{\mathrm{rcb}}/{R}_{\mathrm{core}}\longrightarrow 1$ and confirming numerically:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{M}_{\mathrm{env}}}{{M}_{\mathrm{core}}}\right|}_{\min }=4.4\times {10}^{-5}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{{M}_{\oplus }}\right)}^{0.74}{\left(\displaystyle \frac{a}{0.42\mathrm{au}}\right)}^{0.44},\end{eqnarray} \tag{ 13 }$$

for τ_KH = 100 Myr. The numerical prefactor changes to 7 × 10⁻⁵ for τ_KH = 300 Myr.

The photospheric radius—the observable—is a few scale heights above R_rcb. Correction for the photosphere is made using

$$\begin{eqnarray}&&{R}_{\mathrm{phot}}={R}_{\mathrm{rcb}}+\mathrm{ln}\left(\displaystyle \frac{{\rho }_{\mathrm{rcb}}}{{\rho }_{\mathrm{ph}}}\right)\displaystyle \frac{{{kT}}_{\mathrm{eq}}}{\mu {{\rm{m}}}_{H}g},\end{eqnarray} \tag{ 14 }$$

where ρ_ph = (2/3)μm_Hg/kT_eqκ is the density at the photosphere and $g\equiv {{GM}}_{\mathrm{core}}/{R}_{\mathrm{rcb}}^{2}$ is the surface gravity.

Once the disk gas dissipates and the planets are laid bare to stellar insolation, those that are closest to the star are expected to lose their gaseous envelopes by hydrodynamic winds, either by photoevaporation (e.g., Owen & Wu 2013) or by internal heat (e.g., Ikoma & Hori 2012; Owen & Wu 2016; Ginzburg et al. 2018). The key difference between the two mechanisms that bear on observations is the source of insolation: whereas the former depends on the high-energy flux, the latter depends on the bolometric flux. Using this difference to validate one process over another remains to be carried out. There is a discernible shift in the position of the gap toward larger radius around more massive host stars (Fulton & Petigura 2018; Cloutier & Menou 2020; Berger et al. 2020). To reproduce this feature, photoevaporative model requires stellar-mass-dependent core mass distribution (Wu 2019), whereas this is a natural prediction of Parker wind, core-powered envelope mass-loss model (Gupta & Schlichting 2020). For solar-type stars, the two mechanisms predict similar location and shape of the gap in the radius–period distribution. Since the goal of this paper is to assess the likelihood of broad sub-Neptune core mass functions for a fixed mass of the host star, we limit our analysis to photoevaporative mass loss for simplicity. We discuss potential effects of varying stellar masses in Section 4.

Following Owen & Wu (2017), we evolve the envelope mass over 1 Gyr according to the energy-limited mass loss (e.g., Lopez & Fortney 2013)

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{env}}=-\eta \displaystyle \frac{{L}_{\mathrm{HE}}{R}_{\mathrm{phot}}^{3}}{4{a}^{2}G({M}_{\mathrm{core}}+{M}_{\mathrm{env}})},\end{eqnarray} \tag{ 15 }$$

where η = 0.1 is the mass-loss efficiency factor, and L_HE is the high-energy luminosity of the star (e.g., Ribas et al. 2005; Jackson et al. 2012)

$$\begin{eqnarray}{L}_{\mathrm{HE}}=\left\{\begin{array}{ll}{10}^{-3.5}\,{L}_{\odot } & t\lt 100\,\mathrm{Myr},\\ {10}^{-3.5}\,{L}_{\odot }{\left(\displaystyle \frac{t}{100\mathrm{Myr}}\right)}^{-1.5} & t\geqslant 100\,\mathrm{Myr}.\end{array}\right.\end{eqnarray} \tag{ 16 }$$

Orbital periods are drawn from the empirical distribution following Petigura et al. (2018)

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\mathrm{log}P}=0.52\,{P}^{-0.1}\left[1-\mathrm{Exp}\left(-{\left(\displaystyle \frac{P}{11.9\mathrm{days}}\right)}^{2.4}\right)\right]\end{eqnarray} \tag{ 17 }$$

and then converted to orbital distance assuming a solar mass host star.

We first show that late-time gas accretion alone produces a gap in the radius–period distribution (see Figure 2). The amount of gaseous envelope a core can accrete drops sharply below ∼2M_⊕ as their gas masses are limited by the maximally cooled isothermal state. The exponential dependence of this isothermal envelope mass on the core mass (Equation (5)) implies a bimodal distribution of envelope mass fractions and therefore a bimodal distribution of radii, for a smooth, underlying core mass function (see Figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 demonstrates that the location of the primordial "radius valley" shifts to smaller radii farther from the star. As the disk gets colder, planet's Bondi radius increases and so the isothermal limit rises. Figure 3 illustrates this behavior where the rocky-to-enveloped transition shifts to smaller core masses at longer orbital periods. This negative slope of the valley in the radius–period space is reminiscent of that observed (Fulton et al. 2017; Van Eylen et al. 2018). We see a larger separation between the rocky and the enveloped planetary population for dust-free gas accretion. As Figure 3 shows, this difference arises from both the generally more rapid accretion and weaker dependence on core mass for dust-free envelopes.

As we will show in the next section, gas accretion needs to be dust-free in order for the primordial radius gap (and the post-evaporation gap) to align with the observation. From numerically fitting the rocky-to-enveloped transition mass (i.e., where M_iso, Equation (5), computed numerically, crosses M_env, Equation (4)) versus orbital period, we find the transition mass ∝P^−0.31. Since ${M}_{\mathrm{core}}\propto {R}_{\mathrm{core}}^{4}$, we find the radius valley R_valley ∝ P^−0.08, in good agreement with Van Eylen et al. (2018) (within 1σ error) and Martinez et al. (2019) (within 1.5σ error). Following the same exercise but using temperature scaling expected from active disks T_disk ∝ a^−3/4, we obtain R_valley ∝ P^−0.15, which agrees with both Van Eylen et al. (2018) and Martinez et al. (2019) within 1.5σ error.

In the previous section, we showed how late-time gas accretion alone can produce a gap in the radius distribution and how the shape of this gap in radius–period space is in agreement with the observations. Next, we consider how envelope mass loss further molds the final radius distribution. Figure 4 demonstrates that the location of the radius valley carved out by photoevaporative mass loss is robustly situated at ∼1.6 R_⊕ regardless of the primordial population. As Owen & Wu (2017) cogently explain, gas-enveloped planets transform to bare rocky cores by photoevaporation when their envelope mass-loss timescale ≲100 Myr, the typical saturation timescale of high-energy luminosity of host stars. For our choice of parameters, this transition occurs for M_core ∼ 5–6 M_⊕ whose core radius is ∼1.5–1.6R_⊕ for Earth composition.

Zoom In Zoom Out Reset image size

Where the initial conditions make a difference is in the depth and the width of the gap. As illustrated in Figure 4, the narrow valley and peak in the distribution of radii are more likely to appear in dust-free envelopes (blue lines) with smaller outer radius (smaller f_R) that are assembled in hotter disks (higher f_T), and built from top-heavy core mass functions (smaller ξ). Cores that are slightly less dense than the Earth (ρ_c = 0.8ρ_⊕) reproduce better the observed location of the radius valley as the radii of bare cores puff up slightly pushing the location of the valley to ∼1.7–1.8R_⊕.

The narrowness of the radius peak for dust-free envelopes as opposed to dusty envelopes can be understood from the weaker dependence of M_env on M_core (see Equations (3) and (4) as well as Figure 3). For a given range of M_core, the confines of possible envelope mass fractions and therefore radii are more limited.

Smaller f_R reduces the maximum M_env/M_core and so keeps the primordial radius peak closer to the valley. Since photoevaporative mass loss effectively carves out the large radii peak and adds them to the lower radii, observations are better reproduced when the initial radius valley is narrower.

In hotter disks, the isothermal maximal M_env/M_core shrinks so that the rocky-to-enveloped transition appears at higher core masses. The result is a positive shift in the location of the primordial radius valley. The gas accretion rate for dust-free envelopes also reduces (see Equation (4)) and so the sub-Neptune peak shifts closer to ∼2R_⊕, bringing the primordial distribution of radii in better agreement with the observation (see the faint blue line in the top third panel of Figure 4). Since the locations of the valley are coincident with that expected from photoevaporative mass loss, we only observe a slight reduction in the peak at ∼2R_⊕ and a slight shallowing of the valley at ∼1.6R_⊕.

We observe a loss of a peak in the radius distribution when the underlying core mass function is too bottom-heavy (ξ = 1.3). While we defer detailed formal fitting of models to the data for future analyses, it is already apparent that the allowed range of ξ appears tightly constrained, under the ansatz that the core mass distribution follows a power law. It may be possible to restore the radius peak even with ξ = 1.3 with sufficiently high f_T but we judge f_T > 3 to be unlikely as it implies the disk is hot enough to melt iron at ∼0.1 au.

The combinations of parameters that provide the model radius distribution agreeing best with the observation are highlighted in Figure 5, corresponding to dust-free envelopes and f_R, f_T, ξ, ρ_c/ρ_⊕ = (0.1, 2, 0.7, 0.8). Between the primordial and evaporated population, we see a closing of the gap at ≲30 days as sub-Neptunes are whittled down to bare rocky objects by photoevaporation. This transformation also fills in the valley in a one-dimensional radius histogram. In fact, planet population that evolved from late-time gas accretion without any subsequent mass loss (right column) appears to reproduce best the observed radius valley.

Zoom In Zoom Out Reset image size

The observed radius valley closes at ∼10 days and widens toward ∼100 days (Fulton & Petigura 2018). In photoevaporation models that assume all cores have started with ≳0.01% by mass envelope, this triangular delta is hard to reproduce if the underlying core mass function is assumed to be flat (see, e.g., Owen & Wu 2013). Accounting for core mass-dependent initial envelope mass fraction alleviates the need for peaked core mass distribution as the radius gap is already carved out and photoevaporative mass loss preferentially transforms short-period sub-Neptunes into super-Earths, closing the gap inside ∼30 days (see the middle column of Figure 5).

The required disk temperatures may be uncomfortably high. We note that for accretion disks, the midplane temperature at ∼0.1 au can be as high as ∼2000K (see D'Alessio et al. 1998, their Figure 3), consistent with our f_T = 2. D'Alessio et al. (1998) accounted for both accretional and irradiation heating where the former tends to dominate in the inner disk. The expected midplane temperature from accretion depends on the accretion rate to the power of 1/4. D'Alessio et al. (1998) adopted the typical T Tauri accretion rate $\dot{M}={10}^{-8}{M}_{\odot }\,{\mathrm{yr}}^{-1}$. With f_ρ = 0.001 with respect to MMEN that we adopted here, this accretion rate can be attained with the Shakura–Sunyaev viscous parameter α ∼ 8 × 10⁻³; it is larger than what is usually inferred from studies of turbulence in protoplanetary disks such as HL tau (Pinte et al. 2016), but below the largest detected level of turbulence in, e.g., DM Tau (Flaherty et al. 2020). For optically thick disks, the midplane temperature further enhances by a factor τ^1/4 where τ is the vertical optical depth. D'Alessio et al. (1998) compute τ ∼ 400 at 0.1 au (see their Figure 2), which can be accommodated for with f_ρ = 0.001 assuming a well-mixed dust–gas mixture within the disk (see Lee et al. 2018, their Figure 1). (This enhancement of midplane temperature by τ^1/4 applies to passive disks as well.) Gas accretion onto planets can still proceed in a dust-free manner as long as the dust grains coagulate and settle within the planetary atmosphere (e.g., Ormel 2014). Figure 6 demonstrates that a qualitative result—that late-time gas accretion alone can reproduce the radius valley at the observed location—holds when we explicitly adopt the disk temperature scaling that is appropriate for a disk heated by viscous accretion (i.e., active disks).

Zoom In Zoom Out Reset image size

Matching the location of the primordial radius valley with that observed requires that the sub-Neptune cores are slightly less dense, e.g., ∼80% of the Earth composition, in rough agreement with what is reported for short-period super-Earths by Dorn et al. (2019) and more generally by Rogers & Owen (2020). We note that some of the model parameters that produce a broad peak may be narrowed by taking into account core-envelope interaction, in particular, the dissolution of gas into the magma core (Kite et al. 2019). Assessing the effect of core-envelope mixing at formation is a subject of our ongoing work.

We note that the primordial radius valley is expected to shift toward larger sizes around higher-mass stars, assuming their disks are hotter. For disks heated by stellar irradiation, ${T}_{\mathrm{disk}}\propto {T}_{\mathrm{eff}}{({R}_{\star }/a)}^{1/2}$, where ${T}_{\mathrm{eff}}\propto {M}_{\star }^{1/7}$ is the effective temperature of the star and ${R}_{\star }\propto {M}_{\star }^{1/2}$ is the radius of the star, all evaluated for fully convective, pre-main-sequence stars (e.g., Kippenhahn et al. 2012).⁸ For these passive disks, ${T}_{\mathrm{disk}}\propto {M}_{\star }^{11/28}$. From numerical fitting, we find the rocky-to-gas-enveloped transition core mass ${M}_{{\rm{c}},\mathrm{trans}}\propto {T}_{\mathrm{disk}}^{1.09}$ and ${R}_{c}\,\propto {M}_{c}^{1/4}$ so the radius valley ${R}_{\mathrm{val}}\propto {M}_{\star }^{0.11}$. For disks heated by accretion, ${T}_{\mathrm{disk}}\propto {({M}_{\star }\dot{M}/{R}_{\star }^{3})}^{1/4}{(a/{R}_{\star })}^{-3/4}$ (Clarke & Carswell 2007). Taking $\dot{M}\propto {M}_{\star }^{1.95}$ (Calvet et al. 2004), we find ${T}_{\mathrm{disk}}\,\propto {M}_{\star }^{0.7}$. Again, from numerical fitting but for active disks, ${M}_{{\rm{c}},\mathrm{trans}}\propto {T}_{\mathrm{disk}}^{2.10}$ and so ${R}_{\mathrm{val}}\propto {M}_{\star }^{0.37}$. Both estimates are within the 1σ error bar estimate by Berger et al. (2020).

We conclude this section by noting that computing disk temperature depends on uncertain factors including dust grain size distribution that can vary both radially and vertically. The disk temperatures we adopted in this paper are in the realm of possibility but more accurate comparison with observations will require better understanding of the thermal structure of the protoplanetary disks and their dependence on the host stellar mass.

Late-time gas accretion alone can carve out a valley in radius–period distribution of exoplanets. In fact, in many of the cases we explore, the agreement with the data becomes worse after taking into account photoevaporative mass loss, as the radius gap fills up. Intriguingly, David et al. (2020) report the radius gap appears more filled in around stars older than ∼2 Gyr, potentially consistent with our results as long as the dominant mass loss process operates over 1–2 Gyr, which can be either photoevaporation (e.g., King & Wheatley 2021) or core-powered mass loss (e.g., Ginzburg et al. 2018) (or both). We note however that Berger et al. (2020) report no discernible change in the depth of the gap across stellar ages of ∼1 Gyr. What the two studies have in common is that the relative number fraction of super-Earth rises around older stars, suggesting long-term mass-loss processes are likely in effect.

There remain uncertainties in the exact magnitude of the mass loss for both photoevaporation and core-powered mass-loss models. In the picture of photoevaporation, the unknown strength of planetary magnetic fields can shield against high-energy stellar photons (Owen & Adams 2019). Furthermore, there is an order of magnitude variation in the magnitude and time evolution of stellar EUV and X-ray luminosity (Tu et al. 2015). In the picture of core-powered envelope mass loss, the amount of gas mass that can be lost via wind depends on the structure of the outer envelope subject to uncertain opacity sources. Even if the cores hold enough thermal energy to unbind the entire envelope, the timescale of heat transfer depends on the unknown viscosity and Prandtl number of the magma/rocky core (e.g., Stamenković et al. 2012; Garaud 2018; Fuentes & Cumming 2020). Further advances in both theory and observations should iron out these uncertainties.