GAP OPENING IN 3D: SINGLE-PLANET GAPS

The initial disk profile assumes the following surface density and sound speed profiles:

$$\begin{eqnarray}&&{\rm{\Sigma }}={{\rm{\Sigma }}}_{0}{\left(\displaystyle \frac{R}{{R}_{{\rm{p}}}}\right)}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&c={c}_{0}{\left(\displaystyle \frac{R}{{R}_{{\rm{p}}}}\right)}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 5 }$$

where we set ${{\rm{\Sigma }}}_{0}=1$,³ and we choose ${c}_{0}=0.05$. The disk scale height h is therefore $c/{{\rm{\Omega }}}_{{\rm{k}}}=0.05R$ with ${{\rm{\Omega }}}_{{\rm{k}}}$ evaluated at the midplane; in other words, the disk aspect ratio $h/R=0.05$ is constant. This sound speed profile is fixed in time, corresponding to an irradiated disk that has a cooling time much shorter than the dynamical time. In 3D hydrostatic equilibrium, the initial density structure is

$$\begin{eqnarray}&&\rho ={\rho }_{0}{\left(\displaystyle \frac{R}{{R}_{{\rm{p}}}}\right)}^{-\tfrac{3}{2}}\exp \left(\displaystyle \frac{{GM}}{{c}^{2}}\left[\displaystyle \frac{1}{r}-\displaystyle \frac{1}{R}\right]\right),\end{eqnarray} \tag{ 6 }$$

where ${\rho }_{0}={{\rm{\Sigma }}}_{0}/\sqrt{2\pi {h}_{0}^{2}}$ is the initial gas density at the location of the planet.

The initial velocity field assumes hydrostatic equilibrium, taking into account gas pressure but not viscosity. As a result, the initial disk has an orbital frequency of

$$\begin{eqnarray}&&{\rm{\Omega }}=\sqrt{{{\rm{\Omega }}}_{{\rm{k}}}^{2}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial p}{\partial r}}\end{eqnarray} \tag{ 7 }$$

and zero radial and polar velocity. Given our isothermal equation of state and sound speed profile, our disk model is subject to the vertical shear instability (e.g., Urpin & Brandenburg 1998; Lin & Youdin 2015). However, our disk viscosity of $\nu =2.5\times {10}^{-6}$ is sufficient to prevent the instability from growing, as shown by Nelson et al. (2013).

Both the inner and outer radial boundaries are fixed at their initial values. In the polar direction, we simulate only the upper half of the disk to reduce computational cost. We use reflecting boundaries both at the top, to prevent mass from entering or leaving the domain, and at the midplane. Imposing reflection symmetry across the midplane (polar velocity ${v}_{\theta }=0$) means that certain dynamics are not captured by our simulation, such as flows that cross the midplane, and bending waves near the midplane. It remains to be determined whether these restrictions affect the results of our study.

Our simulation domain extends radially from 0.4R_p to 2${R}_{{\rm{p}}}$ and spans the full $2\pi$ in azimuth. The polar angle in our 3D simulations spans 0°–12°, equivalent to a vertical extent of 4 scale heights.

We use 270 logarithmically spaced cells in the radial direction, 810 uniformly spaced cells in the azimuthal direction, and 45 uniformly spaced cells in the polar direction. This resolution translates to about 8, 6, and 11 cells per scale height in each of the $\{r,\phi ,\theta \}$ directions, respectively. This choice is motivated by Figure 2 of FSC14.

As in FSC14, we define the surface density in the gap, ${{\rm{\Sigma }}}_{\mathrm{gap}}$, as a spatial average over the annulus spanning $R={R}_{{\rm{p}}}-{\rm{\Delta }}$ to ${R}_{{\rm{p}}}+{\rm{\Delta }}$ with ${\rm{\Delta }}\equiv 2\,{r}_{{\rm{H}}}$, excised from $\phi ={\phi }_{{\rm{p}}}-{\rm{\Delta }}/{R}_{{\rm{p}}}$ to ${\phi }_{{\rm{p}}}+{\rm{\Delta }}/{R}_{{\rm{p}}}$. Throughout this paper, we will also show velocity plots and maps that are temporal and/or spatial averages. In all instances, velocities are averages weighted by gas density.

We perform four 3D simulations with q varying between 0.0005 and 0.004 (0.5–4 ${M}_{{\rm{J}}}$). For our parameters, the planet's Hill radius ${r}_{{\rm{H}}}$ ranges from 1.1 to 2.2 times the local disk scale height h₀. We also run four 2D simulations using the same parameters for comparison; these are nearly exact replicates of the simulations from FSC14. Figure 1 plots the gap depth contrast ${{\rm{\Sigma }}}_{\mathrm{gap}}/{{\rm{\Sigma }}}_{0}$ versus time for our 2D and 3D simulations. Because of the costly nature of 3D simulations, we were only able to simulate up to 10³ orbits. By comparison, FSC14 ran up to ∼10⁴ orbits to achieve steady state. Comparing our Figure 1 with Figure 5 of FSC14 indicates that our results for ${{\rm{\Sigma }}}_{\mathrm{gap}}/{{\rm{\Sigma }}}_{0}$ have converged to within a factor of 2 of theirs. We consider such convergence sufficient to identify differences between 2D and 3D gaps. Our results show some time variability, particularly for high-mass runs, that we will investigate in Section 3.3.

Zoom In Zoom Out Reset image size

Figure 1 shows clearly that ${{\rm{\Sigma }}}_{\mathrm{gap}}$ behaves in 3D similarly to how it does at lower dimensionality. In particular, the final values of ${{\rm{\Sigma }}}_{\mathrm{gap}}$ in 3D scale as ${q}^{-2}$, as predicted by the 0D models of FSC14 and Kanagawa et al. (2015). The agreement between 2D and 3D simulations is astonishingly close for ${M}_{{\rm{p}}}\leqslant {M}_{{\rm{J}}}$. For ${M}_{{\rm{p}}}=2{M}_{{\rm{J}}}$ and $4{M}_{{\rm{J}}}$, ${{\rm{\Sigma }}}_{\mathrm{gap}}$ decreases initially more rapidly in 3D than in 2D, but ultimately reaches a slightly larger value (about ∼80% larger in the case ${M}_{{\rm{p}}}=2{M}_{{\rm{J}}}$).

For 3D simulations, to verify that our fiducial choice for smoothing length ${r}_{{\rm{s}}}=0.1{r}_{{\rm{H}}}$ is sufficiently small, we perform a test case at ${M}_{{\rm{p}}}=4{M}_{{\rm{J}}}$ with ${r}_{{\rm{s}}}=0.05{r}_{{\rm{H}}}$ and find that ${{\rm{\Sigma }}}_{\mathrm{gap}}$ changes negligibly. For 2D simulations, the correct choice of ${r}_{{\rm{s}}}$ is less clear; although any value of order h₀ is suitable in principle (Müller et al. 2012), the value of ${{\rm{\Sigma }}}_{\mathrm{gap}}$ is sensitive to ${r}_{{\rm{s}}}$. For example, we find that using ${r}_{{\rm{s}}}=1\,{h}_{0}$ instead of our fiducial $0.5\,{h}_{0}$ results in a ∼40% larger ${{\rm{\Sigma }}}_{\mathrm{gap}}$. In this regard, we should not overinterpret the exceptionally close agreement between 2D and 3D results.

Figure 2 plots the azimuthally averaged surface density profiles, excluding a small $\pm {\rm{\Delta }}/{R}_{{\rm{p}}}$ azimuthal range around the planet. In terms of both gap width and the sharpness of gap edges, the 3D gaps are remarkably similar to 2D gaps. The resemblance is confirmed in face-on views of surface density, as presented in Figure 3. The 2D gaps show more distinct streamers within the gaps. These filaments appear to be related to the Rayleigh instability operating at gap edges, as we discuss in Section 3.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The flow pattern in 3D may be expected to deviate significantly from that in 2D. In 3D, not only are vertical velocities allowed, but radial velocities can be larger since planetary and viscous torques do not have to cancel locally; the former is mostly exerted in the midplane, while the latter can be exerted anywhere.

Figures 4 and 5 illustrate the meridional flow in our $1{M}_{{\rm{J}}}$ and $4{M}_{{\rm{J}}}$ simulations. In agreement with the results of Szulágyi et al. (2014) and Morbidelli et al. (2014), we find that within the gap, there exists a radial flow directed away from the planet's orbit near the midplane, as well as another flow delivering gas back into the gap at higher altitudes (Figure 4). Whereas Morbidelli et al. (2014) describe a closed, localized circulation, where the gas is trapped in a cyclical flow near the gap, the flow pattern that we uncover is less localized. We find vigorous motion outside the gap as well as inside with multiple eddy-like structures (Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 is a schematic drawing of the overall meridional flow field. Outside the gap, steep density gradients drive a viscous flow into the gap ("gradient-driven" flow shown in red). At the same time, the planet drives a fast, denser midplane flow directed away from the gap ("planet-driven" flow shown in blue). These two flows collide near the gap edges, where more complex flow patterns emerge. Because midplane disk pressure decreases outward, the planet-driven flow encounters little resistance moving beyond the outer gap edge and is able to reach a radial distance of ∼4${r}_{{\rm{H}}}$ away from the planet's orbit before being deflected. By comparison, on the opposite, interior half of the gap, the planet-driven flow is deflected at ∼2${r}_{{\rm{H}}}$ because of the higher gas pressure characterizing the inner gap edge. Where the planet-driven flow collides with the gradient-driven flow, the local pressure and density are enhanced. The buildup of pressure at the midplane drives gas upward (the "combined" flow shown in magenta), which eventually falls back down toward the midplane. Meridional eddies are created whose rotation poles point in the $-\phi$ direction, just outside either gap edge. The eddies are strongest near the planet and persist downstream for ∼1–2 rad.

Zoom In Zoom Out Reset image size

The meridional eddies extend as high as our vertical simulation boundaries allow. Gas is carried all the way from the disk midplane to a few scale heights above and back. This vertical circulation enables the system to reach a near-steady state: the angular momentum acquired by a gas element at low altitude, where the planetary torque dominates, nearly cancels with the angular momentum acquired at high altitude, where the viscous torque dominates. This change in the sign of the net torque is illustrated by the left panel of Figure 7, which shows that the planetary torque, integrated over a portion of the outer disk for the $1{M}_{{\rm{J}}}$ case, dominates near the midplane, but falls rapidly with altitude and is eventually overtaken by the viscous torque. Note that the precise altitude where this sign change occurs is radially dependent, because the planetary torque is strongest close to the planet, while the viscous torque is strongest near the gap edge. The right panel of Figure 7 plots the angular momentum flux in the polar direction, ${F}_{\theta }$, integrated over the same portion of the outer disk:

$$\begin{eqnarray}&&{F}_{\theta }={\int }_{0}^{2\pi }{\int }_{1.1}^{1.8}\langle \rho {v}_{\phi }{v}_{\theta }r\sin \theta \rangle r\,{dr}\,d\phi ,\end{eqnarray} \tag{ 8 }$$

where the brackets indicate a time average from 1000 to 1010 orbits. In Figure 7 we separate ${F}_{\theta }$ into its upward (${v}_{\theta }\lt 0$) and downward (${v}_{\theta }\gt 0$) components. The difference between the upward and downward fluxes is small but significant: the net residual is similar in magnitude to the planetary and viscous torques plotted in the left panel of Figure 7. This result is consistent with our interpretation that torque balance is achieved through vertical transport. A more explicit verification of the balance is difficult because the flow varies strongly in time and space, and flow patterns are not strictly closed.

Zoom In Zoom Out Reset image size

Figure 8 plots the vertically and azimuthally averaged meridional speeds of our 3D simulations. The speeds come within an order of magnitude of the sound speed over distances comparable to the planet's orbital radius and are supersonic closest to the planet (e.g., there is a Mach ∼5 shock deep inside the $4{M}_{{\rm{J}}}$ planet's Hill sphere). The planet's influence extends well beyond Hill sphere scales; for example, the most distant principal Lindblad resonances are located at $0.63\,{R}_{{\rm{p}}}$ in the inner disk (2:1 resonance) and $1.59\,{R}_{{\rm{p}}}$ in the outer disk (1:2 resonance).

Zoom In Zoom Out Reset image size

Figure 9 plots 2D maps of the vertical velocity ${v}_{{\rm{z}}}$, vertically averaged over the upper half of the disk above the midplane (recall that our simulations assume reflection symmetry about the midplane), and gives a sense of the flow's azimuthal variations. The fastest downward flows in the gap are concentrated within ∼1${r}_{{\rm{H}}}$ of the planet. The velocities are comparable to the sound speed, as befits gas that has freely fallen a distance $\sim {r}_{{\rm{H}}}$ over a timescale ${{\rm{\Omega }}}^{-1}$. Farther away azimuthally, the flow speed reduces by factors of a few and fluctuates spatially and temporally.

Zoom In Zoom Out Reset image size

Vertical velocities might be used one day to distinguish planet-opened gaps from other kinds of gaps (e.g., Zhang et al. 2015; Béthune et al. 2016), provided that velocity signatures above and below the disk midplane along the observer's line of sight do not cancel to zero. A net velocity of zero would be obtained if the disk were exactly reflection symmetric about the midplane, optically thin in the line used to measure Doppler shifts, and viewed face-on. Rotational transitions from CO may be optically thick, even in the low-density environment of a gap. For example, at a temperature of 100 K and a CO/H₂ number abundance of 10⁻⁴, the J = 4–3 CO line is optically thick at line center down to Σ ∼ 10⁻⁴ g/cm².

Here we offer an explanation for the origin of gap streamers, sometimes called filaments, that have been observed in gap-opening simulations (e.g., de Val-Borro et al. 2006). These streamers arise from unsteady gap edges. In our simulations, they are seen in both 2D and 3D, but more so in 2D, and when ${M}_{{\rm{p}}}\geqslant 2{M}_{{\rm{J}}}$.

The Rayleigh instability can operate at gap edges where the rotation profile is significantly modified by gas pressure gradients. This dynamical instability occurs wherever the specific angular momentum decreases outward:

$$\begin{eqnarray}&&\displaystyle \frac{\partial l}{\partial R}\lt 0,\end{eqnarray} \tag{ 9 }$$

where $l={R}^{2}{\rm{\Omega }}$ is the specific angular momentum. Combining this with Equation (7) imposes a condition on the second radial derivative of the gas pressure, and by extension the sharpness of gap edges (e.g., Yang & Menou 2010). Kanagawa et al. (2015) have used this condition to model gap profiles.

Figure 10 plots the azimuthally averaged values of $\partial l/\partial R$ for our 2D and 3D simulations. For 3D, we plot the midplane values. We find that when ${M}_{{\rm{p}}}\geqslant 2{M}_{{\rm{J}}}$ in 2D, or ${M}_{{\rm{p}}}=4{M}_{{\rm{J}}}$ in 3D, the gap edges are Rayleigh unstable. These cases coincide with the appearance of unsteady gap edges and enhanced streamers. We have continued our 2D, $4{M}_{{\rm{J}}}$ simulation to 10⁴ orbits to verify that the gap edges remain unstable even at late times.

Zoom In Zoom Out Reset image size

There is a slight difference between the 3D and 2D angular momentum profiles, in that the former violates the Rayleigh criterion to a lesser degree. Coinciding with this difference—and arguably causally related—is the fact that the streamers leaking into the gap are less prominent in 3D than in 2D (Figure 3). In general, as demonstrated in Figure 11, the 3D torques are somewhat weaker than the 2D torques—by up to ∼50% in the case of the $4{M}_{{\rm{J}}}$ run. We note that this implies that a larger smoothing length ${r}_{{\rm{s}}}$ in 2D may produce a better match to 3D results for higher planet masses.

Zoom In Zoom Out Reset image size

The Rayleigh instability may limit the extent to which gaps can be emptied. As the explicit disk viscosity decreases, the gap deepens, its edges become more prone to instability, and more streamers leak into the gap (see also Section 4.2 of FSC14). Thus, the Rayleigh instability serves to provide an effective viscosity, transporting angular momentum outward to counter the planet's torque and establishing a "floor" on the gap surface density.

Disk flows are directed onto the planet above its poles and away from the planet in the equator plane (Kley et al. 2001; Klahr & Kley 2006; Tanigawa et al. 2012; Morbidelli et al. 2014; Szulágyi et al. 2014; Fung et al. 2015; Ormel et al. 2015). The giant planets simulated here should be in the final "runaway" phase of their mass accretion history, when the self-gravity of their gas envelopes is significant. Accretion rates for planets in this regime, such as those simulated by D'Angelo et al. (2003) and Machida et al. (2010), and calculated analytically by Tanigawa & Tanaka (2016), show a clear dependence on the ambient gas density. The gap depth should therefore matter for planet accretion rates and luminosities. The converse should also be true. To assess the effect of planetary accretion on the gap density, we would need to adopt a different planet boundary condition than our present smoothing length prescription.

The large-scale vertical and radial flows induced by giant planets may affect a number of disk processes: the radial transport of grains (see, e.g., Ogliore et al. 2009, for possible implications of the Stardust mission); dust filtration at gap edges, which can impact the appearance of transitional disks (e.g., Zhu et al. 2014); the ability of dust to settle vertically, which can affect the disk's spectral energy distribution (e.g., Chiang et al. 2001); and disk chemistry (e.g., Bergin et al. 2007, p. 751; Dutrey et al. 2014, p. 317).

The wide radial extents of the cavities in transitional disks have led some to invoke gap opening by multiple planets (Zhu et al. 2011). Simulations in 2D of multiplanet gaps have found that common gaps (the merged gaps of more than one planet) have significantly higher gas density than single-planet gaps (Duffell & Dong 2015), making it more difficult to reproduce the transparency of transition disk holes. We plan to test this result in 3D in a forthcoming paper.