Limits on Protoplanet Growth by Accretion of Small Solids

TL;DR

This study models the growth limits of small planetary cores by pebble accretion, revealing that high accretion rates lead to turbulent, dusty envelopes that hinder further core growth due to fragmentation and opacity effects.

Contribution

It introduces detailed time-dependent envelope models including particle evolution and turbulence, providing new insights into pebble accretion constraints and envelope dynamics.

Findings

Pebble fragmentation dominates over dust re-sticking.

Envelope turbulence limits core growth rate to ~10^{-7} M_earth/yr.

Envelope metallicity and opacity increase with accretion rate.

Abstract

This paper identifies constraints on the growth of a small planetary core (0.3 M$_{\oplus}$) that accretes millimeter-sized pebbles from a gaseous disk. We construct time-dependent spherical envelope models that capture physical processes not included in existing global hydrodynamic simulations, including particle size evolution, dust transport, and realistic gas equations of state. We assume a low enough disk density that pebbles are marginally coupled to the gas and are trapped efficiently near the core Bondi radius. Pebbles then drift rapidly enough to experience erosion by sandblasting, mutual collisions, and sublimation of water ice. We find that pebble fragmentation is more efficient than dust re-sticking. Therefore the large pebble accretion rate $\dot{M}_p$ needed to build a core of mass > M$_{\oplus}$ leads to a high envelope metallicity and grain opacity. Above $\dot{M}_p$ $\sim 10^{-7}$ M$_{\oplus}$/yr, and without other luminosity sources, convective motions expand near the Bondi radius. The warm, dusty, and turbulent envelope buffers the inward drift of pebble debris: given a turbulent concentration factor $f_{turb} \geq$ 1 near the lower convective boundary, the core growth rate is limited to $1 \times 10^{-7} f_{turb}$ M$_{\oplus}$/yr and the e-folding time 3$f_{turb}^{-1}$ Myr. The remainder of the solid debris is expelled as highly processed silicates. Pebble ice never reaches the core, and the envelope contains comparable amounts of H2/He and metals. We interpret our results using simpler steady models and semi-analytical estimates. Future global simulations incorporating the processes modelled here are needed to understand the influence of rotation and vertical disk structure.

Figures

Tables

Table 1. Table 1: Gas Opacity Coefficients

For all T			T <800 K	T >800 K
c1	10.602	c8	-14.051	82.241
c2	2.882	c9	3.055	-55.456
c3	6.09$\times {10}^{-15}$	c10	0.024	8.754
c4	2.954	c11	1.877	0.7048
c5	-2.526	c12	-0.445	-0.0414
c6	0.843	c13	0.8321	0.8321
c7	-5.490	—	—	—

Table 2. Table 2: Adjustable Parameters of the Model

Parameter	value	note
${\dot{M}}_p$	${10}^{-7}$, ${10}^{-6}$ $M_{\oplus }$ yr-1	pebble accretion rate
${\dot{M}}_{\mathrm{plan}}$	0, ${10}^{-7}$ $M_{\oplus }$ yr-1	planetesimal accretion rate
$M_c$	$0.3M_{\oplus }$	initial protoplanet (core) mass
$R_c$	5200 km	core radius (inner boundary)
${\rho }_{g,\mathrm{disk}}$	10-11 g cm-3	disk gas density
$T_{\mathrm{disk}}$	100 K	disk temperature
$R_B$	$3.3\times {10}^{10}$ cm	Bondi radius ($\sim$ outer boundary)
$a_p$	0.1 cm	radius of injected pebbles
$a_{d,0}$	1 $\mu$m	dust radius at $t=0$
$X_{d,\mathrm{disk}}$	$0.01$	disk metallicity
$X_{\mathrm{ice}}=X_{\mathrm{sil}}$		mass fractions in injected dust/pebbles
$V_f$	${10}^2/{10}^3$ cm s-1	dust/ice fragmentation speed

Table 3. Table 3: Derived Parameters of the Model

${\dot{M}}_p$ [$M_{\oplus }$ yr-1]	M${}_{}{}^{\mathrm{sat}}{}_{H_2}{}^{}$ [$M_{\oplus }$]	M${}_{}{}^{\mathrm{sat}}{}_Z{}^{}$ [$M_{\oplus }$]	$E_{\mathrm{bind}}^{\mathrm{sat}}$ [erg]	${\dot{M}}_{\mathrm{sett},\mathrm{core}}/{\dot{M}}_p$	RRCB/RB
${10}^{-7}$	3.28$\times {10}^{-6}$	2.05$\times {10}^{-6}$	6.65$\times {10}^{33}$	1	0.45
$1.5\times {10}^{-7}$	3.58$\times {10}^{-6}$	2.54$\times {10}^{-6}$	7.63$\times {10}^{33}$	0.93	$>1$
$2\times {10}^{-7}$	3.72$\times {10}^{-6}$	2.75$\times {10}^{-6}$	8.07$\times {10}^{33}$	0.79	$>1$
$3\times {10}^{-7}$	3.67$\times {10}^{-6}$	2.68$\times {10}^{-6}$	7.92$\times {10}^{33}$	0.48	$>1$
${10}^{-6}$	3.48$\times {10}^{-6}$	2.42$\times {10}^{-6}$	7.36$\times {10}^{33}$	0.11	$>1$
${10}^{-5}$	3.41$\times {10}^{-6}$	2.32$\times {10}^{-6}$	7.14$\times {10}^{33}$	0.013	$>1$
${10}^{-7}$ (${\dot{M}}_{\mathrm{plan}}={10}^{-7}$)	2.72$\times {10}^{-6}$	1.44$\times {10}^{-6}$	5.19$\times {10}^{33}$	0.61	$>1$
${10}^{-6}$ (${\dot{M}}_{\mathrm{plan}}={10}^{-7}$)	2.42$\times {10}^{-6}$	1.09$\times {10}^{-6}$	4.37$\times {10}^{33}$	0.072	$>1$

Equations

$$R_{\mathrm{B}}$$

$$\text{mbox}\frac{R_{\mathrm{B}}}{R_{\mathrm{H}}}$$

$${\rho }_{g,\mathrm{disk}}\lesssim 1\times 10^{-11}\frac{T_{\mathrm{disk},2}{a}_{p,-1}{\rho }_{s,0}}{{\epsilon }_{V,-1}}{\left(\frac{M_c}{0.3M_{\oplus }}\right)}^{-1}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}.$$

$${\dot{M}}_p\sim 4\pi r^2{\tau }_{\mathrm{stop}}\frac{G{M}_c}{r^2}\cdot n_p\frac{4\pi }{3}{\rho }_s{a}_p^3,$$

$${\tau }_{pp}(r\lesssim R_{\mathrm{B}})\sim n_{p}4\pi a_p^{2}r\sim 2\left(\frac{{\rho }_g}{{\rho }_{g,\mathrm{disk}}}\right)\frac{{\rho }_{g,-11}}{T_2^{1/2}(a_{p,-1}{\rho }_{s,0}{)}^2}\left(\frac{{\dot{M}}_p}{10^{-6}{M}_{\oplus }{\mathrm{yr}}^{-1}}\right)\left(\frac{r}{R_{\mathrm{B}}}\right).$$

$$a_{d,\mathrm{frag}}\sim \frac{9V_{f}G{M}_c{\rho }_{g,\mathrm{B}}}{4{\rho }_s{c}_{g,\mathrm{B}}}=4\frac{V_{f,2}{\rho }_{g,\mathrm{B},-11}}{T_{\mathrm{B},2}^{3/2}}\left(\frac{M_c}{0.3M_{\oplus }}\right)\mu \mathrm{m}.$$

$$\frac{V_{\mathrm{con}}}{c_g}\sim {\left(\frac{G{M}_c^2}{R_c{t}_{\mathrm{acc},c}}\right)}^{1/3}{\left(4\pi R_{\mathrm{B}}^2{\rho }_g{c}_g^3\right)}^{-1/3}\sim 0.3\frac{T_2^{1/6}}{{\rho }_{g,-11}^{1/3}{R}_{c,9}^{1/3}}{\left(\frac{t_{\mathrm{acc},c}}{10\mathrm{Myr}}\right)}^{-1/3}$$

$$F_{\mathrm{D}}=-\frac{4\pi }{3}{\rho }_g{a}^2{v}_{\mathrm{th}}\Delta V(\mathrm{Epstein}),$$

$$F_{\mathrm{D}}=-6\pi \eta a\Delta V(\mathrm{Stokes}),$$

$${\tau }_{\mathrm{stop}}=\frac{{\rho }_{s}a}{{\rho }_g{c}_g}\times \max \left(\frac{4}{9},\frac{a}{{\ell }_g}\right).$$

$$\frac{d{a}_p}{dt}=-0.63{\left(\frac{{\mu }_{{\mathrm{H}}_2\mathrm{O}}}{2\pi kT}\right)}^{1/2}\frac{P_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{sat}}}{{\rho }_s}.$$

$$\frac{d}{dt}\left(\frac{4\pi }{3}{\rho }_s{a}_p^3\right)=-Y(\Delta V)\times \pi a_p^2{\rho }_d\Delta V.$$

$$Y(\Delta V)=0.7{\left(\frac{\Delta V}{V_f}\right)}^{\alpha },$$

$${\dot{n}}_d=n_d{n}_p\cdot Y\pi a_p^2\Delta V(a_p)-n_d^2\cdot 4\pi a_d^2\Delta V(a_d).$$

$${\dot{n}}_d=\frac{3Y{n}_d{\dot{M}}_p}{16\pi r^2{\rho }_s{a}_p}\left(1-\frac{X_d}{X_{d,\mathrm{col}}}\right),$$

$$X_{d,\mathrm{col}}=\frac{9Y{\dot{M}}_p{c}_g}{64\pi G{M}_c{\rho }_s{a}_p}=0.47\frac{T_2^{1/2}}{a_{p,-1}}\left(\frac{Y}{10}\right){\left(\frac{t_{\mathrm{acc},p}}{0.3\mathrm{Myr}}\right)}^{-1}$$

$${\tau }_{\mathrm{stop}}(a_d)\sim \frac{4\pi V_f^2{r}^3{\rho }_g}{L_{\mathrm{con}}}.$$

$${\tau }_{\mathrm{stop}}(a_d)\sim \frac{V_f{r}^2}{GM(r)}.$$

$${\tau }_{\mathrm{col}}=\frac{1}{n(a_d)4\pi a_d^2{V}_t(a_d)}<\frac{r}{V_{\mathrm{con}}},$$

$$\frac{a_d}{a_d+4{\ell }_g/9}<4X_d^2\frac{{\rho }_g}{{\rho }_s}{\mathcal{M}}_{\mathrm{con}}\frac{r}{{\ell }_g}.$$

$$\frac{{\tau }_{\mathrm{col}}{V}_{\mathrm{con}}}{r}\sim \frac{3}{4X_d}\frac{c_g^{2}r}{GM(r)}\frac{{\mathcal{M}}_{\mathrm{con}}}{1+9a_d/4{\ell }_g}<1$$

$$\frac{1}{a}\frac{da}{dt}=-f_v\frac{D_{\mathrm{sil}}{\mu }_{\mathrm{sil}}}{{\rho }_s{a}^2{k}_{\mathrm{B}}T}(P_{\mathrm{sil}}^{\mathrm{sat}}-P_{\mathrm{sil}}).$$

$$D_{\mathrm{sil}}\simeq 0.6\frac{(kT/{\mu }_{\mathrm{H}}{)}^{1/2}}{(a_{{\mathrm{H}}_2}+a_{\mathrm{sil}}{)}^2{n}_{\mathrm{H}}}.$$

$$\frac{1}{a}\frac{da}{dt}\simeq -3\times 10^{-2}{f}_v{\left(\frac{a}{\mathrm{mm}}\right)}^{-2}{\left(\frac{T}{2000\mathrm{K}}\right)}^{1/2}\frac{X_{\mathrm{sil}}}{X_{\mathrm{H}}}\frac{P_{\mathrm{sil}}^{\mathrm{sat}}-P_{\mathrm{sil}}}{P_{\mathrm{sil}}}{\mathrm{s}}^{-1}.$$

$$\frac{d{P}_g}{dr}$$

$$\text{mbox}\frac{L_{\mathrm{rad}}}{4\pi r^2}$$

$$\frac{d({\tilde{\rho }}_g{\tilde{T}}_g)}{d\tilde{r}}$$

$$\text{mbox}\frac{{\tilde{T}}_g^{2-\beta }}{{\tilde{\rho }}_g}\frac{d{\tilde{T}}_g}{d\tilde{r}}$$

$${\epsilon }_{\mathrm{rad}}=\frac{3{\tau }_{\mathrm{B}}}{64\pi }\frac{L_{\mathrm{rad}}{c}_{g,\mathrm{disk}}^4}{(G{M}_c{)}^2{\sigma }_{\mathrm{SB}}{T}_{\mathrm{disk}}^4}=\frac{3{\tau }_{\mathrm{B}}}{64\pi }\frac{c_{g,\mathrm{disk}}^4}{G{R}_c{t}_{\mathrm{acc},c}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{disk}}^4}$$

$$L_{\mathrm{rad}}=\frac{G{M}_c^2}{R_c{t}_{\mathrm{acc},c}}.$$

$${\tau }_{\mathrm{B}}\equiv \frac{X_d}{1-X_d}{\kappa }_d(T_{\mathrm{disk}}){\rho }_{g,\mathrm{disk}}\frac{G{M}_c}{c_{g,\mathrm{disk}}^2}=0.77{\left(\frac{X_d}{1-X_d}\right)}_{-2}{\rho }_{g,\mathrm{disk},-11}{T}_{\mathrm{disk},2}^{\beta }{\left(\frac{M_c}{0.3M_{\oplus }}\right)}^{2/3}{\left(\frac{t_{\mathrm{acc},c}}{\mathrm{Myr}}\right)}^{-1},$$

$${\epsilon }_{\mathrm{rad}}=0.024{\left(\frac{X_d}{1-X_d}\right)}_{-2}\frac{{\rho }_{g,\mathrm{disk},-11}}{T_{\mathrm{disk},2}^{2-\beta }}{\left(\frac{M_c}{0.3M_{\oplus }}\right)}^{2/3}{\left(\frac{t_{\mathrm{acc},c}}{\mathrm{Myr}}\right)}^{-1}.$$

Full text

Astrophysical dust processes (99), Planetary atmospheres (1244), Radiative transfer (1335), Planet formation (1241)

1 Introduction

The accretion of small, collisional “pebbles” has been suggested as an efficient channel for planet growth, in the first instance in order to explain the assembly of Uranus and Neptune in our solar system. Pebble accretion may be enhanced either by collisional cooling in a thin, pebble-dominated disk (Goldreich et al., 2004a) or by gas drag (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012). Here we focus on pebble accretion from a gaseous protoplanetary disk (PPD), examining how collisional effects may limit core growth from pebble debris.

A key constraint on this accretion channel is that the mm-cm sized pebbles must only be mildly coupled to the gas near the planetary Bondi radius, as defined by the temperature of the ambient gas disk: if the coupling is tight then pebble inflow is limited by the residency of gas outside the Bondi radius. Since the flow speed near the Bondi radius is typically a modest fraction of the sound speed (e.g. Popovas et al. 2018), accreting pebbles must drift at relatively high speeds ($\sim 100$ m s*-1*) with respect to the gas. Pebbles formed by the adhesion of micron-sized grains through van der Waals forces are fragile, and binary collisions at such speeds are destructive (Beitz et al., 2011; Schräpler et al., 2012; Hill et al., 2015). These collisions occur with significant optical depth if the pebble accretion rate is high enough to grow a core above an Earth mass over an interval of 1 Myr.

Even faster abrasion of pebbles occurs, as we show, through collisions with dust grains embedded in the gas (see Jacquet & Thompson 2014 in the context of solid-gas shocks). Small target grains have a much longer residency time in the envelope than do the pebbles, meaning that their abundance is not fixed and can substantially exceed the abundance in the ambient disk. Pebble abrasion combined with mutual fragmenting collisions between the small particles causes an exponential rise in the envelope dust abundance.

Icy pebbles are also destroyed by rapid sublimation as the temperature rises above 100 K. Pebble destruction inhibits the direct incorporation of water or silicates into the core, but does allow enrichment in dust and the extraction of heat from the gaseous component through ice sublimation (Stevenson, 1984; Hori & Ikoma, 2011; Chambers, 2017). The net result is that gas around the Bondi radius can become sufficiently dusty to suppress radiative transport, well outside the inner hydrostatic gaseous envelope where pebble destruction has previously been considered (Pollack et al., 1986; Podolak, Pollack, & Reynolds, 1988; Venturini, Alibert, & Benz, 2016; Brouwers et al., 2018).

The retention of gas and solids depends on efficient radiative cooling (e.g. Lee & Chiang 2015). When cooling is weak, Ormel, Shi, & Kuiper (2015) and Alibert (2017) showed using three-dimensional isothermal simulations that the accreted envelope surrounding a core behaves like an open system, with the envelope material continually being recycled back into the gas disk. As a result, the accretion of pebbles from the neighboring disk into a bound, quasi-spherical envelope (concentrated inside the Bondi radius of the planetary core) may be intrinsically inefficient (Kuwahara et al., 2019). Lambrechts & Lega (2017) included heating by planetesimal accretion onto the core surface and energy exchange between the gas and diffusing radiation. Kurokawa & Tanigawa (2018) also considered the case where gas accreted close to the core is able to cool radiatively, but ignored the potentially strong heating effects associated with the accretion of a modest number of large planetesimals. The refluxing of envelope material away from higher-mass cores, driven by planetesimal accretion, has been considered by Rafikov (2006) in an analytic spherical envelope model with parameterized opacity, but the effects of dust production and ice sublimation were not included.

Here we demonstrate (i) an even stronger suppression of direct pebble accretion onto the core than has previously been obtained; and (ii) a limitation on envelope growth due to the reflux of gas and dust by convective motions back into the disk. A key result is that fragmentation can become a runaway process, so that the density of small particles settling toward the core is strongly enhanced as their size and inward drift speed is reduced. Previous models of dust deposition in the envelope typically include growth by adhesion but not collisional fragmentation, and imply inward drift speeds exceeding $\sim 10^{2}$ m s*-1* (Mordasini, 2014; Ormel, 2014).

We find that the limitation on core accretion is strongest for the water ice carried by pebbles. In addition, the H2/He mass accumulated in the atmosphere is strongly tied to the onset of convection near the Bondi radius. We develop a time-dependent model of the structure of the gas envelope accreted from a PPD that includes the feedback of pebble destruction on radiative opacity and the possibility of convection near and beyond the Bondi radius. The fiducial core mass is $0.3\,M_{\oplus}$, which is significantly heavier than the bound planetesimals which have been found to collapse directly via the streaming instability in pebble-rich gas (Johansen & Lambrechts, 2017). Cores of around this fiducial mass were found in hydrodynamic simulations to accrete pebbles vigorously (Lambrechts & Johansen, 2012, 2014).

The capture rate of pebbles by a planetary core could cover a range of values, and a goal of this study is to determine the critical rate above which the envelope is heated rapidly enough, and its opacity is high enough, that much of the pebble debris is returned to the disk. The limiting accretion rate $\dot{M}_{p}$ so obtained can be compared with the rate needed to grow the core above a mass $M_{c}\sim(1-10)M_{\oplus}$ in the lifetime of a PPD. In a variant on our default luminosity model, we also consider an additional source of heating by the accretion of large planetesimals that easily penetrate the convective envelope; this can compete with (or dominate) the gravitational heating by infalling pebble debris.

An important consideration is the rate at which pebble debris circulating in the lower envelope will be incorporated into the core. Popovas et al. (2018, 2019) have calculated pebble trajectories in a three-dimensional simulation of a core embedded in a shearing gas disk, finding rapid downward advection of pebbles by convective flows. Pebbles reaching the computational core boundary are allowed to be accreted, releasing gravitational energy and driving convection around the core. Although we also find that turbulent diffusion of grains through the convective layer is faster than secular drift in the central gravitational field, particles drifting at the silicate fragmentation speed $\sim 1$ m s*-1* are very tightly coupled to the gas. In fact, the simulations of Popovas et al. (2018, 2019) do not resolve the dynamics of small particles with short stopping times near the core boundary. We evaluate the core growth rate in terms of the drift rate of embedded grains across the lower convective-radiative boundary, which is proportional to the dust concentration. The turbulent intensity decreases toward the convective boundary, meaning that this concentration may be enhanced with respect to the interior of the convection zone, by a factor $f_{\rm turb}$. The core accretion rate is evaluated for $f_{\rm turb}=1$, but the possibility of an enhancement due to turbulent pumping should be kept in mind.

The plan of this paper is as follows. Our treatment of the interactions of solid particles with gas, and with each other, is described in Section 2. Section 3 presents a simple analytical steady-state model showing how convection expands within a hydrostatic envelope that is progressively loaded with small particles. Section 4 describes our one-dimensional (spherical) hydrostatic numerical model of the gas envelope, including heat exchange between sublimating solids and gas and the convective transport of embedded particles calculated in the mixing-length approximation. Results are presented in Section 5 for various pebble accretion rates in the $\dot{M}_{p}=10^{-7}-10^{-6}\,M_{\oplus}$ yr*-1* range (corresponding to a growth times $\sim 10^{6}-10^{7}$ yr to reach $10\,M_{\oplus}$ starting from $0.3\,M_{\oplus}$), and accretion rates $\dot{M}_{\rm plan}=0$, $10^{-7}\,M_{\oplus}$ yr*-1* of larger planetesimals. The effect of pebble composition (icy vs. dry) is also considered. The implications of our work and some outstanding issues are summarized in Section 6. The Appendix details a simple analytic model which is consistent with the detailed numerical results.

2 Interactions of Solid Particles

Destructive processes acting on pebbles are the main subject of this section. We begin by summarizing why fast pebble drift is characteristic of efficient pebble accretion. A high pebble accretion rate (needed to build a substantial core over a Myr interval) is connected to a significant optical depth to destructive binary collisions between pebbles. We explain the relative importance of sandblasting by small grains embedded in the gas, and show that dust production is not limited by resticking of grains before the dust/gas ratio reaches of the order unity.

Next we quantify our treatment of pebble destruction in the numerical model, defining a self-consistent one-size model for pebbles embedded in a turbulent, convective atmosphere. Although our default envelope model does not include the effect of silicate dust sublimation on core accretion, we present additional arguments supporting the suppression of silicate ‘rainout’ in the inner envelope, as a consequence of the large atmospheric scale height and silicate mass fraction.

In what follows, we use ‘pebble’ to denote 0.1-1 cm sized objects, and ‘planetesimal’ to denote 1-10 km sized objects. We focus on the regime $\dot{M}_{p}\gg\dot{M}_{\rm plan}$: pebbles are the main potential source of planetary material, but because the pebbles disintegrate in the outer envelope, planetesimals can supplement the heat flux through the envelope.

2.1 Some Introductory Estimates

The fiducial core mass is111This is also comparable to the planetesimal isolation mass near 5 AU in the minimum-mass solar nebula of Hayashi (1981) (Goldreich et al., 2004b). $0.3\,M_{\oplus}$, which is light enough that gas pressure gradient forces play a significant role in determining the gas profile within the core’s tidal radius: in other words, the Bondi radius is smaller than the Hill radius222Throughout this paper, we use the shorthand $X=X_{n}\times 10^{n}$, where quantity $X$ is presented in c.g.s. units. For example $\rho_{s,0}$ is the pebble density in units of $1=10^{0}$ g cm*-3*. ($R_{\rm H}=a_{\rm orb}(M_{c}/3M_{\odot})^{1/3}$),

[TABLE]

where $a_{\rm orb}$ is the orbital semi-major axis around a Solar-mass star and $c_{g,\rm disk}=(k_{\rm B}T_{\rm disk}/\mu_{g})^{1/2}$ is the isothermal gas sound speed in the ambient disk of temperature $T_{\rm disk}$ and mean molecular weight $\mu_{g}$. (In what follows, “disk” refers to the medium surrounding the core.) Such a core is just massive enough to sublimate silicate grains embedded in a high-entropy H2/He atmosphere, but its gravity is too weak for the impacts of marginally-bound large planetesimals to convert significant silicate mass to vapor.

Efficient accretion of pebbles inside the Bondi radius implies an upper bound on the PPD gas density $\rho_{g,\rm disk}$. Writing the gas flow speed around the core outside the Bondi radius as $V_{g}=\varepsilon_{V}c_{g,\rm disk}$, we require that the pebble drift speed $\Delta V$ with respect to the gas is $\Delta V\gtrsim V_{g}$ for the pebble to be trapped. At some radius around $R_{\rm B}$, the gas flow bifurcates from inflow-outflow to a hydrostatic, rotational motions (e.g. Kurokawa & Tanigawa 2018). Near this bifurcation point, $\Delta V\sim\tau_{\rm stop}GM_{c}/R_{\rm B}^{2}\sim\tau_{\rm stop}c_{g}^{4}/GM_{c}$. The radius $a_{p}$ of a trapped pebble is typically smaller than the gas mean-free path, and so the stopping time $\tau_{\rm stop}\sim\rho_{s}a_{p}/\rho_{g,\rm disk}c_{g,\rm disk}$, where $\rho_{s}$ is the density of the pebble material. Then the condition $\Delta V\gtrsim V_{g}$ implies

[TABLE]

This density is typical of the zone outside the water ice line (semi-major axis $\gtrsim 3$ AU) when about one Jupiter mass of gas remains in the PPD. For an orbit typical of an extra-Solar super-Earth (0.1 AU), it implies a low gas surface density $\sim 1$ g cm*-2* (corresponding to inefficient planetary migration by planet-disk torques).

Our model is schematically illustrated in Figure 1. Pebbles experience destructive collisions passing through dust-rich gas near the outer boundary of the bound envelope. When solid particles are subject to a strong, transient acceleration, they drift differentially through the gas with respect to each other, and especially with respect to small grains stuck in the gas. High-speed collisions with dust grains embedded in the gas effectively sandblast away the pebble material (see Schräpler & Blum 2011; Schräpler et al. 2018 for experimental constraints and the application to PPDs and Jacquet & Thompson 2014 for solid-gas shocks). This destruction channel, which builds a strongly bimodal size distribution, is more effective than collisions between similarly-sized particles, which are the dominant process in a mildly turbulent PPD (e.g. Krijt et al. 2015). For a typical drift speed $10^{2}$ m s*-1*, a grain is expected to liberate $Y=10-100$ times the grain mass from the pebble surface (Jacquet & Thompson, 2014; Schräpler et al., 2018).

When the seed dust density is small, binary collisions between pebbles provide an additional population of small targets. Around the Bondi radius, the pebble accretion rate

[TABLE]

where $n_{p}$ is the space density of pebbles. We can write the accretion rate in terms of an accretion time, $\dot{M}_{p}=M_{c}/t_{{\rm acc},p}$. Then the optical depth for binary collisions within an outer radiative layer of the envelope ($\rho_{g}>\rho_{g,\rm disk}$, $T\simeq T_{\rm disk}$) is

[TABLE]

Here we have restored a numerical factor in the gas drag law (Equation (9)). This expression can also be expressed in terms of the gas flow speed $\varepsilon_{V}$ through Equation (2). It applies in a situation preceding the build-up of small grains by sandblasting and ice sublimation, in which case the outer envelope is nearly isothermal but the gas density increases exponentially inward (Section 4.1). Binary pebble collisions are not included in our time-dependent envelope model due to the relative effectiveness of ice sublimation and sandblasting.

Turning now to pebble-dust interactions, each infalling pebble intersects a column $X_{d}\rho_{g}GM_{c}/c_{g}^{2}\sim X_{d}\rho_{g,-11}T_{2}^{-1}(M_{c}/M_{\oplus})$ g cm*-2* of grains of mass fraction $X_{d}$. Small grains have a long residency time in the gas, and their density exponentiates within the bound envelope in response to the pebble ablation. We find that the mass fraction of particles smaller than the wavelength peak of the Planck function ($\lambda_{\rm max}(T_{\rm disk})\sim 30\,\mu$m; Equation (70)) must exceed $X_{d}\sim 0.3$ before the outer envelope switches to a convective state (Section 3.1; Figure 12). As a result, the column of intercepted grains is high enough that, combined with a large mass-loss multiplicity, each mm-cm sized pebble is destroyed already near the outer boundary of the envelope. Our numerical model, described in Section 4, self-consistently calculates the grain density in response to the sourcing from grain-pebble collisions, and the loss by advection out beyond the Bondi radius.

Given the relatively high abundance of the target particles, one must also consider interactions between them. The particle size as limited by binary fragmenting collisions is quite small. Setting the radial drift speed through the static gas equal to the fragmentation speed, $\Delta V(a_{d})=\tau_{\rm stop}GM_{c}/r^{2}=V_{f}$, one finds near the Bondi radius,

[TABLE]

Secondary collisions between similarly-sized collision fragments larger than $a_{d,\rm frag}$ will also be destructive (e.g. Güttler et al. 2010), and have a significant optical depth when the inward mass flow carried by the fragments rises to a fraction $\sim(a_{d}/a_{p})^{2}$ of $\dot{M}_{p}$ (consider Equation (4) with $a_{p}$ replaced by $a_{d}\ll a_{p}$). This means that the particle fragments with long residency times that are liberated by pebble-grain collisions will, as their abundance rises, break down to a near-monomer size. The pebble debris is therefore an effective source of radiative opacity. The interactions of even smaller (e.g. sub-micron sized) grains are considered in Section 2.3.

Eventually, the evolution of grains in the inner bound gas envelope is determined by convectively forced collisions. At high pebble accretion rates, the convective energy flux is dominated by a small fraction of the pebble material that settles onto the core, with a corresponding accretion time $t_{{\rm acc},c}=M_{c}/\dot{M}_{\rm sett}(R_{c})$ and luminosity $GM_{c}^{2}/R_{c}t_{{\rm acc},c}$. Then the convective Mach number reaches reaches

[TABLE]

at the Bondi radius. The turbulent acceleration of particles within these dynamic gas flows is high enough to prevent fragmented particles from re-sticking, thereby maintaining a dense population of micron-sized grains.

2.2 Pebble Dynamics and Destruction

Solid particles in the envelope are divided into two components: (i) large pebbles (of radius $a_{p}\sim$ mm), which are injected by accretion inside the Bondi radius; and (ii) dust grains, which are produced mainly by ‘sandblasting’ of pebbles during collisions with other grains embedded in the gas, and whose size deeper in the envelope is regulated by binary collisions.

Pebbles are less tightly coupled to the gas than the smaller grains, and so drift more rapidly with respect to the local gas flow. Nonetheless, mm-sized pebbles are at least marginally coupled to the convective motions in the outer envelope, once these motions develop in response to the growing radiative opacity. The mean free path between gas molecules is $\ell_{g}\equiv\mu_{\rm H}/\sigma_{\rm H}\rho_{g}\sim 10^{2}\rho_{g,-11}^{-1}$ cm in the outer envelope. Here $\sigma_{\rm H}\sim 3\times 10^{-15}$ cm*-2* is the mutual collision cross section of H2 molecules. As a result, $\ell_{g}>a$ for both grains and pebbles, and the drag force is

[TABLE]

where $v_{\rm th}\sim 1.6c_{g}=1.6(kT_{g}/\mu_{g})^{1/2}$ is the mean thermal speed of a gas molecule. In the inner envelope, fragmenting grains can have radii larger than $\ell_{g}$, whereas $\Delta V$ is much smaller than the gas sound speed. Then the drag is in the Stokes regime,

[TABLE]

where $\eta\simeq(1/3)\rho_{g}v_{\rm th}\ell_{g}$ is the viscosity. For the stopping time $\tau_{\rm stop}=(4\pi/3)\rho_{s}a^{3}\Delta V/F_{\rm D}$ we use the interpolation formula

[TABLE]

We take into account the erosion of pebbles by two effects: sublimation and sandblasting. Pebbles will sublimate as they drift inward into the hotter parts of the envelope. Where the water vapor pressure is below the saturation value, icy pebbles shrink at the rate (Podolak, Pollack, & Reynolds, 1988)

[TABLE]

Here, the numerical factor matches the measured sublimation rate at low pressure (Haynes et al., 1992). In addition, as pebbles drift into the envelope, they will be sandblasted by small dust particles coupled to the gas. This will result in the mass loss rate:

[TABLE]

In the above $Y(\Delta V)$ is a yield factor controlling the pebble mass ejected per collision

[TABLE]

where $V_{f}$ is the fragmentation velocity, and the coefficient and the index $\alpha\sim 1$ are determined by laboratory experiments (Schräpler et al., 2018). In our atmosphere models $Y=O(10-30)$ in the outer pebble destruction zone. In practice, we find that sandblasting is the dominant destruction channel.

2.3 Limitation on Very Small Grains in Outer Radiative Layer

The accumulation of very small dust grains (of a radius smaller than Equation (5)) may eventually be limited by mutual sticking. Here we focus on an outer radiative layer of the envelope. Sticking becomes important when the dust mass fraction $X_{d}$ reaches a level which we now estimate. This is proportional to the pebble accretion rate $\dot{M}_{p}$ and also to the multiplication factor $Y$ appearing in the mass loss formula (11). Given that the grains are effectively stuck in the gas, the net rate of change of the grain volume density is

[TABLE]

Here both grains and pebbles are in the Epstein drag regime. Relating the pebble volume density to the accretion rate via Equation (3), we obtain

[TABLE]

where for pebbles of icy composition,

[TABLE]

and the pebble accretion time $t_{{\rm acc},p}\equiv M_{c}/\dot{M}_{p}$. This dust mass fraction is comparable to the equilibrium metallicities observed in our time-evolved envelope models, for a pebble accretion rate $\dot{M}_{p}=10^{-6}\,M_{\oplus}$ yr*-1*.

2.4 Equilibrium Dust Size in the Convective Atmosphere

Interior to the pebble erosion layer, we represent the size of the embedded dust particles as being locally peaked at a value $a_{d}$ that is regulated by fragmentation and the rate of mutual collisions. At high accretion luminosities, dust grains in the upper envelope are stirred rapidly by turbulence, and broken down to sizes smaller than (5). In the inner envelope, the collision speed between dust particles is determined by the direct gravitational acceleration.

First consider collisions between grains that are exposed to convectively driven turbulence. Neglecting the effects of rotation, the eddy speed at a scale $\ell\ll r$ is $V_{t}\sim V_{\rm con}(\ell/r)^{1/3}$. Small particles with stopping time $\tau_{\rm stop}$ decouple from the convective eddies at a scale where $\ell/V_{t}\sim\tau_{\rm stop}$, and collide with each other at a speed $\Delta V\sim V_{t}\sim V_{\rm con}(V_{\rm con}\tau_{\rm stop}/r)^{1/2}$. This collision speed decreases as the particles lose mass, until the threshold for fragmentation is reached. Equating $\Delta V$ with the fragmentation speed $V_{f}$ gives $\tau_{\rm stop}V_{\rm con}/r\sim(V_{f}/V_{\rm con})^{2}$. The fragmentation velocity depends on chemical composition, size, porosity, and other factors. Here we simply use 1 m s*-1* and 10 m s*-1*, respectively, for dust and ice aggregates (Blum & Wurm, 2008; Wada et al., 2008, 2009). The size-dependent stopping time is related to $V_{f}$ and $L_{\rm con}$ via

[TABLE]

When the direct gravitational acceleration dominates turbulent stirring, the particle drift speed is instead $\Delta V=\tau_{\rm stop}GM(r)/r^{2}$; hence at the threshold for fragmentation,

[TABLE]

The maximum particle radius $a_{f}(r)$ is obtained by taking the minimum of the right-hand sides of Equations (16) and (17). Typically we find that turbulent stirring dominates near the Bondi radius; hence $a_{f}$ increases inward until Equation (17) dominates.

The convective motions transport embedded particles inward from the Bondi radius toward the core. After an initial phase of collisional breakdown in the outermost envelope, the threshold size for fragmentation increases inward (see Section 3.2 for examples). Inflowing particles reach this threshold only if two conditions are satisfied: first, the mean free time for collisions between similarly sized particles $\tau_{\rm col}$ must remain smaller than the convective time; and, second, the particles must not be so large and compact that that they bounce rather than stick. The second condition is a complicated one, in that it depends on the porosity of the particles (Blum & Wurm, 2008; Okuzumi et al., 2012). We therefore only apply the first constraint,

[TABLE]

where $n(a_{d})=3X_{d}\rho_{g}/4\rho_{s}a_{d}^{3}$ is the particle number density and we approximate $X_{d}\equiv\bar{\rho}_{d}/\rho_{g}$. In the case where convective stirring dominates, one requires

[TABLE]

Alternatively, when the particle drift is dominated by the central gravity, then one requires

[TABLE]

for the particles to reach the local fragmentation threshold.

To determine the local value of the particle size, we first determine $a_{f}$ from Equations (16) and (17). We set $a_{d}=a_{f}$ if the corresponding collision time is short enough for the local fragmentation bound to be reached, as determined by either Equation (19) or (20). Otherwise, if collisional equilibrium is lost at some radius, then we freeze the particle size interior to this radius. This overall method is analogous to the model that Birnstiel, Klahr, & Ercolano (2012) developed for PPDs.

2.5 Suppression of Silicate Rain

Our default model sets aside the effects of silicate sublimation in the inner gasesous envelope. This represents a conservative choice of core-envelope boundary condition for embedded grains. In Section 5.3, we consider the alternative situation where rain-out of silicate particles is suppressed by re-evaporation. This is motivated by the following argument that these particles (of partially solid and liquid composition) will efficiently sublimate during their descent to the core, as a result of the large partial pressure of silicate vapor and large atmospheric scale height.

The silicate saturation vapor pressure increases rapidly inward, so that if the core mass is high enough, a falling rain particle must encounter an inner layer of the envelope where the vapor pressure is well below saturation. Within this layer, particles whose size is limited by fragmentation ($a\sim 0.1$ mm) are in the Stokes drag regime. Then we must take into account the diffusion of silicate vapor molecules away from a sublimating particle, and Equation (10) must be replaced by (Pruppacher & Rasmussen, 1979)

[TABLE]

Here $f_{v}>1$ is the venting factor and the diffusion coefficient of silicate molecules through the lighter gas is (Lifshitz & Pitaevskii, 1981)

[TABLE]

In this expression, $n_{\rm H}\equiv X_{\rm H}\rho_{g}/\mu_{\rm H}$ is the number density of ${\rm H}_{2}$ and He molecules and $a_{\rm H_{2}}\sim a_{\rm sil}\sim 10^{-8}$ cm are the effective collisional radii of ${\rm H}_{2}$, SiO and ${\rm O}_{2}$. Then one finds

[TABLE]

The particle would take $\sim 10^{6}$ seconds to fall at the fragmentation speed ($\sim 1$ m/s) through the inner $10^{3}$ km of the envelope, meaning that it must sublimate well before reaching the core-envelope boundary. The relative efficiency of sublimation in this context, compared with the incidence of ‘dry rain’ in the Earth’s atmosphere, can be ascribed to a combination of two factors: the relatively high silicate vapor mass fraction (typically $X_{\rm sil}\sim(0.3-1)X_{\rm H}$ when refluxing of pebble debris from the envelope begins to balance pebble accretion) and an infall time time a factor $\sim 10^{3}$ times larger than for raindrops. Of course, much larger planetesimals that survive transport through the envelope may still inject solid/liquid silicates into the core.

3 Envelope Expansion in Response to

Planetesimal Heating and Dust Loading: Semi-Analytic Treatment

As a prelude to the numerical investigation described in Sections 4 and 5, we now present a semi-analytic model of an envelope with a low dust density and outer radiative structure. This model connects smoothly to a uniform medium of a fixed temperature $T_{\rm disk}$ and gas density $\rho_{g,\rm disk}$. As the dust abundance rises, the inner convective layer grows in size, and eventually overwhelms the Bondi radius. Assuming that the dust abundance in the outer radiative layer is tied directly to the solid accretion rate, we can work out the critical value of $\dot{M}_{p}$ above which convection extends beyond the Bondi radius, some refluxing of solids will occur, and the self-consistency of the model breaks down. We next quantify the reduction in dust size due to convective stirring, and demonstrate how a convective layer that traps accreted solids may expand in response to deep heating by large planesimals.

3.1 Two-Layer Convective-Radiative Envelope

Even a small amount of embedded dust acts as a powerful coolant, in the sense that any excess of gas temperature over ambient radiation temperature is rapidly converted to thermal radiation via collisions between gas molecules and grains (e.g., Chiang & Goldreich 1997). The presence of an ambient blackbody radiation field of temperature $T_{\rm disk}$ is therefore assumed, with an internal energy amounting to a small fraction $\sim 10^{-5}T_{\rm disk,2}^{3}\rho_{g,\rm disk,-11}^{-1}$ of the ambient gas energy density. The temperature is buffered near $T_{\rm disk}$ in the outer bound envelope, around the Bondi radius, and grows in the inner, optically thick envelope.

In a state of spherical hydrostatic and radiaive equilibrium at local temperature $T_{g}$ and uniform luminosity $L_{\rm rad}$,

[TABLE]

The opacity $\kappa_{d}(T)\propto T^{1+\beta}$ is per unit mass of dust. We first consider the case of uniform dust mass fraction $X_{d}$ (although the equations presented here do not depend on this assumption). We rescale radius to $R_{\rm B}$, temperature to $T_{\rm disk}$, and gas density to $\rho_{g,\rm disk}$. Then the preceding equations can be written

[TABLE]

where

[TABLE]

and the luminosity is expressed in terms of a core accretion time,

[TABLE]

Taking $\kappa_{d}=230(T/100~{}{\rm K})^{1+\beta}$ cm2 g*-1* (see Appendix A), the characterisic optical depth at the Bondi radius is

[TABLE]

and

[TABLE]

The solution to Equations (3.1) is shown in Figure 2 for $\beta=0$ and a range of values of the parameter $\varepsilon_{\rm rad}$ ($10^{-4}$, $10^{-3.5}$, $10^{-3}$, … $10^{0.5}$, 10). When $\varepsilon_{\rm rad}\ll 1$, the temperature remains flat some distance inside $R_{\rm B}$, whereas the density exponentiates inward,

[TABLE]

The entropy grows with radius in this outer zone (red solid curves in Figure 2). The outer envelope forms a stable radiative layer, and represents a “buoyancy barrier” such as is seen in the simulations of Kurokawa & Tanigawa (2018) (which however did not include planetesimal heating or dust production by pebble fragmentation). The density in the inner power-law core of the envelope is enhanced by a factor $\varepsilon_{\rm rad}^{-1}$.

The inner power-law profile obtained from Equations (3.1), which is shown as the dotted red lines in Figure 2, is convectively unstable. The adiabatic extension of the outer radiative profiles appears as the blue curves in Figure 2. Convection sets in when $d\ln T/dr<(\gamma-1)d\ln\rho/dr$, corresponding to $(\gamma-1)(2-\beta)<1$ and $\beta>-{1\over 2}$ for an adiabatic index $\gamma\simeq 1.4$. Then the inner density profiles scales as $\rho(r)\propto r^{-1/(\gamma-1)}$. The inner and outer solutions cross at

[TABLE]

A general feature of this simplified solution, which survives in the more complete numerical treatment, is that a rising dust abundance increases $\varepsilon_{\rm rad}$, reduces the density in the inner envelope, and pushes the radiative-convective boundary $R_{\rm rad-con}$ outward. One observes in Figure 2 that the profile is entirely convective for $\varepsilon_{\rm rad}\gtrsim 0.3$.

In particular, to determine whether the profile is convective at the Bondi radius, we substitute $\widetilde{\rho}_{g}$, $\widetilde{T}\simeq 1$ into Equations (3.1) and find the condition

[TABLE]

Since $\varepsilon_{\rm rad}$ is proportional to $X_{d}/(1-X_{d})$, e.g. $\varepsilon_{\rm rad}=\varepsilon_{\rm rad}^{*}X_{d}/(1-X_{d})$, this condition is equivalent to

[TABLE]

Equivalently, for an ambient temperature of $\sim 100$ K and density $\sim 10^{-11}$ g cm*-3*, convective refluxing across the Bondi radius can occur if the accretion rate rises above a modest threshold corresponding to $t_{{\rm acc},c}\sim 10$ Myr. This is only a few percent of the accretion rate that is needed to build the core up to $\sim 10\,M_{\oplus}$ over a Myr.

This model can be easily extended to allow for a steady, inward radial drift of solid dust particles through the outer radiative layer of the envelope. The fragmentation of pebbles in the outer envelope sources dust grains, which must drift inward in the absence of convective motions. A simple first estimate of the dust density is obtained by balancing the dust mass flux with the imposed pebble mass flux, and assuming that the dust speed is equal to the fragmentation speed,

[TABLE]

As a result, the dust size is typically $a_{d}<\lambda_{\rm max}(T)$ (Equation (5)), and so the dust opacity is independent of $a_{d}$ (see Appendix A) and depends mainly on $\dot{M}_{p}$ and $V_{f}$. The dust mass fraction solves

[TABLE]

Figure 3 shows how the boundary $R_{\rm rad-con}$ between the inner convection zone and the outer radiative layer depends on the imposed solid accretion rate, for $M_{c}=0.3\,M_{\oplus}$, $T_{\rm disk}=10^{2}$ K, and $\rho_{g,\rm disk}=10^{-11}$ g cm*-3*. The core luminosity is taken to be $GM_{c}\dot{M}_{p}/R_{c}$, which is self-consistent as long as the convection zone remains confined approximately to the Bondi sphere. One observes that refluxing sets in when the core e-folding time $t_{{\rm acc},p}\lesssim 4$ Myr; otherwise all the trapped pebble material should be accreted. This constraint is in good agreement with the time-dependent envelope calculations presented in Section 5. By contrast, the envelope models of Ormel (2014) and Mordasini (2014) find more compact convective zones because they allow for adhesive but not fragmenting collisions between grains, and therefore incorporate a much larger grain size and inward drift speed.

3.2 Dust Size Reduction by Turbulent Stirring and Secular Drift

To proceed further, one needs a model for the evolution of the dust size, which is described in Section 2.4 and implemented in the numerical model of Sections 4 and 5. Here we estimate analytically the equilibrium dust size in a fully developed convective zone, both in the outer parts where dust particles are in the Epstein drag regime ($a_{d}<\ell_{g}$) and the turbulent acceleration dominates, and in the inner envelope where the Stokes drag regime. Our estimates agree well with the model results.

In a fully convective state, the spherical temperature and density profiles around the core take the simple form

[TABLE]

These profiles connect smoothly to a uniform ambient medium. In this section, we choose a ratio of specific heats $\gamma=1.4$ and mean molecular weight $\mu=2.3\,m_{u}$, as appropriate for an atmosphere dominated by H2/He.

The luminosity is normalized in terms of the core accretion time $t_{{\rm acc},c}$,

[TABLE]

(This expression assumes high convective efficiency, the justification for which is given in Appendix B.) The convection remains subsonic everywhere in the atmosphere as long as $L_{\rm core}$ ($t_{{\rm acc},c}$) remain below (above) a critical value:

[TABLE]

This is maximized at $r=2R_{\rm B}/7$, where ${\cal M}_{\rm con}<1$ as long as

[TABLE]

and

[TABLE]

To determine the marginally fragmenting grain size in the outer, turbulent convective envelope, we make use of Equation (16). Substituting for $L_{\rm con}$ in terms of the accretion time (Equation (37)) gives

[TABLE]

At $r\sim 2R_{\rm B}/7$ the coefficient works out to $\sim 0.03$ $\mu$m for silicate grains with fragmentation speed $\sim 1$ m s*-1*, increasing by a factor $\sim 10^{2}$ for icy grains with fragmentation speed 10 times higher.

In the inner envelope, the dominant acceleration acting on grains can either be due to turbulence or the central core gravity, with the latter dominating when the planetesimal accretion rate is relatively low ($\lesssim 3\times 10^{-8}\,M_{\oplus}$ yr*-1*). In the latter case, taking $V_{f}=1$ m s*-1* as appropriate for the silicate grains that would be present at temperatures well above 100 K, one has from Equations (9) and (17)

[TABLE]

This shifts to a scaling $a_{d}(r)\propto r$ when the grain size drops below $\ell_{g}\propto r^{5/2}$ (Epstein regime).

Figures 4 and 5 give examples of how the fragmentation model of Section 2.4 is implemented in our numerical model. The solid curves in Figure 4 show the size as limited by fragmentation in response to turbulent stirring, and the dashed curves in response to radial drift in the gravitational field. The large break at $\sim 10^{10}$ cm represents a transition from amorphous silicate-rich conglomerates with fragmentation speed $V_{f}\sim 1$ m $s^{-1}$ to ice-rich grains which stick at higher speeds, $V_{f}\sim 10$ m s*-1*. The innermost breaks in both red and blue curves represent a transition from an inner Stokes drag regime to the outer Epstein drag regime.

3.3 Expansion of a Dust-Loaded Convective Layer Driven by Deep Accretion Heating

Starting from an initial envelope profile with an outer radiative layer, our numerical simulations show that most of the pebble mass is deposited in this outer layer. The differential growth in mean molecular weight $\mu$ in the radiative layer, relative to the interior convective zone, creates an additional source of instability near the radiative-convective boundary. This enhances the growth of the convective layer both by reducing the radiative energy flux and through salt-finger mixing across the boundary. Note also that if a majority of the pebble material were deposited in the deep convective layer, it would rapidly mix within a much larger mass as compared with the outer radiative shell, thereby limiting the relative growth of $\mu$ in the convective layer.

Such salt-fingering effects are not taken into account in our calculations. We use a simplified transport model in which the temperature gradient is the minimum of the adiabatic gradient and the radiative gradient (Section 4.1).

The opposing situation, where the mean molecular weight grows fastest in the convective layer, leading to the appearance of a stabilizing mean molecular weight gradient at its outer boundary, is typically not realized in our simulations. Even in such a circumstance, the rapid rise in dust opacity in the convective layer, combined with a small continuing injection of heat, would force the convective layer to expand. Supposing a fixed mass $\Delta M$ of lower mean molecular weight material to sit outside a radius $R_{\rm con}$, a fractional change $\delta s/s$ in the entropy function $s=P/\rho^{\gamma}$ drives a relatively large expansion $\delta R_{\rm con}/R_{\rm con}\sim 5\delta s/s$, where $M_{\rm con}$ is the mass of the convective layer.

To show this, consider a situation where the convective layer is sandwiched between a dense core (radius $R_{c}$), and thin outer radiative layer positioned at $r>R_{\rm con}\gg R_{c}$. The total gravitating mass is dominated by the core mass $M_{c}$, and the masses $M_{\rm con}$ and $M_{\rm rad}$ of the convective and radiative layers are taken to be fixed. The convective layer is adiabatic, with a pressure profile $P(r)=P(R_{\rm con})[\rho/\rho(R_{\rm con})]^{\gamma}$. Substituting this equation of state into the equation of hydrostatic equilibrium, and defining $s=P/\rho^{\gamma}$, one finds

[TABLE]

The convective mass is

[TABLE]

where $f(\gamma)\equiv\int_{0}^{1}x^{2}(x^{-1}-1)^{1/(\gamma-1)}dx$. Applying pressure balance across the outer convective-radiative boundary gives $s\cdot[\rho(R_{\rm con})]^{\gamma}=GM_{c}M_{\rm rad}/4\pi R_{\rm con}^{4}$. One sees that both terms in Equation (44) are functions of $s\cdot R_{\rm con}^{4-3\gamma}$. Taking $\delta M_{\rm con}=0$ then gives

[TABLE]

This implies a strong dependence of $R_{\rm con}$ on $s$, namely $\delta R_{\rm con}/R_{\rm con}\sim 5\delta s/s$, for an equation of state dominated by diatomic gases.

For all of the preceding reasons, we neglect the influence of a mean molecular weight gradient on the radiative-convective transition (explored for example by Théado & Vauclair 2012; Leconte & Chabrier 2012; Vazan et al. 2016; Lozovsky et al. 2017), except through its indirect effect on the opacity.

4 Numerical Model

We now introduce a spherical hydrostatic model of the gaseous envelope surrounding a growing planetary core, which is embedded in a PPD. A spherical approximation to the envelope structure is well motivated inside the core Bondi radius. The envelope metallicity generally exceeds the Solar value, being dominated in the outer parts by solid grains and in the inner parts by water and silicate vapor derived by pebble accretion. A luminosity $L_{\rm tot}(r)$ is transported through the envelope by a combination of convection and radiative diffusion.

4.1 Atmospheric Structure

The evolution of the envelope is calculated using a sequence of static models, governed by the equations of hydrostatic equilibrium and thermal diffusion,

[TABLE]

and

[TABLE]

Here $P$ and $T$ are total pressure and temperature, and $M(r)$ the enclosed mass. The temperature gradient $\nabla\equiv{d\log T}/{d\log P}$ is taken to be the minimum of the radiative and adiabatic gradients,

[TABLE]

as defined below. The gas mass density is determined from $P$ and $T$,

[TABLE]

The mean dust mass density $\bar{\rho}_{d}=X_{d}\rho_{g}/(1-X_{d})$ is evolved separately, as described below.

4.2 Effects of Ice and Silicate Sublimation

We numerically split the envelope into two parts. (i) An interior zone where water ice is absent, where we use tabulated equations of state (EOSs). For water vapor we use the NIST/STEAM V3.0 EOS (Harvey & Lemmon, 2013), and for H2/He the EOS of Militzer & Hubbard (2013); Militzer (2013). (ii) The outer envelope, where ice is present but the gases are dilute, and a simpler ideal gas approximation is adopted for both water vapor and H2/He.

The density of a hydrogen/helium/water vapor mixture is calculated using the additive-volume rule (e.g. Fontaine et al. 1977),

[TABLE]

Here $X_{\rm H}$ and $X_{\rm H_{2}O,v}=1-X_{\rm H}-X_{\rm ice}-X_{\rm sil}$ are the mass fractions of H2/He and water vapor; $X_{\rm ice}$ and $X_{\rm sil}$ are the mass fractions in water ice and silicates. Turbulent mixing adjusts the total H2O mass fraction in vapor and ice components to a nearly uniform value in the convection zone, $X_{\rm H_{2}O}=X_{\rm ice}+X_{\rm H_{2}O,v}\simeq$ const.

The sublimation of water ice extracts heat from the gas, according to the heat equation

[TABLE]

where $s$ is the specific entropy, $l_{\rm H_{2}O}=2.3\times 10^{10}$ erg g*-1* is the latent heat of sublimation and $C_{p}$ is the specific heat at constant pressure,

[TABLE]

The coefficients here correspond to 2 rotational degrees of freedom for H2 and 3 for H2O. We take $\mu_{\rm H}=2.3\,m_{u}$ to represent a H2/He mixture, and $\mu_{\rm H_{2}O}$ is the mass of a water molecule. Substituting this into the heat equation above and inserting the equation of hydrostatic equilibrium gives a fourth identity to supplement Equations (46), (47), (48) and (50),

[TABLE]

The adiabatic gradient $\nabla_{\text{ad}}$ is now defined as follows. The total pressure in the water sublimation zone is

[TABLE]

where $P_{\rm H_{2}O}^{\rm sat}=P_{\rm H_{2}O}^{0}e^{-l_{\rm H_{2}O}\mu_{\rm H_{2}O}/kT}$ is the saturation pressure, $P_{\rm H_{2}O}^{0}=3.6\times 10^{13}$ dyne cm*-2*. Differentiating Equation (54) with respect to radius and combining with the equation of hydrostatic equilibrium gives

[TABLE]

Combining this with Equation (53) gives

[TABLE]

where

[TABLE]

and $C_{p}$ is given by Equation (52). Finally,

[TABLE]

These equations are equivalent to those derived by Ingersoll (1969) for a wet adiabat. We emphasize that Equation (56) is only valid in the water ice sublimation zone, where H2/He gas and water vapor are approximated as ideal gases. In the inner envelope, where ice has entirely vanished, $\nabla_{\text{ad}}$ is retrieved directly from the tabuled EOSs.

Fully adiabatic profiles solving Equations (46), (47), (50), (52), (56), and (57) are shown in Figure 6 in the case of an ideal gas ($\gamma=1.4$ ratio of specific heats) EOS. A broad plateau in the temperature in the outer envelope corresponds to the ice sublimation layer. The innermost envelope around temperature 1700 K also shows the effects of silicate sublimation. These are computed in an analogous manner by considering pure silica, and considering the reaction ${\rm SiO}_{2}\leftrightarrow{\rm SiO}+{1\over 2}{\rm O}_{2}$, implying a mean molecular weight $\mu_{\rm sil}=40m_{u}$ in the gas phase (Thompson & Stevenson, 1988). The sublimation energy is $\ell_{\rm sil}=1.6\times 10^{11}$ erg g*-1*, and the vapor saturation pressure $P_{\rm sil}^{\rm sat}(T)=3.2\times 10^{14}e^{-(6\times 10^{4}{\rm K})/T}$ (Krieger, 1967).

Our default envelope models do not include the effects of silicate sublimation on the settling of dust across the core-envelope boundary. In Section 5.3, we consider how the imposition of a sublimation barrier to dust settling changes the global mass flow of solids at low core luminosity. One observes in Figure 6 that a relatively weak temperature plateau is associated with silicate sublimation, sitting just above the inner computational boundary. On that basis, we set aside the effects of heat exchange between solid and vapor silicates when considering the effect of dust sublimation on mass transport.

4.3 Radiative Transport

The radiative temperature gradient is computed throughout the envelope as

[TABLE]

where

[TABLE]

is the net gas $+$ dust opacity. The radiative luminosity has several contributions,

[TABLE]

The accretion luminosity in pebbles (mass accretion rate $\dot{M}_{p}$) is approximated as follows. A first contribution represents the disintegration of the pebbles due to water ice sublimation; the second represents the settling of silicate pebble debris onto the core:

[TABLE]

Here, $R_{\rm sub}$ is radius inside of which water ice disappears. The mean radial drift of the dust through the gas is dissipative, releasing the gravitational binding energy of the pebble material,

[TABLE]

This quantity grows in magnitude toward the core, hence our approximation in Equation (62). The accretion luminosity of larger planetesimals which can fully penetrate the envelope is

[TABLE]

Finally $L_{\rm 26Al}$ is the luminosity emitted from the envelope due to the decay of 26Al contained in dust within the envelope,

[TABLE]

Here $\tau_{\rm 26Al}=1.0$ Myr. Both $L_{\rm 26Al}$ and the hydrostatic contraction luminosity $L_{\rm KH}$ of the core and inner envelope are typically negligible compared with the two preceding luminosities, for the relevant density ($\rho\gtrsim 1$ g cm*-3*) and opacity ($\kappa\gtrsim 1$ cm2 g*-1*; Lee & Chiang 2015).

The gas opacity $\kappa_{g}$ is represented following Equations (3)-(5) in Freedman et al. (2014). This prescription depends on the gas metallicity and is thus suitable for a water vapor-enriched atmosphere. One defines

[TABLE]

where

[TABLE]

and

[TABLE]

where $c_{xx}$ are the coefficients shown in Table 1 and “${met}$” is the gas metallicity.

The grain opacity is taken to be the minimum of twice the geometric opacity and the Rosseland mean small-grain opacity $\kappa_{d,\rm R}$ derived in Appendix A,

[TABLE]

Here, $\kappa_{\rm geom}=3X_{d}/4\rho_{s}a_{d}$ for spherical grains of radius $a_{d}$, material density $\rho_{s}$, and mass $X_{d}$ per unit total mass, and

[TABLE]

where $\lambda_{\rm max}(T)\equiv hc/4.95\,k_{\rm B}T$ is the peak wavelength in the Planck function. The coefficient $Q^{\prime}=0.35$ for small silicate-carbon grains, but is a factor $\sim 1/(6-7)$ smaller when these grains are coated with an icy mantle. Dust size $a_{d}$ is calculated self-consistently using Equations (16) and (17).

4.4 Dust Transport

The transport of small dust grains through the envelope is governed by the advection-diffusion equation (Thoul et al., 1994):

[TABLE]

Here, $C=\rho_{d}/(\rho_{d}+\rho_{g})$ is the dust concentration, $D_{\rm tot}$ is the total effective diffusion coefficient, and the second term on the right-hand side represents inward settling of the grains in the central gravitational field, ${\bf V}_{\rm sett}=-\tau_{\rm stop}g(r)\hat{r}=-\tau_{\rm stop}GM(r)\hat{r}/r^{2}$. Although settling can dominate turbulent stirring as a source of differential grain velocities in the inner envelope, the bulk transport of small grains is dominated by convection wherever it is present. The source function $\dot{\rho}_{p}/(\rho_{d}+\rho_{g})$ represents continuous pebble accretion, and therefore is centered at the pebble destruction radius as determined via Equations (10) and (11). We convolve the source term with a narrow normal distribution to enhance numerical stability.

We solve Equation (71) numerically using FiPy: A Finite Volume PDE Solver Using Python333https://www.ctcms.nist.gov/fipy/ (Guyer et al., 2009). Although $D_{\rm tot}$ in principle receives a contribution from molecular diffusion, eddy diffusion completely dominates in convective zones, and the advection term dominates in radiative zones. We set $D_{\rm tot}=D_{\rm con}+D_{\rm mol}$, where

[TABLE]

and

[TABLE]

Here, $H$ is the pressure scale height (corresponding to a mixing length equal to $H$), and ${\rm St_{d}}$ is the dimensionless dust Stokes number ${\rm St_{d}}=V_{\rm con}\tau_{\rm stop}/H$, with $\tau_{\rm stop}$ defined in Equation (9).

The superadiabatic gradient $\nabla-\nabla_{\rm ad}$ is obtained from a mixing length model of convection, following the prescription of Cox & Giuli (1968). One first defines the ratio of convective to radiative conductivity,

[TABLE]

where $Q=1-(\partial\ln\mu/\partial\ln T)_{P}$ with $\mu$ the mean molecular weight, and $C_{p}$ the heat capacity at constant pressure, as derived from the EoS. We then define the quantity

[TABLE]

Finally, the convective efficiency

[TABLE]

is obtained by solving the polynomial equation

[TABLE]

Once $\zeta$ is known, one immediately obtains $\nabla-\nabla_{\rm ad}$ and the convective velocity from Equation (72).

Note that in the mixing length model, $V_{\rm con}\propto\sqrt{\nabla-\nabla_{e}}$, where $\nabla_{e}$ is the temperature gradient of the rising element. This factor can be approximated as $\sqrt{\nabla-\nabla_{\rm ad}}$ when the element experiences minimal heat loss by radiation, corresponding to efficient adiabatic convection. The applicability of this condition, and of the corresponding simple approximation $V_{\rm con}^{3}\simeq L_{\rm con}/\rho$ (implying comparable contributions of the kinetic energy and enthalpy perturbation to the energy flux) is reviewed in Appendix B.

The dust mass fraction in the outer envelope is $X_{d}=X_{\rm ice}+X_{\rm sil}$. The ice mass fraction is easily derived by subtracting the saturation water vapor density from the net solid$+$vapor water mass density, $X_{\rm ice}=X_{\rm H_{2}O}-P_{\rm H_{2}0}^{\rm sat}(T)\mu_{\rm H_{2}O}/k_{\rm B}T\rho$. The analogous procedure is used for silicates when we construct alternative models with dust sublimation in the inner envelope.

4.5 Initial Conditions and Boundary Conditions

The envelope is initialized as a solution to the hydrostatic and radiative Equations (46)-(50) and (56)-(60) with $M_{c}=0.3\,M_{\oplus}$, uniform dust abundance $X_{d}=X_{d,\rm disk}=0.01$, $X_{\rm ice}=X_{\rm sil}=X_{d}/2$, and grain size $a_{d}=1$ $\mu$m (Table 2).

Constant conditions $\rho_{g}=\rho_{g,\rm disk}=10^{-11}$ g cm*-3*, $T=T_{\rm disk}=100$ K and $X_{d}=X_{d,\rm disk}$ are maintained at the outer computational boundary $R_{\rm max}=4\times 10^{10}(M_{c}/0.3~{}M_{\oplus})$ cm $\gtrsim R_{\rm B}$. Grains produced by the destruction of larger pebbles are allowed to flow outward across this boundary and, since the external grain density is held fixed, this outflow represents a permanent loss from the envelope to the disk.

The inner computational boundary sits at the core radius $R_{c}=5.2\times 10^{8}$ cm. Here an open boundary condition is also applied for silicate grains, meaning that settling grains that reach the core are accreted and removed permanently from the envelope. In our default configuration, we do not allow silicates to vaporize in the inner envelope even when this is thermodynamically possible. The effect of an inner silicate sublimation barrier (Section 2.5) on the envelope mass flow is considered as an alternative case.

As regards the outer boundary condition on $X_{d}$, two arguments suggest that the equilibrium metal concentration in the envelope is relatively insensitive to the ambient metallicity. The high dust abundance that is reached in our envelope models (10-30 times Solar) is mainly due to the weak coupling of pebbles to the gas around the Bondi radius, which allows the solids to concentrate with respect to the H2/He gas. An outflow of metal-enriched gas will increase the metallicity of the co-orbital region in the surrounding disk, and so the feedback on the equilibrium envelope metallicity should be considered.

First, we tested the consequences of a free-flow outer boundary condition, by extending the computational domain to well beyond the Bondi radius. Dust transport was tracked in this outer zone, and convection allowed to develop there. We found the results inside radius $R_{\rm B}$ to differ minimally from those obtained from the default boundary conditions. Second, one can estimate analytically the sensitivity of the equilibrium $X_{d,\rm env}$ in the envelope to the ambient $X_{d,\rm disk}$. Given a fixed pebble accretion rate, and a steady state relation between solid accretion and excretion rates,

[TABLE]

one sees that it is difference $X_{d,\rm env}-X_{d,\rm disk}$ in dust abundance between envelope and disk that is constrained. If we were arbitrarily to raise $X_{d,\rm disk}$ from the assumed Solar abundance up to the value that is reached in a given envelope model, then the new envelope metallicity would be a factor of 2 larger. In practice, the dust component of the ejected gas must be diluted by mixing with lower-metallicity disk gas, or by the reassembly of small grains into pebbles.

The entropy of the inner envelope does not initially match the ambient disk entropy, but continuing solid accretion combined with secular dust accumulation in the envelope removes this imbalance in entropy. This result can be seen both semi-analytically and numerically. The steady envelope solution presented in Section 3.1 has a uniform entropy in the inner convection zone that is lower than the entropy of the ambient medium. (In this model, the central radiation source that drives the convection represents deep planetesimal impact heating.) These inner and outer zones are connected by a layer with radially increasing entropy (Figure 2). As the envelope is loading with an increasing amount of dust (here produced by pebble destruction), the entropy of the convection zone rises. Finally the entropy is nearly uniform in the fully convective state with the parameter $\varepsilon_{\rm rad}\gtrsim 0.3$. Only then is the pebble debris recycled back into the disk.

Figure 7 shows the time evolution of the gas ‘entropy’ function $P/\rho_{g}^{7/5}$ for our nominal time-dependent model with icy pebbles (see Section 5) and heating driven by the accretion of pebble debris. The gas entropy is flat in the inner convection zone, and rises through the outer envelope due to a combination of radiative transport and also the extraction of heat by ice sublimation. As time progresses, the radiative layer shrinks and the convective plateau grows. The outer entropy bump does not entirely disappear in the fully convective state, due to the continuing sublimation of ice. See Section 5.3 for further details.

In sum, we do not expect that a “buoyancy barrier” should be present when envelope heating overwhelms cooling. The numerical model constructed by Kurokawa & Tanigawa (2018) does not include planetesimal heating or dust production. A fuller exploration of the interplay between heating, cooling, and hydrodynamic expulsion will require incorporating these effects into global simulations.

5 Results

We now describe the time-dependent behavior of the accreted atmosphere, while it remains in thermal contact with the disk. We ran simulations using the model described in Sections 4.1-4.4 and the boundary conditions summarized in Section 4.5, with a uniform set of parameters summarized in Table 2. Our focus is on the early growth of a core, $M_{c}=0.3\,M_{\oplus}$. The default model assumes pebbles of mixed silicate/ice composition ($X_{\rm ice}=X_{\rm sil}$), which we compare with the case of dry pebbles ($X_{\rm ice}=0$). The chosen core mass is just above the minimum mass that accretes a high-entropy H2/He envelope hot enough to completely sublimate silicates in a finite layer outside the core (when the core itself is silicate-rich). Our default model assumes that the core-envelope boundary is permeable to silicate dust. We also compare the results for envelope and core mass accumulation with an alternative situation where silicate sublimation prevents growth of the core by silicate ‘rain’ (see Section 2.5 for motivation).

The ambient disk temperature and H2/He mass density are taken to be $T_{\rm disk}=100$ K and $\rho_{g,\rm disk}=10^{-11}$ g cm*-3* (Table 2), corresponding to moderate coupling of the pebbles to the gas near the Bondi radius (Equation (2)). The outer and inner computational boundaries sit at $R_{\rm max}=4\times 10^{10}$ cm $\gtrsim R_{\rm B}$ and $R_{c}=5.2\times 10^{8}$ cm. Gas density and temperature are matched to the disk values at radius $R_{\rm max}$, and pebbles are introduced into the envelope in the manner described in Section 4.

A principal goal is to determine the minimum luminosity (core accretion rate) that is needed to drive refluxing of the small particles across the Bondi radius. We fix the pebble accretion rate $\dot{M}_{p}$ at a conservative values of $10^{-7}$, $10^{-6}$ M⊕ yr*-1* (the larger quantity is needed to e-fold the core mass up to $\sim 10\,M_{\oplus}$ in a Myr) and test two values (0 and $10^{-7}$ M⊕ yr*-1*) for the planetesimal accretion rate $\dot{M}_{\rm plan}$. The core is assumed to have a rocky composition; to obtain the same heating rate by accretion onto a rock-ice core with mean density $\widetilde{\rho}_{c}\times 3$ g cm*-3*, the accretion rate would need to rise by a modest factor $1.6(\widetilde{\rho}_{c}/0.5)^{-2/3}$.

5.1 Sandblasting and Sublimation of Pebbles

Infalling pebbles are rapidly eroded by collisions with small ambient grains and by sublimation, as is seen in Figure 8. After $X_{d}$ builds up in the outer envelope, destruction is almost complete a small distance inside the outer boundary of the computational domain. We show two cases, one in which the outer envelope remains radiative ($\dot{M}_{p}$ = $10^{-7}M_{\oplus}$ yr*-1*), and one in which convection extends beyond the Bondi radius ($\dot{M}_{p}=10^{-6}M_{\oplus}$ yr*-1*).

One see that, although destruction is faster in the presence of convection, secular drift in the central gravitational field is still strong enough to eliminate the pebbles quickly in the outer radiative part of a low-luminosity envelope. Pebble debris settles rapidly toward the inner convective layer: at the fragmentation speed $V_{f}\sim 10$ m s*-1* for icy particles, this takes a month or so, and for silicate particles about a year. The dust opacity in the outer envelope saturates over this timescale (Section 5.3).

5.2 Thermal and Hydrostatic Structure of the Envelope

We now analyze the radial structure of the accreted envelope. Two effects produce a temperature plateau in the outer envelope: (i) hydrostatic equilibrium in a radiative atmosphere with outer temperature buffered close to the disk temperature by dust-gas heat exchange (see Figure 2); and (ii) sublimation of water ice, which limits the growth of the gas temperature above the water sublimation temperature.

The simplest case of nearly ice-free pebbles ($X_{\rm ice}=10^{-4}$) is shown in Figure 9, for both low and high planetesimal accretion rates. The right panel shows the settling rate of dust embedded in the atmosphere due to the d.c. component of the radial drift (Equation 63). In the low-$\dot{M}_{p}$ model (which does not produce refluxing of dust across the Bondi radius), this is an appreciable fraction of the imposed pebble accretion rate. At the core boundary there is a balance between pebble accretion and transmission of dust. By contrast, when the accretion rate is high enough to drive dust refluxing (and limit the growth of the envelope), the accretion rate of silicate dust onto the core is about $0.1\,\dot{M}_{p}$.

The core accretion rate $\dot{M}_{\rm sett}(R_{c})$ saturates at $\sim 1\times 10^{-7}\,M_{\oplus}$ yr*-1* for intermediate values of the pebble accretion rate (see Table 3). Convection reaches the Bondi radius when $\dot{M}_{p}$ slightly exceeds this value, with most of the pebble debris continuing to accrete onto the core until $\dot{M}_{p}$ rises above $\sim 2\times 10^{-7}\,M_{\oplus}$ yr*-1*. We therefore derive a core mass e-folding time due to pebble accretion no shorter than $\sim 3$ Myr: trapping of pebbles at higher rates implies higher rates of refluxing but no faster core growth.

Turning next to the case of icy pebbles ($X_{\rm ice}=X_{\rm sil}$), one sees in Figure 10 that the temperature profile shows an outer plateau at both low and high accretion luminosities. The inner density is about a factor 2 higher and the dust settling rate is similar to the case of ice-free pebbles.

The spatial distributions of H2/He, silicate dust, ice, and water vapor mass are shown in Figure 11 for the higher $\dot{M}_{p}$ model. Ice is present beyond a radius $R_{\rm sub}=9\times 10^{9}$ cm, which marks the termination of the temperature plateau seen in Figure 10.

5.3 Evolution of the Envelope Mass and Metallicity in the Baseline Model

The growth of the dust mass embedded in the envelope is shown in Figures 12 and 13 for the default case of icy pebbles. We present a sample of pebble accretion rates, extending between $10^{-7}$ and $10^{-6}\,M_{\oplus}$ yr*-1*. Particles moving near the fragmentation speed settle rapidly through the outer radiative layer. At the two highest accretion rates, $\dot{M}_{p}=(0.3,1)\times 10^{-6}M_{\oplus}$ yr*-1*, the envelope is saturated with dust after a few thousand years, and the H2 mass also saturates (red curves in Figure 13). Beyond this point, convection continuously ejects most of the accreted solids back into the disk in the form of small grains. That is,

[TABLE]

On the other hand, in the most slowly accreting model ($\dot{M}_{p}=1\times 10^{-7}M_{\oplus}$ yr*-1*), the core mass continues to grow at the imposed pebble accretion rate. The masses of dust and gas in the envelope still reach a steady state, showing an almost identical gas/dust ratio near unity (blue curves in Figures 12 and 13).

Here we have made the default assumption that the core-envelope boundary is transmitting to silicate dust. Alternatively, when sublimation of silicate dust does impose a strong barrier to core growth, a very different mass partitioning between core and envelope is obtained. It is the dust mass stored in the envelope which grows as pebbles continue to accumulate; the core mass remains frozen (cyan curves in Figure 13).

The preceding results, which are summarized in Table 3, can be readily understood in terms of the evolving position $R_{\rm rad-con}$ of the outer radiative-convective boundary (RCB), as shown in Figure 14. At the higher pebble accretion rates, the bound envelope begins with a significant radiative layer but quickly captures pebble material. The dust created in this outer layer rapidly settles toward the interior, and the RCB pushes rapidly outward in reseponse to the growing opacity, until convection extends past the Bondi radius. By contrast, in the lowest-luminosity ($\dot{M}_{p}=10^{-7}\,M_{\oplus}$ yr*-1*) case the RCB never approaches the Bondi radius.

Recycling of pebble debris sets in at lower $\dot{M}_{p}$ if a modest fraction of the solid accretion flux is carried by large, penetrating planetesimals. For example, the $\dot{M}_{p}=10^{-7}\,M_{\oplus}$ yr*-1* model becomes fully convective after adding $\dot{M}_{\rm plan}=10^{-7}\,M_{\oplus}$ yr*-1* in planetesimal accretion heating (Table 3).

One observes in the high-luminosity models with refluxing a rapid limit cycle involving the transient appearance of a narrow radiative layer just interior to the computational boundary (Figure 14). This effect is driven by the finite lifetime of the convective eddies: the accretion of low-metallicity gas pushes the RCB inward, but this material then mixes with the dusty envelope, forcing the convection back outward. This cycle is shown in more detailed in Figure 15, which plots the spatial and temporal evolution of the dust mass fraction. Dust deposited in the outer radiative envelope due to sandblasting settles inward quickly, pushing the RCB outwards. When the envelope becomes fully convective, it loses dust in the outer scale height or so, which creates a transient radiative zone. As this cycle continues, the dust mass builds up in the envelope; the radial width of the transient radiative zone decreases with time, until the envelope reaches a steady state around $t\sim 10^{3}$ yr.

This effect is, to some extent, an artifact of a one-dimensional treatment of the convection, combined with the presence of a sharp outer computational boundary. More realistically, the dust abundance is significantly larger in upflowing convective plumes than in gas accreted from beyond the Bondi radius, with the latter providing opacity holes through which much of the radiation escapes. One can expect a more gradual transition in $X_{d}$ when the surrounding disk is included in the computational model.

5.4 Dust Size Distribution

Details of the dust size distribution have already been presented in Figures 4 and 5. The dust size is limited by fragmentation, either resulting from turbulent stirring (Equation (16)) or from secular drift (Equation (17)), with the drag law (Epstein or Stokes) being determined self-consistently. Figure 4 shows that turbulent fragmentation dominates near the Bondi radius, but the secular drift speed is higher inside $\sim 10^{10}$ cm from the core. A substantially larger dust size also results near the Bondi radius in cases where the outer envelope is radiative.

6 Summary and Discussion

We have identified a significant constraint on the growth of planetary cores by the accretion of small (mm-cm sized) pebbles. When ambient conditions are optimal for pebble accretion, pebbles are susceptible to fragmentation (Section 2). Pebble accretion is most efficient when the drift speed of a pebble through the embedding gas toward the local center of gravity is comparable to flow speed of the gas near the the core. In such a circumstance, the pebble drift speed far exceeds the fragmentation speed of an amorphous conglomerate of (sub)-micron sized grains. The water ice component of pebbles is rapidly lost by sublimation, releasing embedded silicate grains and preventing the water ice from reaching the core.

A high pebble accretion rate, as needed to fuel the rapid growth of a small core into a Neptune-mass planet ($t_{{\rm acc},p}=M_{c}/\dot{M}_{p}\lesssim 3\times 10^{5}$ yr), implies (i) a significant dust mass fraction in the medium surrounding the core, and (ii) efficient binary collisions between pebbles and, even more so, between pebbles and dust grains suspended in the gas. The build-up of small particles in the envelope, which do not accrete efficiently onto the core, triggers rapid sandblasting of high-speed pebbles around the Bondi radius. Indeed, convection sets in within the outer envelope only when the mass fraction of particles smaller than $\lambda_{\rm max}(T_{\rm disk})\sim 30\,\mu$m is larger than $\sim 0.3$.

Pebble fragmentation does not, by itself, imply inefficient core growth. Even as their density exponentiates in the outer radiative layer, dust particles moving near the fragmentation speed will settle fairly rapidly reaching the inner convective layer in a year or less. When the imposed pebble accretion rate exceeds $\sim 10^{-7}\,M_{\oplus}$ yr*-1*, we find that convective motions grow and extend beyond the planet’s gravitational sphere of influence as defined by ambient disk conditions (Section 5). Then some of the debris of the accreted pebbles is recycled back into the disk, and the accretion rate of the remaining debris onto a $0.3\,M_{\oplus}$ core stabilizes at $\sim 1\times 10^{-7}\,M_{\oplus}$ yr*-1*, corresponding to an e-folding time of $\sim 3$ Myr. The loss of dust and gas eventually counterbalances accretion, causing the solid mass in the envelope to plateau. The envelope density profile remains shallow and its entropy equilibrates near the ambient disk entropy.

When the pebbles are icy, the water component sublimates in the outer envelope and does not accrete onto the core. Ice sublimation extracts heat from the gas in the outer envelope, which raises the inner density of the envelope. We find, however, that continuing accretion of hydrogen from the disk limits the water vapor mass fraction and maintains a signficant H2/He component in the envelope. This in turn limits the density growth in the inner envelope that is caused by heat extraction from the outer envelope.

Embedded grains do settle toward the core through the inner convective envelope. The rate at which the core gains silicate mass is estimated from the radial drift time (Equation (63)). This accretion rate $\dot{M}_{\rm sett}$ is proportional to the dust concentration in the inner envelope. Allowing for a relative concentration $f_{\rm turb}$ of dust at the base of the envelope as compared to its interior, $\dot{M}_{\rm sett}$ is less than $10\,f_{\rm turb}\%$ of the imposed pebble accretion rate, for $\dot{M}_{p}=10^{-6}\,M_{\oplus}$ yr*-1*. By comparision, a short accretion time $M_{c}/\dot{M}_{p}\sim 3\times 10^{5}$ yr needed to build a Neptune-mass planet in the lifetime of the protoplanetary disk.

The aforementioned effects were demonstrated using a self-consistent, spherical model of the bound envelope (Section 4) which includes the thermodynamics of water phase exchange, pebble destruction, dust collision and size evolution, convection in the mixing length approximation, multiple luminosity sources, dust transport handled using the advection-diffusion equation, analytical prescriptions for dust and gas radiative opacity, and tabulated non-ideal EOSs for water and hydrogen. The steady steady behavior produced by this model is in agreement with the convective/radiative structure seen in a steady-state spherical model (Section 3).

The specific example studied numerically is a core of mass $0.3\,M_{\oplus}$. A simple analytic model of a hydrostatic envelope loaded with dust gives a nearly identical result (Equation (33)) for the critical deep heating rate (due to planetesimal accretion) above which refluxing balances pebble accretion. The growth time of the core due to the settling of pebble debris in the inner envelope scales as $t_{{\rm acc},c}=M_{c}/\dot{M}_{\rm sett}(R_{c})\propto(X_{d,\rm env}\rho_{g,\rm disk})^{-1}M_{c}^{-4/3}$; here $X_{d,\rm env}$ is the envelope metallicity and in Equation (63) we assume that the inward drift speed of the pebble debris is the fragmentation speed $V_{f}$. Given a fixed ratio of pebble drift speed to ambient gas sound speed ($\rho_{g}\propto M_{c}^{-1}$ at the Bondi radius, from Equation (2)), refluxing becomes slightly more difficult as $M_{c}$ increases, $t_{{\rm acc},c}\propto X_{d,\rm env}^{-1}M_{c}^{-1/3}$. Setting aside this constraint of efficient pebble accretion, the minimum $X_{d,\rm env}$ that forces convection at the Bondi radius scales as $X_{d,\rm env}\propto\rho_{g,\rm disk}^{-1}M_{c}^{-2/3}\,t_{{\rm acc},c}$ (Equation (33)). Hence, the core grows moderately faster at the onset of refluxing as its mass increases, $t_{{\rm acc},c}\propto M_{c}^{-1/3}$ (the dependence on $\rho_{g,\rm disk}$ cancels). Finally, we may consider the minimum $X_{d,\rm env}$ needed for dust refluxing across the Bondi radius to balance pebble accretion, $X_{d,\rm env}\rho_{g}V_{\rm con}4\pi R_{\rm B}^{2}\sim M_{c}/t_{{\rm acc},p}$, where $t_{{\rm acc},p}$ is the accretion time associated with the imposed pebble trapping rate. Given the preceding scaling for $t_{{\rm acc},c}$, one finds from Equation (6) that $V_{\rm con}\propto M_{c}^{0}\rho_{g,\rm disk}^{-1/3}$. Hence pebble accretion is cancelled for $X_{d,\rm env}\propto\rho_{g,\rm disk}^{-2/3}M_{c}^{-1}t_{{\rm acc},p}^{-1}$, which once again decreases with increasing $M_{c}$.

6.1 Caveats

1. An important, unresolved issue about pebble accretion is whether it proceeds mainly from a thin, pebble-dominated disk (Goldreich et al., 2004b), or by trapping from a gas-dominated medium (Ormel & Klahr, 2010). Our calculations address the latter possibility, but may still be applicable to the first situation when the bound envelope is convective and mixes embedded solids thoroughly in all three dimensions.

2. The 3D numerical simulations of Popovas et al. (2018) follow the trajectories of solid particles through the gas around a planetary core embedded in a shearing gas disk. Allowing for heating of the gas by transfer of accretion energy from the pebbles, they find evidence for extended cold downflows, which help to advect pebbles to the core. However, they also find that particles drifting at the silicate fragmentation speed of $\sim 1$ m s*-1* can be cycled repeatedly within convective flows; more generally, the dynamics of small particles with short stopping times is not resolved near the core. Our approach to dust accretion suggests that such particles experience slow secular drift toward the core, otherwise being entrained in the rapid convective flows and surviving for a long time in the envelope.

3. The main limitation of our model is the use of a spherical approximation and the neglect of rotational effects near the Bondi radius. The physical processes we investigate need to be incorporated into a global, 3D simulation of a core embedded in a gas disk. Efficient pebble accretion without strong fragmentation might be easier if the core were surrounded by an orbiting gas disk that trapped pebbles for many orbits. Then tightly-coupled pebbles might spiral into the core at reduced speeds relative to the gas. The extent to which such centrifugally-supported structures form around low-mass cores (Bondi radius $\ll$ Hill radius) remains an open question.

4. The recycling of ejected debris back to the disk may enhance the metallicity of the accreted material, but we have argued that this does not significantly perturb the equilibrium envelope state as long as the ejected debris is diluted by mixing with disk gas.

5. The effect of silicate dust sublimation in the inner envelope was not included in our default envelope model, and is not needed to suppress core growth. Nonetheless, we argue that the re-evaporation of silicate “rain” can have a significant additional suppressing effect on solid accretion by the core when its mass exceeds $0.2-0.3\,M_{\oplus}$.

6. The size of a grain near the threshold for fragmentation is typically much smaller than the assumed pebble size of 1 mm, and is small enough to produce significant opacity ($a_{d}<\lambda_{\rm max}(T)$). Furthermore, we have argued that reassembly of spallated monomers into smaller amorphous conglomerates cannot happen fast enough to limit the exponential growth of opacity by the sandblasting process. Nonetheless, the effect of collisions of pebbles with small, marginally-fragmenting conglomerates needs further attention.

6.2 Implications and Future Directions

1. Pebble destruction by sandblasting (Schräpler & Blum, 2011; Jacquet & Thompson, 2014) may be important well outside the quasi-spherical, bound envelope treated here. Hydrodynamic models of pebbles and gas interacting with a core should take into account the feedback of dust-enriched gas back into the disk, and its effect on radiation transfer and gas entropy.

2. When the ambient disk temperature lies not too far below the sublimation temperature of water ice, we find a strong gradient in gas density near the Bondi radius, associated with the absorption of heat by sublimating ice. This means that the gas envelope accreted around a core can maintain a well-defined boundary even while the envelope remains thermodynamically coupled to the disk, that is, while it maintains a similar specific entropy.

3. Our models with a reduced luminosity – which fail to develop convection near the Bondi radius – maintain only a modest metallicity in the envelope in spite of rapid pebble accretion, typically $\sim 20$ times Solar. That is because the accretion of pebbles is accompanied by the continued accretion of hydrogen and helium, as sublimating ice absorbs heat from the envelope and the gas component contracts and cools. This result, if verified at large core masses, would argue against the formation of icy Neptune-mass planets by pebble accretion from a gas-rich disk. There is also a synergy between this result and recent accretion models of Jovian planets, which never develop a substantial core (see, e.g., Lozovsky et al. 2017); and recent measurements of the gravitational moments of Jupiter, which allow for a dilute core in that planet (Wahl et al., 2017).

4. We have not considered processes that might mix core material back into the envelope. In the generic situation considered here, the envelope has a higher entropy and lower metallicity than the core, meaning that the core-envelope boundary is not susceptible to either salt-fingering (higher metallicity and higher entropy on the outside) or semi-convection (lower metallicity and lower entropy on the outside). However, convective dredge-up could supplement the process of dust loading of a gas envelope by pebble accretion.

5. Ongoing pebble accretion has important implications for the retention of planetary atmospheres during a giant impact phase. The saturation in the growth of the envelope implies a strong limit to the gravitational binding energy of the envelope to the core (Table 2). Such a low-mass envelope is more susceptible to expulsion by giant impacts.

6. This work suggests an interesting scenario for the origin of high temperature minerals (crystalline silicates and CAIs) observed in comets (for example Hanner et al. (1994); Hanner, Lynch, & Russell (1994); Brownlee et al. (2006)). In our model, pebbles are processed thermally and mechanically at very high temperatures in the inner envelope before being convectively transported and ejected back into the disk. The combination of high temperature and high envelope metal enrichment satisfies some of the basic constraints on the birth environment of these minerals; we leave detailed modeling to future work.

The broad lesson here is that the structure of planetary atmospheres during the planet assembly phase can be dramatically altered by pebble accretion, and more generally is very sensitive to the mode of solid delivery. In the context of future 3D accretion simulations and time-dependent envelope models, our work implies that it is essential to incorporate a local fragmentation model for embedded grains and to treat their evolution separately from infalling pebbles. The combination of accretion heating by settling solids with the build-up of grain opacity is found to have a powerful feedback effect on gas flows in the outer envelope. By contrast, many current planet formation models allow for the growth of grains that are sourced by infalling planetesimals, while ignoring fragmentation feedback on the grain size and therefore on the mass fraction and opacity of settling grains (e.g. Mordasini 2014; Ormel 2014). It is also essential to simulate the interaction of lower core masses with pebble-loaded gas disks, as our model shows strong feedback effects even for $M_{c}=0.3\,M_{\oplus}$, and implies a significant modification of the solid particle population in the ambient disk through refluxing.

M.A.-D. thanks the Department of Physics at the American University of Beirut where he spent one semester as a visiting faculty. The authors thank T. Birnstiel for discussions on dust size evolution, D. Valencia for discussions on opacities and equations of state, and K. Menou for discussions that helped to motivate this project. Finally, we thank the referee Anders Johansen for comments on our treatment of collisional fragmentation.

Appendix A Small-grain Dust Opacity

Here we motivate our choice of Rosseland mean opacity $\kappa_{d,\rm R}(T)$ for small grains (of radius $a_{d}\ll\lambda_{\rm max}(T)=hc/4.95\,kT$. The grain opacity used in our atmosphere models is the minimum of this value and the geometric opacity $\kappa_{\rm geom}$ (see Equation (69)). Given a wavelength dependence $\kappa(a_{d},\lambda)=Q^{\prime}\cdot(2\pi a_{d}/\lambda)\cdot\kappa_{\rm geom}(a_{d})=Q^{\prime}\cdot 3\pi X_{d}/2\rho_{s}\lambda$ for grains of density $\rho_{s}$, the Rosseland mean works out to

[TABLE]

Here $\nu=c/\lambda$ and $B_{\nu}(\nu,T)$ is the Planck function. The opacity of small silicate-carbon grains, of density $\simeq 2.5$ g cm*-3* is derived by Ossenkopf & Henning (1994) to be $\kappa_{d}=7.0$ cm2 g*-1* at $\lambda=1.3$ mm (see their Figure 9). This corresponds to $Q^{\prime}=0.50$ and a Rosseland mean

[TABLE]

as weighted by gas mass. Equivalently, $Q_{R}(a_{d},T)=0.73\,Q^{\prime}\cdot 2\pi a_{d}/\lambda_{\rm max}(T)=0.36\,\cdot 2\pi a_{d}/\lambda_{\rm max}(T)$.

Ice coatings deposited by condensation onto silicate grain surfaces could reduce the opacity by a factor $1/(6-7)$ (Figure 9 of Ossenkopf & Henning (1994)).

Appendix B Convective Efficiency

The relative importance of convection and radiation in transporting energy through the envelope is encapsulated in the parameter $A$ given by Equation (74). This is plotted as a function of radius in Figure 16 for our highest $\dot{M}_{p}$ model, in which convection reaches the outer computational boundary and envelope growth is quenched. In both cases, $A$ attains a minimum value $O(10^{2})$ in the outer envelope, corresponding to energy transport dominated by convection.

Figure 17 shows that the simple approximation $F_{\rm con}\simeq\rho V_{\rm con}^{3}$ is adequate in the convective regions of the envelope after it has reached a steady state. In the same formulation of mixing length theory used in this paper, one finds $F_{\rm con}\simeq 0.36\alpha\rho V_{\rm con}^{3}$, where $\alpha$ is the ratio of mixing length to pressure scale height. We take $\alpha=1$ in the numerical model.

In our high-luminosity models, the superadiabaticity $\nabla-\nabla_{\rm ad}\sim 8\times 10^{-6}$ near the core, comparable to the value in the interiors of giant planets, and implying almost perfectly adiabatic convection. This quantity rises to $\sim 10^{-3}$ near the Bondi radius, where the assumption of full adiabaticity approaches the limit of applicability.