DENSITY WAVES EXCITED BY LOW-MASS PLANETS IN PROTOPLANETARY DISKS. I. LINEAR REGIME

The discovery of a number of "hot Jupiters"—extrasolar giant planets residing on very tight orbits—calls for understanding the dynamical pathways that brought these objects to their current orbits. Due to the difficulty of forming these planets in situ, it was suggested that potential hot Jupiters form at larger separations from their stars and are subsequently transported inward by some migration mechanism (Papaloizou & Terquem 2006; Papaloizou et al. 2007). In protoplanetary disks, planets act as sources of nonaxisymmetric density waves (Goldreich & Tremaine 1979, 1980, hereafter GT80) and coupling of planetary gravity to these waves acts to transfer angular momentum between the disk and the planet. This ultimately results in the orbital migration of embedded planets. At the same time, density waves generated by planets can affect the structure of the disk (potentially leading to gap opening) once the angular momentum that they carry is transferred to the disk material. Thus, better understanding of the density wave properties and processes of their excitation, propagation, and damping plays a crucial role in the global picture of planet formation.

Evolution of the density wave can be separated into the linear and nonlinear stages. The former applies to density waves excited by low-mass planets, which have not had enough time to propagate very far from their launching site (a more precise condition is formulated in Section 2). The linear theory of disk–planet interaction was pioneered by Goldreich & Tremaine (1979, GT80) who provided a description of the density wave excitation in two-dimensional (2D) gaseous disks. This theory has been subsequently refined in various ways, e.g., by calculating the asymmetry of planetary interaction with inner and outer portions of the disk (Ward 1986), and extending it to fully three-dimensional disks (Tanaka et al. 2002).

These studies concentrated on understanding the behavior of perturbed quantities Fourier decomposed in azimuthal angle (Artymowicz 1993a, 1993b; Korycansky & Pollack 1993), often addressing only the integral wave properties such as the total angular momentum carried by the waves, or the differential Lindblad torque (Ward 1986). Only recently has the linear evolution of the wave properties in physical space, namely the spatial dependence of the perturbed fluid velocity and surface density, been addressed by Goodman & Rafikov (2001, hereafter GR01) in the local shearing sheet approximation. Subsequently, Ogilvie & Lubow (2002) and Rafikov (2002a) studied the shape of the density wake in the global setting. This paved way for subsequent exploration of the nonlinear wave evolution (GR01; Bate et al. 2002), understanding of which requires knowledge of the wave characteristics in physical space.

In parallel with these analytical and semi-analytical studies, numerous attempts to verify linear theory by direct hydrodynamical simulations were made. However, the majority of these studies compared with theoretical predictions only the secondary or derived properties of the planet–disk system related to the density wave excitation and evolution. For example, in their simulations Lin & Papaloizou (1993), Ward (1997), D'Angelo et al. (2003), and D'Angelo & Lubow (2008) measure the rate of planet migration due to its gravitational interaction with the disk—a quantity that is determined by the small difference between the summed (over all azimuthal harmonics) one-sided torques that the planet exerts on the inner and outer portions of the disk. Other studies (Lin & Papaloizou 1993; Muto et al. 2010) tried to measure the planetary mass at which a gap opens in the disk—a process that depends not only on the magnitude of the angular momentum carried by the waves but also on the intricate details of their dissipation. In both situations the verification of the linear theory ends up being rather indirect.

It is certainly much more convincing and more robust to verify linear theory by comparing with its predictions the primary wave properties, such as the spatial distribution of the perturbed fluid properties derived from simulations. Recently, some effort in this direction has been made. In particular, D'Angelo & Lubow (2008) investigated the distribution of the torque imparted by the planet on the disk as a function of radial distance. Muto et al. (2010) showed the spatial structure of the density perturbation associated with the wave obtained in their 2D simulations. Unfortunately, in both cases no quantitative comparison with linear theory predictions was provided.

In this work, we use 2D local inviscid hydrodynamical simulations in the shearing sheet geometry to verify the linear theory. For the first time, we provide direct comparison of the primary physical variables—perturbed surface density profile, torque density in both physical and Fourier spaces, and angular momentum flux (AMF)—derived from simulations with the theoretical predictions. Our verification effort is the most thorough, sophisticated, and quantitative to date, and we achieve agreement with theory at the level of several percent for all physical variables explored. In addition, using our numerical results we address some remaining issues in the linear theory, such as the different expressions for the torque density in Fourier space (GR80; Artymowicz 1993a). We also formulate a set of requirements which should be satisfied by simulations to reach convergence and to reliably reproduce the linear theory results. The way in which we make the comparison and the agreement we achieve can serve as a standard for checking the performance of different codes.

The paper is structured as follows. We briefly summarize the main results of the linear density wave theory in Section 2. A description of the code and our numerical setup is given in Section 3, including a discussion of potential problems associated with implementing orbital advection algorithm in studies of the disk–satellite interaction. We present our numerical results and compare them with theoretical predictions in Section 4, followed by the investigation of the sensitivity of simulation outcomes to the different numerical algorithms in Section 5. In Section 6, we describe and discuss a commonly ignored numerical problem in disk simulations which if not handled appropriately can lead to incorrect results. Finally, a summary is provided in Section 7.

In the linear theory of disk–planet interaction laid out in cylindrical geometry (Goldreich & Tremaine 1979, GT80), the planetary potential is normally decomposed into a series of Fourier harmonics. Assuming all perturbed fluid variables to be proportional to exp [i(m(φ − Ω_pt) + ∫^rk_r(r')dr')] (m is a positive integer) and neglecting the self-gravity of the disk, one can obtain the dispersion relation for the mth wave harmonic in the WKB limit in the form

$$\begin{eqnarray} &&m^2[\Omega (r)-\Omega _p]^2=\kappa ^2+k_r^2c_s^2, \end{eqnarray} \tag{ 1 }$$

where κ is the epicyclic frequency (κ = Ω for a Newtonian potential), k_r is the wavenumber, Ω is the fluid angular frequency, and Ω_p is the pattern speed of the perturbation equal to the angular frequency of the point-mass perturber moving on a circular orbit with semimajor axis r_p, that is, Ω_p = Ω(r_p). As Equation (1) demonstrates, there are two (inner and outer) Lindblad resonances (LR) where k_r → 0 and the mth harmonic of planetary potential couples best to the wave-like perturbation. These resonances are located at radii r_{L, m} given by Ω(r_{L, m}) = Ω_pm/(m ± 1), which in the limit of m ≫ 1 results in

$$\begin{equation} x_{L,m}=\pm \frac{2h}{3\mu }, \end{equation} \tag{ 2 }$$

where x ≡ r − r_p, h = c_s/Ω_p is the vertical disk scale length, and μ ≡ mh/r_p.

The superposition of all the harmonics with different azimuthal wavenumbers m forms inner and outer spiral density waves. In physical space, fluid perturbation associated with the wave is concentrated in narrow wake (with azimuthal extent of order h), the location of which is described in the local limit by (GR01)

$$\begin{eqnarray} &&y_w \approx -\mbox{sgn}(x)\frac{3}{4}\frac{x^2}{h}, \end{eqnarray} \tag{ 3 }$$

where y ≡ r_p(φ − φ_p) and azimuthal angle φ is counted in the prograde sense (with respect to planetary orbital motion). This location corresponds to the stationary phase of different Fourier harmonics composing the wave, as can be easily seen from Equation (1). The wake shape in the global setting accounting for cylindrical geometry of the disk and Keplerian profile of Ω(r) was explored by Ogilvie & Lubow (2002) and Rafikov (2002a).

Each of the perturbation harmonics carries AMF, which for the mth harmonic is given by (GT80)

$$\begin{eqnarray} && F_H(m)=F_H^{\rm WKB}(m)\frac{F_H(m)}{F_H^{\rm WKB}(m)}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} && F_H^{\rm WKB}(m)=\frac{4}{3}\frac{m^2\Sigma _0}{\Omega _p^2} \left(\frac{{\it G M}_{p}}{r_p}\right)^2\left[2K_0(2/3)+K_1(2/3)\right]^2.\quad\ \ \end{eqnarray} \tag{ 5 }$$

Here Σ₀ is the unperturbed surface density of the disk, F^WKB_H(m) is the asymptotic Fourier contribution of the flux reached beyond the launching region on one side of the disk calculated in the WKB approximation (K₀ and K₁ are the modified Bessel functions of zeroth and first order), and the expression is only accurate for m ≫ 1. The ratio F_H(m)/F^WKB_H(m) accounts for the difference between the precise value of F_H(m) and F^WKB_H(m). This difference becomes important for m ≳ r_p/h ≫ 1 since, as demonstrated by GT80, F_H(m) decays exponentially with m at high azimuthal wavenumbers—the so-called torque cutoff phenomenon. As a result, the total AMF carried by the wave,

$$\begin{eqnarray} &&F_H=\sum _{m=1}^\infty F_H(m) \approx 0.93\,({\it G M}_{p})^2\frac{\Sigma _0 r_p\Omega _p}{c_s^3}, \end{eqnarray} \tag{ 6 }$$

is dominated by harmonics with m ∼ r_p/h (or μ ∼ 1), which are excited within a radial separation of order h from the planet; see Equation (2). This local nature of the wave excitation has important implications for the possibility of analytical treatment of the wave evolution.

In real protoplanetary disks, the torques acting on the inner and outer disks are slightly different (fractional difference ∼h/r_p ≪ 1), giving rise to a net torque on the planet and leading to Type I planetary migration (Ward 1986). However, in shearing sheet configuration the inner and outer parts of the disk are identical so there is no net torque acting on the planet and migration does not arise.

There is an important planetary mass scale characterizing the excitation and propagation of the density waves, the so-called thermal mass (GR01), which we define² here as

$$\begin{eqnarray} &&M_{\rm th}=\frac{c_s^3}{G\Omega _{p}}\approx 12\, M_\oplus \left(\frac{c_s}{1 \mbox{ km s}^{-1}}\right)^{3} \left(\frac{M_\odot }{M_\star }\right)^{1/2} \left(\frac{r_p}{\rm AU}\right)^{3/2}.\nonumber\\ \end{eqnarray} \tag{ 7 }$$

There are several reasons for the significance of this mass scale. First, a Hill radius of the planet with M_p = M_th is of order the scale height of the gaseous disk h. Second, the Bondi radius R_B ≡ GM_p/c²_s, which is the size of a region in which planetary gravity strongly perturbs local pressure distribution, also becomes comparable to h for M_p = M_th (Rafikov 2006). Third, and most importantly for our present study, planets with M_p ≳ M_th generate density waves which are nonlinear from the very start (δΣ/Σ₀ ≳ 1).

In the small planet, mass limit of

$$\begin{eqnarray} &&M_{p}\lesssim M_{\rm th}, \end{eqnarray} \tag{ 8 }$$

perturbations induced by the planetary gravity are weak (δΣ/Σ ≲ 1) and excitation of the density wave by the planet can be studied purely in the linear approximation. It is in this limit that the linear theory of wave driving developed by GT80 and all of its results are valid. In the opposite limit of M_p ≳ M_th the nonlinear effects must be accounted for even at the stage of wave excitation, which is not feasible analytically.

As mentioned above, the generation of the wave is largely completed at a distance of order h away from the planet, because this is where the most important Lindblad resonances lie. After that, in the absence of viscosity, the wave travels freely, no longer affected by the planetary gravity. At this stage of propagation nonlinear effects start to accumulate and finally result in wave breaking, formation of a shock, and transfer of wave energy and angular momentum to the disk fluid. Prior to the steepening into a shock the AMF F_H carried by the wave far from the planet (outside of the excitation region |x| ∼ h) is exactly conserved in the absence of linear dissipation. This nonlinear stage has been explored analytically in GR01 in the limit (8). The numerical results on the nonlinear wave evolution will be presented in Dong et al. (2011, hereafter Paper II).

We do not cover this stage in this work, but we mention the result of GR01 for the radial separation away from the planetary orbit at which the shock induced by the nonlinear effects first occurs:

$$\begin{equation} l_{\rm sh}\approx 0.8 \left(\frac{\gamma +1}{12/5}\frac{ M_{p}}{M_{\rm th}}\right)^{-2/5}h, \end{equation} \tag{ 9 }$$

where γ is the adiabatic index of gas. For M_p ≪ M_th, one finds the shock to be well separated from the planet, l_sh ≳ h, so that the wave excitation occurs in the linear regime, as asserted before. However, even for M_p ≪ M_th the linear theory of wave propagation fails for x ≳ l_sh. In this work we are always exploring the low planetary mass limit (8) and consider only |x| ≲ l_sh so that linear theory applies.

Most analytical work on the linear evolution of density waves was carried out in application to 2D disks. To provide a meaningful comparison with these studies, in our current work we employ 2D hydrodynamical simulations in local shearing box configuration. The computational tool for this investigation is Athena, a grid-based code for astrophysical gas dynamics using higher-order Godunov methods. Athena is written in a conservative form, so it conserves mass, momentum, and energy down to machine precision. The mathematical foundations of the algorithms are described in Gardiner & Stone (2005, 2008), and a comprehensive description of the implementation and tests of the algorithms is given in Stone et al. (2008).

The implementation of the shearing box approximation in Athena is described in Stone & Gardiner (2010). This approximation adopts a frame of reference located at radius r_p corotating with the disk at orbital frequency Ω_p = Ω(r_p). In this frame, the 2D hydrodynamics equations are written in a Cartesian coordinate system (x, y, z) that has unit vectors ${\bf \hat{i}}$, ${\bf \hat{j}}$, and ${\bf \hat{k}}$ as (Stone & Gardiner 2010)

$$\begin{eqnarray} && \frac{\partial \Sigma }{\partial t} + {\bf \nabla\, \cdot\, } (\Sigma {\bf v}) = 0, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} && \frac{\partial \Sigma {\bf v}}{\partial t} + {\bf \nabla\, \cdot\, } (\Sigma {\bf vv} + p) = \Sigma \Omega _p^{2}(2qx{\bf \hat{i}}) - 2\Omega _p {\bf \hat{k}} \times \Sigma {\bf v} - \Sigma \nabla \Phi _p, \nonumber\\ \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} && \frac{\partial E}{\partial t} + \nabla \cdot (E{\bf v} + p \cdot {\bf v}) = \Omega _p^{2}\Sigma {\bf v} \cdot (2qx{\bf \hat{i}}) - \Sigma {\bf v}\cdot \nabla \Phi _p,\quad\ \ \end{eqnarray} \tag{ 12 }$$

where Φ_p is the gravitational potential of the planet (whose form we will discuss in the next section), p is the gas pressure, and E is the total energy density (sum of the internal and kinetic energy, but excluding the gravitational energy). The shear parameter q is defined as

$$\begin{equation} q = - \frac{d {\rm ln} \Omega }{d {\rm ln} r}, \end{equation} \tag{ 13 }$$

so that for Keplerian flow q = 3/2. The equation of state (EOS) we use in the simulations is the ideal gas law with

$$\begin{equation} E_{\rm in} = \frac{p}{\gamma -1}, \end{equation} \tag{ 14 }$$

where E_in is the internal energy and γ = 5/3. We also use an isothermal EOS, in which case p = Σc²_s, where c_s is the isothermal sound speed. We do not include explicit viscosity (i.e., there is no explicit linear dissipation in the system) and we do not account for the self-gravity of the disk (thus we concentrate on studying the low-mass disks).

To accurately resolve the flow in the vicinity of a gravitating mass without making the time step too short one usually resorts to softening the potential of the perturber. There are many ways of doing this but all of them come up with the form of the potential that converges to the Newtonian potential Φ_K = −GM_p/ρ at large separations $\rho =\sqrt{x^2+y^2}$, beyond the softening length r_s. At small separations, ρ ≪ r_s, smoothed potentials become finite allowing the time step of simulations to stay finite. To test the sensitivity of our results to the specific method of potential softening we tried three different forms of the softened planetary potential. The second-order potential

$$\begin{eqnarray} &&\Phi _{p}^{(2)}=-{\it G M}_{p}\frac{1}{\big(\rho ^2+r_{\rm s}^2\big)^{1/2}} \end{eqnarray} \tag{ 15 }$$

converges to Φ_K at ρ ≫ r_s as (r_s/ρ)² (which means the fractional error is O((r_s/ρ)²) as r_s/ρ → 0). This is a standard form of the potential used in a majority of numerical hydrodynamic studies. We also studied the fourth-order potential

$$\begin{eqnarray} &&\Phi _{p}^{(4)}=-{\it G M}_{p}\frac{\rho ^2+3r_{\rm s}^2/2}{\big(\rho ^2+r_{\rm s}^2\big)^{3/2}} \end{eqnarray} \tag{ 16 }$$

converging to the point-mass potential as (r_s/ρ)⁴ for ρ ≫ r_s, and the sixth-order potential

$$\begin{eqnarray} &&\Phi _{p}^{(6)}=-{\it G M}_{p}\frac{\rho ^4+5\rho ^2r_{\rm s}^2/2+15r_{\rm s}^4/8}{\big(\rho ^2+r_{\rm s}^2\big)^{5/2}} \end{eqnarray} \tag{ 17 }$$

converging to Φ_K at ρ ≫ r_s as (r_s/ρ)⁶. The difference in accuracy with which these potentials represent Φ_K at large ρ can be quantified by the distance from the planet at which a given potential deviates from Φ_K by 1%. This distance is 7r_s for Φ⁽²⁾_p, 2.3r_s for Φ⁽⁴⁾_p, and 1.5r_s for Φ⁽⁶⁾_p. Note that the first-order potential Φ⁽¹⁾_p = −GM_p/|ρ + r_s|, which we do not use here, gives only 10% accuracy relative to Φ_K even at ρ = 10r_s, and we strongly discourage its use. The effects of different potential prescriptions will be discussed in detail in Section 5.3.

To guarantee accuracy and convergence of our results we carry out an exhaustive exploration of various numerical parameters characterizing our simulations, described in Section 5. Unless noted otherwise, our results shown in Section 4 are obtained with the following numerical parameters: an isothermal EOS, the Roe solver with third-order reconstruction in characteristic variables, and the corner transport upwind unsplit integrator (Stone et al. 2008), a resolution of 256 grid points per scale length h (subsequently denoted 256/h for brevity), and fourth-order accurate potential Φ⁽⁴⁾_p given by Equation (16) with softening length r_s = h/32.

High order of accuracy, high resolution, and an accurate form of planetary potential used in our runs allow us to properly capture the details of wave evolution. High resolution ensures low levels of numerical viscosity and prevents the angular momentum accumulated by the wave from spurious dissipation (see Section 5 for more discussion). We run a series of test simulations with otherwise identical conditions but different explicit Navier–Stokes viscosity, and measure the wave properties. As the explicit viscosity decreases, the simulation results gradually converge to the one with zero explicit viscosity, which indicates that the numerical viscosity dominates the explicit viscosity. For a typical simulation with low M_p = 2.09 × 10⁻²M_th and an isothermal EOS, the effective Shakura–Sunyaev α-parameter (α = ν/hc_s) characterizing our numerical viscosity is found to be below 10⁻⁵. Such small levels of viscosity are expected in dead zones of protoplanetary disks (Gammie 1996), where magnetorotational instability (MRI) may not operate effectively (Fleming & Stone 2003). Also note that even when the MRI operates, MHD turbulence may not act like a Navier–Stokes viscosity.

For the standard simulations in this work to fit the parabolic wake (see Equation (3)) in the box, we use box size 12h × 64h, thus the overall grid resolution in our runs is about 3072 × 16384 (for the standard resolution of 256/h). In a few cases with very small M_p, we extend the simulation box size to 20h × 156h to trace linear wave evolution farther out in x, since smaller M_p delays shock formation (see Equation (9)). Our simulations are run for at least 10, and in some cases up to 50, orbital periods. Figure 1 presents a typical snapshot of one of our simulations showing the density structure and the spiral waves.

Zoom In Zoom Out Reset image size

We use the following boundary conditions (BCs). On x (radial) boundaries, we keep values of all physical variables in ghost zones fixed at their respective unperturbed Keplerian values (i.e., keep the ghost zones as their initial states), and the waves leave through the x boundaries when they reach the edge. In our shearing sheet simulations, experiments show that this x BC has less wave reflection and outflowing of the fluid than the conventional outflow x BC, in which case the ghost zones are copied from the last actively updated column of cells of the simulations at every time step (see also Muto et al. 2010). We do not expect our adopted radial BC to affect wave evolution since significant radial fluid motions are not expected to arise in our simulations anyway.

On the y (azimuthal) boundary, we experimented with two BCs: the conventional outflow BC, as described above, and an inflow/outflow BC. In the latter case, the variables in the ghost zones are fixed at their initial values if they are the physical "inflow" boundaries (the regimes y < 0, x < 0 and y > 0, x > 0), or copied from the last actively updated row of cells if they are the physical "outflow" boundaries (the regimes y > 0, x < 0 and y < 0, x > 0). For the purposes of the current paper, the difference between these two BCs is not significant and resultant density profiles and torque calculation are almost identical (of the level of 10⁻³). We found that with the pure outflow BC, fluid entering the simulation box accumulates some non-zero velocity perturbation on top of the pure linear shear velocity profile. This affects calculation of variables derived from the simulated velocity field, such as potential vorticity, which are analyzed in Paper II. We use the inflow/outflow y BC for our simulations.

At the start of a simulation run we instantaneously turn on the potential of the planet in the center of the linearly sheared fluid flow with uniform surface density. This gravitational perturbation immediately excites an inner and an outer density wave propagating away from the planet, as well as strong transients: vortices that appear near the planet but travel on horseshoe orbits away from it. Before measuring the wave properties, we run simulations for ∼15 orbital periods to let these time-dependent structures move away from the planet. We also tried gradually turning on planetary gravity by linearly increasing strength of the potential from zero to its full value within several orbital periods (Muto et al. 2010; Li et al. 2009), but did not find this trick as effective at removing transient structures as simply waiting for them to move away from the planet.

We now present our main results and compare them with the predictions of linear theory. To mitigate possible nonlinear effects we consider only very low planetary masses in our calculations, ranging between 3.2 × 10⁻²M_th and 2.8 × 10⁻³M_th. According to Equation (7), at 1 AU this mass range corresponds to M_p between 0.4 M_⊕ and 3 lunar masses (note that they are sound speed dependent). As the basis for comparison we primarily use the distribution of the perturbed density (Section 4.1) and the evolution of the AMF carried by the density wave (Section 4.2). We carry out linear theory verification both in physical (coordinate) space and in Fourier space.

4.1. The Density Wave Profile

The most basic and direct comparison between numerical simulations and theory can be performed using the spatial distribution of the perturbed fluid variables. In this work we choose the perturbed surface density δΣ ≡ Σ − Σ₀ as the primary variable for comparison. Since δΣ is spatially concentrated along a narrow wake it makes sense to compare with theory both the overall shape of the wake in x–y coordinates and the density distribution across the wake.

We determine the wake shape in the following way. At each value of x we find the value of y, called y_w(x), at which δΣ reaches its maximum value δΣ_max(x). We then compare the run of y_w(x) with the theoretical prediction (3), which applies in the shearing sheet geometry (GR01). The results are shown in Figure 2. The two curves agree well far from the planet, at |x| ≳ h. The discrepancy at |x| ≲ h is not surprising because in this region (1) the density wave is not yet fully formed and (2) the dispersion relation (1) is significantly affected by the κ² term, which causes the wake profile to deviate from the theoretical prediction (3) valid far from the planet, see Ogilvie & Lubow (2002) for details.

Zoom In Zoom Out Reset image size

Next we investigate the evolution of the wave amplitude as it travels away from the planet by looking at the behavior of the maximum amplitude δΣ_max(x) ≡ δΣ(x, y_w(x)) as a function of x. We plot this quantity scaled by Σ₀(M_p/M_th) (this normalization removes the dependence of δΣ_max on Σ₀ and M_p) in Figure 3(a).

Zoom In Zoom Out Reset image size

At small x ≲ h the value of δΣ_max is large and decreases with x. This behavior has nothing to do with the density wave since the region in the immediate vicinity of the planet represents a quasi-static atmosphere that forms inside the planetary potential well right after the planetary gravity is switched on. This structure is rarely mentioned in the analytical calculations of planet–disk interaction but it shows up prominently in realistic simulations. In hydrostatic balance (in the absence of background velocity caused by differential rotation), the density profile in the vicinity of the planet should be given by (assuming isothermal EOS with constant sound speed c_s)

$$\begin{equation} \Sigma _{\rm atm}(\rho)=\Sigma _0 e^{-\Phi _{p}(\rho)/c_s^2}. \end{equation} \tag{ 18 }$$

In real disks, the background velocity field distorts the atmospheric profile from the circularly symmetric (with respect to planet) shape predicted by this equation. This explains why Σ_atm(x) shown in Figure 3(b) slightly overestimates δΣ_max at small x.

The effect of the atmosphere on δΣ_max rapidly decreases beyond several h from the planet, and δΣ_max starts to rise with increasing x. This is because density wave excitation occurs at |x| ∼ h and from now on the behavior of δΣ_max is determined by the wave-like density perturbation. Conservation of the AMF F_H carried by the wave forces the wave amplitude to increase with distance ∝x^1/2 (GR01). As Figure 3(a) shows this scaling agrees reasonably well with the numerically computed δΣ_max(x) far from the planet.

Note that initially, at separations of several h from the planet, δΣ_max increases with x faster than the theoretical x^1/2 scaling, which is caused by the still ongoing accumulation of the AMF F_H at these separations. In other words, at this location F_H has not yet reached its asymptotic value given by Equation (6); see Section 4.2 for more details.

Far from the planet, beyond x = l_sh (l_sh = 7.5h in the case shown in Figure 3(a)), the peak δΣ rapidly goes down relative to the analytical x^1/2 scaling. This behavior is caused by the appearance of the shock at l_sh and subsequent dissipation of the AMF carried by the wave. These nonlinear processes are explored in detail in Paper II.

Finally, we examine the evolution of the density profile across the wake as a function of radial separation x. We do this by making an azimuthal cut through the density field at fixed x and shifting the resultant one-dimensional density profile in y by y_w given in Equation (3). To eliminate the increase of δΣ resulting simply from the AMF conservation we normalize the density perturbation by x^1/2. We additionally scale δΣ by Σ₀(M_p/M_th), as in Figure 3(a). As mentioned before, in linear theory the normalized density profile should be independent of M_p at a fixed separation x. This is indeed seen in Figure 4(a), where we plot scaled profiles of δΣ computed for several values of M_p at x = 1.33h and find no significant difference between them. The adopted value of x is small enough for the linear theory to apply even for the highest M_p (≈0.03M_th) used in making this figure, i.e., x is always considerably smaller than l_sh (cf. Figure 3(b) and the discussion in Section 4.4).

Zoom In Zoom Out Reset image size

Previously, GR01 computed the evolution of the density profile as a function of x in the linear regime. In Figure 3(b), we provide an analog of Figure 1 of GR01 by showing the scaled density cuts at the same values³ of x—1.33h, 2.67h, 4h, 5.33h—as in GR01. This particular calculation uses a very small planet mass M_p = 3.7 × 10⁻³M_th corresponding to l_s ≈ 8h (see Equation (9)), which puts the values of x used in making this figure well inside the shock position and reduces the impact of nonlinearity on the density profile; see Section 4.4.

The effective width of the density profile in this figure does not significantly vary with x and stays at the level of several h at all times. We note, however, that this property holds only for the density cuts passing through the wake at fixed x, as shown in Figure 3(b). On the contrary, the width of the density profile cutting through the wake at fixed y (not shown here) is shrinking as |x|⁻¹, which follows directly from the dispersion relation (1) since k_r∝|x| for x ≳ h.

Our numerical results agree with linear calculations of GR01 quite well, with quantitative differences at the level of 10% or less. In particular, we are able to reproduce several subtle features of the wake profile evolution found in GR01, which are highlighted by arrows in Figure 3(b). These are (1) the increase of δΣ with x in the density trough in front of the wake, (2) the decrease of δΣ with x in the density trough behind the wake, and (3) the slight drift of density profile with respect to theoretical wake position y = y_w(x) given by Equation (3) as x increases. There are also some minor differences between Figure 3(b) and the results presented in GR01, e.g., the lower (by about 10%) amplitude of δΣ far from the planet, at x ≳ 3h, in our case (note that the vertical axis of Figure 3(b) is different from the vertical axis in Figure 2 of GR01, since we use M_th and h instead of M₁ and l = 2h/3 as mass and length units). We speculate that this may be explained by the gradually accumulating nonlinear effects in our simulations, which are absent in linear calculations of GR01.

Previously, Muto et al. (2010, Figure 6) have shown density profiles resulting from their 2D hydrodynamical simulations in cylindrical geometry, which are similar to those presented in Figure 3. However, their calculations were carried out for planetary masses (minimum M_p = 5 × 10⁻²M_th) larger than those used in our work (maximum M_p = 3.2 × 10⁻²M_th), which does not allow them to isolate linear effects clearly. Also, Muto et al. (2010) did not attempt quantitative comparison of their numerical results with the linear theory predictions.

4.2. The Angular Momentum Flux and Torque Density

We now investigate the behavior of the AMF carried by the wave F_H(x) as a function of x. Linear calculations of this quantity in Fourier space have been first performed in GT80, and later refined in Artymowicz (1993a, 1993b) and Korycansky & Pollack (1993). The torque density dT_H/dx which is a spatial derivative of F_H(x) has been inferred from simulations in Bate et al. (2003) and D'Angelo & Lubow (2008, 2010), although no direct comparison with the linear theory was provided in these studies.

The analytical expression (6) for the integrated AMF scales linearly with the planetary semimajor axis r_p. In the shearing sheet geometry employed in our simulations r_p is ill-defined, and it makes sense to redefine AMF as $\tilde{F}_H\equiv F_H/r_p$. Subsequently, we will drop the tilde for brevity and denote such normalized momentum flux as F_H. We compute F_H(x) from our data according to the following definition:

$$\begin{equation} F_H(x) = \Sigma _0\int _{-\infty }^\infty dy\delta v_y v_x, \end{equation} \tag{ 19 }$$

where v_x is the radial velocity of the fluid and δv_y is the velocity perturbation with respect to the background shear profile in the azimuthal direction.

We compute the torque density (normalized by r_p) exerted by the planet on the fluid at separation x as

$$\begin{equation} \frac{dT_H}{dx} =-\int _{-\infty }^\infty dy\delta \Sigma \frac{\partial \Phi _{p}}{\partial y}, \end{equation} \tag{ 20 }$$

and the integrated torque accumulated by the wave at x is just

$$\begin{equation} T_H(x) = \int _0^x\frac{dT_H}{dx}dx. \end{equation} \tag{ 21 }$$

Conservation of angular momentum ensures that the derivative of Equation (19) (the slope of the AMF) coincides with the torque density (21) in the absence of dissipation.

4.2.1. Comparison with Theory in Fourier Space

We first compare our simulation results with linear theory in Fourier space by studying the behavior of the torque cutoff function F_H(m)/F^WKB_H(m); see Section 2. Previously, GT80 computed a variation of this quantity as a function of m, and the result is shown in Figures 2 and 3 of their paper, providing a basis for comparison. Artymowicz (1993a) subsequently refined this calculation by arguing that evaluating the planetary potential at separations slightly displaced from the classical Lindblad resonance position (2) should lead to more accurate results. Adopting this approach Artymowicz (1993b) provided a simple analytical prescription (Equation (25) of Artymowicz 1993b) to reproduce the behavior of F_H(m)/F^WKB_H(m), which can also be compared with our calculations.

Substituting v_x and δv_y in the form of their Fourier integrals in Equation (19) and manipulating the resulting expression one can write F_H(x) as the integral over the azimuthal wavenumber k_y:

$$\begin{eqnarray} && F_H(x)=\int _0^{\infty }dk_y F_{H,k}(x), \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} && F_{H,k}(x)=4\pi \Sigma _0\,[{\rm Re}\,(v_{x,k})\,{\rm Re}\,(\delta v_{y,k})+{\rm Im}(v_{x,k}){\rm Im}(\delta v_{y,k})],\nonumber\\ \end{eqnarray} \tag{ 23 }$$

where v_{x, k} and δv_{y, k} are the Fourier transforms of v_x and δv_y, which themselves are functions of x. Even though the available analytical calculations of the torque cutoff function were performed only in the limit x ≫ 1, we chose to retain the dependence of F_H(x) on x to explore the evolution of the harmonic content of the AMF with distance from the planet.

Instead of the discrete theoretical WKB Fourier harmonics of the flux F^WKB_H(m) given by Equation (5) we use the following continuous version:⁴

$$\begin{eqnarray} &&F_{H,k}^{{\rm WKB}}=\frac{4}{3}\Sigma _0(k_yh)^2 \left(\frac{{\it G M}_{p}}{c_s}\right)^2\left[2K_0(2/3)+K_1(2/3)\right]^2, \nonumber\\ \end{eqnarray} \tag{ 24 }$$

where we also normalized the final expression by r_p.

We can now compute the real and imaginary components of v_{x, k} and δv_{y, k} using our numerical data and then obtain F_{H, k}(x) from Equation (23). Dividing the result by F^WKB_{H, k} provides us with the ratio F_{H, k}(x)/F^WKB_{H, k}, which can then be compared with the existing theoretical calculations of the torque cutoff. This is done in Figure 5 where we plot the numerical F_{H, k}(x)/F^WKB_{H, k} computed at different separations from the planet against the semi-analytical calculations of the same quantity (in the limit |x| ≫ h) performed by GT80 and Artymowicz (1993b). These calculations were done for rather small M_p = 1.2 × 10⁻²M_th (corresponding to l_sh ≈ 5h) allowing us to see the linear wave evolution out to large separations.

Zoom In Zoom Out Reset image size

One can see in Figure 5(a) that the agreement between theory and simulations is generally quite good. Our results seem to agree better with GT80 torque cutoff prescription than with Artymowicz (1993b). However, at the level of accuracy available to us we are not able to firmly discriminate between these two torque cutoff prescriptions. At intermediate values of μ = k_yh ≈ 0.5–5 our numerical results for F_{H, k}(x)/F^WKB_{H, k} are essentially independent of x and pass between the two aforementioned analytical torque cutoff prescriptions. Only weak sensitivity of F_{H, k}(x)/F^WKB_{H, k} on x is expected since beyond |x| = 2h the power in azimuthal harmonics corresponding to this range of k_yh should have already been fully accumulated by the density wave, and no evolution should arise. Good agreement with the GT80 results is also illustrated in Figure 5(b), which shows the torque cutoff function (multiplied by μ²) in linear space, as in Figure 3 of GT80.

The situation is different for smaller values of k_yh ≲ 0.5, where the numerical curves fall below the analytical asymptotic behavior F_{H, k}/F^WKB_{H, k} → 1 as k_yh → 0. Also, there is noticeable evolution of F_{H, k}(x)/F^WKB_{H, k} with x, with non-zero power extending down to smaller and smaller values of k_y as the wave travels further from the planet. This is because the low-order harmonics corresponding to small k_yh contribute to the wave flux predominantly at separations larger than ∼h (Korycansky & Pollack 1993); see Equation (2). As a result, the farther the wave travels, the more power gets collected by the wave from these lower-order azimuthal harmonics of the planetary potential, but the very small values of k_y ≲ 0.05 never contribute to the wave flux even at |x| = 5h.

Our numerical results also exhibit noticeable disagreement with theory and clear evolution with x at large values of k_yh ≳ 5, and we comment on their origin in Section 4.4.

4.2.2. Comparison with Theory in Physical Space

We now study the behavior of the AMF and torque density in physical space. Understanding the spatial distribution of the latter quantity is important for properly computing the AMF in disks with non-uniform distribution of Σ in the vicinity of the planet, and is thus important for understanding the early stages of gap formation by massive planets.

Here we look at disks with uniform distribution of Σ. In Figure 6, we show the behavior of the AMF F_H(x) computed according to the definition (19) and also of the accumulated torque T_H(x) given by Equations (21). We display our results for four different values of M_p, normalizing F_H(x) and T_H(x) by Σ₀(GM_p)²/(hc²_s)—according to linear theory the shape of the resultant curves should then be independent of M_p. Theoretical scaling of the torque ∝M²_p is indeed largely confirmed by this figure, but see Section 4.4 for more details.

Zoom In Zoom Out Reset image size

To facilitate the comparison with the torque calculation, AMF curves have all been vertically shifted by a small offset to cancel the non-zero value of F_H(0) at x = 0 resulting from the residual gas motion in the horseshoe region and a small intrinsic non-zero starting value of the AMF (R. R. Rafikov & C. Petrovich 2012, in preparation). This brings them in perfect accord with T_H(x) curves close to the planet, which is expected to be the case in the absence of dissipation.⁵ At large separations, beyond l_sh, the F_H(x) curve starts falling below the T_H(x) curve, which is caused by shock formation and dissipation of the wave AMF past this point. Quite naturally, this effect is more pronounced for higher M_p, e.g., M_p = 3.2 × 10⁻³M_th and M_p = 1.2 × 10⁻³M_th, corresponding to smaller l_sh. For smaller masses shocks occur outside the simulation box and essentially no difference between F_H(x) and T_H(x) can be seen. The post-shock behavior of the AMF will be studied in Paper II.

In Figure 7, we plot the torque density dT_H(x)/dx (20) as a function of x. The overall shape of the curve is consistent with previous numerical calculations of the same quantity in cylindrical coordinates performed by Bate et al. (2003) and D'Angelo & Lubow (2008, 2010). As expected from linear theory, beyond ≈2h torque density drops dramatically, and the AMF rises only weakly, with |x|, which is evident from the flattening of AMF curves in Figure 6. This justifies the localization of the wave excitation to the immediate vicinity (within ∼2h) of the planet, as described in Section 2.

To the best of our knowledge, no quantitative analytical description of how F_H(x) should vary with distance from the planet exists in the literature, despite its potential importance for the gap opening problem. To fill this gap we first tried the following theoretical prescription F^LR_H(x) motivated by the calculations of GT80: at each x we compute the lowest azimuthal wavenumber μ_min(x) for which the location of the Lindblad resonance x_{L, m} < x using Equation (2) and then calculate the AMF as the sum of all Fourier contributions F_H(m) given by Equation (4) that correspond to m > (r_p/h)μ_min(x). The result, which we denote by F^LR_H(x), is given by

$$\begin{eqnarray} F_H^{\rm LR}(x)&=&\frac{4}{3}\frac{\Sigma _0({\it G M}_{p})^2}{c_s^2h}\left[2K_0(2/3)+K_1(2/3)\right]^2\nonumber\\ &&\times \int _{\mu _{{\rm min}}(x)}^\infty d\mu \mu ^2 \frac{F_H(\mu)}{F_H^{\rm WKB}(\mu)}, \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} && \mu _{{\rm min}}(x)=\frac{2}{3}\frac{h}{x}, \end{eqnarray} \tag{ 26 }$$

where we adopt the torque cutoff function F_H(μ)/F^WKB_H(μ) calculated in GT80. Note that F^LR_H(x) is in fact independent of the adopted value of r_p since its calculation can be rephrased fully in terms of μ rather than m.

This prescription essentially assumes that the AMF contribution F_H(m) corresponding to a particular azimuthal wavenumber m gets picked up by the density wave in a step-like fashion, solely at a single discrete location corresponding to the position of the mth Lindblad resonance. Previously, this recipe was used in GT80 to compute the asymptotic behavior of the torque density in the limit |x| ≳ h, considering azumuthal harmonics with 1 ≲ m ≲ r_p/h. Here, we essentially extend this prescription to the case of arbitrary m and |x|. Linear calculations by Korycansky & Pollack (1993) show that this approximation is not very accurate, and AMF contribution due to each potential harmonic is in fact accumulated by the wave over an extended range of x. Thus we should not expect F^LR_H(x) given by Equation (25) to accurately represent the real coordinate dependence of the AMF in the linear regime. Nevertheless, this prescription still provides a useful reference point and we plot it in Figure 6(a) (the dotted line, denoted by LR theory). We also plot in Figure 7 the theoretical prescription for the torque density dT^LR_H(x)/dx obtained by differentiating Equation (25) with respect to x.

Zoom In Zoom Out Reset image size

One can see from these figures that the simple-minded theoretical prescription (25) overestimates the torque density at small x and underestimates it at larger x. While the areas under the numerical and theoretical curves in Figure 7, which represent the full accumulated AMF, are close to each other, the profiles of the two curves are quite different. As a result, F^LR_H(x) initially rises faster than the numerical F_H(x) in Figure 6, but eventually the two asymptote to a similar final value. The small remaining difference is caused by both the intrinsic numerical inaccuracy and the fact that unlike the AMF theoretical curve, the numerical AMF curve has been shifted to have a zero starting point.

As the next level of approximation to the behavior of F_H(x) and T_H(x) we used the semi-analytical calculations of these quantities in the linear regime by R. R. Rafikov & C. Petrovich (2012, in preparation). The run of corresponding T_H(x) is displayed in Figure 6 by the dash-dotted curve (the curve has also been shifted downward to cancel the non-zero starting value, and we use the label "Linear theory" to indicate R. R. Rafikov & C. Petrovich's work in all the figures) and clearly demonstrates good agreement between the simulation results and the semi-analytical linear theory. Theoretical torque density based on this T_H(x) is shown in Figure 7 and also agrees well with the numerical results, much better than the derivative of F^LR_H(x). This additionally emphasizes the point that assigning the planetary torque produced by a particular potential harmonic to a single location corresponding to the respective Lindblad resonance (as done in Equation (25)) is not a very accurate procedure.

Finally, we also compare our numerical results with another quantity that can be computed in the linear approximation, namely the "pseudo-AMF" f_J introduced in GR01 as

$$\begin{eqnarray} &&f_J(x)=\frac{2}{3}\frac{c_s^3}{\Omega \Sigma _0|x|}\int _{-\infty }^\infty \left(\delta \Sigma \right)^2 dy. \end{eqnarray} \tag{ 27 }$$

This quantity reduces to F_H(x) at large separations |x| ≳ h from the planet, where δΣ is determined solely by the density wave propagating away from the planet. However, for |x| ≲ h the pseudo-AMF is strongly affected by the presence of the aforementioned atmosphere around the planet (see Section 4.1), which is not a propagating density perturbation and cannot be associated with the density wave (this explains the prefix "pseudo-").

We compute the theoretical value f_J(x) in the linear approximation using unpublished data from GR01 and also derive f_J(x) from our numerical data using the definition (27). The comparison between the two is shown in Figure 8 for very small M_p = 3.7 × 10⁻³M_th. One can indeed see that at small |x| the pseudo-AMF increases with decreasing |x|, analogous to δΣ in Figure 3(a), which is explained by the presence of the atmosphere accumulated in the planetary potential well. Comparison with Figure 6 also shows that f_J(x) indeed reduces to F_H(x) at large |x|. Most importantly, the agreement between the theoretical f_J(x) derived from linear theory and the numerical f_J(x) is very good in the whole range of x explored, which provides additional verification of linear theory of the planetary density wave evolution.

Zoom In Zoom Out Reset image size

4.3. Sensitivity to EOS

4.4. Emergence of the Nonlinear Effects

4.5. "Negative Torque" Phenomenon

We now explore how the results presented in the previous section and their agreement with analytical theory are affected by the variation of purely numerical parameters in our simulations. Table 1 lists the numerical parameters that we varied and the values we explored. Values corresponding to our standard case are indicated in boldface.

Table 1. Parameter Space of the Simulations

Parametersa	Range
Riemann solver (flux function) used	Roe, HLLC
Order of accuracy	2, 3c, 3p
Boundary conditions in y	Outflow, Inflow/Outflow
Resolution of the simulation (cells per h)	64, 128, 256
Planetary potentialb	Φ(2)p, $\bf \Phi _{p}^{\bf (4)}$, Φ(6)p
Softening length	1/8, 1/16, 1/32
Equation of states of the fluid	$\bf \gamma =1$, γ = 5/3
Mass of the planet Mp/Mth	3.2 × 10−2, 1.2 × 10−2, 5.7 × 10−3, 3.7 × 10−3, 2.8 × 10−3

Notes. ^aSee Section 5 for details. ^bSee Section 3.

Download table as:  ASCIITypeset image

5.1. Numerical Solver and Order of Accuracy

In our simulations we compare two different Riemann solvers—Roe's linearized solver (Roe 1981) and Harten–Lax– van Leer–Contact (HLLC; Toro 1999)—with three different algorithms for the spatial reconstruction step (Stone et al. 2008): second order with limiting in the characteristic variables (denoted 2 in this work), which is the predominant choice in literature in this field, and third order with limiting in either the characteristic variables (3c), or in the primitive variables (3p). We find that Roe and HLLC solvers yield nearly identical results both in terms of the density profile and in terms of F_H(x) and T_H(x) (differences are at the 0.1% level).

The effect of using a different order of accuracy on the profile of the density wave is shown in Figure 9(a), where we plot azimuthal density cuts obtained with different numerical settings at x = 1.33h (together with semi-analytical profile computed by GR01 in the framework of linear theory) and x = 4h. One can see that density perturbation is rather insensitive to variation of the order of accuracy.

Zoom In Zoom Out Reset image size

The sensitivity of the AMF calculation to varying order of accuracy is illustrated in Figure 10(a), where the linear theory curve based on the semi-analytical calculation by R. R. Rafikov and C. Petrovich (2012, in preparation) is overplotted for reference. One can see that the higher order of accuracy (3c and 3p) results in a slightly lower asymptotic value of torque and AMF, and the difference in the asymptotic value for the two between the 2 case and the 3p case is about 5%. Furthermore, the 3c and 3p cases agree with each other on the position where the AMF curve starts to deviate from the torque, ∼4h, which is still somewhat smaller than the theoretical prediction l_sh ≈ 5h in this case due to the relatively low resolution. However, the second order of accuracy moves this point inward, which means that the angular momentum carried by the wave starts to dissipate earlier in this case. In particular, we find that using a second order accuracy solver compared with the third order accuracy solver has the similar effect of decreasing the resolution by a factor of two (in terms of advancing the displacement of the AMF-torque separation point; see Section 5.2 and Figure 10(b)).

Zoom In Zoom Out Reset image size

Finally, we also note that calculations with 3p order of accuracy produce considerably noisier velocity field (this explains the noisy corresponding AMF curve in Figure 10(a)) than the other two cases explored. This difference is important for calculation of the velocity-based variables such as potential vorticity, which is used in Paper II as the means of shock detection. For these reasons in our current simulations we use the Roe solver with 3c reconstruction.

5.2. Resolution

5.3. Planetary Potential and Softening Length r_s

5.4. Summary of Convergence Study

We conclude that the overall shape of the density wave profile is generally less sensitive to the variations of the numerical parameters (unless these values are very different from the ranges explored in this work; see below) than the spatial dependence of the AMF or integrated torque. This emphasizes the importance of these integrated quantities as diagnostics of various subtle effects influencing the wave properties.

In the course of this numerical exploration we have also discovered that the "optimal" value of a particular numerical parameter (e.g., softening length r_s) showing closest agreement with the linear theory depends on the values of other numerical parameters (e.g., resolution, type of the potential, etc.) and on M_p. In other words, there is no universal best choice for each numerical parameter of the simulations. In both Figures 9 and 10, only the relative position of the curves representing different parameters matters, not their exact locations and the discrepancies between them and the linear theory. On the other hand, as we demonstrate in Paper II, high order of accuracy, high resolution, and highly accurate form of the potential with small softening length are critical to properly resolve the nonlinear stage of the wave evolution.

Hydrodynamic simulations of the disk–planet interaction must often be run for many orbital periods. This is the case, e.g., in studies of planetary migration or gap opening, in which substantial effects appear on timescales of hundreds to thousands of orbital periods (depending on M_p). This severely restricts the choice of resolution and softening length at which such simulations can run. Typical resolutions in global disk simulations found in the literature (Kley 1999; Bryden et al. 1999; Nelson et al. 2000; D'Angelo et al. 2002, 2003; Li et al. 2009; Yu et al. 2010; Muto et al. 2010) range from 32/h to 1/h, and r_s is usually taken to be between 0.1h and h, in combination with a second order of accuracy solver and the second (see Equation (15)) or even lower order representation of the planetary potential. On the other hand, accurate description of the migration rate and gap opening may only be possible if properties of the density wave excited by the planet are properly captured by the simulation.

To test how reliably the wave structure can be reproduced in global simulations we compare two sets of simulations in our local setting with M_p = 1.2 × 10⁻²M_th. We run one simulation with the typical numerical parameters for local shearing box simulations found in recent literature (typical global simulations in literature use even smaller resolution and r_s): second order of accuracy, resolution 32/h, and Φ⁽²⁾_p potential with r_s = h/4. Another simulation is run with our fiducial parameters—third order of accuracy, resolution 256/h, and Φ⁽⁴⁾_p potential with r_s = h/16—and their results are compared in Figure 11.

Zoom In Zoom Out Reset image size

There are noticeable differences between the two runs. The simulation with typical literature parameters produces a smoother density profile, which deviates from analytical solution by about 20%; see Figure 11(a). At the same time, our fiducial run shows deviations from the analytical profile only at the level of 1%.

The difference is even more pronounced when we compare the AMF and torque behavior in Figure 11(b). Our fiducial simulation yields asymptotic values of F_H(x) and T_H(x) which are within several percent of the expected theoretical value (6). As expected, F_H(x) agrees with T_H(x) all the way until x ≈ 5h, which is precisely the shock location for the M_p used, according to Equation (9). On the contrary, in the typical literature simulation AMF F_H(x) starts deviating from the integrated torque curve T_H(x) already at x ≈ 2.5h, significantly in advance of the expected shock position (a factor of two). This indicates that dissipation and transfer of the angular momentum carried by the wave to the disk fluid start earlier in the typical literature case than they should in reality. The most likely explanation for this behavior is the high level of numerical viscosity in the typical literature run. Moreover, the absolute asymptotic values of F_H(x) and T_H(x) disagree with the theoretical prediction for asymptotic F_H by ∼23% (comparing with ∼2% in our fiducial simulation), which is quite significant. We note that the differences between our fiducial and the typical literature simulations persist also in experiments with larger planet masses, M_p ≲ M_th, while the linear regime of wave excitation still holds.

An underestimate of the wave AMF in the low-resolution simulations may result in an underestimate of the planetary migration rate in global simulations. In our shearing sheet setup we cannot investigate the effect of resolution on the planetary migration rate: by design one-sided torques exerted by the inner and outer parts of the disk on the planet exactly cancel each other, while the migration speed relies on the small imbalance between them. But the very fact that the one-sided torques in the low-resolution case deviate by tens of percent from the high-resolution case suggests that the migration speed should be adversely affected by poor resolution at the same level. Furthermore, the relative imbalance between the one-sided torques can also be a function of resolution, potentially exacerbating the problem.

The discrepancy in the AMF should also affect the fidelity of gap opening by the planet in typical literature global simulations. Lower AMF carried by waves means that the planet is less effective at repelling gas away from its semimajor axis, which would require higher M_p to open a gap. In addition, the gap opening process is sensitive to the spatial distribution of the AMF dissipation (Rafikov 2002b). As a result, spurious dissipation of wave AMF prior to shock formation in low-resolution simulations and the accelerated AMF decay after the shock may introduce artificial effects in the gap opening picture. To summarize, any numerical studies of processes in which density wave production and dissipation plays important roles need to use a high order of accuracy, high resolution, and accurate representation of the planetary potential with small softening length.

In this section, we describe a commonly ignored numerical problem that we discover in disk–planet interaction simulations. Namely, we find that if the time step dt of the simulation is too large, the code cannot resolve the motion of the fluid in regions where the fluid experiences large gravitational acceleration. This issue applies to simulations in general, but it is especially problematic when the orbital advection algorithms are implemented (the Fast Advection in Rotating Gaseous Objects, FARGO, algorithm; Masset 2000).

Masset (2000) introduced a modification of the standard transport algorithm (FARGO) for simulations of differentially rotating systems, which significantly speeds up the calculations. By getting rid of the average azimuthal velocity when applying the Courant condition, this technique results in a much larger time step dt, which is limited by the small perturbed velocity, than the usual procedure, where dt is limited by the full fluid velocity dominated by the differential rotation. FARGO has been implemented in Athena by Stone & Gardiner (2010).

In a shearing box without planets, the dynamical timescale is Ω⁻¹ and is uniform throughout the box. However, when the planet is introduced, the dynamical timescale is spatially varying, and could be characterized by the free-fall time, which is the time it takes a fluid element to fall on the planet assuming a constant acceleration at the grid point. For gravitational potential in the form of Equation (16), the free-fall time is

$$\begin{equation} t_{{\rm ff}}=\sqrt{\frac{2\big(\rho ^2+r_{\rm s}^2\big)^{5/2}}{{\rm G}M_{p}\big(\rho ^2+2.5r_{\rm s}^2\big)}}. \end{equation} \tag{ 28 }$$

Note that t_ff depends on a specific form of the gravitational potential, and a more smoothed potential results in a larger t_ff. The smallest free-fall time in the entire domain, which is the limiting timescale for the simulation, occurs at the grid points that are adjacent to the planet, which have the smallest ρ ($1/\sqrt{2}$ of the cell size r_c in Athena; the separation where the smallest t_ff occurs also depends on the form of potential). In our simulations, we always keep r²_s ≫ r_c² (usually r_s = 8r_c), so the smallest free-fall time is

$$\begin{equation} t_{{\rm ff},s}\approx \sqrt{{\frac{0.8r_{\rm s}^3}{{\rm G}M_{p}}}}. \end{equation} \tag{ 29 }$$

On the other hand, according to the Courant condition when varying the resolution the time step dt∝r_c. So, the ratio t_{ff, s}/dt scales with the numerical parameters and M_p as

$$\begin{equation} \frac{t_{{\rm ff},s}}{dt}\propto \frac{r_{\rm s}^{1.5}}{\sqrt{M_{p}} r_c}. \end{equation} \tag{ 30 }$$

To properly resolve fluid dynamics in the vicinity of the planet, the ratio t_{ff, s}/dt should be kept above a certain threshold.

Our calculations without orbital advection and with dt determined by the Courant condition usually have t_{ff, s}/dt ∼ 150. However, when we turn on the orbital advection algorithm, dt typically increases by a factor of ∼10, and the ratio t_{ff, s}/dt decreases by the same factor. We find that when the time step is too large to properly resolve the fluid motion around the planet, the numerical results will be incorrect, as described below.

Figure 12 shows the azimuthal density profiles δΣ at x = 1.33h for a set of simulations using identical numerical parameters and orbital advection algorithm but with different t_{ff, s}/dt ratios, which we achieve by manually setting dt to be a fraction of the dt set by the Courant condition in FARGO. The two cases with the highest t_{ff, s}/dt yield the density profile in agreement with the theoretical prediction, also demonstrating the convergence of the result at high t_{ff, s}/dt. However, the density profiles in simulations with lower t_{ff, s}/dt clearly deviate from theory, with smaller t_{ff, s}/dt leading to larger discrepancy. In the run with the smallest used t_{ff, s}/dt = 18 which corresponds to the dt set by the Courant condition in FARGO (representing maximum FARGO speedup), the density peak even becomes a density trough. By experimenting with Athena, we generally find that our simulations of the disk–planet interaction produce converged results agreeing with theory only when

$$\begin{equation} t_{{\rm ff},s}/dt\gtrsim 100, \end{equation} \tag{ 31 }$$

while this critical number may depend on the numerical method used in a particular code or some other numerical conditions. We note that in principle all disk–planet simulations have to satisfy condition (31) to ensure correct results, but in experiments we find that simulations with the orbital advection algorithm tend to violate this condition much more easily than simulation without the orbital advection algorithm.

Zoom In Zoom Out Reset image size

Apart from incorrectly reproducing density perturbation, disk–planet simulations violating the t_{ff, s}/dt constraint exhibit other serious problems. In particular, simulations with small t_{ff, s}/dt do not achieve a steady state. We illustrate this in Figure 12 by showing the density profile at half simulation time for two cases. While the t_{ff, s}/dt = 144 case reaches a steady state and develops a time-independent density profile at half-time, the t_{ff, s}/dt = 72 case does not, and its density perturbation decreases with time.

In simulations with t_{ff, s}/dt lower than the critical value, the fluid element approaching the planet at small |x|, instead of moving on a horseshoe orbit, gets trapped by the planetary potential, leading to the formation of a fluid loop around the planet. As simulation progresses the rotational velocity of the fluid loop becomes higher and higher, and eventually the centrifugal force evacuates the region around the planet. This effect, which is a purely numerical artifact of the too small t_{ff, s}/dt, is visually very similar to (and can easily be confused with) the gap opening phenomenon, which is a real physical effect. Simulations exploring the gap opening process with large planet mass and low resolution are likely to have small t_{ff, s}/dt, violating condition (31) (see Equation (30)), which may have detrimental effects on the results. We note that using Athena with orbital advection algorithm, numerical simulations with M_p ∼ M_th and a large r_s = h/4 must have an effective resolution at least 64/h in the vicinity of the planet to avoid this problem (a more smoothed potential may weaken this condition though). However, again we emphasize that the critical value of t_{ff, s}/dt in simulations with the same numerical parameters but using a different code and solver may be different from ours.

The t_{ff, s}/dt threshold severely limits the ability of FARGO to speed up calculations of disk–satellite interaction. For most of our simulations, t_{ff, s}/dt is already rather close to the threshold value without orbital advection algorithm, so there is not much room left to increase dt, which is what the orbital advection algorithm aims to achieve (all our simulations presented outside this section are done without orbital advection technique). However, in many other studies of differentially rotating systems, such as the investigation of MRI in accretion disks, the point-mass objects are absent and the characteristic dynamical timescale is always long enough for the orbital advection algorithm to be an extremely useful tool for speeding up the simulations. At last, we note that both dt and t_{ff, s} may vary during the simulation, and the constraint (31) on their ratio should be satisfied in the entire domain throughout the simulation time.

We have conducted a series of hydrodynamical simulations to study the details of the gravitational interaction between low-mass planets and a protoplanetary disk, and to test predictions of the linear theory of density wave evolution. Our simulations assume local shearing sheet geometry and are carried out in 2D at very high resolution to reduce the effect of numerical viscosity. We focus on both the excitation of the density waves and their propagation away from the planet in the linear regime. To mitigate the effects of nonlinearity we consider very small planetary masses, starting at 0.4 M_⊕ and going down to 3 lunar masses at 1 AU.

We extract the spatial distribution of the density perturbation induced by the planet from our simulations and compare it with the predictions of linear theory. We find good agreement between the two, at the level of several percent when high resolution (typically 256/h) is employed. We also investigate the spatial dependence of the AMF carried by the wave and the distribution of torque induced by the planet on the disk, again finding good agreement with theoretical predictions. In particular, we are able to reproduce the theoretical behavior of the torque cutoff in Fourier space.

We also find various manifestations of nonlinear effects in our simulations even while the waves are formally linear. These include (1) the noticeable steepening of the density profile at larger values of M_p far from the planet; (2) the slight variation of the asymptotic value of the AMF with planetary mass in addition to the expected F_H∝M²_p dependence, causing the discrepancy with linear theory prediction at the level of several percent; and (3) the growth of power in high azimuthal wavenumber harmonics of the AMF in Fourier space.

By carefully studying the spatial distribution of the torque density dT_H(x)/dx in our simulations we discover an interesting "negative torque" phenomenon: dT_H(x)/dx changes sign at some radial separation from the planet (at |x₋| ≈ 3.2h in our simulations), which contradicts the analytical results of GT80. This effect can however be understood in the framework of the linear theory as shown in R. R. Rafikov & C. Petrovich (2012, in preparation). This feature of the torque distribution only has a weak effect on the total accumulated torque in our simulations.

We also carried out a detailed investigation of the effect of different numerical settings on our results in the linear regime. We explored the influence on wave properties of (1) different Riemann solvers with different accuracy, (2) spatial resolution, (3) different forms of the softened planetary potential, and (4) softening length of the planetary potential. We find the spatial distribution of the AMF and torque to be a more sensitive probe of the effects of various numerical parameters on the wave evolution than the distortions of the density wave profile. Based on this study we conclude that a very high resolution (ensuring low numerical viscosity), high order of accuracy, and an accurate prescription of planetary potential with small softening length are critical for accurately reproducing the key features of the wave evolution in linear regime.

We demonstrate that a low order of accuracy, low spatial resolution, and inaccurate potential with large softening length often employed in global disk–planet interaction studies can severely affect the fidelity of their results, especially in applications to planet migration and gap opening. Specifically, we conduct a test run with typical numerical parameters from recent literature. Comparing with analytical theory, the test run produces a lower amplitude of the density wave (by ∼20%), a lower final accumulated torque (by ∼23%), and a largely advanced starting point of wave dissipation (a factor of ∼2).

Most of our own calculations were carried out with third order of accuracy, resolution of 256/h, softening length r_s = h/32, and a potential that rapidly converges to Newtonian (as (r_s/ρ)⁴) at large separation ρ from the planet. This set of numerical parameters allows us to obtain excellent agreement with linear theory in Athena, but is not meant to be universal for other codes. However, the way in which we make the comparison with the theory and the agreement we achieve may serve as a standard framework for future code tests.

We also discover a numerical problem which is largely ignored in previous simulations of disk–planet interactions. To follow the fluid motion correctly, the time step in a simulation has to be smaller than the local dynamical timescale by a certain factor (∼100 in our case using Athena) in the entire domain, including the region where the fluid experiences the largest gravitational acceleration (e.g., the vicinity of the planet). In the context of numerical disk–planet studies, violation of this condition leads to inaccurate calculation of wave properties, lack of convergence on long timescales, and spurious repulsion of gas by a planet similar to gap opening. This timescale issue applies to disk–planet simulations in general, but it particularly affects the ones which use the orbital advection algorithm (FARGO), in which case a significant increase of the time step is usually allowed to speed up calculations.

In Paper II, we will continue our investigation of the disk–planet interaction by looking at the details of the density wave evolution in the nonlinear phase.

We are grateful to Jeremy Goodman, Takayuki Muto, Zhaohuan Zhu, and an anonymous referee for useful discussions and comments. The financial support of this work is provided by NSF grant AST-0908269 and Sloan Fellowship awarded to R.R.R.