The End of Runaway: How Gap Opening Limits the Final Masses of Gas Giants

TL;DR

This paper develops a theory explaining why gas giant planets stop growing at certain masses, showing that circumstellar disc gaps limit their final size, especially at large orbital distances beyond 10 AU.

Contribution

It introduces a self-consistent, time-dependent model of gap formation that predicts the final mass of gas giants based on disc properties, independent of viscosity.

Findings

Final masses are around a few Jupiter masses at 10-100 AU.

Gap opening halts runaway accretion in low-viscosity discs.

Final mass depends on disc lifetime, aspect ratio, and density, not viscosity.

Abstract

Gas giants are thought to form by runaway accretion: an instability driven by the self-gravity of growing atmospheres that causes accretion rates to rise super-linearly with planet mass. Why runaway should stop at a Jupiter or any other mass is unknown. We consider the proposal that final masses are controlled by circumstellar disc gaps (cavities) opened by planetary gravitational torques. We develop a fully time-dependent theory of gap formation and couple it self-consistently to planetary growth rates. When gaps first open, planetary torques overwhelm viscous torques, and gas depletes as if it were inviscid. In low-viscosity discs, of the kind motivated by recent observations and theory, gaps stay predominantly in this inviscid phase and planet masses finalize at $M_{\rm final}/M_\star\sim(\Omega t_{\rm disc})^{0.07}(H/a)^{2.73}(G\rho_0/\Omega^2)^{1/3}$, with $M_\star$ the host stellar mass, $\Omega$ the planet's orbital angular velocity, $t_{\rm disc}$ the gas disc's lifetime, $H/a$ its aspect ratio, and $\rho_0$ its unperturbed density. This final mass is independent of the dimensionless viscosity $\alpha$ and applies to large orbital distances, typically beyond $\sim$10 AU, where disc scale heights exceed planet radii. It evaluates to a few Jupiter masses at 10-100 AU, increasing gradually with distance as gaps become harder to open.

Equations

$$M_{\mathrm{thermal}}\equiv 3(H/a{)}^3{M}_{\star }$$

$$R_{\mathrm{B}}=GM/c_{\mathrm{s}}^2$$

$${\Sigma }_0=70{\text{ g cm}}^{-2}{\left(\frac{a}{10\text{ AU}}\right)}^{-3/2}.$$

$$T_0=100\text{ K}{\left(\frac{a}{10\text{ AU}}\right)}^{-3/7}.$$

$${\rho }_0=\frac{1}{\sqrt{2\text{uppi}}}\frac{{\Sigma }_0}{H}=3\times 10^{-12}{\text{ g cm}}^{-3}{\left(\frac{a}{10\text{ AU}}\right)}^{-39/14}$$

$$H/a=0.08{\left(\frac{a}{10\text{ AU}}\right)}^{2/7}.$$

$${\rho }_{\mathrm{rcb}}=\frac{M}{R_{\mathrm{rcb}}^3}{\left(\frac{k_{\mathrm{B}}{T}_0{R}_{\mathrm{rcb}}}{GM\mu }\right)}^{1/(\gamma -1)}=\frac{M}{R_{\mathrm{rcb}}^3}{\left(\frac{R_{\mathrm{rcb}}}{R_{\mathrm{B}}}\right)}^{1/(\gamma -1)}.$$

$$t_{\mathrm{cool}}\sim \frac{E}{L}\sim {\left(\frac{k_{\mathrm{B}}}{\mu }\right)}^5\frac{\kappa }{\sigma G^4}\frac{T_0}{M^2}{\left(\frac{R_{\mathrm{B}}}{R_{\mathrm{rcb}}}\right)}^{(4\gamma -5)/(\gamma -1)},$$

$$t_{\mathrm{cool}}\sim {\left(\frac{k_{\mathrm{B}}}{\mu }\right)}^5\frac{\kappa }{\sigma G^4}\frac{T_0}{M^2}\approx 10^4\text{ yr }\left(\frac{T_0}{100\text{ K}}\right){\left(\frac{M}{M_{\mathrm{J}}}\right)}^{-2},$$

$$t_{\mathrm{Bondi}}\equiv \frac{M}{\dot{M}}=\frac{c_{\mathrm{s}}^3}{4\text{uppi}{G}^2\rho M},$$

$$t_{\mathrm{Bondi}}=\frac{m}{\dot{m}}=\frac{\Omega c_{\mathrm{s}}^3}{4\text{uppi}{G}^2{\rho }_0{M}_{\star }}\frac{{\rho }_0}{\rho }{m}^{-1}=t_0{h}^3\frac{{\rho }_0}{\rho }{m}^{-1}$$

$$t_0\equiv \frac{{\Omega }^2}{4\text{uppi}G{\rho }_0}$$

$$\frac{\Sigma m^2}{x^3}t={\Sigma }_0{w}^2.$$

$$\frac{{\rho }_0}{\rho }=\frac{{\Sigma }_0}{\Sigma }=\frac{m^2}{x^3{w}^2}t\approx \frac{m^2}{h^5}t$$

$$m_{\mathrm{final}}^{\mathrm{simple}}\equiv \frac{M_{\mathrm{final}}^{\mathrm{simple}}}{M_{\star }}\sim \frac{h^2}{t_0}={\left(\frac{H}{a}\right)}^2\frac{4\text{uppi}G{\rho }_0}{{\Omega }^2}$$

$$M_{\mathrm{final}}^{\mathrm{simple}}\approx 0.4M_{\mathrm{J}}{\left(\frac{a}{10\text{ AU}}\right)}^{11/14},$$

$$\frac{{\rho }_0}{\rho }=\{\begin{array}{l}\begin{array}{rlr}& m^4{h}^{-549/49}{t}^{39/49}& t<t_1\text{(inviscid)}\\ & m^4{h}^{-391/35}{\alpha }^{-4/35}{t}^{5/7}& t_1<t<t_2\text{(partial)}\\ & m^4{h}^{-53/5}{\alpha }^{-7/5}& t>t_2\text{(full viscous)}\end{array}\end{array}$$

$$t(m)=\frac{m}{\dot{m}}=\{\begin{array}{l}\begin{array}{rlr}& t_0^{49/10}{h}^{-201/5}{m}^{147/10}& m<m_1\text{(inviscid)}\\ & t_0^{7/2}{h}^{-143/5}{\alpha }^{-2/5}{m}^{21/2}& m_1<m<m_2\text{(partial)}\\ & t_0{h}^{-38/5}{\alpha }^{-7/5}{m}^3& m>m_2\text{(full viscous)}\end{array}\end{array}$$

$$\frac{M_{\mathrm{final}}^{\mathrm{inviscid}}}{M_{\star }}={\left(\Omega t_{\mathrm{disc}}\right)}^{10/147}{\left(\frac{H}{a}\right)}^{134/49}{\left(\frac{4\text{uppi}G{\rho }_0}{{\Omega }^2}\right)}^{1/3},$$

$$\frac{\Sigma m^2}{h^3}t={\Sigma }_1{w}_1^2.$$

$$w_1={\left(\frac{{\Sigma }_1}{\Sigma }\frac{h^{11}}{m^2}\right)}^{1/5}$$

$$\frac{{\Sigma }_1}{\Sigma }=\{\begin{array}{ll}{m}^2{h}^{-37/7}{t}^{5/7}& t<t_1\\ m^2{h}^{-5}{\alpha }^{-1}& t>t_1\end{array}$$

$$w_1=\{\begin{array}{ll}(h^{8}t{)}^{1/7}& t<t_1\\ (h^6/\alpha {)}^{1/5}& t>t_1\end{array},$$

$$\frac{{\Sigma }_1{m}^2}{w_1^3}t={\Sigma }_2{w}_2^2$$

$$w_2={\left(\frac{{\Sigma }_2}{{\Sigma }_1}\frac{w_1^8{h}^3}{m^2}\right)}^{1/5},$$

$$\frac{{\Sigma }_2}{{\Sigma }_1}=\frac{m^2{t}^{5/7}}{h^{6/7}{w}_1^{31/7}}=\{\begin{array}{ll}{m}^2{h}^{-290/49}{t}^{4/49}& t<t_1\\ m^2{h}^{-216/35}{\alpha }^{31/35}{t}^{5/7}& t_1<t<t_2\\ m^2{h}^{-28/5}{\alpha }^{-2/5}& t>t_2,\end{array}$$

$$\frac{{\Sigma }_2}{\Sigma }=\frac{{\Sigma }_2}{{\Sigma }_1}\frac{{\Sigma }_1}{\Sigma }=\{\begin{array}{ll}{m}^4{h}^{-549/49}{t}^{39/49}& t<t_1\\ m^4{h}^{-391/35}{\alpha }^{-4/35}{t}^{5/7}& t_1<t<t_2\\ m^4{h}^{-53/5}{\alpha }^{-7/5}& t>t_2.\end{array}$$

Full text

††pubyear: 2019††pagerange: The End of Runaway: How Gap Opening Limits the Final Masses of Gas Giants–A

1 Introduction

22footnotetext: 51 Pegasi b Fellow

Perhaps the leading theory for the formation of giant planets is the core accretion model (Perri & Cameron, 1974; Harris, 1978; Mizuno et al., 1978; Mizuno, 1980; Stevenson, 1982; Bodenheimer & Pollack, 1986; Pollack et al., 1996). According to this model, gas giants assemble from the bottom up: a rocky or icy core several times the mass of the Earth accretes gas from the circumstellar disc in which it is embedded. The process must complete within a few million years, the observed disc lifetime (Mamajek, 2009; Williams & Cieza, 2011; Alexander et al., 2014). Initially, gas accretion is regulated by Kelvin–Helmholtz cooling, i.e., how long it takes the nascent atmosphere to radiate away its gravitational potential energy (Piso & Youdin, 2014). When still lighter than the core, the atmosphere grows (read: doubles in mass) on a timescale that increases as roughly the square of its mass (Lee & Chiang, 2015). The situation changes dramatically once the atmosphere exceeds the core in mass. Then the self-gravity of the atmosphere causes it to cool—and therefore accrete mass from the surrounding disc—over ever shorter timescales. Gas accretion is now in the “runaway” phase (Bodenheimer & Pollack, 1986; Pollack et al., 1996; Ikoma et al., 2000; Lee et al., 2014; Piso & Youdin, 2014; Lee & Chiang, 2015; Piso et al., 2015; Berardo et al., 2017, and references therein).

How does runaway accretion stop? In other words, what determines the final mass of a gas giant? The total gas disc mass sets a strict upper bound, but for most discs this is not constraining, as inside an orbital radius of a few hundred AU, discs are estimated to have 10–100 times more gas than Jupiter (see, e.g., Fig. 10 of Tripathi et al., 2017). Even these gas mass estimates are often only lower bounds because they are based on masses inferred from dust emission and gas-to-dust ratios that are assumed to be solar; in reality, inward radial drift of solid particles leads to higher gas-to-dust abundances at large disc radii (e.g., Andrews et al., 2012; Powell et al., 2017). The problem of stopping runaway is made more acute by recent indications that the giant planet mass function is bottom-heavy, with the number of planets per decade in mass decreasing inversely with mass (Nielsen et al., 2019; Wagner et al., 2019). Understanding the endgame of giant planet formation—in particular the accretion luminosities of nascent giants (e.g., Berardo et al., 2017)—is critical for motivating and interpreting observational campaigns to directly image protoplanets.

Gap opening by planets (Goldreich & Tremaine, 1980) has been considered a key mechanism halting a planet’s growth (Lin & Papaloizou, 1993; Bryden et al., 1999; Kley, 1999; Lubow et al., 1999; D’Angelo et al., 2003; Tanigawa & Ikoma, 2007; Lissauer et al., 2009; Machida et al., 2010; Tanigawa & Tanaka, 2016). Gravitational torques exerted by the planet repel gas away from its orbit, carving a gap (annular cavity) in the disc and diminishing the supply of gas to the planet. At the same time, viscous torques intrinsic to the disc diffuse material back into the gap. The essential quantity to compute is the gap depth, i.e., the extent to which the gas density in the planet’s vicinity is lowered relative to the background disc value (Duffell & MacFadyen, 2013; Fung et al., 2014; Kanagawa et al., 2015a).111Kratter et al. (2010; see also Thorngren et al. 2016 and Rosenthal & Murray-Clay 2018) suggested in their subsection 5.2.1 an ad-hoc gap starvation mass based on gap widths. Their discussion motivates a more thorough analysis of accretion rates inside gaps, which is what our study provides. Ginzburg & Sari (2018, hereafter GS18) calculate both gap depths and radial widths by considering in detail how planet-driven waves dissipate within gaps. In addition to providing physical justifications for previous numerical results, they uncover new gap scaling behaviours appropriate to low disc viscosities.

Here we build upon these advancements to re-visit the role of gap opening in determining the final masses of gas giants. We focus on how gaps open in discs having low viscosities ($\alpha\leq 10^{-3}$, where $\alpha$ is the turbulent Mach number introduced by Shakura & Sunyaev 1973), as motivated by Atacama Large Millimeter Array (ALMA) observations that not only reveal gaps in discs but also point to nearly inviscid environments (e.g., Pinte et al., 2016; Flaherty et al., 2017; Dong et al., 2018; Hartmann & Bae, 2018; Zhang et al., 2018). Another novel feature of our analysis is that we account for how gap depths vary (deepen) with time. Our approach is analytic and approximate, to gain intuition and guide future, more precise numerical studies.

The outline of the rest of this paper is as follows. After setting down in the next subsections our model assumptions, we describe in Section 2 how runaway growth proceeds unchecked when gap opening is ignored. The tables turn in Section 3 where we account for gap clearing and provide closed-form expressions for the final masses of gas giants. In Section 4 we summarize and give an outlook.

1.1 Accretion onto sub-thermal planets:

Spherical symmetry at the Bondi radius

We restrict our analysis to planets with masses $M$ below the thermal mass, defined here as that for which the Hill radius $R_{\rm H}$ equals the gas disc’s scale height $H$:

[TABLE]

where $M_{\star}$ is the host star mass and $a$ is the orbital radius. Historically, only planets above $M_{\rm thermal}$ were thought to open gaps and stop growing (e.g., Lin & Papaloizou, 1993; Bryden et al., 1999), the rationale being that planet-driven density waves had to be non-linear in order to shock, dissipate, and transfer angular momentum to the disc. However, Goodman & Rafikov (2001) explained that even linear waves gradually steepen and eventually shock as they propagate away from the planet. Modern numerical calculations (Duffell & MacFadyen, 2013; Fung & Chiang, 2017) confirm that sub-thermal planets are indeed capable of carving out deep gaps (if $\alpha$ is sufficiently low) and can therefore potentially stop growing before reaching $M_{\rm thermal}$.

Sub-thermal planets are easier to model than super-thermal planets in the following sense. In general, the outer radius of a planet embedded in a disc is the smaller of the Hill radius or the Bondi radius

[TABLE]

where $c_{\rm s}$ is the nebular sound speed and $G$ is the gravitational constant. For $M<M_{\rm thermal}$, the length scales order as $R_{\rm B}<R_{\rm H}<H$. This hierarchy simplifies planetary accretion in several ways. First, because sub-thermal planetary atmospheres are on scale $R_{\rm B}<H$, anisotropies in disc density (differences in density between the vertical and in-plane directions) can be ignored; the atmospheres will have a roughly spherical symmetry (Rafikov, 2006; Piso & Youdin, 2014). Second, most of a planet’s repulsive torque is carried by waves launched at radial distances \sim$$H away from the planet (Goldreich & Tremaine, 1980) and deposited at still larger distances (Goodman & Rafikov, 2001). Thus, the gas nearest the planet (at the bottom of the gap, which even for a sub-thermal planet can be deep) has an approximately flat radial density profile, up to distances of at least $H>R_{\rm B}$ (e.g., Kanagawa et al., 2015a). Finally, the fact that $R_{\rm B}<R_{\rm H}$ means that the Keplerian shear velocity at the planet’s outer Bondi radius is subsonic and smaller than the circumplanetary orbital velocity, supporting the classical Bondi picture which neglects angular momentum and takes the accretion velocity to be sonic at $R_{\rm B}$. In other words, a circumplanetary disc does not necessarily form at $R_{\rm B}$, and accretion on that scale may be spherically symmetric (rotation becomes significant on scales $<R_{\rm B}$; see our Section 4.1).

Our restriction to sub-thermal masses and assumption of spherical symmetry mean that we are working in a regime complementary to that considered by Szulágyi et al. (2014), who address the problem of limiting giant planet masses in the context of circumplanetary discs (see also Szulágyi & Mordasini, 2017).

1.2 Background disc model

For most of our calculations we focus on the nominal case of a nascent giant planet located at an orbital separation of $a=10\textrm{ AU}$ from a solar mass star ($M_{\star}=M_{\odot}$). The gas surface density of our disc when unperturbed (i.e., with no planet) resembles that of the minimum-mass solar nebula (MMSN)222 By construction, the MMSN contains just enough gas to form a Jupiter-mass planet. This total mass constraint does not appear in our calculations, as we wish to explore other limitations to runaway growth. (Weidenschilling, 1977; Hayashi, 1981; Chiang & Youdin, 2010):

[TABLE]

The disc temperature is similar to that derived by Chiang & Goldreich (1997):

[TABLE]

The sound speed is given by $c_{\rm s}=\sqrt{\gamma k_{\rm B}T_{0}/\mu}$, where $k_{\rm B}$ is the Boltzmann constant. We take an adiabatic index $\gamma=7/5$ and a molecular weight $\mu=2\textrm{ amu}$ appropriate for molecular hydrogen. The resulting volumetric gas density at the disc midplane is

[TABLE]

(e.g., Frank et al., 2002) with $H=\Omega^{-1}\sqrt{k_{\rm B}T_{0}/\mu}$ the disc’s scale height and $\Omega=(GM_{\star}/a^{3})^{1/2}$ the orbital angular velocity. For our parameter choices, the disc’s aspect ratio equals

[TABLE]

Our nominal disc is thicker (hotter) than in other models (e.g., D’Alessio et al. 1998 compute $H/a\approx 0.06$ in their Fig. 6), possibly causing us to overestimate, by factors of a few, both the range of planet masses that can remain sub-thermal (equation 1) and final planet masses (Section 3.2.3, in particular equation 19). We allow ourselves this latitude as there are other model uncertainties that are also order-unity if not larger.

We assume a nominal disc lifetime of $t_{\rm disc}=3\textrm{ Myr}$, during which the background density $\Sigma_{0}$ remains constant, and after which it vanishes. This time covers only the runaway phase of planetary gas accretion and its aftermath; the assembly of the underlying rocky core, and the initial pre-runaway phase of gas accretion (discussed briefly in Section 2.1) are not covered. What fraction of time rocky core assembly takes of the total disc lifetime (spanning 1–10 Myr; e.g., Mamajek 2009; Pfalzner et al. 2014) is unknown; as such, we may be overestimating the duration of the runaway/post-runaway phase when setting it equal to $t_{\rm disc}=3\textrm{ Myr}$. Note that merely decreasing $t_{\rm disc}$ in order to stop runaway at some desired mass would require unreasonable fine-tuning, as growth timescales become extremely short with increasing mass (it is, after all, called “runaway” for a reason; see, e.g., equation 9). We will see in our full theory of accretion within gaps that final planet masses are insensitive to changes in $t_{\rm disc}$, changing only by several percent when $t_{\rm disc}$ varies by factors of a few (equation 19).

2 Gapless accretion

We first consider how planets accrete gas from discs neglecting gap clearing. This will demonstrate the problem of runaway accretion in both thermodynamic (Section 2.1) and hydrodynamic (Section 2.2) senses.

2.1 Runaway cooling

Pre-runaway gas accretion onto planetary cores has been studied extensively, both in the context of gas giants and short-period sub-Neptunes (Pollack et al., 1996; Ikoma et al., 2000; Papaloizou & Nelson, 2005; Rafikov, 2006; Piso & Youdin, 2014; Lee et al., 2014; Lee & Chiang, 2015; Piso et al., 2015; Ginzburg et al., 2016). These studies find that the planet’s atmosphere is divided into an interior convective region containing most of the mass, and an outer radiative and nearly isothermal envelope which extends to $R_{\rm B}$. Conditions at the radiative–convective boundary dictate the planet’s cooling rate which is the bottleneck for pre-runaway accretion. The atmosphere grows (doubles) on a timescale that increases with its mass as long as it is much lighter than the core. Once the atmosphere and core become comparable in mass, the growth time decreases, initiating the runaway phase.

In this subsection we extend the analytical scaling relations of previous studies to the runaway phase, where the core mass can be neglected. The main difference is that the gravity field is dominated by the gas itself, and not by the core, transforming the atmosphere’s density profile in the convective region into a polytrope akin to those in stars. We keep the following analysis approximate, as we shall see that the cooling rate derived in this subsection is not the ultimate bottleneck for accretion.

The density in the convective region scales with temperature as $\rho(T)/\rho_{\rm rcb}=(T/T_{0})^{1/(\gamma-1)}$, with $\rho_{\rm rcb}$ denoting the density at the radiative–convective boundary (RCB), and $T_{0}\sim T_{\rm rcb}$ since the radiative envelope is isothermal. From hydrostatic equilibrium, the temperature at the planet’s centre is given by $k_{\rm B}T_{\rm c}\sim GM\mu/R_{\rm rcb}$, with $M$ denoting the planet’s mass (assumed dominated by the convective gas layer), and $R_{\rm rcb}$ the radius of the convective polytrope. By substituting $\rho(T_{\rm c})\sim M/R_{\rm rcb}^{3}$ for the central density of the polytrope we find

[TABLE]

The luminosity is calculated by applying the diffusion equation to the radiative layer $L\sim\sigma T_{0}^{4}R_{\rm B}/(\kappa\rho_{\rm rcb})$, where $\sigma$ is the Stefan–Boltzmann constant and $\kappa$ is the opacity at the RCB (see, e.g., Section 2.2 of Ginzburg et al. 2016; we have here omitted order-unity coefficients). The gravitational energy that the planet has to radiate away is $E\sim GM^{2}/R_{\rm rcb}$. The Kelvin–Helmholtz cooling time is therefore

[TABLE]

where we substituted $\rho_{\rm rcb}$ from equation (7). The radius $R_{\rm rcb}$ is related to $R_{\rm B}$ by a logarithmic factor that stems from the exponentially declining density profile in the outer isothermal envelope (Piso & Youdin, 2014; Ginzburg et al., 2016). We omit this logarithmic factor and approximate $R_{\rm rcb}\sim R_{\rm B}$, simplifying equation (8):

[TABLE]

where $M_{\rm J}$ is the mass of Jupiter and we have assumed for simplicity a constant $\kappa=0.1\textrm{ cm}^{2}\textrm{ g}^{-1}$ (from dust; see Appendix C in Piso et al. 2015).

Equation (9) demonstrates how decreasing cooling times lead to runaway growth on a timescale much shorter than the disc’s lifetime. One interesting feature of runaway cooling is its near-insensitivity to the nebula’s density. Pre-runaway cooling is characterized by a similar insensitivity (Lee & Chiang, 2016).

We emphasize that equation (9), which indicates that the cooling/growth timescale decreases with increasing envelope mass, is valid only for atmospheres that outweigh their cores, $M>M_{\rm core}$. By contrast, lighter atmospheres are accreted on timescales that increase with increasing gas mass (Lee & Chiang, 2015). We conclude that atmospheres take the longest time to double their mass when $M\sim M_{\rm core}$. By setting the cooling/growth time in equation (9) to the gas disc lifetime of \sim$$3\times 10^{6} yr, we see that only planets above a critical mass of $M\sim 10$–$20\,M_{\oplus}$ achieve runaway, for our nominal $T_{0}=100\textrm{ K}$ (see Piso & Youdin 2014 and Piso et al. 2015 for a more comprehensive analysis). Our derivation is indicated graphically in Fig. 1. The intersection (just outside and to the left of the plotted range) of the cooling timescale (dashed green line) and the gas disc lifetime (dotted horizontal black line) provides a rough estimate for the critical planet mass (and by extension core mass $M_{\rm core}\sim M$).

2.2 Bondi accretion

As is apparent from equation (9) and Fig. 1, the cooling timescale becomes exceedingly short during runaway. At some stage, the nebula will not be able to replenish gas at a sufficient rate to the planet’s outer boundary ($R_{\rm B}$) to maintain this rapid growth. Accretion ceases to be cooling/thermodynamically limited and becomes instead hydrodynamically limited.

The maximum rate at which gas can be brought to the Bondi radius is given by the free-fall velocity, which is also the sound speed at that radius. The corresponding mass accretion rate is the Bondi rate, $\dot{M}=4\uppi R_{\rm B}^{2}c_{\rm s}\rho$, where $\rho$ is the ambient gas density. The fastest, Bondi-limited growth time is then

[TABLE]

which, like $t_{\rm cool}$, decreases with increasing $M$ and therefore also implies runaway behaviour.

In an attempt to calculate the final masses of gas giants, Tanigawa & Tanaka (2016) utilize an accretion formula that is drawn empirically from two-dimensional simulations by Tanigawa & Watanabe (2002) and scales as $\dot{M}\propto M^{4/3}$ (recently this scaling has been explained analytically by Lee, 2019, in terms of isothermal accretion shocks). These simulations account for hydrodynamical flows that become increasingly complex and rotation-dominated as the planet mass increases. However, in the limit $R_{\rm B}\ll R_{\rm H}$, i.e., $M\ll M_{\rm thermal}$, the geometry and flow are expected to be simple and to match the results of classical Bondi accretion: the envelope boundary at $R_{\rm B}<H$ should be more-or-less spherically symmetric, and the Keplerian shear velocity at $R_{\rm B}$ should be negligible as it is smaller than the sound speed $c_{\rm s}$ (see our Section 1.1). Indeed, as can be appreciated from Fig. 1 of Tanigawa & Tanaka (2016), the numerically simulated three-dimensional accretion rates at the lowest masses match the Bondi scaling $\dot{M}\propto M^{2}$ better than $\dot{M}\propto M^{4/3}$. As we assume in this paper that $M<M_{\rm thermal}$—an assumption that becomes better justified at larger orbital distances—we use the Bondi result (see also the discussion in Sections 1.1 and 4.1).

If we ignore gap clearing, the surrounding density is given by the unperturbed nebular value $\rho=\rho_{0}\approx 3\times 10^{-12}\textrm{ g cm}^{-3}$ at our nominal disc radius of 10 AU (Section 1.2). By comparing equations (9) and (10) using this background (gapless) $\rho$, we find that the Bondi rate limits gas accretion only at unrealistically high masses—above $200\,M_{\rm J}$ for our nominal 100 K disc (and even then, accretion still proceeds in a runaway manner, as explained above). In the next section, we consider how gaps lower $\rho$ and stop runaway at smaller masses.

3 Gap Clearing

Equation (10) indicates that the Bondi-limited growth time depends on the planet’s ambient density $\rho$ (in contrast to the Kelvin–Helmholtz cooling time, which is insensitive to $\rho$; see Section 2.1). Therefore, once a planet opens a gap around its orbit, the decreasing density prolongs the Bondi accretion timescale, conceivably beyond the gas disc lifetime $t_{\rm disc}$.

By balancing gravitational and viscous torques, Fung et al. (2014) derived an analytical scaling for the depth of a gap with respect to the background: $\rho/\rho_{0}\sim\alpha h^{5}/m^{2}$, where $h\equiv H/a$ is the disc’s aspect ratio and $m\equiv M/M_{\star}$ is the planet-to-star mass ratio (see also Duffell & MacFadyen 2013 who found the same scaling empirically). The Shakura & Sunyaev (1973) $\alpha$ parametrizes the kinematic viscosity $\nu\equiv\alpha c_{\rm s}H$. The scaling compares well against some multidimensional hydrodynamical simulations (e.g., Kanagawa et al., 2015b; Fung & Chiang, 2016) and demonstrates that even low-mass planets can open deep gaps, if the viscosity is low enough. By substituting the gap density scaling into equation (10), we find that $t_{\rm Bondi}\propto m$, implying that gap clearing stops accretion from running away. Eventually, the planet will grow to its final mass, for which $t_{\rm Bondi}$ becomes as long as $t_{\rm disc}$. Tanigawa & Tanaka (2016) also333We reiterate that Tanigawa & Tanaka (2016) use an empirical accretion rate $\dot{M}\propto M^{4/3}$, whereas we adopt the Bondi rate $\dot{M}\propto M^{2}$. Our choice is justified in Section 2.2. used this analytical depth scaling to argue that gap opening can limit the final mass of gas giants to about $10\,M_{\rm J}$, a result that depends on the viscosity parameter $\alpha$.

A new analysis by GS18 indicates that the above scaling for $\rho/\rho_{0}$ is altered for very low viscosities. The reason is that the total gravitational torque that the planet exerts is no longer dominated by interaction with gas at a distance $H$ from the planet (as assumed in the derivation of Fung et al., 2014). Once the density there drops sufficiently, the total torque becomes dominated by interaction with gas farther away from the planet, leading to an even wider and deeper gap. At the same time, GS18 also find that the time to open a gap in a low-viscosity disc is long—so long that planets might not have enough time to fully clear their gaps as they grow.

Here we calculate the final mass of a gap-opening planet by taking into account the updated gap profiles from GS18 and the temporal evolution of the gap clearing process. This new analysis leads to qualitatively different results.

3.1 Simplified solution

We first develop intuition by deriving in a simplistic way the final mass of a gap-opening planet in the inviscid ($\alpha=0$) limit. Our aim here is purely pedagogical; the simplified derivation is inaccurate but captures some of the essential arguments. We present it merely as a guide for our more careful analysis in Section 3.2.

We adopt the GS18 dimensionless notation $G=M_{\star}=a=1$. In these units time $t$ is measured in units of the inverse Kepler frequency $\Omega^{-1}=1$. As before, $m\equiv M/M_{\star}$ and $h\equiv H/a$. We rewrite the (now dimensionless) Bondi accretion timescale given by equation (10) as

[TABLE]

where the dimensionless time

[TABLE]

is a function of the background density $\rho_{0}$ and the distance from the star (through $\rho_{0}$ and $\Omega$; $t_{0}$ is related to the unperturbed nebula’s Toomre stability parameter $Q\gtrsim 1$ by an order-unity coefficient).

We now consider how $\rho_{0}/\rho$, i.e., the depth of the gap cleared by the planet, evolves with time. The planet exerts a repulsive torque $\Sigma m^{2}/x^{3}$, where $x$ denotes the radial distance away from the planet where the strongest resonant interactions with the disc occur, inside the gap having characteristic surface density $\Sigma$ (Goldreich & Tremaine, 1980; Goodman & Rafikov, 2001; Fung et al., 2014; Ginzburg & Sari, 2018). The torque displaces gas by means of density waves that travel a radial distance $w\gtrsim x$ before damping and releasing the angular momentum they carry to the disc (Goodman & Rafikov, 2001; Ginzburg & Sari, 2018). Just outside of $w$, the surface density returns to its unperturbed444GS18 use a different notation, with $\Sigma_{0}$ denoting the density at the bottom of the gap and $\Sigma_{\infty}$ the unperturbed density at its top. value $\Sigma_{0}$. The angular momentum required to displace an annulus of width $w\ll 1$ over its own width is $\Sigma_{0}w^{2}$. In the absence of viscosity, the depth and width of the gap grow freely with time as the planet deposits ever more angular momentum to the disc:

[TABLE]

In later sections and the Appendix, we consider the physics underlying $x$ and $w$, in particular how both increase with time, and the relation between $w$ and $\Sigma$. These time dependencies, which are necessary for a self-consistent analysis, are ignored here for simplicity; we assume for now that $x\sim w\sim h$, following previous crude estimates (e.g., Duffell & MacFadyen, 2013; Fung et al., 2014). Then the gap contrast increases as

[TABLE]

from a minimum value of $\rho_{0}/\rho=1$.

Inserting equation (14) into equation (11), we find that the planet’s mass grows logarithmically with time (since $\dot{m}\propto\rho\propto t^{-1}$) and reaches a final mass of

[TABLE]

modulo a logarithmic factor that depends on the gas disc’s lifetime $t_{\rm disc}$.

In the next Section 3.2 we derive a more accurate expression for $M_{\rm final}$ which takes into account the gap’s detailed structure ($x,w\neq h$) and the disc’s finite viscosity. Because $w$ increases slowly, we will find that the gap contrast $\rho_{0}/\rho$ grows sub-linearly with time. This leads to a weak power-law (instead of logarithmic) growth for $m(t)$; see equation (18) and Fig. 1, and equation (19) for the dependence of the final mass on $t_{\rm disc}$.

Nevertheless, the simplified expression in equation (15) exhibits many of the same features of our more carefully derived result in Section 3.2. The final mass depends weakly on $t_{\rm disc}$, has essentially no dependence on $\alpha$ (by construction in this simplified solution), and depends explicitly on $\rho_{0}$. All these features contrast with those of previous estimates of the final mass. For our nominal disc, equation (15) sets a final mass scale of

[TABLE]

differing only by a factor of $\sim$2 from our more accurate equation (19). The scaling of the final mass with distance from the star is also similar (11/14 versus 3/4; see Section 3.2.4). One difference is the scaling with $\rho_{0}$, which is weaker in the full calculation.

3.2 Full solution

3.2.1 Gap depletion

We focus on low but non-zero viscosities for which planet gaps are so deep that the total torque is no longer dominated by interaction with gas at a single radial distance $x\sim h$ away from the planet (as assumed by Goldreich & Tremaine 1980 and Fung et al. 2014). Instead, for $\alpha<m^{5}/h^{14}$, the torque peaks at two locations: one at $h$, and another farther away (GS18). Each interaction contributes a factor of $m^{2}$ to the density contrast, resulting in a “two-step” $m^{4}$ scaling. This is an asymptotic result that appears supported by hydrodynamical simulations by Zhang et al. (2018, see in particular the left column of their Fig. 2). Gaps in discs with higher viscosities obey the classical single-step $m^{2}$ scaling (Duffell & MacFadyen, 2013; Fung et al., 2014). We will assess in Section 3.2.4 the extent to which our assumption that $\alpha<m^{5}/h^{14}$ is satisfied, after we compute final planet masses.

Gaps evolve toward an equilibrium depth that is set by a balance between planetary torques which push gas out, and viscous torques which ooze material back in. While GS18 calculate this equilibrium and the time required to reach it, they do not provide gap profiles at earlier times. In the Appendix we complete the calculation of GS18 and derive pre-equilibrium, time-varying gap densities. The density $\rho$ at the bottom of the gap (where the planet resides) relative to the background density $\rho_{0}$ is given by

[TABLE]

where $t$ is the dimensionless time, $t_{1}=(h^{2}/\alpha^{7})^{1/5}$, $t_{2}=(h^{4}/\alpha^{9})^{1/5}$, and it is understood that the equation applies only for times for which $\rho_{0}/\rho>1$ (this condition turns out to be satisfied for our parameters because Kelvin–Helmholtz cooling is long enough for the planet to clear a significant gap; see Fig. 2).

Equation (17) is the more accurate version of the simplified equation (14). Three gap-clearing stages are now delineated. In the first “inviscid” stage, the planetary torque overwhelms the viscous torque, and so the evolution does not depend on $\alpha$. This is the regime considered in the previous simplified analysis of Section 3.1; whereas before we found $\rho_{0}/\rho\propto m^{2}/h^{5}$ following a one-step scaling, we now find $\rho_{0}/\rho\propto m^{4}/h^{549/49}$ reflecting the two-step nature of low-viscosity, deep gaps. The scaling with time is now slightly sub-linear because we have accounted for how the gap width $w$ slowly increases above $h$. At intermediate times $t_{1}<t<t_{2}$, the first density step (closer to the planet) reaches viscous equilibrium and saturates—hence the dependence on $\alpha$—while the second step continues to steepen. In the final stage $t>t_{2}$, both density steps saturate; this is the GS18 equilibrium result.

3.2.2 Growth rates

We substitute $\rho_{0}/\rho(t)$ from equation (17) into equation (11) to obtain planet growth timescales that account for the gap’s temporal evolution:

[TABLE]

with $m_{1}=t_{0}^{-1/3}h^{58/21}\alpha^{-2/21}$ and $m_{2}=t_{0}^{-1/3}h^{14/5}\alpha^{-2/15}$ marking the transitions to first-step and second-step viscous equilibrium at times $t_{1}$ and $t_{2}$, respectively. As a reminder, the doubling time $t$ is measured in units of $\Omega^{-1}$, $m\equiv M/M_{\star}$, $h\equiv H/a$, and $t_{0}$ is defined in equation (12).

We plot in Fig. 1 the planet growth timescales $t(M)=M/\dot{M}$ (with units restored) for different viscosities $\alpha$. At the smallest masses (earliest times), accretion is described by equation (9): it is cooling-limited and proceeds in runaway fashion (dashed green line). Runaway halts when the planet opens a gap. In the first stage of Bondi-limited accretion within a gap (equation 18, $m<m_{1}$, solid blue line), the gap grows inviscidly since growth timescales are too short for equilibrium to be achieved between planetary torques and the initially weaker viscous torque. Departure from the first inviscid stage occurs at different masses $m_{1}(\alpha)$. In the subsequent second ($m_{1}<m<m_{2}$, dot-dashed red line, steep segment) and third ($m>m_{2}$, less steep segment) stages, viscosity competes more effectively and growth timescales are long enough for the gap to reach partial and finally full viscous equilibrium.

The evolution of the gap depth $\rho_{0}/\rho(M)$ as the planet grows is shown in Fig. 2, obtained by solving equations (11) and (17). Even during the earliest cooling-limited growth phase, a deep gap having $\rho_{0}/\rho\approx 10$–100 forms. At low masses ($m<m_{1}$), the gap deepens on the same timescale as the planet’s mass-doubling time (given first by cooling and then by the inviscid limit); at high masses ($m>m_{2}$), the gap achieves viscous equilibrium (GS18). By the time the planet is done forming (when the growth time is as long as $t_{\rm disc}$, i.e., when a given blue-red curve intersects with the black line), the gap depth $\rho_{0}/\rho\sim 10^{5}$. Note how this final depletion factor can be estimated independently of any gap formation theory: the dotted black line in Fig. 2 is merely given by equating the Bondi-limited growth time in equation (10) to $t_{\rm disc}$. This calculation implies that, for our nominal parameters at 10 AU, a Jupiter-mass planet undergoes its last doubling in mass in a disc whose gas density is depleted locally relative to the MMSN by \sim$$10^{5}.

3.2.3 Final mass

Fig. 1 demonstrates that gap clearing lengthens the Bondi accretion time above the Kelvin–Helmholtz cooling time and ultimately limits a giant planet’s growth. Gap density contrasts increase super-linearly with planet mass ($\rho_{0}/\rho$ scales as $m^{4}$ in the low-viscosity limit, and as $m^{2}$ at higher viscosities); this dependence puts an end to runaway accretion by increasing the mass doubling time with increasing mass. In the inviscid limit, $t$ scales nearly as $M^{15}$ (equation 18; cf. the exponential $t(M)$ in our simplified derivation of Section 3.1). The inviscid doubling time increases so rapidly that it quickly exceeds the disc lifetime $t_{\rm disc}$ (horizontal black line in Fig. 1), finalizing the planet’s mass at

[TABLE]

where $\Omega$ is the planet’s orbital frequency, $H/a$ is the disc’s aspect ratio, and $\rho_{0}$ is its unperturbed background density. For our nominal disc parameters at 10 AU, $M_{\rm final}^{\rm inviscid}\approx 1M_{\rm J}$. Fig. 1 shows that this inviscid final mass applies for $\alpha\lesssim 10^{-4}$, and that the final mass increases to \sim$$2.5M_{\rm J} for $\alpha=10^{-3}$. That equation (19) for the final mass and equation (1) for the thermal mass scale similarly with $H/a$ implies that our assumption that masses stay sub-thermal tends to hold regardless of the disc’s temperature profile.

Equation (19) is the more accurate version of equation (15). Among the features of the final mass as given by equation (19) are that it does not depend on the (uncertain) disc viscosity and depends only weakly on the disc lifetime. While the simpler version reproduces these and other qualitative features of the full calculation, the quantitative agreement between the two is partly coincidental.

3.2.4 Dependence on orbital distance

We present the final planet mass as a function of distance from the star in Fig. 3. The final mass is computed in one of two ways, either in the inviscid limit (equation 19) or when the gap is fully viscously equilibrated (equation 18 for $m>m_{2}$, with $t$ set to $\Omega t_{\rm disc}$). We disregard for simplicity the intermediate stage of partial equilibrium since it spans a small mass range and hardly affects the final mass (see Fig. 1). Fig. 3 indicates that the inviscid limit applies at large separations (e.g., $\gtrsim 5$ AU for $\alpha=10^{-4}$) and the viscous regime applies at small separations; far from the star, diffusion times across gaps are so long that the disc behaves as if it were inviscid.

According to our nominal disc model (Section 1.2), $M_{\rm final}^{\rm inviscid}\propto a^{0.75}$. By comparison, in viscous equilibrium, $M_{\rm final}\propto a^{0.30}$. Also plotted for reference in Fig. 3 is the thermal mass $M_{\rm thermal}\propto a^{6/7}\approx a^{0.86}$ (equation 1). Many of our simplifying approximations (e.g., neglect of rotation; see Section 4.1) break down for $M>M_{\rm thermal}$, and so we see that our final mass estimates are safest at large distances ($a\gtrsim 3$ AU for $\alpha=10^{-4}$ and $a\gtrsim 20$ AU for $\alpha=10^{-3}$), and then only marginally so. The ratio of the thermal mass to the inviscid final mass is insensitive to the disc’s temperature profile, as both vary similarly with $h$.

Finally, we check whether the low-viscosity condition $\alpha<m^{5}/h^{14}$ assumed by our derivation (Section 3.2.1) is satisfied. At our nominal 10 AU, $h\approx 0.08$ and $m_{\rm final}^{\rm inviscid}\approx 0.001$, and so $m^{5}/h^{14}\approx 2$, validating our assumption in the case of the inviscid final mass. Although pre-final masses are smaller and may in principle violate the low-viscosity condition, in practice this is not an issue because the growth time is an extremely steep function of mass (Fig. 1). Furthermore, for the inviscid final mass, $m^{5}/h^{14}$ scales only weakly with distance, as $a^{-0.25}$, implying that the entire inviscid curve plotted in Fig. 3 is computed self-consistently. This implies that the viscous curves in Fig. 3 are also self-consistent, as they yield more massive planets for which $m^{5}/h^{14}$ is only larger. We note that $m^{5}/h^{14}$ is not particularly sensitive to the disc’s temperature profile, as the dependence of $m_{\rm final}^{\rm inviscid}$ on $h$ nearly cancels out the factor of $h^{14}$.

4 Conclusions and discussion

Gap opening by planets is a potential mechanism to stop runaway gas accretion and thereby set final planetary masses. Contrary to some of the older literature (Lin & Papaloizou, 1993; Bryden et al., 1999; Ida & Lin, 2004), this scenario is not restricted to planets exceeding the thermal mass. Sub-thermal planets also carve out gaps, provided their host discs have low viscosities (e.g., Duffell & MacFadyen, 2013; Fung & Chiang, 2017; Ginzburg & Sari, 2018). With this fact in mind, and motivated further by ALMA observations that point to low levels of turbulence in discs, we have revisited the gap starvation hypothesis and calculated post-runaway masses, restricting consideration to sub-thermal planets (whose gas envelopes might still be reliably modelled in 1D) in low-viscosity discs. A new feature of our analysis is an accounting for how gaps may not have enough time to reach equilibrium depletion levels in nearly inviscid environments. We developed a fully time-dependent theory of gap formation and coupled it to planetary growth rates.

We identified a purely inviscid regime in which gap depletion is limited by time and not by viscosity. By “inviscid” we do not mean that the Shakura & Sunyaev (1973) viscosity parameter $\alpha$ is literally zero. Rather, gap formation in all discs, including those with non-zero $\alpha$, initially proceeds as if gas were inviscid. At early times, gravitational torques from planets dominate viscous torques, and gaps deepen and widen freely with no viscous backflow. If $\alpha$ is low enough, gaps never evolve beyond this inviscid phase before the disc dissipates. We derived in this case an explicit expression for the planet’s post-runaway mass—the final “inviscid” mass—that does not depend on $\alpha$. This result is qualitatively different from that of Tanigawa & Tanaka (2016), who assumed viscous equilibrium.

The inviscid limit applies broadly. For $\alpha\leq 10^{-4}$ ($\alpha\leq 10^{-3}$), planets at orbital separations beyond 5 AU (50 AU) that would otherwise undergo runaway Bondi accretion have their growth halted by gaps that clear inviscidly (Fig. 3). For all of the above parameter space, planet masses remain sub-thermal, consistent with our working assumption, albeit marginally so.

At 10 AU, growth stops after a few Myr (when the disc expires) at the final inviscid mass of $1M_{\rm J}$, as given by equation (19). The final inviscid mass scales as $M_{\rm final}^{\rm inviscid}\propto t_{\rm disc}^{10/147}\rho_{0}^{1/3}$, with $t_{\rm disc}$ denoting the disc’s lifetime and $\rho_{0}$ its unperturbed density. The weak sensitivity to $t_{\rm disc}$ and $\rho_{0}$ and the complete lack of dependence on $\alpha$ imply that gas accretion onto giant planets is not sensitive to the details of how gas discs evolve and dissipate. Note that low viscosities do not necessarily conflict with observed disc lifetimes, as modern studies find that discs may have laminar (low $\alpha$) midplanes (e.g., Pinte et al., 2016; Flaherty et al., 2017), with accretion and mass loss restricted to surface layers (e.g., Bai, 2016). This justifies our decoupling of $\alpha$ from $t_{\rm disc}$.

Folding together all the factors that depend on orbital distance $a$ in equation (19), we see that the final inviscid mass for giant planets increases approximately as $a^{0.75}$. Giant planets are expected to be more massive at larger distance where disc aspect ratios are larger and gaps are harder to open (equation 17). The predicted trend of mass with distance can be tested observationally, say with microlensing or direct imaging campaigns. We emphasize that the prediction pertains to giant planet masses and not to planet occurrence rates. Current direct imaging surveys (e.g. Nielsen et al., 2019) indicate that the giant planet frequency declines outside 10 AU, suggesting that the controlling factor for giant planet formation at large distances may not be gas accretion (the subject of this paper) but instead the agglomeration of rocky cores massive enough to undergo runaway (see, e.g., Lin et al., 2018, and references therein). Any observational test must also screen out brown dwarfs, which occur more frequently beyond 10 AU, and which exhibit other demographic differences with giant planets, presumably because the former form by gravitational instability and the latter by core accretion (Nielsen et al., 2019).

We can try, e.g., to compare our predictions against observations of the HR 8799 system. The four known planets are spaced between 14 and 68 AU and appear to have roughly equal (model-dependent) masses of $5-7M_{\rm J}$ (Bowler, 2016). This constant mass seems to fit the shallower viscous curve ($M_{\rm final}\propto a^{0.30}$; see Fig. 3) better than the inviscid limit. For higher viscosities lying outside the inviscid limit, equation (18) and Fig. 1 yield a final mass of several times $M_{\rm J}$, consistent with previous numerical studies (e.g., Lissauer et al., 2009, who studied $\alpha=4\times 10^{-4}$ and $\alpha=4\times 10^{-3}$).

4.1 Unresolved issues

Notwithstanding the improved theory of gap opening from GS18 that we have applied to the problem of stopping runaway, there remain significant uncertainties in our understanding of deep and wide gaps in low-viscosity discs. Some of the theory’s approximations break down when the width of the gap becomes comparable to the orbital radius, as it does in our models. Furthermore, hydrodynamical instabilities that are not accounted for in the 1D theory (e.g., the Rossby wave instability; Lovelace et al. 1999; Li et al. 2000) might in reality limit gap depletion. For more discussion of these and other issues, see Section 5.1 of GS18.

We have assumed that planets, in the act of accreting, do not consume all the material in their vicinity and empty their gaps. That is, we have assumed that whatever gas is locally accreted on the scale of the Bondi radius $R_{\rm B}$ is replenished by radial transport from farther away, at fast enough rates that the gas density at $R_{\rm B}$ is maintained over the planet’s doubling time. Viscous diffusion offers one means of replenishment, and is more effective as the gap transitions out of its inviscid phase, which it does in discs with $\alpha\gtrsim 10^{-4}$ at about the same time that the planet undergoes its last doubling (Fig. 1). Under these conditions, the planet is in viscous communication with more-or-less the entire gap, which retains enough mass at its periphery to supply a Jupiter’s worth of gas at $\sim$10 AU (by construction in the minimum-mass nebula). Another way to replenish gas is via hydrodynamic instabilities, not only multidimensional ones but also those in one dimension such as the Rayleigh instability (Tanigawa & Ikoma, 2007; Yang & Menou, 2010; Ginzburg & Sari, 2018, and references therein). These can smear away sharp density contrasts and thereby fill in voids created by the accreting planet. In general, the problem of determining final giant planet masses is tied to the problem of how discs transport mass and angular momentum; the possibility that planet masses are transport-limited deserves further consideration (see also Tanigawa & Tanaka, 2016, and references therein).

We emphasize that the above issues are germane to any theory that wishes to explain the final masses of gas giants by gap opening. As Fig. 2 indicates, the gap depletion factors necessary to stop accretion at Jupiter-like masses are rather large, on the order of $\rho_{0}/\rho\sim 10^{5}$, a result that follows solely from the Bondi accretion rate and the assumption of a background minimum-mass disc, and not from the specific gap-opening theory we have used. Although we have shown in this paper that such large depletions are achievable within the theory of GS18, they motivate further work on hydrodynamic instabilities and radial transport.

Future investigations can also try to lift our restriction to sub-thermal masses, which as Fig. 3 indicates has prevented us from probing small orbital distances. Super-thermal masses pose a variety of challenges, some of which are discussed in GS18. The perturbations that super-thermal planets induce in the disc are too strong to be described by the linearized equations for wave excitation (Goldreich & Tremaine, 1980) and propagation (Goodman & Rafikov, 2001; Ginzburg & Sari, 2018). Furthermore, as $M$ grows above $M_{\rm thermal}$, and the flow near the planet becomes increasingly rotational and less spherically symmetric, the gas accretion rate should deviate from the classical Bondi formula that we have used. Rotation seems particularly interesting to consider. In general, whether or not the planet is super-thermal, the angular spin velocity of gas accreted by the planet at its outer edge is of order the orbital Keplerian shear $\Omega$. For sub-thermal planets, $\Omega$ is less than $\omega_{\rm breakup}$, the planet’s break-up angular velocity. Under these circumstances, the planet’s envelope is not rotationally supported; it can contract and allow more gas from the nebula to be accreted. But for super-thermal planets whose atmospheres extend to $R_{\rm H}$, $\Omega\sim\omega_{\rm breakup}$; such atmospheres cannot contract without losing angular momentum. If they cannot do that—and see, e.g., Szulágyi et al. (2014) and Batygin (2018) for recent thinking on this issue—then the thermal mass naturally emerges as a limit on the planet’s final mass (based on considerations that are, to leading order, orthogonal to the ideas of gap opening explored in this paper). Hydrodynamical simulations along these lines have been conducted by Fung et al. (2019, in preparation).

Acknowledgements

We thank Jeffrey Fung, Eve Lee, Steve Lubow, Ruth Murray-Clay, Mickey Rosenthal, Re’em Sari, and Takayuki Tanigawa for discussions. We also thank the anonymous reviewer for comments which improved the paper. SG is supported by the Heising-Simons Foundation through a 51 Pegasi b Fellowship.

Appendix A Time-dependent gaps

GS18 calculated the equilibrium gap depth in low-viscosity discs. In this Appendix we generalize their calculation to derive the time-dependent depth, before equilibrium is reached. We adopt their dimensionless notation in which $G=M_{\star}=a=1$.

GS18 found that the gap is composed of two density “steps”. The first step is generated by waves that originate at a distance $h$ from the planet (this step corresponds to the classical calculation of Fung et al., 2014) and deposit their angular momentum at a distance $w_{1}>h$. Waves that are generated at $w_{1}$ can further raise the density profile by depositing angular momentum at $w_{2}>w_{1}$ (see Fig. 2 of GS18, who also show that the third and later steps are insignificant). We denote the surface densities at $h$, $w_{1}$, and $w_{2}$ by $\Sigma$, $\Sigma_{1}$, and $\Sigma_{2}$, respectively. We change the notation slightly with respect to GS18; here we use $\Sigma$ instead of their $\Sigma_{0}$, and $m$ instead of $\mu$ to denote the mass ratio.

We begin by analysing the dynamics of the first step. The torque that the planet generates by interaction with gas at $h$ is given by $\Sigma m^{2}/h^{3}$ (Goldreich & Tremaine, 1980). As long as this torque remains unbalanced by viscosity, gas is displaced from $w_{1}$, which itself expands outward. The angular momentum required to displace an annulus of width $w_{1}$ over its own width is given by $\Sigma_{1}w_{1}^{2}$ (we approximate $h<w_{1}<w_{2}\ll 1$). From these considerations we write

[TABLE]

which relates the wave generation location ($h$) and deposition location ($w_{1}$) according to equation (11) of GS18, generalizing the classical Goodman & Rafikov (2001) result to waves travelling across a deep gap. Equation (20b) also shows that gaps grow wider as they get deeper (as $\Sigma_{1}/\Sigma$ increases). We solve equations (20) for the two unknowns:

[TABLE]

where $t_{1}\equiv(h^{2}/\alpha^{7})^{1/5}$ marks the transition to a viscously balanced regime ($\Sigma m^{2}/h^{3}=\Sigma_{1}\nu=\Sigma_{1}\alpha h^{2}$). The saturated depth in equation (21a) corresponds to the Duffell & MacFadyen (2013) and Fung et al. (2014) scaling. The depth of the gap grows sub-linearly with time because of the increasing width $w_{1}\propto t^{1/7}$. As explained in Section 3.2.1, we disregard the earliest times before a gap has formed (when $\Sigma_{1}/\Sigma=1$).

We write equations analogous to (20) for the second step, which is generated by the planet’s interaction with gas at $w_{1}$ (we integrate the differential torque there; see, e.g., Lubow & Ida, 2010):

[TABLE]

with equation (22b) derived from equation (11) of GS18. We solve equations (22) to obtain the depth of the second step:

[TABLE]

where we substituted $w_{1}$ from equation (21b). At $t>t_{2}\equiv(h^{4}/\alpha^{9})^{1/5}$ the second step reaches viscous equilibrium.

Finally, we combine equations (21a) and (23) to derive the overall depth of the gap

[TABLE]

As indicated in Section 3.2.1, equation (24) applies only when $\alpha<m^{5}/h^{14}$ (GS18). Since $\Sigma=\rho H$ is the density at the bottom of the gap, and $\Sigma_{2}=\Sigma_{0}=\rho_{0}H$ is the unperturbed density at its periphery, we have derived $\rho_{0}/\rho=\Sigma_{2}/\Sigma$ as it appears in equation (17).