PLANET–DISK INTERACTION IN THREE DIMENSIONS: THE IMPORTANCE OF BUOYANCY WAVES

In our calculations, we use both an isothermal (γ = 1 in p∝ρ^γ), and an adiabatic EOS (γ = 5/3 or 7/5). The unit of length in our simulations is defined via the adiabatic scale height $h=c_s/\Omega _p=\sqrt{\gamma }c_{{\rm iso}}/\Omega$, where c_iso ≡ (kT/μ)^1/2 and c_s is the adiabatic sound speed calculated using the midplane temperature. Obviously, c_s = c_iso for an isothermal EOS (γ = 1). This choice keeps sound waves traveling the same distance in a given amount of time in different models so that the wake position is always the same.

The vertical structure of a protoplanetary disk is typically determined by irradiation from the central star so that it is vertically isothermal:

$$\begin{eqnarray} &&\frac{\rho _{0}(z)}{\rho _{00}}={\rm exp} \left(-\frac{z^{2}\Omega ^{2}}{2c_{{\rm iso}}^{2}}\right), \end{eqnarray} \tag{ 4 }$$

where ρ₀₀ is the midplane density. The value of c_iso is determined by the central irradiation flux (Kenyon & Hartmann 1987; Chiang & Goldreich 1997).

We also consider thermally stratified disks in which the disk midplane is hotter than its atmosphere, which is expected if viscous heating is important. Including the fact that such disks will eventually become isothermal at height z, the vertical structure of the disk is better represented as

$$\begin{eqnarray} &&p_{0}(z)=c_{s,s}^{2}\rho _{0}(z)\left[\left(\frac{\Gamma -1}{\Gamma }\right)\left(\frac{\rho _{0}(z)}{\rho _{s}} \right)^{\Gamma -1}+1\right], \end{eqnarray} \tag{ 5 }$$

where c_{s, s} and ρ_s(set as 0.01 ρ_c) represent the sound speed and density at the transition where the disk becomes isothermal (Lin et al. 1990; Bate et al. 2002). However, unlike Bate et al. (2002), we do not necessarily set Γ in Equation (5) the same as γ in the EOS. We derive ρ₀(z) by solving the equation of vertical hydrostatic equilibrium with Equation (5) using the Newton–Raphson method. When Γ → 1 the isothermal disk structure is recovered.

We run a set of five models up to 10 orbits with different disk structure and EOS, summarized below and in Table 1.

• 1.  
  Isothermal disk structure (4) with isothermal EOS (model II).
• 2.  
  Disk structure (5) with adiabatic EOS and Γ = γ = 5/3 (model AA).
• 3.  
  Isothermal disk structure with adiabatic EOS having γ = 5/3 (model IA). This setup is typical for a centrally irradiated disk in which density perturbations are adiabatic.
• 4.  
  Analogous to (4) but with γ = 7/5 (model IA2), designed to illustrate the effect of adiabatic index on the buoyancy wave coupling.
• 5.  
  Disk structure (5) with Γ = 5/3 and adiabatic EOS with γ = 7/5 (model PA), to illustrate viscously heated accretion disks.

In the first two models N = 0 so that buoyancy waves are absent, while the last three models feature non-zero N and support buoyancy waves. All models were run at a resolution of 32/h, but we also run model (3) at higher resolution (64/h) to verify numerical convergence.

The planet mass is M_p = 0.0058 M_th, where M_th is the thermal mass

$$\begin{equation} M_{{\rm th}}\equiv \frac{c_{s}^{3}}{G\Omega _{p}}. \end{equation} \tag{ 6 }$$

Thus, disk–planet coupling is well inside the linear regime. Density wave generated by such a planet in a 2D disk would shock at |x| = 6 h (Goodman & Rafikov 2001), which is outside our box.

First we show the results for models with no buoyancy waves (II and AA). The spatial distributions of the excitation torque density dT/dx (the amount of torque excited per unit radial distance) for these two models are shown in Figures 2(a) and (b) as green and orange curves. One can see that dT/dx rapidly goes to zero for |x| ≲ h, clearly exhibiting the "torque cutoff" phenomenon (Goldreich & Tremaine 1980) in both models.

Zoom In Zoom Out Reset image size

For model II, our simulation agrees with the semi-analytical calculation by Takeuchi & Miyama (1998) remarkably well, and we confirm that Equation (1) holds in this case with q_3D(II) = 0.25q_2D. Although wave structures are different in model AA (Lubow & Pringle 1993), the distribution of dT/dx for this model is quite similar in shape to that of model II, but the amplitude is slightly different.

We next look at the behavior of dT/dx in models with non-zero N. The most important feature of these models, evident in Figure 2, is the weak torque cutoff near the planet—all show a significant torque contribution near x = 0, including model IA2, which provides the best description for an irradiated protoplanetary disk. Since higher-m modes are excited at Lindblad resonances closer to the planet, the strong torque close to the planet also suggests the traditional torque cutoff in Fourier space (Ward 1997) is weaker when buoyancy waves are present.

The one-sided torque T (Figures 2(c) and (d)) in models IA and IA2 is characterized by q_3D(IA) = 0.48 and q_3D(IA2) = 0.41, respectively, which should be compared with q_3D(II) = 0.34 for model II having the same vertical structure as IA and IA2 but different EOS.

In order to identify the origin of this excess torque near the planet, we integrate the volume density of the torque along the y direction

$$\begin{equation} \frac{d^{2}T}{dxdz}=\int _{-\infty }^{\infty }dy\delta \rho (x,y,z)\frac{\partial \phi }{\partial y}. \end{equation} \tag{ 7 }$$

This quantity is shown in Figure 3(a) for models II and IA. We see that the excess torque in model IA comes from the region z ∼ 0.4 h and |x| < h. In Figure 3(d) we plot the density contours in the xy plane at z = 0.4 h. In clear contrast with model II (Figure 3(c)), the density distribution in model IA exhibits density fluctuations/ridges close to the planet which are different from the usual wake structure. They extend to y > 0 and look like rays emanating from the origin. It is these density fluctuations that contribute to the excess torque.

Zoom In Zoom Out Reset image size

These fluctuations also have vertical structure, as shown in Figure 4 where we plot δρ = ρ − ρ₀ and v_z in the yz plane at x = 0.5 h for both models II and IA. Again, significant density fluctuations close to the planet exist only in case IA.

Zoom In Zoom Out Reset image size

Since the only difference between case IA and case II is that the case IA has non-zero Brunt–Väisälä frequency N, it is natural to relate these density fluctuations to buoyancy waves. To confirm this hypothesis, we note that a mode with wavenumber k_y has the buoyancy frequency N whenever the condition

$$\begin{equation} 2Ak_{y}x=N \end{equation} \tag{ 8 }$$

is fulfilled. For comparison, the Lindblad resonance condition is 2Ak_yx = ±κ. In the disk with isothermal structure and adiabatic EOS (models IA and IA2), we have

$$\begin{equation} N=\frac{\Omega z}{h}\sqrt{\frac{\gamma -1}{\gamma }}. \end{equation} \tag{ 9 }$$

With γ = 5/3 and λ_y = 2π/k_y, Equations (8) and (9) can be combined to derive

$$\begin{equation} \frac{\lambda _{y}}{h}=11.543 \frac{x}{z}. \end{equation} \tag{ 10 }$$

If we assume the phase of buoyancy waves is 0 at x = 0, the geometric location of the constant phase 2nπ (n is integer) is given by

$$\begin{equation} \frac{y}{h}=11.543 \frac{x}{z}n\quad {\rm and}\; n=1,2,3,\ldots \,. \end{equation} \tag{ 11 }$$

Curves corresponding to Equation (11) are drawn in the lower right panels of Figures 3 and 4. They agree with the pattern of the density perturbation and v_z quite well, strongly suggesting that these fluctuations and the excess torque are related to buoyancy waves (their dissipation or excitation).

The torque density and integrated torque for model IA2 (isothermal disk with adiabatic EOS and γ = 1.4) plotted as the cyan curves in Figure 2 also exhibit excess torque density due to buoyancy waves close to the planet. However, the integrated torque in this case is characterized by q_3D(IA2) ≈ 0.41, which is lower than q_3D(IA), and the buoyancy wave coupling to the planetary potential is weaker for smaller γ as seen from dT/dx curve.

Blue curves in the right panels of Figure 2 show dT/dx and T(x) for model PA. There is again an apparent torque excess near the planet, which is not surprising since N ≠ 0 in this model, allowing buoyancy waves to be excited.

Previous analytical and semi-analytical studies of buoyancy waves (Lubow & Pringle 1993; Korycansky & Pringle 1995; Lubow & Ogilvie 1998) have been limited to axisymmetric waves or low-m modes excited far away from the planet, at |x| ≫ h. These studies have demonstrated such modes to play rather insignificant role in planet–disk angular momentum exchange.

The main result of this work is that coupling of the planetary potential to the locally excited (|x| ≲ h) buoyancy waves channels a considerable amount of angular momentum into these modes. As a result, the buoyancy torque strongly contributes to the planetary torque, at the level of tens of percent. This result is different from analytical expectations extrapolated from low-m buoyancy modes (Lubow & Ogilvie 1998).

Most previous 3D simulations of planet–disk interaction explored isothermal disk structure with an isothermal EOS (e.g., Bate et al. 2003; D'Angelo & Lubow 2010). An adiabatic EOS has been explored in Paardekooper & Mellema (2006). However, because of the global geometry used in their work, the buoyancy torque may have been hard to distinguish from the corotation torque; see Section 4.3.

The vertical structure of protoplanetary disks around T Tauri stars should be close to isothermal given the dominant role of stellar irradiation in their thermal balance. Existence of buoyancy waves in such disks depends on the thermodynamics of the fluid perturbations on timescales ∼N⁻¹. Density perturbations induced by planets result in temperature fluctuations which tend to be erased by radiative diffusion. If the cooling time t_c of such perturbations is longer than N⁻¹ then modes are adiabatic with γ > 1 (as in models IA and IA2) and buoyancy waves should be present. In the opposite case of t_cN ≲ 1 fluid perturbations behave as isothermal (γ ≈ 1) and buoyancy waves are not excited, as in model II.

We expect that the separation between the two regimes should occur between several to several tens of AU, depending on the disk properties, and that buoyancy waves should be efficiently excited by relatively close-in planets. Further out, t_cN ≲ 1 and excitation of buoyancy waves should be suppressed. Note that both t_c and N are in general functions of z; see Equation (9). This makes cooling considerations dependent not only on the radius but also on the height in the disk.

Whenever buoyancy waves are efficiently excited by planets they can reveal themselves via associated vertical motions (see Figure 1) which should give rise to corrugations of the disk surface aligned with density ridges seen in Figure 3. Near-IR imaging of stellar light scattered by such perturbations at very high resolution (on scales ∼h) can reveal even low-mass planets embedded in protoplanetary disks.

The rate at which a planet migrates is determined by the imbalance of one-sided torques exerted on it by the inner and outer parts of the disk, caused by the global radial gradients of the disk density and temperature (Ward 1986; Tanaka et al. 2002). The net torque can also be caused by the local asymmetries in the disk properties caused by the presence of the planet itself, such as the temperature perturbations due to shadowing near the planet (Jang-Condell & Sasselov 2005), or the planet's own heat output caused by, e.g., the accretion of planetesimals. Such asymmetries generated very close to the planet may not strongly affect the net standard Lindblad torque because of the torque cutoff at such separations. However, their effect on the net buoyancy torque excited primarily at small x may be disproportionately large and affect the planetary migration considerably.

This statement is only true if the buoyancy waves observed in this work can propagate and deposit their angular momentum far from the corotation region. Otherwise, their angular momentum accumulates in the narrow annulus around the planetary orbit and the corresponding torque on the planet may be subject to saturation, like the standard corotation torque (Balmforth & Korycansky 2001; Masset 2001, 2002). This would eliminate the effect of the buoyancy waves on the planetary migration. We will investigate these possibilities in Z. Zhu et al. (in preparation).

Direct measurement of the buoyancy torque effect on the migration speed would require global disk–planet calculations, in which corotation torque appears as well. It may then be difficult to separate the effects of the buoyancy and corotation torques in global simulations. However, the strength of the corotation torque depends on the radial gradients of specific vortensity and entropy across the horseshoe region (Paardekooper & Mellema 2006; Baruteau & Masset 2008; Paardekooper & Papaloizou 2008). By properly setting the disk properties to nullify these gradients one can eliminate the corotation torque in global simulations, thus isolating the contribution due to the buoyancy torque.

Gap opening by massive planets depends not only on the wave excitation but also on wave damping as a means of transferring angular momentum to the disk fluid (Lunine & Stevenson 1982; Rafikov 2002). The global density wake excited by planet is thought to dissipate primarily via shock damping (Goodman & Rafikov 2001; Dong et al. 2011), but buoyancy waves are very distinct from rotation-modified sound waves and should dissipate differently (Lubow & Ogilvie 1998; Bate et al. 2002). Channeling of the wave action (Lubow & Ogilvie 1998) may considerably speed up their nonlinear evolution resulting in more efficient damping. Since we showed that buoyancy waves carry good fraction of the total angular momentum flux, understanding their damping mechanism and spatial pattern of dissipation may be important for clarifying the issue of gap opening by planets (Zhu et al. in prep).

Our work suggests buoyancy waves can play an important role in planet migration and gap opening which demands further studies.