Close-in Super-Earths: The first and the last stages of planet formation in an MRI-accreting disc

TL;DR

This study models the formation and evolution of close-in super-Earths within an MRI-driven protoplanetary disc, highlighting dust dynamics, planetesimal formation challenges, and atmospheric accretion constraints compared to observations.

Contribution

It introduces a detailed model of dust evolution and planet formation in MRI-active discs, addressing the impact of turbulence on dust fragmentation and planetesimal growth.

Findings

Dust fragmentation due to MRI turbulence limits planetesimal formation.

Inner disc conditions prevent gas giant formation, aligning with observed exoplanet distributions.

Predicted planetary envelope masses are larger than those observed, indicating additional processes at play.

Abstract

We explore in situ formation and subsequent evolution of close-in super-Earths and mini-Neptunes. We adopt a steady-state inner protoplanetary gas disc structure that arises from viscous accretion due to the magneto-rotational instability (MRI). We consider the evolution of dust in the inner disc, including growth, radial drift and fragmentation, and find that dust particles that radially drift into the inner disc fragment severely due to the MRI-induced turbulence. This result has two consequences: (1) radial drift of grains within the inner disc is quenched, leading to an enhancement of dust in the inner regions which scales as dust-to-gas-mass-flux-ratio at ~1 AU; (2) however, despite this enhancement, planetesimal formation is impeded by the small grain size. Nevertheless, assuming that planetary cores are present in the inner disc, we then investigate the accretion of atmospheres onto cores and their subsequent photoevaporation. We then compare our results to the observed exoplanet mass-radius relationship. We find that: (1) the low gas surface densities and high temperatures in the inner disc reduce gas accretion onto cores compared to the minimum mass solar nebula, preventing the cores from growing into hot Jupiters, in agreement with the data; (2) however, our predicted envelope masses are still typically larger than observed ones. Finally, we sketch a qualitative picture of how grains may grow and planetesimals form in the inner disc if grain effects on the ionization levels and the MRI and the back-reaction of the dust on the gas (both neglected in our calculations) are accounted for.

Equations

$$\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[r\left({\Sigma }_{\mathrm{d}}\bar{u}-D_{\mathrm{gas}}{\Sigma }_{\mathrm{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\mathrm{d}}}{{\Sigma }_{\mathrm{g}}}\right)\right)\right]=0,$$

$$u_i=\frac{1}{1+{\text{St}}_i^2}{u}_{\mathrm{gas}}+\frac{2}{{\text{St}}_i+{\text{St}}_i^{-1}}{u}_{\mathrm{drift}},$$

$$u_{\mathrm{drift}}=\frac{c_{\mathrm{s}}^2}{2v_{\mathrm{K}}}\frac{\text{dln}P}{\text{dln}r},$$

$$t_{\mathrm{drift}}\approx \frac{1}{2}\frac{d_{\mathrm{trap}}}{u_1(r_{P_{\max }})},$$

$$\begin{array}{rl}{t}_{\mathrm{col}}& =\sqrt{\frac{8}{27\pi }}\frac{{\Sigma }_{\mathrm{g}}}{f_m{\Sigma }_{\mathrm{d}}}\sqrt{\frac{\text{St}}{\alpha }}\frac{1}{\Omega }\\ & \approx 4{\left(\frac{f_m}{0.75}\right)}^{-1}{\left(\frac{{\Sigma }_{\mathrm{d}}}{0.01{\Sigma }_{\mathrm{g}}}\right)}^{-1}{\left(\frac{\text{St}}{10^{-4}}\right)}^{1/2}{\left(\frac{\alpha }{10^{-4}}\right)}^{-1/2}\\ & \times {\left(\frac{\Omega }{10{\text{yr}}^{-1}}\right)}^{-1}\text{yr}\end{array}$$

$$\log Z_{cr2}=0.1{\log }^2\mathrm{St}+0.2\log \mathrm{St}-1.76.$$

$$X(t)$$

$$t_{\dot{X}}=\frac{4\pi d^{2}G{M}_{\mathrm{core}}^{2}X(1+X)}{\eta \pi L_{\mathrm{HE}}}\frac{1}{R_{\mathrm{p}}^3},$$

Full text

††pubyear: 2018††pagerange: Close-in Super-Earths: The first and the last stages of planet formation in an MRI-accreting disc–Close-in Super-Earths: The first and the last stages of planet formation in an MRI-accreting disc

1 Introduction

Recent advances in exoplanet detection, led primarily by the Kepler mission, have uncovered several new classes of exoplanets (e.g. Borucki et al., 2011; Batalha et al., 2013). These are the close-in super-Earths and mini-Neptunes, planets with radii of 1 – 4 R*⊕* and periods of up to 100 days, and they are found to be abundant around solar and sub-solar mass stars (e.g. Fressin et al., 2013; Dressing & Charbonneau, 2015). How (and where in their parent protoplanetary discs) these planets form is a subject undergoing intense study.

One suggestion is that these super-Earths/mini-Neptunes form at larger separations, as more solids are potentially available outside the ice line, and then migrate inwards through the disc (e.g. Ida & Lin, 2008; Kley & Nelson, 2012; Cossou et al., 2013; Cossou et al., 2014). However, the migration scenario predicts that planets in multi-planet systems should typically end up in mean motion resonances, whereas such orbital resonances are rare among the Kepler planets (e.g. Baruteau et al., 2014; Winn & Fabrycky, 2015). Although several mechanisms to either break, or prevent the capture into, resonances have been explored (e.g. Goldreich & Schlichting, 2014; Izidoro et al., 2017; Liu et al., 2017), this discrepancy has not yet been fully resolved. Moreover, the radius distribution of the Kepler planets, shaped by atmospheric photoevaporation, appears consistent with the planetary cores having a rock/iron (Earth-like) composition (Owen & Wu, 2017), implying formation inside the ice line and arguing against significant migration.

An alternative scenario is that these planets form in situ, close to their present orbits. In this case, planetary cores form in the inner protoplanetary disc. They can still be subject to planet-disc interactions and the two are not mutually exclusive. If the cores are subject to the fast type I migration, the innermost planet could stall at the inner disc edge (Ogihara et al., 2015; Masset et al., 2006), explaining why the observed period distribution of the innermost planet peaks around $\sim$ 10 days (Mulders et al. 2018; see also Mulders et al. 2015, Lee & Chiang 2017). Alternatively, type I migration could be suppressed if the surface density profile is flat or has a positive slope in the inner disc (Ogihara et al., 2018), or stalled by the core opening a gap (e.g. Hu et al., 2016). In fact, there is evidence suggesting that the Kepler planets could have been massive enough to open gaps in the inner disc (Wu, 2018).

If the close-in planets do form in situ, a large amount of solids is necessary in the inner protoplanetary disc compared to the amount in the minimum mass solar nebula (Chiang & Laughlin, 2013). These solids may be delivered to the inner disc from the outer disc prior to planet formation (Hansen & Murray, 2012, 2013; Chatterjee & Tan, 2014), through the radial drift of pebbles and rocks (Weidenschilling, 1977; Takeuchi & Lin, 2002; Armitage, 2018). The growth of dust grains in the outer disc and their subsequent radial drift inwards have been confirmed by observations (e.g. Panić et al., 2009; Andrews et al., 2012; Isella et al., 2012; Rosenfeld et al., 2013; Powell et al., 2017).

To create a dust-rich inner disc in which to form planets, the radial drift of dust particles needs to be stopped or slowed down. The radial drift of particles in the Epstein drag regime slows down closer to the star in conventional disc models, and this can concentrate dust in the inner disc to some extent (Youdin & Shu, 2002; Youdin & Chiang, 2004; Birnstiel et al., 2010, 2012; Drążkowska et al., 2016). Another way to halt the radial drift and enrich the inner disc with dust is to trap the solids inside an axisymmetric local gas pressure maximum that is expected to form if the accretion in the inner disc is driven by the magneto-rotational instability (MRI; Kretke et al., 2009; Dzyurkevich et al., 2010; Drążkowska et al., 2013; Chatterjee & Tan, 2014). A gas pressure maximum acts as a trap for the marginally-coupled solids as the gas inwards of the pressure maximum is super-Keplerian, reversing the direction of the radial drift (e.g. Pinilla et al., 2012). Furthermore, in a steady-state disc accreting due to the MRI a pressure maximum forms at the boundary between the thermally ionized innermost disc in which the MRI-induced viscosity is high, and the low-viscosity dead zone in which the MRI is suppressed due to low ionization levels (Gammie, 1996). This local pressure maximum is expected to form at few tenths of AU around solar and sub-solar stars (Chatterjee & Tan, 2014), which is consistent with the orbital distances of the close-in super-Earths and mini-Neptunes.

Mohanty et al. (2018) presented a semi-analytic steady-state inner disc model in which the disc structure, thermal ionization and the viscosity due to the MRI were determined self-consistently. The location of the pressure maximum inferred from this model is similarly in general agreement with the orbital distances of the close-in planets. An important insight from the Mohanty et al. (2018) models of the inner disc is that in a steady state they predict gas surface densities that are considerably lower than those of the minimum mass solar nebula. This is not surprising as the minimum mass solar nebula simply extrapolated the surface density to small separations, whereas in reality the shrinking size of the dead-zone results in more efficient angular momentum transport, hence lowering the surface densities towards smaller separations.

The atmospheres of many of the close-in super-Earths and mini-Neptunes must be H/He dominated (e.g. Jontof-Hutter et al., 2016) and they typically make up 0.1 – 10 % of their total mass (Lopez & Fortney, 2014; Wolfgang & Lopez, 2015). Thus, they are considerably larger than the atmospheres of the planets in the inner Solar system. Outgassing of hydrogen from a rocky core is not sufficient to explain the majority of these atmospheres (Rogers et al., 2011). Thus these atmospheres are most likely composed of gas accreted from the protoplanetary disc after the formation of a solid core. If so, these atmospheres are formed steadily through core accretion.

Lee et al. (2014) (see also Lee & Chiang, 2015, 2016; Lee et al., 2017) argue that core accretion is so efficient that the key concern is how to stop the super-Earth cores from undergoing runaway accretion and becoming gas giants (Mizuno, 1980). This led Lee & Chiang (2016) to suggest that super-Earth/mini-Neptune formation occurred in gas-poor “transition discs”, during the final short-lived phase of disc dispersal. The requirement for a gas-poor inner disc raises the question if planet formation in the gas-poor inner disc arising due to steady-state MRI accretion could be a desirable scenario and a possible alternative to the Lee & Chiang (2016) proposal of atmospheric accretion during disc dispersal.

In this paper, we examine the possibility of the in situ formation of the close-in planets in the inner disc structure arising from MRI-driven accretion, obtained using the self-consistent model of Mohanty et al. (2018). First, in section 2, we examine the evolution of dust in the inner disc and discuss the possibility of planetesimal formation. In section 3, we calculate the atmospheres that super-Earth and mini-Neptune cores can accrete in the gas-poor inner disc implied by the MRI, and then evolve them forward in time, in order to compare our calculations to the data.

An important caveat to our dust calculations is that we evolve the dust assuming the gas profile is fixed in time. Dust grains can act to suppress the MRI by lowering the coupling between the gas and the magnetic field (e.g. Sano et al., 2000; Ilgner & Nelson, 2006), and could thus significantly alter the gas disc structure. Moreover, at high dust-to-gas ratios, which we show can be achieved in the inner disc, dust becomes dynamically important and the dynamical back-reaction on the gas should be taken into account. In section 4, we qualitatively discuss these effects and how they might influence our results. We shall address the self-consistent feedback of dust enhancement on the gas disc structure in subsequent studies.

2 Dust evolution

We consider the evolution of the dust, including growth, fragmentation and radial drift, in a steady-state gas disc which is viscously accreting due to the MRI. The structure of the gas disc (gas surface density, temperature, pressure and viscosity) is obtained from the steady-state models of the inner protoplanetary gas disc calculated by Mohanty et al. (2018). In these models the viscosity ($\alpha$ parameter) is determined self-consistently with the disc structure (Shakura & Sunyaev, 1973), thermal ionization and MRI criteria (Bai & Stone, 2011; Bai, 2011). The parameters of the model are the stellar mass and radius ($M_{\ast}$, $R_{\ast}$); steady state accretion rate ($\dot{M}_{\rm g}$); and minimum viscosity of the gas due to purely hydrodynamical effects ($\alpha_{\rm DZ}$). The minimum viscosity $\alpha_{\rm DZ}$ is the assumed value of $\alpha$ inside the MRI-dead zones; it is a minimum in the sense that such hydrodynamical effects are assumed to dominate over the MRI-induced turbulence if the MRI implies a viscosity lower than $\alpha_{\rm DZ}$. In this work we primarily use $M_{\ast}=1$ M⊙, $R_{\ast}=2.33$ R⊙, $\dot{M}_{\rm g}=10^{-8}$ M⊙ $\rm{yr}^{-1}$ and $\alpha_{\rm DZ}=10^{-4}$. The disc structure for these parameters is shown in Fig. 1. The local gas pressure maximum is at an orbital distance of $\sim$ 0.7 AU, and temperature and surface density at that location are $\sim$ 1000 K and $\sim$ 5000 g cm*-2*, respectively. Outwards from the pressure (and the surface density) maximum the MRI is suppressed and $\alpha=\alpha_{\rm DZ}$.

Throughout this paper we use the above fiducial values for the disc parameters. Here we briefly describe how the disc structure depends on these parameters. For a higher gas accretion rate $\dot{M}_{\rm g}$ the radial $\alpha$ profile is roughly an outward-translated version of the one shown in the top panel of Fig. 1, and inward-translated for a smaller $\dot{M}_{\rm g}$. The radial location of the local gas pressure maximum scales with the accretion rate approximately as $\dot{M}_{\rm g}^{1/2}$. For a higher $\alpha_{\rm DZ}$ the $\alpha$ falls to this value closer to the star, and vice versa, but the value of $\alpha$ as a function of orbital distance remains almost the same otherwise; as a result, the radial location of the local pressure maximum scales with the minimum dead-zone viscosity as $\alpha_{\rm DZ}^{-1/4}$. Furthermore, inwards of the pressure maximum the temperature has to be sufficiently high for thermal ionization of potassium to support the MRI and so it is always larger than 1000 K, regardless of the exact choice of $\dot{M}_{\rm g}$ and $\alpha_{\rm DZ}$. For a steady-state vertically-isothermal $\alpha$-disc the accretion rate is $\dot{M}_{\rm g}=3\pi c_{\rm s}^{2}\alpha\Sigma_{\rm g}/\Omega$ (ignoring an additional factor which depends on the boundary condition at the inner disc edge, and which becomes unimportant far away from the edge). Further assuming (as is approximately the case) that the temperature at the location of the pressure maximum is constant regardless of the disc parameters, it follows that the maximum gas surface density approximately depends on the disc parameters111Derivation of the scalings from the Shakura-Sunyaev equations takes into account the small correction due to the dependence of the temperature on the disc parameters and yields that the maximum surface density depends on the disc parameters as $\dot{M}_{\rm g}^{3/10}$ and $\alpha_{\rm DZ}^{-13/20}$. These small corrections are omitted here for simplicity.as $\dot{M}_{\rm g}^{1/4}$ and $\alpha_{\rm DZ}^{-5/8}$.

As we neglect the dynamical back-reaction of the dust on the gas and the effect of the dust on the MRI, the local evolutionary timescales are considerably shorter than the Myr timescale on which the accretion rate will evolve. Thus, the structure of the gas disc is held fixed, i.e. not evolved in time. Note also that our gas disc model considers MRI-driven accretion only, and magnetic winds, if present, can also affect the inner gas disc structure (e.g. reduce the gas surface density compared to the minimum mass solar nebula even in the absence of the MRI, Suzuki et al., 2016).

2.1 Methods

The dust particle size distribution is evolved using the two-population model of Birnstiel et al. (2012). The dust surface density $\Sigma_{\rm d}$ is evolved using the advection-diffusion equation

[TABLE]

where $r$ is the cylindrical radius, $\bar{u}$ is the dust advection velocity, $D_{\rm gas}$ is the gas diffusivity and $\Sigma_{\rm g}$ is the gas surface density.

The dust advection velocity is a sum of the velocities due to advection with the accreting gas and radial drift. For particles with Stokes number St${}_{i}=\pi\rho_{\rm s}a_{i}/2\Sigma_{\rm g}$ (with $\rho_{\rm s}$ the internal density of the dust and $a_{i}$ the particle size), and adopting the terminal velocity approximation (e.g. Takeuchi & Lin, 2002), the dust velocity is given by

[TABLE]

where

[TABLE]

with $c_{\rm s}$ the speed of sound, $v_{\rm K}$ the Keplerian velocity and $P$ the midplane gas pressure. Small particles (St $\ll$ 1) move with the gas, and larger particles can move faster or slower than the gas, depending on the sign of the pressure gradient.

In the Birnstiel et al. (2012) model the dust surface density $\Sigma_{\rm d}$ and dust advection velocity $\bar{u}$ are the sum and the mass-weighted average, respectively, of the surface density and velocity of two populations of particles: small monomer-sized particles ($a_{\rm 0}=1$ $\mu$m) and large particles ($a_{\rm 1}$). The size of the large particles evolves in time and space. At first, small dust grains grow. Then, at each radius the size of the large particles is set by whichever process yields the smallest grain size: radial drift ($a_{\rm drift}$), where grains larger than $a_{\rm drift}$ radially drift more quickly than they can grow; drift-fragmentation ($a_{\rm df}$), where grains larger than $a_{\rm df}$ fragment due to relative radial drift velocities; or turbulent fragmentation ($a_{\rm frag}$), where grains larger than $a_{\rm frag}$ fragment due to relative velocities induced by turbulence.

As we are interested in the innermost protoplanetary disc, inside the ice line, we set the bulk density of particles to $\rho_{\rm s}=3$ g cm*-3*, and the critical fragmentation velocity to $u_{\rm f}=1$ m s*-1*, based on experiments on collisions of silicate grains (Blum & Münch, 1993; Beitz et al., 2011; Schräpler et al., 2012; Bukhari Syed et al., 2017) of similar size (the regime applicable in the Birnstiel et al. 2012 model used here).222We note that simulations of grain collisions (Meru et al., 2013) indicate that the critical fragmentation velocity could be significantly higher for porous grains than for compact ones, for a range of porosities that is not robustly covered by the above experiments. Now, in the fragmentation-limited regime particle size depends quadratically on $u_{\rm f}$. Thus, if porosity is important, and hence $u_{\rm f}$ is set to e.g. $10$ m s*-1* instead, particle sizes (Stokes number) would be larger by a factor of 100, strongly affecting how coupled a particle is to the gas flow and how susceptible to radial drift. However, since particles also become less porous (to the point of becoming compact, and fragile) in a wide variety of conditions – e.g., in collisions that result in coagulation (simulations by Meru et al. 2013, experiments by Kothe et al. 2010), collisions that result in bouncing (Weidling et al., 2009) and collisions of larger grains with monomers (Schräpler & Blum, 2011) – we opt to use the compact grain value of $u_{\rm f}=1$ m s*-1*.

The viscosity parameter $\alpha$ due to the MRI-turbulence from our gas disc model is determined as a vertical average at each radius. We assume that this vertically-averaged $\alpha$ signifies the strength of turbulence that particles feel, which in turn determines the particle size due to turbulent fragmentation ($a_{\rm frag}$) and the radial turbulent mixing (diffusivity $D_{\rm gas}$). However, the viscosity (and the level of turbulence) can be different at the disc midplane compared to the upper layers of the disc, depending on where the non-ideal magnetohydrodynamic effects suppress the MRI. Dust tends to settle towards midplane, its scale-height being determined by the balance between gravitational settling and turbulent stirring (e.g. Youdin & Lithwick, 2007). The use of a vertically-averaged $\alpha$ could thus be invalid in weakly turbulent regions. Nevertheless, we proceed with this assumption for ease of computation; we do check the robustness of our results by swapping the vertically-averaged $\alpha$ parameter for the midplane value in one run, and recover qualitatively the same results.

2.2 Numerical procedure

The advection-diffusion equation (1) is integrated using an explicit first order in time and second order in space finite element method. The advection term is integrated with an upwind scheme that adopts a van Leer flux limiter; the numerical scheme is described in detail by Owen (2014). The only modification we have made is the inclusion of the Birnstiel et al. (2012) dust evolution algorithm. Our simulations use 262 cells in the radial direction with a 0.002 AU spacing inwards of 0.4 AU and a 0.01 AU spacing outwards. The inner boundary is set to 0.016 AU, and the outer boundary is set to 1 AU (i.e. outside the pressure maximum, but inside the ice line). The time-step is set with respect to the spatial resolution, the advection speed and the diffusion coefficient, so that it obeys the Courant-Friedrichs-Lewy condition. Following Birnstiel et al. (2012), at each time step the size of the large particles is updated to the smallest of the four size limits ($a_{\rm growth}$, $a_{\rm drift}$, $a_{\rm df}$ or $a_{\rm frag}$). At the outer boundary we impose a constant dust accretion rate, which we vary in different runs. This is to mimic the fact that the dust flux is not a fixed quantity and can vary with time due to radial drift of dust from the outer disc (e.g. Birnstiel et al., 2012). Dust particle size at the outer boundary is calculated self-consistently.

Furthermore, we neglect the effects of dust sublimation, and note that the temperature in our gas disc model exceeds the dust sublimation value (\sim$$1500 K) only inwards of $0.1$ AU.

2.3 Results

Our simulations proceed as follows: initially, the dust-to-gas ratio is $0.01$ at all radii and all dust grains are monomers ($a=1$ $\mu$m). We evolve the dust for $0.2$ Myr, by which time it has reached steady state. The steady-state dust-to-gas ratio is shown in Fig. 2 as a function of radius for three different values of dust accretion rates $\dot{M}_{\rm d}$ at the outer boundary condition. For any given dust accretion rate, the steady-state dust-to-gas ratio is roughly constant inwards of the pressure maximum (indicated by the vertical dashed line), and it decreases outwards from the pressure maximum. There is only a moderate accumulation of dust at the pressure maximum, compared to the rest of the inner disc. This implies that the pressure maximum does not efficiently trap dust particles.

Essentially, the particles do not feel significant gas drag, and thus do not significantly feel the effect of the change in the sign of the pressure gradient inwards of the pressure maximum. This happens because, the MRI-induced turbulence causes fragmentation, resulting in small particles. Fig. 3 shows the three dust size limits (due to radial drift, drift-fragmentation and turbulent fragmentation) as functions of radius, calculated in steady state. The smallest of the three ($a_{\rm frag}$), due to turbulent fragmentation, sets the size of the population of large particles in these simulations, which dictates the evolution of dust overall. Particle size is thus limited to only a few millimeters near the pressure maximum (indicated by the vertical dashed line), and the particles are monomer-sized in the innermost disc.

The particle size determines, through the Stokes number, how coupled the dust is to the gas. Thus, the particle size determines to what extent the particles move with the accreting gas towards the star and also by how much they are slowed down or sped up by the gas drag. In this case, inwards of the pressure maximum (and outwards from the pressure minimum at $\sim$ 0.2 AU) the gas drag acts outwards ($u_{\rm drift}>0$). However, as dust particles are small inside the pressure maximum (St $\sim$ $4\times 10^{-4}$ even for the large particles) and their size further decreases inwards, the dust advection velocity is outwards only in a very narrow region. Consequently, after accounting for the diffusivity (i.e., the radial turbulent mixing of dust, which limits the radial gradient of the dust-to-gas ratio), the mass build-up inside the pressure maximum is moderate compared to the rest of the inner disc.

The dust advection velocity used in our method is a mass-weighted average of the velocity of the monomer-sized particles and the large particles (of size $a_{\rm frag}$, Fig. 3). The monomer-sized particles are just advected by the gas through the pressure maximum, but it can be shown that an individual large particle will also not be trapped. This is due to dust particles being in the fragmentation limit, in which the particles are fragmented faster than they drift. The drift timescale for the large particles inside the pressure trap (the region inwards of the pressure maximum where their advection velocity is outwards) can be estimated by

[TABLE]

where $d_{\rm trap}\approx 0.06$ AU is the radial width of the trap, and the particle velocity is $u_{\rm 1}\approx 2$ cm s*-1* (see eq. 2) is evaluated just inwards of the pressure maximum. The velocity $u_{\rm 1}$ decreases inwards, and so this estimate, $t_{\rm drift}\approx 7400$ yr, is a lower limit. The collisional (i.e., fragmentation) timescale for the large particles, $t_{\rm col}=(n\sigma\Delta v)^{-1}$, is much shorter in comparison. Here $n=f_{m}\rho_{\rm d}/m_{\rm d}$ is the number density of large particles, $f_{m}$ is the mass fraction of the large particles ($f_{m}=0.75$ in the fragmentation limit, Birnstiel et al., 2012), $\rho_{\rm d}\approx\Sigma_{\rm d}\Omega/(\sqrt{2\pi}c_{\rm s})$ is the midplane mass density of particles, $m_{\rm d}$ is mass of a single particle, $\sigma$ is the collisional cross section and $\Delta v\approx\sqrt{3\alpha\textrm{St}}c_{\rm s}$ is the typical relative velocity between the particles due to turbulence (Ormel & Cuzzi, 2007). This yields the collisional timescale of:

[TABLE]

where we have expressed the particle size $a_{\rm frag}$ in terms of the Stokes number St. At the pressure maximum $\alpha=10^{-4}$ and so $t_{\rm col}\approx 0.08\Sigma_{\rm g}/\Sigma_{\rm d}$ yr $\ll t_{\rm drift}$ for dust-to-gas ratios $\Sigma_{\rm d}/\Sigma_{\rm g}\gtrsim 0.01$. Thus, instead of becoming trapped in the pressure maximum, dust particles fragment and flow inwards.

In the innermost disc turbulent fragmentation yields monomer-sized particles that are entrained with the gas ($\bar{u}$ $\sim$ $u_{\rm gas}$). And so the radial drift inwards is slowed by particles becoming well coupled to the gas. This is also why, in steady state, the dust-to-gas ratio is roughly constant inwards of the pressure maximum.

Finally, despite the pressure maximum not trapping the inflowing particles, the dust-to-gas ratio is enhanced. Because the dust moves with the gas in the innermost disc, the steady state dust-to-gas ratio there is directly proportional to the ratio of dust-to-gas accretion rates at the outer boundary, $\dot{M}_{\rm d}/\dot{M}_{\rm g}$. Now, the initial dust-to-gas ratio everywhere is $\Sigma_{\rm d}/\Sigma_{\rm g}=10^{-2}$. Preserving this ratio over time in the inner disc would thus require $\dot{M}_{\rm d}/\dot{M}_{\rm g}=10^{-2}$ at the outer boundary. However, the growth of dust grains in the outer disc, and the attendant inward radial drift of grains there, means that $\dot{M}_{\rm d}/\dot{M}_{\rm g}>10^{-2}$ at the outer boundary (i.e., dust accretes inwards preferentially compared to gas). In this case, the inner disc in steady-state will also have $\Sigma_{\rm d}/\Sigma_{\rm g}>10^{-2}$. In other words, as Fig. 2 shows, radial drift of grains from the outer disc leads to an enrichment of solids in the inner disc.

What level of the enrichment is attainable depends on the ratio of dust and gas accretion rates $\dot{M}_{\rm d}/\dot{M}_{\rm g}$, i.e. how quickly the grains drift from the outer disc relative to gas accretion. Since the grain growth in the outer disc is limited by radial drift rather than fragmentation (Birnstiel et al., 2012), high grain drift rates are possible. For example, assuming that $\Sigma_{\rm d}/\Sigma_{\rm g}=10^{-2}$ in the outer disc, achieving $\dot{M}_{\rm d}/\dot{M}_{\rm g}=1$ requires the radial drift velocity of grains ($\approx 2\textrm{ St }u_{\rm drift}$) to reach $10^{2}u_{\rm gas}$. For the standard $\alpha$-disc model and $\alpha=10^{-4}$ (extrapolation of the disc model shown in Fig. 1) this is satisfied if grains grow to St $\sim$ $10^{-2}$. This corresponds to a particle size an order of magnitude below the radial drift limit ($a_{\rm drift}$) throughout the outer disc. Therefore, the grains easily grow large enough to achieve accretion rates of $\dot{M}_{\rm d}/\dot{M}_{\rm g}\gtrsim 1$. These grains in the outer disc will contain ices which will evaporate as grains drift across the ice lines, towards the inner disc. However, even in the outermost disc silicates account for a considerable portion of the total solid mass (e.g., adopting the abundances of oxygen and carbon in their main molecular carriers from Öberg et al., 2011, 23 % of the total oxygen and carbon mass in the outer disc is in silicates and other refractories). Hence, Fig. 2 features dust accretion rates up to $\dot{M}_{\rm d}/\dot{M}_{\rm g}=1$ (which, as argued above, are likely), in which case the dust-to-gas ratio in the inner disc also approaches unity.

2.4 Implications for planetesimal formation

The above results show that the MRI yields a dust-enhanced inner disc, although at the expense of the dust particle size. Furthermore, as there is no trap for the dust particles, the accumulation of dust is limited by the dust inflow rate from the outer disc, and does not increase indefinitely. In this section we explore if further concentration of particles via the streaming instability (SI; Youdin & Goodman, 2005) and subsequent gravitational collapse into planetesimals could be the next step towards forming the close-in super-Earths and mini-Neptunes.

The SI can greatly concentrate dust particles if the ratio of dust-to-gas bulk densities is $\rho_{\rm d}/\rho_{\rm g}$ $\gtrsim$ 1 (Youdin & Goodman, 2005; Johansen & Youdin, 2007). This is most likely to be attained in the disc midplane, as dust particles gravitationally settle. The settling is balanced by turbulent stirring. One source of turbulence is the MRI. To reach $\rho_{\rm d}/\rho_{\rm g}\geq 1$ in the midplane in the presence of such turbulence, the dust-to-gas surface density ratio $\Sigma_{\rm d}/\Sigma_{\rm g}$ needs to be greater than or equal to $Z_{\rm cr1}=\sqrt{\alpha/(\rm{St}+\alpha)}$ (Carrera et al., 2017).

Even in discs that are weakly turbulent or completely laminar, as dust settles the dust-gas interaction leads to turbulence (self-stirring) which can prevent clumping by the SI. In this case, SI can only successfully concentrate dust particles if the dust-to-gas ratio $\Sigma_{\rm d}/\Sigma_{\rm g}$ is greater than a critical value $Z_{\rm cr2}$ that depends on the particle Stokes number (Johansen et al., 2009; Carrera et al., 2015). For St $<0.1$ (relevant to our simulations) this critical value has been most recently revised by Yang et al. (2017), who find

[TABLE]

Small dust grains that are entrained with the gas do not participate in the SI. Hence, to compare our results against the SI criteria, we use the dust-to-gas ratio $f_{\rm m}\Sigma_{\rm d}/\Sigma_{\rm g}$ of the large grain population only (where $f_{\rm m}=0.75$ is the mass fraction of large particles; Birnstiel et al., 2012).

In the top panel of Fig. 4 we compare this dust-to-gas ratio, for the outer boundary condition $\dot{M}_{\rm d}/\dot{M}_{\rm g}=1$, with the above two criteria for the onset of the SI. We find that in the inner disc, turbulence due to the MRI is generally more prohibitive to dust settling than the turbulence due to dust-gas interactions. Both conditions are fulfilled only near the pressure (and density) maximum.

Provided that the SI successfully concentrates dust particles in the disc midplane, the dust bulk density there may reach up to 100 – 1000 times the local gas density (Johansen & Youdin, 2007). Gravitational collapse of such particle concentrations will occur if the dust density exceeds the local Roche density (below which the star can tidally disrupt the fragment), $\rho_{\rm Roche}=9M_{\ast}/4\pi r^{3}$. Comparing $\rho_{\rm Roche}$ with the midplane gas densities in our steady-state MRI disc (Fig. 4, bottom panel), we find that the condition for gravitational collapse is only satisfied near the pressure (and density) maximum, or at larger separations. Importantly, the bottom panel of Fig. 4 shows that the possibility of gravitational collapse of solids in the inner disc is severely limited by the steep gradient of the Roche density.

Overall, given the stellar and disc parameters used here, we find that the SI and the gravitational collapse pathway to planetesimals is viable in the inner disc only in a very narrow region near the pressure maximum. For gas accretion rates larger than the one used here ($\dot{M}_{\rm g}>10^{-8}$ M⊙ yr*-1*, that could be expected in the early phase of disc evolution, e.g. Manara et al., 2012) this conclusion will hold, while for sufficiently smaller accretion rates planetesimals would not form in this way anywhere near or inwards of the pressure maximum. Firstly, in both cases dust evolution is expected to yield roughly the same steady-state dust-to-gas ratios given the same ratio of dust and gas accretion rates at the outer boundary. To confirm this conclusion, we repeat the dust evolution calculations for the gas accretion rate of $\dot{M}_{\rm g}=10^{-9}$ M⊙ yr*-1*, obtaining results very similar to those above. In addition to similar steady-state dust-to-gas ratios, the particle Stokes number sharply drops inwards of the pressure maximum, and the SI is similarly triggered only around the pressure maximum. This is because the slope of the increase in $\alpha$ inwards of the pressure maximum is roughly the same for different $\dot{M}_{\rm g}$. Hence the SI criteria is expected to be fulfilled only near the pressure maximum for higher gas accretion rates as well.

Secondly, even if SI is successfully triggered, to form planetesimals the peak dust density needs to be above the Roche density. The peak dust density scales with the midplane gas density at the pressure maximum, so it scales with the gas accretion rate approximately as $\dot{M}_{\rm g}^{-1/2}$, and with the radial location of the pressure (and density) maximum as $r_{P_{\rm max}}^{-1}$. Since the peak dust density is larger than the Roche density at the pressure maximum for $\dot{M}_{\rm g}=10^{-8}$ M⊙ yr*-1*, and the Roche density $\rho_{\rm Roche}\propto r^{-3}$, for $\dot{M}_{\rm g}>10^{-8}$ M⊙ yr*-1* the gravitational collapse criterion will also be fulfilled near the pressure maximum, while for a sufficiently smaller accretion rate (including $\dot{M}_{\rm g}=10^{-9}$ M⊙ yr*-1*) peak dust density will be too low.

Moreover, it is important to note that, in order for the SI to operate, there must be a relative azimuthal velocity between the dust and gas (e.g. Squire & Hopkins, 2018), in addition to the density criteria invoked above. Such a relative velocity necessarily disappears at the pressure maximum itself, further constricting the region where the SI is viable in our disc.

However, we reiterate that our calculations do not include the effects of dust on the gas dynamics and on the MRI. The latter effect in particular may relax some of the constraints on the SI, and we discuss this possibility in section 4. Overall, if the close-in sub-Neptunes form in situ, in an MRI-accreting inner disc (as suggested by Chatterjee & Tan, 2014), any theory needs to explain how the dust grains can grow larger to either start the SI in a wider vicinity of the pressure maximum, or for the pressure maximum to become an efficient dust trap that can concentrate grains to densities needed for the gravitational collapse, or find an alternative pathway to planetesimal and core formation.

3 Accretion and evolution of planetary atmospheres

We showed above that planetesimal formation through streaming instabilities and gravitational collapse is challenging despite the dust enhancement in the inner disc. Although these conclusions could change (see section 4), it is not presently clear how exactly planetesimals or cores would arise in the inner disc, and we are unable to predict properties of solid cores formed in situ. Nevertheless, orbital distances, radii and masses of many close-in planets have been well determined observationally. Furthermore, their radius distribution has been shown to be consistent with Earth-like composition and thus formation inside the ice line (Owen & Wu, 2017), possibly in the inner disc. Thus, in this section we use observational results to perform a separate test of in situ formation by considering accretion of planetary atmospheres in the inner disc and their subsequent evolution.

To follow the accretion of planetary atmospheres, we assume that solid super-Earth sized cores are present in the inner regions of our gas-poor inner disc structure taken from Mohanty et al. (2018). Furthermore, following the arguments of Lee & Chiang (2016), we ignore additional heating arising from further accretion of solids, as any amount of accretion providing significant heating typically results in the core rapidly reaching run-away masses.

After disc dispersal (which we do not model), we also account for atmospheric mass lost due to photoevaporation. This is an important addition, as several theoretical (e.g. Owen & Wu, 2013; Lopez & Fortney, 2013) and observational (e.g. Lundkvist et al., 2016; Fulton et al., 2017; Fulton & Petigura, 2018) studies have shown that photoevaporation significantly sculpts the exoplanet population after formation.

3.1 Methods

3.1.1 Accretion of planetary atmospheres

We assume that solid super-Earth-sized planet cores accrete their gaseous envelopes in a gas-poor inner disc that is viscously accreting due to the MRI; this disc structure does not evolve in time. For a quasi-steady state envelope the accreted envelope mass fraction after time $t$ can be estimated by the scaling relations from Lee & Chiang (2015) (with an additional factor that accounts for varying gas surface density from Lee et al., 2017; Fung & Lee, 2018):

[TABLE]

for dusty atmospheres, and

[TABLE]

for dust-free atmospheres. Here $Z$ is the metallicity of the atmosphere; $\mu=1/(0.5W+0.25Y+0.06Z)$ is the mean molecular weight, with $W=(1-Z)/1.4$, $Y=0.4(1-Z)/1.4$; $T_{\rm rcb}$ is the temperature at the radiative-convective boundary inside the atmosphere; $f_{\Sigma}=\Sigma_{\rm g}/\Sigma_{\rm MMSN}$ is the ratio of the gas surface density ($\Sigma_{\rm g}$ from our inner disc model, Fig. 1) and the gas surface density profile of the minimum mass solar nebula ($\Sigma_{\rm MMSN}=1700(d/1\textrm{\,AU})^{-3/2}$ g cm*-2*, where $d$ is the orbital radius, Hayashi, 1981). Furthermore, we have assumed the gas adiabatic index is $\gamma=1.2$, and that in dusty atmospheres $T_{\rm rcb}=2500$ K which arises from the disassociation of Hydrogen (see Lee & Chiang (2015), their section 2.1).

We use the expressions (7,8) to calculate how much gas a planet accretes in 1 Myr as a function of core mass, for various metallicities $Z$ and gas surface density factors $f_{\Sigma}$ in the case of dusty atmospheres (with $T_{\rm rcb}$ = 2500 K), and various $Z$, $f_{\Sigma}$ and $T_{\rm rcb}$ in the case of dust-free atmospheres.

3.1.2 Photoevaporation of planetary atmospheres

These accreted atmospheres are then subjected to photoevaporation following disc dispersal. We use a simplified estimate of the photoevaporative mass loss. First, for a given planet core mass $M_{\rm core}$ and (accreted) envelope mass fraction $X$ we find the photospheric radius of the planet $R_{\rm p}$. For this we use a simple model of an atmosphere at hydrostatic equilibrium (Owen & Wu, 2017), in which the solid core is surrounded by an adiabatic convective envelope, topped by an isothermal radiative photosphere.

Next, the mass-loss timescale due to high-energy stellar irradiation is (Owen & Wu, 2017)

[TABLE]

where $d$ is the orbital radius of the planet and $L_{\rm HE}$ is the stellar high-energy flux. We consider a Sun-like star, as in our disc model above.

To determine the final envelope mass fractions, we do not explicitly evolve the atmospheres in time. Instead, we use the fact that most of the mass loss happens during the first $\sim$ 100 Myr after disc dispersal, since during this period the stellar high-energy flux $L_{\rm HE}$ is saturated ($L_{\rm HE}=L_{\rm sat}\sim 10^{-3.5}L_{\odot}$ for a Sun-like star) while after this time it quickly decays (Jackson et al., 2012; Tu et al., 2015).

Thus, if a planet’s mass-loss timescale $t_{\dot{X}}$ is initially (i.e., at the time of disc dispersal) longer than 100 Myr, the planet will not suffer significant mass loss. Here we assume that such a planet remains unchanged by the photoevaporation.

On the other hand, a planet with initial $t_{\dot{X}}<100$ Myr will lose mass. Now, for a given core mass $M_{\rm core}$ and orbital distance $d$, the mass-loss timescale as a function of the envelope mass fraction, $t_{\dot{X}}(X)$, peaks at $X\equiv X_{\rm peak}$ of a few percent, decreasing for both smaller and larger $X$ (Owen & Wu, 2017). Thus, for a planet with an initially small accreted envelope mass fraction ($X<X_{\rm peak}$), the mass loss further shortens the loss timescale. Hence, if such a planet’s initial mass-loss timescale is less than 100 Myr it is subject to runaway mass loss, and we assume it is completely stripped of its atmosphere. Conversely, the mass-loss timescale of a large atmosphere ($X>X_{\rm peak}$) increases as it loses mass, tending towards the peak value of $t_{\dot{X}}(X_{\rm peak})$. If $M_{\rm core}$ and $d$ are such that the latter peak timescale is $\geq$100 Myr, we assume that a planet with such a large initial atmosphere will stall at an envelope mass fraction $X$ corresponding to a loss timescale of $t_{\dot{X}}$ = 100 Myr. However, if the peak timescale is $<$100 Myr, the mass-loss timescale can increase to this peak value and then descend into the runaway regime on the other side while the stellar activity is still saturated; thus, such a planet will lose its entire atmosphere regardless of the accreted $X$. This simple prescription adequately captures the basic physics of atmospheric photoevaporation (see far left panel of Fig. 6 in Owen & Wu, 2017).

3.2 Results

3.2.1 Accretion of planetary atmospheres

Using the scaling relations (7,8) we calculate the envelope mass fractions that planetary cores of various masses accrete from the gas disc in 1 Myr. Results are shown in Fig. 5 for both dusty and dust-free atmospheres of various metallicities, ranging from solar ($Z=0.02$) to the metallicity of Neptune’s atmosphere ($Z=0.2$; Karkoschka & Tomasko, 2011), gas surface density factors $f_{\Sigma}=10^{-4}$ – 1.88, and radiative-convective boundary temperatures $T_{\rm rcb}=1600$ – 2500 K.

The dependence on metallicity $Z$ is non-monotonous, as the accreted envelope mass fraction depends separately on the metallicity and on the mean molecular weight (which is set by the metallicity). The smallest accreted atmospheres, with the rest of the parameters fixed, have $Z$ $\sim$ $0.1$.

The radiative-convective boundary temperature $T_{\rm rcb}$ is expected to be roughly constant in dusty atmospheres, so we only explore the effect of this parameter in dust-free atmospheres. In the latter, $T_{\rm rcb}$ is related to the temperature of the environment. Additionally, accretion of both dust-free and dusty atmospheres depends the density of the environment. Here we are interested in atmospheres that are accreted in the inner disc, near or inwards of the pressure maximum. The location of the pressure maximum is determined by the extent of thermal ionization of potassium in our disc model, and so this corresponds to disc temperatures of $T$ $\gtrsim$ 1000 K, regardless of the exact disc parameters (e.g. gas accretion rate). For a disc temperature of $T$ $\sim$ 1000 K, numerical models of the accreting atmospheres give $T_{\rm rcb}$ $\sim$ $1600$ K (Lee & Chiang, 2015), which thus sets a lower bound on $T_{\rm rcb}$ for our calculations. Moreover, the location of the pressure maximum is also where the gas surface density is highest (see Fig. 1). Colder atmospheres in a more dense environment accrete more. So, to show the maximum accreted atmospheres in an MRI-accreting disc, we plot a set of dusty and dust-free atmospheres (of various metallicities $Z$) for the maximum $f_{\Sigma}=1.88$, and the minimum $T_{\rm rcb}=1600$ K (the latter refers only to the dust-free atmospheres). Conversely, the maximum temperature at which equation (8) is valid (due to the limitations of the opacity tables used by Lee et al., 2014) is $T_{\rm rcb}=2500$ K333Note that, assuming $T_{\rm rcb}$ is directly proportional to the disc temperature $T$, and scaling from the numerical models’ result that $T_{\rm rcb}$ $\sim$ 1600 K corresponds to $T$ $\sim$ 1000 K, yields $T$ $\sim$ 1500 K for $T_{\rm rcb}$ $\sim$ 2500 K. The disc temperature in our model only exceeds 1500 K at radii $<$0.1 AU, so a maximum $T_{\rm rcb}$ of 2500 K is indeed roughly valid over most of our inner disc., and the minimum gas surface density in our inner disc model with respect to the minimum mass solar nebula is $f_{\Sigma}\approx 10^{-4}$ (corresponding to the inner disc edge in Fig. 1). Hotter atmospheres in lower-density environments accrete less, and so to show the smallest accreted dust-free atmospheres, we plot the $f_{\Sigma}=10^{-4}$ atmospheres, with $T_{\rm rcb}=2500$ K for the dust-free atmospheres, and with metallicity $Z=0.1$ (since, as noted above, $Z$ $\sim$ 0.1 yields the smallest atmosphere for any given set of other parameters).

Finally, atmospheres that grow above a threshold of $X=0.5$ undergo runaway accretion and end up as gas giants (Rafikov, 2006). The scaling relations (7,8) are not applicable in this case, neither are we interested in the much rarer close-in Jupiters. Therefore, Fig. 5 is cut off at $X=0.5$, and the grey region indicates how small or large super-Earth/mini-Neptune atmospheres may be at the time of disc dispersal. Overall we see that, if the cores are formed 1 Myr before the dispersal, runaway accretion is avoided for the majority of relevant core masses, but they do accrete significant gaseous envelopes of up to a few $\times$ 10% of core mass.

The envelopes shown in Fig. 5 have been calculated assuming that the accretion lasts for 1 Myr. Since disc lifetimes can be longer (Mamajek, 2009), these envelopes could be conservative estimates if planets form sooner than 1 Myr before disc dispersal. If, for example, the envelopes are accreted for 5 Myr, the envelope mass will double. We do not expect the results to be very sensitive to the exact disc parameters, as long as the cores accrete their atmospheres in a thermally-ionized MRI-active inner disc. We note, however, that our disc model implies the extent of such inner disc does not encompass all observed sub-Neptunes for all relevant accretion rates; e.g. the gas pressure maximum is at an orbital period longer than 100 days only for gas accretion rates of $\dot{M}_{\rm g}$ $\gtrsim$ $3\times 10^{-9}$ M⊙ yr*-1*. Thus, a planet with a longer orbital period might spend at least some time in a colder MRI-dead zone, which we do not take into account.

3.2.2 Photoevaporation of planetary atmospheres

To further check the consistency of core accretion in the MRI-accreting inner disc with observations, we need to consider whether these accreted atmospheres survive photoevaporation. We calculate the final (remaining) envelope mass fraction of the minimum and maximum possible accreted atmospheres (corresponding respectively to dusty atmospheres with $f_{\Sigma}=10^{-4}$ and $Z=0.1$, and dust-free atmospheres with $Z=0.02$, $f_{\Sigma}=1.88$ and $T_{\rm rcb}=1600$ K) for each core mass and as a function of orbital period. Results are shown in Fig. 6. In the case of the maximum accreted atmospheres (top panel), the atmospheres would undergo runaway accretion for core masses $\gtrsim$ 8 M⊕ (indicated by the hatched region), which are thus excluded here.

The figures show that the orbital period at which the atmosphere can be completely evaporated decreases with increasing core mass, and cores that retain their atmospheres generally evolve towards a 1% envelope mass fraction as expected from theory. At 100 days the atmospheres are unaffected by photoevaporation, and at periods shorter than 1 day all planets are predicted to end up as bare cores. Massive cores are predicted to keep their 1 – 50% atmospheres at the majority of orbital periods, and planets with Earth-mass cores are safe from complete mass loss at periods larger than 50 days.

Note that here the orbital period determines the level of high-energy flux that planet experiences and planet equilibrium temperature (and thus planet radius), but does not directly reflect variations in temperature and density of the protoplanetary disc inside which the atmospheres were accreted. As discussed above, the effect of the disc temperature on the accreted envelope mass fraction is negligible for dusty atmospheres. For dust-free atmospheres the dependence is monotonous and the extent of the effect is explored by considering the minimum $T_{\rm rcb}$ we expect in the inner disc, and the maximum $T_{\rm rcb}$ for which the scaling relations (7,8) are valid. Similarly, the dependence on the disc density is explored by considering the smallest and largest values of the ratio of the MRI-disc model and the minimum mass solar nebula surface densities. Thus, by calculating the effect of photoevaporation on both the minimum and maximum accreted atmospheres shown in Fig. 5 for each core mass, we also encompass the possible range of disc densities and temperatures.

3.3 Comparison to observations

Overall, Fig. 6 shows that the envelopes formed in a gas-poor inner disc due to the MRI survive photoevaporation for a large range of orbital periods, and the low gas surface densities are not a hindrance to the formation of mini-Neptunes. On the contrary, the final envelope mass fractions of the planets that do keep their atmospheres are typically overestimated. The planets with core mass larger than 2 M⊕ are predicted to either have a $>1\%$ atmosphere or to be completely evaporated. On the other hand, from the observations, the typical envelope mass fraction of mini-Neptunes that hold onto their atmospheres is 1% (Wolfgang & Lopez, 2015). To look into this further, we compare the predictions of our calculations against the observed mass-radius relationship for sub-Neptune planets in Fig. 7, and against measured masses and radii of individual sub-Neptune planets in Fig. 8.

For the observations in Fig. 7, we show the probabilistic best-fitting mass-radius relationship of Wolfgang et al. (2016): a power law $M/\textrm{M}_{\oplus}=2.7(R/\textrm{R}_{\oplus})^{1.3}$ (indicated by the solid grey line) with a standard deviation of $\pm 1.9$ M⊕ due to an intrinsic scatter in planet mass (the grey region), and an upper limit constraint on the planet density corresponding to a mass-radius relationship for solid cores of Earth-like composition $M/\textrm{M}_{\oplus}=(R/\textrm{R}_{\oplus})^{4}$ (dash-dot-dash line; Valencia et al., 2010). Additionally, the above mass-radius relationship does not capture a significant feature of the observed radius distribution of sub-Neptunes, a decrease in occurrence rates of planets with radii of 1.5 – 2 R⊕ (indicated here by the sheer grey region; Fulton et al., 2017).

To show the predictions of our atmospheric calculations in the mass-radius plane, we take the calculated envelope mass fraction as a function of core mass and period and re-calculate the planet radius at the planet age of 5 Gyr as a function of core mass and period, using the same simple atmospheric evolution model of Owen & Wu (2017). We show the results for the minimum and maximum accreted atmospheres (Fig. 6, and excluding the completely evaporated planets) in Fig. 7 (the yellow and the light yellow region respectively). Note that the light yellow region has a cut-off at about 8 M⊕ because we exclude massive cores that, given the parameters of the maximum accreted atmospheres, would be subject to runaway accretion. The solid-line contours show how the planet mass and radii change as a function of period for the minimum accreted atmospheres. At the orbital period of 100 days the planets are largely unaffected by the atmospheric loss, and closer to the star the photoevaporation removes atmospheres of the lower-mass planets entirely. For the planets that keep their atmospheres at large periods a decrease in period means little to no change in planet mass. Consequently, for these planets a decrease in period results in an increase in planet radius as atmospheres are hotter and more expanded closer to the star due to stronger stellar irradiation. At small periods the atmospheric loss is significant for all planets, and the trend is reversed.

It is clear from Fig. 7 that for planets with radii $R\lesssim 2.3$ R⊕ the core accretion of atmospheres in the inner disc predicts larger planet radii than those observed, due to the overestimated envelope mass fractions. The predicted atmospheres are massive enough to populate the region corresponding to planet radii of 1.5 – 2 R⊕, which is inconsistent with the observed decrease in planet occurrence rates at those radii (sheer grey region). For planets with $R\gtrsim 2.3$ R⊕ there is a region in which the observed (grey) and the predicted (yellow) mass-radius relationships overlap. This overlap corresponds to the (minimum accreted) predicted atmospheres for orbital periods between 20 and 100 days, and a narrow range of short orbital periods (2 – 5 days). Notably, even for the minimum accreted atmospheres, the planet radii, at fixed planet mass, are smaller than those observed only for significant high-energy fluxes at orbital periods of less than about 2 days. Taking into account the full range of accreted atmospheres (up to the maximum accreted atmospheres shown in light yellow) further suggests that the predicted atmospheres are typically larger than the atmospheres of the observed sub-Neptunes.

We further compare the predictions of our calculations to sub-Neptune planets with measured masses and radii (taken from Wolfgang et al. 2016, excluding the planets where only the upper limit on the mass was known). The observed and the predicted radii and masses are shown in Fig. 8 in four panels corresponding to four orbital period bins. As in Fig. 7, the yellow and light yellow regions correspond to the predictions from the minimum and maximum atmosphere mass models respectively. To facilitate comparison against the planets that are bare solid cores in each period bin, here we also show the core masses that are predicted to lose their entire atmospheres in a given period bin (yellow lines shown below the dash-dot-dash lines that represent the Earth-like composition mass-radius relationship). Fig. 8 shows that the masses of the predicted bare cores and the period at which photoevaporation can strip them are largely consistent with those observed. That is, there are no observed planets consistent with the Earth-like composition that are (significantly) more massive than the largest core that the photoevaporation can strip (the upper limit of the yellow line, the predicted bare cores) in each period bin. Fig. 8 also explicitly demonstrates that for the planets that maintain their atmospheres against the photoevaporation, the predicted planet radii are consistent with or larger than the those observed for the majority of the planets. At long orbital periods (20 – 100 days) all planets except one are consistent, within the observational uncertainties, with the predictions (the minimum accreted atmospheres in yellow, the maximum accreted atmospheres in light yellow, and the region in between). At intermediate periods (5 – 20 days), about a third of planets that are not bare cores have radii smaller than the predicted radii at the same mass. Finally, at short periods of less than 5 days, there are noticeably 5 planets with radii of $\sim$ $1.8-2$ R⊕ that are neither consistent with the mass-radius relationship of rocky cores, nor with the presence H/He envelopes. This suggests, potentially, that the cores of these planets could contain significant amounts of ice. Still, majority of the short-period planets are consistent with the predictions. Additionally, we reiterate that, while there might be exceptions, the radius distribution of sub-Neptunes is consistent with cores being largely rocky (Owen & Wu, 2017). Therefore, these results confirm our inference that typically our planets accrete too much gas. In section 4 we suggest possible explanations for this result, such as incorrect assumptions of quasi-hydrostatic accretion and negligible heating from planetesimal accretion, or missing mass-loss mechanisms that might act during, or after, disc dispersal.

Overall, the atmospheres accreted in the inner disc are typically in agreement with or larger than those observed, with the exception of planets with significant high-energy fluxes within a very narrow range. This is because core accretion is so efficient that considerable atmospheres can be accreted in the hot and low-density MRI-accreting inner disc and also maintained against photoevaporation.

4 Discussion and conclusions

In this paper we have investigated two aspects of planet formation in the inner disc that is viscously accreting due to the MRI. We considered the earliest phase of planet formation that is evolution of dust particles, and the final phase that is the shaping of the planetary atmospheres. We present simple calculations that include both the accretion and subsequent photoevaporation of close-in super-Earth/mini-Neptune planets. By coupling this processes to perform end-to-end calculations we are able to asses the viability of the in situ formation model for close-in planets.

In section 2, we consider the evolution of dust grains in the inner disc that is viscously accreting due to the MRI (Mohanty et al., 2018), that features a local gas pressure maximum at the orbital distance of a few tenths of AU. Taking into account the effect of the MRI-induced turbulence on the dust grain size, we find that fragmentation of particles due to turbulent relative velocities limits the particle size to below few millimeters. As a result, the particles are not efficiently accumulated inside the pressure maximum as hypothesised by Chatterjee & Tan (2014); Hu et al. (2018). Regardless of that, as the particles become well coupled to the gas, the radial drift is negated in the inner disc, and the dust-to-gas ratio is enhanced throughout the inner disc. Thus, the local gas pressure maximum might play a lesser role in the in situ planet formation than previously thought.

The pressure maximum is, however, still the location of a local density maximum in both gas and dust. We explored if the resulting inner disc structure that is enriched in dust could be susceptible to the onset of the streaming instabilities. This pathway to planetesimals seems to be viable only in a narrow region near the pressure (and density) maximum, for the chosen disc parameters.

The gas is not, however, evolved in this work and effects of the growing amounts of dust onto the MRI have not been taken into account. Dust grains lower the gas ionization levels by absorbing free charges and enhancing recombination rates, as ions recombine on the grains (Draine & Sutin, 1987; Ilgner & Nelson, 2006). Charged grains are not themselves well coupled to the magnetic field as they are too massive, and so their presence promotes the non-ideal MHD effects which can suppress the MRI (Sano et al., 2000; Ilgner & Nelson, 2006; Wardle, 2007; Salmeron & Wardle, 2008; Bai & Goodman, 2009). The likely result of taking dust effects into account is thus weakened turbulence, and the change of the disc structure in the longer term. The consequences can only be investigated by modelling both the gas and the dust self-consistently. Here, we sketch out a potential scenario by considering the relevant timescales.

Assuming that the steady-state solution of the gas structure (Fig. 1) is reached before dust starts affecting the MRI, we expect the dust enhancement of the inner disc to ensue. At a certain dust-to-gas ratio the dust will suppress the MRI, and we expect the levels of turbulence to adapt almost instantly, as the timescale of the magnetic field regeneration is the orbital timescale $t_{\rm orb}$ (e.g. Balbus & Hawley, 1991). With the decreasing levels of turbulence, the dust particle size will rapidly grow due to particle coagulation. The growth due to coagulation happens on the timescales of $\Sigma_{\rm g}/\Sigma_{\rm d}t_{\rm orb}$ (e.g. Brauer et al., 2008), so faster than $10^{2}t_{\rm orb}$ if the inner disc is indeed enriched in dust.

The gas disc structure would not change over such short timescales. But, in the absence of the viscous heating due to the MRI, the disc will cool at a timescale of $10^{2}$ – $10^{4}t_{\rm orb}$ (the equilibrium thermal timescale from roughly the pressure minimum to the pressure maximum in the model considered here; Mohanty et al., 2018). The pressure profile would also follow this timescale, as the vertical hydrostatic equilibrium is established quickly on the orbital timescale ($t_{\rm orb}$). This would likely result in the pressure maximum moving radially inwards.

Concurrently, due to larger particle size and lower turbulent stirring the particles would vertically settle towards the midplane and radially towards the pressure maximum, increasing both the surface density and the midplane bulk density dust-to-gas ratio there. Such formation of a ring of solids could potentially trigger formation of larger bodies, such as planetesimals (as hypothesized by e.g. Chatterjee & Tan, 2014). The larger particle size and the settling towards the midplane would likely trigger the streaming instability (inwards of the pressure maximum; see Fig. 4). However, it is unclear if this could lead to the formation of planetesimals, as a gravitational collapse is unlikely due to the low bulk dust densities and high Roche density in the inner disc. Moreover, as discussed above, the pressure maximum is expected to move radially inwards, and so will the accumulated dust, whereas the Roche density steeply increases inwards. Thus, over time the pressure maximum would need to accumulate significantly more dust to cross the gravitational collapse threshold.

The gas accretion rate and the gas surface density will change slowly in comparison to the above processes, on the long viscous timescale, $\sim$ $10^{3}$ – $10^{5}t_{\rm orb}$ from the pressure minimum to the pressure maximum in the model considered here (Mohanty et al., 2018). Suppression of the MRI would lead to increased amounts of gas in the inner disc on this timescale. However, if planetesimals are formed, this would clear the inner disc of the dust grains and the MRI could be induced again, decreasing the gas surface density and moving the pressure maximum outwards. At this stage it is unclear whether these competing processes are balanced in a steady state, or the behaviour of the inner disc is dynamic and quasi-periodic. Such a determination can only be investigated through self-consistent modelling of dust, gas and the MRI.

At high dust-to-gas ratios the dust also becomes dynamically important, and affects the gas disc structure through the drag backreaction (Nakagawa et al., 1986). The gas rotation profile is then driven towards Keplerian, and as a result the radial gas pressure profile flattens. This, in turn, slows down the radial drift of dust particles. If dust already pile ups in the inner disc due to radial drift being slower than in the outer disc, the dust backreaction amplifies the effect (Drążkowska et al., 2016). In this work, dust enhancement is driven by the dust grains already being completely coupled to the gas in the innermost disc, and thus the effect of dust backreaction would be limited. However, the backreaction would become important if the dust grains grow (e.g. due to the supression of the MRI-induced turbulence discussed above), especially near the pressure maximum. If the dust grains grow in the innermost disc where the pressure gradient is negative, the backreaction would slow down the loss of dust to the star. However, the backreaction would also limit the concentration of dust that can be achieved at the pressure maximum, since it acts to flatten the overall gas pressure profile (Taki et al., 2016).

Furthermore, if super-Earth and mini-Neptune cores indeed form in situ, would the inferred low gas surface densities due to the MRI allow them to acquire the observed $0.1-10$ % envelope mass fractions? We find that they would. In fact, even after accounting for atmospheric evaporation, the calculated atmospheres tend to overestimate the observed ones.

Could the atmospheric accretion in the MRI-implied disc and the observations be brought into agreement, without invoking an assumption that cores form just before the beginning of disc dispersal (e.g. Ikoma & Hori, 2012; Lee & Chiang, 2016)? The calculations shown here do not include several effects which could contribute.

First of all, for core masses smaller than $10$ M⊕, the discrepancy could be explained by the “boil-off” or core powered mass-loss (Owen & Wu, 2016; Ginzburg et al., 2018; see also Ikoma & Hori, 2012), a process in which a planet atmosphere that had not cooled and contracted before the disc dispersal loses its mass. Upon the dispersal the stellar continuum radiation illuminates the planet and launches a Parker wind. The mass loss causes rapid contraction of the atmosphere, and the contraction in turn shuts off the mass loss. Planets that start out with few 10s of percent atmospheres, may be left with 1 % after the boil-off. This process precedes the mass loss caused by the stellar high-energy flux considered above, and can operate at larger distances from the star.

Secondly, the scaling relations we use to calculate the accreted atmospheres are derived assuming no sources of heating due to planetesimal accretion, or due to heat deposited in the hypothesized final stage of giant mergers of planetary embryos. The latter could be released for several kyr (e.g. Inamdar & Schlichting, 2015), lowering the cooling rate of the atmosphere, and thus allowing less gas to be accreted. Furthermore, the scaling relations assume that the gas inside the planet’s Hill sphere is bound and static. Three-dimensional simulations suggest this may not be true and that high-entropy disc material is recycled between the envelope and disc (e.g. Ormel et al., 2015; Fung et al., 2015; Cimerman et al., 2017) potentially modifying the atmospheric cooling rate.

Finally, if the giant mergers happen between planets, after the disc has fully dispersed, they would likely result in significant atmospheric mass loss. Head-on collisions between Earth/super-Earth-sized planets with few-percent atmospheres can remove tens of percent of the total atmospheric mass (Liu et al., 2015; Inamdar & Schlichting, 2016).

Nevertheless, to avoid the runaway accretion for more massive cores, the low gas surface densities the MRI provides are favourable compared to the MMSN environment. Furthermore, gas-poor conditions in this case are provided in a long-lived state, and not in a transient phase (e.g. a transition disc, as proposed by Lee & Chiang, 2016).

In summary, our results support the hypothesis that the MRI-driven accretion in the inner protoplanetary disc could lead to in situ planet formation. However, there are several avenues that need to be explored in more detail until we can make quantitative predictions. In particular, the feedback between the enhancement of the dust in the inner disc and the suppression of the MRI, the feedback between the dust enhancement and the gas dynamics, and the role of disc dispersal and boil-off/core-powered mass-loss in shaping the final envelope masses of super-Earths/mini-Neptunes are issues that deserve closer study.

Acknowledgements

We thank the reviewer for helpful suggestions that improved the manuscript. We thank Jonathan Tan, Lauren Weiss and Eve Lee for helpful discussions. MRJ is funded by the President’s PhD scholarship of the Imperial College London and the "Dositeja" stipend from the Fund for Young Talents of the Serbian Ministry for Youth and Sport. JEO is supported by a Royal Society University Research Fellowship.