DISK–SATELLITE INTERACTION IN DISKS WITH DENSITY GAPS

We study the tidal coupling of a planet with a non-uniform disk in the shearing sheet geometry (Goldreich & Lynden-Bell 1965), which allows us to neglect geometric curvature effects while preserving the main qualitative features of the system. We also neglect the vertical dimension and assume the disk to be two-dimensional. The dynamics of fluid is then governed by the following equations⁵ (e.g., Narayan et al. 1987, hereafter NGG):

$$\begin{equation} \frac{\partial {\bf v}}{\partial t}+ \left({\bf v}\cdot \nabla \right){\bf v}+ 2{\bf \Omega }\times {\bf v} +4A\Omega x{\bf e}_x =-\frac{\nabla P}{\Sigma }-\nabla \Phi _p, \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{\partial \Sigma }{\partial t}+\nabla \cdot \left({\bf v}\Sigma \right)=0, \end{equation} \tag{ 6 }$$

where v = (u, v) is the fluid velocity with components in the x ≡ r − r_p (radial) and y (azimuthal) directions correspondingly, P is the gas pressure, Φ_p is the planetary potential, Ω_p ≡ Ω_pe_z is the Keplerian rotation frequency at the planet position (r = r_p; we will subsequently drop the subscript "p" for brevity), and A ≡ (r/2)(dΩ/dr) is the shear rate at the same location. Additionally, in our notation c_s is the sound speed and B = Ω  +  A is Oort's B constant.

For an isothermal disk with no planet present (ϕ ≡ 0) and for an arbitrary background density profile Σ = Σ₀(x), Equations (5) and (6) have the following exact steady-state solution:

$$\begin{eqnarray} &&u_{0}=0, \quad v_{0}=2Ax+\frac{c_s^{2}}{2\Omega } \frac{\partial \ln \Sigma _{0}}{\partial x}. \end{eqnarray} \tag{ 7 }$$

After having specified our steady-state solution {Σ₀, u₀, v₀}, we introduce the planetary potential as a perturbation and study the linear behavior of the perturbed density Σ₁ and velocity field $\delta \textbf {\em v}=(u_1,v_1)$. By ignoring the quadratic perturbed quantities that appear in Equations (5) and (6), one can get the following general system:

$$\begin{eqnarray} \frac{\partial u_1}{\partial t}&+&\left(2Ax+\frac{c_s^{2}}{2\Omega } \frac{\partial \ln \Sigma _{0}}{\partial x}\right)\frac{\partial u_1}{\partial y} -2\Omega v_1 \nonumber \\ & +& c_s^{2}\frac{\partial \left(\Sigma _1/\Sigma _0\right)}{\partial x}=- \frac{\partial \Phi _p}{\partial x} \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \frac{\partial v_1}{\partial t}&+&\left(2Ax+\frac{c_s^{2}}{2\Omega } \frac{\partial \ln \Sigma _{0}}{\partial x}\right)\frac{\partial v_1}{\partial y} \nonumber \\ & +&\left(2B+\frac{c_s^{2}}{2\Omega } \frac{\partial ^{2} \ln \Sigma _{0}}{\partial x^{2}} \right) u_1 + c_s^{2}\frac{\partial \left(\Sigma _1/\Sigma _0\right)}{\partial y} =-\frac{\partial \Phi _p}{\partial y} \quad\quad \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \frac{\partial \left(\Sigma _1/\Sigma _0\right)}{\partial t}&+& \left(2Ax+\frac{c_s^{2}}{2\Omega } \frac{\partial \ln \Sigma _{0}}{\partial x}\right)\frac{\partial \left(\Sigma _1/\Sigma _0\right)}{\partial y} \nonumber \\ & +& u_1\frac{\partial \ln \Sigma _{0}}{\partial x}+ \frac{\partial u_1}{\partial x}+\frac{\partial v_1}{\partial y}=0. \end{eqnarray} \tag{ 10 }$$

We then represent all perturbed fluid variables via Fourier integrals as {Σ₁, u₁, v₁, Φ_p} = (2π)⁻¹∫^∞_∞exp (ik_yy){δΣ, u, v, ϕ}, which takes care of the y-dependence, leaving us with the following set of equations in the x-coordinate only. Thus, by defining

$$\begin{equation} \sigma (x)=-2Ak_{y}x-\frac{c_s^{2}k_{y}}{2\Omega } \frac{\partial \ln \Sigma _{0}}{\partial x}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \tilde{B} = B +\frac{c_s^{2}}{4\Omega } \frac{\partial ^{2} \ln \Sigma _{0}}{\partial x^{2}}, \end{equation} \tag{ 12 }$$

we get the linearized system in steady state (∂/∂t ≡ 0)

$$\begin{equation} {-}i\sigma u-2\Omega v +c_s^{2}\frac{\partial \left(\delta \Sigma /\Sigma _0\right)}{\partial x} =- \frac{\partial \phi }{\partial x} \end{equation} \tag{ 13 }$$

$$\begin{equation} {-}i\sigma v+2 \tilde{B} u +i k_{y}c_s^{2}\frac{\delta \Sigma }{\Sigma _{0}} = - ik_{y} \phi \end{equation} \tag{ 14 }$$

$$\begin{equation} {-}i\sigma \frac{\delta \Sigma }{\Sigma _{0}}+u\frac{\partial \ln \Sigma _{0}}{\partial x}+ \frac{\partial u}{\partial x}+ik_{y}v = 0, \end{equation} \tag{ 15 }$$

where the Fourier component of the planetary potential produced by a point mass at the origin and its derivative are given by

$$\begin{equation} \phi (k_y,x) = -\frac{GM_p}{\pi }K_0(|k_yx|), \end{equation} \tag{ 16 }$$

$$\begin{equation} \frac{\partial \phi }{\partial x} = \mbox{sgn}(x)\frac{GM_p}{\pi }k_y K_1(|k_yx|). \end{equation} \tag{ 17 }$$

Solving Equations (13)–(15) for the azimuthal velocity perturbation v, we obtain the second-order ordinary differential equation

$$\begin{eqnarray} \frac{\partial ^{2} v}{\partial x^{2}} &+& \frac{\partial v}{\partial x}\frac{\partial \ln \left(\Sigma _{0}/\Delta \right)}{\partial x} +v\left[\! \frac{\sigma ^{2}}{c_s^{2}}-\frac{\Omega }{c^{2}\tilde{B}}\Delta + \frac{k_{y}\sigma \Lambda }{2\tilde{B}} \frac{\partial \ln \left(\sigma \Lambda /\Delta \right)}{\partial x} \!\right] \nonumber \\ &&= {-}\frac{\partial \phi }{\partial x} \frac{2\tilde{B}}{c_s^{2}} +\phi \left[ \frac{k_{y}\sigma }{c_s^{2}} +\frac{k_{y}^{2}\Lambda }{2\tilde{B}}\frac{\partial \ln \left(\Lambda /\Delta \right)}{\partial x}\right], \end{eqnarray} \tag{ 18 }$$

where we have defined

$$\begin{equation} \Lambda = 1-\frac{2}{k_{y}\sigma }\left(\frac{\partial \tilde{B}}{\partial x}- \tilde{B}\frac{\partial \ln \Sigma _{0}}{\partial x}\right), \end{equation} \tag{ 19 }$$

$$\begin{equation} \Delta = c_s^2k_{y}^{2}\Lambda +4\tilde{B}^{2}. \end{equation} \tag{ 20 }$$

The perturbed radial velocity u and the surface density perturbation δΣ are directly expressed in terms of v via the following relations:

$$\begin{eqnarray} && u=-\frac{i}{\Delta }\left(k_y c_s^2\frac{\partial v}{\partial x}-2 \tilde{B} \sigma v + 2\tilde{B}k_y\phi \right), \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} && \delta \Sigma =\frac{\Sigma _0}{\Delta }\left(2\tilde{B}\frac{\partial v}{\partial x} +k_y\sigma \Lambda v - k_y^2 \Lambda \phi \right). \end{eqnarray} \tag{ 22 }$$

One can easily check that in the limit ∂Σ₀/∂x → 0 Equation (18) reduces to Equation (5) of RP12, which was derived for a uniform disk. Note also that the formalism presented here can be easily extended to incorporate temperature gradients in the disk by adding corresponding terms to the unperturbed azimuthal velocity v₀ in Equation (7).

We model the density structure of an axisymmetric gap around the satellite orbit by a particular three-parameter density profile (symmetric with respect to x = 0)

$$\begin{eqnarray} &&\frac{\Sigma _{0}(x)}{\Sigma _\infty }=1-\frac{1-\Sigma _{\rm min}}{1+(x/\delta)^{n} \exp [{-}(2.5\delta /x)^{2}]}, \end{eqnarray} \tag{ 23 }$$

where the surface density is normalized to its value at infinity Σ_∞. In this profile the dimensionless parameter Σ_min controls the depth of the gap, δ regulates its width, while an even positive integer n sets the steepness of the density gradient at the gap edge. The density variation between Σ_minΣ_∞ inside the gap and Σ_∞ outside of it is predominantly due to the power-law dependence on x in the denominator of the second term on the right-hand side of Equation (23). The exponential dependence in the denominator is introduced only to ensure a flat bottom at x ≲ δ and allows us to avoid large density gradients close to the origin. This simplifies our subsequent analysis and helps us get important contributions to the torque from just outside the gap, exactly where we want to test our theory. The factor 2.5 in the exponential is introduced for the parameter δ to better correspond to the width of the inner, flat part of the gap.

In the left panel of Figure 1 we plot the density profile for the disk model with parameters δ = 5h, Σ_min = 3 × 10⁻⁴, and n = 2, 4, 6, respectively. These profiles represent rather broad gaps: for the models with δ = 5h shown in Figure 1 one finds that the density Σ_∞/2 is reached at |x| ≈ 10h, |x| ≈ 8.5h, and |x| ≈ 7.7h for n = 2, 4, and 6, respectively. Evidently, as n is increased, the profile becomes boxier and the gradients at the gap edge (they are largest for 4h ≲ |x| ≲ 15h) grow.

Zoom In Zoom Out Reset image size

Non-zero pressure gradients intrinsic to inhomogeneous disks affect the unperturbed azimuthal velocity (7) and change the locations of the LRs, at which different harmonics of planetary potential couple to the disk, compared to the purely Keplerian case. In Appendix A we demonstrate that in the shearing sheet approximation a given k_y harmonic has the LR condition satisfied at x_L(k_y), satisfying the following implicit relation:

$$\begin{eqnarray} &&k_y(x_L)=\pm \frac{2}{3}\frac{\kappa (x_L)}{\Omega _p} \left(x_L-\frac{h^2}{3}\frac{\partial \ln \Sigma _{0}}{\partial x}\Big |_{x_L}\right)^{-1}, \end{eqnarray} \tag{ 24 }$$

where the modified value of the local epicyclic frequency κ for the inhomogeneous disk is given by (see Appendix A)

$$\begin{eqnarray} &&\kappa (x)=\Omega _p\left(1+ h^2\frac{\partial ^{2} \ln \Sigma _{0}}{\partial x^{2}}\right)^{1/2}. \end{eqnarray} \tag{ 25 }$$

In a uniform disk (∂ln Σ₀/∂x = 0) one recovers the usual expressions κ = Ω_p and x_L(k_y) = 2/(3k_y).

For a class of density profiles (23) considered in this work the modification to the x_L(k_y) relation (24) compared to the purely Keplerian case comes predominantly from the change of the local epicyclic frequency as given by Equation (25). For this reason we explore the dependence of κ on x in some detail and plot it in the right panel of Figure 1 for the density profiles shown in the left panel of the same figure.

Deep inside the gap Σ₀ is almost constant because of the assumed exponential dependence in Equation (23) and κ(x) → Ω_p. As the gap edge is approached, at |x| ≳ 4h, ∂²ln Σ₀/∂x² increases and κ does the same: in all displayed models κ goes up to ≳ 2Ω_p. Beyond |x| ≳ 5h the density profiles change their curvature from positive to negative and κ becomes smaller than Ω_p. At even larger |x| the density profile flattens out and κ converges to Ω_p at |x| ≳ 15h. The deviations of the epicyclic frequency compared to its value in a uniform disk become more extreme as we increase n (increase boxiness of Σ₀(x)), with the differences between models being more dramatic in regions where κ < Ω_p. The minimum values of κ are 0.36Ω_p, 0.27Ω_p, and 0.03Ω_p for profiles with n = 2, 4, and 6, respectively; see Figure 1.

In fact, one can see from Equation (25) that, depending on the curvature of Σ₀, the gap profile could have κ² < 0 at the gap edge, which would mean that the disk is Rayleigh unstable. For this reason in this work we ignore profiles with extremely deep gaps and strong density gradients at their edges since these are likely to result in κ² < 0 at some place. Our boxiest model with n = 6 might be regarded as a limiting profile (for assumed values of δ and Σ_min) since the epicyclic frequency closely approaches zero at x ≈ 6h.

Armed with this knowledge, we plot in Figure 2 the dependence (24) for the density model given by Equation (23) with the same set of parameters as in Figure 1: δ = 5h, Σ_min = 3 × 10⁻⁴, and n = 2, 4, 6. By comparing these two figures, one can see that the resonances are displaced outward relative to the uniform density case in regions where κ > Ω_p, while the opposite happens when κ < Ω_p. As a result, there is a concentration of the resonances at the gap edge where this transition occurs.

Zoom In Zoom Out Reset image size

Indeed, in a uniform disk with κ = Ω_p the radial interval 4h < |x| < 6h contains LRs only for the modes with k_yh in the range (0.11, 0.17). At the same time, for the gap in the form (23) with n = 2 the same radial interval contains LRs for modes satisfying 0.05 ≲ k_yh ≲ 0.37; for n = 6 model modes with 0.009 ≲ k_yh ≲ 0.4 are being excited in the same radial interval. Thus, concentration of resonances at the gap edge is more pronounced for density profiles with larger gradients, meaning larger deviations of κ from Ω_p.

Another striking feature of the x_L(k_y) dependence clearly visible in Figure 2 is that even for not very sharp gaps (for rather low values of n) there are regions close to the gap edge where a single k_y mode can be excited at two or three different locations (for a fixed sign of x), as happens, e.g., for k_yh = 0.3 in this figure. This splitting of resonances in inhomogeneous disks has important implications for the torque behavior in Fourier space; see Section 5.

Having specified the gap density profile, the solution for the azimuthal velocity perturbation v is obtained by direct numerical integration of Equations (11), (12), and (16)–(20); knowing v, we find u and δΣ using Equations (21) and (22).

Numerical integration uses the same method as was employed in RP12 for uniform disks: we shoot numerical solutions from the origin and match them to the WKB outgoing waves far from the planet (see also Korycansky & Pollack 1993). As shown in NGG, the exact homogeneous solution to Equation (18) for a homogenous disk, i.e., the parabolic cylinder functions, can be well approximated by the WKB outgoing waves

$$\begin{eqnarray} &&v(x\rightarrow \pm \infty)\sim \sqrt{\frac{2}{3k_{y}x}}e^{\pm i3k_{y}x^{2}/4h}, \end{eqnarray} \tag{ 26 }$$

when x/h ≫ (8/9)[1 + (k_yh)²]/(k_yh)². Moreover, for this approximation to accurately represent an inhomogeneous solution of Equation (18) one requires xk_y ≫ 1, so that ϕ ≪ 1 and ∂ϕ/∂x ≪ 1.

To additionally account for the presence of a density gap in the disk, we use the fact that for every profile ∂Σ₀/∂x → 0 when x is several times the characteristic width of the gap δ and, therefore, the asymptotic solution will still be given by Equation (26). All conditions together imply that x/h ≫ min {1/(k_yh), 1, δ/h}, which is the criterion we use to match our numerical solutions to the outgoing waves.

Our numerical calculations use a fourth-order Runge–Kutta integrator with spatial resolution of h/400 and x ∈ [ − 220h, 220h] for 720 uniformly log-spaced values of k_yh between 0.01 and 15 and potential softening length of 10⁻⁴h (its value is not so important for us since we are considering relatively broad gaps).

In Figure 3 we show one of the outcomes of our calculations—the behavior of the surface density perturbation Σ₁(x, y) in physical space. We obtain it by first calculating δΣ using Equation (22) with the numerically determined v and then performing the inverse Fourier transform. In Figure 3 azimuthal cuts of Σ₁(x, y) at different values of x are normalized by (Σ₀(x)x)^1/2 to ensure similar wake amplitude at different radial separations from the perturber; this scaling is expected⁶ from the conservation of the AMF carried by the density wave (Goodman & Rafikov 2001; Rafikov 2002a).

Zoom In Zoom Out Reset image size

In panel (a) we display our results for a shallow gap with a minimum surface density at the satellite position of 0.1Σ_∞ (corresponding density profile is shown in panel (b)). For comparison we also display the wake profiles for a uniform density disk. Deep inside the gap, at x = h, where the Σ₀ is essentially constant in our model, we find our Σ₁ to exactly coincide with the uniform disk calculation, as expected (our normalization of Σ₁ takes care of the difference in amplitude). However, farther out from the planet the results start to differ. In particular, at x = 10h the azimuthal location of the wake in a non-uniform disk gets shifted by ≈h relative to the uniform case, but the differences in the shape of the density perturbation remain at the modest level.

The situation is a bit different for deep gaps, as illustrated in panel (c), where we plot azimuthal density cuts through the wake for a gap profile with Σ_min = 0.0003Σ_∞ (shown in panel (d)). Again, at x = h density profiles coincide for both uniform and non-uniform cases. However, at the gap edge (at x = 5h) there are now substantial differences between the profiles not only in the wake position but also in shape. The differences become more pronounced as the wave propagates even further: at x = 10h the peak value of Σ₁ is shifted by ∼4h compared to the uniform case and the overall wake shape is considerably distorted. This is not surprising since the surface density eigenfunctions at the edge of the gap in a non-uniform disk are considerably distorted compared to their uniform disk analogs; see Section 5.1.

Note that the amplitude of scaled Σ₁ in Figure 3 does not vary too much with x. Since the AMF conservation predicts that Σ₁∝[Σ₀(x)x]^1/2, it is clear that the amplitude of Σ₁ grows as the wake propagates out from the bottom of the gap. At the same time, the relative density perturbation Σ₁/Σ₀(x) scales as [x/Σ₀(x)]^1/2 and decreases as the wake climbs up the density contrast at the gap edge. This may have important implications for the wake damping, as discussed in Section 8.

To derive the asymptotic torque density behavior (4) in physical space, GT80 have explicitly assumed that the angular momentum associated with a particular low-order (k_yh ≲ 1) harmonic of the perturbing potential gets added to the density wave exactly at the LR corresponding to that mode, i.e., x_L(k_y) = 2/(3k_y). Dong et al. (2011b) have extended this assumption to modes with arbitrary values of k_y to obtain a so-called LR prescription for the AMF F^LR_H(x), which they used as a benchmark for comparing with their numerical results for uniform disks. This prescription is given by the following formula (Dong et al. 2011b):

$$\begin{equation} F_H^{\rm LR}(x) = \frac{1}{h}\int _{\mu _{{\rm min}}(x)}^\infty d\mu \frac{F_H(\mu)}{F_H^{\rm WKB}(\mu)}F_H^{\rm WKB}(\mu), \end{equation} \tag{ 30 }$$

$$\begin{eqnarray} F_H^{\rm WKB}(\mu)&=&\frac{4}{3}\Sigma _\infty \mu ^2 \left(\frac{G M_{\rm p}}{c_s}\right)^2 \nonumber \\ && \times \left[2K_0(2/3)+K_1(2/3)\right]^2, \end{eqnarray} \tag{ 31 }$$

where F_H(μ)/F^WKB_H(μ) is the ratio of the actual AMF contribution carried by a particular harmonic with μ ≡ k_yh far from the perturber to the value of the same quantity obtained using WKB approximation. This ratio has been computed numerically by GT80 for arbitrary μ; in the limit μ → 0 it tends to unity.

In Equation (30) μ_min(x) = (2/3)h/x is the minimum value of k_yh for which the position of a corresponding LR (in a uniform disk) satisfies x_L(k_y) < x. Thus, F^LR_H(x) represents the sum of the AMFs (measured at x → ∞) carried by all individual Fourier modes having their LRs lying at x_L(k_y) < x. Note that in computing μ_min(x) we do not use the positions of LRs accounting for disk non-uniformity (Equation (24))—we found this prescription for x_L(k_y) to provide F^LR_H inconsistent with numerical results; see Section 6.1.

The prescription (30) is formulated for a disk with constant density Σ_∞. We extend it to disks with gaps by first defining an LR torque density in a uniform disk as dT^LR(x)/dx|_u ≡ dF^LR_H(x)/dx and then using this torque density in Equation (3):

$$\begin{eqnarray} &&\frac{dT^{\rm LR}(x)}{dx}=\frac{\Sigma _0(x)}{\Sigma _\infty }\cdot \frac{dF_H^{\rm LR}(x)}{dx}. \end{eqnarray} \tag{ 32 }$$

Note that for |x| ≳ h one has μ_min ≲ 1 and dT^LR(x)/dx|_u reduces to the standard GT80 expression (4). Most studies of tidal coupling in non-uniform disks have used this asymptotic form of the prescription (32), i.e., Equation (3) with dT/dr|_u given by Equation (4), to characterize the torque density; use of Equation (30) simply allows us to extend this prescription to arbitrary values of x.

The prescription (30) is equivalent to assuming that each potential harmonic exerts torque only in the immediate vicinity of its LR, a notion that has been shown by RP12 to be incorrect in the shearing sheet approximation even for μ ≲ 1. In reality LRs have finite width, which leads to the nontrivial interference between them, resulting in the negative torque density phenomenon at |x| ≳ 3.2h (Dong et al. 2011b; RP12) not captured by the GT80 analysis. Despite this failure of the assumed discreteness of the LR, we still use F^LR_H(x) in this work to compare with our new results as it provides an interesting reference point.

In Figure 4 we show the torque density for individual Fourier harmonics $(dT/dx)_{k_y}(x)$ computed using Equation (28), at different radial separations x. These calculations assume our fiducial density profile (23) with parameters δ = 5h, Σ_min = 3 × 10⁻⁴, and n = 2. We show both the torque density per unit local surface density (i.e., $(dT/dx)_{k_y}(x)$ normalized by Σ₀(x), left panels) and the pure $(dT/dx)_{k_y}(x)$ (divided by constant Σ_∞, right panels). We also display $(dT/dx)_{k_y}|_u$ computed for a uniform disk (see Figure 3 of RP12), normalizing it by Σ_∞ in the left panels (which yields the torque density per unit local surface density in a homogeneous disk) and also multiplying by Σ₀(x)/Σ²_∞ (to incorporate the local density correction) in the right panels. In all cases torque density is properly normalized to make it dimensionless. In the left panels of this figure we indicate the position of the LRs with arrows, both for an inhomogeneous (see Equation (24)) and a uniform disk.

Zoom In Zoom Out Reset image size

In panels (a) and (b), we plot $(dT/dx)_{k_y}(x)$ for k_yh = 1. The corresponding LR lies at x/h = 2/3, in the flat bottom of the gap, and is thus the same for both uniform and fully non-uniform disks. Moreover, the waveform inside the gap perfectly matches that of the uniform disk since the profile has constant density there. Note that the torque density is not sharply localized at the LR position x_L(k_y) but is extended over a reasonably broad radial range around x_L, with quite pronounced (even though damping) oscillations of $(dT/dx)_{k_y}(x)$ clearly visible even for x ≫ x_L. This additionally confirms the finite and non-negligible width of resonances previously pointed out by Artymowicz (1993) and RP12.

At x > 4h, where Σ₀ starts to increase, we observe that the torque density decreases its amplitude quite abruptly compared to the uniform disk calculation; see panel (b). By looking at the shaded regions, one can note that the two calculations start differing at Σ₀ ∼ 10⁻³Σ_∞, where the density gradients become significant. This damping of the amplitude compared to the homogeneous disk case is the consequence of incorporating the disk non-uniformity in the fluid equations, also clearly visible for other values of k_yh at the gap edge (see below).

In panels (c) and (d) we show the torque density for k_yh = 0.32 mode, for which there are multiple LRs at x/h = 2.1, 3.8, and 4.5; see Figure 2 (the uniform disk has, of course, just a single one at x/h = 2/(3k_yh) ≈ 2.1). These are indicated with arrows in panel (c), and it can be seen from both panels that the extra resonances located at x ∼ 4h give rise to positive torque density contribution and overall significant phase shift compared to the uniform disk case. There is a hint of the resonances alternating sign between positive at x/h = 2.1, negative at 3.8, and then positive again at 4.5, consistent with the analytical theory that predicts that the sign of the torque for a given mode is determined by the derivative of D(r) = κ(r)² − m²[Ω(r) − Ω_p]², which changes signs at the extra resonances (GT80). However, the interference of different resonances does not allow us to disentangle their separate contributions and makes ambiguous any interpretation of $(dT/dx)_{k_y}(x)$ behavior via forcing produced at a set of discrete locations. On the other hand, the amplitude of the torque density oscillations far from the resonances is greatly reduced in the self-consistent calculation of $(dT/dx)_{k_y}(x)$ so that in panel (d) we even multiply the uniform disk curve by 0.05 for it to have an amplitude comparable to that of the fully non-uniform disk calculation.

The mode with k_yh ≈ 0.14 shown in Figures 4(e) and (f) has a single LR at x_L ≈ 5h, which coincides with the resonance position in a uniform disk (see Figure 2). Despite this coincidence, the waveforms for the uniform and non-uniform disk calculations are very different from each other in terms of both shape and the overall amplitude (note that the uniform disk curve in panel (f) has been multiplied by 0.1).

The mode with k_yh = 0.08 shown in Figures 4(g) and (h) has its (single) LR at x ≈ 5.1h, considerably shifted inward from the corresponding uniform disk resonance position at x = 8.3h. The torque density normalized by Σ₀(x) (see panel (g)) has a rather broad spatial distribution extending from x ∼ 5h so x ∼ 8h before starting to oscillate at an amplitude reduced roughly by a factor of two compared to the uniform disk calculation.

Finally, in panels (i) and (j) we show the mode k_yh = 0.04 for which x_L ≈ 16h lies where the density gradients are negligible. As a result, calculations in uniform and non-uniform disk cases perfectly coincide, demonstrating that far from gap torque excitation is well described by the uniform disk theory of RP12.

To summarize, the spatial distribution of the torque density for individual Fourier harmonics can be substantially modified in non-uniform disks in regions where the density gradients are large, compared to the case of a constant density disk. The differences arise both because of the displacement of the LRs and because of the different mathematical structure of the fluid equations in the non-uniform disk case.

We now study the behavior of the AMF as a function of k_y for x → ∞. Since, as mentioned in Section 4, $F_{H,k_y}(x)$ is identical (up to a constant offset) to the integrated torque, we look at the behavior of

$$\begin{eqnarray} &&T(k_y)=\int _0^\infty \left(\frac{dT}{dx}\right)_{k_y} dx \end{eqnarray} \tag{ 33 }$$

in Fourier space, where $(dT/dx)_{k_y}$ is given by our numerical calculations.

As a point of comparison we also use the variation of $F_H(k_y)\equiv F_{H,k_y}(x\rightarrow \infty)$ with k_y in a homogeneous disk (previously computed by GT80), which we additionally multiply by Σ₀(x_L(k_y)) to account for the disk non-uniformity (where we take x_L(k_y) = 2/(3k_y)). This re-scaling implies that each mode carries the AMF from a uniform disk corrected by the local density at the position where the excitation of this mode occurs in a uniform disk. This procedure is a version of the LR theory described in Section 4.1 in Fourier space, as it incorporates the two major ingredients from that prescription: the localized mode excitations at an LR and the AMF computed for a uniform density disk.

In Figure 5, we plot the torque from our calculations and the AMF from the LR theory for our fiducial model with δ = 5h, Σ_min = 3 × 10⁻⁴, and n = 2. The first thing we note is that both expressions coincide in the limits of long wavelengths k_yh < 0.05 and short wavelengths k_yh > 2. This is expected to happen since in both limits the satellite excites waves in the flat parts of the disk—far outside the gap for k_yh < 0.05 and deep inside the gap, in the flat part of the density profile for k_yh > 2. For that reason the positions of the resonances and the waveforms coincide in our non-uniform calculation and in a uniform disk. But in the intermediate range of k_y there are appreciable differences between the LR theory and our calculations, mainly due to effect of the shifted LRs as we discuss next.

Zoom In Zoom Out Reset image size

First, for k_yh ∼ 0.05–0.1 Figure 2 demonstrates that LRs are shifted toward the perturber and the excitation occurs in the lower density region of the disk compared to the case when resonances lie at x = 2/(3k_y) (as appropriate for homogeneous disk). Because of this reduction of the background surface density, the torque produced by modes in this range of k_y ends up being lower than the AMF computed according to the LR theory (which assumes x_L = 2/(3k_y)), even though the coupling to the perturbation potential is stronger at smaller |x|. Conversely, modes with k_yh ∼ 0.1–0.2 have their resonances shifted outward, into the higher density part of the gap profile, which explains their higher values of T(k_y) compared to the LR theory AMF.

At higher values of k_y there is a point where resonances split and take place at multiple locations, as discussed in Section 2.3. According to Figure 2, this is the case for k_yh ≳ 0.2, which explains the intricate shape of T(k_y) for these modes—a peak at k_yh = 0.2–0.3, where the LRs are split into three (see Figures 4(c) and (d) for illustration of the torque density behavior for these modes). The innermost resonance almost coincides with that in a uniform disk, but the other two are located further out, in a higher density region of the gap, explaining the spike of T(k_y) for these values of k_y.

Complicated structure of the T(k_y) curve around the peak at k_yh ∼ 0.32 does not have a straightforward explanation and is probably related to the intricate behavior of the torque density when the excitation is taking place at multiple locations at the gap edge.

In the top panel of Figure 6, we plot the torque density dT/dx (integrated over all k_y; see Equation (28)) per unit of the local surface density Σ₀(x) for our fiducial model and make a comparison with the LR theory behavior given by Equation (32). In the lower panel of Figure 6 we directly compare dT/dx and dT^LR/dx. One can clearly see a significant discrepancy between the two.

Zoom In Zoom Out Reset image size

The disagreement at |x| ≲ 5h, where the surface density is close to constant (the bottom part of the gap), is expected based on the RP12 results—one can see that our dT/dx is negative for x ∼ (3–4)h, which is just the manifestation of the negative torque density phenomenon (Dong et al. 2011b) not captured by the dT^LR/dx. Also, the overall shape of dT/dx and dT^LR/dx is different for |x| ≲ 3h—again well known from the uniform disk studies of Dong et al. (2011b) and RP12.

The discrepancies at |x| ≳ 5h where the density gradients are important can only be understood by fully accounting for the disk non-uniformity. According to GT80, for |x| ≳ h the LR theory predicts Σ⁻¹₀dT/dx ≈ 2.5(h/x)⁴ (in this section we refer to the torque density in units of Σ_∞G²M²_p/(hc_s)²); see Equation (4) and Section 4.1. However, our self-consistent non-uniform disk calculations suggest that the torque density decays exponentially ∝exp (− 0.8x/h) for |x| ≳ 5h, as shown in the upper panel of Figure 6. This result casts serious doubt on all torque prescriptions used in the literature that adopt the GT80 functional form to account for the tidal torque. The exponential falloff is clearly related to the presence of strong density gradients at the gap edge, but we could not find a simple qualitative explanation for this particular form of the torque decay. We just note here that this behavior of dT/dx has nothing to do with the presence of the exponential term in Equation (23), which is introduced only to guarantee the flat density profile in the bottom part of the gap. We provide more details on the exponential decay of dT/dx in Appendix B.

Apart from the discrepancy in the functional form, the self-consistent specific torque density is also more concentrated toward the perturber at the gap edge and has a larger amplitude there. This effect can be understood from our previous analysis of the shifted LRs that tend to accumulate at the gap edge (see Figures 2 and 4), increasing excitation torque density in this region and reducing it outside of it, at |x| ≳ 10h. This effect has not been taken into account in previous studies, but it does drastically change the overall shape of the torque density.

We tried to improve the performance of the LR theory (Section 4.1) by assigning torque contributions of individual modes not to x_L = 2/(3k_y) (which is appropriate only for a uniform disk), but to x_L(k_y) given by the self-consistent calculation in a non-uniform disk; see Equation (24). Unfortunately, the agreement with our fully self-consistent calculation has gotten even worse (for this reason we did not pursue this idea further), which suggests that properly accounting for the resonance overlap is also very important for getting the correct spatial distribution of the torque.

Far enough from the gap edge, where Σ₀(x) flattens out, the excitation torque density becomes negative, which is a direct manifestation of the negative torque density phenomenon (Dong et al. 2011b). For the adopted gap profile this happens at x ≈ 15h, and dT/dx stays negative out to infinity and agrees with the asymptotic behavior −0.63(h/x)⁴ derived in RP12.

Additionally, in the bottom panel of Figure 6 (dot-dashed curve) we plot dT/dx calculated according to the simplified prescription (3) but using the correct excitation torque density for the uniform disk dT/dx|_u calculated in RP12 instead of the GT80 prescription (4). One can see that this calculation fails miserably compared to the correct dT/dx: it predicts negative excitation torque density for x > 3.2h, which is the direct consequence of the negative torque density phenomenon in the uniform disks. Because of this, as we will later see in Figure 7, the use of this prescription results in the negative integrated torque deposited in the density wave by the perturber, which cannot be true in a real disk.

Zoom In Zoom Out Reset image size

The point of this last exercise was to illustrate the failure of the simple prescription (3), which yields unphysical results even when it uses physically motivated ingredients—the spatial behavior of dT/dx|_u derived in RP12 as opposed to the asymptotic formula (4).

In Figure 7, we plot the integrated torque T(x) for our fiducial model and the AMF F_H(x) calculated using Equation (29). The two have small relative vertical offset but are otherwise identical, as expected from the discussion in Section 4. We also display T^LR(x) obtained by integrating dT^LR/dx given by Equation (32) from zero to x.

All these expressions almost coincide within the gap. This happens because in the flat-bottomed gap different methods of computing integrated torque should yield the same result: a total torque (or AMF) of ≈0.93 (GT80) in units of the plot times Σ_min = 3 × 10⁻⁴. This equivalence of the full torque calculations in the Fourier space (T^LR(x)) and physical space (T(x)) has been previously mentioned in RP12.

At the gap edge, for |x| ≳ 5h, one finds that T(x) increases faster than T^LR(x), which is a direct consequence of the torque density dT/dx in a self-consistent non-uniform disk calculation being rather concentrated toward the perturber; see Section 6.1. Quite remarkably, despite this difference at the gap edge, the accumulated torque at infinity nearly coincides with that given by the simple theoretical prescription T^LR(x): our self-consistent calculation yields full torque ≈1.2 × 10⁻³, while the LR theory predicts ≈1.3 × 10⁻³ (see Figure 7 at x = 30h for reference). This means that the total torque exerted on the disk (but not the spatial distribution of the torque density!) is well described by the simple theoretical prescription (essentially T^LR(x) based on asymptotic scaling (4)) that has been broadly used in the literature.

We do not have a fully satisfactory explanation for this coincidence except for the following observation. In the lower panel of Figure 5 we display the torque in Fourier space (on linear scale) for k_yh < 0.4. There it can be seen that our non-uniform disk calculation generates almost the same amount of torque as the theoretical prescription (integrals under the curves are roughly the same), meaning that the torque deficit with respect to the LR theory visible for k_yh ≲ 0.1 is almost exactly compensated by the torque excess for k_yh ≳ 0.1. This means that the effect of resonances, which are being shifted to lower density regions resulting in lower torque contribution (smaller k_yh), is closely counterbalanced by the torque enhancement due to the resonances displaced to higher density region (larger k_yh). This is an interesting result, which only weakly depends on the exact gap profile as we empirically show next.

As discussed in Section 2.2 and shown in Figure 1, the parameter n controls the steepness of the density gradients at the gap edge. In panels (a) and (d) of Figure 8 we plot the excitation torque density and the integrated torque, respectively, for n = 2, 4, 6, while fixing δ = 5h and Σ_min = 3 × 10⁻⁴.

Zoom In Zoom Out Reset image size

As n is increased, the gradients become larger and, therefore, the amplitude of the variation of κ² around the gap edge becomes larger. This leads to higher resonance accumulation closer to the perturber, and Σ₀ gets higher there as well. As a result, the torque density gets more concentrated toward the gap center with increasing n (see Figure 8(a)): for n = 6 almost all the torque is excited within the range 5h ≲ x ≲ 10h, while for n = 2 this region extends out to x ∼ 15h.

The integrated torque also increases with n since closer to the gap center excitation by the perturber is more effective, depositing more angular momentum into the density waves. We see from Figure 8(d) that the density profile with n = 2 results in a total torque of ∼1.2 × 10⁻³ in units of (GM_p)²Σ_∞/hc²_s, while using n = 6 we obtain ∼2.1 × 10⁻³.

This figure also demonstrates that for all values of n the LR theory predicts more extended torque density distribution than we find in self-consistent calculations; see Section 6.1. Nevertheless, the integrated torque at |x| → ∞ is still well predicted by the LR theory, the point we previously made in Section 6.2. The largest relative difference of ∼13% is found for n = 6, for which we obtain the integrated torque of ∼2.1 × 10⁻³, while the LR prescription gives ∼2.4 × 10⁻³.

In Figures 8(b) and (e) we vary the width of the gap by slightly changing δ, but fixing Σ_min = 3 × 10⁻⁴ and n = 2. As expected, for narrower gaps the torque density peaks closer to the planet and the integrated torque reaches higher values, simply because the strength of the tidal coupling is higher closer to the perturber.

Figure 8(e) shows that the integrated torque predicted by the LR theory tends to better match the self-consistent calculation as the gap width is increased. Also, the dependence of the T(x → ∞) on δ is not linear: decreasing δ from 6h to 5.5h results in ∼14% increase of the integrated torque, while going from 5.5h to 5h results in the increase of ∼25%. Using the fact that the LR theory gives a reasonably good approximation to the integrated torque behavior, and adopting the prescription (3) with dT/dr|_u given by Equation (4), one can estimate that the total torque excited outside the density gap scales as T(x > δ)∝(δ/h)³. We verified this simple scaling by calculating integrated torque for density profiles with δ = (4–8)h. This behavior is in agreement with Lin & Papaloizou (1984), who present the same scaling for the one-sided torque.

In Figures 8(c) and (f) we vary the dimensionless gap depth Σ_min from 10⁻⁴ to 10⁻³, while fixing n = 2 and δ = 5h. One can see that changing Σ_min only affects the torque density and the full torque accumulated inside the gap (in a linear fashion); outside the gap (for x > 4h) torque structure remains unchanged.

As Σ_min is increased, there is a point at which the integrated torque starts being dominated by the excitation inside the gap, in the immediate vicinity of the perturber. One can easily predict when this happens. For a flat-bottomed gap profile like the one we consider in this work the torque accumulated inside the gap is T(x < δ) ∼ 0.9Σ_minΣ_∞G²M²_p/hc²_s (RP12 and GT80) if the gap width δ is larger than about 2h—the extent or the region where most of the torque excitation occurs in a uniform disk. This exceeds the torque T(x > δ) excited outside the gap for

$$\begin{eqnarray} &&\Sigma _{\rm min}\gtrsim \frac{T(x>\delta)hc_s^2}{0.9\Sigma _\infty G^2M_p^2}\approx 2.8h^3 \int _\delta ^\infty \frac{\Sigma _0(x)}{x^4}dx, \end{eqnarray} \tag{ 34 }$$

since the value of T(x > δ) is well approximated by the LR theory using the prescription (3) with dT/dr|_u given by Equation (4). In the case shown in Figure 8(c) and (f) the critical value of Σ_min is ≈10⁻², i.e., the integrated torque in this case is dominated by the excitation at the edge of the disk only when the density inside the gap is depleted below 1% of Σ_∞.

Finally, after trying different combinations of parameters for the density profile, we observe that the torque density always decays exponentially outside the gap in regions where the density gradients are important, as discussed in Section 6.1. In Appendix B we construct an empirical fit for the torque density in this region represented by Equation (B1).

We noticed in Section 3.2 that the conservation of the angular momentum forces the relative amplitude of the density wake Σ₁/Σ₀ to decrease as the wave climbs up the density gradient at the gap edge. Interestingly, even though the normalized Σ₁/Σ₀ does not vary much as the wave travels through the gap edge, the angular momentum density varies a lot; this is because the latter is sensitive to the overall wake shape and not just its amplitude.

The variation of the relative density perturbation has interesting implications for the strength of the density waves outside the gap. In particular, from the fact that the normalized wake amplitude in Figure 3 does not vary much with x it follows that the wave nonlinearity measured by the ratio Σ₁/Σ₀∝(x/Σ₀(x))^1/2 should decrease by a factor of ∼[(δ/h)/Σ_min]^1/2 between the bottom of the gap and the part of the disk just outside of the gap, where Σ₀ ∼ Σ_∞. For clean gaps with small Σ_min this factor can be quite significant, ≈0.2 for a gap with the width of δ = 5h and density contrast at the bottom Σ_min = 10⁻² (this factor is ∼0.05 for the situation shown in Figure 3(c)). Thus, even if the density wave starts out highly nonlinear (Σ₁/Σ₀ ≳ 1) inside the gap near the planet, it may become linear upon propagating out of the gap. The initial wave nonlinearity in the flat part of the density profile near the planet is Σ₁/Σ₀ ∼ M_p/M_th, where M_th ≡ c³_s/(ΩG) is the characteristic planetary mass at which the Hill radius becomes comparable to the disk scale height. This mass is about 10 M_⊕ at 1 AU and about 60 M_⊕ at 10 AU. Thus, a Jupiter mass planet at 10 AU would generate a highly nonlinear wave with Σ₁/Σ₀ ∼ 5 at the bottom of the gap, but the nonlinearity should drop to Σ₁/Σ₀ ∼ 1 outside the gap for δ = 5h and Σ_min = 10⁻², even if we neglect the wave damping (not captured by our linear theory). This demonstrates that even rather massive planets capable of opening gaps may still have their density waves in the linear regime outside the gap, depending mainly on the density contrast Σ_min.

This point is rather important for the possibility of the nonlinear wave damping. Goodman & Rafikov (2001) have studied the nonlinear density wave evolution in a uniform disk and demonstrated that shock formation resulting from nonlinear effects can lead to effective damping of the wave. This study has been extended by Rafikov (2002a) to include the possibility of radial variation of Σ₀. Our linear results suggest, in agreement with the latter work, that the nonlinear wake distortion (which ultimately results in shock formation) should slow down as the wave gets out of the gap, shock formation gets postponed, and the nonlinear dissipation becomes less effective at transferring wave angular momentum to the disk material.

For massive planets the wave might shock close to the planet, while it is still propagating through the roughly uniform bottom part of the gap profile: Goodman & Rafikov (2001) estimate the shocking length of just 2h for a 2M_⊕ planet at 1 AU, and it is even smaller for more massive planets. Then the efficient density wave dissipation will start inside the gap but then will appreciably slow down as the wave climbs out of the gap (even though the shock will still persist). This complication should certainly affect the picture of the nonlinear wave damping and may produce important feedback for the self-consistent calculations of the gap density profile.

Theory of disk–satellite interaction has been extensively used for understanding the interaction of protoplanetary disks with embedded planets and the evolution of SMBH binaries surrounded by circumbinary accretion disks, among other astrophysical settings. Our extensions of this theory have direct impact on our understanding of these systems.

In particular, our finding of the torque density being quite concentrated near the perturber suggests that regardless of the details of the wave damping, the angular momentum carried by the density waves should be deposited in the disk closer to the perturber than the conventional theory would suggest. This certainly facilitates the formation of density gaps or cavities since for the same mass of the perturber a narrower annulus of the disk needs to be cleaned by the tide. This may make Type II migration possible for lower mass planets than was previously thought.

To fully understand the impact of our results on the Type II migration and SMBH inspiral, one needs to know a self-consistent gap shape, which can be computed only when the wave damping mechanism is specified. In this regard we note that the majority of studies of the gap or cavity opening simply ignore the issue of wave damping process and take dT/dr|_d = dT/dr (Lin & Papaloizou 1986; Chang et al. 2010), i.e., assume immediate deposition of the excitation torque into the disk. However, the validity of this approximation has never been demonstrated, and calculations that do properly take wave damping into account (Ward 1997; Rafikov 2002b) result in a rather different picture of the gap opening. Our results on the wave amplitude evolution and its implications for the nonlinear wave damping in Sections 3.2 and 8.1 should be useful for understanding the spatial distribution of dT/dr|_d and calculating the gap shape (and orbital evolution of the perturber) fully self-consistently.

Our present work is cast in terms of the shearing sheet approximation, which limits the direct application of its results to systems in which gaps are narrow compared to the semimajor axis of the perturber. This is the case for, e.g., planets in protoplanetary disks (Lin & Papaloizou 1986; Takeuchi et al. 1996) and SMBH binaries with extreme mass ratio surrounded by an accretion disk (so-called extreme mass-ratio inspirals; Armitage & Natarajan 2002; Chang et al. 2010; Kocsis et al. 2011).

SMBH binaries of similar mass do not fall in this category since they are known to clear out a large and clean cavity in circumbinary gaseous disk, with inner radius comparable to the binary separation (e.g., MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Shi et al. 2012; Farris et al. 2011). However, even though the shearing sheet approximation breaks down in these systems, our results can still be used to gain qualitative insight into the dynamics of cavity clearing. In particular, concentration of the torque density toward the perturber means that the cavity should be opened by lower mass ratio SMBHs than was thought before.

Also, our results can help understand the excitation torque density patterns derived from simulations, most of which (MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Shi et al. 2012; Farris et al. 2011; Roedig et al. 2012) agree that outside the cavity dT/dr exhibits a complicated oscillatory pattern (some of these studies include additional physics such as MHD or general relativity). This behavior was not recognized before, and some earlier studies (e.g., Armitage & Natarajan 2002) tried modeling numerically derived dT/dr in the form similar to GT80's result (4) but with the amplitude reduced by a factor of ∼10⁻². Such low amplitude most likely results from averaging the actual spatially rapidly oscillating pattern of dT/dr and then applying the ∝|r − r_p|⁻⁴ scaling to reproduce the total integrated torque in the simulations. Clearly, usage of this prescription is likely to lead to misleading results (Liu & Shapiro 2010; Chang et al. 2010) in modeling gap shape.

MacFadyen & Milosavljević (2008) suggested that the oscillatory behavior of the torque distribution in their simulations is consistent with forcing at the 2:3 (m = 2) outer LR and tried fitting linear waveforms from Meyer-Vernet & Sicardy (1987) to spatial variation of the perturbed fluid variables. Even though they were able to reproduce qualitatively the oscillatory behavior of the numerical dT/dr, the amplitude and phase of the torque density were substantially different between the numerical calculation and the analytical waveform, qualitatively consistent with the LR shift and reduction of perturbation amplitude in our calculations; see Figure 4(f). This is certainly not surprising since not only were the Meyer-Vernet & Sicardy (1987) waveforms derived for a uniform disk, but they are also valid only in the close vicinity of the LR. The latter limitation was removed in RP12, who came up with a fully global analytical approximation for the waveforms, but again only for a uniform disk case. Our present work is free from these constraints since it provides a way of computing $(dT/dr)_{k_y}$ for arbitrary disk surface density profile at any distance from the resonance; see Section 5.1. Future extensions of our approach to full cylindrical geometry will be directly applicable to the astrophysical systems such as SMBH binaries in which broad gaps or cavities are typical.

Finally, another suitable application of our results is the Type I migration of planets near sharp density transitions, where the linear theory is directly applicable (Masset et al. 2006; Hasegawa & Pudritz 2011; Masset 2011). In particular, we confirm that the commonly used LR prescription for the integrated one-sided gravitational torque (a quantity that determines the migration speed) coincides with our results within 10% for all the density gap profiles we studied and may be accurate enough to determine the migration speed of small planets.