HOW EMPTY ARE DISK GAPS OPENED BY GIANT PLANETS?

Although our study is primarily numerical, we derive here an approximate analytic relation for Σ_gap that we will use to put our numerical results in context. We admit at the outset that our derivation can hardly be called such, as our reasoning will ignore many details and make assumptions not carefully justified. But as the rest of our paper will show, the simple relation we now present will yield results surprisingly close to those of detailed numerical simulations.

First examine the outer disk, exterior to the planet's orbit. The outer Lindblad torque exerted by the planet transmits angular momentum outward at a rate

$$\begin{equation} T_{\rm L} \sim q^2 \left(\frac{r}{h} \right)^3 \Sigma _{\rm gap} \,\Omega ^2 r^4, \end{equation} \tag{ 1 }$$

where Ω is the disk angular frequency (Goldreich & Tremaine 1980). This is the total torque from linear perturbations to one side of the planet, integrated over all resonances up to the torque cutoff at azimuthal wavenumber m ∼ r/h (see, e.g., Equation (2) of Crida et al. 2006, and references therein). Note that we have used Σ_gap, the surface density averaged over the bottom of the gap, in our evaluation of the integrated Lindblad torque. Our justification for this choice is that the integrated torque is dominated by resonances at the torque cutoff, i.e., at distances ∼h from the planet (see Figures 2 and 3 of Goldreich & Tremaine 1980), and bottoms of gaps opened by giant planets in gas disks typically extend this far (see, e.g., the simulated gaps of Crida et al. 2006, or our Figure 2). One caveat is that the linear theory from which Equation (1) originates is formally valid only for low-mass planets for which q ≲ (h/r)³; the highest-mass planets we simulate will violate this condition, and indeed for such super-Jovian objects our numerical simulations will reveal deviations from the simple-minded scaling law we derive in this section.

In steady state, the outward transmission of angular momentum by the (outer) Lindblad torque must be balanced by the angular momentum transmitted inward by the viscous torque:⁵

$$\begin{equation} T_{\rm v} \sim \Sigma _0 \,\nu \,\Omega \,r^2, \end{equation} \tag{ 2 }$$

where ν = α(h/r)²Ωr² is the kinematic viscosity (e.g., Frank et al. 2002). Here we have used Σ₀, the surface density at the gap periphery—or more conveniently, the surface density at the planet's location if the planet were massless—to evaluate T_v. The viscous torque depends on the gradient of Σ, and this gradient is larger outside the flat-bottomed gap than inside it (see, e.g., the gap profile shown in Figure 2; by "gap periphery" we mean a location like r ≈ 1.4, where the gradient might reasonably be approximated as Σ₀/r, which is what Equation (2) essentially assumes). Conscripting Σ₀ in this way is a gross simplification, but alternatives would require that we actually compute the precise shape of the gap, which is what we are trying to avoid with our order-of-magnitude derivation.

Usually one thinks of viscous torques as transmitting angular momentum outward, but in the gap edge of the outer disk, the direction of viscous transport is inward because the surface density there has a sharp and positive gradient (dΣ/dr > 0).

Setting T_L = T_v yields

$$\begin{equation} \frac{\Sigma _{\rm gap}}{\Sigma _0} \sim \frac{\alpha (h/r)^5}{q^2} \,. \end{equation} \tag{ 3 }$$

Exactly the same scaling relation applies to the inner disk, interior to the planet's orbit. The signs in the inner disk are reversed from those in the outer disk: the inner Lindblad torque transmits angular momentum inward, while the local viscous torque transmits angular momentum outward.

As we were completing our numerical tests of Equation (3) and preparing our manuscript for publication, we became aware of the study by Duffell & MacFadyen (2013), who found the same scaling relation on purely empirical grounds (although these authors did not explicitly vary the Mach number r/h). Duffell & MacFadyen (2013) concentrated on the low-mass q ≲ 10⁻⁴ regime. Our study complements theirs by studying the high-mass q ≳ 10⁻⁴ regime; we will see to what extent Equation (3) also holds true for giant planets. See also our Section 4.3, where we discuss to what extent the short derivation given in this subsection captures the whole story.

Section 2 contains our numerical methods and simulation parameters. Section 3 presents our results for Σ_gap(q, α, h/r). Section 4 concludes and charts directions for future work.

PEnGUIn (Piecewise Parabolic Hydro-code Enhanced with Graphics Processing Unit Implementation) is a Lagrangian, dimensionally split, shock-capturing hydrodynamics code. We defer a detailed description to a future paper and here just mention a few salient points. The code is written in CUDA-c and runs on multiple GPUs for accelerated performance. It uses the piecewise parabolic method (PPM; Colella & Woodward 1984), and its main solver is modeled after VH-1 (Blondin & Lufkin 1993) with two main differences: (1) PEnGUIn explicitly conserves angular momentum by replacing, as one of the quantities a fluid element carries in the Lagrangian frame, the angular speed with the specific angular momentum, and (2) PEnGUIn uses a non-iterative Riemann solver for isothermal flows; see Balsara (1994), and note that their Lagrangian re-map formulation applies not only to strictly isothermal flows, but also to locally isothermal flows under the two-shock approximation by allowing for different sound speeds for left- and right-moving waves. PEnGUIn also contains a module that computes the divergence of the stress tensor through piecewise parabolic interpolations.

Technically, PEnGUIn's reference frame is a barycentric frame that rotates at Ω_p; thus, the planet's position is fixed in time, but the Coriolis force is not computed as an explicit source term. Rather, it is absorbed into the conservative form of the angular momentum equation (Kley 1998).

GPU acceleration is optimal for programs that are massively parallel; PEnGUIn achieves parallelization through memory management and domain splitting. The hardware used for this paper are three GTX-Titan graphics cards all connected to a single node. Running in double precision on all three cards simultaneously, PEnGUIn takes 12 s to run per planetary orbit for (q, α, h/r) = (10⁻³, 10⁻², 0.05).

For comparison with PEnGUIn, we also carried out simulations with ZEUS90: a modern version of ZEUS (Stone & Norman 1992; Hawley & Stone 1995) written in FORTRAN 90. It is a three-dimensional (3D), operator-split, time-explicit, Eulerian finite-differencing magnetohydrodynamics code, widely used to simulate a variety of systems, including magnetorotationally unstable circumbinary disks (Shi et al. 2012) and warped disks (Sorathia et al. 2013). For our application, we suppress the vertical dimension and magnetic fields; implement the Navier–Stokes viscosity module; and add a planetary potential having a specified time dependence. The von Neumann–Richtmyer artificial viscosity, commonly used to capture shock waves, is switched off in the presence of an explicit viscosity. The reference frame for ZEUS90 is a non-rotating frame centered on the star, and so ordinarily there is an extra term in Φ due to the indirect potential: $GM_{\rm p} r \cos (\phi -\phi _{\rm p})/r_{\rm p}^{2}$. However, we find in practice that the indirect term results in a "wobbling" of the disk that is difficult to reconcile with our fixed boundary conditions on fixed circles (see Equations (9)–(11) below and surrounding discussion). The wobbling generates spurious time variability that increases with increasing q; therefore, we drop the indirect potential in all ZEUS90 simulations with q ⩾ 1 × 10⁻³ (apparently Lubow et al. 1999 also dropped the indirect potential in their simulations with ZEUS; see also Zhu et al. 2011, who found that their planet-disk simulations were not sensitive to the indirect term). Running in double precision on 128 cores in parallel, ZEUS90 takes 26 s to run per planetary orbit for (q, α, h/r) = (10⁻³, 10⁻², 0.05).

For our parameter study we vary:

• 1.  
  q from 10⁻⁴ to 10⁻²,
• 2.  
  α from 10⁻³ to 10⁻¹, and
• 3.  
  h/r from 0.04 to 0.1.

Other properties of our simulations are as follows.

2.3.1. Initial and Boundary Conditions

Our simulation domain spans 0 to 2π in azimuth, and from r_in = 0.4 to r_out = 2.5 in radius (in units where the planet-star separation r_p = 1). Initial conditions correspond to a steady-state accretion disk having constant α, constant h/r, and a rotation curve modified by the radial pressure gradient:

$$\begin{eqnarray} \Sigma &= \Sigma _0 \, (r/r_{\rm p})^{-1/2} \,, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} v_{r} &= -\frac{3}{2}\, \alpha \, \left(\frac{h}{r}\right)^2\sqrt{\frac{G(M_\ast +M_{\rm p})}{r}}\,, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \Omega &= \sqrt{1-\frac{1}{2}\left(\frac{h}{r}\right)^2}\sqrt{\frac{G(M_\ast +M_{\rm p})}{r^3}}\,, \end{eqnarray} \tag{ 11 }$$

with Σ₀ = 1 (we could have chosen any value for Σ₀ because we do not calculate the gravity of the disk—the disk exerts no gravitational backreaction on the planet, nor does it self-gravitate). At both inner and outer radial boundaries, we fix Σ, v_r, and Ω to their values determined by the above equations. These fixed boundary conditions ensure a steady inflow of mass across the simulation domain—as is appropriate for real accretion disks.⁶ Figure 1 illustrates how our boundary conditions enable a planet-less disk to relax over a viscous diffusion timescale to the equilibrium profile described by Equations (9)–(11).

Zoom In Zoom Out Reset image size

Ideally the radial boundaries should be placed far enough away that waves launched from the planet damp before they reach the edges of the domain. Goodman & Rafikov (2001) calculated that nonlinear steepening of waves causes them to dissipate over length scales of ∼3h (for q = 10⁻³ and h/r = 0.05; the damping length scales as q^−0.4 and (h/r)^2.2). Our outer radial boundary of r_out = 2.5 (∼15–40 h away from the planet, depending on our choice for h/r) is distant enough that outward-propagating waves largely dissipate within the domain. For our inner radial boundary of r_in = 0.4 (∼6–15 h away from the planet), the situation is more marginal; depending on the simulation, waves are still present at our disk inner edge. However, the main focus of our paper is the surface density deep within the gap, and we have verified that this quantity changes by no more than ∼10% as we shrink r_in from 0.4 to 0.2. Thus, we opt for the larger boundary to keep code timesteps longer. For simplicity we eschew wave-killing zones (cf. de Val-Borro et al. 2006). In practice, the Godunov-type scheme used by PEnGUIn is effective at absorbing waves at fixed boundaries, even more so than using wave-killing zones (Z. Zhu 2013, private communication).

To avoid strong shocks at the beginning of the simulation, the planet mass is ramped from zero to its assigned value over an initial "warm-up" phase. In PEnGUIn, M_p increases according to M_p(t)/M_* = qsin ²[(Ω_pt/20)(10⁻³/q)]. For q = 10⁻³, this takes five orbits. In ZEUS90, the planet mass grows linearly from zero to its desired value in one orbit. Both warm-up schemes proved stable.

2.3.2. Grid Resolution

Our grid spacings are logarithmic in radius and uniform in azimuth. For h/r = 0.05, the resolution is 270 (r) × 810 (ϕ) for PEnGUIn and 256 × 864 for ZEUS90 (the latter choice yields square grid cells). Figure 2 attests that gap surface densities have largely converged at our standard resolution. We scale our grid cell size with h/r, i.e., with sound speed c_s, so that sound waves of a given frequency are equally well resolved between simulations. Cold disks with small h/r are especially costly, which is why we do not vary h/r below 0.04. Code timesteps scale with grid cell sizes according to the Courant–Friedrichs–Lewy condition, with the Courant number chosen to be 0.5 for PEnGUIn and 0.4 for ZEUS90.

Zoom In Zoom Out Reset image size

2.3.3. Smoothing Length r_s

2.3.4. Σ_gap and Convergence with Time

Our metric for gap depth is the space- and time-averaged surface density Σ_gap in the planet's co-rotation region, normalized to Σ₀ = 1 (the surface density at r = r_p = 1 in the absence of the planet). As judged from snapshots like those shown in Figures 3 and 4, the annulus spanning r = r_p − Δ to r_p + Δ with Δ ≡ 2max (R_H, h), excised from ϕ = ϕ_p − Δ/r_p to ϕ_p + Δ/r_p, is visibly depleted and reasonably uniform. For most simulations, this is the area over which we average Σ to calculate Σ_gap.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In a few cases the outer edge of the gap is visibly eccentric (see Section 3.2 for more discussion), and the circular annulus we have defined above becomes contaminated with non-gap material and is no longer suitable for measuring Σ_gap. In these cases, we keep the circular inner gap edge and the azimuthal excision as defined above, but approximate the outer gap edge with an ellipse having semimajor axis r_p + Δ, and an eccentricity and apsidal orientation estimated by eye from snapshots (for a sampling, jump to Figures 10 and 11).

Each simulation runs until Σ_gap appears to have converged in time; see Figure 5 for an example. The time required to reach convergence scales approximately as the viscous timescale $r_{\rm p}^{2}/\nu$, shortening somewhat with larger q. Each value of Σ_gap that we report in Table 1 is averaged over a time interval that starts at time t_conv near the end of the simulation and that lasts for duration Δt. For many models, there is actually no need to time-average because the time variability in Σ_gap is fractionally small (less than a few percent). However, some models exhibit greater variability, particularly for q approaching 10⁻² or h/r = 0.1. The fluctuations appear periodic, with periods ranging from 0.5–1 P_p and amplitudes up to order unity. For these more strongly time-variable cases, we also record in Table 1 the maximum and minimum values of Σ_gap that occur during the averaging interval.

Zoom In Zoom Out Reset image size

Table 1. Simulated Gap Depths

q	α	h/r	PEnGUIn	ZEUS90	PEnGUIn		ZEUS90		PEnGUIn	ZEUS90
			Σgapa (Σ0)		tconvb (Pp)	Δtb (Pp)	tconv(Pp)	Δt(Pp)	Comments
1 × 10−4	10−3	0.05	4.6 × 10−1	4.9 × 10−1	20000	10	10000	10		d
2 × 10−4	10−3	0.05	1.9 × 10−1	2.0 × 10−1	20000	10	10000	10		d
5 × 10−4	10−3	0.05	3.0 × 10−2	3.1 × 10−2	20000	10	10000	10		d
1 × 10−3	10−3	0.05	4.0 × 10−3	3.9 × 10−3	20000	10	10000	10
2 × 10−3	10−3	0.05	$8.5^{+0.4}_{-0.3}\times 10^{-4}$	$1.2^{+0.1}_{-0.2}\times 10^{-3}$	6000	10	6000	10		e
5 × 10−3	10−3	0.05	$2.1^{+0.2}_{-0.2}\times 10^{-4}$	2.6$^{+0.5}_{-0.8}\times 10^{-4}$	6000	10	6000	10	f	e
1 × 10−2	10−3	0.05	$1.7^{+0.4}_{-0.5}\times 10^{-4}$	$1.4^{+0.7}_{-0.5}\times 10^{-4}$	6000	10	6000	10	f
1 × 10−3	10−2	0.05	1.4 × 10−1	1.5 × 10−1	2000	10	2000	10
2 × 10−3	10−2	0.05	2.7 × 10−2	2.8 × 10−2	2000	10	2000	10
5 × 10−3	10−2	0.05	$5.5^{+0.9}_{-1.0}\times 10^{-3}$	$4.6^{+2.6}_{-1.8}\times 10^{-3}$	2000	10	2000	10
1 × 10−2	10−2	0.05	$2.0^{+0.6}_{-0.8}\times 10^{-3}$	$2.7^{+0.6}_{-0.4}\times 10^{-3}$	2000	10	2000	10	e
1 × 10−3	10−1	0.05	6.1 × 10−1	7.2 × 10−1	300	10	300	10
2 × 10−3	10−1	0.05	3.5 × 10−1	4.9 × 10−1	300	10	300	10
5 × 10−3	10−1	0.05	9.8 × 10−2	1.7 × 10−1	300	10	300	10
1 × 10−2	10−1	0.05	3.8 × 10−2	$4.9^{+0.2}_{-0.2}\times 10^{-2}$	300	10	300	10
1 × 10−3	10−3	0.1	2.2 × 10−1	$2.4^{+0.1}_{-0.1}\times 10^{-1}$	7000	10	7000	10
5 × 10−3	10−3	0.1	$1.5^{+0.5}_{-0.1}\times 10^{-2}$	$1.6^{+0.1}_{-0.1}\times 10^{-2}$	6000	1000	7000	10	c
1 × 10−2	10−3	0.1	$8.4^{+0.7}_{-0.6}\times 10^{-3}$	$1.0^{+0.1}_{-0.1}\times 10^{-2}$	6900	100	7000	10	c
5 × 10−4	10−3	0.04	6.0 × 10−3	6.5 × 10−3	17000	10	10000	10		d
2 × 10−3	10−3	0.04	$2.7^{+0.2}_{-0.2}\times 10^{-4}$	$2.9^{+0.6}_{-0.6}\times 10^{-4}$	10000	10	5000	10		e
2 × 10−3	10−2	0.04	$7.9^{+1.1}_{-0.9}\times 10^{-3}$	$7.0^{+0.4}_{-0.5}\times 10^{-3}$	4000	10	4000	10

Notes. ^aAveraged over time and over a partial annulus centered on the planet, as defined in Section 2.3.4 and delineated in Figures 3 and 4. For visibly eccentric outer disks (see Comments column), the outer edge of the measurement annulus is made eccentric to conform to the gap shape (e.g., Figures 10 and 11). Maximum and minimum values of Σ_gap are given for runs for which these values deviate from the time-averaged value by more than a percent. All surface densities are in units where Σ(r = 1) = 1 in a steadily accreting, planet-less disk. ^bThe time t_conv is taken near the end of a simulation, when Σ_gap appears to have nearly converged to its steady-state value (see Figure 5). Formally t_conv marks the beginning of the time interval, of duration Δt, over which Σ_gap is averaged. All times listed are in units of the planetary orbital period, P_p. ^cHighly unsteady outer gap edge (e.g., Figure 9) and long-term time variability. A longer Δt is chosen to capture the variability. ^dIndirect potential included (Section 2.2). ^eEccentric outer disk (e.g., Figures 10 and 11). Measured eccentricities are ∼0.10–0.15 and apsidal precession periods are ∼300–600P_p. ^fAn eccentric outer disk is observed at t ∼ 1000P_p, but the eccentricity damps away by t_conv. The damping is probably artificial (Section 3.2).

Download table as:  ASCIITypeset image

Gap depths Σ_gap(Σ₀) are recorded in Table 1 and plotted against each of the parameters q, α, and h/r in Figures 6, 7, and 8, respectively. Overall the agreement between the two codes, which utilize completely different algorithms, is remarkably good.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 attests that for q ≲ 5 × 10⁻³, gap depths scale roughly as q⁻²—as our analytic scaling (3) predicts. For q ≳ 5 × 10⁻³, the curves flatten somewhat (more on the behavior at large q in Section 3.2). In Figure 7, Σ_gap appears to scale with α to a power between 1 and 1.5. For comparison, our analytic scaling (3) predicts Σ_gap∝α¹. The empirical dependence on h/r is similarly steeper than the analytic dependence: Equation (3) predicts that Σ_gap∝(h/r)⁵ whereas Figure 8 shows that the power-law indices vary between 5 and 7.

We obtain a "best-fit" power-law relation by minimizing the function y = ∑[ln Dq^Aα^B(h/r)^C − ln Σ_gap]²/N over the parameters (A, B, C, D). The sum is performed over N data points, excluding the discrepant runs at q = 0.01 and runs for which Σ_gap/Σ₀ > 0.2 (i.e., runs for which gaps hardly open). With these exclusions, there are N = 13 data points from PEnGUIn, best fitted by

$$\begin{eqnarray} && 10^{-4} \le \, q \le 5 \times 10^{-3} \,: \nonumber\\ &&\quad\Sigma _{\rm gap}/\Sigma _0 = \, 0.14 \left(\frac{q}{10^{-3}} \right)^{-2.16}\left(\frac{\alpha }{10^{-2}} \right)^{1.41} \left(\frac{h/r}{0.05} \right)^{6.61} \,.\qquad \end{eqnarray} \tag{ 12 }$$

The N = 13 points from ZEUS90 are best described by a very similar formula:

$$\begin{eqnarray} && 10^{-4} \le \, q \le 5 \times 10^{-3} \,: \nonumber\\ &&\quad\Sigma _{\rm gap}/\Sigma _0 = \, 0.15 \left(\frac{q}{10^{-3}} \right)^{-2.12}\left(\frac{\alpha }{10^{-2}} \right)^{1.42} \left(\frac{h/r}{0.05} \right)^{6.45} \,.\qquad \end{eqnarray} \tag{ 13 }$$

Both of these relations fit their respective data to typically better than 20%; the largest deviation in the PEnGUIn fit is 40%, corresponding to (q, α, h/r) = (2 × 10⁻³, 10⁻³, 0.05), and for ZEUS90 the largest deviation is 50%, corresponding to (q, α, h/r) = (10⁻³, 10⁻³, 0.05).

Because our two codes agree so well, and because the fits are good, we are confident the deviations between our empirical scaling (say, Equation (12) from PEnGUIn) and our analytic scaling (3) are real and reflect physical effects not captured by our analytic scaling. And because our analytic scaling (3) matches exactly the scalings found numerically by Duffell & MacFadyen (2013) at low q ≲ 10⁻⁴—whereas our empirical relation (12) applies to high q ≳ 10⁻⁴—these physical effects manifest for giant (Jupiter-like) planets, not lower-mass (Neptune-like) planets.⁷ We have not elucidated what this physics is, although there might be some clues from the gap behavior at the very highest values of q we tested, as discussed in Section 3.2.

At the same time, we emphasize that the deviations between Equations (3) and (12), though (probably) real, are not large. If we insist on fitting the PEnGUIn data using Equation (3)—i.e., if we fix (A, B, C) = (− 2, 1, 5) and allow only the coefficient D to float—then the data deviate from Equation (3) by typically a factor of two, and at most a factor of three. Thus, the physical effects not captured by our analytic scaling, whatever they are, do not lead to order-of-magnitude changes in gap depth, at least over the range of parameters tested.

The dependence of Σ_gap on q flattens at the highest values of q considered (Figure 6). Fitting the N = 8 points from PEnGUIn for which q ⩾ 5 × 10⁻³ and Σ_gap/Σ₀ < 0.2 yields

$$\begin{eqnarray} && 10^{-2} \ge \, q \ge 5 \times 10^{-3} \,:\, \Sigma _{\rm gap}/\Sigma _0 =\nonumber\\ \ms &&\quad 4.7 \times 10^{-3} \left(\frac{q}{5 \times 10^{-3}} \right)^{-1.00}\left(\frac{\alpha }{10^{-2}} \right)^{1.26} \left(\frac{h/r}{0.05} \right)^{6.12} \,.\qquad \end{eqnarray} \tag{ 14 }$$

Similarly for the N = 8 points from ZEUS90 we obtain

$$\begin{eqnarray} 10^{-2}& \ge & q \ge 5 \times 10^{-3} \,:\,\Sigma _{\rm gap}/\Sigma _0 = 5.6 \times 10^{-3} \nonumber\\ \ms &&\times \left(\frac{q}{5 \times 10^{-3}} \right)^{-1.02}\left(\frac{\alpha }{10^{-2}} \right)^{1.34} \left(\frac{h/r}{0.05} \right)^{6.12}.\qquad \end{eqnarray} \tag{ 15 }$$

Although at first glance one might attribute the flattening of the trend of Σ_gap with q to the onset of strong shocks, we do not believe this is the correct interpretation. In the strong-shock regime, where disturbances excited by the planet are non-linear at launch, the torque exerted on the disk by the planet scales as q¹(h/r)⁰ (e.g., Hopkins & Quataert 2011, their Section 2.3).⁸ Then the same arguments in Section 1.1 yield Σ_gap∝q⁻¹(h/r)². This analytic relation reproduces the scaling index for q given by our empirical relations (14) and (15), but fails to reproduce the empirical scaling index for h/r. Furthermore, the flattening begins at an apparently "universal" q-value of ∼5 × 10⁻³ that is independent of h/r, whereas in the strong-shock interpretation, the critical q-value should scale as (h/r)³ (i.e., the expected critical q is given by the so-called thermal mass).

We do not have an explanation for the flatter slope of −1 at high q. We speculate that it might be caused by the most massive companions causing material at the gap edge to "leak" into the gap. The most massive planets disturb the gap edge so strongly that local instabilities send streamers of gas into the gap. These streamers, which de Val-Borro et al. (2006) called "filaments," are prominent in the high-q snapshots in Figure 9 (and can be seen even at q = 10⁻³ in the PEnGUIn snapshot in Figure 3). The filaments appear to originate from unsteady structures along gap edges; similar structures were seen by, e.g., Kley & Dirksen (2006, their Figures 1 and 7).

Zoom In Zoom Out Reset image size

We also observe evidence for an eccentric outer disk at high q; see Figure 10 and the "Comments" column in Table 1. The outward transport of angular momentum by waves launched at the 1:3 outer eccentric Lindblad resonance pumps the eccentricity of the outer disk (Lubow 1991; Papaloizou et al. 2001). The q-value for which disks become eccentric depends on α, h/r, and disk mass (Lubow 1991; Papaloizou et al. 2001; Kley & Dirksen 2006; D'Angelo et al. 2006; Dunhill et al. 2013). For a planet held on a fixed circular orbit embedded in a non-gravitating disk for which h/r = 0.05 and α ≈ 0.005, Kley & Dirksen (2006) found that q ≳ 0.003 led to eccentric disks; for α ≈ 0.01, the required q ≳ 0.005. Their findings are in line with ours.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In two runs with PEnGUIn, an eccentricity appears in the outer disk at early times but damps away by the time Σ_gap converges (see Table 1). The eccentricity damping is probably an artifact of our outer circular boundary at r_out = 2.5r_p, which according to Kley & Dirksen (2006) is too close to the planet to properly simulate eccentric disks. The danger posed by the outer boundary should lessen as the α-viscosity increases and disturbances excited by the planet are more localized; this may be why circularization occurs only for our lowest α = 10⁻³ runs at high q.

We note that the fine-structure filaments threading the gaps are seen in most of our q ⩾ 0.001 simulations, while only a few of these cases evince eccentric outer disks. Moreover, the filaments observed by de Val-Borro et al. (2006) do not appear associated with disk eccentricity. It is unclear to us whether the filaments and disk eccentricity are directly related.

Generally PEnGUIn and ZEUS90 agree very well on Σ_gap (e.g., Figures 6–8), typically differing by no more than a few tens of percent, and often much better. Because the codes rely on fundamentally different algorithms—one is a shock-capturing Lagrangian-remap code, while the other is an Eulerian code using the upwind method—their agreement lends confidence in the accuracy of our results. Some minor, systematic differences include: (1) Σ_gap for ZEUS90 is larger than for PEnGUIn; (2) PEnGUIn usually resolves a higher density peak near the planet; (3) PEnGUIn typically requires a longer time to converge; and (4) near the outer boundary where the resolution is coarser, ZEUS90 has difficulty relaxing to the steadily accreting solution described by Equations (9)–(11), deviating from the correct Σ by up to 50%. All of these differences may be attributed to the fact that PEnGUIn uses PPM, which is a fourth-order method for uniform grids (third-order for non-uniform grids), whereas ZEUS90's algorithm is only second-order in space. Thus, PEnGUIn tends to be more accurate and less numerically diffusive than ZEUS90, at the cost of taking a longer time to resolve sharp features.

A key innovation of PEnGUIn is its use of GPU technology to accelerate computations. PEnGUIn's speed on a single GTX-Titan graphics card can rival that of a traditional CPU cluster having ∼100 cores. For this paper we ran PEnGUIn on a desktop computer housing three graphics cards. With specialized hardware we can connect up to eight GPUs to a motherboard. PEnGUIn's scalability with the number of cards approaches linear as the resolution increases; our speed on three cards is 2.33 × that of a single card at our standard resolution; if the resolution is doubled (quadrupled), the speed enhancement factor increases to 2.64 (2.92) as PEnGUIn takes more full advantage of GPU's parallelism. Currently PEnGUIn can run on a single node only, and would need to be modified to run on multiple nodes. In a multi-node cluster of GPUs, the scaling of speed with the number of cards per node is unlikely to be linear because communication between nodes is significantly slower than between cards on a single node.

Are the gaps empty enough to reproduce the low optical depths characterizing the cavities of transitional and pre-transitional disks? Surface density contrasts from models of disks like PDS 70 (Dong et al. 2012) and GM Aur (Calvet et al. 2005) are 10³ or more. We have found that certain sets of planet-disk parameters can achieve such contrasts. For example, (q, α, h/r) = (5 × 10⁻³, 10⁻³, 0.05) produces contrasts of ∼5000 (Table 1). Lower mass planets could also be made to work with lower α-viscosities and/or cooler disks with lower h/r; as a further example, (q, α, h/r) = (2 × 10⁻³, 10⁻³, 0.04) generates a contrast of ∼3000. The dependence on disk temperature is especially sensitive: Σ_gap∝(h/r)^6.6∝T^3.3.

The surface density contrasts reported in this paper are all underestimates insofar as we have neglected accretion onto the planet, but arguably the disk accretion rate cannot be reduced by more than a factor of order unity, lest the planet starve the host star and violate observed stellar accretion rates (Zhu et al. 2011; see also Lubow & D'Angelo 2006). We would argue further that our gas surface density contrasts are also underestimates of dust surface density (i.e., optical depth) contrasts, to the extent that mechanisms like dust filtration at outer gap edges (e.g., Zhu et al. 2012) deplete dust relative to gas in gaps.

Given these findings, we feel that when it comes to transition disks, the problem is not so much gap depth, but gap width. A single planet embedded in an accreting disk generates a gap too narrow in radial width to explain the expansive cavities observed in transition disks. To connect to observations would seem to require that we expand our study to include multiple planets or brown dwarfs within a viscous, gravitating disk, as has been done by Zhu et al. (2011). These authors discounted α ≲ 0.002—and were therefore compelled to invoke additional channels of opacity reduction (e.g., grain growth) to explain transition disks—because multiple planets were found to be dynamically unstable at low α/high Σ (see their page 8). The incompatibility of multiple giant planets with low α is a result that we would like to see confirmed independently and further explored.

All our empirical scalings for Σ_gap suggest that arbitrarily low values of α generate arbitrarily clean gaps. We expect, however, the scalings to break down for small enough α. Without an intrinsic disk viscosity, a planet may stir the gap edge in such a way as to trigger local instabilities and turbulent diffusion. For example, with α = 0, the outer gap edge might be so sharp as to be Rayleigh unstable. There should therefore be a minimum value or "floor" on Σ_gap caused by a minimum planet-driven viscosity. Just such a floor has been reported by Duffell & MacFadyen (2013); see their Figure 7. Other simulations of planets in inviscid disks, concentrating on orbital migration and not gap depth, have been carried out by Li et al. (2009) and Yu et al. (2010).

Similar arguments suggest that there might also be a floor on Σ_gap at high q. The streamers/filaments that invade the gap are densest for the highest q-values we tested.

In Section 1.1, we presented an analytic derivation of gap depth Σ_gap. We discovered that the power-law scalings in our analytic relation (3) match precisely those reported from numerical experiments for low-mass planets by Duffell & MacFadyen (2013). Our analytic relation can even reproduce approximately (deviating systematically by factors of a few) the empirical results we found for gaps carved by Jupiter-like planets.

The success of our breezy analytic derivation for Σ_gap is surprising. Our derivation is "zero-dimensional" ("0D") because it considers only the total rates of angular momentum transport—a.k.a. the total "angular momentum luminosity" or "angular momentum current"—integrated over all azimuth and radius; it ignores the complicated details of how the torques are actually applied differentially in space. We have already noted in Section 1.1 how it is not completely obvious that Σ₀ characterizes the viscous torques in gap edges. Furthermore, our 0D treatment considers only the wave and viscous contributions to the total angular momentum current and neglects the contribution from advection (i.e., the transport of angular momentum associated with a non-zero radial velocity v_r) and the contribution from azimuthal pressure variations (Crida et al. 2006).

We can address some of these problems and make the leap from 0D to 1D by considering the azimuthally averaged, radially dependent torque balance equation for a steady-state accretion disk perturbed by a planet (see, e.g., Equation (3) of Lubow & D'Angelo 2006):

$$\begin{equation} \frac{d(3\nu \Sigma \Omega r^2)}{dr} = -\Sigma \Omega r^2 v_r + 2 \Sigma r \Lambda (r), \end{equation} \tag{ 16 }$$

where the Lindblad torque per unit mass is

$$\begin{equation} \Lambda (r) = {\rm sgn} \, (r - r_{\rm p}) \frac{f GM_\ast q^2}{2r} \,\left(\frac{r}{r-r_{\rm p}} \right)^4 \end{equation} \tag{ 17 }$$

and f is an order-unity constant. The left-hand side of Equation (16) accounts for the viscous torque, while the first term on the right-hand side accounts for angular momentum transport by advection. We define x ≡ r − r_p and $\dot{M} \equiv 2\pi \Sigma v_r r =$ constant  <  0 and approximate Ω = Ω_p(r/r_p)^−3/2 so that νΩr = ν_pΩ_pr_p = constant (variables subscripted by p take their values at the planet's orbital radius). For the outer disk, Equation (16) simplifies to

$$\begin{equation} \alpha \left(\frac{h}{r} \right)^2 \frac{d(r\Sigma)}{dr} = \frac{|\dot{M}| \Omega _{\rm p} r_{\rm p}}{6 \pi GM_\ast } \left(\frac{r}{r_{\rm p}} \right)^{-1/2} + \frac{f}{3} \,q^2 \, \left(\frac{r}{x} \right)^4 \, \Sigma \,. \end{equation} \tag{ 18 }$$

From this equation it becomes apparent how the outer accretion disk responds when repelled outward by a planet: for q > 0, the surface density gradient dΣ/dr steepens, just enough that the viscous torque can exceed the Lindblad torque and maintain a steady flow of mass inward (i.e., carry the $\dot{M}$ imposed at infinity across the planet's orbit).

Setting $\dot{M}=0$, and working in the WKB limit where d/dr ≫ 1/r, gives the standard zero-inflow solution: an exponential profile for Σ(r), commonly used in the literature (e.g., Lubow et al. 1999; de Val-Borro et al. 2007; Mulders et al. 2013). Keeping $\dot{M} \ne 0$ alters Σ(r): it still resembles an exponential but is shifted outward (for fixed f), as Figure 12 demonstrates. We find for the parameters chosen that Equation (18) describes well the gap profile from our 2D simulations, down to a distance of ∼3–4 Hill radii away from the planet. But inside this cutoff distance, the 1D solution fails critically—it falls much too steeply to recover the actual flat-bottomed gap.

Zoom In Zoom Out Reset image size

The problem of determining the gap depth analytically in 1D appears tantamount to the problem of understanding what happens inside this cutoff distance. Lindblad torques shut off here, the tidal gravitational field of the planet is especially strong, and circulating streamlines give way to horseshoe orbits. One-dimensional analytic treatments may be inadequate to the task of modeling how gas navigates from the outer disk to the inner disk through a series of "horseshoe turns" (Lubow et al. 1999; Kley & Nelson 2012). As far as analytic treatments go, it may be that to do better than 0D requires at least 2D.