Envelopes of embedded super-Earths II. Three-dimensional isothermal simulations

TL;DR

This study uses three-dimensional isothermal hydrodynamic simulations to explore how the properties of gaseous envelopes around super-Earth cores depend on core mass and disc interactions, revealing complex flow patterns and recycling processes.

Contribution

It provides new insights into the flow dynamics and envelope properties of embedded super-Earths across different core masses using detailed 3D simulations.

Findings

High-mass cores experience supersonic polar inflows and shock-driven turbulence.

Envelope properties are sensitive to the ratio of Bondi radius to core size.

Gas recycling influences runaway gas accretion and atmosphere retention.

Abstract

Massive planetary cores embedded in protoplanetary discs are believed to accrete extended atmospheres, providing a pathway to forming gas giants and gas-rich super-Earths. The properties of these atmospheres strongly depend on the nature of the coupling between the atmosphere and the surrounding disc. We examine the formation of gaseous envelopes around massive planetary cores via three-dimensional inviscid and isothermal hydrodynamic simulations. We focus the changes in the envelope properties as the core mass varies from low (sub-thermal) to high (super-thermal) values, a regime relevant to close-in super-Earths. We show that global envelope properties such as the amount of rotational support or turbulent mixing are mostly sensitive to the ratio of the Bondi radius of the core to its physical size. High-mass cores are fed by supersonic inflows arriving along the polar axis and shocking on the densest parts of the envelope, driving turbulence and mass accretion. Gas flows out of the core's Hill sphere in the equatorial plane, describing a global mass circulation through the envelope. The shell of shocked gas atop the core surface delimits regions of slow (inside) and fast (outside) material recycling by gas from the surrounding disc. While recycling hinders the runaway growth towards gas giants, the inner regions of protoplanetary atmospheres, more immune to mixing, may remain bound to the planet.

Tables

Table 1. Table 1: Simulations presented in this paper: label, pressure scale height H 𝐻 H , Bondi radius B 𝐵 B (both normalised by r c subscript 𝑟 𝑐 r_{c} ), maximal percentage of rotational support v φ / v K subscript 𝑣 𝜑 subscript 𝑣 K v_{\varphi}/v_{\mathrm{K}} reached in the equatorial plane, maximal Mach number v z / c s subscript 𝑣 𝑧 subscript 𝑐 𝑠 v_{z}/c_{s} reached by the time-averaged flow in a cone of half-opening angle π / 16 𝜋 16 \pi/16 about the polar axis.

Label	$H$	$B$	$\%$ equatorial $v_{\phi }/v_K$	polar $v_z/c_s$
H16B8	$16$	$8$	$1.20\pm 0.01$	$0.14\pm 0.01$
H16B16	$16$	$16$	$38\pm 1$	$2.4\pm 0.1$
H16B32	$16$	$32$	$86\pm 2$	$5.8\pm 0.1$
H16B64	$16$	$64$	$105\pm 1$	$9.4\pm 0.1$
H32B16	$32$	$16$	$31.3\pm 0.3$	$1.20\pm 0.08$
H32B32	$32$	$32$	$85.4\pm 0.6$	$5.92\pm 0.03$

Equations

$${\partial }_t\rho$$

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

$$\Phi =-\frac{G{m}_c}{r}+q{\Omega }^2{x}^2+\frac{1}{2}{\Omega }^2{z}^2,$$

$$r_{\mathrm{H}}\equiv {\left(\frac{G{m}_c}{3{\Omega }^2}\right)}^{1/3}=a{\left(\frac{m_c}{3m_{\star }}\right)}^{1/3}$$

$$v_x=v_z=0,v_y=q\Omega x,\rho ={\rho }_0\exp \left(-z^2/2h^2\right),$$

$$\frac{r_{\mathrm{B}}}{h}=\frac{m_c}{m_{\mathrm{th}}}=3{\left(\frac{r_{\mathrm{H}}}{h}\right)}^3,$$

$$m_c>m_{\mathrm{B}}\equiv m_{\star }{\left(\frac{h}{a}\right)}^3{\left(\frac{3}{4\pi }\frac{m_{\star }}{{\rho }_c{a}^3}\right)}^{1/2}$$

$$=0.007m_{\oplus }{\left(\frac{m_{\star }}{m_{\odot }}\right)}^{3/2}{\left(\frac{a}{0.1\text{mbox}AU}\right)}^{-1}{\left(\frac{h/a}{5\%}\right)}^3{\left(\frac{{\rho }_c}{5\text{mbox}{gcm}^{-3}}\right)}^{-1/2},$$

$$B\equiv \frac{r_{\mathrm{B}}}{r_c},H\equiv \frac{h}{r_c}.$$

$$\left[\rho ,v_r,v_{\theta },v_{\phi }\right](r,-\vartheta ,\phi )=\left[+\rho ,+v_r,-v_{\theta },-v_{\phi }\right](r,\vartheta ,\phi +\pi ),$$

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

$$\frac{\partial m}{\partial t}=-4\pi r^2\left(\overline{\rho }\cdot \overline{v_r}+\overline{{\rho }^{\prime }{v}_r^{\prime }}\right).$$

$${\partial }_t\left(n\rho \right)+\nabla \cdot \left(n\rho v\right)=0$$

$$\dot{M}\approx 0.1{\rho }_0{r}_{\mathrm{B}}h{c}_s\approx 0.1{\Omega }^{-1}{\rho }_0{a}^3\frac{m_c}{m_{\star }},$$

$$\frac{\partial }{\partial t}{\int }_{\mathcal{V}(r)}n\rho \mathrm{d}V$$

$$\displaystyle=-4\pi r^{2}\bigl{(}\left\langle n\right\rangle\cdot\left\langle\rho v_{r}\right\rangle+\tilde{c}\cdot\|\tilde{n}\|\cdot\left\|\widetilde{\rho v_{r}}\right\|\bigr{)}.$$

$$\tau (r)\equiv 2\frac{{\int }_{\mathcal{V}(r)}\rho \mathrm{d}V}{{\oint }_{\mathcal{S}(r)}\Vert \widetilde{\rho v_r}\Vert \mathrm{d}S}=2\frac{(4/3)\pi r^3{⟨\rho ⟩}_{\mathcal{V}}}{4\pi r^2\Vert \widetilde{\rho v_r}\Vert }\gtrsim \frac{{⟨\rho ⟩}_{\mathcal{V}}}{{⟨\rho ⟩}_{\mathcal{S}}}\frac{r}{c_s},$$

$$\frac{\mathrm{d}\mathbf{\varpi }}{\mathrm{d}t}\equiv {\partial }_t\mathbf{\varpi }+(\mathbf{v}\cdot \nabla )\mathbf{\varpi }=(\mathbf{\varpi }\cdot \nabla )\mathbf{v}.$$

Full text

††pubyear: 2018††pagerange: Envelopes of embedded super-Earths II. Three-dimensional isothermal simulations–Envelopes of embedded super-Earths II. Three-dimensional isothermal simulations

1 Introduction

Over $4000$ planetary candidates have been discovered by the Kepler mission (Borucki et al., 2010; Batalha et al., 2013). Simultaneous measurements of masses and radii available for some planets give insights about their composition, and reveal a population of ‘super-Earths’ — planets with a mass in the range $(1-10)~{}M_{\oplus}$ and a radius in the range $(1-4)~{}R_{\oplus}$. Given the low average density of many such planets (Weiss & Marcy, 2014; Rogers, 2015), up to tens of percent of the mass of super-Earths could reside in their gaseous atmospheres (Lopez & Fortney, 2014; Wolfgang & Lopez, 2015).

In one-dimensional models of embedded planetary envelopes, a planet can accrete as much gas as radiative cooling permits (Lee & Chiang, 2015). Eventually, no static equilibrium can be maintained beyond a critical core mass $\sim 5-20~{}M_{\oplus}$ (Perri & Cameron, 1974; Mizuno et al., 1978; Rafikov, 2006), leading to a phase of runaway gas accretion. This scenario is successful in forming gas giant planets (Mizuno, 1980; Pollack et al., 1996), but it also suggests that, unless forming late in the life of the disc (Lee et al., 2014; Lee & Chiang, 2016), many super-Earths should have turned into gas giants over the lifetime of the disc.

One way to counteract the radiative cooling of planetary envelopes and to prevent runaway core accretion is to continuously replace their material by high-entropy gas from the disc (Ormel et al., 2015b). The efficiency of this ‘recycling’ process motivated dedicated studies of the flow structure around planetary cores (Fung et al., 2015; Kurokawa & Tanigawa, 2018), which confirm the potential importance of this mechanism. In Béthune & Rafikov (2019, hereafter Paper I) we explored the properties of embedded planetary envelopes using two-dimensional (2D) simulations. In this second paper of a series, we extend the 2D results of Paper I to three dimensions by running a suite of 3D simulations of the flow near the core. We explore how the flow properties change as the core mass varies across the so-called thermal mass scale (Rafikov, 2006). We restrict ourselves to isothermal simulations and postpone a more realistic treatment of the thermodynamics to a later study.

As part of describing 3D envelopes, we will primarily focus on differences in flow morphology, mass accretion, the formation of rotationally supported envelopes, and the efficiency of gas recycling as a function of core mass. We will also systematically explore the differences with the 2D results brought about by extending the flow in the vertical dimension. Compared to Bate et al. (2003), we will focus on scales smaller than the pressure scale height of the disc, and will consider gas accretion in an inviscid context.

Another subject of our study is the role of the finite size of the core in determining the flow properties around the core. To handle the large contrast between the core radius and the stratification scale of the disc, many previous studies used a softened gravitational potential and/or sink cells in place of the core (e.g., Kley et al., 2001; Bate et al., 2003; Machida et al., 2010; Tanigawa et al., 2012). Ormel et al. (2015a) and Fung et al. (2014) noticed the importance of properly resolving the core scale, and its influence on the entire envelope, as we verified in 2D in Paper I. We follow the same approach here and fully account for the boundary condition imposed by the core surface.

The paper is organised as follows. We present our physical setup and numerical model in section 2. The results are described starting with a reference case in section 3 and then sampling different values for the core mass and size in section 4. We discuss these results focusing on the differences between the 2D and 3D situations in section 5. Our main findings are summarised in section 6.

2 Method

The method employed to address the current problem is similar to that in Paper I, with the obvious distinction that the flow can now evolve in three dimensions. We will therefore focus on the main distinctions from the 2D case; Paper I should be consulted for more details.

2.1 Physical model

2.1.1 Notations and main equations

We model the 3D flow of the disc fluid in the vicinity of a planetary core with mass $m_{c}$ and radius $r_{c}$ orbiting the central star $m_{\star}$ on a circular orbit with a semi-major axis $a$, so that the corresponding angular frequency is $\Omega=(Gm_{\star}/a^{3})^{1/2}$. Let $\rho$ be the gas density, $P$ its pressure and $v$ its velocity. We adopt an isothermal equation of state $P(\rho)=\rho c_{s}^{2}$ with a constant sound speed $c_{s}$.

We study the dynamics of the disc fluid in a rotating frame centered on the planetary core in the ‘local’ approximation of a Keplerian disc (Hill, 1878). Let $(x,y,z)$ be Cartesian coordinates measured from the center of the core, with $x$ along the radius from the star to the core, $y$ along the orbit of the core and $z$ normal to the disc midplane. Let $(r,\theta,\varphi)$ be spherical coordinates centered on the core, with the polar axis $\theta=0$ along the $z$ direction and $\varphi=0$ along the $x$ direction. We also make use of cylindrical coordinates $(R,\varphi,z)$, with $R=r\sin\theta$ being the cylindrical radius measured from the polar ($z$) axis.

The gas density and momentum evolve according to

[TABLE]

where $\Phi$ is the full gravitational potential detailed below.

2.1.2 Gravitational potential

The gravitational potential of the star is expanded to second order about the orbital radius of the core. It is also expanded to second order in the vertical direction, so the disc will be stratified on a pressure scale $h$. The total gravitational potential $\Phi$ is

[TABLE]

where we assume a Keplerian shear $q=-3/2$. Unlike Paper I, in this work we use the exact Newtonian potential without any softening near the core surface. The transition from the star’s tidal potential to the core-dominated potential occurs at one Hill radius

[TABLE]

away from the core in the equatorial plane. Rotational support of the gas in the envelope around the core will be measured via the Keplerian acceleration $a_{\mathrm{K}}\equiv Gm_{c}/r^{2}$ and the Keplerian velocity $v_{\mathrm{K}}\equiv\sqrt{Gm_{c}/r}$ due to the core only.

In the absence of the core (i.e. with $m_{c}=0$) the potential in the form (3) allows for a steady shear flow in the disc with

[TABLE]

assuming that the midplane density $\rho_{0}$ is constant in space. As in Paper I, we neglect the headwind arising in sub-Keplerian discs and caused by global pressure gradients; its effect has been considered by, e.g., Ormel et al. (2015b); Kurokawa & Tanigawa (2018).

2.1.3 Characteristic space and time scales

The sound speed sets the characteristic pressure scale height of the disc $h\equiv c_{s}/\Omega$ as well as the Bondi radius of the core $r_{\mathrm{B}}\equiv Gm_{c}/c_{s}^{2}$. In the absence of other characteristic length scales, the properties of the flow around a planetary core are expected to be determined by $r_{\mathrm{B}}/h$. This ratio coincides with the ratio of the core mass to its ‘thermal mass’ $m_{\rm th}$ and is closely related to the Hill radius (Rafikov, 2006):

[TABLE]

see definition (4).

The physical size of the core $r_{c}$ is an additional free length scale of the problem, which we will always assume being small compared to the pressure scale, $r_{c}\ll h$. In our simulations the core is represented by a spherical boundary at $r=r_{c}$, with the requirement that matter should not flow through it, see section 2.2.4. The core is capable of inducing substantial perturbation in the disc flow via its gravity whenever $r_{\mathrm{B}}\gtrsim r_{c}$, which sets a lower limit on the core mass of

[TABLE]

where $\rho_{c}$ is the bulk density of the core. In this work we will always be in the regime $m_{c}>m_{\rm B}$, focusing primarily on $m_{c}\sim m_{\rm th}$

As in Paper I, we find it convenient to introduce dimensionless length scales by normalizing $r_{\mathrm{B}}$ and $h$ by $r_{c}$:

[TABLE]

By our assumptions $B>1$ (as $m_{c}>m_{\rm B}$) and $H\gg 1$. Numerical estimates of these parameters typical for Kepler super-Earths can be found in Paper I.

The time required to restore a quasi-static equilibrium through the envelope is the sound crossing-time across one pressure scale $\sim\Omega^{-1}$ (Miki, 1982). However, deep near the core hydrodynamic fluctuations evolve on time scales $\sim r_{c}/c_{s}=(\Omega H)^{-1}$. This scale separation makes it computationally expensive to integrate (1) and (2) over more than a few tens of orbital times ($2\pi/\Omega$) of the core around the central star. Such timescales are short relative to planet migration or disc clearing, so our results should be interpreted as a quasi-instantaneous snapshot in the planet formation history.

2.2 Numerical method

Our numerical approach is to evolve equations (1) and (2) in spherical $(r,\theta,\varphi)$ geometry, naturally extending the calculations done in Paper I to a third dimension. Modifications introduced by this extension are discussed next.

2.2.1 Integration scheme

The numerical method is largely the same as in Paper I. We use the finite-volume code PLUTO 4.0 (Mignone et al., 2007) to integrate (1) and (2) in conservative form. Primitive variables are estimated at cell interfaces by linear reconstruction with VanLeer’s slope limiter (Van Leer, 1979). Godunov fluxes are computed via the Roe approximate Riemann solver (Roe, 1981). The equations are integrated in time via an explicit second order Runge-Kutta scheme; the Courant-Friedrichs-Lewy (CFL) stability criterion is satisfied by timesteps $\Delta t\leq 0.3\Delta x/U$ for characteristic velocities $U$ over a grid spacing $\Delta x$. The Coriolis acceleration is included in a conservative fashion by integrating (2) in a frame rotating with angular frequency $\Omega$ about the polar axis (Kley, 1998; Mignone et al., 2012).

2.2.2 Computational domain

The computational domain is defined in spherical geometry by $\left(r,\theta,\varphi\right)\in\left[r_{c},128r_{c}\right]\times\left[0,\pi\right]\times\left[0,2\pi\right]$. It is meshed by $128\times 80\times 160$ grid cells with a logarithmic spacing in the radial direction and constant spacings in the other dimensions. This moderate resolution will allow us to perform a series of simulation runs at a reasonable computational cost. The dependence of most diagnostics on resolution was examined in Appendix A of Paper I, where convergence was found starting at $128$ cells in the radial direction. This represents approximately $18$ cells between $[r,2r]$ for any $r$.

2.2.3 Stabilizing procedures

The numerical scheme does not intrinsically preserve the positivity of density, so we provide additional dissipation to stabilize the solver in extreme conditions. If the relative pressure variations exceed a factor of $10$ between neighboring cells, or if the local Mach number $v/c_{s}\geq 10$, we switch to the MINMOD slope limiter and to the HLL approximate Riemann solver in the concerned cells (Van Leer, 1997).

For simulations with $H=16$, the computational domain goes as high as eight pressure scale heights above the midplane. In hydrostatic equilibrium, the density given by (5) becomes extremely low at such altitudes. Small momentum fluctuations easily result in high velocities, turbulence and shocks in the upper layers of the disc. Because our interest is on the dense parts of the disc, the gravitational potential of the star $\Omega^{2}z^{2}/2$ is set to a constant value $8\Omega^{2}h^{2}$ for $z>4h$. In the absence of a core, the hydrostatic disc density would therefore be constant above four pressure scale heights.

2.2.4 Initial and boundary conditions

The initial conditions consist of the background, unperturbed shear flow of the disc in the form (5), where we additionally extrapolate $\rho$ to a constant for $z\geq 4h$. To prevent a violent relaxation of the flow toward the core starting from this initial state, the core mass $m_{c}$ is linearly increased from zero to its nominal value over the first $\Omega t/2\pi=2$ orbital times. As soon as the core potential is fully introduced, the properties of the envelope evolve on a much longer timescale (see Figure A.1 of Paper I).

The azimuthal boundary conditions are periodic. The conditions in the polar dimension respect the spherical topology: about the axis $\theta=0\pm\vartheta$ in one hemisphere,

[TABLE]

and similarly in the opposite hemisphere about $\theta=\pi\pm\vartheta$.

At the outer radial boundary, we prescribe the initial, unperturbed state in the ghost cells. Because this condition does not smoothly match the conditions inside the computational domain, a discontinuity appears in the outermost active grid cells. We find no noticeable effect of this discontinuity on the flow near the core; to avoid possible effects of the outer boundary on our results, the outer $R>64r_{c}$ will always be excluded from our analysis. We emphasize that disc-scale processes such as the formation of a density gap cannot be captured in our model without ad-hoc prescriptions at the outer radial boundary.

At the inner radial boundary, we prevent mass inflows by the adequate symmetrization in the ghost cells

[TABLE]

Enforcing an even symmetry for $\Phi_{c}$ helps cancelling the gravitational acceleration at the interface. We monitor the mass losses through this boundary by integrating cell-averaged mass fluxes. The mass losses through the inner radial boundary are negligible to better than $10^{-2}$ accuracy compared to the mass accretion rates through the envelope, so the core surface is indeed impermeable.

The choice of $\left(v_{\theta},v_{\varphi}\right)$ in the inner radial ghost cells will likely influence the long-term evolution of the envelope. In the absence of explicit viscosity, the spin of the core can only affect its envelope via turbulent or numerical diffusion. In Paper I, we used a two-dimensional conservation argument to prescribe $v_{\varphi}$ in a stress-free fashion. There is no simple and equivalent prescription in three dimensions. We opt for simplicity to make boundary effects easier to identify. The meridional velocity is set to $v_{\theta}=0$ in the inner radial ghost cells. The azimuthal velocity $v_{\varphi}$ could be influenced by a number of torques resulting, for example, from the angular momentum of accreted planetesimals or from tidal synchronization with the star. We set $v_{\varphi}=(q/2)\,r\sin\left(\theta\right)$, i.e. a solid-body rotation with the same vorticity as in the Keplerian flow. This corresponds to the steady state of a solid sphere under the friction of a viscous shear flow.

2.3 Units and conventions

The orbital frequency $\Omega$ and the core radius $r_{c}$ are taken as frequency and distance units. The isothermal sound speed of the disc is $c_{s}=\Omega h=H$ and we take the gravitational constant $G=1$, so the core mass becomes $m_{c}=BH^{2}$ in these units. We take the midplane density of the background disc $\rho_{0}$ as density unit. We label each simulation run by its isothermal pressure scale height H# and its Bondi radius B#.

In our analysis we will sometimes separate a given quantity $X$ into a ‘secular’ (moving-average) component $\overline{X}$ and a fluctuating component $X^{\prime}$. The secular term $\overline{X}$ is estimated by averaging over a set of simulation snapshots; it is assumed to vary only moderately on the time interval considered. By default, we use six snapshots spanning three orbits of the core when doing this averaging. The fluctuations are computed by subtracting the average $\overline{X}$ in each snapshot. Spatial averages are denoted with brackets $\left\langle\cdot\right\rangle$. Depending on the context, they will refer to azimuthal averages or to averages on the sphere $\left\langle X\right\rangle_{\mathcal{S}}=\oint_{\mathcal{S}}X\,\mathrm{d}S/4\pi r^{2}$.

The parameters and main diagnostics of the simulations presented in this paper are listed in Table 1. The highest $H=32$ considered in this study is still smaller by about a factor of 2 compared to what is expected for super-Earths at 0.1 AU from the star (Rogers, 2015). This restriction comes from the CFL constraint on sound waves at the grid scale near the surface of the core. The confirmation of our findings at larger scale separations $H>32$ and longer timescales $\Omega t/2\pi>100$ will require additional computational resources.

We refer to the runs using their values of $H$ and $B$. We take H16B16 as our reference (fiducial) case, against which the other simulations will be compared.

3 Results: fiducial simulation

We start by presenting the results of the fiducial simulation run H16B16. As it has $B/H=1$, the core mass equals the thermal mass $m_{c}=m_{\rm th}$ (Rafikov, 2006), corresponding to the transition from low to high-mass cores. In this regime, gravitational torques caused by the core are expected to trigger a non-linear response in the disc close to the planet (Korycansky & Papaloizou, 1996), i.e. density perturbations at $r\sim h$ away from the planet should be of order $\rho_{0}$.

3.1 Flow structure

We start by discussing the time-averaged characteristics of the flow in run H16B16. Figure 1 shows the density and velocity distributions in the equatorial plane after averaging over three orbital periods. The flow can be naturally separated into distinct regions (Fung et al., 2015): the background shear flow (left and right sides), the co-orbital flow completing U-turns (top and bottom), and circular streamlines near the core. The flow in these different regions roughly follows the contours of constant gravitational potential. Density wakes are launched from about one pressure scale away from the core, extending into spiral density waves on larger scales.

In comparison with the 2D case, the gas density ratio between the atmosphere on top of the core and the background disc is several orders of magnitude larger in 3D. This difference is discussed in more detail in section 4.3. The closed streamlines orbit the core with a prograde orientation ($v_{\varphi}>0$), but they span a smaller area than the equivalent 2D flow.

Figure 2 shows the time-averaged flow in the poloidal plane $(R,z)$ after averaging in the azimuthal direction $\varphi$. Deviations from spherical symmetry clearly appear in the density contours. A key feature obvious in this figure is that the mass flow is oriented toward the core at high latitudes and away from the core in the equatorial plane. Every simulation presented in this paper displays such a global circulation, with the same orientation — outflow near the equatorial plane — as found by, e.g, Bate et al. (2003), Fung et al. (2015), Lambrechts & Lega (2017), Kuwahara et al. (2019). The isodensity contours are pinched toward the core near the poles and outward in the equatorial plane. The polar inflow reaches a Mach number $v_{z}/c_{s}\approx 2.4$ on the polar axis before shocking on the dense inner envelope, at $z\approx 2.4r_{c}$. This transition from supersonic to subsonic poloidal flow happens at the lower part of the sonic surface on Figure 2. The maximum Mach number reached in the time-averaged polar inflow of each run is indicated in Table 1.

3.2 Quasi-steady state

At the beginning of a simulation, every diagnostic smoothly evolves as the mass of the core is progressively increased. Once the gravitational potential of the core is fully set, the evolution becomes much slower. Momentum balance is established on sound crossing timescale $h/c_{s}=\Omega^{-1}$, so a quasi-steady state is effectively achieved after four orbits. Even though the mass of the envelope keeps slowly evolving in time, the mass accretion rate $\dot{M}$ is roughly constant over the limited integration time of each simulation (after the initial build-up phase during which the core mass is increased).

On Figure 3, the radial acceleration $a_{R}\equiv\partial_{t}v_{R}$ is compared to the Keplerian acceleration $a_{\mathrm{K}}\equiv m_{c}/R^{2}$. We decompose $a_{R}$ into several components arising in the equation of motion: a gravitational term $-\nabla\Phi$ including the tidal potential of the star, a pressure term $-\nabla P$, and the inertial term including $(v\cdot\nabla)v$ and the Coriolis acceleration. Each term is evaluated from the time-averaged flow variables and then azimuthally averaged in the equatorial plane.

The sum of these different terms is close to zero: the flow maintains a radial momentum balance dominated by pressure against gravity. The pressure jump marks the boundary of the envelope at one pressure scale away from the core. An important feature of the momentum balance in 3D is that the rotational support is weaker than in the equivalent 2D case (cf. Figure 3 of Paper I). As a result, pressure support dominates over rotational support in the equatorial plane of this 3D run.

The upper panel of Figure 4 shows that there is a net mass flux accreting through the envelope. As $\partial_{r}\left[r^{2}\rho v_{r}\right]\approx 0$, this accretion flow is steady in time except in a narrow shell on top of the core where mass accumulates. Sampling the density and velocity distributions over time, we can separate the secular and fluctuating components of the radial mass flux. The mass $m(r,t)$ inside a sphere of radius $r$ evolves on long (secular) timescales as

[TABLE]

The mass inside $r<h$ increases almost linearly in time, with a volume averaged density growing as $\left\langle\rho\right\rangle/\rho_{0}\sim 1.77\,\Omega t/2\pi$. At this rate, it would take approximately $300$ orbits to reach the mass of the corresponding hydrostatic envelope. As our simulations are substantially shorter, we do not reach the stage of the envelope being in nearly hydrostatic equilibrium, although we do detect a reduction of the inflow velocity in time. We recall that without self-gravity or radiative transfer, there is no characteristic density scale, so the envelope mass is proportional to the background disc density $\rho_{0}$ regardless of the core mass. The rate at which the envelope grows matches the surface-integrated mass flux to better than $10^{-2}$ accuracy, so there is no noticeable mass loss through the inner radial boundary (surface of the core).

The mass flux $\left\langle\overline{\rho^{\prime}v_{r}^{\prime}}\right\rangle$ estimated from the density and velocity fluctuations is positive at every radius. This fluctuating mass flux does not only come from turbulent motions, but also from the fact that the ‘secular’ mass flux evolves over the time-averaging interval. Indeed, as the envelope grows more massive ($\left\langle\rho^{\prime}\right\rangle>0$), the accretion velocity tends to decrease ($\left\langle v_{r}^{\prime}\right\rangle>0$), so the correlation $\left\langle\overline{\rho^{\prime}v_{r}^{\prime}}\right\rangle$ has a non-zero secular component that is generally positive. By narrowing the time-averaging interval, we can reduce the amplitude of this secular contribution and still estimate the amplitude of turbulent motions. We found that restricting the averaging interval to three orbits is sufficient to extract a statistically meaningful turbulent mass flux. Doing so, we find that the turbulent mass flux is smaller than the secular component by more than a factor of $10$ at every radius, i.e. negligible with respect to mass accretion. This is similar to what was found in 2D case in the analogous run with $B=H=16$ (see Paper I).

We measure the contribution of rotational support in the envelope of run H16B16 by computing the maximal value of $v_{\varphi}/v_{\mathrm{K}}\approx 38$ per cent in the equatorial plane. To become more massive (more pressure-supported), the envelope must lose momentum. The lower panel of Figure 4 shows the net torque exerted on the fluid inside the envelope. The net torque oscillates around zero; its maximal amplitude corresponds to a momentum removal timescale of $\sim 10^{2}$ orbits. Mass accumulates only near the core surface, so accretion is likely allowed by momentum losses in its boundary layer. In the midplane, the core surface can exchange angular momentum with the surrounding fluid by both turbulent and numerical diffusion. Dedicated boundary layer simulations would be required to disentangle these two effects (e.g., Belyaev et al., 2013). In the polar regions, the supersonic inflows terminate in shocks (see Figure 2), and these shocks accurately delimit the region where mass accumulates. Because the gas is taken to be isothermal, there is no heat generation at shocks, so energy is effectively dissipated there. Because the shocked gas comes from high latitudes, it carries virtually no angular momentum. The isothermal shocks thus allow mass to accumulate deep inside the envelope without any need for angular momentum extraction.

3.3 Time variability

Velocity fluctuations can enhance the mixing of the gas from different regions of the envelope, as well as the exchange of gas between the disc and the envelope. Even without inducing a net mass flux, turbulent diffusion could affect the thermal structure of an envelope by homogenizing its entropy toward the background disc value. We start adressing this issue by examining the amplitude of velocity fluctuations in our isothermal simulations.

The envelope of run H16B16 features small-scale variability on short timescales. We estimate the turbulent Mach number $\mathrm{M}\equiv\tilde{v}/c_{s}$ from the velocity fluctuations $\tilde{v}^{2}\equiv\overline{v^{\prime}\cdot v^{\prime}}$, and represent its equatorial and poloidal distributions on Figure 5 and Figure 6 respectively.

In the equatorial plane, Figure 5 reveals variability in the entire envelope. The turbulent Mach number is maximal inside $r/r_{\mathrm{H}}<1/2$ with $\mathrm{M}\approx 30$ per cent. There is a thin shell directly on top of the core where the turbulent Mach number drops below $10$ per cent. None of our isothermal simulations show supersonic turbulence in the equatorial plane: the turbulent Mach number always saturates below $50$ per cent.

Figure 6 shows that the velocity fluctuations are stronger away from the midplane. The variability at $z/h\approx\pm 0.2$ is excited by helical motions of the fluid as shown on Figure 10 of Fung et al. (2015); these sources of variability are not axisymmetric. Variability is mostly excited in the polar regions, with velocity fluctuations reaching supersonic amplitudes. The turbulent Mach number goes as high as $\mathrm{M}\approx 1.2$ at the polar shocks, half the amplitude of the time-averaged inflow Mach number. The sonic surface marked on Figure 2 actually fluctuates on short time and spatial scales near the core. In the post-shock region surrounding the core, the turbulent Mach drops from $\mathrm{M}\approx 1$ at $z=2~{}r_{c}$ down to $\mathrm{M}\approx 10^{-2}$ at the core surface. Numerical dissipation might damp the turbulent motions at the core surface, but most of the momentum dissipation has already occured through the polar shocks. The only simulation run maintaining subsonic variability in the entire envelope is H16B8, with a turbulent Mach number $\mathrm{M}\lesssim 6$ per cent on the polar axis.

3.4 Envelope recycling

We now look at the issue of the envelope recycling by exchange of gas between the disc and the envelope interior (Ormel et al., 2015b). In Paper I we used passive tracers to diagnose fluid exchanges between the different parts of the flow; this is a method we resort to in this work as well. In our fiducial run the flow can be considered as being in a quasi-steady state after ten orbits of integration time. At this point we inject a tracer fluid with constant concentration $n=1$ inside the Bondi sphere. The tracer is passively advected by the flow according to $\partial_{t}n+v\cdot\nabla n=0$, or equivalently

[TABLE]

for the traced mass. We follow the evolution of the traced fluid over ten additional orbits since tracer injection.

Figure 7 shows the tracer distribution in the equatorial plane of run H16B16 ten orbits after tracer injection. The remaining tracer is concentrated inside $r\lesssim r_{\mathrm{H}}/2$, where the time-averaged flow features circular streamlines orbiting the core (see Figure 1). The tracer concentration displays a spiral pattern offset with respect to the spiral density waves. These spiral streams do not result from gas compression: they reveal how the traced fluid is channelled out of the envelope by gas flowing from high latitudes. The tracer concentration is still of order unity at the surface of the core.

Figure 8 shows the average tracer distribution in the poloidal plane of run H16B16. The iso-concentration contours form lobes restricted to $\left|z/h\right|\lesssim 0.3$ and pinched toward the core on the polar axis. The traced fluid is never transported high above the midplane, it is only depleted there. The tracer concentration has barely decreased from its initial value near the polar caps of the core. In this region where $R\lesssim r_{c}$, recycling is inefficient up to the polar shocks at $\left|z\right|\approx 2.4~{}r_{c}$, inside of which radial motions are largely suppressed.

Figure 9 shows the evolution over time of the shell-averaged tracer concentration in run H16B16. It takes approximately two orbits to reach a quasi-steady concentration profile. Compared to two-dimensional simulations, there is no plateau in the concentration profile: the tracer concentration is substantially reduced by envelope mixing down to the core surface at $r\approx r_{c}$. The average tracer concentration keeps slowly decreasing inside $r\lesssim 2~{}r_{c}$. The radially cumulated tracer mass decreases with an e-folding decay timescale going from $10$ to $30$ orbits over the duration of the simulation. This tracer dispersal timescale is not adequately defined because the geometry of the tracer distribution keeps evolving in time. Different regions of the envelope are recycled on different timescales: material in the midplane is evacuated faster than material near the poles, and the flow does not efficiently homogenize the tracer distribution near the core. We discuss the use of another dispersion timescale in section 4.5.

4 Results: parameter exploration

We now explore the sensitivity of the different diagnostics presented in the previous section when varying the core mass and size via the parameters $H$ and $B$.

According to equation (6), $B/H=m_{c}/m_{\rm th}$, so that increasing (decreasing) $B$ and $H$ in the same proportions corresponds to decreasing (increasing) the core radius without changing its mass relative to the thermal mass or the mass of its host star. On the other hand, changing $B$ alone while keeping $H$ fixed corresponds to varying the mass of the core without changing its size relative to the disc scale height $h$.

Figure 10 gathers the poloidal maps of time and azimuthally-averaged density and mass flux in our series of simulations listed in Table 1. The top four panel have $H=16$ and increasing values of $B$, while the bottom two panels have $H=32$. We proceed to compare them according to the mass and size of the core.

4.1 Varying the core mass

We start by comparing simulations having the same $H$ to examine the influence of the core mass on the envelope properties. A series of runs with $H=16$ and $B=8$, $16$, $32$ and $64$ shown in Figure 10a-d spans $m_{c}/m_{\rm th}=0.5$, $1$, $2$, and $4$, thus exploring the transition from low mass cores (previously studied by, e.g., Ormel et al. 2015b) to high mass cores (e.g., Bate et al., 2003; Szulágyi et al., 2016) using a single numerical setup.

First, we look at the density distribution around the core. At the lowest mass sampled (H16B8, $m_{c}=0.5m_{\rm th}$), the iso-density contours are rather close to being spherical with only a small amount of oblateness that increases with radius, see Figure 10a. This deviation from spherical symmetry is caused by a weak rotational support in the envelope, which will be examined in more detail in section 4.3. As $m_{c}$ is increased, the density distribution becomes compressed toward the midplane, see the fiducial run H16B16 with $m_{c}=m_{\rm th}$ on Figure 10b. As $m_{c}$ increases further to $m_{c}=2~{}m_{\rm th}$ (H16B32, Figure 10c), more mass becomes concentrated in an oblate overdensity near the midplane. Eventually, in the run H16B64, where the core is four times more massive than $m_{\rm th}$, the density contours are no longer focused on the core: there is a local density maximum at $r\approx 1.8r_{c}$ in the midplane, disjoint from the core surface. This is only possible if rotational support dominates the radial momentum balance, i.e. $v_{\varphi}/v_{\mathrm{K}}\gtrsim 1$. These are unambiguous signatures of a forming circumplanetary disc.

Next, we look at the poloidal flow pattern in this set of simulations. As mentioned in section 3.1, the gas inflow towards the core occurs at high latitudes regardless of the core mass. Some of this flow gets accreted by the core, while a fraction gets diverted towards the midplane and is then funneled outwards, returning to the bulk of the disc for all cores with $m_{c}\lesssim 2m_{\rm th}$. This circulation pattern featuring an outflow of gas in the midplane is very different from the picture of equatorial accretion established in the 2D case (e.g., Kley, 1999). The midplane gas outflow present in the circumplanetary region for all planetary masses is likely to oppose or even prevent the accretion of small particles (pebbles) by the core (Kuwahara et al., 2019).

The situation is a little different only for the most massive core of run H16B64 ($m_{c}=4m_{\rm th}$), as can be see in Figure 10d. Here the inflow has too much angular momentum to reach the core and is deflected toward the midplane by the ‘centrifugal barrier’ even at intermediate latitudes. The resultant pileup of mass in a disc-like structure prevents the high-latitude flow from being reoriented outward in the equatorial plane. Instead, the gas keeps piling up in this circumplanetary disc whose mass increases on orbital timescales. We find that in this simulation the volume-averaged density (inside $r<h$) grows linearly in time as $\left\langle\rho\right\rangle/\rho_{0}\sim 7.0\,\Omega t/2\pi$. Although this is four times faster than in the fiducial run H16B16 (see section 3.2), the time required for the envelope to reach the mass corresponding to hydrostatic equilibrium is now much larger111The mass of a hydrostatic envelope grows as $\int_{r_{c}}^{h}4\pi\rho_{0}r^{2}e^{r_{\mathrm{B}}/r}\mathrm{d}r$, which is faster than exponential with the core mass..

One can also see that variation of $m_{c}$ strongly affects the characteristics of the near polar gas inflow. In the low-mass case of run H16B8 (Figure 10a) the flow remains subsonic in the whole domain. However, as pointed out in section 3.1, a transsonic surface appears around polar axis in the fiducial run H16B16. The inflow along this axis gets eventually arrested in a standing shock not too far from the core surface, and mass accumulates there. As $m_{c}$ increases, the sonic surface extends to lower latitudes, and the standing shock at the poles moves closer to the core surface. Finally, in the run H16B64, supersonic accretion streams span a wide cone about the polar axis. At high latitudes, the inflow reaches Mach $\mathrm{M}\approx 9.4$ and shocks at $z\approx 1.3~{}r_{c}$, nearly on the core surface. The post-shock medium is only barely resolved with seven grid cells in this case. The flow is subsonic only very near the core surface and also in the (detached) thick torus around the forming circumplanetary disc.

The strength of the polar shocks, i.e. the jump in Mach number $\mathrm{M}$, can roughly be estimated by assuming that the flow encounters the shock after falling freely from infinity. At the surface of the core, the corresponding Mach number would be $\mathrm{M}=(2Gm_{c}/r_{s})^{1/2}/c_{s}=(2B)^{1/2}$. This estimate works reasonably well when correcting for the distance between the core surface and the shock front (see section 4.2). The maximum values of the polar Mach number are listed in Table 1.

4.2 Varying the core size

We now compare runs having the same ratio of $B/H$ to isolate the effect of the core size $r_{c}$ on the envelope properties.

Figure 10e represents the sub-thermal mass regime for a small core with $H=32$; it should be compared to Figure 10a, also featuring a core with $m_{c}=0.5~{}m_{\rm th}$ but with $H=16$. Despite the general similarities of the flow pattern (e.g., the orientation of the mass flux circulation), one notices several differences. First, a transsonic surface emerges in the polar regions of run H32B16 even though it is in the sub-thermal mass regime. Second, the iso-density contours are mildly pinched near the midplane. This is explained by the significant level of rotational support in this case, with a maximum $v_{\varphi}/v_{\rm K}\approx 31$ per cent; in run H16B8 this ratio does not exceed two percent. Third, the near-midplane outflow is now confined to a narrower range of latitudes around the equatorial plane of the simulation.

Figure 10f shows the poloidal density and mass flux maps of run H32B32, for which the core radius is twice smaller than in the reference run H16B16 (see Figure 10b), while $m_{c}=m_{\rm th}$ for both. We notice several differences between the two cases. The density gradient, as measured by $\left|\partial\log\rho/\partial\log r\right|$ in the equatorial plane, is lower than in H16B16; to maintain radial momentum balance, the (maximum) fraction of rotational support $v_{\varphi}/v_{\mathrm{K}}\approx 85$ per cent is consequently larger in run H32B32 compared to the $38$ per cent of H16B16. Also, the sonic surface encompasses a wider cone around the core in run H32B32, i.e. a larger fraction of the inflow becomes supersonic before reaching the core.

On the polar axis itself, the time-averaged flow of run H32B32 reaches a Mach number $\mathrm{M}\approx 5.9$ before shocking on the densest parts of the envelope. The stationary polar shocks are located at $z/r_{c}\approx 1.7$ (14 grid cells), closer to the core than in the reference $H=16$ case. The free-fall velocity at this height is $\sqrt{2Gm_{c}/(1.7r_{c})}$, corresponding to $\mathrm{M}\approx 6.1$ for a particle initially at rest away from the core. The agreement with the Mach number measured in simulation, $\mathrm{M}\approx 5.9$ (see Table 1), suggests that the inflow starts feeling a pressure gradient only near the shock, in the deepest layers of the envelope.

Despite these differences, the overall progression of the envelope properties with $m_{c}$ described in section 4.1 for the case $H=16$ remains qualitatively the same for $H=32$, i.e. for twice smaller cores. Higher masses still result in stronger polar shocks, a higher degree of the rotational support of the envelope and more disc-like density distributions.

4.3 Rotational support

We now provide a more systematic overview of the role played by the rotational support in determining the structure of the core envelope. Figure 11 shows the relative measure of rotational support $v_{\varphi}/v_{\mathrm{K}}$ in the equatorial plane of our series of our isothermal simulations; the maximum values of $v_{\varphi}/v_{\mathrm{K}}$ are gathered in Table 1. This ratio $\approx 1$ per cent is negligible in H16B8, for which the envelope is almost entirely supported by pressure gradients against gravity. The fraction of rotational support grows to approximately $35$ per cent in runs H16B16 and H32B16. It reaches $85$ per cent in runs H16B32 and H32B32, with the two curves being almost superimposed. As expected from the density contours of Figure 10 (d), H16B64 achieves full rotational support close to the core. Even in this case, the ratio $v_{\varphi}/v_{\mathrm{K}}$ decreases with radius in the envelope.

This figure shows that the degree of rotational support in the midplane depends primarily on $B$ and only very weakly on $H$. The definition (8) then implies that at a given orbital frequency ($\Omega$) and for a given disc temperature ($c_{s}$), the degree of rotational support of the entire envelope is a function of $m_{c}/r_{c}$. One might expect that for $r\gg r_{c}$ the envelope structure should not depend on the core size anymore as long as $r_{c}$ is smaller than any other characteristic scales of the problem, such as $r_{\rm B}$ and $h$. In the two-dimensional simulations of Paper I, we examined the dependence of $v_{\varphi}/v_{\mathrm{K}}$ on the size of the core and demonstrated that it should indeed vanish in the limit of $H\gg 1$ (i.e., $r_{c}\ll h$). In this limit it is natural to expect $v_{\varphi}/v_{\mathrm{K}}$ to depend only on the core mass, in which case it would be sensitive to $B/H=m_{c}/m_{\rm th}$ and not on $B$ alone. However, runs H16B8 and H32B16 with the same $B/H=0.5$, as well as H16B16 and H32B32 with $B/H=1$, have very different $v_{\varphi}/v_{\mathrm{K}}$ profiles in Figure 11. The core mass is therefore not the only parameter controlling rotational support in these simulations.

On the other hand, we found in Paper I that the convergence rate of $v_{\varphi}/v_{\mathrm{K}}$ with the core size is rather slow. In practice, none of the simulations presented in Paper I or in the current paper effectively reaches this regime with $H\leq 32$. Figure 11 demonstrates that $v_{\varphi}/v_{\mathrm{K}}$ depends on both $m_{c}$ and $r_{c}$ (through $B\equiv Gm_{c}/c_{s}^{2}r_{c}$ mainly), so we must conclude that the cores considered in this study are still large enough to have an effect on their envelope properties. This also implies that the rotational support of the atmospheres of super-Earths discovered by Kepler could be noticeably affected by the finite size of the core.

4.4 Mass accretion

Here we systematically examine how the mass accretion rate onto the core changes between the different simulation runs. Mass accretion, i.e. the permanent trapping of gas from the disc in the vicinity of the core, is regulated by the energy dissipation of the incoming flow near the core, allowing it to become gravitationally bound. The mass of an isothermal envelope is limited from above only by its hydrostatic value, for which rotational support vanishes. In practice, the steady state mass of the envelope is never reached in our simulations because of their relatively short integration time.

We display the mass accretion rates measured in our isothermal simulations in Figure 12. These accretion rates are calculated by integrating the mass flux over spherical shells of different radii centered on the core; as such they account for both the inflow at the poles as well as the outflow in the equatorial regions. The accretion rates measured this way decrease by less than $0.5$ per cent over the duration of each simulation.

Among the runs presented in this paper, only H16B8 has a negligible accretion rate. Incidentally, it is also the only run where the polar inflows do not become supersonic (see Table 1 and Figure 10). This supports the idea that the polar shocks are the primary cause of energy dissipation and mass trapping in the envelope, as opposed to torques in the equatorial plane (cf. Paper I). The polar inflows of run H16B8 are channelled toward the midplane with no physical dissipation that would allow the incoming gas to become bound to the core. As a result, the mass accretion rate is compatible with zero in this simulation. An analogous (non-accreting) situation was presented by Ormel et al. (2015b) with $B=10$ and $H\gg B$.

Moving on to runs with supersonic gas inflows, we find that during the initial stage of envelope accumulation (captured in this study) all such simulations maintain a universal accretion rate

[TABLE]

irrespective of the parameters $B$ or $H$ as long as $B\gtrsim 16$, see Figure 12. Equation (13) shows that in runs with supersonic inflows, the accretion rate $\dot{M}$ scales linearly with the mass of the core and is independent of the thermodynamic conditions in the disc ($\dot{M}$ is essentially independent of $c_{s}$). This is different from the standard Bondi accretion rate, for which $\dot{M}\propto m_{c}^{2}c_{s}^{-3}$ (Bondi, 1952), because in our case the background state away from the core (the sheared flow of the disc) is different from the static and homogeneous background of Bondi (1952).

To better understand our $\dot{M}$ behavior we use the results of Krumholz et al. (2005), who studied Bondi-like accretion in a 3D isothermal flow with non-zero vorticity — a setup very similar to ours except for the lack of density stratification in the vertical direction, and the absence of the central star (which introduces a Hill sphere and horseshoe orbits into the problem). By demanding the streamline deflection in the vicinity of the accretor to be strong enough to result in shocks and binding of the flow, they arrive at an estimate $\dot{M}=C\rho_{0}r_{\mathrm{B}}hc_{s}$ in the limit of the flow vorticity $\omega\gtrsim c_{s}/r_{\mathrm{B}}$, see their equations (28)-(29). In our case $\omega\sim\Omega$, so this limit corresponds to $r_{\mathrm{B}}\gtrsim h$, i.e., close to the regime in which the scaling (13) was obtained. The coefficient $C$ in Krumholz et al. (2005) involves additional logarithmic dependence on $\omega r_{\mathrm{B}}/c_{s}\gtrsim 1$ and is generally larger than the proportionality constant in (13). We attribute this difference to a somewhat different setup: unlike Krumholz et al. (2005), in our case the flow can reach the core only within a limited azimuthal range due to the barrier imposed by the tidal potential.

If the $\dot{M}\propto m_{c}$ behavior given by (13) persists on longer timescales, and if the core is massive enough for its hydrostatic envelope mass to exceed $m_{c}$, then an exponential growth of the envelope may be possible when it becomes self-gravitating. However, this conclusion applies only to the (rather unrealistic) case of a purely isothermal envelope for which the behavior (13) was established. Moreover, complications such as gap formation are likely to limit mass supply from the disc into the feeding zone of the core on long timescales, slowing down the exponential growth eventually (Ginzburg & Chiang, 2019).

4.5 Envelope recycling

We also look at the evolution of the recycling characteristics of the flow within the envelope as a function of $m_{c}$ and $r_{c}$. As in section 3.4, we let every simulation settle to a quasi-steady state over ten orbits, and then we inject a passive tracer with constant concentration $n=1$ inside the Bondi sphere. We let the flow evolve for ten additional orbits and draw the resulting concentration profiles on Figure 13. Larger final concentrations can be interpreted as longer recycling timescales.

Figure 13 shows a pattern of $\langle n(r)\rangle$ very similar to that in Figure 11. After ten orbits, run H16B8 has an average tracer concentration larger than $10$ per cent inside $r/r_{\mathrm{H}}<1/4$, corresponding to the longest recycling timescale. The average concentration profiles of runs H16B16 and H32B16 are nearly superimposed and decrease smoothly with radius, down to $10^{-5}-10^{-4}$ at large radii $r\sim h$. In the remaining runs, the tracer concentration drops by a factor $10^{3}-10^{4}$ outside the dense shell surrounding the core, with runs H16B32 and H32B32 having very similar $\langle n(r)\rangle$ profiles. The polar shocks seem to delimit regions of low (inside) and high (outside) recycling efficiency. The absence of a concentration plateau with $\langle n(r)\rangle=1$ at small radii means that in these runs recycling operates down to the core surface. This progression clearly demonstrates that the profiles of $\langle n(r)\rangle$ depend primarily on the value of $B$ and only weakly on $H$ (as in Figure 11). Thus, the finite size of the core, and not only its mass, has an important effect on the recycling properties of the flow.

We mentioned in section 3.4 that the e-folding decay timescale of the cumulative tracer mass may not be well suited to characterize the recycling efficiency of the flow. It is worth looking for additional diagnostics based on the time-averaged flow variables alone (e.g., via streamlines integration, Ormel et al., 2015b). Here we construct one-dimensional recycling diagnostics similar to those of Kurokawa & Tanigawa (2018) in that they do not rely on a tracer fluid. Let $\tilde{X}\equiv X-\left\langle X\right\rangle_{\mathcal{S}}$ be the deviation of a quantity $X$ relative to its shell average at a given radius, and $\|X\|^{2}\equiv\oint_{\mathcal{S}}X^{2}\,\mathrm{d}S/4\pi r^{2}$ be the squared norm of $X$ on the sphere. Integrating (12) over the volume $\mathcal{V}(r)$ of a sphere $\mathcal{S}(r)$ of radius $r$ centered on the core, the enclosed tracer mass evolves according to

[TABLE]

The first product in parenthesis involves the net mass flux across the sphere $\left\langle\rho v_{r}\right\rangle$. In the second term, $\|\tilde{n}\|$ characterizes the amplitude of variation of the tracer concentration at the surface of the sphere, while $\widetilde{\rho v_{r}}$ is the dispersion of the radial mass flux on the sphere relative to its shell average (cf. Figure 9 of Fung et al., 2015). The correlation between $\tilde{n}$ and $\widetilde{\rho v_{r}}$ on the sphere is measured by $-1\leq\tilde{c}\leq+1$. This correlation remains positive at every radius, but it evolves in time with the tracer distribution in space.

We define recycling as the replacement of the mass inside a given volume by velocity streams passing through this volume. Recycling conserves the total mass inside the volume and can take place in stationary flows ($\left\langle\rho v_{r}\right\rangle=0$). For these reasons, we discard the first product in the parenthesis of (14). The tracer concentration allows us to track fluid motions, but the recycling properties of the flow should be accessible from the velocity field alone. The only term independent of the tracer distribution is $\left\|\widetilde{\rho v_{r}}\right\|$; this term measures the circulation of mass in and out of the sphere, including both the polar inflows and equatorial outflows.

Figure 14 shows the radial profiles of the mean mass flux dispersion $\left\|\widetilde{\rho v_{r}}\right\|$ in our simulations. The ordering of these curves is the inverse of Figure 13: H16B8 is much lower than the other curves, H16B16 and H32B16 are superimposed, and the three remaining curves feature a strong mass flux dropping only near the core surface. This pattern of dependence primarily on $B$ and not $H$ is again reminiscent of Figures 11 & 13. Increasing $B$ enhances the mean radial mass flux dispersion deeper in the envelope, down to the edge of the tracer shells visible in Figure 13. For the most massive cores, the mass flux $\left\|\widetilde{\rho v_{r}}\right\|$ saturates at a sonic accretion/decretion value $\sim\left\langle\rho\right\rangle_{\mathcal{S}}c_{s}$, where $\left\langle\rho\right\rangle_{\mathcal{S}}$ is the shell-averaged density. The correspondence between Figure 13 and Figure 14, in terms of ordering and transition radii, supports that $\left\|\widetilde{\rho v_{r}}\right\|$ essentially captures the recycling efficiency of the flow.

Assuming that the entire envelope is recycled by the mass flux across its boundary, a recycling timescale can be defined as222The prefactor $2$ makes this definition exact for a velocity flux tube having a constant mass flux across its section.

[TABLE]

where we use the norm $\left\|\widetilde{\rho v_{r}}\right\|$ (and its aforementioned saturation at $\left\langle\rho\right\rangle_{\mathcal{S}}c_{s}$) instead of the absolute value $\left|\widetilde{\rho v_{r}}\right|$ to make a link with (14) and Figure 14.

We found that the recycling timescale $\tau(r)$ determined this way varies from $30$ to $200$ orbits inside $r<r_{\mathrm{H}}$ in run H16B8, in relatively good agreement with the tracer dispersion timescale determined from Figure 13. However, in the fiducial run H16B16 we find $\tau>300$ orbits for $r<r_{\mathrm{H}}$. This is more than ten times longer than the tracer dispersion timescale (see Figure 13), since in every simulation run except H16B8, the volume-averaged concentration has decreased by $\gtrsim 50$ per cent over ten orbits. This discrepancy results from the inefficient recycling of the shocked shell surrounding the core, whence the numerator of (15) is largely over-estimated. The volume actually recycled should not include the quasi-static, tracer-rich shells on top of the core. The quantitative agreement between $\tau(r)$ and the tracer dispersion timescale should therefore be limited to low-mass cores, for which the entire envelope is recycled by the mass flux crossing its boundary.

5 Discussion

5.1 Comparison of 2D and 3D simulations

In Paper I, we studied the characteristics of the envelopes of embedded planetary cores in a 2D framework, which has also been done before by Miki (1982), Kley (1999), and Ormel et al. (2015a). Based on our current results, we now provide a systematic comparison between the 2D and 3D simulation outcomes. We discuss different envelope characteristics, highlighting both the similarities and differences of the 2D and 3D setting, and referring to existing literature.

The easiest comparison one can make between 2D and 3D setups is in the equatorial plane of the flow. In 3D, the equatorial density distribution has a morphology similar to 2D, although with weaker shocks (Bate et al., 2003). The equatorial velocity field can still be partitioned into distinct regions — horseshoe, bulk of the disc, and the envelope itself (Fung et al., 2015).

For a given core size and mass, 3D envelopes are more pressure-supported in the equatorial plane than their 2D equivalent. For example, run H16B8 has maximum $v_{\varphi}/v_{\mathrm{K}}\approx 70$ per cent in 2D, while in 3D we find maximum $v_{\varphi}/v_{\mathrm{K}}\approx 38$ per cent. This is due to vertical motions bringing mass from the polar cones without accumulating much vorticity (Ormel et al., 2015b). Indeed, azimuthal averaging around the core reveals large vertical motions (see Figure 10) as seen in global simulations of embedded cores (Kley et al., 2001; Bate et al., 2003). The mean flow is oriented toward the core at high latitudes and radially outward in the midplane. In 2D vertical motions are naturally absent, and the flow in the equatorial plane around the core can flow inward under gravitational torques, opposite to the 3D isothermal case.

In 2D mass accretion onto the core must necessarily be accompanied by the loss of angular momentum of the fluid orbiting the core; otherwise, the “centrifugal barrier” (Ormel et al., 2015a) would not allow accretion to happen. As shown in Paper I, the angular momentum transport is typically effected by the spiral shocks emerging in the circumplanetary disc, although other possibilities exist as well (e.g., by viscous torques, Kley, 1999; Lubow et al., 1999).

In our 3D runs the situation is very different. There we find that polar material falls nearly freely towards the core. Depending on the mass of the core, these vertical motions can reach supersonic velocities before shocking on the densest parts of the envelope (Ayliffe & Bate, 2009b; Szulágyi & Mordasini, 2017). Mass accretion is allowed by the dissipation at the standing isothermal shocks in the polar regions (Tanigawa et al., 2012), which allow the incoming mass to become gravitationally bound to the core. In this way the accreted material thus bypasses the centrifugal barrier, or the “bottleneck” of a circumplanetary disc (Rivier et al., 2012).

For these reasons, the mass accretion rate $\dot{M}$ that we find in the 3D case ends up being quite different from its 2D analog. Indeed, we can rewrite (13) as $\dot{M}\sim\Omega\Sigma a^{2}(m_{c}/m_{\star})(a/h)$, where $\Sigma\sim\rho_{0}h$ is the disc surface density. This form highlights the difference with the behavior suggested in Tanigawa & Tanaka (2016) and Lee (2019) who, based on the 2D isothermal simulations of Tanigawa & Watanabe (2002), proposed the following scaling: $\dot{M}\sim\Omega\Sigma a^{2}(m_{c}/m_{\star})^{4/3}(a/h)^{2}$. Such a difference with the $\dot{M}$ behavior found in our 3D isothermal simulations (and interpreted using the analysis of Krumholz et al. 2005) may have affected the conclusions of these studies.

We also find a correlation between the strength of shocks and turbulent variability, mostly excited near the polar axis. Although the time-averaged flow features a smooth sonic surface, the shock surface is prone to variability on short space and time scales. Our 3D isothermal simulations confirm the presence of significant turbulent variability on top of a secular evolution of the envelope Ormel et al. (2015b); Nelson & Ruffert (2013).

All of our 3D simulations feature prograde inner envelopes, in agreement with previous numerical studies (e.g., Bate et al., 2003; Wang et al., 2014). In 2D setup this property naturally results from the conservation of vortencity (see Miki, 1982, and Paper I). However, as we show next, strict vortensity conservation is apparently not necessary to allow formation of the prograde envelopes in a 3D setup.

5.2 Vortensity conservation

Let $\bm{\omega}\equiv\nabla\times\mathbf{v}$ be the flow vorticity; the potential vorticity, alias vortensity, is defined in a co-rotating frame as $\bm{\varpi}\equiv\left(\bm{\omega}+2\mathbf{\Omega}\right)/\rho$. The background vortensity in the midplane of a Keplerian shear flow (5) is $\varpi_{0}\equiv\Omega/(2\rho_{0})$.

Vortensity is closely related to rotational support in the radial momentum balance (see Ormel et al., 2015a, and Paper I). Its evolution in barotropic flows obeys

[TABLE]

In the inviscid 2D case, the vertical component of vortensity $\varpi_{z}$ is conserved in the Lagrangian sense as the right hand side of equation (16) vanishes identically. This conservation of $\varpi_{z}$ is violated when shocks are present in the flow. In Paper I we demonstrated vortensity production at the spiral shocks launched by the core gravity in a 2D disc. The vortensity perturbation $\delta\varpi=\varpi_{z}-\varpi_{0}$ can be positive or negative depending on the geometry of the shock front. However, away from the shocks $\varpi_{z}$ is strictly conserved in the 2D simulations of Paper I.

The situation is different in 3D. Figure 15 shows the relative deviations of the vertical vortensity $\delta\varpi_{z}/\varpi_{0}\equiv(\varpi_{z}-\varpi_{0})/\varpi_{0}$ in the equatorial plane of our fiducial run H16B16 ($m_{c}=m_{\rm th}$). The vortensity itself remains positive everywhere despite order-unity variations, which explains the prograde rotation of the flow around the core in 3D. However, the vortensity deviation $\delta\varpi_{z}$ exhibits sign variations in different parts of the domain.

$\delta\varpi_{z}$ is negative and roughly constant in the innermost region $r\lesssim 2r_{c}$; every simulation presented here has such an inner envelope with a reduced vortensity. This feature arises because the polar inflows coming from high latitudes bring mass along the direction of unperturbed vortex tubes. Such mass accretion occurs without substantial accumulation of vorticity, reducing the vortensity in this region and causing $\delta\varpi<0$. Mathematically, the negative $\delta\varpi_{z}$ arises because the right hand side of (16) is non-zero in 3D and is dominated by the $\varpi_{z}(\partial v_{z}/\partial z)<0$ term in this part of the flow.

Further away from the core, but still inside the Hill sphere and extending out of the envelope, positive vortensity deviations appear in a non-axisymmetric spiral pattern, with $\delta\varpi_{z}$ alternating sign at a given radius. Note that these vortensity deviations are not superimposed on the spiral density waves of Figure 1. They again emerge due to the non-zero right-hand side of (16). Because the (vertical) inflow velocity $|v_{z}|$ decreases with the cylindrical radius $R$, the toroidal vorticity $\omega_{\varphi}\simeq-(\partial v_{z}/\partial R)>0$ in a wide cone about the $\theta=0$ axis, and the vortensity $\varpi_{\varphi}>0$ as well. As shown on Figure 1, the flow is not axisymmetric around the core; in particular, the vertical velocity $v_{z}$ varies with the azimuthal angle $\varphi$. As a result, the toroidal vorticity $\omega_{\varphi}$ is sheared into a vertical one by the azimuthally-varying vertical flow. Mathematically, the right hand side $(\varpi_{\varphi}/R)(\partial v_{z}/\partial\varphi)$ is non-zero and alternates sign as a function of the angle $\varphi$.

This vertical vortensity produced in the Hill sphere is then passively advected by the outward midplane flow into the bulk of the disc. This rather non-trivial sequence of processes explains the variations of $\delta\varpi_{z}$ in our 3D isothermal simulations.

5.3 Finite core radius effects

Under the simplifications of this study, the problem admits two free length scales: the pressure scale $h$ and the Bondi radius $r_{\mathrm{B}}$. It is commonly assumed that the properties of the planet-disc interaction are asymptotically independent of the core size in the limit of small cores ($H\equiv h/r_{c}\gg 1$). Physically, this means that the flow properties should be set by the core mass alone. Because of the relation (6) one would then expect $B/H=r_{\mathrm{B}}/h$, and not $B$ or $H$ individually, to determine the envelope characteristics.

Following Ormel et al. (2015a), in Paper I we examined this assumption via 2D simulations and found that the core radius controls the properties of the entire envelope at least for $H\leq 32$. Using a simplified 1D model based on the assumption of vortensity conservation, we then found that rotating envelopes around smaller cores ($H\sim 10^{3}$, $B/H\sim 10^{-2}$) could still have a non-vanishing sensitivity to the core radius.

We make the same observation in our 3D isothermal simulations having $H\leq 32$. Simulations with the same ratio of $B/H$ produce different outcomes depending on the size of the core relative to its Bondi radius. In particular, the degree of rotational support $v_{\varphi}/v_{\mathrm{K}}$, which determines whether a circumplanetary disc forms, depends primarily on $B$ and only moderately on $H$, see section 4.3. Similarly, the Mach number of the polar inflows depends primarily on $B$ (section 4.1). By controlling the strength of the polar shocks, the core radius also influences turbulent variability and mixing properties of the flow (section 4.5). This sensitivity of the envelope characteristics to $r_{c}$ is analogous to what we find in 2D, even though arguments based on vortensity conservation no longer apply (see section 5.2).

On the other hand, we find that the mass accretion rate onto the core depends only on the mass of the core to a good accuracy, see Figure 12 and Equation 13. We emphasize that this scaling applies only in our twenty orbit-long simulations, when the mass of the envelope is much lower than its equivalent hydrostatic value (see section 4.4). Nevertheless, this dependence supports the idea that mass accretion rates can be reliably measured in inviscid isothermal simulations even when the core radius is not spatially resolved (see Section 5 in Machida et al., 2010), as long as accretion is mediated by polar shocks.

5.4 Recycling properties

The poloidal circulation in 3D allows material from the envelope to be more efficiently replaced by fresh disc material than in 2D. Indeed, the radial distribution of passive tracer that we find in 2D in Paper I typically exhibits a rather extended plateau near the core with unperturbed tracer distribution, $\langle n\rangle=1$. This plateau can be identified with the rotationally-supported region around the core, which forms a closed system.

In the 3D case, such a concentration plateau, if present, almost always features $\langle n\rangle$ below unity. This means that recycling operates at some level down to the very surface of the core. The overall morphology of the recycling flows connecting the envelope and the disc is, obviously, also different between 2D (where everything is confined to an equatorial plane) and 3D. In the latter case, the midplane is recycled on orbital timescales, whereas recycling takes several tens of orbits in the polar regions. For sufficiently massive cores inducing polar shocks, the post-shock medium is isolated from the rest of the envelope. Turbulent mixing is inefficient in this innermost shell, so the fluid undergoing the most dissipation is also the least efficiently recycled.

In 2D, the radial distribution of a passive tracer fluid can, to zeroth order, be considered as depending primarily on $m_{c}/m_{\rm th}$ — the non-recycled region is quite extended for $m_{c}\lesssim m_{\rm th}$, but shrinks substantially for $m_{c}>m_{\rm th}$. In 3D, as mentioned in section 5.3, the recycling characteristics of the flow depend not only on the core mass but also on its radius.

As opposed to the non-isothermal simulations of Cimerman et al. (2017), Lambrechts & Lega (2017), and Kurokawa & Tanigawa (2018), in our runs envelope recycling operates efficiently down to the core surface in the absence of shocks (run H16B8) and down to the shocks otherwise. The reduced recycling efficiency reported by Kurokawa & Tanigawa (2018) reflects the formation of a circumplanetary disc in the limit of short cooling time (see their Figure 4). In the simulations of Cimerman et al. (2017) and Lambrechts & Lega (2017), an isolated inner envelope appears due to the efficient gas cooling near the core. The resulting contraction should eventually stop when recycling balances cooling in the entropy budget of the envelope (Cimerman et al., 2017).

5.5 Transition from sub-thermal to super-thermal mass regimes

Both Paper I and our present study explored the evolution of envelope characteristics as the core mass was gradually varied from the sub-thermal to the super-thermal regime. Previous studies typically focused on each regime separately (Ayliffe & Bate, 2009a; Ormel et al., 2015a, b; Szulágyi et al., 2016) and employed different setups. This complicates the interpretation of the changes that occur as $m_{c}$ is varied, something that we do naturally here and in Paper I using a single setup for all values of $m_{c}$.

A notable exception is the recent study by Kuwahara et al. (2019), who similarly explored the divide between sub-thermal and super-thermal mass cores using a single 3D isothermal numerical setup. The focus of their study was on exploring the evolution of the characteristics of the midplane outflow as $m_{c}$ was varied; in particular, they demonstrated a significant role played by the midplane outflow in preventing the accretion of solids by the core. Our focus is somewhat different as we concentrate on the gas mass accretion rate, envelope recycling properties, and detailed comparison with the 2D simulations.

We find that both in 2D and 3D the transition from low to high mass cores is accompanied by a considerable change in the envelope properties. Rotational support is insignificant for $m_{c}<m_{\rm th}$, but becomes very important for $m_{c}>m_{\rm th}$ both in 2D and 3D, with a rotationally-supported disc forming around the cores with $m_{c}\gtrsim m_{\rm th}$. In 3D the flow becomes supersonic as the core mass reaches $m_{\rm th}$, starting with polar regions. Radial density distribution, mass accretion, turbulent variability, mixing properties of the flow all change as $m_{c}$ crosses the threshold at $m_{\rm th}$. At the same time, many characteristics of the flow end up depending not only on $m_{c}/m_{\rm th}$ but also on the core size, at least for (relatively low) values of $H$ used in this work and Paper I, see section 5.3.

6 Summary

We have studied the gaseous envelopes surrounding embedded planetary cores via inviscid, isothermal, 3D hydrodynamic simulations. Our numerical setup was designed to explicitly include the core as a spatially-resolved impermeable boundary. The simulations were evolved over twenty orbits, allowing us to study the quasi-instantaneous state of embedded planetary atmospheres. Using this setup we probed the transition from low to high-mass cores ranging from $0.5~{}m_{\rm th}$ to $4~{}m_{\rm th}$, as well as explored the effect of varying the size of the core. Our main conclusions are as follows.

In agreement with other studies we find that a core embedded in a 3D disc drives a meridional circulation pattern in its Hill sphere: the azimuthally-averaged flow moves toward the core at high latitudes and away from the core in the equatorial plane. The equatorial flow circulates with a prograde orientation around the core, despite the lack of strict conservation of the vertical component of the vortencity. 2. 2.

Similar to the 2D case, we find that for the moderate values of $H\leq 32$ explored in this 3D study the finite size of the core, and not just its mass, affects many properties of the flow. The envelope becomes rotationally supported as $B\equiv r_{\mathrm{B}}/r_{c}$ increases, with only a weak dependence on $H\equiv h/r_{c}$. Full rotational support is achieved for $r_{\mathrm{B}}/r_{c}=64$, but the outer envelope is pressure suppported even then. 3. 3.

The polar inflows fall nearly freely towards high-mass cores. The inflow velocity increases with $B$ and becomes supersonic before shocking on the deep envelope for cores with $B\gtrsim 16$. The dissipation associated with these isothermal shocks allows mass accretion in the vertical direction to occur on orbital timescales; the mass accretion rate scales linearly with the core mass but is independent of the core radius or sound speed in the disc; sonic turbulence appears in the same regime, with turbulent Mach numbers of order unity near the shocks. 4. 4.

Different regions of the envelope are recycled on different timescales, ranging from just a couple of orbits in the midplane region to several tens of orbits near the poles. The shocked gas on top of the core remains dynamically bound to the core while the rest of the envelope is recycled on orbital times for $B\geq 32$.

The main simplification of this study is the isothermal assumption for the entire flow. Choosing a different equation of state appears to affect the circulation of the flow in inviscid simulations where the embedded core is spatially resolved (Cimerman et al., 2017). In particular, Kurokawa & Tanigawa (2018) report that recycling can be largely suppressed in adiabatic flows due to a buoyancy barrier. The inclusion of a more sophisticated treatment of the gas thermodynamics is the natural extension of the current work, which will be presented in a future study.

Acknowledgements

Financial support of this work by the Isaac Newton Trust, Department of Applied Mathematics and Theoretical Physics and STFC through grant ST/P000673/1 is gratefully acknowledged. We thank the referee for making a number of useful comments and suggestions to improve this paper. W.B. thanks Richard Nelson and Pablo Benitez-Llambay for sharing their insight on this topic.