Circumplanetary Disk Dynamics in the Isothermal and Adiabatic Limits

Machida et al. (2008) demonstrated that the local, shearing sheet model of a planet embedded in a Keplerian disk can be described by a set of nondimensional equations that is characterized by a single parameter—namely, the ratio between the planet mass and the "disk thermal mass," which we write as:

$$\begin{eqnarray}&&{q}_{\mathrm{thermal}}=\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{* }\,{H}_{{\rm{p}}}^{3}}=\displaystyle \frac{q}{{H}_{{\rm{p}}}^{3}},\end{eqnarray} \tag{ 1 }$$

where ${M}_{{\rm{p}}}$ is the planet's mass, ${M}_{* }$ is the star's mass, q is the planet-to-star mass ratio, and ${H}_{{\rm{p}}}$ is the disk aspect ratio evaluated at the planet's position. Here, ${q}_{\mathrm{thermal}}$ is related to the fundamental length scales in CPD dynamics, which include the Hill radius ${r}_{{\rm{H}}}$, the Bondi radius ${r}_{{\rm{B}}}$, the scale height of the background disk ${h}_{{\rm{p}}}={H}_{{\rm{p}}}{R}_{{\rm{p}}}$, and the half width of the co-orbital horseshoe region ${x}_{{\rm{s}}}$. The Hill and Bondi radii are:

$$\begin{eqnarray}&&{r}_{{\rm{H}}}={R}_{{\rm{p}}}{\left(\displaystyle \frac{q}{3}\right)}^{1/3},\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{r}_{{\rm{B}}}=\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{c}_{{\rm{s}}}^{2}}={R}_{{\rm{p}}}\displaystyle \frac{q}{{H}_{{\rm{p}}}^{2}},\end{eqnarray} \tag{ 3 }$$

where G is the gravitational constant, ${R}_{{\rm{p}}}$ is the radial location of the planet, and ${c}_{{\rm{s}}}$ is the sound speed of the gas. To provide a sense of scale, we can also write:

$$\begin{eqnarray}&&{q}_{\mathrm{thermal}}=0.7\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10\,{M}_{\oplus }}\right){\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{H}_{{\rm{p}}}}{0.035}\right)}^{-3},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{r}_{{\rm{H}}}=0.031\,\mathrm{au}\left(\displaystyle \frac{{R}_{{\rm{p}}}}{\mathrm{au}}\right){\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{\oplus }}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-1/3},\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{r}_{{\rm{B}}}=0.027\,\mathrm{au}\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10\,{M}_{\oplus }}\right){\left(\displaystyle \frac{{c}_{{\rm{s}}}}{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-2}.\end{eqnarray} \tag{ 6 }$$

For ${x}_{{\rm{s}}}$, Masset et al. (2006) found that it can be separated into two regimes:

$$\begin{eqnarray}{x}_{{\rm{s}}}=\left\{\begin{array}{ll}{R}_{{\rm{p}}}\sqrt{q/{H}_{{\rm{p}}}} & {\rm{if}}{q}_{\mathrm{thermal}}\lesssim 1,\\ 2.5{r}_{{\rm{H}}} & {\rm{if}}{q}_{\mathrm{thermal}}\gtrsim 1.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

The ratio between any two of these four length scales, ${r}_{{\rm{H}}}$, ${r}_{{\rm{B}}}$, ${h}_{{\rm{p}}}$, and ${x}_{{\rm{s}}}$, can be expressed in terms of ${q}_{\mathrm{thermal}}$ and ${q}_{\mathrm{thermal}}$ only. For instance, ${r}_{{\rm{B}}}/{r}_{{\rm{H}}}={3}^{1/3}{q}_{\mathrm{thermal}}^{2/3}$, ${h}_{{\rm{p}}}/{r}_{{\rm{H}}}={3}^{1/3}{q}_{\mathrm{thermal}}^{-1/3}$, and ${x}_{{\rm{s}}}/{r}_{{\rm{H}}}={3}^{1/3}{q}_{\mathrm{thermal}}^{1/6}$ in the ${q}_{\mathrm{thermal}}\lesssim 1$ regime. Hence, ${q}_{\mathrm{thermal}}$ alone should be sufficient to determine the dynamics of the CPD, as long as we restrict ourselves to considering only gravity and hydrodynamics. While it is beyond the scope of this work, we note that other physical parameters such as thermal diffusivity and optical thickness will add extra dimensions to this problem. We focus here on how CPD dynamics depends on ${q}_{\mathrm{thermal}}$, and will verify that simulations with the same ${q}_{\mathrm{thermal}}$ but different q and ${H}_{{\rm{p}}}$ will yield the same results.

We use the graphics processing unit–accelerated hydrodynamics code PEnGUIn (Fung 2015) to simulate planets embedded in disks. PEnGUIn solves the Euler equations:

$$\begin{eqnarray}&&\displaystyle \frac{D\rho }{{Dt}}=-\rho \left({\rm{\nabla }}\cdot {\boldsymbol{v}}\right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{v}}}{{Dt}}=-\displaystyle \frac{1}{\rho }{\rm{\nabla }}p-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{De}}{{Dt}}=-\displaystyle \frac{1}{\rho }{\rm{\nabla }}\cdot p{\boldsymbol{v}}-{\rm{\nabla }}{\rm{\Phi }}\cdot {\boldsymbol{v}},\end{eqnarray} \tag{ 10 }$$

where ρ is the gas density, p is the gas pressure, ${\boldsymbol{v}}$ is the gas velocity, and $e=u+| {\boldsymbol{v}}{| }^{2}/2$ is the specific total (internal + kinetic) energy. In isothermal simulations, we discard Equation (10) and set $p/\rho$ to be globally constant. Here, Φ is the combined gravitational potential of the star and the planet. In spherical coordinates centered on the star, where {R, ψ, θ} represent the radial, azimuthal, and polar coordinates, respectively, Φ can be written as:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Phi }} & = & -\displaystyle \frac{{{GM}}_{* }}{R}-\displaystyle \frac{{{GM}}_{{\rm{p}}}}{\sqrt{{R}^{2}+{R}_{{\rm{p}}}^{2}-{{RR}}_{{\rm{p}}}\sin \theta \cos \psi ^{\prime} +{r}_{{\rm{s}}}^{2}}}\\ & & +\displaystyle \frac{{{GM}}_{{\rm{p}}}R\sin \theta \cos \psi ^{\prime} }{{R}_{{\rm{p}}}^{2}},\end{array}\end{eqnarray} \tag{ 11 }$$

where $\psi ^{\prime} =\psi -{\psi }_{{\rm{p}}}$ is the angular distance from the planet. The third term on the right is the indirect potential due to our frame centering on the star rather than the center of mass. We have assumed that the planet is on a circular orbit in the midplane. The smoothing length ${r}_{{\rm{s}}}$ is included to prevent numerical instability very close to the planet, and could represent the size of the planet's solid core. For example, the radius of the Earth at 1 au is about 1.6% of the Earth's Bondi radius. In this work, we do not assume a size for the core, and instead aim to quantify how our results depend on ${r}_{{\rm{s}}}$. Following Fung et al. (2017), we set ${r}_{{\rm{s}}}$ to be the length of three grid cells (resolution is described in Section 2.1.2; also see Table 1 for values of ${r}_{{\rm{s}}}$).

Table 1.  Model Parameters

Model #	${H}_{{\rm{p}}}$	EOS	${M}_{{\rm{p}}}$	${q}_{\mathrm{thermal}}$	Resolution	${r}_{{\rm{s}}}$ a	${t}_{\mathrm{sim}}$ b
			(${M}_{\oplus }$)		(cells/$\min [{r}_{{\rm{B}}},h]$)	(${r}_{{\rm{B}}}$)	($2\pi {{\rm{\Omega }}}_{{\rm{p}}}^{-1}$)
1	0.035	isothermal	1.4	0.1	64	0.047	21
2	0.035	isothermal	3.5	0.25	64	0.047	21
3	0.035	isothermal	7	0.5	64	0.047	21
4	0.035	isothermal	14	1	64	0.047	21
5	0.035	isothermal	28	2	64	0.023	21
6	0.035	isothermal	56	4	64	0.012	21
7	0.035	adiabatic	1.4	0.1	64	0.047	11
8	0.035	adiabatic	3.5	0.25	64	0.047	11
9	0.035	adiabatic	7	0.5	64	0.047	11
10	0.035	adiabatic	14	1	64	0.047	11
11	0.035	adiabatic	28	2	64	0.023	11
12	0.035	adiabatic	56	4	64	0.012	11
13	0.1	isothermal	333	1	64	0.023	21
14	0.1	isothermal	1333	4	64	0.012	100c
15	0.1	adiabatic	333	1	64	0.023	11
16	0.1	adiabatic	1333	4	64	0.012	11
17d	0.035	isothermal	1.4	0.1	512e	0.006	3f
18	0.035	isothermal	14	1	512e	0.006	3f

Notes.

^aThe smoothing length ${r}_{{\rm{s}}}$ is equal to three times the size of the smallest cells. ^bUnless otherwise stated, results are time-averaged over the last orbit. ^cResults are time-averaged between the 20th to 21st orbit, the same as other isothermal runs, but we then extend it to 100 orbits to study the effects of gap opening (see Section 3.5). ^dModel #17 is also simulated with Athena++ using the same physical parameters but a different numerical setup. See Section 2.2. ^eUnlike other models where the resolution is nearly uniform inside all of ${r}_{{\rm{B}}}$, in models #17 and #18, only within $\sim 0.1\,{r}_{{\rm{B}}}$ is the resolution equivalent to 512 cells $/{r}_{{\rm{B}}}$. ^fIn models #17 and #18, we do not perform time averaging over the last orbit.

Download table as:  ASCIITypeset image

Simulations are performed in a spherical grid centered on the star, as well as in a rotating frame at the planet's orbital frequency ${{\rm{\Omega }}}_{{\rm{p}}}$, fixing the planet in space. This is advantageous because it reduces numerical diffusion caused by advection. The Coriolis force due to frame rotation is not explicitly computed; rather, it is absorbed into the conservative form of the angular momentum equation (Kley 1998).

For the equation of state (EOS), we use the ideal gas law, such that:

$$\begin{eqnarray}&&p={c}_{\mathrm{iso}}^{2}\rho \end{eqnarray} \tag{ 12 }$$

in isothermal simulations, where ${c}_{\mathrm{iso}}$ is the isothermal sound speed, and

$$\begin{eqnarray}&&p=u\rho (\gamma -1)\end{eqnarray} \tag{ 13 }$$

in adiabatic simulations, where γ is the ratio of specific heats (or adiabatic index). We choose $\gamma =7/5$.

Mass and momentum fluxes across simulation cells are conservative, but the total energy is not. In PPDs, the orbital speed is highly supersonic, making the kinetic energy dominant over the internal energy by orders of magnitude. Thus, a conservative scheme in total energy would produce significant noise in the internal energy of the gas. We therefore opt to conserve the internal energy instead. Testing shows that this leads to significantly more stable flows.

In most cases, we run our simulations to 11 or 21 planetary orbits (${t}_{\mathrm{sim}}$ in Table 1). In fact, they all appear to reach quasi-steady states even after just ∼2 orbits, in the sense that our results would not be significantly different even if we had terminated them at just ~2 orbits. Nonetheless, we run them much longer to confirm our results are robust. There are noticeable temporal fluctuations, particularly in the isothermal cases. These fluctuations are not qualitatively important, but can affect our quantitative results. Therefore, unless otherwise stated, all our results are time-averaged over the last orbit.

Even though our simulation grid centers on the star, for analysis it is convenient to use a cylindrical coordinate centered on the planet. Thus, we will also be using {r, ϕ, z} to denote the radial, azimuthal, and vertical coordinates where the planet is at $r=z=0$.

2.1.1. Initial and Boundary Conditions

We assume the initial disk is axisymmetric, is in hydrostatic equilibrium with the star's gravity, and has a power-law profile in the radial direction:

$$\begin{eqnarray}&&\rho ={\rho }_{0}{\left(\displaystyle \frac{R\sin \theta }{{R}_{{\rm{p}}}}\right)}^{-3}\exp \left[-\displaystyle \frac{{{GM}}_{* }}{{{Rc}}_{\mathrm{iso}}^{2}}\left(\displaystyle \frac{1}{\sin \theta }-1\right)\right],\end{eqnarray} \tag{ 14 }$$

where ${\rho }_{0}$ is a background normalization (i.e., the ambient disk density at the planet's orbital radius, not including perturbations by the planet). In the code, ${\rho }_{0}$ is set to 1; the exact value is immaterial because gas self-gravity is neglected.

The power of −3 is chosen such that the gas surface (vertically integrated) density Σ scales as ${R}^{-3/2}$. We also denote ${{\rm{\Sigma }}}_{0}$ as the surface density at the planet's location. We choose ${c}_{\mathrm{iso}}$ to be either $0.035{R}_{{\rm{p}}}{{\rm{\Omega }}}_{{\rm{p}}}$ or $0.1{R}_{{\rm{p}}}{{\rm{\Omega }}}_{{\rm{p}}}$ (see Table 1), which respectively correspond to disk aspect ratios ${H}_{{\rm{p}}}$ of 0.035 or 0.1 at the planet's location. To make the comparison as direct as possible, Equation (14) is used for both isothermal and adiabatic simulations, and gas pressure is also initialized as $p={c}_{\mathrm{iso}}^{2}\rho$ in both cases. In other words, the initial conditions in both the isothermal and adiabatic simulations are completely identical. We note that the sound speed for adiabatic gas is not ${c}_{\mathrm{iso}}$, but rather ${c}_{{\rm{s}}}=\sqrt{\gamma p/\rho }={\gamma }^{1/2}{c}_{\mathrm{iso}}$. To keep our notation simple, the Bondi radius is evaluated as ${r}_{{\rm{B}}}={{GM}}_{{\rm{p}}}/{c}_{\mathrm{iso}}^{2}$ regardless of the EOS.

To establish a hydrostatic disk, initially there is no radial or polar motion, and the azimuthal rotation frequency is:

$$\begin{eqnarray}&&{\rm{\Omega }}=\sqrt{\displaystyle \frac{{{GM}}_{* }}{{R}^{3}}+\displaystyle \frac{1}{R\rho }\displaystyle \frac{\partial p}{\partial R}}.\end{eqnarray} \tag{ 15 }$$

Our simulation domain spans ${R}_{{\rm{p}}}-10\,{r}_{{\rm{B}}}$ to ${R}_{{\rm{p}}}+10\,{r}_{{\rm{B}}}$ radially, the full 2π azimuthally, and $\pi /2-3\,{H}_{{\rm{p}}}$ to $\pi /2$ in the polar direction, which is from the disk midplane to about three scale heights above.

We impose periodic boundary conditions in the azimuthal direction, and reflective boundaries in the polar direction. Reflective boundaries are used in the midplane to enforce symmetry, and at the top to prevent gas from flowing in or out of the simulation box. The radial boundaries are fixed to the initial values. Additionally, we place wave-killing zones next to the radial boundaries to help reduce wave reflections. They are prescribed as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial X}{\partial t}=\displaystyle \frac{X(t=0)-X}{{t}_{\mathrm{kill}}}{\left(1-\displaystyle \frac{| r-{r}_{\mathrm{bound}}| }{{L}_{\mathrm{kill}}}\right)}^{2},\end{eqnarray} \tag{ 16 }$$

where X corresponds to the fluid properties ρ, p, and each component in ${\boldsymbol{v}};$ ${t}_{\mathrm{kill}}$ is the damping timescale; ${r}_{\mathrm{bound}}$ is the position of either the inner or the outer radial boundary; and ${L}_{\mathrm{kill}}\gt | r-{r}_{\mathrm{bound}}|$ is the width of the kill zone. We choose ${t}_{\mathrm{kill}}=0.2\,\pi \,{{\rm{\Omega }}}_{{\rm{p}}}^{-1}$, which is one-tenth of the planet's orbital period, and ${L}_{\mathrm{kill}}=h$, except when ${q}_{\mathrm{thermal}}=0.1$ (models #1, 7, and 17), where we have ${L}_{\mathrm{kill}}=0.1\,h$ instead.

2.1.2. Resolution

We use a nonuniform grid to concentrate resolution near the planet. If we denote L as the maximum distance away from the planet along one of the three coordinates, N as the number of cells within L, and $i$ as the $i\mathrm{th}$ cell away from the planet, then ${x}_{{\rm{i}}}$, the distance from the $i\mathrm{th}$ cell to the planet, is:

$$\begin{eqnarray}&&{x}_{{\rm{i}}}=i{\rm{\Delta }}{x}_{\min }+\left(L-N{\rm{\Delta }}{x}_{\min }\right){\left(\displaystyle \frac{i}{N}\right)}^{a},\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&a=\displaystyle \frac{\mathrm{ln}\left(1-\tfrac{{\rm{\Delta }}{x}_{\max }-{\rm{\Delta }}{x}_{\min }}{L-N{\rm{\Delta }}{x}_{\min }}\right)}{\mathrm{ln}\left(1-\tfrac{1}{N}\right)}.\end{eqnarray} \tag{ 18 }$$

Here, ${\rm{\Delta }}{x}_{\min }$ determines the resolution near the planet, while ${\rm{\Delta }}{x}_{\max }\gg {\rm{\Delta }}{x}_{\min }$ is the cell size farthest away from the planet. When ${q}_{\mathrm{thermal}}\leqslant 1$ and ${H}_{{\rm{p}}}=0.035$, along each of the three directions {R, ψ, θ}, we have $L=\{10{r}_{{\rm{B}}},\,\pi ,\,3{H}_{{\rm{p}}}\}$ and ${\rm{\Delta }}{x}_{\max }=\{{r}_{{\rm{B}}}/4,\,0.1,\,{H}_{{\rm{p}}}/8\}$, and ${\rm{\Delta }}{x}_{\min }$ is either ${r}_{{\rm{B}}}/64$ or ${r}_{{\rm{B}}}/512$ (see Table 1) but is the same in all directions. When ${H}_{{\rm{p}}}=0.1$, L in the radial direction is 0.6 ${R}_{{\rm{p}}}$ instead. When ${q}_{\mathrm{thermal}}={r}_{{\rm{B}}}/h\gt 1$, we use ${\rm{\Delta }}{x}_{\max }=\{h/4,\,0.1,\,{H}_{{\rm{p}}}/8\}$ and ${\rm{\Delta }}{x}_{\min }=h/64$. We use $N=\{192,\,240,\,128\}$ in models #1–16, and $N=\{224,\,272,\,192\}$ in models #17–18; in terms of the total number of cells in the grid, they correspond to $384\times 480\times 128$ and $448\times 544\times 192$, respectively.

A main goal of this work is to resolve and analyze the rotation in the gas around protoplanets. To confirm we can achieve this goal, we take our fiducial model where ${q}_{\mathrm{thermal}}=1$ and ${H}_{{\rm{p}}}=0.035$ (the same as models #4 and #18), and simulate it under different resolutions ranging from 16 to $128\,\mathrm{cells}/{r}_{{\rm{B}}}$. The smoothing length ${r}_{{\rm{s}}}$ is always set to be three times the cell size. Figure 1 plots the rotation curves around the planet from these simulations. At $64\,\mathrm{cells}/{r}_{{\rm{B}}}$, we find that we have reached numerical convergence to within a percent level for the bulk of the Bondi sphere—whereas very close to the planet, $\lesssim 0.1{r}_{{\rm{B}}}$, we start seeing the effects of the smoothed gravitational potential and find slower rotation speeds. Using this test as a guide, we use $64\,\mathrm{cells}/{r}_{{\rm{B}}}$ as our fiducial resolution ($64\,\mathrm{cells}/h$ when ${r}_{{\rm{B}}}\gt h$, as in the cases with models #5, 6, 11, 12, 14, and 16) to study dynamics on the ${r}_{{\rm{B}}}$ scale, and enhance our resolutions by factors up to 8 in some models to study dynamics deep within ${r}_{{\rm{B}}}$.

Zoom In Zoom Out Reset image size

For one particular model, Model 17 in Table 1, we carry out a similar simulation but with a quite different numerical setup using the Athena++ code (Stone et al. 2008). In contrast to the PEnGUIn setup, we adopt a spherical-polar grid centered on the planet. In the radial direction, the grid spacing is uniform in logarithmic space, with 256 cells from ${r}_{\min }=3\times {10}^{-5}$ ${R}_{{\rm{p}}}$ to ${r}_{\max }=0.35$ ${R}_{{\rm{p}}}$. At the inner boundary, a reflecting boundary condition is adopted. At the outer boundary, the fluid variables are fixed to their initial values so that the gas flow is orbiting around the central star. The grid is uniform in the polar direction, with 32 grids from 0 to $\pi /4$. The polar boundary condition (Zhu & Stone 2018) is adopted at the pole, while a reflecting boundary condition is adopted at the midplane. The grid is uniform in the azimuthal direction, with 128 grids from 0 to 2π. With this grid structure, each cell at the disk midplane has the same length in all three directions. At the inner boundary, where resolution is the highest, the edge of each cell has a length of ${r}_{{\rm{B}}}/2378$.

The planet's gravitational acceleration is smoothed by the following function:

$$\begin{eqnarray}&&f=\displaystyle \frac{{\left(r-{r}_{\min }\right)}^{2}}{{\left(r-{r}_{\min }\right)}^{2}+{r}_{\mathrm{smooth}}^{2}},\end{eqnarray} \tag{ 19 }$$

where ${r}_{\mathrm{smooth}}=4\times {10}^{-6}$ ${R}_{{\rm{p}}}$. A density floor of 10⁻¹⁰ times the initial midplane density at R_p is adopted. All other planet and disk setups are the same as the PEnGUIn setup as described in Section 2.1.1.

We run the simulation for 1.1 planetary orbits. If the CPD were rotating at Keplerian speed, it would correspond to $1.37\times {10}^{4}$ orbits at ${r}_{\min }$.

Determining which part of the gas is bound to the planet is nontrivial. For example, the specific energy of gas is not a good indicator, because it is not a conserved quantity. One way to infer boundedness is from kinematics. Gas that flows away from the planet is unbound; otherwise, it is bound and should be considered a part of the CPD.

Figures 8 and 9 show two examples of the azimuthally averaged meridional flow pattern, one each for the isothermal and adiabatic CPDs. Generally, gas flows vertically toward the planet along the poles, and away from the planet near the midplane, as has been shown in numerous previous studies (e.g., Fung et al. 2015; Ormel et al. 2015; Béthune & Rafikov 2019; Kuwahara et al. 2019). In isothermal cases, this inflow is supersonic and can reach $4\sim 5\,{c}_{\mathrm{iso}}$. When adiabatic, on the other hand, the inflow is subsonic and generally does not exceed $0.1\,{c}_{\mathrm{iso}}$. In the midplane, the flow is directed away from the planet at larger distances, but toward the planet closer in. We can therefore use the location where the sign of the midplane radial velocity changes to characterize the sizes of our CPDs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The inward flow occurs at around $\sim 0.1{r}_{{\rm{B}}}$ for isothermal CPDs, and $\sim 0.2{r}_{{\rm{B}}}$ for adiabatic ones. This measurement has some uncertainty because there are substantial temporal fluctuations in the velocity field close to the planet. In particular, the radial velocity can frequently change sign. Nonetheless, when averaged over time, we consistently find these specific locations to be where the radial velocity (centered on the planet) changes direction in all our models with ${q}_{\mathrm{thermal}}\leqslant 1$. For the largest ${q}_{\mathrm{thermal}}$ tested, ${q}_{\mathrm{thermal}}=4$, the isothermal CPD size is still about $0.1\,{r}_{{\rm{B}}}$, but the adiabatic CPD size shrinks and becomes closer to $0.1\,{r}_{{\rm{B}}}$. We expect the CPD size to eventually scale with ${r}_{{\rm{H}}}$ instead of ${r}_{{\rm{B}}}$ as ${q}_{\mathrm{thermal}}$ increases, but our parameter space does not extend far enough to quantify that regime.

Another indicator of boundedness is the vertical density structure. If the gas is bound and has no significant vertical motion, it should settle into vertical hydrostatic equilibrium. The isothermal vertical hydrostatic density profile is:

$$\begin{eqnarray}&&\rho (r,z)={\rho }_{0}(r)\exp \left[\displaystyle \frac{{r}_{{\rm{B}}}}{\sqrt{{r}^{2}+{z}^{2}}}-\displaystyle \frac{{r}_{{\rm{B}}}}{r}\right]\exp \left[-\displaystyle \frac{{z}^{2}}{2{h}_{\mathrm{PPD}}^{2}}\right].\end{eqnarray} \tag{ 22 }$$

In the limit of $z\ll r$, it can be written as:

$$\begin{eqnarray}&&\rho (r,z)\approx {\rho }_{0}\exp \left[-\displaystyle \frac{{z}^{2}}{2}\left(\displaystyle \frac{1}{{h}_{\mathrm{CPD}}^{2}}+\displaystyle \frac{1}{{h}_{\mathrm{PPD}}^{2}}\right)\right],\end{eqnarray} \tag{ 23 }$$

where ${h}_{\mathrm{CPD}}=\sqrt{{r}^{3}/{r}_{{\rm{B}}}}$ is the expected scale height of the CPD if it is vertically settled.

The disk scale height as a function of distance from the planet is shown in Figure 10 for our high-resolution isothermal models (#17 and #18); the same results but at fiducial resolution are shown in Figure 11, with isothermal runs (#1–6) on the left and adiabatic runs (#7–12) on the right. We measure the CPD scale height h using the following definition:

$$\begin{eqnarray}&&0.68=\displaystyle \frac{{\int }_{0}^{h}\rho \,{dz}}{{\int }_{0}^{{z}_{\max }}\rho \,{dz}},\end{eqnarray} \tag{ 24 }$$

where ${z}_{\max }$ marks the top vertical boundary of the simulation. In other words, 68% of the gas mass lies below h. This definition does not explicitly depend on the local sound speed and so can be used consistently in both the isothermal and adiabatic runs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For isothermal CPDs, their vertical profiles follow the hydrostatic solution within $\sim 0.1{r}_{{\rm{B}}}$. This holds for all values of ${q}_{\mathrm{thermal}}$ tested, as shown in the left panel of Figures 11 and 10. Beyond this distance, the disk expands vertically and follows the background scale height ${h}_{\mathrm{PPD}}$ instead. The gas outside $0.1{r}_{{\rm{B}}}$ must therefore be unbound to the planet and passing through the Bondi sphere so rapidly that it barely reacts to the planet's gravity. This is consistent with the CPD sizes inferred from kinematics.

For the adiabatic runs, we measure small dips in h directly above the planets, ranging from $\sim 30 \%$ for ${q}_{\mathrm{thermal}}=0.1$ to a factor of ∼4 for ${q}_{\mathrm{thermal}}=4$. These dips do not extend beyond $\sim 0.2{r}_{{\rm{B}}}$, which is again consistent with our interpretation that gas is unbound beyond that point.

When the local PPD aspect ratio is about 0.1, ${q}_{\mathrm{thermal}}=0.1$ corresponds roughly to 0.1 Jupiter mass. Analyses of gaps in PPDs suggest that planets around this size may be common between 10 and 100 au (Zhang et al. 2018). If these planets are present, our results indicate that their signature on the disk surface should be small if their CPDs are close to adiabatic. This may explain why they are not observed directly, despite their prominent gaps.

In the classical 2D picture, the background Keplerian shear provides the source of angular momentum for CPDs. Gas is accreted by the planet through the L1 and L2 Lagrange points; material enters the Hill sphere with an angular momentum of roughly ${r}_{{\rm{H}}}^{2}{{\rm{\Omega }}}_{{\rm{p}}}$. Setting this equal to the Keplerian angular momentum around the planet $\sqrt{{{GM}}_{{\rm{p}}}r}$, one finds $r={r}_{{\rm{H}}}/3$, implying that the CPD is rotationally supported within $\sim {r}_{{\rm{H}}}/3$ (Quillen & Trilling 1998). Two-dimensional calculations focusing on the effects of tidal truncation produced a similar disk size (Martin & Lubow 2011).

This picture is modified significantly in 3D. Previous studies have found that, in 3D, planets accrete gas from the vertical direction instead. That gas originates directly above the planet and could be co-orbiting with it. A small orbital velocity difference between the gas and the planet would mean the gas has a lower angular momentum than in 2D CPDs.

How low might this angular momentum be? The lower it is, the smaller the rotationally supported region is, and the higher the required resolution becomes. This presents a challenge to our ability to resolve the CPD. In Section 2.1.2, we have demonstrated that our fiducial resolution is converged for the rotation speed on scale $\sim {r}_{{\rm{B}}}$, but it may not be sufficient if the rotationally supported region turns out to be much smaller than ${r}_{{\rm{B}}}$. We shall bear this in mind as we proceed.

The Keplerian rotation around a planet can be expressed in terms of ${c}_{\mathrm{iso}}$ and ${r}_{{\rm{B}}}$:

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{{\rm{K}}}}{{c}_{\mathrm{iso}}}={\left(\displaystyle \frac{r}{{r}_{{\rm{B}}}}\right)}^{-1/2},\end{eqnarray} \tag{ 25 }$$

as long as we use ${c}_{\mathrm{iso}}$ when defining ${r}_{{\rm{B}}}$ (Equation (3)), and the specific angular momentum profile is similarly:

$$\begin{eqnarray}&&\displaystyle \frac{{l}_{{\rm{K}}}}{{r}_{{\rm{B}}}{c}_{\mathrm{iso}}}={\left(\displaystyle \frac{r}{{r}_{{\rm{B}}}}\right)}^{1/2}.\end{eqnarray} \tag{ 26 }$$

We will normalize the speeds and distances of our results by ${c}_{\mathrm{iso}}$ and ${r}_{{\rm{B}}}$ to directly compare simulations with different planet and disk parameters. Figure 12 plots the midplane rotation profiles from the isothermal simulations (models #1–6) in the left panel, and the adiabatic results (models #7–12) are on the right. These profiles are azimuthally averaged in the planet-centered frame, but we note that within ${r}_{{\rm{B}}}$ there generally is little azimuthal variation.

Zoom In Zoom Out Reset image size

3.2.1. Isothermal Disks

Strikingly, within a distance of $\sim 0.2,{r}_{{\rm{B}}}$, all angular momentum profiles from various ${q}_{\mathrm{thermal}}$ converge to a single value in units of ${r}_{{\rm{B}}}{c}_{\mathrm{iso}}$. For $0.1\leqslant {q}_{\mathrm{thermal}}\leqslant 1$, the profiles nearly lie on top of each other. The ${q}_{\mathrm{thermal}}=2$ and 4 simulations have a shorter normalized smoothing length ${r}_{{\rm{s}}}/{r}_{{\rm{B}}}$ (see Table 1), and so they reach higher speeds at very short distances.

The fact that we can achieve higher rotation speeds by decreasing ${r}_{{\rm{s}}}/{r}_{{\rm{B}}}$ suggests that ${r}_{{\rm{s}}}$ still has too large of an effect on the rotation at our fiducial resolution. We therefore increase the resolution inside $\sim 0.2{r}_{{\rm{B}}}$ by a factor of eight, correspondingly reducing ${r}_{{\rm{s}}}$ by a factor of eight, to produce models #17 and #18. Figure 13 plots the angular momentum profiles from those models. They confirm that rotationally supported, Keplerian disks indeed exist around these planets.

Zoom In Zoom Out Reset image size

The angular momentum profiles around subthermal planets can be approximated as a superposition of a constant value and the background shear that one would obtain in the absence of the planet. A formal fit gives:

$$\begin{eqnarray}&&\displaystyle \frac{l}{{r}_{{\rm{B}}}{c}_{\mathrm{iso}}}=\displaystyle \frac{{l}_{\max }}{{r}_{{\rm{B}}}{c}_{\mathrm{iso}}}\,-\,\displaystyle \frac{3{q}_{\mathrm{thermal}}}{4}{\left(\displaystyle \frac{r}{{r}_{{\rm{B}}}}\right)}^{2},\end{eqnarray} \tag{ 27 }$$

where ${l}_{\max }=0.23\,{r}_{{\rm{B}}}{c}_{\mathrm{iso}}$ and the second term on the right corresponds to the background shear. We overplot this profile in Figure 13 to show that it compares well with our empirical profiles.

The value $0.23\,{r}_{{\rm{B}}}{c}_{\mathrm{iso}}$ corresponds to the Keplerian angular momentum at about $0.05\,{r}_{{\rm{B}}}$. This is the size of the Keplerian disk, which we denote as ${r}_{{\rm{K}}}$. The scaling ${r}_{{\rm{B}}}{c}_{\mathrm{iso}}$ can alternatively be expressed as ${x}_{{\rm{s}}}^{2}{{\rm{\Omega }}}_{{\rm{p}}}$, where ${x}_{{\rm{s}}}$ takes the subthermal branch in Equation (7). Therefore, ${l}_{\max }$ may be related to the incoming momentum in the horseshoe orbits, which is consistent with that idea that the CPD is fed by the horseshoe flow (Fung et al. 2015). Even though our parameter space only covers down to ${q}_{\mathrm{thermal}}=0.1$, given the lack of dependence on ${q}_{\mathrm{thermal}}$ in the rotation profile, we expect our results to apply to all subthermal planets. This implies that, under isothermal conditions, even smaller planets (such as the Earth) should have rotationally supported CPDs inside ${r}_{{\rm{K}}}\approx 0.05{r}_{{\rm{B}}}$.

For superthermal planets, the maximum angular momentum in their CPDs is also about $0.23{r}_{{\rm{B}}}{c}_{\mathrm{iso}}$, but the size of the region with this specific angular momentum rapidly shrinks with increasing ${q}_{\mathrm{thermal}}$, to the point that the overall profile significantly deviates from Equation (27). The left panel of Figure 12 shows that when ${q}_{\mathrm{thermal}}=4$, ${l}_{\max }$ is reached at just about $0.15{r}_{{\rm{B}}}$.

In this regime, it is likely that we are beginning to see the transition of the CPD from being limited by the Bondi radius to the Hill radius. Clearly, one should not expect ${l}_{\max }$ to scale with ${r}_{{\rm{B}}}$ indefinitely, or else the size of the ${r}_{{\rm{K}}}$ disk will eventually exceed ${r}_{{\rm{H}}}$. One commonly suggested scaling for superthermal planets is ${r}_{{\rm{K}}}\sim {r}_{{\rm{H}}}/3$. If we take that scaling, then the transition would occur near ${r}_{{\rm{H}}}/3={r}_{{\rm{B}}}/20$, corresponding to ${q}_{\mathrm{thermal}}\sim 10$.

To address the possibility of code bias, we also compare our inferred values of ${l}_{\max }$ and ${r}_{{\rm{K}}}$ to the simulations by Wang et al. (2014), who used the Antares code and its static mesh refinement to attain resolutions comparable to our models #17 and #18. Figure 14 plots the angular momentum profiles from two of their models with subthermal planet masses of ${q}_{\mathrm{thermal}}=0.2$ and 0.4. Their results agree with ours; the maximum angular momentum from their profiles differ from ours by about 10%, and their Keplerian disk sizes are also similar. We emphasize that the three codes we have used for comparison, PEnGUIn, Athena++, and Antares, all used different setups: PEnGUIn uses a nonuniform spherical grid centered on the star, Athena++ uses a logarithmic spherical grid centered on the planet, and Antares uses a cylindrical grid with mesh refinement centered on the star. Additionally, a similar isothermal simulation carried out by Tanigawa et al. (2012) also found a maximum specific angular momentum of about $0.2\,{r}_{{\rm{B}}}{c}_{\mathrm{iso}}$ (expressed in their units as $0.7\,{r}_{{\rm{H}}}^{2}{{\rm{\Omega }}}_{{\rm{p}}}$) for a ${q}_{\mathrm{thermal}}=3$ planet. The agreement between all these results lends confidence to our findings.

Zoom In Zoom Out Reset image size

Similar experiments have been performed by Ormel et al. (2015). They simulated planets with ${q}_{\mathrm{thermal}}=0.01$ and report that there is negligible rotation in the CPD. Their simulation domain extends as close to the planet as about $0.055\,{r}_{{\rm{B}}}$. Because this is similar to ${r}_{{\rm{K}}}$, the region where we expect to see Keplerian rotation is cut out from their domain. Their Figure 4, top panel, suggests that if they had set their inner boundary smaller, they would have seen faster rotation.

We also inspect how the rotation rate changes vertically. Figure 15 plots the azimuthally averaged r–z rotation profiles from models #1, 2, and 4. The speed is slower at higher altitudes, but the overall prograde rotation pattern does extend vertically to about $1\,{r}_{{\rm{B}}}$. At higher planet masses, rotation in the CPD appears to become more columnar.

Zoom In Zoom Out Reset image size

3.2.2. Adiabatic Envelopes

The right panel of Figure 12 tells a much different story for adiabatic CPDs. Unlike the isothermal cases, we do not find speeds close to the Keplerian value in any of our simulations. Even at ${q}_{\mathrm{thermal}}=4$, the rotation speed reaches only about one-third of the Keplerian speed.

Also unlike the isothermal CPDs, where rotation profiles follow a similar pattern regardless of ${q}_{\mathrm{thermal}}$, adiabatic CPDs increase in rotation speed as ${q}_{\mathrm{thermal}}$ increases. This trend can be understood as an effect of the Coriolis force. Because the adiabatic gas around the planet is, as we have seen in Figure 3, roughly in hydrostatic equilibrium, vertical gas flow toward the planet must be deflected and turned to planar motion as it flows over this ball of hydrostatic gas. In the planet's frame, the deflected gas must then be torqued by the Coriolis force into prograde rotation. Assuming the bound atmosphere has a size of ${{Cr}}_{{\rm{B}}}$, where C is a scaling coefficient, the rotation speed by the time the gas reaches the equator plane is approximately:

$$\begin{eqnarray}&&{v}_{\phi }\approx {v}_{\mathrm{gas}}{{\rm{\Omega }}}_{{\rm{p}}}\displaystyle \frac{{{Cr}}_{{\rm{B}}}}{{v}_{\mathrm{gas}}},\end{eqnarray} \tag{ 28 }$$

where ${v}_{\mathrm{gas}}$ is the meridional speed of the gas as it gets deflected around and flows over the planet's atmosphere, ${v}_{\mathrm{gas}}{{\rm{\Omega }}}_{{\rm{p}}}$ is the Coriolis force (dropping the factor of two), and ${{Cr}}_{{\rm{B}}}/{v}_{\mathrm{gas}}$ is the timescale of the deflection (also dropping order-unity prefactors). We can combine this expression with Equation (25) to scale it with the Keplerian speed. At the atmosphere's boundary, $r={{Cr}}_{{\rm{B}}}$, the rotation speed as a fraction of the Keplerian speed is then:

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\phi }}{{v}_{{\rm{K}}}}\approx {C}^{3/2}\displaystyle \frac{{r}_{{\rm{B}}}{{\rm{\Omega }}}_{{\rm{p}}}}{{c}_{\mathrm{iso}}}={C}^{3/2}{q}_{\mathrm{thermal}}.\end{eqnarray} \tag{ 29 }$$

This fraction scales linearly with ${q}_{\mathrm{thermal}}$, in rough agreement with our results in the right panel of Figure 12. We have seen that adiabatic CPDs are approximately bound within $0.2\,{r}_{{\rm{B}}}$. Plugging this into Equation (29), we get ${v}_{\phi }/{v}_{{\rm{K}}}\sim 0.1{q}_{\mathrm{thermal}}$, which is in good quantitative agreement with our results. Furthermore, this suggests that adiabatic CPDs can potentially become rotationally supported if ${q}_{\mathrm{thermal}}\gtrsim 10$.

Our adiabatic results can be roughly compared to simulations where effects of radiative transfer are included, such as those by Szulágyi et al. (2016), Szulágyi (2017), Cimerman et al. (2017), and Lambrechts & Lega (2017). Before the gas can cool significantly, it is roughly adiabatic and is comparable to our simulations. Our results are similar to those by Cimerman et al. (2017) and Lambrechts & Lega (2017), who also found little-to-no rotation in their circumplanetary gas. We agree qualitatively with Szulágyi et al. (2016) and Szulágyi (2017), in that they measured slower rotation when the gas cools more slowly, but we note that their simulations use large values of ${q}_{\mathrm{thermal}}$ ranging from 8 to 80, significantly different from our parameter space.

Isothermal models are commonly interpreted to represent the final state of the protoplanet's atmosphere, after it has fully cooled to the background nebular temperature (e.g., Lee & Chiang 2015, their Figure 4 and related discussion). In 1D, spherically symmetric models, such a final state would be described by the hydrostatic profile of Equation (20). In 3D hydrodynamical simulations, we have seen that the density profile differs (Figure 4), having much lower densities. It seems clear, then, that 1D models overestimate the gas mass in isothermal planetary atmospheres, and to this extent may overpredict the likelihood of giant planet formation. Here, we measure the total gas mass, ${M}_{\mathrm{CPD}}$, in our CPDs and compute the gas-to-core mass ratios, $\mu \equiv {M}_{\mathrm{CPD}}/{M}_{{\rm{p}}}$.

We compute ${M}_{\mathrm{CPD}}$ by summing the gas mass within a sphere of $0.1\,{r}_{{\rm{B}}}$ for the isothermal cases or $0.2\,{r}_{{\rm{B}}}$ for the adiabatic cases, following the CPD sizes measured in Section 3.1. Because we only simulate half the disk and assume midplane symmetry, we multiply the sum total by two to get the full ${M}_{\mathrm{CPD}}$. A note about units: the mass so computed is scaled to the ambient nebular gas density ${\rho }_{0}$, and has units of ${M}_{* }$. Because the code takes ${\rho }_{0}=1$ in units of ${M}_{* }/{R}_{{\rm{p}}}^{3}$ (the exact value is immaterial because gas self-gravity is neglected), to scale to any other nebular density ${\rho }_{0}$, we multiply by ${\rho }_{0}/({M}_{* }/{R}_{{\rm{p}}}^{3})$. To then convert into physical units, we multiply by ${M}_{* }$. In sum, to convert ${M}_{\mathrm{CPD}}\,(\mathrm{code}\,\mathrm{units})$ into ${M}_{\mathrm{CPD}}$ in physical units, we compute ${M}_{\mathrm{CPD}}={M}_{\mathrm{CPD}}\,(\mathrm{code}\,\mathrm{units})\times {\rho }_{0}{R}_{{\rm{p}}}^{3}/{M}_{* }$ $\times \,{M}_{* }\,={M}_{\mathrm{CPD}}\,(\mathrm{code}\,\mathrm{units})\times {\rho }_{0}{R}_{{\rm{p}}}^{3}$.

We find the gas-to-core mass ratios μ in both our isothermal and adiabatic simulations to scale the same way with ${q}_{\mathrm{thermal}}$. For subthermal planets, not surprisingly, we find ${M}_{\mathrm{CPD}}$ to scale with the volume of the Bondi sphere, ${r}_{{\rm{B}}}^{3}$, which implies $\mu \propto {r}_{{\rm{B}}}^{3}/{M}_{{\rm{p}}}\propto {M}_{{\rm{p}}}^{2}{H}_{{\rm{p}}}^{-6}\propto {q}_{\mathrm{thermal}}^{2}$. This scaling breaks down when we reach superthermal masses, where μ starts to scale linearly with ${q}_{\mathrm{thermal}}$ instead; this can be understood as ${M}_{\mathrm{CPD}}\propto {r}_{{\rm{B}}}^{2}{h}_{\mathrm{PPD}}$, where ${h}_{\mathrm{PPD}}={R}_{{\rm{p}}}{H}_{{\rm{p}}}$ is the protoplanetary (circumstellar) disk scale height. The left panel of Figure 16 shows these scalings match well with our measurements.

Zoom In Zoom Out Reset image size

For isothermal CPDs, $\mu ={\mu }_{\mathrm{iso}}$ is about three times higher in our high-resolution models compared to our fiducial ones. Because we trust the high-resolution results more—but also because the lower-resolution fiducial simulations better sample parameter space—we use the fiducial models to guide our scaling and the high-resolution models to normalize these scalings. Our final, empirical measurement of ${\mu }_{\mathrm{iso}}$ is:

$$\begin{eqnarray}{\mu }_{\mathrm{iso}}=\left\{\begin{array}{ll}24\,{q}_{\mathrm{thermal}}^{2}\times {\rho }_{0}{R}_{{\rm{p}}}^{3}/{M}_{* }, & {q}_{\mathrm{thermal}}\leqslant 0.5;\\ 12\,{q}_{\mathrm{thermal}}\times {\rho }_{0}{R}_{{\rm{p}}}^{3}/{M}_{* }, & {q}_{\mathrm{thermal}}\gt 0.5,\end{array}\right.\end{eqnarray} \tag{ 30 }$$

where the multiplicative factor allows us to scale to any desired background nebular density ${\rho }_{0}$ (see above note about units). The Athena++ simulation for model #17 produces a higher ${\mu }_{\mathrm{iso}}$ than PEnGUIn does by one order of magnitude. This discrepancy can also be seen in Figure 4. As mentioned in the beginning of Section 3, this is likely due to mass accumulating in front of the reflecting boundary used by Athena++. The effects of different boundary conditions at the planetary core need to be investigated further in the future. PEnGUIn has no boundary at the planet's location, so its results are easier to interpret.

Figure 17 further illustrates this difference between PEnGUIn and Athena++ and investigates the steadiness of our CPDs. It plots the mass flux, $\dot{M}$, across the sphere at a given distance r centered on the planet. PEnGUIn shows a close balance between the in (toward the planet) and out (away from the planet) fluxes, while the influx dominates in Athena++, resulting in a higher accretion rate. Despite the difference, the net $\dot{M}$ values are small in both cases. We get roughly ${M}_{\mathrm{CPD}}/\dot{M}\sim {10}^{3}\mbox{--}{10}^{4}\,{{\rm{\Omega }}}_{{\rm{p}}}^{-1}$, which we consider nearly steady. A similar level of steadiness is found in all of our models.

Zoom In Zoom Out Reset image size

We emphasize that our measurements of ${\mu }_{\mathrm{iso}}$ are many orders of magnitude below what they would be if the gas were to follow the 1D hydrostatic profile described by Equation (20). This holds true for both PEnGUIn and Athena++ results, and represents one of the most important differences between 1D models and 3D hydrodynamics simulations.

For adiabatic CPDs, resolution is less of a concern. Because adiabatic density profiles closely follow the 1D hydrostatic solution at fiducial resolution, we believe the numerical solution to have converged. A fit to our measured values of ${\mu }_{\mathrm{ad}}$ is given by:

$$\begin{eqnarray}{\mu }_{\mathrm{ad}}=\left\{\begin{array}{ll}0.5\,{q}_{\mathrm{thermal}}^{2}\times {\rho }_{0}{R}_{{\rm{p}}}^{3}/{M}_{* }, & {q}_{\mathrm{thermal}}\leqslant 1.6;\\ 0.8\,{q}_{\mathrm{thermal}}\times {\rho }_{0}{R}_{{\rm{p}}}^{3}/{M}_{* }, & {q}_{\mathrm{thermal}}\gt 1.6.\end{array}\right.\end{eqnarray} \tag{ 31 }$$

This result is more sensitive to the choice of CPD size, because adiabatic CPDs are much less centrally concentrated than isothermal ones. If we had used $1\,{r}_{{\rm{B}}}$, for example, instead of $0.2\,{r}_{{\rm{B}}}$, then ${\mu }_{\mathrm{ad}}$ would increase by about an order of magnitude. This is not concerning because we have seen evidence that the gas beyond $\sim 0.2\,{r}_{{\rm{B}}}$ is unbound (Section 3.1).

We can estimate numerical values for ${\mu }_{\mathrm{iso}}$ for real-world applications. For this, we use the same disk density profile as in Equation (14) and set ${\rho }_{0}={{\rm{\Sigma }}}_{\mathrm{MMSN}}/\sqrt{2\pi {h}_{\mathrm{PPD}}^{2}}$, where ${{\rm{\Sigma }}}_{\mathrm{MMSN}}=1700\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ is the surface density of the minimum-mass solar nebula (MMSN) at 1 au. For the temperature profile, we choose one such that ${h}_{\mathrm{PPD}}/R=0.035\,{(R/1\mathrm{au})}^{1/4}$. For ${M}_{* }$, we use one solar mass. Plugging these values into Equation (30), we show in the right panel of Figure 16 the gas-to-core mass ratios of planets with core masses ranging from 5 to 20 ${M}_{\oplus }$. Within a few au, the cores are superthermal and ${\mu }_{\mathrm{iso}}$ is on the order of a few percent. Outside a few au, the cores are subthermal and have gas-to-core mass ratios that decrease with distance to values less than a percent.

As mentioned at the beginning of this subsection, the final state of an atmosphere that cools and concomitantly accretes is an isothermal one. As such, isothermal simulations might be expected to yield maximum gas-to-core mass ratios. If we interpret ${\mu }_{\mathrm{iso}}$ as shown in Figure 16 along these lines, then we might conclude that gas mass fractions of super-Earth cores always remain much less than unity, even when such cores are embedded in a gas-rich disk like the MMSN. This finding agrees with inferred gas-to-core mass ratios of observed super-Earths (e.g., Wu 2019), but it would also imply that gas giants cannot form at distances of a few au (where most are actually found; see Nielsen et al. 2019) unless the background disk were at least an order of magnitude more massive than the MMSN. On the other hand, it also seems possible that ${\mu }_{\mathrm{iso}}$, though representing a maximally cooled state, does not necessarily equal the maximum μ possible. We return to this possibility in Section 4.1.

The flow structure in a CPD is close to axisymmetric within $\sim 0.2{r}_{{\rm{B}}}$, but becomes asymmetric as it merges with the background Keplerian flow. Figures 18 and 19 illustrate the flow patterns in the midplane for isothermal and adiabatic runs, respectively. For isothermal gas, the outflow (traced by magenta lines) is directed predominantly through the L1 and L2 Lagrange points when ${q}_{\mathrm{thermal}}\gt 1$ (Figure 18, right panel); by comparison, when ${q}_{\mathrm{thermal}}\lt 1$, these channels widen (Figure 18, left panel). This is expected because outflow speeds are generally subsonic, so if the gravitational potential at ${r}_{{\rm{H}}}$ is much larger than the internal energy of the gas (so-called "cold" flows), then matter can only exit the Hill sphere near the Lagrange points. If instead the internal energy dominates, then it becomes possible to overflow the Hill sphere in all directions. This criterion to open up the outflow channel can be written as:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{c}_{\mathrm{iso}}^{2}}\displaystyle \frac{{GM}}{{r}_{{\rm{H}}}}={3}^{\tfrac{1}{3}}{q}_{\mathrm{thermal}}^{\tfrac{2}{3}}\gt 1.\end{eqnarray} \tag{ 32 }$$

This translates to ${q}_{\mathrm{thermal}}\gt 0.6$, which we find to be consistent with our results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The story is similar with the adiabatic cases, although there are some differences. In both panels of Figure 19, we find the outflow to be focused toward the outer disk. The fact that this outflow connects to horseshoe orbits (Fung et al. 2015) implies the outward horseshoe turns (blue streamlines that turn radially outward) are wider than the inward turns (red streamlines that turn radially inward). This asymmetry, which is strongest for subthermal planets, has been seen in previous studies (e.g., Paardekooper & Mellema 2008; Masset & Casoli 2009; Jiménez & Masset 2017) and is related to the entropy gradient in the PPD. Our setup introduces a positive radial entropy gradient when the gas is adiabatic, which is indeed expected to widen the outward horseshoe turns.

Interestingly, as planet mass increases and becomes superthermal, isothermal and adiabatic results seem to converge, as seen in the right panels of Figures 18 and 19. This implies that when the planet is superthermal, CPD dynamics is dictated by gravity and the EOS is relegated to a more minor role.

We also look into vertical variations in the horseshoe orbits. Figure 20 plots the streamlines of the widest horseshoe orbits at different altitudes. Fung et al. (2015) and Masset & Benítez-Llambay (2016) showed that, for isothermal disks, horseshoe orbits should align into columns. We confirm that this remains true in our isothermal simulations, as shown in the left panel of Figure 20. Fung et al. (2015) suggested that it is an effect similar to Taylor–Proudman columns, and therefore might not apply to non-isothermal disks, where baroclinicity can alter the vorticity of the gas. The right panel of Figure 20 shows our results for an adiabatic case. The width of the horseshoe column gradually shrinks as altitude increases, and becomes about half its midplane value at three scale heights. We conclude that, while baroclinicity does introduce some variations, horseshoe orbits are still mostly columnar in adiabatic disks.

Zoom In Zoom Out Reset image size

Planets torque the gas in their orbits and open gaps in the PPD. In this section, we look into how gap opening affects the CPD. We have seen that the CPD mass—or equivalently, the gas-to-core mass ratio μ—is proportional to the background density, as described by Equations (30) and (31). From that, we can naively expect the CPD mass to decrease as the gap forms. To test this, we extend model #14, which has the highest planet mass and exerts the strongest planetary torque, to 100 orbits.

In the left panel of Figure 21, we compare how the mean surface densities of the CPD and of the PPD gap evolve with time. For the CPD, we compute the surface density by integrating the total mass within a cylinder of radius of $1\,{r}_{{\rm{H}}}$ around the planet and a vertical length equal to our simulation domain, and divide that by the surface area $2\pi {r}_{{\rm{H}}}^{2}$. For the PPD gap density, we do the same for the region between $R=\{{R}_{{\rm{p}}}-{r}_{{\rm{H}}},{R}_{{\rm{p}}}+{r}_{{\rm{H}}}\}$, ${\rm{\Psi }}=\{{{\rm{\Psi }}}_{{\rm{p}}}-0.5,{{\rm{\Psi }}}_{{\rm{p}}}+0.5\}$, and $\theta =\{\pi /2-0.3,\pi /2\}$, with the CPD region excised.

Zoom In Zoom Out Reset image size

The two surface densities evolve over time following a similar pattern. Because the sound crossing time in the CPD, $\sim {r}_{{\rm{H}}}/{c}_{\mathrm{iso}}$, is about 1 ${{\rm{\Omega }}}_{{\rm{p}}}^{-1}$, it is not surprising that the CPD reacts quickly to the emptying gap. We therefore conclude that Equations (30) and (31) can also be used for gap-opening planets, as long as ${\rho }_{0}$ accounts for gap depletion.

Observationally, we know that PPDs dissipate over a few million years. Taken at face value, our results imply that as PPDs dissipate, CPDs should dissipate with them. In reality, however, we do not expect CPD evolution to play out so simply—in part because we have neglected cooling of the bound gas, which enables them to contract and survive the loss of external pressure from the dissipating nebula. In 1D cooling models, the evolution of the atmosphere is controlled by its radiative–convective boundary, whose properties are insensitive to the nebular density at large (e.g., Lee & Chiang 2015). Massive planetary envelopes can accrete and survive even in nearly gas-free disks, at least in 1D (Lee et al. 2018). We will return to this tension between 1D cooling models and 3D hydrodynamic models in Section 4.1.

We note that our exploration of how PPD gaps influence CPDs is also limited because our simulations are optimized for smaller-scale CPDs and not for larger-scale phenomena. Our spatial resolution is poor far from the planet, with attendant problems in numerical diffusion. Moreover, 100 orbits are far from sufficient to evolve the gap to a steady state.

We found that the gas-to-core mass ratio μ remains below 10% even for a 20 ${M}_{\oplus }$ core surrounded by fully cooled, isothermal gas (right panel of Figure 16). Theoretically, the adiabatic and isothermal cases should bracket a planet's thermal (i.e., accretion) history—the atmosphere/CPD starts off behaving adiabatically on timescales shorter than the cooling time, and becomes isothermal on timescales longer than the cooling time (e.g., Lee et al. 2014; Lee & Chiang 2015; Ginzburg et al. 2016; Coleman et al. 2017). It is expected that μ will evolve from the adiabatic to the isothermal state monotonically.⁶ Thus, given our result that ${\mu }_{\mathrm{iso}}\lt 10 \%$, it would seem unlikely that envelope self-gravity would ever become significant enough to trigger "runaway accretion" and gas giant formation (e.g., Pollack et al. 1996; Ikoma et al. 2000).

Are there ways out of this conclusion? Is it possible for ${\mu }_{\mathrm{iso}}$ to be larger than 10%? In 3D, we have seen that isothermal CPDs are rotationally supported within $\sim 0.05{r}_{{\rm{B}}}$. Rotationally supported envelopes can have arbitrary masses and density profiles (within the bounds of gravitational and hydrodynamic, e.g., Rayleigh stability). Indeed, the densities given by PEnGUIn and Athena++ do not agree inside $0.05{r}_{{\rm{B}}}$ (Figure 4). This leaves much room for speculation. An isothermal CPD could potentially become more massive if there are angular momentum transport mechanisms that allow it to accrete. Zhu et al. (2016), for example, reported shock-driven and vortex-driven accretion in their 2D simulations. Reality might be a mixture of the 1D and 3D models. The outer parts of the atmosphere (still within $0.05{r}_{{\rm{B}}}$) may become radiative, nearly isothermal, and disky, with complex 3D flow structures like what we have seen in this work, while the inner parts may be convective, nearly adiabatic, and spherically symmetric. Gas might accrete across the isothermal disk and pile on top of the adiabatic envelope, in a fashion similar to the way circumstellar disks feed protostars.

We illustrate these ideas in Figure 22 by drawing some schematic evolutionary paths. One-dimensional models predict a cooling phase followed by a runaway phase after μ reaches unity (orange path). Our 3D simulations, taken at face value, indicate that cooling alone leads to much smaller values of μ (cyan path). However, if one combines cooling, 3D hydrodynamics, and disk accretion physics, then one might produce an evolutionary path resembling the green path. Gas giants may form if disk accretion (and eventually, self-gravity) push μ above unity.

Zoom In Zoom Out Reset image size

The presence of prograde, low-inclination, low-eccentricity "regular" satellites around the giant planets in our solar system suggests that there once existed rotationally supported CPDs around them, much like the ones we discover in our isothermal simulations. Canup & Ward (2006) found that the total mass of each satellite system is lower than its host's mass by a factor of ∼10⁻⁴, which implies, if one assumes a gas-to-solid mass ratio of 100, a gaseous-CPD-to-planet mass ratio of ∼10⁻². This is encouraging, as it is of the same order as our measured values for ${\mu }_{\mathrm{iso}}$.

If satellite systems are formed in CPDs like the ones we simulated, we expect the former to have sizes comparable to or smaller than (if inward migration of solids is significant) $0.05\,{r}_{{\rm{B}}}$ (for ${q}_{\mathrm{thermal}}\lt 1$). We test this expectation here. To evaluate ${r}_{{\rm{B}}}$, one needs to estimate ${c}_{\mathrm{iso}}$ where the planets formed. For simplicity, we will assume they formed near their current positions, and use the same temperature profile as the one used in Section 3.3, where the temperature is about 300 K at R = 1 au and scales as ${R}^{-1/2}$. We approximate the sizes of the regular satellite systems using the semimajor axes of their outermost members.

4.2.1. Jupiter

4.2.2. Saturn

4.2.3. Uranus

4.2.4. Neptune