Envelopes of embedded super-Earths I. Two-dimensional simulations

TL;DR

This paper uses two-dimensional hydrodynamical simulations to study the structure and dynamics of gaseous envelopes around super-Earth cores embedded in protoplanetary disks, revealing rotational support and shock phenomena.

Contribution

It provides new insights into envelope structure and flow dynamics around super-Earths, highlighting the importance of rotation, core size, and shocks, which are not captured by one-dimensional models.

Findings

Envelopes gain rotational support with increasing core mass.

High-mass cores develop non-axisymmetric shocks in their envelopes.

Gas self-gravity weakly affects envelope structure but influences disk density waves.

Abstract

Measurements of exoplanetary masses and radii have revealed a population of massive super-Earths --- planets sufficiently large that, according to one dimensional models, they should have turned into gas giants. To better understand the origin of these objects, we carry out hydrodynamical simulations of planetary cores embedded in a nascent protoplanetary disk. In this first paper of a series, to gain intuition as well as to develop useful diagnostics, we focus on two-dimensional simulations of the flow around protoplanetary cores. We use the pluto code to study isothermal and adiabatic envelopes around cores of sub- to super-thermal masses, fully resolving the envelope properties down to the core surface. Owing to the conservation of vortensity, envelopes acquire a substantial degree of rotational support when the core mass increases beyond the thermal mass, suggesting a limited applicability of one-dimensional models for describing the envelope structure. The finite size of the core (relatively large for super-Earths) also controls the amount of rotational support in the entire envelope. Steady non-axisymmetric shocks develop in the supersonic envelopes of high-mass cores, triggering mass accretion and turbulent mixing in their interiors. We also examine the influence of the gas self-gravity on the envelope structure. Although it only weakly alters the properties of the envelopes, the gas gravity has significant effect on the properties of the density waves triggered by the core in the protoplanetary disk.

Tables

Table 1. Table 1: List of isothermal self-gravitating simulations with their label, pressure scale H 𝐻 H , Bondi radius B 𝐵 B , Toomre parameter Q 𝑄 Q , and Bondi mass m B subscript 𝑚 B m_{\mathrm{B}} relative to the core mass.

Label	$H$	$B$	${\log }_{10}Q$	$m_B/m_c$
IH16B16Q100	$16$	$16$	$2$	$6.7\times {10}^{-4}$
IH16B16Q10	$16$	$16$	$1$	$7.2\times {10}^{-3}$
IH16B16Q3	$16$	$16$	$1/2$	$3.9\times {10}^{-2}$
IH32B4Q3	$32$	$4$	$1/2$	$3.2\times {10}^{-3}$
IH32B8Q100	$32$	$8$	$2$	$1.9\times {10}^{-4}$
IH32B8Q10	$32$	$8$	$1$	$2.1\times {10}^{-3}$
IH32B8Q3	$32$	$8$	$1/2$	$7.0\times {10}^{-3}$

Equations

$$H$$

$$B$$

$$\frac{B}{H}=\frac{m_c}{m_{\mathrm{th}}}$$

$$m_{\mathrm{th}}=\frac{c_s^3}{\Omega G}=m_{\star }{\left(\frac{h}{a}\right)}^3\approx 40m_{\oplus }\left(\frac{m_{\star }}{m_{\odot }}\right){\left(\frac{h/a}{5\%}\right)}^3$$

$$\Phi ={\Phi }_c+{\Phi }_g+q{\Omega }^2{x}^2,$$

$${\partial }_t\rho$$

$${\partial }_t\left[\rho v\right]$$

$${\partial }_{t}E$$

$$\left[\rho ,v_r,P,{\Phi }_c\right](r_c-x)=\left[+\rho ,-v_r,+P,+{\Phi }_c\right](r_c+x).$$

$${\Phi }_c(r)=\frac{G{m}_c}{2ϵ}\log \left(\frac{\sqrt{r^2+{ϵ}^2}-ϵ}{\sqrt{r^2+{ϵ}^2}+ϵ}\right),$$

$$-\underset{\text{gravity}}{\underbrace{{\partial }_r\Phi }}-\underset{\text{pressure}}{\underbrace{{\partial }_{r}P/\rho }}+\underset{\text{inertial support}}{\underbrace{\left(2\Omega +v_{\phi }/r\right)v_{\phi }}}=0.$$

$$r_{\mathrm{H}}\equiv {\left(\frac{m_c}{3{\Omega }^2}\right)}^{1/3},$$

$${\rho }_{\mathrm{hs}}(r)={\rho }_0\exp \left[-\frac{\Phi (r)}{c_s^2}\right],$$

$${\partial }_t\varpi +v\cdot \nabla \varpi =\varpi \cdot \nabla v.$$

$${\partial }_t⟨\bar{\rho }⟩=-\frac{1}{r}{\partial }_r\left[r⟨\bar{\rho }\bar{v}+\overline{{\rho }^{\prime }{v}_r^{\prime }}⟩\right],$$

$${\partial }_t⟨r\rho v_{\phi }⟩=-⟨\rho {\partial }_{\phi }\Phi +2\Omega r\rho v_r⟩-\frac{1}{r}{\partial }_r⟨r\left(r\rho v_{\phi }\right)v_r⟩.$$

$$\overline{\rho v_r{v}_{\phi }}=\bar{\rho }{\bar{v}}_r{\bar{v}}_{\phi }+\bar{\rho }\overline{v_r^{\prime }{v}_{\phi }^{\prime }}+\overline{{\rho }^{\prime }{v}_r^{\prime }}{\bar{v}}_{\phi }+\overline{{\rho }^{\prime }{v}_{\phi }^{\prime }}{\bar{v}}_r.$$

$$({\partial }_t+v\cdot \nabla )n=0.$$

$$\mathcal{B}\equiv \frac{v^2}{2}+\mathcal{H}+\Phi -{\Phi }_{\mathrm{H}},$$

$$\Theta \equiv T{\left(\frac{P_0}{P}\right)}^{\left(\gamma -1\right)/\gamma }\sim \exp \left(s\right).$$

$$Q$$

$$\approx 2.83\times 10^2\left(\frac{h/a}{5\%}\right)\left(\frac{m_{\star }}{m_{\odot }}\right){\left(\frac{\Sigma }{500\mathrm{g}{\mathrm{cm}}^{-2}}\right)}^{-1}{\left(\frac{a}{1\mathrm{AU}}\right)}^{-2}$$

$$\Gamma$$

$$={\oint }_{\mathcal{C}(r)}v\cdot \mathrm{d}\ell +2\Omega \pi (r^2-r_c^2)\simeq 2\pi r{⟨v_{\phi }+\Omega r⟩}_{\mathcal{C}}.$$

$$\frac{\partial \log \rho }{\partial \log r}$$

$$\frac{1}{r}{\partial }_r\left[r{v}_{\phi }\right]$$

$$\Delta {\Phi }_g=4\pi \rho$$

Full text

††pubyear: 2019††pagerange: Envelopes of embedded super-Earths I. Two-dimensional simulations–B.3

1 Introduction

The Kepler survey has revealed that about $50\%$ of Sun-like stars host at least one planet with radius between $1.25-4\mathrm{R}_{\oplus}$ and orbital period less than $145$ days — the so-called super-Earths (Batalha et al., 2013; Fressin et al., 2013). Measurements of both the mass and radius of these exoplanets show that their density decreases as their size increases beyond $1.5\mathrm{R}_{\oplus}$ (Weiss & Marcy, 2014). This density decrease is associated with the presence of gaseous atmospheres (Rogers, 2015), with volatile masses reaching several percent of the mass of the solid core (Lopez & Fortney, 2014; Wolfgang & Lopez, 2015).

In the standard paradigm of planet formation (Safronov, 1969; Kusaka et al., 1970), solid cores of protoplanets grow by accreting planetesimals (Nakagawa et al., 1983; Wetherill & Stewart, 1989) and, possibly, pebble-sized grains (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012). This scenario is supported by the observed distributions of exoplanet masses and radii (Matsuo et al., 2007; Howard et al., 2012; Mordasini et al., 2012). While embedded in a circumstellar disk, the growing core is surrounded by a dense and extended envelope (Perri & Cameron, 1974; Rafikov, 2006). The observed atmospheres are most likely remnants of such primitive envelopes, modified during later stages of photoevaporation (Owen & Jackson, 2012; Owen & Wu, 2013), impact erosion (Inamdar & Schlichting, 2016; Yalinewich & Schlichting, 2018) and outgassing (Chachan & Stevenson, 2018).

As the core grows by accretion of solids, the mass of its hydrostatic envelope increases faster than the mass of the core. Quasi-static envelopes contract and accrete gas as fast as radiative cooling permits (Lee & Chiang, 2015). Hydrostatic and thermal equilibrium cannot be sustained beyond a critical core mass, triggering a runaway cooling and accreting phase (Mizuno et al., 1978; Pollack et al., 1996). At orbital separations $0.1-1\mathrm{AU}$, the “runaway gas accretion” phase is triggered above a critical core mass $5-20\mathrm{M}_{\oplus}$ (Rafikov, 2006), even for vigorously accreting cores. Higher-mass cores are expected to turn into gas giants over the lifetime of the disk (Rafikov, 2011), according to one-dimensional models. Owing to their masses, many super-Earths should thus have turned into gas giants (Lee et al., 2014). That this did not happen is a puzzle.

Various mechanisms have been proposed to moderate or prevent the runaway growth of gaseous envelopes. Examples include the ambient pressure drop following disk clearing (Owen & Wu, 2013; Ginzburg et al., 2016), the tidal (or internal) heating of close-in atmospheres (Ginzburg & Sari, 2017; Ginzburg et al., 2018) or the enhanced opacity of dust-enriched disks (Lee et al., 2014). Hydrodynamic mechanisms have also been studied, primarily under simple thermodynamic assumptions. Due to conservation of the potential vorticity (alias vortensity), two-dimensional envelopes should become rotationally supported (Miki, 1982). If rotational support becomes dominant, a circumplanetary disk should form and regulate gas accretion (Rivier et al., 2012; Tanigawa et al., 2012). The formation of such disks depends on the efficiency of radiative cooling against heating processes, either viscous (Ayliffe & Bate, 2009) or shock-induced (Szulágyi et al., 2016).

Ormel et al. (2015b) pointed out the possible importance of atmospheric recycling — dynamic exchange of mass between the envelope and the surrounding disk — for preventing the transition of the super-Earths into a gas giant. A notable feature of their model is the inclusion of a spatially resolved inner boundary in the simulations domain, mimicking the presence of a solid core. Most previous studies have represented the core via a smoothed gravitational potential (e.g., D’Angelo & Bodenheimer, 2013) and/or a mass sink (e.g., Machida et al., 2010), focusing on the properties of the flow on the scale of the disk thickness. However, Ormel et al. (2015a) showed the influence of this boundary on the global properties of two-dimensional planetary envelopes.

This paper is the first in a series of numerical studies of super-Earth atmospheres, with increasing degree of realism and complexity. Its goal is to go beyond the work of Ormel et al. (2015a) and to provide detailed 2D exploration of the atmospheres of massive (super-thermal, see section 2) cores, which would be capable of triggering runaway gas accretion according to 1D models (Mizuno, 1980). We investigate the role of envelope recycling, with a specific care taken of the core boundary. In doing this we pay particular attention to the balance of pressure versus rotational support in the envelope, and to the mixing properties of the flow.

The plan of the paper is the following. We present our model and its numerical implementation in section 2. Two appendices complete this section with a series of tests. Section 3 contains the results of isothermal, adiabatic, and finally self-gravitating simulations. Section 7 provides a discussion of our results in light of previous studies, and summarize our main conclusions in section 8.

2 Method

2.1 Physical setup

We consider a solid core orbiting a star on a circular trajectory in the midplane of a protoplanetary disk. In this paper, we focus on the two-dimensional properties of the flow around the core, with simple assumptions about the thermodynamics of the gas. The gas self-gravity is included in some cases, see section 6.

Let $r_{c}$ be the radius of the core, $m_{c}$ its mass, $a$ its semi-major axis, and $\Omega=(Gm_{\star}/a^{3})^{1/2}$ its Keplerian frequency around the central star of mass $m_{\star}$. The sound speed of the gas $c_{s}$ is linked to the hydrostatic pressure scale height of the disk $h=c_{s}/\Omega$. Another important lengthscale of the problem is the Bondi radius111The Bondi radius is often defined with the Keplerian escape velocity $r_{\mathrm{B}}=2Gm_{c}/c_{s}^{2}$, twice larger than our definition. $r_{\mathrm{B}}=Gm_{c}/c_{s}^{2}$.

We define the dimensionless ratios $H\equiv h/r_{c}$ and $B\equiv r_{\mathrm{B}}/r_{c}$, and evaluate their typical values for a super-Earth core orbiting a Solar-mass star:

[TABLE]

For an Earth-like planet at $1\mathrm{AU}$, both $H$ and $B$ would be much larger than these numerical estimates. However, for characteristic super-Earth-like core masses and separations from the star, $B$ can be of order unity and $r_{c}$ can reach several percent of $h$. As a result, the physical core size might be able to influence the flow structure on scales $\sim h$. This is an important characteristic of super-Earths that distinguishes them from e.g. the cores of gas giants at 5-10 AU.

From (2) and (2) one can form the dimensionless ratio

[TABLE]

which is independent of the core radius. This ratio is equivalent to the ratio of the core mass $m_{c}$ to its ‘thermal mass’ (Rafikov, 2006):

[TABLE]

At this important mass scale the perturbations induced by the core in the surrounding disc (e.g., density waves) become non-linear. Also, $m_{c}=m_{\rm th}$ corresponds to $r_{\mathrm{B}}=h$, and both scales are approximately equal to the core’s Hill radius defined as $r_{\mathrm{H}}=a\left(m_{c}/3m_{\star}\right)^{1/3}$ (up to a factor $3^{1/3}$).

Note that $h$ and $r_{\mathrm{B}}$ are defined with respect to the temperature of the disk. If the disc is non-isothermal, then any heating process occuring near the core will increase the effective pressure scale and decrease the effective Bondi radius. In this work we use either isothermal ($P=\rho c_{s}^{2}$) or adiabatic ($P\propto\rho^{\gamma}$) equation of state for the gas. The isothermal limit would correspond to instantaneous radiative heating/cooling, while the adiabatic limit assumes that no radiative losses occur. In the adiabatic case, we prescribe the exponent $\gamma=7/5$, as would be appropriate for a diatomic gas in three dimensions (even though we use a 2D setup in this work). For simplicity, chemical effects such as thermal dissociation and reactions between different species are not taken into account. Also, unlike Kley (1999), we consider only inviscid flows in which momentum dissipation should be limited to shocks (Lubow et al., 1999).

2.2 Numerical setup

We aim at resolving scales ranging from a fraction of the core radius to several pressure scales $h$. Since $H\gg 1$, the envelope should isolate the core from the global dynamics of the disk. We follow the core along its orbit around the star in ‘the local approximation’ (Hill, 1878). Let $(x,y)$ be Cartesian coordinates with the origin at the center of the core, with $x$ axis pointing from the star to the core and $y$ along the orbit of the core. The total gravitational potential $\Phi$ is split into three components: the core potential $\Phi_{c}$, the gas self-gravity $\Phi_{g}$, and the tidal potential of the star in the co-orbiting frame. In the local approximation, this decomposition reads

[TABLE]

where $q\equiv\Omega^{-1}\partial_{x}v_{y,0}$ is the dimensionless shear rate of the background flow $v_{y,0}$ (unperturbed by the core). We set the shear rate to its Keplerian value $q=-3/2$ throughout this paper.

For simplicity, we neglect the head-wind experienced by the core in sub-Keplerian disk (its effect was previously examined by Ormel et al. (2015a)). If the midplane pressure $P\propto r^{-n}$, then the core faces a head-wind with Mach number $v/c_{s}\simeq(n/2)(h/a)$, i.e. a few per cent for typical values of the opening angle $h/a\approx 5\%$.

A key aspect of our model is the inclusion of the core as a physical boundary. Matter should not be able to cross the boundary, allowing the accumulation of mass and angular momentum on top of the core. However, resolving the core comes at a cost. At spatial scales smaller than the core size, hydrodynamic fluctuations evolve on time scales shorter than $r_{c}/c_{s}=(\Omega H)^{-1}\ll\Omega^{-1}$. Following the dynamics on scales $\sim r_{c}$ over several orbits is computationally affordable only for moderate values of $H$. Luckily, the time required to restore a quasi-static equilibrium is the sound crossing-time $\sim\Omega^{-1}$ (see Miki, 1982, and Appendix A). For this reason, we focus on time intervals of a few tens of orbits only, short relative to disk clearing or planet migration timescales (Gorti et al., 2016; Fung & Lee, 2018). In this sense, we look at the quasi-instantaneous state of the envelope in this work.

We call $\rho$ the gas (surface) density, $v$ the velocity, $e$ the internal energy per unit mass, $P$ the thermal pressure, $\Phi$ the gravitational potential from (5) and $E\equiv\rho\left(e+v^{2}/2+\Phi\right)$ the total energy density. With these notations, the equations describing the dynamics of mass, momentum and energy are:

[TABLE]

We close this system with either an isothermal ($P=\rho c_{s}^{2}$) or adiabatic ($P=\left(\gamma-1\right)\rho E$) equation of state. In the isothermal case, only (6) and (7) are actually integrated.

2.2.1 Integration scheme

We use the finite-volume code pluto (Mignone et al., 2007) to integrate (6)$-$(8) in conservative form. The primitive variables $\left(\rho,v,P\right)$ are estimated at cell interfaces by linear reconstruction with Van Leer’s slope limiter (Van Leer, 1979). Godunov fluxes are then computed via the Roe approximate Riemann solver (Roe, 1981). A shock-flattening strategy is set to stabilize the solver in regions of strong pressure contrast. If the relative pressure variations between neighboring cells exceed a factor of $5$, we locally switch to the MINMOD slope limiter and to the HLL approximate Riemann solver (Van Leer, ). The time-stepping is performed via an explicit second-order Runge-Kutta scheme; a Courant number of $0.3$ is used for CFL stability.

Specific care is taken of the Coriolis acceleration: to ensure the local conservation of angular momentum, the equations are discretized in a frame rotating with angular velocity $\Omega$ around the core axis (Kley, 1998; Mignone et al., 2012). When taking the gas self-gravity into account, the gas potential $\Phi_{g}$ is obtained by solving Poisson’s equation as described in Appendix B.

2.2.2 Computational domain

We use polar coordinates $(r,\varphi)$ centered on the core, with $\varphi=0$ along the axis from the star to the core. The computational domain is $\left(r,\varphi\right)\in\left[r_{c},128r_{c}\right]\times\left[0,2\pi\right]$. The radial interval is meshed with $512$ logarithmically spaced grid cells, and the azimuthal interval is uniformly meshed by $640$ cells. This corresponds to $73$ cells over $\left[r,2r\right]$ for any $r$. The convergence of most diagnostics was tested for resolutions ranging from $64$ to $512$ radial cells (see Appendix A).

The logarithmic grid spacing allows a fine spatial resolution near the core at the expense of computational time. The CFL constraint is dominated by sound waves in the innermost grid cells. For moderate values of $H=16$, approximately $5\times 10^{4}$ time steps are required to integrate one orbit of the core; this number increases linearly with $H$ in isothermal simulations. The timestep constraint is even more demanding in adiabatic simulations because the sound speed increases near the core. We therefore focus our analysis on $H\leq 32$ simulations; according to (2), this parameter range is most relevant to large cores at small orbital separations in thin disks.

2.2.3 Initial and boundary conditions

The initial conditions consist of the unperturbed shear flow of the disk $\left(\rho,v_{x},v_{y},P\right)=\left(\rho_{0},0,-3\Omega x/2,\rho_{0}c_{s}^{2}\right)$. To avoid a violent relaxation towards the core, the potential of the core $\Phi_{c}$ is gradually introduced, increasing linearly over the first two orbits (see Appendix A).

We impose periodic boundary conditions in the azimuthal direction. In our local setup we cannot capture certain features of the global dynamics of the disk (e.g., the opening of a gap by the core) without ad-hoc prescriptions at the outer radial boundary. The shear flow is supersonic with respect to the boundary for $\left|x\right|>2h/3$, so unless filtering out the characteristic waves that propagate inward, any density/velocity perturbation reaching the boundary should generate discontinuities inside the computational domain. For simplicity, we impose the initial conditions of an unperturbed shear flow at the outer radial boundary. Because of the sustained perturbations induced by the planet on the flow, a velocity discontinuity appears at the outermost radial grid cells in the portion of the flow leaving the computational domain. We expect no influence of this discontinuity in the inner regions close to the core. As a precaution, the outer regions $r>64r_{c}$ are always be excluded from our analysis.

At the inner radial boundary, we require that mass should not flow through it. This is achieved by the adequate symmetrization

[TABLE]

The even symmetry of $\Phi_{c}$ across the boundary effectively cancels the gravitational acceleration at the interface. The residual mass flux through the boundary is always monitored in our simulations.

There is not a unique acceptable boundary condition for the azimuthal velocity $v_{\varphi}$ at the core surface. This choice will likely affect the long term behavior of our simulations by injecting or extracting angular momentum from the envelope. Let $\omega\equiv\nabla\times v$ be the vorticity, and $\varpi_{z}\equiv\left(\omega_{z}+2\Omega\right)/\rho$ be the vertical vortensity of the flow. In the inviscid case this component should remain constant in the absence of shocks (see Ormel et al., 2015a, and section 7.2). We use this property to design a stress-free inner boundary. Since $v_{r}$ is already set by the zero-mass-flux condition, we can impose $\omega_{z}$ (and therefore $\varpi_{z}$) by setting the derivative $\partial_{r}v_{\varphi}$. The velocity $v_{\varphi}(r)$ in the ghost cells is estimated by first-order integration from the active domain, followed by $8$ implicit Jacobi iterations. The desired vortensity $\varpi_{z}$ is thus maintained to better than $1\%$ accuracy. Since vortensity is not conserved around high-mass cores (see section 3.3), we impose the azimuthally-averaged value of the vortensity on top of the core instead of the background $\varpi_{0}$.

We emphasize that the system (6)-(8), along with the prescribed boundary conditions, admits an infinity of steady two-dimensional solutions. In the absence of mixing, whether viscous or turbulent, isothermal solutions can be parametrized by a radial profile of vortensity near the core, and adiabatic solutions admit the entropy profile as an additional degree of freedom. By slowly introducing the potential of the core, we are studying the class of solutions having a constant vortensity and entropy away from shocks.

2.2.4 Gravitational potential of the core

In this simplified study, we neglect the vertical dimension of the protoplanetary disk: the flow is constrained to evolve in the midplane of the disk, where the core lies. If this model is meant to represent a three-dimensional disk in some average sense, then the gravitational potential of the core $\Phi_{c}$ should be averaged correspondingly.

It is customary to smooth the gravitational potential of the core in order to avoid singularities when the core size is unresolved at the grid scale (Müller et al., 2012). Let $\Phi_{\mathrm{Newton}}=-Gm_{c}/r$ be the Newtonian potential of a point mass $m_{c}$; the most common smoothing technique (Plummer, 1911) includes a smoothing length $\epsilon$ via Plummer potential $\Phi_{\mathrm{Plummer}}(r)=-Gm_{c}/\sqrt{r^{2}+\epsilon^{2}}$. To avoid spurious mass fluxes through the core boundary of their three-dimensional simulations, Ormel et al. (2015b) used a potential $\Phi_{\mathrm{Ormel}}$ that is force-free ($\partial_{r}\Phi_{\mathrm{Ormel}}=0$) near the core surface (see their equation 2). As indicated above, we prevented such mass leaks by full symmetrization of the ghost cells at the inner boundary.

In this work we opted to integrate the Newtonian potential $-1/\sqrt{r^{2}+z^{2}}$ over the height of the core $z\in\left[\pm r_{c}\right]$ to obtain the softened potential

[TABLE]

with the smoothing length $\epsilon=2r_{c}$. Figure 1 illustrates these different smoothing techniques. In all cases, the smoothing becomes apparent only near the core radius. Our potential $\Phi_{c}$ nearly coincides with $\Phi_{\mathrm{Plummer}}(\epsilon=r_{c})$ for $r>r_{c}$ so we expect no significant difference with a Plummer type of smoothing.

2.2.5 Units and conventions

From now on, we set the gravitational constant $G=1$. We take the orbital frequency $\Omega$ and the core radius $r_{c}$ as frequency and distance units. With this choice, the isothermal sound speed of the disk $c_{s}=\Omega h=H$ and the core mass $m_{c}=r_{\mathrm{B}}c_{s}^{2}=B\,H^{2}$. Without self-gravity, we take the background disk density $\rho_{0}$ as density unit; this value is tuned in the self-gravitating case, see section 6.

We label each simulation run by its equation of state (I for isothermal, A for adiabatic), its pressure scale height H# and its Bondi radius B#. Spatial averages are represented by brackets $\left\langle\cdot\right\rangle$, taken in the azimuthal dimension $\varphi$ by default.

3 Results for the fiducial isothermal simulation

We present our results by first examining in detail the outcomes of a single (fiducial) run. We take the run IH16B16 as such a reference isothermal simulation. With $B=H$, the core mass equals its thermal mass (Rafikov, 2006), delimiting the transition from low-mass to high-mass cores. This simulation is integrated for a total of twenty orbital periods of the core.

Figure 2 shows the steady flow obtained in IH16B16 after time-averaging over $5$ orbits (ten snapshots). As described by Miki (1982), the flow can be partitionned into several distinct regions: the background shear flow (leftmost and rightmost parts), the co-orbiting gas on horse-shoe orbits (lower and upper parts), and closed streamlines around the core. The fluid on closed streamlines has a prograde orientation $v_{\varphi}>0$. The pressure scale height $h$ demarcates the dense envelope from the background shear flow. Spiral density waves are launched near the edge of the envelope and saturate in stationary shock waves. These shock waves extend out to the outer boundary of the computational domain. They maintain a near-parabolic shape due to the background shear; the amplitude of the density jump increases at first as a result of angular momentum injection by the planetary torque (Goodman & Rafikov, 2001) before saturating at $\Delta\rho/\rho_{0}\approx 2$ away from the core due to numerical dissipation.

3.1 Radial momentum balance

In a steady and non-accreting state, the azimuthally-averaged velocity $\left\langle v_{r}\right\rangle\simeq 0$ (as supported by our numerical results), so that we can express radial momentum balance as

[TABLE]

These contributions are drawn on Figure 3 and compared to the Keplerian acceleration $a_{\mathrm{K}}\equiv m_{c}/r^{2}$. The gravitational acceleration is sub-Keplerian at small radii due to the smoothing of the core potential $\Phi_{c}$. The gravitational acceleration changes sign at the Hill radius

[TABLE]

separating the core-dominated potential from the tidally-dominated potential. The shock waves form near this transition, causing the steep increase in pressure support. The total radial acceleration is negligible inside $r<r_{\mathrm{H}}$, so the envelope is indeed in steady state. The gravitational acceleration is balanced by inertial (Coriolis + centrifugal) and pressure support in roughly equal proportions, so the envelope deviates significantly from a hydrostatic solution.

3.2 Rotational support

The departure from a hydrostatic equilibrium can be measured by comparing the density to its hydrostatic value. For isothermal flows with azimuthal symmetry, the hydrostatic density profile is

[TABLE]

depending only on the local potential $\Phi$ and sound speed $c_{s}$. In the opposite limit of no presure support, a flat density profile can be sustained by fully rotational support if the azimuthal velocity equals222We neglect the contribution of the Coriolis acceleration to the total inertial support; this term is always small near the core of rotationally supported envelopes, see section 7.2. $v_{\mathrm{K}}\equiv\sqrt{r\partial_{r}\Phi}$.

The density and angular velocity profiles of run IH16B16 are drawn on Figure 4, respectively normalized by $\rho_{\mathrm{hs}}$ and $v_{\mathrm{K}}$. The density always increases near the core surface, but not as fast as predicted by (13); the density at the core surface reaches only $4\times 10^{-3}$ of its hydrostatic value in this run. The importance of rotational support is signified by $v_{\varphi}/v_{\mathrm{K}}$ reaching up to $70\%$ in the envelope. The azimuthally-averaged gravitational acceleration is outward outside approximately one Hill radius (see Figure 3) because of Keplerian shear in the disk, so there can be no rotational support beyond.

3.3 Vortensity conservation and generation

Rotationally supported envelopes form in our simulations as a result of vortensity conservation. In inviscid barotropic flows, the vortensity $\varpi\equiv\left(\nabla\times v+2\Omega\right)/\rho$ evolves according to

[TABLE]

In two dimensions, the axial component $\varpi_{z}$ can vary only at shocks (e.g., Dong et al., 2011); it is otherwise conserved along streamlines. This quantity determines the amount of rotational support achieved in the envelope (see Ormel et al., 2015a, and section 7.2).

In the regime $B\gtrsim H$, the stationary shock waves visible in Figure 2 allow vortensity generation. Figure 5 reveals substantial deviations from the Keplerian vortensity value $\varpi_{0}$ in run IH16B16. Vortensity is generated at the shocks and advected downstream. The shock fronts extend upstream of the core into the co-orbital region (lower left and upper right). This material gains a significant amount of vorticity despite the shocks being weak there ($\Delta\rho/\rho\rightarrow 0$, see Figure 2).

The flow obtained in this simulation is not strictly steady, and Figure 5 reveals filamentary patterns along the shock front in the post-shock medium. These filaments are only resolved by a few grid cells, but they do not vanish upon time averaging. They appear in simulations with $B\geq H$ and always develop along the shock front, with no clear signature in the density distribution. With a local shear rate comparable to the rotation rate of the flow, it is tempting to interpret these filamentary patterns as marginally stable sheared waves; a detailed stability analysis is however beyond the focus of this paper. No such features appear in the low-mass regime $B/H<1$, for which the flow is stationary in agreement with Ormel et al. (2015a).

3.4 Mass and momentum transport

In our simulations the flow features small-scale fluctuations over a quasi-steady state. This scale separation allows us to represent a flow variable $X$ as a laminar component $\overline{X}$ (moving average over time) plus short-timescale fluctuations $X^{\prime}$. Under this decomposition, the radial mass transport obeys

[TABLE]

where the second term in the right-hand side comes from correlations of density and velocity fluctuations. These two contributions to mass flux are respectively labeled as ‘laminar’ and ‘turbulent’. We estimate the time-averages $\overline{X}$ by stacking $21$ snapshots spanning $10$ orbits in the quasi-steady state of run IH16B16. We compute the fluctuations by subtracting $X^{\prime}=X-\overline{X}$ in each snapshot; we then compute the correlations between fluctuating terms and average them over the same set of snapshots.

The laminar and turbulent mass fluxes are drawn on the upper panel of Figure 6. The turbulent mass flux is negligible inside the envelope. The laminar mass flux is oriented toward the core and causes mass accumulation near the core surface. The mass contained inside the Bondi radius increases by only $0.2\%$ over ten orbits, so the flow is quasi-steady to a high degree. This increment corresponds to the mass flux passing through the Bondi radius, there is no measurable mass loss into the core through the inner radial boundary.

The same analysis can be performed for the radial transport of angular momentum, as seen in the non-inertial frame:

[TABLE]

The first bracket includes gravitational tides and the Coriolis acceleration. The second bracket corresponds to the radial flux of momentum advected by the flow. We decompose this flux into one laminar and three turbulent components:

[TABLE]

The different terms of (16) are shown on the lower panel of Figure 6. Inside the envelope, the derivative $\partial_{t}\left\langle r\rho v_{\varphi}\right\rangle$ fluctuates around zero with amplitudes small relative to $\left\langle\Omega r\rho v_{\varphi}\right\rangle$. The angular momentum of the envelope is thus conserved to a good approximation. The turbulent momentum flux is small at all radii. At the envelope boundary $r\approx h$, the gravitational torque induced by the star excites density waves that carry angular momentum away from the core.

Since the net torque inside the envelope is negligible, mass accretion is driven by momentum losses at the envelope boundary. Numerical dissipation induces negligible accretion rates at the current resolution (cf. Appendix A), so the isothermal shocks are responsible for momentum dissipation. This accretion mechanism differs from viscously-driven accretion (e.g., Kley, 1999) and primarily concerns massive cores in 2D (see Lubow et al., 1999, and section 4.3).

3.5 Envelope recycling

One of the goals of this study is to characterize the recycling of the envelope material for different core masses. As recycling represents the exchange of mass between the envelope near the core and the adjacent parts of the disc, unbound to the core, we need to come up with ways of characterizing this exchange. Here we describe two such methods, utilizing a passive tracer section 3.5.1 and gravitational binding of the envelope section 3.5.2, and illustrate their application in our fiducial run.

3.5.1 Tracer fluid

To examine the mixing properties of the flow, we follow the transport of a passive tracer in run IH16B16. We let the flow evolve toward a quasi-steady state over eight orbits after the core potential has been fully introduced; we then inject a tracer with concentration $n=1$ inside the Bondi disk. The tracer concentration is then advected by the flow, obeying

[TABLE]

We set the tracer concentration to zero in the radial ghost cells, so the boundaries can only absorb the traced fluid.

Figure 7 shows the evolution of the tracer concentration over $10$ orbits since injection in run IH16B16. At small radii $r/r_{\mathrm{B}}\lesssim 1/2$, the tracer concentration remains unity, meaning that this region does not exchange matter with the outer envelope. This region is defined by having only closed streamlines in the time-averaged (‘laminar’) flow, confirming that turbulent velocity fluctuations do not contribute to gas mixing in the inner envelope. The tracer concentration drops by more than three orders of magnitude near $r/r_{\mathrm{B}}\approx 1/2$, beyond which it is advected away from the core in less than two orbits. Recycling thus appears to be either switched on or off, depending on the topology of the time-averaged flow alone. In principle, the fluid orbiting on closed streamlines around the core could radiate its internal energy, and eventually remain bound to the core by gravity. To determine the extent of such a gravitationally bound envelope, we move on to an energetic arguments below.

3.5.2 Gravitationally bound envelopes

If the nebula in which the core is embedded were to disperse with time, the envelope would expand to match the reduced ambient pressure (Owen & Jackson, 2012). For the core to keep a gravitationally bound atmosphere, the energetic content of the gas must not allow it to escape the potential well of the core. The critical energy barrier $\Phi_{\mathrm{H}}$ is the total potential (5) evaluated at one Hill radius. Whether the envelope is gravitationally bound can be assessed via the Bernoulli number

[TABLE]

where the enthalpy is $\mathcal{H}\equiv c_{s}^{2}\log\rho$ in isothermal flows. We verified that $\mathcal{B}$ is accurately conserved along streamlines in the time-averaged flow of our simulations. A gravitationally bound streamline should therefore satisfy $\mathcal{B}<0$.

Figure 8 represents the radial distribution of Bernoulli number and its different contributions in run IH16B16. The kinetic term is smaller than the enthalpy contribution despite the flow orbiting the core at sonic velocity $v_{\varphi}/c_{s}\approx 1$. The Bernoulli number becomes negative inside $r/r_{\mathrm{B}}\lesssim 1/5$, corresponding to approximately $3r_{c}$ in this case. Since the flow is immediately maintained to its initial temperature, neither adiabatic compression nor shocks can heat the envelope. Isothermal envelopes represent the coldest solutions given a background disk temperature, so they provide an upper bound on the size of the gravitationally bound region. As long as the kinetic contribution is small, the Bernoulli number is only a function of $r_{\mathrm{B}}/r$ inside the Hill radius; the radius at which $\mathcal{B}=0$ should therefore be a fraction of the Bondi radius. We consistently find this radius near $r_{\mathrm{B}}/5$ in every isothermal simulation, regardless of $B/H$. If the disk were to disperse, we would expect the outer parts of the envelope $r/r_{\mathrm{B}}\gtrsim 1/5$ to expand and escape the gravity of the core. This region includes closed streamlines, unaffected by recycling under the current pressure confinement from the disk onto the envelope.

4 Parameter exploration: isothermal simulations

We now examine how the outcomes found in a fiducial case change as we vary $B$ and $H$ in a series of isothermal simulations. This corresponds to changing the mass of the core $m_{c}=B\,H^{2}$ (in units of $m_{\rm th}$) and its radius relative to the disk scale height $r_{c}/h=H^{-1}$.

4.1 Pressure versus rotational support

Figure 9 displays the radial profiles of density relative to the hydrostatic value $\rho_{\rm hs}$ in a series of isothermal simulations. For a fixed $H$ (upper or lower panel), increasing $B$ (the planet mass) leads to progressively higher densities at the same radius. However, the density does not increase as fast as predicted for a hydrostatic profile (13); the ratio $\rho/\rho_{\mathrm{hs}}$ decreases with $B$, as Figure 9 shows. Comparing the two panels, we see that curves with the same $B$ but different $H$ are essentially identical but offset by a factor $2$ horizontally. The amount of pressure support $\rho/\rho_{\mathrm{hs}}$ thus depends on $B$ and on the distance to the core surface, but only weakly on $H$. In other words, the ratio of $B/H$, independent of $r_{c}$, does not uniquely determine the amount of pressure support as a function of $r/h$.

As $\rho/\rho_{\mathrm{hs}}$ decreases, the amount of rotational support $v_{\varphi}/v_{\mathrm{K}}$ increases toward unity. Therefore, the amount of rotational support in the entire envelope also depends on $r_{c}/h$. We show in section 7.2 that this dependence on the core radius is expected to vanish in the limit $H\rightarrow\infty$.

4.2 Tracer fluid

After the potential of the core is introduced, we let the flow settle to a quasi-steady state over eight orbits in simulations with different values of $H$ and $B$. We then inject a passive tracer fluid with concentration $n=1$ inside the Bondi disk, and let it evolve over ten orbits. The final tracer distributions are drawn on Figure 10.

As in section 3.5.1, every simulation features a concentration plateau below a critical radius, so tracer mixing is inefficient in the direct vicinity of the core. Comparing both panels one notes that, for a fixed value of $H$, this critical radius remains near $r_{\mathrm{H}}/3$ as long as $B\leq H$, i.e. $m_{c}\leq m_{\rm th}$. For larger values of $B>H$ the critical radius moves much closer to the core, meaning that most of the envelope was recycled by the disk flow over the duration of the simulation. Closed streamlines are still present in the time-averaged flow of every simulation, but a fraction of these streamlines have a reduced tracer concentration. We therefore attribute recycling to time-dependent fluctuations. Figure 10 suggests the onset of turbulent tracer transport beyond a critical core mass $B/H>1$.

4.3 Mass and tracer transport around massive cores

We recall that in run IH16B16 the envelope displays an extended tracer plateau (Figure 10) and essentially no turbulent mass flux (Figure 6). However, with a core twice more massive, run IH16B32 features tracer mixing much deeper in the envelope. We examine the radial mass fluxes in run IH16B32 on Figure 11.

The net mass flux is oriented towards the core, with an amplitude $\sim 10^{3}$ larger than in the fiducial case IH16B16. The mass contained within the Bondi disk increases by $40\%$ over ten orbits, corresponding to the mass flux passing through the Bondi radius. The mass flux through the inner radial boundary fluctuates around zero with amplitudes at the level of $\pm 10^{-6}m_{c}$ per orbit, negligible compared to mass accretion rate of the envelope. Compared to run IH16B16, there is now a significant turbulent mass flux inside the envelope, reaching amplitudes larger than $10\%$ of the laminar one. The turbulent mass flux of (15) allows mass to flow across the time-averaged velocity field. It allows the mixing of the envelope with the background flow, regardless of momentum transport.

Figure 12 shows the time-averaged vortensity distribution in run IH16B32. Compared with IH16B16 on Figure 5, we see that stationary shocks emerge inside the Hill radius. These shocks alter the vortensity distribution in the rotating envelope, the sign of the vortensity jump being sensitive to the geometry of the shock (Kevlahan, 1997). The core is surrounded by a low-vortensity disk extending to $\lesssim 4r_{c}$, and delimited by a high-vortensity ring near $5r_{c}$. The low-vortensity region encloses the remaining tracer fluid visible in Figure 10. The shock locations are also correlated with the presence of a turbulent mass flux in Figure 11. We therefore attribute the envelope recycling to the formation of steady shocks in the circulating flow within the Hill sphere of the core.

Momentum dissipation at these shocks also drives the laminar mass flux toward the core. As vortensity decreases near the core, the role of rotational support decreases (see Ormel et al., 2015a, and section 7.2). To maintain radial momentum balance, the envelope must increase its pressure support by accreting mass. Without including the gas gravity, this accretion process cannot produce envelopes more massive than the hydrostatic limit given by Equation 13. On long time scales, the asymptotic distribution of vortensity (rotational support) will be sensitive to the viscous processes at play (e.g., turbulence) and cannot be uniquely determined a priori in an inviscid context.

5 Parameter exploration: non-isothermal simulations

Having explored the characteristics of the flow in isothermal setup, we now change our thermodynamic assumptions and simulate the flow with a adiabatic equation of state $P\propto\rho^{\gamma}$. The energy equation (8) is integrated conservatively, and the entropy per unit mass $s$ is conserved along streamlines as long as the flow does not shock. This situation can be relevant in the optically thick limit, when the heat due to adiabatic compression is not efficiently radiated away. The exponent is $\gamma=7/5$ everywhere and at all time, neglecting chemical and ionization effects. Fresh gas with the entropy of the background disk is continually supplied at the outer radial boundary. We still expect convergence to quasi-steady states, as no energy losses are allowed and entropy production at the shocks is rather slow. The control parameters $H$ and $B$ (which depend on temperature) now represent the thermodynamic properties of the background disk flow.

5.1 Adiabatic versus irreversible heating

As the core mass is progressively introduced in our adiabatic runs, the envelope contracts and heats adiabatically, the temperature of the gas $T\equiv P/\rho$ increasing as $\rho^{\gamma-1}$. Once the gravitational potential of the core is set, the gas still heats and cools adiabatically as it passes by the core. In addition to adiabatic heating, the shocks forming at the envelope boundary induce irreversible heating, changing entropy $s$ of the gas (Rafikov, 2016). To measure the amount of irreversible heating at shocks, we compute the potential temperature

[TABLE]

This quantity represents the temperature that a fluid element would have if its pressure were adiabatically relaxed to the background disk pressure $P_{0}$. Being only a function of the local entropy, $\Theta(s)$ is also conserved along streamlines and increases at shocks. We normalize it to the background disk temperature $T_{0}=H^{2}$.

Figure 13 shows the distribution of potential temperature in the quasi-steady state of run AH16B32. Compared to the isothermal equivalent (see Figure 12), the incoming gas shocks further away from the core in the adiabatic case, so the envelope spans a larger area. The potential temperature increases by up to $10\%$ at the envelope boundary; it is then passively advected inside and out of the envelope. The absence of irreversible heating near the core implies that the shocked gas never blends into the inner envelope.

The upper panel of Figure 14 shows the radial profiles of temperature in a series of adiabatic simulations. The temperature always decreases with radius and increases with the mass of the core. The temperature at the surface of the core is up to seven times larger than its background value in run AH32B32. The temperature profiles of runs AH16B16 and AH32B32 are nearly superimposed as a function of $r/h$. Unlike isothermal envelopes, the radial structure of adiabatic envelopes apparently depends only on $B/H$ and not on the core radius. This behavior is expected if rotational support is small and if the entropy remains close to its background value (e.g., equation 52 of Rafikov, 2006).

The lower panel of Figure 14 shows the radial profiles of potential temperature in the same series of adiabatic simulations. We observe irreversible heating for $B\gtrsim H$, but the potential temperature variations remain below $10\%$ for the range of parameters considered. The shock-induced heating is small compared to the total heating (see the upper panel of Figure 14), so these envelopes can be considered as isentropic to a good approximation. In runs AH16B8 and AH16B16, the spiral shocks do not extend upstream of the core into the co-orbital flow, as is the case on Figure 13. The shocked material is advected away without mixing into the envelope, so the heat deposition is localized near $r\approx h$. In runs AH16B32 and AH32B32, the shocked gas blends into the envelope and mildly increases its entropy down to $r\gtrsim h/8$. Inside $r\lesssim h/8$, the potential temperature keeps its initial value, implying that the shocked gas never enters the innermost regions of these envelopes.

5.2 Envelope recycling

To examine the efficiency of envelope recycling, we inject a tracer fluid in the quasi-steady state of each adiabatic simulation. The tracer concentration $n=1$ inside the Bondi disk and zero outside, and it is passively advected over ten orbits. The final concentration profiles are drawn on Figure 15. Recycling affects the envelope inside the Hill radius in every case. Every run except AH16B8 features a concentration plateau near the core. The radial extent of this plateau shows no obvious scaling with $B$ or $H$. Both runs AH16B16 and AH32B32 feature a concentration plateau extending to $\approx r_{\mathrm{H}}/4$. Their concentration profiles are also similar to the isothermal equivalent IH16B16 (see Figure 10). Contrarily to their isothermal analogues, there is no concentration plateau near the core in run AH16B8, and there is an extended plateau in run AH16B32.

In run AH16B8, the envelope is rotationally supported at only $v_{\varphi}/v_{\mathrm{K}}\leq 3\%$, instead of $30\%$ in the equivalent isothermal run IH16B8. Only the innermost streamlines circulate around the core, so most of the envelope is recycled by the shear flow on orbital timescales. In run AH16B32, the envelope is rotationally supported at $v_{\varphi}/v_{\mathrm{K}}\approx 30\%$, instead of $70\%$ in run IH16B32. Because of the temperature increase, the circulating flow remains subsonic in run AH16B32, with $v_{\varphi}/c_{s}\lesssim 60\%$ instead of $v_{\varphi}/c_{s}\in\left[1,3\right]$ in run IH16B32. We showed that stationary shocks develop in the circulating flow of run IH16B32 (see Figure 12), and we correlated these shocks to turbulent mixing (see Figure 11). Because the envelope of run AH16B32 remains subsonic, no shocks develop and turbulent tracer transport does not operate. Laminar and circular streamlines extend up to $r\lesssim r_{\mathrm{H}}/3$, corresponding to the extent of the tracer plateau visible on Figure 15.

6 Parameter exploration: self-gravitating disks

As the gas density near the core increases, the mass of the envelope could become comparable to the mass of the core. The gravitational acceleration induced by the gas could then alter the structure of the envelope even when background disk is not susceptible to a gravitational instability. We test the sensitivity of our previous results to inclusion of the gas self-gravity in a series of isothermal simulations.

Gravitational instability is expected to occur in massive disks, when the Toomre parameter

[TABLE]

is below unity. In an isothermal, hydrostatic and non-self-gravitating disk, the midplane density $\rho_{0}$ is linked to the surface density $\Sigma$ via $\rho_{0}=\Sigma/h\sqrt{2\pi}$. We sample values of $Q\in\left[10^{1/2},10^{2}\right]$ by tuning the background density $\rho_{0}$. Whether realistic or not in the context of super-Earths formation (Chiang & Laughlin, 2013), such low values of $Q$ are required to reveal gas-gravity effects in the envelope.

Given the gas density distribution, we solve Poisson’s equation to obtain its gravitational potential. The specific method and its validation are presented in Appendix B. Characteristics of self-gravitating simulations are listed in Table 1, where the labels are appended by the corresponding value of $Q$.

6.1 Self-gravitating envelopes

In the limit $Q\rightarrow\infty$, the gravity of the gas is negligible compared to the gravity of the core, so the density distribution should scale linearly with the background density $\rho_{0}$. If we decrease $Q$ (increase $\rho_{0}$) while maintaining the distribution of $\rho/\rho_{0}$, then the mass of the envelope increases relative to the mass of the core and adds up to it. Self-gravitating envelopes are therefore expected to pull more gas toward the core and become more massive.

We integrate the mass $m_{\mathrm{B}}$ inside the Bondi disk and find that it remains below $4\%$ of the mass of the core at the end of every two-dimensional simulation presented here (see Table 1). To measure the contribution of the gas gravity, we compare $m_{\mathrm{B}}$ to the mass $m_{0}\equiv\pi(r_{\mathrm{B}}^{2}-r_{c}^{2})\rho_{0}$ of the background disk over the same area. The ratio $m_{\mathrm{B}}/m_{0}$ increases by less than $30\%$ when decreasing $Q$ in our sample of simulations. Run IH32B8Q3 has a density corresponding to $10^{-1/2}\approx 32\%$ of the critical density for gravitational fragmentation, which is $10^{3/2}\approx 32$ times larger than the density of run IH32B8Q100, yet the ratio of $m_{\mathrm{B}}/m_{0}$ is only $16\%$ larger. As in their non-self-gravitating analogues, runs IH32B32Q100 and IH32B32Q10 feature mass accretion onto the core due to steady shocks inside their rotating envelope. The mass accretion rate normalized by $\Omega m_{0}$ is $25\%$ larger in run IH32B32Q10 compared to the less massive case IH32B32Q100.

We previously demonstrated the importance of rotational support for the envelope dynamics. We now examine how the radial momentum balance of the envelope is affected by the gas self-gravity. We draw on Figure 16 the radial profiles of $v_{\varphi}/v_{\mathrm{K}}$ measured in the time-averaged flow of a series of self-gravitating simulations. Note that $v_{\mathrm{K}}\equiv\sqrt{r\partial_{r}\Phi}$ includes the mass of the enclosed gas. We find that the degree of rotational support $v_{\varphi}/v_{\mathrm{K}}$ decreases with $Q$, i.e. it increases with the mass of the envelope. The gas gravity affects the envelopes of low and high-mass cores similarly with respect to their rotational support. The ratio $v_{\varphi}/v_{\mathrm{K}}$ increases by less than $10\%$ in every case considered. This is to be expected given the $<4\%$ envelope mass increase measured inside the Bondi radius. We conclude that for the range of core masses considered here, two-dimensional envelopes are only moderately affected by their self-gravity.

6.2 Self-gravitating density waves

We find that in our self-gravitating simulations the mass in the spiral density waves can become comparable to the mass enclosed inside two-dimensional envelopes. We illustrate this by showing the density distribution of run IH32B4Q3 on Figure 17. With $B=4$, the density at the surface of the core is only $26$ times larger than the background disk density, corresponding to $29\%$ of its hydrostatic value from (13). As previously, stationary shock waves are launched into the disk around one pressure scale $h$ away from the core. These shocks extend in the co-orbiting (horse-shoe) region beyond $2h$ upstream of the core. Such extended shocks appeared only around high-mass cores in non-self-gravitating simulations; see Figure 12 for example. The core of run IH32B4Q3 is embedded in a homogeneously shocked medium, with circulating streamlines restricted to the innermost $r\lesssim r_{\mathrm{H}}/2$. The rest of the shocked gas is recycled on orbital timescales.

Given the density distribution of Figure 17, the gravitational potential $\Phi_{g}$ solution of Poisson’s equation is represented on Figure 18. The iso-potential contours are elongated along the axis of the shocked and dense region. The potential of the gas remains deeper than $50\%$ of $\Phi_{c}(r_{c})$ over an area several times larger than the Hill disk. The potential of the gas effectively dominates over the potential of the core outside the Hill radius, competing with the tidal potential of the star. Run IH32B4Q3 thus illustrates how a low-mass core can launch self-gravitating density waves into a massive protoplanetary disk. The density perturbations, although not as pronounced as around massive cores, are amplified by their self-gravity beyond one pressure scale away from the core. With stronger shocks spanning a larger area, self-gravity could significantly alter the observational signatures of even the relatively low-mass planetary cores (Dong et al., 2015).

7 Discussion

7.1 Comparison with previous studies

Despite the difference in computational power, our simulations reproduce most of the flow properties already identified by Miki (1982) in two-dimensional inviscid simulations. The flow pattern can be divided into three distinct regions: the background shear, the co-orbital (horse-shoe) flow, and an inner region where streamlines circulate around the core. As a result of vortensity conservation, the innermost regions become rotationally supported (with a prograde orientation) when the envelope mass increases. The radial momentum balance of envelopes forming around massive cores may therefore be far from hydrostatic, contrary to the usual assumption of one-dimensional models (e.g., Mizuno et al., 1978; Pollack et al., 1996; Rafikov, 2006).

Contrarily to Kley (1999) and Nelson et al. (2000), we do not include any explicit viscosity and we do not remove mass at the location of the core. Our setup allows steady envelopes to form with no mass accretion toward the core. However, we do observe mass accretion in a subset of simulations with massive cores, satisfying $r_{\mathrm{H}}/h\gtrsim 1$ ($m_{c}\gtrsim m_{\rm th}$). Shocks form in this regime (Korycansky & Papaloizou, 1996; Lubow et al., 1999) on scales previously under-resolved. Shocks break vortensity conservation. The sign of the vortensity jump depends on the geometry of the shock (Kevlahan, 1997); we find that the inner envelope generally loses vortensity. With a reduced vortensity, the envelope is less rotationally supported. To maintain radial momentum balance, the envelope must increase its pressure support by accreting mass, so shocks effectively lower the ‘centrifugal barrier’ around sufficiently massive cores in two dimensions (Ormel et al., 2015a). In the low-mass regime $B/H<1$, our results agree with those of Ormel et al. (2015a) in terms of vortensity conservation and resulting rotation profiles (see section 7.2 below).

In addition to momentum dissipation, shocks drive temporal variability in the envelope. The ‘turbulent’ mass flux is always small compared to the laminar one, but it allows the mixing of material from the different regions of the flow.

Recycling the envelope material with fresh gas from the disk is a possible way to counteract the radiative cooling of the envelope (Ormel et al., 2015b). We postpone the discussion of recycling to later three-dimensional results, noting that the strength of the shocks should be significantly reduced in three-dimensions (e.g., Bate et al., 2003).

The main caveat of our 2D model is the local approximation. Planets in the parameter range $B/H\gtrsim 1$ are expected to open a gap in the disk (Kley & Nelson, 2012), which we cannot capture without ad-hoc boundary conditions. Global models are also necessary to capture the accumulation of vortensity in the co-orbiting flow, which could affect the migration rate of the planet (Ward, 1991; Koller et al., 2003; Paardekooper & Papaloizou, 2009) and that of dust grains through the orbit of the planet (Weidenschilling, 1977).

7.2 Global sensitivity to the core radius

Results of section 4 demonstrate the sensitivity of the flow structure to the core size $r_{c}$. To gain insight into the origin of this dependence, we analyze a simplified model for the envelope structure to examine its dependence on the size of the core. As mentioned in section 3.3, in barotropic flows the vortensity $\varpi_{z}\equiv\left(\omega_{z}+2\Omega\right)/\rho$ is constant along streamlines. Assuming that the $\varpi_{z}$ has kept its initial value $\varpi_{0}$ everywhere, we integrate the vorticity flux over the disk of radius $r$ around the core and apply Stokes theorem:

[TABLE]

We have neglected the contribution from the inner radial boundary at $r=r_{c}$ and assumed a purely circular flow in (22b). This relation links the velocity on a contour to the mass enclosed inside this contour. It is accurately satisfied in our simulations, from the inner radius up to the envelope boundary at $r\simeq h$. In the hydrostatic limit, the mass contained in the isothermal envelope $\pi r^{2}\left\langle\rho\right\rangle_{\mathcal{D}}$ diverges when the inner radius $r_{c}\rightarrow 0$. The amount of rotational support in such an envelope should therefore depend on the size of the core.

The azimuthal velocity $v_{\varphi}$ must increase with the mass of the envelope, eventually becoming positive (prograde, Miki, 1982). Assuming axisymmetry, we can compute the structure of an isothermal envelope in radial momentum balance:

[TABLE]

where $\varpi_{z}$ is the $z$-component of the vortensity, assumed to be fixed. Ormel et al. (2015a) examined the same system omitting the Coriolis acceleration in (23). The Coriolis term is always negligible around rotationally supported envelopes, but we keep it for the sake of completeness.

We solve this system numerically by imposing $v_{\varphi}(r_{c})=0$, $\rho(h)=1$, and the Keplerian vortensity $\varpi_{z}=\varpi_{0}\equiv\Omega/2\rho_{0}$ everywhere. The solutions are in good agreement with those obtained from direct simulations (see Figure 4). Our purpose here is not to reproduce exactly the two-dimensional results, but to capture the main features of the equilibrium, which do not depend on the exact choice of $\rho(h)$. Given the density and velocity profiles, we can compare the pressure to the inertial support against gravity at every radius in (11). We focus on the profiles of $v_{\varphi}/v_{\mathrm{K}}$ in the prograde ($v_{\varphi}>0$) portion of the envelope; the angular velocity $v_{\varphi}$ smoothly takes negative values outside of this inner region.

First, we set the core size via $H\equiv h/r_{c}=64$ according to (2) and we compute solutions for different core masses by varying $r_{\mathrm{B}}/h$. This family of solutions is represented on Figure 19. The amount of rotational support, as measured by $v_{\varphi}/v_{\mathrm{K}}$, increases in the entire envelope as a function of $r_{\mathrm{B}}/h$. In this case, the threshold $v_{\varphi}/v_{\mathrm{K}}=50\%$ is attained for $B\equiv r_{\mathrm{B}}/r_{c}\gtrsim 8$. The radial extent of the prograde envelope also increases with $r_{\mathrm{B}}/h$ and approaches the outer boundary $r=h$ as $r_{\mathrm{B}}/h$ increases. Solutions with $r_{\mathrm{B}}/h>1$ are not represented because they become excessively stiff and require higher resolutions.

Next, we set the ratio of $B/H=2.4\times 10^{-2}$ as estimated in (3), and we vary the core radius relative to the local pressure scale via $H$. With this parametrization, the ratio $r_{\mathrm{H}}/h$ remains constant, the mass of the core increases as $m_{c}\propto H^{3}$ so the Keplerian velocity $v_{\mathrm{K}}\propto H^{3/2}$. This family of solutions is represented on Figure 20. The envelope becomes more rotationally supported as $H$ increases, i.e. as the core radius $r_{c}$ decreases with respect to $h$. All the simulations presented in this paper fall in this regime sensitive to the core radius. However, the curves appear to converge to an asymptotic profile of $v_{\varphi}/v_{\mathrm{K}}$ over the whole envelope. We conclude that for sufficiently small cores, the amount of rotational support in the envelope should eventually become independent of the core radius.

8 Summary

We performed two-dimensional inviscid hydrodynamic simulations of embedded protoplanetary cores in a local model of Keplerian disk. We focused on the properties of the dense envelope surrounding the core in a regime where its Bondi radius spans across the hydrostatic pressure scale of the disk (i.e. as the core mass is varied from sub-thermal to super-thermal). The core was included as a spatially resolved boundary, and gas was treated as either isothermal or adiabatic. We also implemented a Poisson solver in the pluto code in order to include the gas self-gravity in a subset of simulations. Our main conclusions are the following.

Vortensity is conserved in two-dimensional flows around low (sub-thermal) mass cores, so the amount of rotational support in the envelope depends on the mass inside the envelope and, a priori, on the size of the core. This dependence is expected to vanish in the limit of cores small relative to the pressure scale of the disk, which is computationally more challenging to simulate. 2. 2.

Stationary shocks form inside the envelope of massive cores, when the circulating flow becomes supersonic. By altering the vortensity distribution, these shocks allow mass accretion and effectively break the centrifugal barrier. They also drive a turbulent mixing of the envelope material with the background disk. Stationary shocks also form in the co-orbital flow, far upstream of massive cores. 3. 3.

Adiabatic envelopes are more pressure-supported than their isothermal analogues; they are more efficiently recycled around low-mass cores but less susceptible to shock-induced mixing in the high-mass regime. The shock-induced (irreversible) heating is small compared to the adiabatic heating due to gas compression in every case studied. 4. 4.

The self-gravitating potential of the gas is most sensitive to the large-scale density waves launched into the disk. Self-gravity affects primarily the extent and intensity of these large-scale shock waves, whereas the properties of the inner envelope only weakly depend on the Toomre parameter $Q$ of the disk.

The extension of this study to three-dimensions, with a focus on the recycling properties of the flow, will be the subject of the subsequent work.

Acknowledgements

We thank the anonymous referee for providing constructive comments that improved the clarity of this paper. Financial support of this work by Isaac Newton Trust, Department of Applied Mathematics and Theoretical Physics and STFC through grant ST/P000673/1 is gratefully acknowledged.

Appendix A Convergence study

We performed a series of isothermal 2D simulations with $H=16$ and $B=8$ to evaluate the spatial resolution required for convergence of most diagnostics. The radial domain $r/r_{c}\in\left[1,128\right]$ is meshed by either $64$, $128$, $256$ or $512$ logarithmically spaced cells.

Figure 21 shows the temporal evolution of the surface-averaged density inside the Bondi disk. The average density increases during the first two orbits as the potential of the core is progressively introduced. It takes two more orbits for the envelope to adjust to the final potential and reach a quasi-steady state. The envelope mass varies by less than $10\%$ depending on the spatial resolution. After that, the envelope mass evolves on much longer timescales: the doubling time ranges from $10^{4}$ to over $10^{5}$ orbits from the lowest to the highest resolution. It corresponds to the timescale of kinetic energy dissipation by numerical viscosity in the envelope.

On Figure 22, we examine how the mass accretion flux depends on the grid resolution. This flux is estimated from the time and cell-averaged primitive variables $(\rho,v_{r})$. There is a clear convergence with resolution toward a flat profile with zero mass flux in the envelope: the mass flux in the envelope goes down by almost a factor ten every time the resolution is doubled. We choose to use $512$ radial cells to make such mass fluxes insignificant. The measured mass flux is oriented inward at the core surface, outward in the envelope and inward out of the envelope. This suggests mass depletion near the core surface and accumulation at the envelope boundary. In reality, the inner radial boundary condition cancels the mass flux through the core surface.

We emphasize that the pluto code uses a reconstruction scheme to estimate the primitive variables at the cell interfaces before converting to conservative variables such as the mass flux. Unfortunately, the mass flux actually used by the code is not easily accessible during computations. To confirm the absence of mass losses through the inner radial boundary, we measured the increase rate of the mass contained inside $r<h$. Independently, we averaged the integrated mass flux $2\pi r\rho v_{r}$ inside $r/h\in[1/4,1]$. These two diagnostics match to $10^{-2}$ accuracy when there is a net mass flux through the envelope as in Figure 6 or Figure 11.

Appendix B Self-gravity in pluto

B.1 Implementation

The gravitational potential $\Phi_{g}$ satisfies Poisson’s equation

[TABLE]

in the entire computational domain with appropriate boundary conditions. The structures and routines used for the representation and resolution of the parallel problem of solving this equation come from the PETSc library (Balay et al., 1997, 2018). The Laplacian operator is discretized via second order finite difference. Boundary conditions are applied in the appropriate rows of the operator matrix. Given a density distribution, the linear problem for $\Phi_{c}$ is solved iteratively via a biconjugate gradient method (Yang & Brent, 2002). By construction, the potential $\Phi_{g}$ is defined at cell centers; the interface values used by the hydrodynamic solver are computed by linear interpolation. Note that the boundary conditions are applied in the first active cell of the domain, and not at the boundary nor in the ghost zones.

We implemented and tested this Poisson solver in cartesian, cylindrical and spherical geometries as described in the following sections. In the case of two-dimensional cylindrical geometry, the Green’s function of the Laplacian operator is the potential of an infinite line $\Phi_{g}(r)\sim\log(r)$. Instead, we use the spherical representation of the Laplacian operator over the cylindrical coordinates $(r,\varphi)$, so that the potential of a point mass varies as the three-dimensional Newtonian potential $\sim 1/r$.

The boundary conditions in the angular directions respect the topology of the domain (periodic in the azimuthal angle $\varphi$). The potential is set to $\Phi_{g}=0$ at the outer radial boundary, so as to fix a reference value. By doing so, we enforce an axial symmetry of the gas potential at large distances from the core. We keep this choice for simplicity, since we cannot know a priori the boundary values $\Phi_{g}\left(r_{\mathrm{out}},\varphi\right)$ for non-axisymmetric density distributions. The condition $\partial_{r}\Phi_{g}=0$ is imposed at the inner radial boundary, consistent with the absence of gas below the core radius $r_{c}$. Any other choice would also favor an unwanted mass flux through the inner boundary.

The parallel Poisson solver can significantly affect the performances of the code depending on the spatial resolution and dimensionality of the problem. In 3D spherical geometry, with a resolution of $[128\times 80\times 160]$ on $(r,\theta,\varphi)\in\left[1,128\right]\times\left[0,\pi\right]\times\left[0,2\pi\right]$, the Poisson solver takes the equivalent of $2.3$ hydrodynamic timesteps to determine $\Phi_{g}$ for a constant density distribution and Dirichlet radial boundary conditions. This computational overhead can reach a factor $25$ in high-resolution cylindrical (2D) simulations. Since the large-scale structures of the density distribution do not evolve on short timescales near the core, it is reasonable to solve (25) for $\Phi_{g}$ every $n=4$ hydrodynamic timesteps, which we opted to do in our simulations (see section B.3).

B.2 Static tests

The Poisson solver was tested in static configurations for cartesian, cylindrical and spherical geometries, in two and three dimensions. In these static tests, the mass distribution $\rho$ is not allowed to evolve in time. Simple prescriptions for the source term $\rho$ allow comparing $\Phi_{g}$ with analytical solutions of Poisson’s equation.

In 1D, imposing a constant density $\rho=1$ on the interval $x\in\left[0,1\right]$ with homogeneous Dirichlet conditions, the exact solution is $\Phi(x)=2\pi x(x-1)$. Because the boundary conditions are imposed at the center of the first active cells and not at their boundaries, the approximate $\Phi_{g}$ is offset with respect to the exact one by a constant. After subtracting this constant, the residual error oscillates between $\pm 10^{-7}$, so the first and second derivatives of $\Phi_{g}$ are accurately captured.

We perform a similar test within our standard cylindrical setup by prescribing $\rho(r)=1$ for $r<4r_{c}$ and $\rho(r)=0$ beyond. The numerical and theoretical profiles of $|\Phi_{g}(r)|$ are in excellent agreement as shown on Figure 23. The relative error on $\Phi_{g}$ is less than $1\%$ inside $r\leq 32$ and reaches $2\%$ at $r=92r_{c}$.

To test the proper implementation of the boundary conditions, additional static tests were performed with an asymmetric source term. For fully periodic domains, the compatibility condition333In Fourier space, the zeroth-order (constant) component of (25) must reduce to zero. $\int\rho=0$ over the whole domain is enforced by subtracting the volume-averaged density in (25).

B.3 Dynamic test: Jeans instability

We verified that our implementation performs well in dynamic situations by reproducing Jeans’ instability in Cartesian geometry with periodic boundary conditions. The domain $\left(x,y\right)\in\left[0,1\right]\times\left[0,2\right]$ is meshed with $64\times 128$ cells, the fluid is initialized with zero velocity, and the density $\rho=1+\epsilon$ is flat with a white noise of amplitude $\epsilon=10^{-6}$. The Jeans length is set to $\ell_{\mathrm{J}}=3/2$, so harmonic perturbations are unstable only in the $y$ direction.

We start by computing $\Phi_{g}$ at every timestep of the Runge-Kutta 2 integrator (every two substeps). The linear mode with length $L=2$ has an expected growth rate $s=2\pi\left(c_{s}/L\right)\sqrt{|1-\left(L/\ell_{\mathrm{J}}\right)^{2}|}\approx 5.54\left(c_{s}/L\right)$. As demonstrated on Figure 24, the measured growth rate $5.520\pm 0.001$ matches the theoretical one to $4\times 10^{-3}$ accuracy in this configuration. We reproduced Jeans’ instability by solving (25) every $n=4$ and $n=10$ timesteps, and measured relative errors of $2\%$ and $5\%$ on the growth rates respectively.