Gaseous drag on a gravitational perturber in Modified Newtonian Dynamics and the structure of the wake

We calculate the structure of a wake generated by, and the dynamical friction force on, a gravitational perturber travelling through a gaseous medium of uniform density and constant background acceleration g_ext, in the context of Modified Newtonian Dynamics (MOND). The wake is described as a linear superposition of two terms. The dominant part displays the same structure as the wake generated in Newtonian gravity scaled up by a factor mu^{-1}(g_ext/a_0), where a_{0} is the constant MOND acceleration and mu the interpolating function. The structure of the second term depends greatly on the angle between g_{ext} and and the velocity of the perturber. We evaluate the dynamical drag force numerically and compare our MOND results with the Newtonian case. We mention the relevance of our calculations to orbit evolution of globular clusters and satellites in a gaseous proto-galaxy. Potential differences in the X-ray emission of gravitational galactic wakes in MOND and in Newtonian gravity with a dark halo are highlighted.


INTRODUCTION
According to the Brownian theory, a 'macroscopic' particle moving in a fluid experiences a fluctuating force as a consequence of the graininess of the fluid, resulting in dynamical friction (DF) and diffusion in the velocity-space. Understanding the nature of the DF force experienced by a gravitational object that moves against a mass density background is of great importance for describing the evolution of gravitational systems and the exchange of angular and linear momentum between gravitational subsystems. There are two limits of interest in astrophysics: pure collisionless media and fully collisional gas. If the pertuber only interacts with the surrounding particles gravitationally, either the medium is collisionless or collisional, the DF force is the result of the formation of an overdensity wake behind the perturber.
The DF depends on the nature of the force between the perturber and the background particles. Departures from Newtonian gravity have been proposed as an alternative to dark matter in galaxies (e.g., Modified Newtonian Dynamics; Milgrom 1983). Ciotti & Binney (2004) found that in collisionless backgrounds, the characteristic timescale for DF in Modified Newtonian Dynamics (MOND) is much shorter than in the Newtonian equivalent system with dark mat-⋆ E-mail:jsanchez@astroscu.unam.mx ter. If confirmed, this result may have strong astrophysical implications and may be used to distinguish between dark matter and modified gravity. Using the MOND scaling found in Ciotti & Binney (2004), the globular clusters in Fornax dwarf spheroidal galaxy should have spiraled into the centre of Fornax in ∼ 0.5 Gyr, forming a visible galactic nucleus (Sánchez-Salcedo et al. 2006). Ciotti & Binney (2004) derived the MOND DF timescale in a restrictive plane-parallel geometry, which is a too specialized case, and in a non-standard statistical formulation. These caveats led Nipoti et al. (2008) to carry out fully N-body simulations of the dynamics of bars in MOND systems. Their simulations confirmed the scaling of the DF timescale predicted in Ciotti & Binney (2004). Still, it is not entirely clear how the background responds to the perturbation, and the description of the structure of the density wake induced by a small and massive perturber is lacking. Since MOND is a nonlinear theory, the gravitational acceleration induced by the perturber is affected by the presence of an external gravitational field. In particular, the DF force is expected to depend on the angle between the direction of the mean field and the velocity of the perturber. Our aim is to understand the nature of DF in MOND.
In this paper we consider gaseous media and study the structure of the density wake induced by a gravitational perturber in MOND. Gaseous DF has found many applications in astrophysics, ranging from protoplanet accretion all the way to the motion of galaxies in clusters (e.g., Ostriker 1999;Kim 2007;Conroy & Ostriker 2008). In a pure baryonic universe as in MOND theory, the relative role of gas becomes even more important because all the contribution to the mass is accounted for by the gas and the stars. Moreover, the complementary view of hydrodynamics may provide new insights into the differences and analogies between Newtonian and modified gravities. Although the details on the internal density structure of the induced wake may depend on whether the system is fully collisional or collisionless, the fluid limit is useful for understanding the main conceptual features introduced by the change of the gravity law.
The linear response of the gaseous medium to a gravitational perturber in Newtonian gravity is well-documented (Dokuchaev 1964;Ruderman & Spiegel 1971;Just & Kegel 1990;Ostriker 1999;Sánchez-Salcedo & Brandenburg 1999, 2001Kim & Kim 2007;Kim et al. 2008). In all these works, a minimum radius rmin in the Coulomb logarithm was introduced in order to regularise the gravitational potential of a point mass. Although the formulae were derived for rectilinear orbits in homogeneous and infinity media, simple 'local' extensions have been proven very successful in more realistic situations, e.g., smoothly decaying density backgrounds or when the perturber is moving on a circular orbit (Sánchez-Salcedo & Brandenburg 2001;Kim & Kim 2007).
The paper is organized as follows. In §2, we discuss the basic concepts on the ideal problem of a point particle moving at constant speed within a uniform gas in MOND (the Bondi-Hoyle problem). In §3, we outline the linear derivation of the basic equations for calculating the steady-state density wake generated by an extended body. In §4 we describe the structure of the resulting wake and evaluate the DF force exerted on the perturber. We then discuss some implications of our results in §5. Finally, we conclude in §6.

THE BONDY-HOYLE PROBLEM IN MOND
We consider a gravitational point particle at the origin of our coordinate system, surrounded by a gas whose velocity far from the particle is where c∞ is the sound speed of the gas at infinity and M is the Mach number. We are interested in the wake produced by the gravitational interaction with the perturber in the context of MOND. Bekenstein & Milgrom (1984) suggested a Lagrangian theory and also a nonlinear differential equation for the nonrelativistic MOND gravitational potential produced by a mass density distribution ρ: where a0 ∼ 10 −8 cm s −2 is an universal acceleration and µ(x) is a monotonic and continuous function with the property that µ(x) ≃ x for x ≪ 1 (deep MOND regime) and µ(x) ≃ 1 for x ≫ 1 (Newtonian regime). Only for very special configurations (one-dimensional symmetry -spherical, cylindrical or plane symmetric systems-or Kuzmin discs) the MOND acceleration g is related to the Newtonian acceleration, g N by the algebraic relation µ(|g|/a0)g = g N (Brada & Milgrom 1995). At variance with the Poisson equation, MOND is a non-linear theory. This implies that the gravitational potential generated by the perturber depends in a complex way on the external field g ext in which it is immersed. The external field is to be thought of as the field of an enveloping system in the absence of the perturber. For example, the perturber can be a galaxy in the external field of a cluster of galaxies, or a dwarf spheroidal galaxy in the field of its parent galaxy. In MOND, it is crucial to include the external field g ext to describe DF on a perturber. The case g ext = 0 is not very relevant when considering DF because it corresponds to a situation where all the surrounding gas is bound to the perturber, whereas DF concerns the interaction with unbound distant particles having large impact parameters. MOND permits to have unbound particles if g ext = 0. For these reasons, we will assume that g ext is a no-null vector field. Along this paper we will assume that the unperturbed gas density is constant all over the space and that the body moves in a uniform rectilinear orbit, which seems to be in contradiction with the inherent assumption that there is an external gravitational field. Note, however, that this approximation is amply used even in Newtonian dynamics where in most of the cases (e.g., a galaxy inside the galaxy cluster) the gravitational body is immersed in a stratified medium. The 'local' approximation, that is, estimating the drag force at the present location of the perturber as if the medium were homogeneous but taking appropriately the Coulomb logarithm, has been proven very successful in both collisionless and collisional systems (e.g., Sánchez-Salcedo & Brandenburg 2001, and references therein). So far we are not interested in including the complications of density gradients in the unperturbed surrounding medium.

The far-field equation: external dominated field
Far away enough from the object, the change in the external gravitational potential can be treated as a linear perturbation. This is valid as soon as the physical size of the perturber is supposed to be small as compared to the characteristic length-scale of the medium. Using the subscript 0 for referring to unperturbed quantities, the linear field equation is where Φ1 is the change in Φ produced by an increment ρ1 in ρ (Milgrom 1986). Here µ0 ≡ µ(gext/a0), L0 is the logarithmic derivative of µ0 (in the unperturbed system) and E0 represents a 3 × 3 matrix with elementsẼ0,i,j =ê0,iê0,j , whereê0 is the unit vector in the direction of g ext , which is r dependent. As Milgrom (1986) pointed out, Φ1 satisfies an equation analogous to the electrostatic field equation with an inhomogeneous and anisotropic dielectric tensor.
Since we are not interested in the distortions in the wake by tidal forces, we will assume that g ext is a constant vector all over the space. This assumption also facilitates a comparison with the more familiar Newtonian case. If so, µ0, L0 andê0 are also constant and the linear field equation for a constant external field can be written as (Milgrom 1986): It can be seen that the potential becomes Newtonian (but with a larger effective gravitational constant) and anisotropic. We must note that, although the above equation was derived for |∇Φ1| ≪ gext, it is also valid when gext ≫ a0, regardless the value of |∇Φ1|. In this limit µ0 ≈ 1 and L0 ≈ 0; the conventional Poisson equation is recovered.

The Bondi-Hoyle radius
In the classical problem of a slow Brownian particle in a fluid, the field particles are assumed to form a heat bath.
That is, they all stay close to thermal equilibrium, despite the presence of the Brownian particle and, hence, the equation of motion is solved by perturbation. In the case of a point gravitational perturber immersed in a perfect gaseous medium, the Bondi-Hoyle radius defines the region where the response of the gas is linear. In Newtonian dynamics the Bondi-Hoyle radius is rBH ≡ GM/(c 2 ∞ (1 + M 2 )). Streamlines whose impact parameter is less than 2rBH , will bend significantly and pass through a shock. Hence, within 2rBH it is not any longer a small perturbation. In order to regularise the gravitational potential of a point mass, one has to introduce a minimum radius in the formulae of the DF drag (rmin ≈ 2rBH ).
In the Appendix A, we estimate the Bondi-Hoyle radius for a point mass in MOND. It is shown that the MOND Bondi-Hoyle radius is larger than in Newtonian gravity by a factor between µ −1 0 and µ −1 0 (1 + L0) −1/2 , depending on the angle between the velocity of the particle and the external field. To get a sense of values of rBH in typical cases, consider a galaxy of 5×10 11 M⊙ orbiting supersonically M ≈ 1.5 in a cluster of galaxies with a sound speed of intracluster gas of ∼ 1500 km s −1 . The Bondi-Hoyle radius in this case is rBH µ −1 0 kpc. In a typical galactic cluster µ0 = 0.3-1 (see, e.g., fig. 7 in Sanders & McGaugh 2002), hence rBH 1-3 kpc, implying that, if the interaction with the intracluster gas is merely gravitational, the linear approximation is satisfactory for studying the gaseous wake even quite close to the galaxy. For extended perturbers with characteristic size much larger than 2rBH , the flow is essentially laminar at any location. In the remainder of the paper we will describe the perturbation on the gas using linear theory.

Modelling the perturber
Our aim is to study the large-scale gravitational perturbation induced by a small perturber travelling through a much larger system. We obtain that, beyond a certain distance from the perturber, |∇Φ1| ≪ gext. In other words, the farfield wake is expected to be in the external field-dominated regime. In order to highlight the differences between genuine MOND and Newtonian gravity, we restrict our considerations to situations in which the potential is dominated by the external field everywhere, i.e., we will consider an extended perturber of mass M and characteristic size rp, where the internal acceleration gint ≈ GMp/(µ0r 2 p ) ≪ gext 1 . The following density-potential pair, which corresponds to a "modified" Plummer model, is an exact solution of Eq. (4) where rp andrp ≡ (1 + L0) 1/2 rp are the characteristic radii. Note that for L0 = 0 and µ0 = 1 (Newtonian limit), it corresponds to the classical spherical Plummer model. Sometimes it is useful to express ρp in terms of the central density of the object, ρc, as follows where The central density can be written in terms of the central density for the spherical Plummer model in classical Newtonian gravity,ρc, as ρc =ρc/(1 + L0) 1/2 . The selection of a Plummer model was for analytical purposes. As long as the size of the system is much larger than the perturber, the structure of the wake and the drag force experienced by the perturber are not expected to be sensitive to the details of the potential close to the body.
In order to approach this problem analytically, we study the simplest case: an externally dominated perturber. In real life, there are some Galactic dwarf spheroidal galaxies that are known to be in this regime (e.g., Milgrom 1995;Sánchez-Salcedo & Hernandez 2007). It is likely that some subclusters and groups of galaxies, with low internal accelerations, embedded in a main massive galaxy cluster (such as the NGC 4911 group in the Coma Cluster) lie also in this regime (e.g., Sanders & McGaugh 2002).
In Eqs (5) and (6), the external field was taken along z and, therefore, in the same direction as the incident flow (see Eq. 1). In this section we will focus on this axisymmetric case. The derivation of the equations when g ext is perpendicular to v∞ is postponed up to §4.1.3. These two situations brackets a general case where the external field has an arbitrary angle with respect to the velocity of the flow at infinity.

Linear equations in an external dominated field
In the following, we give the linear derivation of the wake in a medium with unperturbed density ρ0 and adiabatic sound speed c∞, ignoring gas self-gravity and any magnetic fields. As stated in Eq. (1), it is assumed that the gravitational perturber is seated at the origin of our coordinate system and the gas velocity far from the perturber is v∞ = Mc∞ẑ.
In the axisymmetric case, v∞ and g ext are parallel. We are interested in the steady-state density enhancement D(r) = (ρ − ρ0)/ρ0 produced by the gravitational interaction with the perturber. The steady-state linearized basic dynamical equations for adiabatic perturbations ρ = ρ0 + ρ ′ and v = v∞ + v ′ are: and Our strategy is to eliminate v ′ everywhere. By substituting equation (9) in the divergence of equation (10), we obtain that D satisfies the differential equation where L is the linear differential operator The operator L arises frequently in fluid dynamics (e.g., Landau & Lifshitz 1959). In Eqs (9)-(12) we have not specified the law of gravity. In MOND, the Laplacian of the potential in Eq. (11) can be expressed in terms of ρp using the MOND field equation (4). The equation for D becomes From now on, it will be convenient to use the following dimensionless variables:x = x/rp,ŷ = y/rp and z = z/(1 + L0) 1/2 rp. The hat symbol over a certain variable χ will be used to denote that χ is written withx,ŷ andẑ as the arguments 2 , e.g.,D(x,ŷ,ẑ) = D(x, y, z). The secondorder derivative in Eq. (13) can be performed as soon as the potential is known. Evaluating the second-order derivative of the potential in Eq. (13) using the potential given in Eq. (6), and rearranging and grouping the terms, the so-lutionD can be expressed as a linear superposition of two contributionsD =D1 +D2. Each one satisfies the following differential equations: and wherer 2 =x 2 +ŷ 2 +ẑ 2 , 2 At this stage it is obvious that this transformation does not conserve mass in the sense that if ρ is a density field, then R ρdxdydz = Rρ dxdŷdẑ. and the operatorL iŝ with According to Eq. (18) We recall that the equation forD in the Newtonian case iŝ which is naturally recovered from the above equations just by taking µ0 = 1 and L0 = 0 (so that M eff = M). The com-ponentD1 obeys a differential equation similar to the Newtonian case and hence is identical to the wake induced by a perturber with mass densityρ1 in conventional Newtonian gravity, once G has been replaced by a larger effective value G/µ0. The profileρ1 corresponds to the classical (spherical) Plummer model in the transformed coordinates (x,ŷ,ẑ), multiplied by a form factor between 0.59 (deep-MOND limit) and 1 (Newtonian limit). In analogy to the Newtonian case (Eq. 21), a fictitious Newtonian perturber with a pseudo-density mass density distri-butionρ2 = T0ρcĝ2 generates a component identical toD2.
We refer toρ2 as pseudo-density because it may take negative values. Interestingly, the total mass associated with this distribution iŝ Our natural choice was to adopt the Cauchy principal value of the integral, which we denote by the symbol PV (e.g., Mathews & Walker 1970). Note thatρ2 decays more slowly with radius thanρ1. Hence, the pseudo-densityρ2 becomes larger thanρ1 at large radii. The importance of the contribution of D2 to the wake and to the gravitational drag is difficult to forsee without a quantitative study.

Formal solution
The solution of Eqs (14) and (15) can be obtained using the retarded Green's function, which may be derived following different paths (e.g., Just & Kegel 1990;Ostriker 1999).
In particular, the Green function can be found after Fourier transforming and imposing causality when choosing the contour of integration in the complex plane for M eff > 1 (e.g., Just & Kegel 1990;Furlanetto & Loeb 2002). In the steadystate, the perturbed density fieldsD1 andD2 arê where for i = 1, 2, with β 2 eff ≡ M 2 eff − 1 and  Figure 1 shows the integralÎ1, which is proportional to the perturbed densityD1, in theẑ-R plane, of a M eff = 0.8, 1.13, 1.5 body. In the deep MOND regime, they correspond to physical Mach numbers of 0.53, 1.25 and 1.9, respectively. So far we are only interested in the 'far-field' perturbed density, hence we will not delve into details regarding the near-field (within a few core radius from the perturber). A subsonic perturber generates a density distribution with contours of constant density corresponding to similar ellipses with eccentricity M eff . For supersonic motions, however, the region of perturbed density is confined within the rear Mach cone, dragged by its apex by the perturber. The surfaces of constant density within the wake correspond to hyperbolae in theẑ-R plane, with eccentricity e = M eff . This is expected because, as we show in Section 3.2, the equation for D1 has the same form as in the Newtonian case with M eff , once the density distribution of the perturber is rescaled by a factor (1 − T0)(1 + L0) −1/2 , and G is replaced by G/µ0. Using this analogy, we can take advantage of the analytical results of Ostriker (1999) to find the 'far-field' perturbed density: where ξ ′′ = 2 for supersonic perturbers and 1 in the subsonic regime. In Eq. (27) we used thatM1 =ρc Rĝ 1dxdŷdẑ = M/r 3 p . From the equation above, we see that the isodensity contours are z 2 + R 2 (1 − M 2 ) =const, i.e. ellipses or hyperbolae with eccentricity e = M. In the Appendix B we reconsider the perturbation as a time-dependent rather than a steady state problem. Figure 2 showsÎ2 for the same three effective Mach numbers as previously considered (M eff = 0.8, 1.13 and 1.5). In the subsonic regime,Î2, with a bipolar structure, is positive along the z axis and negative in the perpendicular plane. For supersonic perturbers, an overdense bump at the head of the perturber is generated. Interestingly,Î2 displays a drop in density along the surface of the Mach cone. This negative jump in the surface of the Mach cone is very remarkable at M eff = 1.13 and dilutes at larger effective Mach numbers. Our calculations show that at M eff = 1.75 there is no jump in the Mach cone surface. For perturbers with M eff > 1.75 the density jump in the Mach cone becomes positive. We see that, in general,Î2 may take values comparable toÎ1. We note, however, that the real densityD2 is related toÎ2 through a factor T0 ≤ 1/6 (see Eq. 24). Figure 3 contains the superposition (1 − T0)Î1 + T0Î2, which is proportional to D, in the deep MOND regime (i.e. L0 = 1 and T0 = 1/6). In the subsonic regime, the inclusion of the component D2 plays an important role. The resulting contours of isodensity can be fitted by ellipsoids defined by the equation R 2 + z 2 /q 2 =const, where q is the flattening parameter of the distribution in the physical R-z plane. According to our discussion in the preceding section, for M eff = 0.8 and L0 = 1 (i.e. M = 0.53), we know that the isodensity contours of D1 have q = 0.85. The superposition of components D1 plus D2 generates a density distribution with q = 1.23 (in the deep-MOND regime). In order to have an axis ratio of q = 0.85, the body should be moving at M = 0.85. For M = 0.2, the flattening parameter is 0.98 if only D1 is considered, and becomes 1.38 when D2 is also taken into account. As a general conclusion, the wake for a subsonic body is more flattened along the direction of motion in MOND than in Newtonian dynamics.

The component D2 in the axisymmetric case
In the supersonic deep-MOND case, the overdensity head is still visible when both components are added. The inclusion of D2, however, does not change significantly the structure of the wake within the Mach cone for M eff ≤ 1.5. At larger Mach number the contribution of D2 relative to D1 becomes less and less important. Since the structure of D1 is a scaled version of the wake generated in the Newtonian case, the difference between the structure of a wake generated by a small highly-supersonic perturber in MOND, would be likely too subtle to be distinguished from a perturber of fictitious mass (1 − T0)µ −1 0 in Newtonian gravity. Differences only appear in the vicinity of the perturber.

External field orthogonal to the velocity flow
Suppose now that the external field is along the x-axis. By denoting nowx = x/rp,ŷ = y/rp,ẑ = z/rp, and by reasoning entirely analogous to that leading to Eqs (14)-(20), one finds that when the external field is perpendicular to the velocity flow, M eff ,ĝ1 andĝ2 are given by and In this case, the dependence on L0 does not factorize so that one needs to calculate the integrals (25) for each pair (M eff , L0). We found numerically that D1 is very axisymmetric around the x axis; the effect of the gravitational dilation on D1 is small. For instance, when adopting L0 = 1, the angular variations of D1 are less than 2%, 2%, 8% and 13% for M = 0.4, M = 0.8, M = 1.25 and M = 1.8, respectively. D1 is almost undistinguishable (differences of ∼ 5 percent or less) from its counterpart in the axisymmetric case, and thus they are not shown. In particular, the opening angle of the Mach cone for supersonic perturbers in physical coordinates (z, R), is the same as in the axisymmetric case.
Unlike D1, D2 is expected to depart from axisymmetry about the z axis. We will focus again on the deep-MOND limit (L0 = 1). Figure 4 showsÎ2, which is proportional to D2, in three perpendicular planes: y = 0, x = 0 and z = 10rp, for M = 0.53 and M = 1.25. When the perturber moves subsonically, the structure of D2 in the y = 0 plane looks pretty much like in the axisymmetric case after a rotation of π/2 (compare Fig. 2 and Fig. 4). In the x = 0 plane, however, the density map is notoriously different. It clearly shows that the configurationD2 has not an axial symmetry around the z axis.
For subsonic perturbers, D2 may be able to change the flattening parameter of the wake. As an example, consider M = 0.8 and L0 = 1. In this situation, D1 has isodensity contours with q = 0.6. If the contribution of D2 is added, the composed wake exhibits q = 0.43 in the y = 0 plane, and q = 0.54 in the x = 0 plane (not shown). For a supersonic perturber, the overdensity regions in D2 are not located any longer at the head of the body; regions at the front exhibit a decrease in density. Zones with density depletion, that is D2 < 0, can be found at both downstream and upstream. The maximum density enhancement in D2 appears smaller, by a factor ∼ 4 at M = 1.25, than in the axisymmetric configuration. Figure 4 shows the complexity of the structure in the (x, y) plane (lower panels). The isodensity contours in that plane turn up to be elongated along the x-axis.

The gravitational drag on the body
Once the structure of the wake D is constructed, it is straightforward to evaluate the drag force exerted on the perturber by its wake. By symmetry, the only non-vanishing component lies along the z axis. In particular, the drag force in the axisymmetric case is If FDF,1 and FDF,2 denote the contribution to the drag force by components D1 and D2, respectively, we have: where We first estimate the relative contribution of FDF,2 as compared to FDF,1. For supersonic perturbers, F2 and F1 were calculated by carrying out the integration in Eq. (35) over all our box domain. Hence, theẑ-integral has upper/lower limitsẑ = ±25 and theR-integral has limitŝ R = 0, 25. This implies that our steady-state wakes have an extent along z of zmax ∼ 25rp to 35rp, depending on L0, and a Coulomb logarithm ln Λ ∼ ln zmax/2rp ∼ 3 (Sánchez-Salcedo & Brandenburg 1999).
For a purely steady-state density perturbation, the forward-backward symmetry in the subsonic case argues that zero net force acts on the perturber. However, Ostriker (1999) noticed that this conclusion is misleading because, although complete ellipsoids exert no net force on the perturber, there are always cut-off ones within the sonic sphere that exert a gravitational drag. Thus, we estimate the gravitational drag on a subsonic perturber by integrating Eq. (35) over the largest sonic sphere contained in our computational domain.
The ratio F2/F1 depends on the angle between v∞ and g ext . For illustration, Fig. 5 shows its behaviour as a function of M eff for supersonic bodies moving along the external gravitational field. In such a case, F2 is positive at large Mach numbers, implying that it contributes to drag the body. It becomes zero at M eff ≃ 1.75 and negative for lower values. The relative contribution of D2 becomes more important when approaching the transonic motion. At large Mach numbers, the ratio F2/F1 increases monotonically but very slowly.
In order to visualize the importance of including FDF,2, Figure 6 shows the total drag force exerted on the perturber as a function of M, together with FDF,1, in the deep-MOND limit (L0 = 1, T0 = 1/6). As anticipated, the sign of FDF,2 depends on the angle between v∞ and g ext . In contrast to the axisymmetric case, when v∞ and g ext are perpendicular, FDF,2 is positive for transonic Mach numbers and becomes negative (reduces friction) at high Mach numbers. The contribution to the drag by the component D2 is more important in the axisymmetric case but it is only noticeable (> 20%) at 1 < M < 1.5. Our numerical calculations show that FDF,1 scales with the size of the box domain as ∝ ln zmax, whereas FDF,2 increases somewhat slower with zmax. Therefore, the relative importance of FDF,2 is expected to be less for larger zmax. Now, we wish to compare the drag force in deep MOND (L0 = 1) and in Newtonian gravity (L0 = 0). Figure 6 also shows the Newtonian drag force experienced by a body of mass M , travelling at Mach number M when the wake has the same extent as in the MOND case, i.e. zmax = 25 √ 2rp. The drag force in MOND is a factor αµ −2 0 larger than in Newton, where α is a form factor that depends on the Mach number and on the angle between v∞ and g ext . In the axisymmetric case, α ≃ 0.6 at low Mach numbers (M 0.5), and becomes α ≃ 0.5 at 0.5 < M 1.0. For supersonic Mach numbers, α ≃ 0.4 at 1 < M 1.5 and increases monotonically with Mach number up to ≃ 0.8 at high Mach numbers. The explanation for α = 0.8 at high Mach numbers is covered in detail in the Appendix C. When the direction of motion is perpendicular to the external field, α is very similar to its value in the axisymmetric case at M < 1. At supersonic velocities, α ≃ 0.6 at 1 < M 1.5, and falls monotonically down to 0.5 at high Mach numbers. Hence, the drag force may vary with the angle between v∞ and g ext by as much as 50 percent.

Wakes by galaxies and falling groups in clusters
In this section, we discuss the implications of X-ray observations of the morphology of wakes in clusters of galaxies for modified gravities. As in Furlanetto & Loeb (2002), let us assume that the galaxy or group of galaxies moves supersonically through a constant density cluster core surrounded by a isothermal envelope: where ρ0 is the density in the core and rc is the core radius. The emitted surface brightness S = R (ǫ f f /4π)dl, where ǫ f f is the bremsstrahlung free-free volume emissivity, will be enhanced in the wake by: where Σ is the column density, and χ ≈ 4 if the gas is isothermal, or ≈ 5 if it varies adiabatically (Furlanetto & Loeb 2002). Maps of X-ray for a M = 1.25 body with v∞ perpendicular to the external gravitational field are shown in Fig. 7. The structure of the emission in the wake is very similar to that formed by a "compact" gravitational perturber under Newton gravity. The projected X-ray emission along line-of-sights perpendicular to the direction of motion is roughly independent of the location in the wake, except near the edges of the cone, and is given by As a consequence, in searching for the wake, one may expect an abrupt jump in surface brightness at the edge of its cone. We derived the gravitational wake induced by an extended body with a Plummer profile having a fast density decay at large radii (ρ ∝ r −5 ) in an attempt to model the baryonic mass of a certain bound object. Wakes in the cold dark matter (CDM) scenario are expected to be different than in purely baryonic MOND because the dark matter component in the halo of galaxies decays as r −2 , much more slowly than the baryonic mass density (≤ r −3 ). To illustrate this, Figure 7 also shows the X-ray emission generated by a pseudo-isothermal perturber with core radius rp in the Newtonian case. We see that the X-ray emission in the standard CDM scenario is more cuspy. A MOND wake can be distinguished from a CDM wake by detecting a sharp X-ray enhancement along the Mach cone. Furlanetto & Loeb (2002) made a detailed analysis of the wake morphology in CDM models for collisionless and fully collisional (fluid) dark matter (FDM) in the supersonic case. They found that, because in the collisional case the dark halo is truncated by ram pressure stripping, the X-ray emission of the wake is rather flat, similar to that we find in MOND. Due to the quantitative similarity between the wake in MOND and in FDM, many of the observational suggestions raised by Furlanetto & Loeb Figure 5. Ratio between F 2 and F 1 , as defined in Eq. (35), versus the effective Mach number, for the axisymmetric case. At M eff < 1, the ratio depends on the adopted value for L 0 and hence is not shown.
to distinguish between FDM and CDM can be used to distinguish between dark matter or MOND. We can repeat the reasoning of Furlanetto & Loeb (2002) and argue that the observations of the wake of the elliptical galaxy NGC 1404 in the core of the Fornax group marginally support CDM against MOND, but the evidence is very weak. In the last decade, this type of observations has improved considerably (e.g., Drake et al. 2000;Neumann et al. 2001;Machacek et al. 2005Machacek et al. , 2007Sun et al. 2006). Still, it is difficult to draw some firm conclusions because of the difficulty to isolate the structure of the gravitational wake from the hydrodynamical wake, that is the mass in the wake stripped from the own galaxy by ram pressure. We must warn that MOND and FDM predict the same structure of the gravitational wake past a galaxy but, in many other astrophysical aspects, they must give different predictions since they are not equivalent.
In principle, observations of the tails of subclusters of galaxies are a potential route to distinguish betweeen collisionless CDM and MOND. However, it is a well established issue that MOND still requires dark matter at cluster scales (e.g., The & White 1988;Sanders 1999). The inclusion of an isothermal dark matter component in MOND erases somewhat the abovementioned differences between the MOND wake and the CDM wake. The observed displacement between the X-ray peaks and the associated mass distribution, as derived from lensing data in the Bullet Cluster, basically rules out FDM (Markevitch et al. 2004) but not necessarily MOND (Angus et al. 2007).

DF timescale in a spherical system
In §4.2 we derived the DF force experienced by a body of mass M travelling on a rectilinear orbit through a homogeneous fluid medium in deep-MOND and found that: where α depends on the Mach number and on the angle that makes the velocity of the perturber and the direction of the external field. How does it compare with the DF force in a collisionless medium? The formula for the collisionless MOND DF force was derived by Ciotti & Binney (2004). They show that it is similar to the Newtonian case but replacing G → Ga0/gext plus an extra factor of √ 2: where with σ the velocity dispersion of the Maxwellian distribution of velocities of background particles (e.g., Binney & Tremaine 1987;Sánchez-Salcedo et al. 2006). By comparing Eqs (38) and (39), we can generalize the conclusion of Ostriker (1999) but for MOND: since the functional form of the gaseous DF drag is much more sharply peaked near M = 1, the drag force near M = 1 is larger in a gaseous medium than in a collisionless medium with the same density and σ = c∞. At M ≈ 1, using α ≈ 0.5 (see §4.2) there is factor of ∼ 3 difference in the MOND force between the two cases. For M < 1, the drag force is generally larger in a collisionless medium than in a gaseous medium, because pressure forces in a gaseous medium create a much more symmetric density perturbation in the background. For many problems of astrophysical interest, it is convenient to have the DF inspiraling timescale for a body that is initially on a circular orbit. For circular orbits, the above formula is correct as long as the maximum impact parameter in the Coulomb logarithm is taken as ∼ 2Rp, where Rp is the instantaneous orbital radius of the perturber (Sánchez-Salcedo & Brandenburg 2001;Kim & Kim 2007). Consider a massive body embedded in the gaseous outer spherical envelope of a galaxy of mass MG, with circular speed vc = (GMGa0) 1/4 . Suppose that the gas is isothermal and in hydrostatic equilibrium in the gravitational potential of the parent galaxy. In a typical galaxy, the outer parts are expected to be in the deep-MOND regime, hence gext = √ GMGa0/r and µ0 = (GMG/a0) 1/2 r −1 . For the sake of clarity, let us assume that the gas has the virial temperature: c∞ ≃ vc/ √ 2. Imposing hydrostatic equilibrium, the unperturbed gas in the envelope will pursue the following profile ρg = ρs " r rs where ρs is the gas density at the radius of reference rs. By substituting the values for µ0 and ρg into Eq. (38), and equating the rate of decrease of angular momentum d(M vcr)/dt to the torque rFDF , we find the deep-MOND evolution of the orbital decay of a massive perturber's nearcircular orbit in a gaseous isothermal spherical distribution where rinit is the orbital radius at t = tinit and with µs = µ0(rs). Consider now its equivalent Newtonian system (ENS) that is, the Newtonian system (dark matter plus gas) in which the baryonic gas has exactly the same density distribution as in the MOND system. The density distribution of dark matter in the ENS is: In the ENS, the dark matter component added to have the same "dynamics" also contributes to the DF experienced by the body. The frictional force in the ENS is The factor 0.428 arises because the dark matter has been considered collisionless (see, e.g., eq. [7-25] of Binney & Tremaine 1987). Equating again angular momentum loss with the torque, we find τENS, the time that the body takes to reduce its distance a factor e in the ENS: where Rg ≡ ρ dm /ρg =const. The ratio between the DF timescale in MOND and in ENS is: τM τENS = 4.6(1 + 0.428Rg )µ 2 init .
Here µinit ≡ µ0(rinit) and we used that α ≃ 0.5 for M = 1.4 (see §4.2). In terms of Rt,init defined as the ratio between the mass in dark matter and the baryonic mass (stars plus gas) inside rinit, we have τM τENS = 4.6(1 + 0.428Rg ) In dwarf protogalaxies, before the gas is turned into stars, Rg ≃ Rt, and for typical values of dark matter contents in these systems (Rt ≃ 20), the gaseous DF timescale in MOND is 10 times shorter than in the ENS. Condensed objects that form early (e.g., globular clusters) could spiral into the centre of their host galaxy more rapidly than would be predicted in the standard CDM haloes.

SUMMARY
After a discussion about the Bondi-Hoyle problem, we have calculated the gravitational density wake of a perturber moving through a uniform gaseous medium in MOND. A hydrodynamical treatment provides useful insight into the problem of DF. In order to describe the response of the farfield medium to a small perturber, we focused on the case when the perturbation is dominated by a constant external field. The analytical knowledge of the wake structure in this simple case is useful to interpret fully non-linear simulations of the MOND dynamics of realistic systems. The structure of the wake depends on the angle between the velocity vector and the external gravitational field. We have considered two cases: the velocity of the perturber being either parallel or orthogonal to the external field. For an intermediate case, the structure of the wake will lay somewhere between the extremes described here. The density wake is decomposed into two linear contributions. The morphology of the dominant contribution is a scaled version of the Newtonian wake, Figure 7. Normalized X-ray surface brightness maps for a 'compact' body moving at M = 1.25 in the deep-MOND nonaxisymetric case for a line-of-sight along the y-axis (upper panel) and along the x-axis (middle panel), along with the X-ray map for an 'extended' pseudo-isothermal body in Newtonian gravity (bottom). The background field lays along the x-axis. The observer's line of sight makes an angle π/2 with respect the velocity of the perturber.
being the MOND perturbed density a factor µ −2 0 larger than in Newtonian gravity. The second contribution to the wake depends greatly on the angle which the direction of motion makes with the external field.
The MOND DF force on the perturber induced by its own wake, as a function of the Mach number has been also derived and compared with the drag force in the Newtonian case. It is important to know the dependence of the drag force strength on Mach number to study the circularization of orbits. Our results confirm earlier analyses suggesting that the DF force is higher in MOND than in Newtonian gravity. The DF drag is larger when the motion of the perturber is along the external field, especially at Mach numbers between 1-1.5. In the context of the sinking satellite problem, the recent claim that the existence of an extended globular cluster population in Fornax is problematic for MOND gains strength (Sánchez-Salcedo et al. 2006;Nipoti et al. 2008). We find that MOND predicts an even faster orbital decay, especially when the satellite arrives to the halo of a galaxy before the gas has turned into stars.
In MOND we show that a supersonic perturber generates a wake in the gas with a well-defined Mach cone in which the the surface density increases substantially in a narrow region and then flattens. Because CDM haloes extend to larger radii than the baryonic mass, the wakes in the two models have significant morphological differences.
Inherent to the Chandrasekhar treatment of DF is the assumption that the perturber interacts with unbound particles. In MOND, this requirement immediately demands the inclusion of the external field. Once the background gravitational field is included, two-body orbits of distant particles are unbound, the interaction of the perturber with the medium is quasi-Newtonian, and the effects of distant encounters can be simply added. This topic is closely related to studies of relaxation processes in self-gravitating media. In a forthcoming paper we will present a standard derivation of the DF drag in collisionless systems in MOND.