NONLINEAR DYNAMICAL FRICTION IN A GASEOUS MEDIUM

A massive object in orbital motion suffers from orbital decay due to a negative torque caused by gravitational interaction with its own gravitationally induced wake created in the background medium. This process, commonly referred to as dynamical friction (DF), occurs in not only collisionless environments (e.g., Chandrasekhar 1943) but also collisional gaseous backgrounds (e.g., Ostriker 1999). The DF in a gaseous medium is of great importance in understanding the formation and evolution of planets, binary stars, supermassive black holes, etc. For instance, gravitational interaction between a protoplanet and its environmental disk causes the former to migrate toward a central star, naturally explaining the presence of "hot Jupiters" found from radial velocity surveys (Butler et al. 2006, and references therein). The migration also helps a planet grow faster in mass by providing an expanded gas-feeding zone at an enhanced accretion rate, which may overcome the failure of in situ core-accretion scenario in building a giant planet within a typical disk lifetime (see, e.g., Alibert et al. 2005). In the case of nuclear black holes in merging galaxies, they are expected to first experience the DF to form a binary and then coalesce into a supermassive black hole by emitting gravitational waves. Friction of nuclear black holes against the collisionless stellar background appears to be inefficient due to scattering and depletion of stars near the black holes, which is known as the "final-parsec problem" (see Milosavljević & Merritt 2003, and references therein). However, recent numerical N-body/smoothed particle hydrodynamics (SPH) simulations show that the gravitational drag from the gaseous background is sufficient to form a black hole binary in a relatively short time (∼1–10 Myr; e.g., Escala et al. 2004, 2005; Dotti et al. 2006, 2007; Mayer et al. 2007; Cuadra et al. 2009).

Thanks to a seminal paper of Ostriker (1999), DF in a gaseous medium is well understood as long as density wakes have small amplitudes. Earlier theoretical work by Dokuchaev (1964), Ruderman & Spiegel (1971), and Rephaeli & Salpeter (1980) considered density wakes in a steady state and found that the drag force vanishes for a subsonic perturber, while it becomes remarkably similar to the collisionless counterpart for supersonic cases. Using a time-dependent linear perturbation theory, on the other hand, Ostriker (1999) found that the gravitational drag force on a point-mass perturber with mass M_p moving at velocity V_p on a straight-line trajectory through a uniform gaseous medium with density ρ_∞ and sound speed a_∞ is given by

$$\begin{equation} F_{\rm lin}= \frac{4\pi \rho _\infty (GM_p)^2}{V_p^2} \left\lbrace \begin{array}{@{}ll} \frac{1}{2}\ln \left(\frac{1+\mathcal {M}}{1-\mathcal {M}}\right)-\mathcal {M}, & \mathcal {M}<1,\nonumber\\ \frac{1}{2}\ln \left(1-\frac{1}{\mathcal {M}^2}\right)+\ln \left(\frac{V_p t}{r_{\rm min}}\right), & \mathcal {M}>1, \end{array}\right.\nonumber\\ \end{equation}\vspace*{2pt} \tag{ 1 }$$

where $\mathcal {M}\equiv V_p/a_{\rm \infty }$ is the Mach number, t is the time elapsed after the introduction of the perturber, and r_min is the minimum radius introduced to avoid the singularity in the force evaluation. Equation (1) shows that the gaseous DF force becomes identical to the Chandrasekhar (1943) formula for the collisionless drag for $\mathcal {M}\gg 1$, and is, albeit small, nonzero even for subsonic perturbers because the time dependence breaks the symmetry in the density wakes (see also Just & Kegel 1990). Equation (1) has been applied to various astrophysical situations including orbital decay of compact objects in accretion disks (e.g., Narayan 2000; Karas & Subr 2001) and heating of an intracluster medium by supersonically moving galaxies in clusters (e.g., El-Zant et al. 2004; Faltenbacher et al. 2005; Kim et al. 2005; Kim 2007; Conroy & Ostriker 2008).

While the result of Ostriker (1999) is valid in a strict sense only for a linear-trajectory perturber in a uniform medium, it has proven to be applicable to more general cases. For example, Kim & Kim (2007) considered a circular-orbit perturber with orbital radius r_p in a uniform gaseous medium and found that Equation (1) is a reasonable approximation to the gaseous drag force on it, provided V_pt = 2r_p. Sánchez-Salcedo & Brandenburg (2001) numerically found that the orbital decay of a Plummer sphere with a softening radius r_s in a radially stratified medium is consistent with the prediction of Equation (1), if V_pt/r_min = (0.35r_p/r_s)^2.34. Also, Barausse (2007) showed that Equation (1) remains valid even for a perturber with relativistic speed if the relativistic correction factors are included. While Dotti et al. (2006) found that the orbital decay of black hole binaries took longer than the prediction of Equation (1) for a single perturber, the discrepancy between the numerical and analytical results can be reconciled, at least partly, by taking allowance for the fact that an object in a binary experiences not only a negative torque due to its own wake but also a positive torque from the companion wake. For an equal-mass binary, Kim et al. (2008) found that the positive torque is on average about 40% of the negative torque.

Since the results of Ostriker are based on the assumption that density wakes remain in the linear regime, the validity of Equation (1) for very massive perturbers has yet to be seen. The strength of gravitational perturbations due to a body with mass M_p can be measured by the dimensionless parameter

$$\begin{equation} \mathcal {A}=\frac{GM_p}{a_{\rm \infty }^2r_s}, \end{equation} \tag{ 2 }$$

which roughly corresponds to the perturbed density at a distance r_s from the perturber relative to the background density (e.g., Just & Kegel 1990; Ostriker 1999), and is equal to the Bondi radius r_B = GM_p/a²_∞ relative to r_s. For systems with $\mathcal {A}\gg 1$, the density wakes are clearly in the nonlinear regime and the linear perturbation analyses are likely to fail. Identifying r_s with the gravitational softening radius of a perturber (or, equivalently, its size), $\mathcal {A}$ is in the ranges of ∼0.1–1 for galaxies embedded in typical intracluster media, ∼10–100 for protoplanets in protostellar disks, and ∼10⁶–10⁸ for supermassive black holes near galaxy centers, suggesting that the wakes of massive compact objects can readily be nonlinear. Indeed, Escala et al. (2005) reported that the orbital decay time of supermassive black hole binaries with $\mathcal {A}\sim 10\hbox{--}100$ depends on M_p much less sensitively than the results of the linear theory, which may be caused primarily by the nonlinear effects.

In this paper, we investigate nonlinear DF of a very massive perturber in a gaseous medium using numerical hydrodynamic simulations. In order to isolate the effects of the perturber mass and its velocity on the DF force, we consider a perturber following a straight-line trajectory in a uniform gaseous medium, similar to that in Ostriker (1999). We model the perturber as a Plummer sphere that does not possess any solid surface and merely provides gravitational potential perturbations to the background medium that would otherwise remain static and uniform; to make contact with the results of the linear theory, we ignore the accretion of gas onto the perturber in the current work. Our primary objectives are to find the changes in distributions of density wakes with $\mathcal {A}$ and $\mathcal {M}$, and to quantify the resulting gravitational drag forces in comparison with the linear cases.

Nonlinear responses of a background to a massive perturber moving at a supersonic speed have been extensively studied in the context of the Bondi–Hoyle–Lyttleton (BHL) accretion (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952; see also review of Edgar 2004, and references therein). Unlike our models where mass accretion to a perturber is prohibited, however, the BHL accretion problem considered a perturber containing a defined surface through which gas is either accreted or reflected. It was Hunt (1971, 1979) who first solved for BHL accretion flows numerically, finding that the collisional nature of gas supports a bow shock in front of a supersonic perturber. Later studies found that the BHL accretion flows exhibit unstable behaviors such as flip-flop motions of the accretion shocks and vortex shedding when the condition of axisymmetry is relaxed (e.g.,Matsuda et al. 1987; Fryxell & Taam 1988; Taam & Fryxell 1988; Matsuda et al. 1989; Ruffert 1994; Foglizzo & Ruffert 1997, 1999; Foglizzo et al. 2005). While these numerical works on the BHL accretion explored temporal evolution and distribution of density wakes as well as nonlinear features in some great detail, they mainly concentrated on the gravitational focusing and resulting accretion rate of gas onto the perturber. Although some authors (e.g., Shima et al. 1985; Shankar et al. 1993; Ruffert 1996) presented values for aerodynamic and gravitational drag forces, only a limited range of $\mathcal {A}$ was covered. In this work, we run a number of numerical simulations by varying $\mathcal {A}$ and $\mathcal {M}$ systematically in order to quantify the dependences of the gravitational drag force on these parameters. A brief comparison between our results with those from the BHL accretion studies will be presented.

This paper is organized as follow: in Section 2, we describe numerical methods we employ for nonlinear simulations. In Section 3, as a code test we revisit the cases with a spatially extended, linear perturber with $\mathcal {A}=0.01$. We compare the resulting distributions of density and velocity wakes with those of analytical results and provide a way to handle the effect on the DF force of a softening radius which is necessary for numerical simulations. Evolution and quasi-steady distributions of fully nonlinear density wakes and the associated drag forces are presented in Section 4. Finally, in Section 5 we summarize our findings and discuss their astrophysical implications.

In this paper, we consider an unmagnetized, inviscid, non-radiating, non-self-gravitating gaseous medium and study its gravitational responses to a massive perturber moving along a straight-line trajectory using numerical simulations. The background gaseous medium is initially static and uniform with density ρ_∞ and adiabatic sound speed a_∞. We adopt an adiabatic equation of state with an index γ = 5/3 throughout simulations. The simulations are carried out on a two-dimensional (R, z)-plane in cylindrical symmetry, where R and z denote the distances from and along the axis of symmetry, respectively. We assume that the perturber exerts only gravity to the surrounding medium and does not possess any surface, so that neither accretion nor reflection of the gas is allowed. The perturber is modeled as a Plummer sphere with mass M_p that is moving at a constant speed V_p along the R = 0 axis toward the positive z-direction, with the gravitational potential

$$\begin{equation} \Phi _P(R,z,t)=-\frac{GM_p}{\big(R^2+r_s^2+[z-V_pt]^2\big)^{1/2}}, \;\;\;{\rm for} \;t\ge 0, \end{equation} \tag{ 3 }$$

where r_s is the softening radius.

We take r_s, a_∞, and t_cross = r_s/a_∞ as the units of length, velocity, and time, respectively, in our simulations. Then, our models are completely parameterized by $\mathcal {A}$ and $\mathcal {M}$. We run a total of 58 models with $\mathcal {A}$ varying from 0.01 to 600 and $\mathcal {M}$ in the range of 0.5–4.0. We solve the basic equations of ideal hydrodynamics using FLASH3 (Fryxell et al. 2000), an Eulerian hydrodynamics code that implements a direct piecewise parabolic method solver of Colella & Woodward (1984) for high-order spatial reconstruction. Although the FLASH3 code is capable of both uniform grid and adaptive mesh refinement calculations, we adopt the uniform grid method since the accurate evaluation of the drag force requires the whole computational domain to be well resolved. As we will show below, we find that it is necessary to have at least five zones per r_s to obtain converged results for the drag forces. Our largest grid models have 3072 × 12,288 zones in (R, z); we make sure that our computational domain is taken to be large enough to contain the whole density wake in a given model. The simulations are typically carried out until t/t_cross = 600 when most of the wakes are well resolved and reach a quasi-steady state.

Time-dependent linear perturbation theories for the DF drag force in a gaseous medium usually study the responses of gas to a low-mass, point-mass perturber corresponding to r_s = 0 in the Plummer potential (e.g., Just & Kegel 1990; Ostriker 1999), which requires to introduce the cutoff radius r_min in the linear force formula (Equation (1)). In numerical simulations, on the other hand, one needs to assign a nonzero value to the softening radius, which in turn makes it unnecessary to use the cutoff radius in the force evaluation. Since our goal is to compare the nonlinear drag force on a massive perturber with the linear prediction, we have to first find a proper relationship between r_s and r_min that makes the numerical and analytical results consistent with each other when A ≪ 1. Motivated by this consideration, in this section we briefly present the results of numerical simulations for a low-mass perturber with $\mathcal {A}=0.01$, and compare the resulting distributions of density and velocity wakes and the drag forces with those from the linear theories. This will also allow us to check the accuracy of our numerical experiments.

Regarding the linear wakes with which numerical results will be compared, it is worth mentioning that there are several analytical methods for finding solutions for the perturbed density and velocity fields. Just & Kegel (1990) utilized Fourier transform for the space variables and Laplace transform for the time variable, finding expressions both for the density and velocity wakes. Instead, Ostriker (1999) used a retarded Green's function technique and found an expression only for the density wake that is identical to the result of Just & Kegel (1990). We found that while the analytical formula for the density wake agrees well with our numerical results for a low-mass perturber, the simulated velocity field differs from the expression given by Just & Kegel (1990). In the Appendix, we revisit the time-dependent linear theory using Fourier transforms both for the space and time variables. As we will show below, our expressions (A19) and (A20) for the perturbed velocities are in good agreement with the numerical results, confirming that Equation (47) of Just & Kegel (1990) contains a typographical mistake.

Figure 1 shows as color-scale images the snapshots on the z–R plane of the perturbed density α = ρ/ρ_∞ − 1 (top), the parallel velocity v_z (middle), and the perpendicular velocity v_R (bottom) to the line of motion for a model with $\mathcal {A}=0.01$ and $\mathcal {M}=1.5$. Note that the coordinates are normalized by a_∞t. The perturber initially introduced at (z, R) = (0, 0) has moved to $(\mathcal {M}a_{\rm \infty }t, 0)$ at time t. The characteristic features of a supersonic wake consisting of a sonic sphere with radius a_∞t centered at the initial perturber location and a Mach cone bounded by Mach waves located at $R=-(\mathcal {M}^2-1)^{-1/2}(z-\mathcal {M}a_{\rm \infty }t)$ for $z>a_{\rm \infty }t/\mathcal {M}$ are apparent in the top panel. Also plotted as black solid contours are the results of the linear perturbation theory (Equations (A10), (A19), and (A20)), which are overall in good agreement with the numerical results. A careful inspection of Figure 1, however, reveals that the numerical results deviate slightly from the analytical ones especially near the sonic radius and the Mach waves. This is due to the fact that the perturber in our numerical models is modeled as an extended Plummer sphere rather than a point mass.

Zoom In Zoom Out Reset image size

One can semianalytically construct the density wake α_ext of an extended perturber by convolving the density wake α due to the corresponding point mass with the extended mass distribution ρ_ext of the perturber (e.g., Just & Kegel 1990; Furlanetto & Loeb 2002). For a Plummer sphere we use, the convolution theorem gives

$$\begin{equation} \alpha _{\rm ext}(\mbox{\boldmath {$x$}},t)= \frac{1}{M_p} \int \alpha (\mbox{\boldmath {$x$}}-{ {\bm x}^\prime },t)\rho _{\rm ext}(\mbox{\boldmath {$x^\prime $}},t) \;d^3\!{ {\bm x}^\prime }, \end{equation} \tag{ 4 }$$

where ρ_ext(x, t) = 3M_pr²_s(R² + r²_s + (z − V_pt)²)^−5/2/(4π). Figure 2 plots as solid lines the profiles of α_ext for $\mathcal {A}=0.01$ and $\mathcal {M}=1.5$ along the cuts at R/a_∞t = 0.20 and z/a_∞t = 0.92 marked as dotted lines in the top panel of Figure 1, which are in excellent agreement with the simulation outcomes (open circles). Compared with the point-mass results (dashed lines), the extended mass distribution tends to smear out the discontinuities at the boundary of the sonic sphere and the Mach cone. This makes sense since a perturbed density at one location is a superposition of sonic perturbations with various strengths launched by all the mass elements comprising the extended body.

Zoom In Zoom Out Reset image size

Allowing for the extended mass distribution, the gravitational drag force exerted on the perturber in the negative z-direction can be obtained by directly evaluating the integral

$$\begin{equation} F(t)=\int \!\!\!\!\int \frac{G\rho _{\rm ext}(\mbox{\boldmath {$x^\prime $}},t)\rho (\mbox{\boldmath {$x$}},t) (z-z^\prime)}{|\mbox{\boldmath {$x$}}-\mbox{\boldmath {$x^\prime $}}|^3}\ d^3\mbox{\boldmath {$x$}}\ d^3\mbox{\boldmath {$x^\prime $}}, \end{equation} \tag{ 5 }$$

where ρ(x, t) is the wake distribution. Although Equation (5) requires the $\mbox{\boldmath {$x^\prime $}}$-integration to be performed over the entire Plummer sphere, we empirically found that the drag force on the region with distance from the perturber larger than 10r_s has a negligible contribution to the total. Thus, in practice, we limit the integration to the region within 10r_s that contains about 98.5% of the total perturber mass. Figure 3 plots as open squares the numerical drag forces on a low-mass perturber with $\mathcal {A}=0.01$ but differing $\mathcal {M}$. The solid line corresponds to Equation (1) with

$$\begin{equation} r_\mathrm{min}=0.35\mathcal {M}^{0.6}r_s, \end{equation} \tag{ 6 }$$

which gives the best fit to the supersonic results of our adiabatic simulations. Equation (6) is our prescription for the cutoff radius when we compare the numerical drag forces with the analytic results. Note that Sánchez-Salcedo & Brandenburg (1999) suggested r_min = 2.25r_s based on their isothermal simulations. The difference between the two prescriptions may be due in part to using different equations of state and in part to low resolution (1 cell per r_s) in Sánchez-Salcedo & Brandenburg (1999).

Zoom In Zoom Out Reset image size

4.1. Wake Evolution

4.1.1. Supersonic Cases

We begin by describing the temporal evolution of our fiducial model with $\mathcal {A}=20$ and $\mathcal {M}=1.5$; the evolution of other supersonic models is qualitatively similar. Figure 4 illustrates the density and velocity structures of this model in a comoving frame with the perturber located at (s ≡ z − V_pt, R) = (0, 0). Only the region with |s|/r_s ⩽ 20 and 0 ⩽ R/r_s ⩽ 20 is shown. Unlike the linear cases with $\mathcal {A}\ll 1$ where the density wakes have too small amplitudes to launch shock waves, the perturber with $\mathcal {A}=20$ emits strong perturbations that quickly develop into a bow shock. The upstream gas moving toward the perturber along the symmetry axis is first accelerated by the gravity of the perturber, is shocked to subsonic speed, and then piles up near the center of the perturber. This creates a steep pressure gradient in between the perturber and the shock, tending to push the shock away from the perturber. On the other hand, the gas flowing above (not far away from) the symmetry axis is deflected toward the perturber even before entering the shock and decreases its speed after the shock. This gas thus has a longer time to be exposed to the gravity of the perturber as it moves toward the symmetry axis. The gravitational potential is so deep that the material arriving at the rear side of the perturber can be pulled back toward the perturber, creating a stagnation point just as in the BHL accretion flows (e.g., Matsuda et al. 1987; Fryxell et al. 1987). This produces a strong counterstream that moves into the upstream direction along the symmetry axis, as well as a primary vortex in the s–R plane (Figure 4(a)).¹ The counterstream combined with the pressure gradient in the front side pushes the shock front away from the perturber.

Zoom In Zoom Out Reset image size

The advance of the shock front in the upstream direction allows more time for the shocked gas to be affected by the gravity, strengthening the counterstream (as well as the primary vortex) and thus increasing the detached distance δ of the shock measured along the symmetry axis. At the same time, the center of the vortex with a lower density than the surrounding buoyantly rises toward the high-R regions. This decreases the ram pressure of the counterstream exerted on the shock front and the shock advance slows down. Figure 5 plots the time evolution of δ for some selected models, while Figure 6 traces the trajectory of the center of the primary vortex in our standard model ($\mathcal {A}=20$ and $\mathcal {M}=1.5$) on the s–R plane. For our standard model, δ keeps increasing to 13r_s at t/t_cross = 100. At this time, the vortex is located near at (s, R) ≃ (−0.2r_s, 5r_s) (Figures 4(b) and 6). The shock front soon overshoots a potential equilibrium position where the thermal pressure of the postshock gas supports it against the perturber gravity, and begins oscillating around the equilibrium position.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shock oscillation changes the velocity field in the postshock region, causing the vortex to move in the counterclockwise direction around its mean position (−0.5r_s, 5r_s) following the background flow. Kelvin's circulation theorem implies that the vortex becomes stronger as it moves toward the symmetry axis, amplifying the strength of the counterstream near the perturber. The velocity of the counterstream becomes largest when the vortex arrives closest to the perturber on its oscillation path (Figure 4(d)). The shock front that was moving toward the perturber is pushed by the strong counterstream, and reverses its motion. The vortex rises to high R as the shock moves away from the perturber, weakening the counterstream, and the oscillation cycle repeats quasi-periodically.

Figures 5 and 6 show that the amplitudes of the detached shock oscillations as well as the vortex movements grow secularly with time as the shock oscillation continues. This can be understood as follows. When the shock front is displaced from the equilibrium position away from the perturber (e.g., Figure 4(c)), the shocked subsonic gas can acquire an extended time to be gravitationally influenced, similarly to the early situation when the shock advances away from the perturber. With the stagnation point moving away from the perturber, more mass and momentum are added to the counterstream, amplifying the vortex oscillation. In addition, when the strong counterstream collides directly with the shocked subsonic flow near the symmetry axis, the material at the interface is injected upward in the lateral direction and then carried in the negative s-direction by the background flow (see, e.g., the region at s/r_s ∼ 12 and R/r_s ∼ 0–3 in Figure 4(e)). This produces small vortices near the interface that move downstream and merge with the primary vortex, again strengthening the latter. Consequently, the primary vortex is able to move closer to the symmetry axis and amplify the counterstream in the next cycle, making the shock oscillation overstable.

As the center of the primary vortex moves away from the symmetry axis during its last phase of the overstable oscillation (t/t_cross ∼ 465–492), it becomes less gravitationally bound and can thus be more easily influenced by the background flow. For all the supersonic models we run in this paper, we find that whenever the vortex rises above a half of the accretion radius, also known as the Hoyle–Lyttleton radius, defined by r_A ≡ 2GM_p/V²_p, it is swept away by the background flow in the negative s-direction (Figure 6). When the primary vortex is located outside 0.5r_A, the counterstream associated with the vortex is weak and occurs in the far rear side of the perturber, unable to pass through the center of the perturber. The counterstream is instead resisted by the flow moving downstream right across the perturber and pushed up in the lateral direction (Figure 4(g)). This develops another vortex with an opposite sense of rotation to the primary one, capable of pushing the latter away from the perturber (Figure 4(h)). With the primary vortex carried away downstream, the flow near the perturber reaches a quasi-steady state in which shocked subsonic gas moves almost parallel to the symmetry axis. Since the associated kinetic energy is much smaller than the thermal and gravitational energies of the gas, the density distribution around the perturber becomes nearly hydrostatic, as will be shown below.

4.1.2. Subsonic Cases

Figure 7 displays density and velocity structures of a nonlinear subsonic model with $\mathcal {A}=20$ and $\mathcal {M}=0.5$ in a comoving frame with the perturber. Snapshots at t/t_cross = 60, 150, and 600 are shown. Sound waves launched from the perturber at t = 0 propagate radially outward into the surrounding medium, forming a spherical causal region within which the medium is affected by the sonic perturbations. Since the perturber is spatially extended, however, the boundary of the causal region is not as sharp as in the case of a point mass, although the most dominant perturbations still come from the perturber center. Unlike supersonic cases, this model always involves subsonic flows and never produces a shock. Nevertheless, the overall flow pattern and late-time density structure near the perturber of this model is very similar to those in the postshock subsonic regions of supersonic models. First of all, the strong gravitational pull forms a counterstream and an associated vortex ring near the symmetry axis (Figure 7(a)). The counterstream moving in the upstream direction interacts with the incident flow (Figure 7(b)). The gas at the interface is pushed up toward the high-R regions and then carried downstream, creating small vortices with low density (Figure 7(c)).

Zoom In Zoom Out Reset image size

The primary vortex slowly rises in the R-direction due to buoyancy, and merges with the small vortices. Since this model does not contain a shock that would confine the region of influence and since the causal region keeps expanding at a sonic speed, the flows in the high-R regions are almost parallel to the symmetry axis. As a result, the primary vortex keeps rising as there is no momentum input in the background flow capable of pushing it back toward the symmetry axis. At the end of the run (t/t_cross = 600), the primary vortex in this model arrives at (s, R) = (3r_s, 28r_s), corresponding to 0.18r_A. It is uncertain whether the primary vortex will be carried downstream when it goes beyond r = 0.5r_A in a manner similar to supersonic cases. At any event, the density distribution close to the perturber is well described by the condition of hydrostatic equilibrium at late time (see Section 4.2).

4.2. Quasi-steady Density Wakes

Figure 8 illustrates changes in the quasi-steady density wakes with varying $\mathcal {A}$ on the s–R plane for supersonic models with fixed $\mathcal {M}=1.5$ at t/t_cross = 600. The left panels show large-scale views of the wakes at −1400 ⩽ s/r_s ⩽ 100 and 0 ⩽ R/r_s ⩽ 600, while the region near the perturber with −60 ⩽ s/r_s ⩽ 20 and 0 ⩽ R/r_s ⩽ 60 is enlarged in the right panels. In each panel, the perturber is located at s = R = 0, and the black line connecting the points (s, R)/r_s = (0, 0) and (−60, 54) marks the boundary of the Mach cone characteristics of the linear density wake for $\mathcal {M}=1.5$. When $\mathcal {A}\ll 1$, the sonic perturbations are too weak to produce a shock, and the high-density ridge of the wake follows the Mach cone fairly well, although it is broadened due to the extended mass distribution of the perturber. As $\mathcal {A}$ increases to unity, sonic perturbations even outside the Mach cone attain substantial amplitudes enough to induce a bow shock that is attached to the center of the perturber (within the resolution limit). Since the density wake is effectively shifted toward the perturber compared with the linear case and still located preferentially at the rear side of the perturber, the resulting drag force will be larger than the linear counterpart (see Section 4.4).

Zoom In Zoom Out Reset image size

As $\mathcal {A}$ increases further, the shocked material gathered around the perturber begins to build up a strong pressure barrier which the incident flow cannot easily penetrate. This naturally makes the shock detached. Figures 8(c) and (d) show that the postshock density distribution around the perturber with $\mathcal {A}=10$ or 20 is almost spherically symmetric, indicating that the kinetic energy of the gas there is much smaller than the thermal and gravitational potential energies. To check if this is indeed the case, we plot in Figure 9 the density profiles along the symmetry axis and the R-axis from the center of a massive perturber with $\mathcal {A}=20$; both supersonic and subsonic models with $\mathcal {M}=1.5$ and 0.5 at t/t_cross = 600 are presented. Also shown as a dotted line is the density distribution under the assumption of hydrostatic equilibrium

$$\begin{equation} \rho =\rho _0 \left\lbrace 1+\frac{(\gamma -1)\mathcal {A}a_{\rm \infty }^2}{a_0^2} \left[\frac{r_s}{\big(r^2+r_s^2\big)^{1/2}}-1\right] \right\rbrace ^{1/(\gamma -1)}, \end{equation} \tag{ 7 }$$

where r = (s² + R²)^1/2 denotes the distance from the perturber center, and ρ₀ and a₀ are the density and the adiabatic sound speed at r = 0, respectively.

Zoom In Zoom Out Reset image size

For both supersonic and subsonic models, the density distributions at r/r_s ≲ 10 along the z⁺-, z⁻-, and R-cuts are virtually identical to each other, and are in remarkable agreement with the predictions of Equation (7) with ρ₀/ρ_∞ = 40 and a₀/a_∞ = 3.8 for $\mathcal {M}=1.5$, and ρ₀/ρ_∞ = 52 and a₀/a_∞ = 3.8 for $\mathcal {M}=0.5$. The sharp drop-offs of the density in the supersonic model near r/r_s = 12 and 32 along the z⁺- and R-directions, respectively, are of course due to the bow shock. In the subsonic model, the presence of the primary vortex makes the local density decreased at R/r_s = 28. Note that the quantity a²/ρ^γ−1 measures the specific entropy and thus is conserved in an adiabatic flow without involving a shock, as is the case in the subsonic model. For the supersonic model, however, a²/ρ^γ−1 is increased from unity to 1.34 because of a shock jump. This corresponds to the shock Mach number of 2.3, which is larger than $\mathcal {M}=1.5$ because the flow is accelerated by the perturber even before experiencing the shock. In the regions with r/r_s ≳ 10, the density is overall larger along the z⁻- than z⁺-directions, providing non-vanishing drag forces. Nonetheless, the presence of hydrostatic cores in the density wakes of massive supersonic perturbers makes the drag forces smaller than the linear results.

While the density wake near a massive perturber with $\mathcal {A}\ge 1$ is quite different from the linear counterpart, the distant part of the wake is more or less the same. Figure 10 plots the exemplary profiles of the normalized perturbed density, $\alpha /\mathcal {A}$, along the cuts with R = 0 and R/r_s = 200 for models shown in Figure 8. Note that $\alpha /\mathcal {A}$ is nearly independent of $\mathcal {A}$ in regions far away from the perturber (e.g., s/r_s < −900 region in Figure 10(a) and s/r_s < −500 region in Figure 10(b)); in these regions, the gravitational potential perturbations are weaker by more than two orders of magnitudes than at the perturber location and thus locally in the linear regime. Even with low amplitudes of local perturbations, however, the region with −900 < s/r_s < −100 along the symmetry axis behind the perturber has the nonlinear density wake that deviates considerably from the linear case. This is because the gas flowing in this region was already affected by strong gravitational potential in the upstream region, and has a diverging velocity field that reduces the perturbed density. Small fluctuations of nonlinear density wakes apparent in Figures 8 and 10 near the R = 0 axis are thought of as arising from sonic perturbations induced by the primary vortex and its oscillations discussed in Section 4.1.

Zoom In Zoom Out Reset image size

4.3. Detached Shock Distance

We have shown in the previous subsection that a sufficiently massive supersonic perturber generates a density wake that is characterized by a bow shock standing ahead of the perturber and a surrounding hydrostatic envelope. Figure 5 shows that the quasi-steady value of the detached shock distance δ is larger for models with larger $\mathcal {A}$ or smaller $\mathcal {M}(>1)$. To quantify the dependences of δ upon $\mathcal {A}$ and $\mathcal {M}$, we introduce the nonlinearity parameter

$$\begin{equation} \mathcal {\eta }\equiv \frac{\mathcal {A}}{\mathcal {M}^2-1}, \end{equation} \tag{ 8 }$$

and plot in Figure 11 the normalized shock distance against $\mathcal {\eta }$. The various symbols give the mean values of δ temporally averaged over t/t_cross > 50 that ignores the initial wake-development phase. The standard deviations of δ are also indicated by error bars. The numerical results are remarkably well described by the two simple power laws: $\delta /r_s=2(\mathcal {\eta }/2)^{2.8}$ for $0.7\lesssim \mathcal {\eta }\lesssim 2$ and $\delta /r_s=\mathcal {\eta }$ for $\mathcal {\eta }>2$.

Zoom In Zoom Out Reset image size

The behavior of δ with η can be qualitatively understood as follows. For a very massive perturber, the kinetic energy of the incident flow along the symmetry axis is almost entirely converted to the thermal energy of the postshock flow that supports a hydrostatic envelope against the gravitational potential of the perturber. When the shock is strong $(\mathcal {M}\gg 1)$, the postshock thermal energy proportional to $a_{\rm \infty }^2\mathcal {M}^2$ balances the gravitational potential energy −GM_p/δ at the shock location, resulting in $\delta /r_s\propto \mathcal {A}/\mathcal {M}^2$. In the limit of $\mathcal {M}\rightarrow 1$, the flow does not in principle produce a shock, corresponding to δ → ∞. In practice, the gravitational acceleration is able to turn an incident nearly transonic gas into a weakly supersonic flow, but the shock that may form is located very far away from the perturber anyway. In view of the detached shock distance, η may be a better indicator of nonlinearity than $\mathcal {A}$; for fixed $\mathcal {A}$, the shape of a density wake due to a supersonic perturber becomes similar to the linear counterpart as $\mathcal {M}$ increases.

It is well known from hydrodynamic experiments and corresponding theories that a supersonic flow over a solid object with a blunt nose develops a detached bow shock when the nose angle is larger than the maximum angle allowed by the postshock flow (e.g., Liepmann & Roshko 1957, and references therein). The shocked gas becomes subsonic and slowly adjusts its velocity as it flows downstream to meet the boundary conditions at the surface of the object. In our simulations, the incoming flow toward a massive perturber recognizes the hydrostatic envelope as a spherical obstacle. Since the nose angle of a spherical body with respect to the incident flow is 90°, the shock must be detached.

Even though near-hydrostatic envelopes that form in our simulations are not entirely impenetrable, we want to measure their effective sizes as perturbing obstacles. This can be achieved by comparing our numerical results for the detached shock distances with those from non-gravitating hydrodynamic theories (or lab experiments). Assuming that a bow shock ahead of a spherical body with radius R_s has a spherical shape near the symmetry axis, Guy (1974) showed that the standoff distance of the shock is approximately given by

$$\begin{equation} \frac{\delta }{R_s} = \left[\frac{(\gamma -1)\mathcal {M}^2+2}{4(\mathcal {M}^2-1)} + 1\right]^{1/(2K)} -1, \end{equation} \tag{ 9 }$$

where

$$\begin{equation} \frac{1}{K(\mathcal {M})} = \frac{1}{2} \left[ 1 + \frac{2}{\gamma +1}\frac{1-\mu ^2}{\mu }\right] \left[ 2\mu + 1 + \frac{1}{\mathcal {M}^2}\right], \end{equation} \tag{ 10 }$$

with $\mu ^2\equiv ((\gamma -1)\mathcal {M}^2 +2)/ (2\gamma \mathcal {M}^2 -\gamma +1)$. Equation (9) has proven to explain the experimental data quite well (e.g., Heberle et al. 1950; Schwartz & Eckerman 1956; van Dyke & Milton 1958; see also Schreier 1982). Figure 12 plots as a solid line δ/R_s from Equation (9) with γ = 5/3 as a function of $\mathcal {M}$. Also plotted as various symbols are our numerical results for the ratio of δ to the BHL radius, r_BHL ≡ GM_p/(V²_p + a²_∞) for models with $\mathcal {\eta }>2$, averaged over t/t_cross > 50. Again, error bars indicate the standard deviations. A rough agreement between the two results suggests that the BHL radius can be a useful measure of the effective size of the hydrostatic sphere.

Zoom In Zoom Out Reset image size

4.4. Gravitational Drag Force

For a given wake distribution ρ(x, t) at time t of a perturber, it is straightforward to calculate the DF force on it by performing integration in Equation (5). Figure 13 plots temporal changes of the DF forces normalized by 4πρ_∞(GM_p/a_∞)² for models with differing $\mathcal {A}$, but with fixed $\mathcal {M}=1.5$, over the course of the wake evolution. The dotted line corresponding to the linear DF force (Equation (1)), with r_min given in Equation (6), closely follows the numerical results for $\mathcal {A}=0.01$, showing that the linear drag force increases logarithmically with time. The drag forces for nonlinear cases with high $\mathcal {A}$ also have a similar logarithmic time dependence, although they fluctuate for a while in response to the oscillations of primary vortices as well as detached bow shocks before a quasi-steady state is attained; the fluctuation amplitudes are typically ∼4–16%, with a smaller value corresponding to larger $\mathcal {A}$ and $\mathcal {M}$. Note that for $\mathcal {M}=1.5$, the normalized drag force decreases with increasing $\mathcal {A}$, indicating that the nonlinear effect makes the DF force smaller than the linear estimate. This is because a higher value of $\mathcal {A}$ implies a correspondingly larger detached shock distance, and a hydrostatic envelope with front–back symmetry near the perturber contributes a negligible amount to the net DF force.

Zoom In Zoom Out Reset image size

The dimensionless drag forces at t/t_cross = 600 when the wakes are in a quasi-steady state are given in Figure 14 for various models with different $\mathcal {A}$ and $\mathcal {M}$. For all the supersonic models, the DF force is a decreasing function of $\mathcal {A}$. The reduction of the normalized DF force is larger for models with $\mathcal {M}\sim 1$ than highly supersonic models. The Mach number corresponding to the maximum drag force shifts from unity to ∼1.5 as the wake becomes highly nonlinear. For subsonic models, on the other hand, the nonlinear drag forces for $\mathcal {A}=20$ and 50 show some fluctuations (represented by error bars) associated with slowly evolving vortices present in the wakes, but their respective time-averaged values are very close to the drag forces in the linear regime. In fact, the similarity between the linear and nonlinear drag forces on subsonic perturbers is expected since the linear density wakes also possess front–back symmetry in the vicinity of the perturber (Ostriker 1999). Since the subsonic DF forces are dominated by the far field where perturbations are weak regardless of the perturber mass, the linear results should be valid even for very massive perturbers.

Zoom In Zoom Out Reset image size

The gravitational drag force certainly depends on both $\mathcal {A}$ and $\mathcal {M}$, but the discussions given above suggest that it may be through the nonlinearity parameter $\mathcal {\eta }$. To check this, we plot in Figure 15 the ratio of the nonlinear DF force F to the linear prediction F_lin as a function of $\mathcal {\eta }$. Again, various symbols give temporal averages of F/F_lin over t/t_cross > 50, and their standard deviations are indicated by error bars. For $\mathcal {\eta }\lesssim 0.7$ with which a bow shock that barely forms is attached to a perturber, F/F_lin ≈ 1. When $\mathcal {\eta }$ is increased to ∼0.7–2, the shock becomes detached, but its standoff distance is not so large. In this case, most of the material in the wake is still located behind of, but closer to, the perturber in comparison with the linear wake (see, e.g., Figure 8(b)), resulting in the DF force slightly larger (by less than 20%) than the linear counterpart. In highly nonlinear cases with $\mathcal {\eta }>2$, however, the presence of a large hydrostatic envelope makes the drag force reduced considerably. For $100 > \mathcal {\eta }>2$, the numerical results are well fitted by

Zoom In Zoom Out Reset image size

$$\begin{equation} F=F_{\rm lin}\left(\frac{\mathcal {\eta }}{2}\right)^{-0.45}. \end{equation} \tag{ 11 }$$

Figure 15 also plots the gravitational drag force from the hydrodynamic simulations of Shima et al. (1985) for BHL accretion flows, which are consistent with our numerical results. We defer to Section 5 a more detailed discussion of our results in connection with the BHL flows.

To ascertain that the reduction of the DF force in highly nonlinear supersonic cases is really caused by the presence of spherically symmetric hydrostatic envelopes near the perturbers, we calculate the drag force by imposing a cutoff radius r_min such that only the regions in the wake with r > r_min participate in the force evaluation (Equation (5)). Figure 16 plots the resulting dependence of F upon r_min for the model with $\mathcal {A}=20$ and $\mathcal {M}=1.5$ at t/t_cross = 600. The vertical dashed line marks the location of the bow shock along the symmetry axis, while the dotted line indicates a slope of −1. Note that the drag force is independent of r_min for r_min < δ, clearly demonstrating that the hydrostatic sphere surrounding the perturber has a negligible contribution to the net drag force. Figure 16 shows F ∝ −ln(r_min/r_s) for large r_min, analogous to the linear cases (see Equation (1)). From the study of BHL accretion flows to large gravitating bodies, Shankar et al. (1993) similarly found that the drag force declines logarithmically as the size of the accretor increases.

Zoom In Zoom Out Reset image size

4.5. Resolution Dependence

Finally, we remark on the effect of numerical resolution on the DF force. Figure 17 plots the time evolution of the detached shock distance as well as the drag force on a perturber with $\mathcal {A}=10$ and $\mathcal {M}=1.5$ from models with different resolution: Δz/r_s = 0.1, 0.2, 0.4, and 0.8, where Δz is the grid spacing. In models with Δz/r_s = 0.4 and 0.8, the shock has a larger standoff distance without noticeable oscillations, and the DF force is correspondingly smaller than in models with higher resolution. Compared with the Δz/r_s = 0.2 model, the model with Δz/r_s = 0.1 presents the shock oscillations with higher amplitudes and arrives at quasi-steady equilibrium earlier. This resolution dependence of the shock oscillations was also noted by Matsuda et al. (1989) for adiabatic BHL flows onto an absorbing perturber. The resolution study shown in Figure 17 suggests that our numerical results for the DF forces are reliable as long as five or more grids per r_s are taken.

Zoom In Zoom Out Reset image size

DF of bodies orbiting in a gaseous medium is of great importance in various astronomical systems ranging from protoplanetary disks to galaxy clusters. In previous analytic studies of DF, it has been assumed that the mass of a moving object is small enough for the induced density wake to have low amplitudes and be thus in the linear regime. However, there are many astronomical situations such as in a merger of black holes near a galaxy center and migration of protoplanets, where a perturber is so massive that the induced wakes are well in the nonlinear regime. In this paper, we use numerical hydrodynamic simulations to explore nonlinear gravitational responses of the gas to, and the resulting drag force on, a perturber with mass M_p moving straight at velocity V_p in an initially static, uniform background with density ρ_∞, and adiabatic sound speed a_∞. The perturber is represented by a Plummer sphere with softening radius r_s. Unlike in the BHL problems, the perturber in our models does not possess a defined surface through which gas is accreted or reflected. Assuming an axial symmetry, we solve the basic equations on the (R, z)-plane, where R and z denote the directions perpendicular and parallel to the motion of the perturber, respectively. Our numerical models are completely characterized by two dimensionless parameters: $\mathcal {A}=GM_p/(a_{\rm \infty }^2r_s)$ and $\mathcal {M}=V_p/a_{\rm \infty }$. To study DF in various situations, we run as many as 58 models that differ in $\mathcal {A}$ and $\mathcal {M}$. Our standard models have five zones per r_s, but we also run models with different resolutions to ensure that the density and velocity wakes are fully resolved.

For supersonic models $(\mathcal {M}\,\,{>}\,\,1)$, we find that a massive perturber with $\mathcal {A}\,\,{\gg}\,\, 1$ produces a bow shock ahead of it through which the incident supersonic flow becomes subsonic. In the beginning, the postshock flow in the nonlinear cases develops transient features such as a primary vortex ring and an associated counterstream near the symmetry axis, which causes the shock to oscillate around its equilibrium position. The shock oscillation added by the counterstream is overstable, amplifying the amplitudes of its oscillation and the vortex movement in the (R, z)-plane (see Section 4.1). The vortex ring is eventually shed downstream from the perturber, leaving the nearly stationary bow shock and a quasi-hydrostatic envelope that surrounds the perturber. On the other hand, subsonic models with $\mathcal {A}\,\,{\gg}\,\, 1$ without involving a shock retain a primary vortex and many other small-scale structures until the end of runs. Nevertheless, strong gravity makes the density distribution near the perturber still very close to that under hydrostatic equilibrium.

By comparing the numerical results from various supersonic models with differing $\mathcal {A}$ and $\mathcal {M}$, we find that the simulation outcomes such as the detached shock distance and the gravitational DF force are very well quantified by a single parameter η defined in Equation (8). When $\mathcal {\eta }\lesssim 0.7$, the system is in the linear regime where a Mach cone is bounded by low-amplitude Mach waves, and the drag force is just the same as the analytic linear value F_lin of Ostriker (1999). When $\mathcal {\eta }$ is moderate at $0.7\lesssim \mathcal {\eta }\lesssim 2$, the Mach waves turn into a bow shock that is weakly detached, with the standoff distance $\delta \approx 2(\mathcal {\eta }/2)^{2.8}r_s$ from the perturber. In this case, the density wake is slightly shifted toward the perturber compared with the linear counterpart, causing the drag force to be larger than F_lin by a factor of less than 1.2. In the highly nonlinear regime with $100 > \mathcal {\eta }>2$, however, the detached shock distance behaves as δ = ηr_s, and the nonlinear drag force F is given by F/F_lin = (η/2)^−0.45. The reduction of the drag force compared with the linear value is because the density wake close to the perturber is in near-hydrostatic equilibrium and thus contributes a negligible amount to the net DF force. Since front–back symmetry notable for nonlinear wakes exists also in the linear wakes with $\mathcal {M}\,\,{<}\,\,1$, the nonlinear drag force on a massive subsonic perturber is similar to the linear prediction.

As mentioned in Section 1, our model simulations differ from those of the BHL accretion flows in terms of the boundary conditions. In our models, a perturber simply provides a gravitational potential for the background gas and does not hold any surface, while models for the BHL accretion considered a defined surface through which the gas is accreted or reflected. The different boundary conditions might lead to different evolutions and structures of wakes. For instance, the BHL accretion flows around a perturber present nonlinear features including a detached bow shock, vortex oscillations, and counterstreams in the forward/backward directions, just as in our models, if the perturber is not perfectly absorbing (Shima et al. 1985; Fryxell et al. 1987; Matsuda et al. 1989). When the perturber has a totally absorbing surface, on the other hand, the wakes are relatively quiescent, flowing nearly spherically into the perturber that absorbs the angular momentum carried by the accreting gas (Shima et al. 1985; Fryxell et al. 1987; Koide et al. 1991; Ruffert 1994). In the case of an extremely small perturber compared with the accretion radius, even an absorbing body produces violent features since it is not able to accept the whole accreting gas (Koide et al. 1991).

Despite these differences, the detached distances of bow shocks and the resulting DF forces appear to be not so sensitive to the boundary conditions adopted. Figure 18 plots the detached shock distances inferred from the results of two-dimensional axisymmetric simulations (Shima et al. 1985; Fryxell et al. 1987; Matsuda et al. 1989; Koide et al. 1991) and the tabulated values from three-dimensional nonaxisymmetric runs (Ruffert 1994) for the BHL accretion flows, by taking r_s equal to the radius of the perturber. The BHL results under the reflecting boundary conditions are larger than our results by a factor of only ∼1.2, while those under the absorbing boundary conditions are smaller by a factor of ∼0.4. Since the work on the BHL accretion mostly focused on mass accretion rates, only a few of them evaluated the gravitational drag forces. In Figure 15, we plot using star symbols the drag forces from the adiabatic runs with γ = 5/3 in Shima et al. (1985), which are the only data we acquired from the literature that allow a reliable comparison with our results.² Note that the drag force is larger on a purely absorbing perturber for which the detached shock distance is smaller. The DF forces from Shima et al. (1985) are roughly consistent with, and follow a similar trend to, our results. This is probably because even though density wakes are nowhere close to being hydrostatic in BHL accretion flows, they somehow maintain front–back symmetry that makes the net DF force smaller.

Zoom In Zoom Out Reset image size

While we have considered in this paper only an adiabatic gas with index γ = 5/3, we note that the DF force may depends sensitively on γ. Simulations of the BHL accretion flows reported that as γ decreases, the detached shock distance decreases and the density wake becomes shaped increasingly into a cone-like structure similar to the linear wake (Hunt 1979; Shima et al. 1985; Ruffert 1995, 1996). This is presumably because bow shocks in lower-γ models should be stronger in order to compensate for the diminished pressure in supporting the postshock gas against the perturber's gravity. Consequently, drag forces are larger for the gas with smaller γ (Shima et al. 1985).

Escala et al. (2005) carried out numerical simulations of a DF-induced merger of supermassive binary black holes due to an adiabatic gas with γ = 5/3. Their Figure 9 presents the effect of the black hole mass on the evolution of the binary separation. From this figure, we infer that the orbital decay time scales with the black hole mass as τ_decay ∝ M^−0.3_p in their simulations, which is somewhat shallower than the linear expectation of τ_decay ∝ M⁻¹_p (e.g., Ostriker 1999). On the other hand, our results predict τ_decay ∝ M^−0.55_p if the perturber is sufficiently massive, suggesting that the delayed orbital decay of supermassive black holes in Escala et al. (2005) is partly due to the nonlinear effect. Of course, there are many other factors that may change τ_decay. While we consider quite an ideal situation in which a perturber is moving straight in a static, uniform medium, the gas in Escala et al. (2005) is distributed in a rotating, radially stratified, self-gravitating disk and the black holes follow curvilinear, possibly eccentric, orbits. In addition, the Mach number of the black holes in their models is likely to vary during the decay toward the orbital center, which may also modify the decay time. Since these compact objects have a very small size, it is also an issue whether numerical models resolve detached bow shocks, which is crucial in evaluating the drag force accurately. We will discuss these in subsequent work.

We acknowledge a helpful report from an anonymous referee and Eve C. Ostriker for useful suggestions. The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. This work was supported by KICOS through the grant K20702020016-07E0200-01610 provided by MOST. Simulations were performed by using the supercomputing resource of the Korea Institute of Science and Technology Information through the grant KSC-2009-S02-0008.

Using the Fourier transform for the space variables and the Laplace transform for the time variable, Just & Kegel (1990) derived an analytic solution for the velocity field in the linear wake, but unfortunately the resulting expression (their Equation (47)) contains an typographical mistake. Here, we rederive the expression using Fourier transforms for both space and time variables. We consider a point-mass perturber with mass M_p on a straight-line trajectory through a uniform medium with density ρ_∞ and adiabatic sound speed a_∞. The basic equations of hydrodynamics for an inviscid, nonmagnetic, non-self-gravitating gas are

$$\begin{equation} \frac{\partial \rho }{\partial t}+\mbox{\boldmath {$\nabla $}}\cdot (\rho \mbox{\boldmath {$v$}})=0 \end{equation} \tag{ A1 }$$

and

$$\begin{equation} \frac{\partial \mbox{\boldmath {$v$}}}{\partial t}+\mbox{\boldmath {$v$}}\,{\bm \cdot}\, \mbox{\boldmath {$\nabla $}}\mbox{\boldmath {$v$}} = -\frac{1}{\rho }\mbox{\boldmath {$\nabla $}}P-\mbox{\boldmath {$\nabla $}}\Phi _{\rm ext}, \end{equation} \tag{ A2 }$$

where Φ_ext is the gravitational potential of the perturber with density ρ_ext, satisfying the Poisson equation $\mbox{\boldmath {$\nabla $}}^2\Phi _{\rm ext}=4\pi G\rho _{\rm ext}$. Assuming that the density and velocity fields induced by the perturber have small amplitudes, we linearize Equations (A1) and (A2) to obtain

$$\begin{equation} \frac{\partial \alpha }{\partial t}+\mbox{\boldmath {$\nabla $}}\,{\bm \cdot}\, \mbox{\boldmath {$v$}}=0 \end{equation} \tag{ A3 }$$

and

$$\begin{equation} \frac{\partial \mbox{\boldmath {$v$}}}{\partial t} + a_{\rm \infty }^2\mbox{\boldmath {$\nabla $}}\alpha = -\mbox{\boldmath {$\nabla $}}\Phi _{\rm ext}, \end{equation} \tag{ A4 }$$

where α ≡ ρ/ρ_∞ − 1 is the dimensionless perturbed density.

Taking Fourier transformation of Equations (A3) and (A4), we have

$$\begin{equation} \left(k^2-\frac{\omega ^2}{a_{\rm \infty }^2}\right)\widehat{\alpha } =\frac{4\pi G}{a_{\rm \infty }^2}\widehat{\rho }_{\rm ext} \end{equation} \tag{ A5 }$$

and

$$\begin{equation} \widehat{\mbox{\boldmath {$v$}}}=\frac{\omega \mbox{\boldmath {$k$}}}{k^2} \widehat{\alpha }, \end{equation} \tag{ A6 }$$

where symbols with and without hat denote a Fourier pair defined by

$$\begin{equation} \Psi (\mbox{\boldmath {$x$}},t)=\frac{1}{(2\pi)^4}\int \!\!\!\!\int \!\!\!\!\int \!\!\!\!\int \widehat{\Psi }(\mbox{\boldmath {$k$}},\omega)e^{i(\mbox{\boldmath {$k$}}\cdot \mbox{\boldmath {$x$}}-\omega t)} \,d^3\mbox{\boldmath {$k$}}\,d\omega \end{equation} \tag{ A7 }$$

and

$$\begin{equation} \widehat{\Psi }(\mbox{\boldmath {$k$}},\omega)=\int \!\!\!\!\int \!\!\!\!\int \!\!\!\!\int \Psi (\mbox{\boldmath {$x$}},t)e^{-i(\mbox{\boldmath {$k$}}\cdot \mbox{\boldmath {$x$}}-\omega t)}\,d^3\mbox{\boldmath {$x$}}\,dt, \end{equation} \tag{ A8 }$$

with Ψ representing any physical quantities.

For a perturber introduced at t = 0 and moving at a constant speed V_p along the positive z-direction, $\rho _{\rm ext}=M_p\,\delta (\mbox{\boldmath {$x$}}-V_pt\mbox{\boldmath {$\mathrm{e_z}$}})\,\mathcal {H}(t)$, where $\delta (\mbox{\boldmath {$x$}})$ is the Kronecker delta and $\mathcal {H}(t)$ is a Heaviside step function. This results in

$$\begin{equation} \widehat{\rho }_{\rm ext}= M_p\left[2\pi \delta (\omega -k_z V_p)+\frac{i}{\omega -k_z V_p}\right]. \end{equation} \tag{ A9 }$$

The second term in Equation (A9) is originated from the step function, which vanishes if a steady state is assumed (e.g., Furlanetto & Loeb 2002).

Substituting Equation (A9) into Equation (A5), and taking an inverse Fourier transform of the resulting expression for $\widehat{\alpha }$, one can obtain the solution of the perturbed density for finite-time perturbations. The detailed procedure is similar to that of the Fourier–Laplace transform method of Just & Kegel (1990), so that we simply write the result

$$\begin{equation} \alpha =\frac{GM_p}{a_{\rm \infty }^2}\frac{\rm\Upsilon_1}{(s^2+R^2(1-\mathcal {M}^2))^{1/2}}, \end{equation} \tag{ A10 }$$

where s ≡ z − V_pt and ϒ₁ = 2, 1, and 0 in region I ($R^2+z^2>a_{\rm \infty }^2t^2,\ R^2<-sz,\ s^2+R^2(1-\mathcal {M}^2)>0$, $\mathcal {M}>1$), region II (R² + z² < a²_∞t²), and any other region, respectively. Note that Equation (A10) is identical to Equation (10) of Ostriker (1999) based on the retarded Green's function technique.

Combining Equations (A5), (A6), and (A9) together in favor of $\widehat{\mbox{\boldmath {$v$}}}$ and taking its inverse Fourier transform, we have

$$\begin{eqnarray} \mbox{\boldmath {$v$}}&=&\frac{GM_p}{4\pi ^3a_{\rm \infty }^2}\int \!\!\!\!\int \!\!\!\!\int \!\!\!\!\int \left[\frac{2\pi \delta (\omega -k_z V_p)}{k^2-\omega ^2/a_{\rm \infty }^2}\right.\nonumber\\ &&\left.\;+\,\frac{i}{\big(k^2-\omega ^2/a_{\rm \infty }^2\big)(\omega -k_z V_p)}\right] \frac{\omega \mbox{\boldmath {$k$}}}{k^2} e^{i(\mbox{\boldmath {$k$}}\cdot \mbox{\boldmath {$x$}}-\omega t)}\,d^3\mbox{\boldmath {$k$}}\,d\omega.\nonumber\\ \end{eqnarray} \tag{ A11 }$$

The second term in the integrand of Equation (A11) has three simple poles on the real ω axis at ω₁ = k_zV_p and ω_2,3 = ±ka_∞. To ensure the proper analyticity, we must displace them by an infinitesimal amount into the upper or lower half ω-plane. Since our problem requires v = 0 everywhere for t < 0, the proper choice of the pole displacements should be ω₁ into the upper half ω-plane and ω_2,3 to the lower half ω-plane. For t < 0, then, a contour of integration consisting of the real ω axis and an infinitely large semicircle in the upper half ω-plane encloses only ω₁, and the calculus of residues guarantees that after the integration over ω, the first and second terms in the integrand of Equation (A11) become identical with each other with opposite sign, resulting in no velocity perturbations before the introduction of the perturber.

For t > 0, the contour that goes along an infinitely large semicircle in the lower half ω-plane instead contains the poles at ω_2,3. The remaining steps in the integration with respect to ω are quite straightforward to yield

$$\begin{eqnarray} \mbox{\boldmath {$v$}}&=&\frac{GM_p}{4\pi ^2a_{\rm \infty }}\int \!\!\!\!\int \!\!\!\!\int \left(\frac{e^{ika_{\rm \infty }t}-e^{-ik_zV_pt}}{k+k_z\mathcal {M}}\right.\nonumber\\ &&\left.-\, \frac{e^{-ika_{\rm \infty }t}-e^{-ik_zV_pt}}{k-k_z\mathcal {M}} \right) \frac{\mbox{\boldmath {$k$}}}{k^2} e^{i\mbox{\boldmath {$k$}}\cdot \mbox{\boldmath {$x$}}}\,d^3\mbox{\boldmath {$k$}}. \end{eqnarray} \tag{ A12 }$$

To carry out the integration in Equation (A12), it is convenient to use spherical coordinates in the k-space as

$$\begin{equation} \mbox{\boldmath {$k$}}=k (\sin \theta \cos (\varphi +\zeta), \sin \theta \sin (\varphi +\zeta), \cos \theta), \end{equation} \tag{ A13 }$$

which is oriented relative to the space coordinates

$$\begin{equation} \mbox{\boldmath {$x$}}=(R\cos \zeta,R\sin \zeta,z), \end{equation} \tag{ A14 }$$

satisfying $\mbox{\boldmath {$k$}}\cdot \mbox{\boldmath {$x$}}=k(z\cos \theta +R\sin \theta \cos \varphi)$. We then obtain

$$\begin{eqnarray} \mbox{\boldmath {$v$}}&=&\frac{GM_p}{4\pi ^2a_{\rm \infty }}\int \!\!\!\!\int \!\!\!\!\int \left [ \frac{e^{ik(z\cos \theta +a_{\rm \infty }t)}- e^{iks\cos \theta }}{1+\mathcal {M}\cos \theta }\right.\nonumber\\ &&\left. -\,\frac{e^{ik(z\cos \theta -a_{\rm \infty }t)} - e^{iks\cos \theta }}{1-\mathcal {M}\cos \theta } \right ] e^{ikR\sin \theta \cos \varphi }\nonumber \\ &&\times (\sin ^2\theta \cos (\varphi +\zeta)\mbox{\boldmath {$\mathrm{e_x}$}}+\sin ^2\theta \sin(\varphi +\zeta)\mbox{\boldmath {$\mathrm{e_y}$}}\nonumber\\ &&+\,\sin \theta \cos \theta \mbox{\boldmath {$\mathrm{e_z}$}}) \,\,dk\,d\theta \,d\varphi, \end{eqnarray} \tag{ A15 }$$

where $\mbox{\boldmath {$\mathrm{e_x}$}}$, $\mbox{\boldmath {$\mathrm{e_y}$}}$, and $\mbox{\boldmath {$\mathrm{e_z}$}}$ are the unit vectors in the x-, y-, and z-directions, respectively. Without loss of generality, we take ζ = 0, in which case $\mbox{\boldmath {$\mathrm{e_x}$}}$ and $\mbox{\boldmath {$\mathrm{e_y}$}}$ correspond to the unit vectors in the radial $(\mbox{\boldmath {$\mathrm{e_R}$}})$ and azimuthal $(\mbox{\boldmath {$\mathrm{e_\varphi }$}})$ directions, respectively. It then follows that the azimuthal velocity vanishes (i.e., v_φ = 0) since this component in the integrand is an odd function of φ.

For the velocity component in the direction of motion, we substitute θ = π − θ' to the second term inside the square bracket of Equation (A15) to obtain

$$\begin{eqnarray} v_z&=&\frac{GM_p}{\pi a_{\rm \infty }} \int ^\pi _0\!\!\!d\theta \!\int ^\infty _0\!\!\!dk\, \frac{\sin \theta \cos \theta J_0(kR\sin \theta)}{1+\mathcal {M}\cos \theta }\nonumber\\ &&\;\left\lbrace \cos [k(z\cos \theta +a_{\rm \infty }t)] -\cos (ks\cos \theta) \right\rbrace. \end{eqnarray} \tag{ A16 }$$

Using [G6.671.8],³ inserting the substitutions

$$\begin{equation} \sin \eta =\frac{(R^2+z^2)\cos \theta +za_{\rm \infty }t}{R\big(R^2+z^2-a_{\rm \infty }^2 t^2\big)^{1/2}} \end{equation} \tag{ A17 }$$

and

$$\begin{equation} \sin \eta =\frac{(R^2+s^2)^{1/2}}{R}\cos \theta, \end{equation} \tag{ A18 }$$

for the first and second terms of Equation (A16), respectively, and applying [G2.551.3] for the η-integrations, we obtain

$$\begin{eqnarray} v_z&=&{-}\frac{GM_p}{V_p} \bigg[ \frac{1}{(R^2+s^2)^{1/2}} -\frac{\rm\Upsilon_1}{(s^2+R^2(1-\mathcal {M}^2))^{1/2}} \nonumber\\ &&-\frac{\rm\Upsilon_2}{(R^2+z^2)^{1/2}} \bigg], \end{eqnarray} \tag{ A19 }$$

where $\rm\Upsilon_1=$2, 1, 0 and $\rm\Upsilon_2=1$, 0, 1 for region I ($R^2+z^2>a_{\rm \infty }^2t^2,\ R^2<-sz,\ s^2+R^2(1-\mathcal {M}^2)>0$, $\mathcal {M}>1$), region II (R² + z² < a²_∞t²), and any other region, respectively.

For the velocity component in the radial direction, the use of [G6.671.1] similarly leads to

$$\begin{eqnarray} v_R&=&\frac{GM_p}{RV_p} \bigg[ \frac{s}{(R^2+s^2)^{1/2}} -\frac{s\rm\Upsilon_1}{(s^2+R^2(1-\mathcal {M}^2))^{1/2}}\nonumber\\[-2pt] && -\frac{z\rm\Upsilon_2}{(R^2+z^2)^{1/2}} \bigg]. \end{eqnarray} \tag{ A20 }$$

The results of our time-dependent linear analyses for density and velocity fields are confirmed by direct numerical simulations for a low-mass perturber with $\mathcal {A}=0.01$ presented in Section 3.

In terms of our notation,⁴ Equation (47) of Just & Kegel (1990) can be rewritten as

$$\begin{eqnarray} (v_R, v_z)&=&\frac{GM_p}{RV_p} \bigg[ \frac{(s, -R)}{(R^2+s^2)^{1/2}} -\frac{(s, -R)\rm\Upsilon_1}{(s^2+R^2(1-\mathcal {M}^2))^{1/2}}\nonumber\\[-2pt] && -\frac{(z, -R)}{(R^2+z^2)^{1/2}} \bigg]. \end{eqnarray} \tag{ A21 }$$

Comparison of Equation (A21) with Equations (A19) and (A20) shows that Equation (47) of Just & Kegel (1990) takes ϒ₂ = 1 everywhere, which is clearly incorrect.