Issue 
A&A
Volume 589, May 2016



Article Number  A10  
Number of page(s)  17  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201527629  
Published online  05 April 2016 
Dynamical friction for supersonic motion in a homogeneous gaseous medium
Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: daniel.thun@unituebingen.de; rolf.kuiper@unituebingen.de; schmidtstudent@unituebingen.de; wilhelm.kley@unituebingen.de
Received: 23 October 2015
Accepted: 27 January 2016
Context. The supersonic motion of gravitating objects through a gaseous ambient medium constitutes a classical problem in theoretical astrophysics. Its application covers a broad range of objects and scales from planetesimals, planets, and all kind of stars up to galaxies and black holes. In particular, the dynamical friction caused by the wake that forms behind the object plays an important role for the dynamics of the system. To calculate the dynamical friction for a particular system, standard formulae based on linear theory are often used.
Aims. It is our goal to check the general validity of these formulae and provide suitable expressions for the dynamical friction acting on the moving object, based on the basic physical parameters of the problem: first, the mass, radius, and velocity of the perturber; second, the gas mass density, soundspeed, and adiabatic index of the gaseous medium; and finally, the size of the forming wake.
Methods. We perform dedicated sequences of highresolution numerical studies of rigid bodies moving supersonically through a homogeneous ambient medium and calculate the total drag acting on the object, which is the sum of gravitational and hydrodynamical drag. We study cases without gravity with purely hydrodynamical drag, as well as gravitating objects. In various numerical experiments, we determine the drag force acting on the moving body and its dependence on the basic physical parameters of the problem, as given above. From the final equilibrium state of the simulations, for gravitating objects we compute the dynamical friction by direct numerical integration of the gravitational pull acting on the embedded object.
Results. The numerical experiments confirm the known scaling laws for the dependence of the dynamical friction on the basic physical parameters as derived in earlier semianalytical studies. As a new important result we find that the shock’s standoff distance is revealed as the minimum spatial interaction scale of dynamical friction. Below this radius, the gas settles into a hydrostatic state, which – owing to its spherical symmetry – causes no net gravitational pull onto the moving body. Finally, we derive an analytic estimate for the standoff distance that can easily be used when calculating the dynamical friction force.
Key words: gravitation / hydrodynamics / shock waves / waves / methods: numerical
© ESO, 2016
1. Introduction
The motion of gravitating objects through a gaseous ambient medium is a classical problem in theoretical astrophysics. Of special interest is the force acting on the object as it determines typical evolution or equilibration time scales. This force is caused by two main mechanisms, purely hydrodynamic drag and dynamical friction. Even without gravity, the flow of the gas will be perturbed by the body because it acts as an obstacle and blocks the gas in front of it. This causes an exchange of momentum between the embedded body and the surroundings, leading to a hydrodynamic drag force on the object. On Earth, for example, this gives rise to the air resistance of all flying objects, a drag force that tends to oppose the motion of the body. For a moving massive body in an astrophysical setting, its gravitational attraction onto the surrounding leads to the formation of a wake of higherthanaverage density behind the moving object, and, as a result, to a gravitational pull of the dense region onto the body. This additional force is again directed opposite to the motion of the body and leads to a slow down. This force is usually called gravitational drag force or dynamical friction (Chandrasekhar 1943; Ostriker 1999). Hence, the total drag force acting on a gravitating body is the sum of the hydrodynamic drag and dynamical friction.
Its astrophysical application covers a wide range of fields, as pointed out recently by Lee & Stahler (2014). In the common envelope phase, the decay time of the secondary companion is driven by dynamical friction (Ricker & Taam 2008, 2012), as is the important survival rate of planets around evolved stars (Villaver & Livio 2009). Further applications include the settling of massive stars in molecular clouds (Chavarría et al. 2010), the drag on a star in the accretion flow around black holes (Narayan 2000), the coalescence of black holes (Armitage & Natarajan 2005), and the migration of planetesimals (Grishin & Perets 2015). In the framework of planet evolution, dynamical friction plays a role in the change of planetary inclination owing to the interaction with the disk (Rein 2012; Teyssandier et al. 2013). The dynamical interactions of evolved stars that move through the interstellar medium (ISM) determine the observational properties of the dust in the envelope (Slavin et al. 2004; van Marle et al. 2011; Meyer et al. 2014a,b, 2015).
This incomplete list of applications shows that the phenomenon of dynamical friction is widespread and of general importance. However, there is still no satisfying derivation of the force on the object despite many years of research (see e.g. Lee & Stahler 2014). In this study, we tackle this important problem of dynamical friction via direct numerical modelling. We perform hydrodynamical simulations of a gravitating body moving supersonically through a homogeneous medium and extract the dynamical friction on the object. We analyse in detail the scaling of the drag with important physical parameters of the problem such as Mach number and mass of the object and others, and compare this to existing formulae for dynamical friction.
The paper is organized as follows. In Sect. 2, we summarize the status on formulae for dynamical friction. In Sect. 3, we describe the hydrodynamics equations and the physical and numerical setup of the simulations. In Sect. 4, we present comparison simulations to laboratory experiments for nongravitating objects moving at supersonic speed through air. In Sect. 5, we present the simulations of gravitating objects for various physical parameters and discuss the implications on the dynamical friction. In Sect. 6, we present a new formula for the drag force. In Sect. 7, we give a brief summary of the study.
2. Analytical estimates of dynamical friction
A first account of dynamical friction was given by Chandrasekhar (1943) who calculated the dynamical friction of a star embedded in a stellar cluster. His results were later adapted to the motion of stars through the ISM or galaxies through the intergalactic medium. Here, the ambient gas is typically treated as a static hydrodynamic continuum which is traversed by a moving body of mass M and the velocity V_{∞}. Using the impulse approximation, where the moving mass perturbs the surroundings only for a finite time, Ruderman & Spiegel (1971) derived a general formula for the dynamical friction force F_{DF} acting on the mass M: (1)In Eq. (1), ρ_{∞} denotes the ambient constant density, G is the gravitational constant, and C_{A} stands for the Coulomb logarithm (2)where s_{min} and s_{max} correspond to the minimum and maximum impact parameter of the interaction. A very similar expression for the drag on a star moving in a zero temperature medium was derived by Bondi & Hoyle (1944) who analysed the change in momentum by the gas that is collected in a very thin line behind the object. Additionally, they derived relations for the gas accretion rate onto the body.
To apply Eq. (1)to various physical problems, the relevant length scales s_{min} and s_{max} need to be known. Typically, for s_{max} the maximum extension of the medium is assumed. For the minimum length, s_{min}, the situation is not so clear and no general accepted recipe exists. Often, either the physical radius R of the object is chosen or, when the body becomes very small, s_{min} ≈ max(R,R_{A}). Here R_{A} stands for the accretion radius (3)Recently, Cantó et al. (2011) derived an approximation for s_{min} using ballistic orbit theory; however, it does not account for the pressure. Obviously, the question of the correct minimum radius of the object is still not clarified and remains ambiguous.
While Eq. (1)was derived for a pressureless medium (dust), some corrections have to be applied if pressure effects play a role. In the supersonic case – with V_{∞}>c_{∞}, where c_{∞} denotes the soundspeed of the unperturbed medium – a bow shock forms in front of the object where additional energy can be dissipated. Using the linearized fluid equations Ruderman & Spiegel (1971) and Rephaeli & Salpeter (1980) derive a correction to the above equation that depends on the Mach number of the object.
In an important work, Ostriker (1999) extended the above analysis and studied the timedependent dynamical friction force acting on a massive object embedded into a gaseous medium in the sub and supersonic regime. Her analysis resulted in an expression for the coefficient C_{A}, very similar to that of Rephaeli & Salpeter (1980), which includes a timedependent maximum radius, s_{max} ~ V_{∞}t, where t is the time the object has interacted with the gas. Her analysis is based on linear perturbation theory and is useful for studying timedependent phenomena. The expression for the supersonic case reads (4)where ℳ = V_{∞}/c_{∞} denotes the Mach number of the problem. Obviously, if we set s_{max} = V_{∞}t, then Eqs. (2)and (4)agree with each other in the limit of large ℳ because pressure effects become less important for highly supersonic flows. Here, we are not interested in the timedependent process, but focus on the final stationary state instead.
In addition to these analytic estimates there have been many numerical studies of moving gravitating objects embedded in a gas. The first to address this issue numerically was Hunt (1971) who used a special shock fitting method to calculate the flow around the body and the accretion rate onto it. Later, Shima et al. (1985) were the first to use modern fluid dynamical methods on this problem and performed axisymmetric simulations in spherical polar coordinates. In their now classic paper they studied the drag force and the accretion rate onto the object as a function of the velocity and found very rough qualitative agreement with Eqs. (1)and (2), when using for s_{min} the inner radius of the computational domain. Most of the subsequent simulations dealt with open inner boundary conditions and focused on the mass accretion rate onto the object. An introductory summary to this BondiHoyleLittleton accretion process is given in the review article by Edgar (2004). In this paper we do not study accretion onto the object and only study rigid bodies. The main focus lies on the accurate computation of the dynamical friction and – as a prerequisite – on the determination of the minimum integration radius s_{min}.
The validity of the formula by Ostriker has been demonstrated by SánchezSalcedo & Brandenburg (2001) who showed agreement for moderate Mach numbers in nonlinear, isothermal simulations. However, they used a Plummertype potential with a smoothing length much larger than R_{A}. Later, Kim & Kim (2009) reanalysed dynamical friction for extended bodies through numerical simulations. However, they did not consider objects with a rigid surface (like stars), but used again a Plummertype potential with a smoothing length, and hence their results may be more applicable to galaxies moving in the intergalactic medium. They compare their results to the formula by Ostriker (1999) and provide a new fitting formula involving R_{A}. They used an axisymmetric cylindrical coordinate system, which suffers in terms of limited resolution close to the centre.
Starting with Ruffert (1994) and Ruffert & Arnett (1994) there have been many simulations considering the threedimensional (3D) flow around an accreting object. An interesting feature discovered in these first studies is the onset of unstable flow when breaking the axial symmetry, an effect that is most pronounced in planar twodimensional (2D) flows that are nonphysical, however. A comprehensive summary of the 2D and 3D simulations that have been carried out has been given by Foglizzo et al. (2005). They point out that the origin of the instability, i.e. under what conditions it is expected and its effect on the drag force, is still not understood. In this paper we avoid the question of instability and focus on purely axisymmetric flows.
It is useful to compare the drag formula Eq. (1) to standard hydrodynamic drag formula. Dimensional analyses (Landau & Lifshitz 1966) yields the hydrodynamical drag force acting on an object with cross section A: (5)The dimensionless drag coefficient C_{D} contains the details of how the body and the medium physically interact with each other, e.g. through surface effects and the dependence on laminar vs. turbulent flows. If we use the accretion radius to set the interaction cross section , the ratio of the hydrodynamical drag and the dynamical friction is given by (6)In this study, we determine the dependence of dynamical friction acting on a rigid, gravitating body moving supersonically through a homogeneous medium on the basic physical parameters of the system. First, we demonstrate that the dynamical friction scales like Eq. (1) and we derive an analytic relation for the as yet undetermined minimum length scale s_{min} of the problem, and finally give a convenient expression for the dynamical friction in general.
3. Physics and numerics
3.1. Problem setup
We model gravitating, spherical objects with a given mass M and radius R moving with supersonic speed V_{∞} through a homogeneous medium that is characterized by its mass density ρ_{∞}, pressure p_{∞}, and temperature T_{∞} or soundspeed c_{∞}. As an initial condition, the object is placed instantaneously into the ambient medium and the evolution of the gaseous system is followed using timedependent hydrodynamical simulations. After the system reaches a steady state we determine the gravitational pull of the gas onto the object (the dynamical friction) as well as the hydrodynamic drag.
3.2. Equations
We study the motion of the ideal gas by solving the Euler equations Here, ρ denotes the gas mass density, u the velocity, p the gas pressure, and e = e_{kin} + e_{th} the total (kinetic and thermal) energy density of the gas. The acceleration due to external forces a_{ext} is given by the gravitational attraction of the moving object (10)where r is the distance from the centre of the object to the position under consideration. We close the equations of hydrodynamics with the ideal gas equation of state: (11)Gas pressure, density, and temperature are related via (12)where μ is the mean molecular weight, k_{B} the Boltzmann constant, and m_{H} the mass of the hydrogen atom. The soundspeed of the gas is given by (13)In our simulations we use a constant value of the adiabatic exponent γ. We study the influence of different γvalues on the dynamical friction in dedicated parameter series, see Sect. 5.7.
From the given values of the density ρ_{∞} and pressure p_{∞} of the initially homogeneous medium, its adiabatic index γ, and the speed V_{∞} of the moving object, the associated Mach number of the system is given as (14)
3.3. Numerics
Simulations are carried out using the opensource code PLUTO (Mignone et al. 2007), version 4. Part of the simulations were carried out with a new inhouse developed CUDA version for usage on GPUs. We use a RungeKutta time stepping scheme of second order with a secondorder reconstruction of states in space, a van Leer limiter in the reconstruction step, and a HartenLaxvan Leer solver for the Riemann problem.
Equations (7)to (9)are solved in the comoving frame of the object, i.e. the object is at rest in the modelling frame and the initial velocity of the surrounding homogeneous gas corresponds to the physical speed of the solid body. We use a 2D spherical grid (r,θ) assuming axial symmetry. The object is fixed at the origin of the coordinate system and the gas flows into the negative zdirection. The solid surface of the moving spherical object is represented as boundary conditions at the radial inner boundary of the computational domain. Here, we use reflecting boundary conditions, corresponding to a nonaccreting object. The computational domain extends in the radial direction from the radius R of the object up to R_{domain} = 100or1000 R. The symmetry axis of the domain is aligned with the trajectory of the object. In the polar direction, the computational domain extends from 0 to π using a grid spacing uniform in angle. To obtain high resolution at the area of interest around the object and to ensure an approximately quadratic grid spacing in the radial and the polar direction of each grid cell, we use a logarithmic grid spacing in the radial direction.
At the outer radial boundary R_{domain}, we set the boundary conditions according to the flow of the surrounding gas. In the upper hemisphere (for θ ∈ [0,π/ 2]), we implemented an inflow boundary condition, i.e. all ghost cells are set to the unperturbed values of the surrounding medium. In the lower hemisphere (for θ ∈ [π/ 2,π]), we implemented a zerogradient boundary condition so that the outflowing material can leave the computational domain without reflections. As already mentioned, we perform simulations in axial symmetry and therefore at θ = 0 (positive zaxis) and at θ = π (negative zaxis) we make use of axisymmetric boundary conditions.
Simulations with a size of the body of R = 0.1 R_{⊙} are performed on a grid consisting of 700 × 480 grid cells, simulations with a smaller size of the body of R = 0.01 R_{⊙} are performed on a grid consisting of 860 × 400 grid cells. We present a corresponding convergence study in detail in Appendix C.
3.4. Assumptions and simplifications
The numerical experiments are performed to determine the dynamical friction in a general astrophysical context. Hence, we do not take into account effects that depend on a specific system under investigation such as radiative heating and cooling. The flow is assumed to be inviscid. Selfgravity of the gas is not taken into account. Effects of magnetic fields are not investigated. The longevolution problem of a timedependent slowdown of the moving body by the acting dynamical friction is not included in this study (see SánchezSalcedo & Brandenburg 2001 for a discussion of this issue). The formula for dynamical friction is derived within the comoving frame of the moving body. The object is treated as a rigid body where no accretion through its surface is allowed. The medium is considered to behave like an ideal gas.
3.5. Calculation of the drag force
The total drag force acting on an object moving through an ambient medium can be calculated from the momentum balance in the final equilibrium state (Landau & Lifshitz 1966). In our case, the object is moving along the zdirection and we have to consider the zcomponent of Eq. (8), which reads in equilibrium (15)where a_{z} denotes the zcomponent of the gravitational acceleration due to the moving object. Integrating over the whole volume of the computational domain we can convert the first two parts into surface integrals and obtain (16)where df is the surface element and dV the volume element, and the subscripts in and out refer to the inner and outer boundary of the computational domain. The first term of Eq. (16) represents the total hydrodynamical force acting on the surface of the body, which is the sum of a momentum transport, , through the body’s surface (e.g. for porous objects or open, accreting bodies) and a pressure force, , acting directly on the body. The second term is the dynamical friction force on the body, F_{DF}, which is obtained by integrating over the whole volume of the domain. The sum of these two contributions must be balanced by the corresponding momentum transport and pressure terms at the outer boundary, and , respectively. In our case, the first surface integral in Eq. (16) is taken at the surface of the moving object – here the radius, R, of the spherical body – and the last surface integral at the outer boundary of the domain, R_{domain}. In the case of an impermeable rigid body, , and we obtain the final force balance as (17)We calculate all these forces for our simulations and evaluate the importance of the different contributions. Equation (17) shows that the drag acting on an object can either be evaluated at the inner boundary (plus gravity) or solely at the outer boundary.
4. Comparison to laboratory experiments
In computational astrophysics, it is only rarely possible to test numerical algorithms against laboratory experiments. However, in the case of the problem on hand – a sphere moving supersonically through a gaseous medium – experimental data is available for nongravitating moving bodies. To check the validity of our numerical ansatz and prove the accuracy of the code, including setup, boundary conditions, and the calculation of the drag force, we perform comparison simulations of a body moving supersonically through a homogeneous gaseous medium (air), which can be compared to existing data from laboratory experiments, namely by van Dyke (1982) and Billig (1967).
Fig. 1 Visualization of a comparison simulation result of a sphere embedded in air flow (γ = 1.4) at a Mach number ℳ = 1.53. The homogeneous gas arrives from the top and flows in the negative zdirection. Shown is the density distribution for the final equilibrium state in the central part of the computational domain around the sphere. Upper panel: numerical results, black lines denote isodensity contours. Bottom panel: laboratory data from van Dyke (1982, Fig. 266). 

Open with DEXTER 
4.1. Morphology
The blackandwhite photograph in van Dyke (1982), Fig. 266, shows the result of a laboratory experiment of a spherical nongravitating body moving with supersonic velocity through air. We present the right half of the original image in Fig. 1, bottom panel. The morphology of the system is characterized by a standing shock front ahead of the moving sphere, a clearly visible shock boundary between the unperturbed and the perturbed gaseous medium, and a lowdensity region past the moving object.
We model the same experiment numerically within our framework described in the previous section. Here, we switch off the gravitational force of the moving body and use an adiabatic exponent for air of γ = 1.4. The body is moving into the positive zdirection; actually, simulations are performed in the comoving frame of the body, hence, the gas flow is initialized into the negative zdirection. The homogeneous gas is initialized with a constant density ρ_{∞}. The Mach number of the flow is set to ℳ = V_{∞}/c_{∞} = 1.53. In the whole computational domain the initial condition is given by the unperturbed flow. At the start of the simulation the rigid sphere is embedded and perturbs the flow. We run the model until an equilibrium state has been reached.
In Fig. 1, upper panel, we show the morphology of the gas density around the sphere in the numerical experiment. Clearly visible is the bow shock in front of the sphere where the gas flow changes from supersonic to subsonic velocities, the shock front dividing the perturbed from the unperturbed gas, and the lowdensity region behind the moving body. As expected, the gas density reaches a maximum directly in front of the sphere. A visual comparison with the experimental data from van Dyke (1982), Fig. 266, as presented in Fig. 1, shows excellent agreement between the numerical experiment and the laboratory experiment in terms of their morphological characteristics. Furthermore, the magnitude of the density jump at the shock front from the numerical simulation agrees closely with the wellknown analytical RankineHugoniot jump conditions. In contrast to the axisymmetric and inviscid numerical result, the laboratory experiment shows that the lowdensity region past the moving object is subject to weak turbulence that cannot be captured in the idealized simulation setting.
The following comparison checks for the shock front physics and its dependence on the Mach number.
4.2. Standoff distance of the shock front
As we point out below, the shock’s standoff distance, R_{SO}, plays a crucial role in determining the dynamical friction on a gravitating moving body. Here, we measure R_{SO} from the centre of the moving spherical object. To test the dependence of the shock front on the problem’s Mach number, ℳ, we compute a sequence of models for the same setup as in the previous section but with a variety of different inflow velocities V_{∞}. Figure 2 shows the resulting standoff distance of the shock, measured from the centre of the sphere, as a function of ℳ.
Fig. 2 Quantitative comparison of the shock’s standoff distance as function of Mach number. The solid black line denotes the relationship derived from laboratory experiments as given in the main text. Small red dots represent results from numerical experiments for the same setup. Blue dots denote early numerical results by van Dyke & Gordon (1959). 

Open with DEXTER 
The R_{SO} decreases with increasing ℳ and approaches asymptotically a value of R_{SO} → 1.14 in the limit of highly supersonic flow. This outcome of the numerical experiments can be directly and quantitatively compared to the laboratory experiments by Billig (1967), using a derived fit to the laboratory data given by R_{SO}/R ≃ 1 + 0.143exp(3.24/ℳ), where R denotes the geometrical radius of the moving sphere. The results from the numerical and laboratory experiments are quantitatively in very good agreement, as shown in Fig. 2. They are consistent with each other at both extremes, weak shocks (ℳ = 1.53) and strong shocks (ℳ → ∞). In the regime of moderate shocks (ℳ ~ 2...4), the numerical experiments show slightly larger standoff distances than the laboratory experiments. Additionally, in the numerical framework, the position of the shock front is associated with a certain grid cell, more specifically its centre; no further interpolation has been applied. Hence, the visible steps in the numerical data visualize the finite spatial resolution of the numerical experiments.
4.3. Hydrodynamical drag force
Fig. 3 Force contributions in the equilibrium state according to Eq. (17) for the nongravitating, rigid sphere in air, given as a function of the Mach number. denotes the pressure force measured at the surface of the sphere, the pressure force measured at the outer boundary of the computational domain, and the momentum transport through the outer boundary. The sum of these three forces is denoted by F_{tot}. The curve describes a parabola with F ∝ ℳ^{2}. 

Open with DEXTER 
In Fig. 3 we show the contributions of the different forces acting on the gas and the sphere. Shown are the pressure force acting on the body, , together with the pressure and momentum flux at the outer boundary. As is expected from Eq. (17), the sum of all three force contributions, F_{tot}, adds up to zero. The error for the total force increases slightly for higher Mach numbers, but it is always below 2% for all simulations. In this case the hydrodynamic drag acting on the body, given by , is negative which indicates a force opposing the direction of motion of the body. The standard drag force, as given in Eq. (5), depends quadratically on the velocity V_{∞} of the object. Our numerical results clearly confirm this expected quadratic scaling with V_{∞}, see added solid line in Fig. 3.
5. Simulations of dynamical friction
Fig. 4 Density and velocity structure of the final quasistationary state around a gravitating object moving with supersonic motion through a gaseous homogeneous medium. The body is moving towards the positive zdirection. Density is colourcoded and black lines denote isodensity contours. Velocities are given as arrows, scaled by the speed. 

Open with DEXTER 
Now we turn to astrophysical applications and consider the motion of a gravitating body through a homogeneous gaseous medium. For supersonic speeds a shock front ahead of the moving body is produced that is very similar to the morphology found in the laboratory experiments for the nongravitating bodies. Figure 4 shows the steady state solution for the gas mass density and the velocity field for one of our simulations. In front of the object a bow shock forms where the material is decelerated from supersonic to subsonic speeds. After passing the shock front, matter close to the object settles into a hydrostatic envelope. The major difference to the nongravitating case is that behind the moving object a wake of higherthanaverage density instead of lower density is formed, which is a direct consequence of the gravitational attraction of the body. In turn, this wake of higher density yields a gravitational pull onto the object, which slows it down. This is the phenomenon of dynamical friction.
Overview of the series of simulations performed.
As shown above in Eq. (17) the total drag acting on a gravitating rigid body is the sum of the hydrodynamic drag (pressure force on the body, ) and the dynamical friction, F_{DF}. In our case, behind the shock front a spherical hydrostatic shell forms around the object such that the total pressure force on the object is negligible. Hence, the total drag on the object is given solely by dynamical friction, and in the following we concentrate on this part only (see Fig. 6 below). For very high Mach numbers the separation of the shock from the object’s surface becomes very small such that no hydrostatic envelope can form, and pressure effects will become important again.
The dynamical friction of a body moving with supersonic speed through a gaseous homogeneous medium denotes a welldefined problem, which involves a manageable amount of dependencies. In the course of this section, we first give an overview of the relevant problem parameters and how they are expected to affect the amount of dynamical friction (Sect. 5.1). Afterwards, we present our realized simulation series, each dedicated to investigating the impact of a specific parameter, and we discuss the simulation results in terms of scaling laws (Sects. 5.5 to 5.7). Most importantly, the simulations’ outcome reveal the standoff distance as the minimum spatial scale of the forming anisotropic density structure around the moving body. We analyse and discuss this new aspect indepth and derive a semianalytical relation between the standoff distance and the accretion radius from fits to the numerical data in Sect. 5.3. Finally, we combine these findings to derive a convenient semianalytical expression for the dynamical friction (Sect. 6).
We would like to point out again that we do not assume an a priori validity of the classic drag formula as stated in Eq. (1). Our procedure is to check the scaling of each parameter in Eq. (1) individually and confirm the existence of the logarithmic term. After this we demonstrate that the minimum distance s_{min} in Eq. (2) is closely related to the standoff distance of the shock from the object. We then present a new formula for calculating s_{min} from the basic physical parameter of the problem.
5.1. Relevant parameters of the problem
First of all, the moving object is characterized by its mass M and velocity V_{∞}. The impact of these parameters on the dynamical friction is analysed and discussed in Sects. 5.5 and 5.6. For high enough masses or small enough velocities, the smallest scale of interaction with the gas is set by the gravity of the moving object rather than its geometrical radius. In these cases, the moving body can also be treated as a point mass, as is usually done in semianalytical approaches, e.g. in the derivation by Ostriker (1999). Its geometrical radius R becomes important for the interaction with the surrounding gas in cases of either low mass or high velocities. We investigate the impact of the geometrical radius in Sect. 5.3.
The gaseous homogeneous medium is characterized by two of the four quantities: mass density ρ_{∞}, pressure p_{∞}, temperature T_{∞}, or soundspeed c_{∞}. The remaining two are then determined according to Eqs. (12)and (13). In the following, we will choose the gas mass density ρ_{∞} and soundspeed c_{∞} as the independent variables; the mass density enters directly the force term of the gravitational pull of the wake onto the moving body and the soundspeed of the medium sets the Mach number ℳ = V_{∞}/c_{∞} of the shock.
The thermodynamics of the gas is controlled by its caloric and thermal equations of state Eqs. (11)and (12)with the adiabatic index γ as the only free parameter. Most importantly for the shock physics, the adiabatic index controls how the gaseous medium reacts in case of compression and expansion.
While the object is moving through the gaseous medium, a wake of higher density will form behind the object, which grows in time. For a known interaction time t, the extent of the wake is constrained by s_{max} = V_{∞}t. We discuss the dependence of the dynamical friction on the maximum extent of the wake in more detail in Sect. 5.4.
In summary, the dynamical friction of a massive body moving supersonically through a gaseous homogeneous medium depends only on the mass of the body M, the velocity of the body V_{∞}, its geometrical radius R, the mass density ρ_{∞} of the gas, its soundspeed c_{∞}, the adiabatic index γ, and the maximum extent of the wake s_{max}.
In Ruderman & Spiegel (1971), the authors give an expression for the dynamical friction as given in Eq. (1). In this formula, the dynamical friction also depends on the minimum length scale s_{min} of the interaction, a result of the spatial integration limits of the analytical derivation. From the analysis of the relevant parameters above, it follows that the parameter s_{min} is actually not a free parameter, but has to be a function of the parameters given above. In Sect. 5.3, we compute the minimum length scale of the interaction from our numerical solutions, associate this length scale with the shock’s standoff distance, and derive its dependence on the relevant problem parameters.
In the following, we present the numerical experiments performed. The dependence of the dynamical friction on each of the relevant parameters is determined in a single or multiple dedicated simulation series. Table 1 gives an overview of the simulation series and their physical and numerical parameters. Owing to the scalefreedom in the equations of the problem, the gas mass density can be chosen arbitrarily; we use a value of ρ_{∞} = 1.5 × 10^{4} gcm^{3} in all the simulations performed. In the table, the velocity of the body is given in units of Mach. Simulations with varying soundspeed of the medium (series “S”) use the same velocity of the perturber, hence yield different Mach numbers as well. Simulations with varying adiabatic index of the medium (series “G”) use varying initial values of the gas pressure to keep the soundspeed the same in all simulations of the series.
5.2. The gas mass density
The dynamical friction should scale linearly with the density of the environment: (18)This scaling behaviour is a direct consequence of the fact that the total dynamical friction is given by the sum of all the gravitational pulls from the environment onto the body. Each gravitational pull scales linearly with the density of the environment. Additionally, the hydrodynamical equations (see Sect. 3.2) are scalefree in density, and so the flow morphology is independent of the initial density of the medium.
5.3. Minimum spatial scale of interaction
The total dynamical friction acting on the moving body is given by the spatial integral over the gravitational pull of the surrounding gaseous medium. According to Eqs. (1)and (2), it is expected to scale with the value of a minimum spatial interaction scale s_{min} according to the Coulomb logarithm: (19)Below this minimum interaction scale, the gaseous medium is assumed to exert no net gravitational pull on the body. What actually sets the minimum interaction scale is an open discussion in the literature.
5.3.1. Spatial force analysis
In our study, the approach of direct numerical experiments allows us to properly determine the minimum spatial interaction scale in the case of gaseous media quantitatively. In order to derive the impact of each spatial scale individually, we first calculate the gravitational drag force acting on the object from the gas inside a specific shell at radius r_{i} with a radial shell thickness of r_{i + 1/2}−r_{i−1/2}, (20)where the volume of shell S(r_{i}) is given by (21)The total gravitational force acting on the object is given by the sum over all shells: (22)Even though the above Eq. (20)is written in vector form, only its zcomponent is nonzero due to the axial symmetry of the problem. The net gravitational drag will lead to a slowdown of the body, which is moving along the positive zdirection.
In Fig. 5, we show the fractions of the gravitational force f(r_{i}) as a function of radius or distance to the moving body, respectively.
Fig. 5 Fractions of the gravitational drag force f(r_{i}) to the total gravitational drag force acting on the moving object as a function of distance to the object. The vertical solid lines mark the standoff distance R_{SO} of the shock ahead of the object, and the dashed vertical lines denote the accretion radius R_{A} according to Eq. (3). Results are shown for simulations with three different Mach numbers ℳ as labelled in the bottom left corners. 

Open with DEXTER 
The spatial analysis of the force is shown for three simulations with Mach numbers ℳ = 2,4,and8. Simulation parameters are given in Table 1, series “V4”.
The spatial analysis clearly reveals the socalled standoff distance R_{SO} of the shock front as the searched for minimum spatial interaction scale s_{min}: (23)The associated bow shock forming in front of the body is visible in the density morphology depicted in Fig. 4. The shock front R_{SO} denotes the radius, where the flow changes from supersonic to subsonic velocities. Moreover, the density morphology around the body is characterized by the formation of a hydrostatic envelope extending up to R_{SO}. This hydrostatic envelope is very close to spherical symmetry, as also depicted in Fig. A.1. Hence, although this region marks the highest gas mass density, and in principle might have the strongest gravitational impact owing to its closeness, the net gravitational pull within this hydrostatic envelope (r<R_{SO}) turns out to be negligible. The stronger compression in the forward direction of the trajectory of the moving body yields a slight deviation from spherical symmetry, which actually results in a small gravitational acceleration instead of a drag, but this effect is negligible compared to the total dynamical friction and is only marginally visible in Fig. 5 by eye for the case of highest Mach number; here the fraction of the dynamical friction force close to but still below the standoff distance has a positive sign, indicating positive acceleration of the body.
Fig. 6 Force contributions for the equilibrium state according to Eq. (17) for the gravitating, rigid sphere, given as a function of the Mach number. denotes the pressure force measured at the surface of the sphere, denotes the total dynamical friction calculated according to Eq. (22), the pressure force measured at the outer boundary of the computational domain, and the momentum transport through the outer boundary. The sum of these four forces is denoted by F_{tot}. 

Open with DEXTER 
Clearly, the shells inside the sphere around the object with radius R_{SO} do not contribute to the total dynamical friction force, because around the object a spherically symmetrical hydrostatic envelope forms. This is confirmed in Fig. 6 where we plot the individual contributions to the total drag force on the spherical body as a function of Mach number. As shown, the contribution of the pressure force, , acting on the object is negligible. From this we can conclude that i) to calculate the drag acting on the object it is sufficient to consider the dynamical friction alone; and that ii) the relevant quantity for the minimum interaction scale s_{min} is given by R_{SO}. Additional simulations with fixed Mach number ℳ but different object masses M give the same result for the spatial analysis of the force.
5.3.2. Convenient expression for the standoff distance
It is the aim of this study to derive a general expression for the dynamical friction that allows the acting gravitational drag to be computed from the relevant problem parameters without the need of direct numerical simulations. As revealed in the previous section, this expression for the dynamical friction will include a scaling with (24)In the next step of our investigation, we derive an expression for the standoff distance, which is based on the basic problem parameters only. Therefore, we make use of the apparent relationship of the standoff distance R_{SO} to the accretion radius R_{A} (see Eq. (3)for the definition of the accretion radius).
Figure 5 seems to indicate that R_{SO} is directly proportional to the accretion radius R_{A}. Since R_{A} scales as we check these scaling laws for the standoff distance as well. We use a series of simulations where we systematically vary the mass M of the object and its velocity V_{∞}. Simulation parameters are given in Table 1, series “M2” and “V4”. The results are shown in Fig. 7.
Fig. 7 Correlation of the shock’s standoff distance R_{SO} with object velocity or Mach number ℳ (top panel) and with object mass M (bottom panel). Solid lines denote linear fits in the loglog plane. 

Open with DEXTER 
Both panels show leastsquare fit lines that support the fact that R_{SO} ∝ M/ ℳ^{2}. For larger ℳ the standoff distance becomes smaller according to R_{SO} ∝ ℳ^{2} as shown in the top panel of Fig. 7. For fixed ℳ, lowering the object mass yields a smaller standoff distance in agreement with R_{SO} ∝ M as shown in the bottom panel of Fig. 7. This outcome confirms the scaling laws as long as the geometrical extent of the body R is small compared to R_{SO}. For high Mach shocks and for low body masses, respectively, the standoff distance cannot decrease to smaller radii and approaches the size of the body R instead. In all other cases, R_{SO} is directly proportional to R_{A}: (25)An important difference between the two radii is given by the fact that the standoff distance depends on the thermodynamics of the system. As a consequence, we study the impact of the adiabatic index γ on the standoff distance and the final dynamical friction term in Sect. 5.7 below.
Fig. 8 Correlation of the shock’s standoff distance R_{SO} with the nonlinearity parameter η, as given in Eq. (26), for various model sequences. In the top panel, R_{SO} is normalized by the object radius R. For η> 1, a dashed line representing the relation R_{SO}/R = η is superimposed. In the bottom panel, R_{SO} is normalized by the corresponding accretion radius. 

Open with DEXTER 
To investigate further the proportionality between the two radii, we perform additional simulations varying the mass, radius, and velocity of the object as well as the soundspeed of the medium. Detailed simulation parameters are given in Table 1. To find a useful general expression for the standoff distance R_{SO} we utilize the nonlinearity parameter (26)as introduced by Kim & Kim (2009). In Fig. 8 we show the standoff distance R_{SO} and the ratio of R_{SO}/R_{A} as a function of the nonlinearity parameter η. As seen in the top panel, for η> 1 (larger M, lower ℳ) we find the scaling R_{SO}/R ∝ η with an unknown proportional constant close to unity as indicated by the added dashed line. This is in good agreement with the results of Kim & Kim (2009). The largest deviation from this trend is seen in the rightmost green triangle, which belongs to a simulation in the weakly supersonic regime ℳ = 1.25. A spatial analysis of the final configuration of this model shows that behind the shock front and around the object the gaseous medium is not in hydrostatic equilibrium – as for the supersonic cases – but displays small residual motions that did not disappear even in longterm runs. So, we expect deviations from the following relation (27)for nonsupersonic flows, a regime that is not part of the present study and that will have to be explored in future studies. For comparison, we denote that the rightmost red square belongs to a simulation with ℳ = 2, but even larger η than the green triangle run, and matches the found relationship much more closely.
For η< 1 (larger ℳ, lower M) we enter the regime where the standoff distance is limited by the size of the object R, and therefore the ratio R_{SO}/R remains constant. For this regime of low η, Kim & Kim (2009) find a different scaling because they used a smoothed gravitational potential to model their perturber instead of the solid surface used in our study. Using the definition of η (Eq. (26)), the proportionality R_{SO}/R ∝ η implies (27)Relation (27)represents a suitable expression for the standoff distance. As pointed out above, the proportionality constant for the simulations performed so far seems to be very close to unity, but the relation still lacks the effect of different equations of state, i.e. any dependence on the adiabatic index γ. We investigate this dependence in detail in the following section.
The bottom panel of Fig. 8 confirms again the proportionality between the standoff distance and the accretion radius for different Mach numbers, object masses, and soundspeeds in the regime R_{SO} ≫ R. We also checked that a change in the size R of the moving body does not influence the standoff distance as long as R remains substantially smaller than R_{SO}. Hence, the curves for R = 0.1 and R = 0.01 are identical.
5.4. Maximum extent of the wake
The maximum extent of the wake past the body determines the size of the perturbed region, which will contribute to the total dynamical friction via its gravitational pull. In Eq. (1), the maximum extent of the wake enters as the outer radius of the integral s_{max} in the Coulomb logarithm C_{A}: (28)For an infinite medium, the size of the wake s_{max} is given by the duration of the movement of the perturbing object. For a body starting its journey at t = 0 with a constant velocity V_{∞}, the timedependent size of the wake s_{max} at a later time t is given by the distance travelled s_{max}(t) = V_{∞}t.
Relation (28)implies that the dynamical friction increases as long as the body is moving through the medium, leading to a larger and larger wake. In principal, the dynamical friction seems to approach an infinitely high force for an infinitely long travelling object. In practice, the dynamical friction cannot grow infinitely, not even for an infinitely large medium, because the dynamical friction acting on the body leads to a slowdown of the object on a dynamical timescale of the gasbody interaction (SánchezSalcedo & Brandenburg 2001). The dynamical timescale of dynamical friction is given as the ratio of the initial momentum of the moving object and the acting force t_{DF} = MV_{∞}/F_{DF}. Hence, the velocity of the perturber monotonically decreases in time until the body is at rest. As a consequence, the dynamical friction first increases (owing to a growing extent of the wake) to a maximum value and then decreases (owing to the slowdown of the body).
Here, we model the evolution of the system on timescales much smaller than the timescale for the slowdown of the body to extract the acting dynamical friction as a function of the relevant initial system parameters. As a next step, we check the scaling of the dynamical friction force with the extent of the wake. We compute the resulting dynamical friction by numerically integrating the gravitational pull of the forming wake onto the body up to different radii within the computational domain. As a result, we obtain the individual contributions of the wake at each radius. For further analysis, we extended the size of the computational domain to R_{domain} = 25 R_{⊙} (model C); other simulation parameters correspond to the fiducial case (model F) given in Table 1.
As shown in Fig. 9, the scaling of the numerically determined drag force with the extent of the wake follows the expected logarithmic relation.
Fig. 9 Scaling of dynamical friction with the size of the forming wake for model C in Table 1. 

Open with DEXTER 
In the following, simulations make use of an arbitrary size of the computational domain, which should be much larger than the standoff distance of the shock. Finally, the acting dynamical friction can then be determined using the scaling relation Eq. (28). We use an outer radius of the computational domain of either 100 R or 1000 R, depending on the specified radius R of the moving body.
5.5. Mass of the perturber
The dimensional part of the dynamical friction formula – as written in Eqs. (1)and (5)– is expected to scale with the mass of the object squared: (29)The physical explanation for this scaling relation is that the mass of the object first causes the accumulation of high density in the wake (a gravitational backreaction, which scales linearly with the mass of the object) and secondly, the gravitational pull of the wake onto the object denotes a gravitational interaction, which again scales linearly with the mass of the object. This argument is based on a linearization of the equations to relate the forming wake to a linear scaling with M.
Within our numerical framework, we can compute the total dependence of the dynamical friction on the mass, and can further check, if the dependence on M can be properly split into a scaling with M^{2} in the dimensional part and its further dependence within the dimensionless Coulomb logarithm C_{A}. We have performed numerical simulations for various masses M of the perturber. Simulation parameters are given in Table 1, series “M2”. The resulting scaling of the ratio of dynamical friction and the Coulomb logarithm is shown in Fig. 10.
Fig. 10 Correlation of the ratio of dynamical friction and Coulomb logarithm with the mass of the perturber. 

Open with DEXTER 
The numerical experiments confirm the expected scaling very accurately.
5.6. Velocity of the perturber
The dimensional part of the dynamical friction formula – as written in Eqs. (1)and (5)– is expected to scale with the inverse of the velocity of the object squared: (30)This scaling can only be true for supersonic motion and actually denotes a peculiarity of the dynamical friction because usually hydrodynamical drag forces (see Eq. (5)) scale with the square of the velocity, . In the case of dynamical friction, however, a faster body is already further away from the highdensity wake once it has formed. The distance of the body to the forming wake scales linearly with its velocity and the gravitational pull scales with the inverse of the distance squared.
As shown above, this inverse scaling with of the dynamical friction force can be understood in terms of an effective cross section that is given by the accretion radius, R_{acc}, that decreases with higher velocities of the body. This gives rise to the scaling of the dynamical friction with the inverse of the velocity of the object squared.
We have performed numerical simulations for various velocities V_{∞} of the perturber. Simulation parameters are given in Table 1, series “V4”. The resulting scaling of the ratio of dynamical friction and drag coefficient is shown in Fig. 11.
Fig. 11 Correlation of the ratio of dynamical friction and Coulomb logarithm with the velocity of the perturber. 

Open with DEXTER 
Again, the numerical experiments confirm the expected scaling. Deviations from the scaling law occur for very high Mach numbers because the standoff distance cannot shrink to arbitrarily small values and approaches the geometrical radius of the body instead (compare discussion in Sect. 5.3 for details, especially Fig. 7).
5.7. Adiabatic index
All the simulations presented in the previous sections use an adiabatic index of γ = 5/3. This was also the one and only value for the adiabatic index investigated in Kim & Kim (2009). However, the compression in the shocked region around the object, as well as the compression in the wake behind the moving body, depends on the equation of state in use and so on the value of the adiabatic index. In an extreme case (e.g. for very large values of the adiabatic index), the gaseous medium becomes incompressible because any compression (i.e. increase in density) causes the pressure to increase infinitely, which causes an infinitely strong pressure force, and the medium directly relaxes towards an isodensity morphology again.
To further investigate the generality of our results so far we perform a variety of additional model sequences using different adiabatic indices. Detailed setup parameters of these simulation series are given in Table 1.
The correlation of the standoff distance with the velocity of the perturber is shown in Fig. 12, but now for a variety of different values of the adiabatic index γ.
Fig. 12 Correlation of the shock’s standoff distance with the Mach number for different values of the adiabatic index γ. Dots denote results from the numerical simulations. The coloured lines represent the scaling of R_{SO}/R ∝ ℳ^{2}. 

Open with DEXTER 
For large Mach numbers ℳ, the standoff distance approaches the geometrical radius of the object. For lower ℳ, as indicated by the coloured lines, all simulation series confirm the scaling law of R_{SO}/R ∝ ℳ^{2} as previously determined in Fig. 7.
In Fig. 13 top panel, the correlation of the standoff distance with the nonlinearity parameter η is shown for different values of γ.
Fig. 13 Correlation of the shock’s standoff distance with the nonlinearity parameter η given in Eq. (26)(top panel) and the Mach number (bottom panel). The top panel is analogous to the top panel of Fig. 8, but now with results for different values of the adiabatic index γ. Dots denote results from the numerical simulations. Coloured lines represent the relations discussed in the main text. 

Open with DEXTER 
Coloured lines indicate linear fits in the loglog plane of the plot. Since we found R_{SO}/R ∝ η for γ = 5/3 in agreement with Kim & Kim (2009), we assume the following more general relation, (31)where the factor g(γ) describes the dependence on the adiabatic index γ. Specific values of this function are obtained according to the linear fits to the numerical experiments, namely
For large values of the adiabatic index, namely γ = 3 and 4, the accuracy of the fit values has to be taken with care, because the fits only rely on three or four data points.
In the bottom panel of Fig. 13, the standoff distance is shown in units of the accretion radius as a function of the Mach number. Here too the previously determined correlation between the standoff distance and the accretion radius (Eq. (27)) can be generalized to include the dependence on the adiabatic index γ via (32)This function is presented as coloured lines in the bottom panel of Fig. 13. The values for g(γ) are as given above.
In Appendix B, we derive an analytical estimate of the γdependence of the standoff distance: (33)This relation is shown in comparison to a simulation series of varying adiabatic index in Fig. 14.
Fig. 14 Standoff distance in units of accretion radius R_{SO}/R_{A} as a function of the adiabatic index γ. The black line denotes the analytically estimated upper limit given in Eq. (33). 

Open with DEXTER 
The estimate does not allow an exact determination, but it denotes an upper limit of the standoff distance. More importantly, the estimate follows the trend of the γdependence when compared to the numerical experiments.
Owing to the similarity to the γdependence, we choose the analytical estimate as the basis function of a fit to the numerical data, but allow now for a shift in γ → γ + a and in R_{SO}/R_{A} → R_{SO}/R_{A} + b(34)with a and b as free fitting parameters. A leastsquare fit of this function to the numerical data yields a = 0.1 and b = 0.18.
Using Eq. (32)for the ℳ = 4.0 numerical data in combination with the fitted function Eq. (34)gives an approximate function of g(γ)(35)with a = 0.1 and b = 0.18.
As a further check of the approximate function, we can compare this relation with the tabular values above derived via Eq. (31), see Fig. 15.
Fig. 15 Comparison of the semianalytical approximation of the function g(γ), see Eq. (35), shown as a black line with numerically obtained values. Small red dots denote the simulation series “G” with varying adiabatic index; each red dot represents a single simulation. Larger blue dots denote values derived by fitting the scaling of the standoff distance with the nonlinearity parameter for a variety of different simulation series shown in the top panel of Fig. 13; each blue dot represents such a fit value to a full simulation series. See main text for details of the derivation. 

Open with DEXTER 
The approximate function g(γ) for the overall γdependence gives reasonable results in comparison to the numerical outcome. As stated previously, for large values of the adiabatic index, namely γ = 3 and 4, the accuracy of the fit values (big blue dots) has to be taken with care, because the fits rely only on three or four data points.
6. New formula of dynamical friction
As a final step, we combine our obtained scaling law for the standoff distance, R_{SO}, and the dependence on the adiabatic index, γ, to formulate a new general expression for the dynamical friction of a gravitating body of mass M and radius R moving supersonically with velocity V_{∞} through a homogeneous gaseous medium of mass density ρ_{∞} and soundspeed c_{∞}: (36)The determination of the drag force from the basic problem parameters only requires an a priori knowledge of the standoff distance. As shown above, the shock’s standoff distance R_{SO} can be approximated by (37)with g(γ) given by Eq. (35).
Strictly speaking, the formula above for the standoff distance is only valid for a standoff distance larger than the geometrical radius of the object. In the regime of large, lowmass bodies moving with high velocity, the standoff distance approaches the value of the geometrical radius instead, and s_{min} can be replaced by the object size R.
We check the overall accuracy of this approximate determination of dynamical friction in direct comparison to the various numerical experiments performed here, covering a broad parameter space in terms of the dimensionless nonlinearity parameter η = [2 (ℳ^{2}−1)/ℳ^{2}R/R_{A}] ^{1}. The result of this comparison is shown as the ratio of the force in the experiments and the approximate formula in Fig. 16.
Fig. 16 Comparison of the newly derived approximate formula for dynamical friction F_{DF} with data from various numerical experiments . The corresponding ratio is given as function of the nonlinearity parameter η. The top panel shows results from a variety of simulations with adiabatic index γ = 5/3, varying the Mach number, soundspeed, and object mass. The bottom panel shows results from a variety of simulations with varying Mach number for different adiabatic indices as labelled. 

Open with DEXTER 
The approximate formula can be used as a convenient tool to estimate the drag force from the basic parameters in reasonable accuracy. Especially in the regime of 1 ≤ η ≤ 50 the estimated values match the numerical experiments. The largest deviations are observed in the marginally supersonic regime (ℳ close to unity, denoted by the rightmost greenish triangle) and in cases where the shock’s standoff distance approaches the size of the object. This effect is visible when comparing the numerical outcome in the regime of large η for simulations with larger (green diamonds) and smaller object radii (red squares) in the top panel of Fig. 16.
Fig. 17 Normalized drag force as a function of the Mach number. The solid line resembles our newly derived formula as given in Eq. (36). The numerical data are given by the black dots, for which we we used simulation series V4. 

Open with DEXTER 
To compare how well our newly derived formula approximates the drag force, F_{DF}, we plot in Fig. 17 the numerically calculated drag force divided by the factor 4πρ_{∞}^{(}GM/c_{∞}^{)}^{2} for varying Mach numbers of the perturber as black dots. The relation (38)is shown as a solid black line, where we use C_{A} = ln(s_{max}/R_{SO}) with R_{SO} according to Eq. (32). For the maximum radius in C_{A} we used r_{max} = 1000R. As can be seen, our approximative function for C_{A} reflects the numerical data very well.
7. Summary
We derived a new formula for dynamical friction of a body moving with supersonic speed through a homogeneous gaseous medium. This formula was obtained by following an ansatz of direct numerical modelling of the dynamical problem.
In 11 simulation series that include a total of slightly more than 100 individual simulations, we scanned the parameter space of the problem and derived the scaling relations of the dynamical friction with the mass, velocity, and radius of the moving body; the gas mass density, its soundspeed, and the value of the adiabatic index; and the maximum extent of the forming wake.
7.1. Scaling relations
As expected, the dimensional part of the dynamical friction formula as given in Eq. (36)scales proportionally to the mass of the moving body squared, and to the inverse of its velocity squared. Furthermore, the dynamical friction F_{DF} scales linearly with the gas mass density and is proportional to the logarithm of the ratio of the spatial extents of the forming wake and the shock’s standoff distance.
7.2. Minimum spatial interaction scale
The numerical experiments allow us to compute the total dynamical friction acting on the body and also to investigate the individual contributions to the dynamical friction induced as a function of the distance to the perturber. The outcome of this analysis reveals that the region within the standoff distance R_{SO} of the shock does not contribute to the dynamical friction of a gaseous medium (see Fig. 5 and associated main text for details). At radii smaller than the standoff distance, a spherically symmetric stratified atmosphere forms around the moving body, which (because of its symmetry) induces no net gravitational pull or pressure force. Only at radii above the standoff distance does a wake of higher density form behind the moving body and induce a gravitational drag force.
Hence, in the wellknown formulae by Ruderman & Spiegel (1971) or Ostriker (1999), it is the standoff distance R_{SO} of the shock that determines the minimum length scale s_{min} required in the socalled Coulomb logarithm C_{A} = ln(s_{max}/s_{min}): (39)The formation of a spherically symmetric atmosphere (with radius R_{SO}) around a gravitating object implies that the hydrodynamical drag force vanishes in contrast to the nongravitating case, but for very high Mach numbers the shock moves very close to the object so that R_{SO} ≈ R. In this case, no spherically symmetric envelope around the object forms and the hydrodynamical drag becomes important again.
7.3. Relating the standoff distance to the accretion radius
In a second step, to allow an easy computation of the standoff distance for general dynamical friction problems, we derived the relation of the shock’s standoff distance to the common definition of the accretion radius, (40)depending only on the mass M and velocity V_{∞} of the perturber. The standoff distance is directly proportional to the accretion radius and the proportionality constant only depends on the Mach number and the adiabatic index of the gas: (41)Additionally, the geometrical radius R of the moving body yields a lower limit for the standoff distance (R_{SO} ≥ R).
Via an analytical estimate of the standoff distance and further fitting of the remaining free parameters, we gave an approximate relation of the dependence of the shock’s standoff distance on the adiabatic index γ of the gaseous medium (42)with γ′ = γ + a and values of a = 0.1 and b = 0.18.
7.4. Updated formula of dynamical friction
Finally, we combined the derived scaling relations and the finding on the standoff distance to give an update of the wellknown formula of dynamical friction acting on a body of mass M moving at supersonic speed V_{∞} through a gaseous homogeneous medium of mass density ρ_{∞}, (43)where s_{max} denotes the extent of the wake and R_{SO} can be obtained using the derived formulae above.
Acknowledgments
We thank Chris Ormel for very helpful discussions. R.K. acknowledges financial support from the EmmyNoetherProgram of the German Research Foundation (DFG) under grant No. KU 2849/31. Part of the numerical simulations were performed on the bwGRiD cluster in Tübingen, which is funded by the Ministry for Education and Research of Germany and the Ministry for Science, Research and Arts of the state BadenWürttemberg, and the cluster of the Forschergruppe FOR 759 funded by the DFG. All plots in this paper have been made with the Python library matplotlib (Hunter 2007).
References
 Armitage, P. J., & Natarajan, P. 2005, ApJ, 634, 921 [NASA ADS] [CrossRef] [Google Scholar]
 Billig, F. S. 1967, J. Spacecraft Rockets, 4, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Cantó, J., Raga, A. C., Esquivel, A., & SánchezSalcedo, F. J. 2011, MNRAS, 418, 1238 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, ApJ, 97, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Chavarría, L., Mardones, D., Garay, G., et al. 2010, ApJ, 710, 583 [NASA ADS] [CrossRef] [Google Scholar]
 Edgar, R. 2004, New Astron. Rev., 48, 843 [NASA ADS] [CrossRef] [Google Scholar]
 Foglizzo, T., Galletti, P., & Ruffert, M. 2005, A&A, 435, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grishin, E., & Perets, H. B. 2015, ApJ, 811, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Hunt, R. 1971, MNRAS, 154, 141 [NASA ADS] [Google Scholar]
 Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, H., & Kim, W.T. 2009, ApJ, 703, 1278 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1966, Hydrodynamik (Berlin: AkademieVerlag) [Google Scholar]
 Lee, A. T., & Stahler, S. W. 2014, A&A, 561, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meyer, D. M.A., Gvaramadze, V. V., Langer, N., et al. 2014a, MNRAS, 439, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, D. M.A., Mackey, J., Langer, N., et al. 2014b, MNRAS, 444, 2754 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, D. M.A., Langer, N., Mackey, J., Velázquez, P. F., & Gusdorf, A. 2015, MNRAS, 450, 3080 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R. 2000, ApJ, 536, 663 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, E. C. 1999, ApJ, 513, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Rein, H. 2012, MNRAS, 422, 3611 [NASA ADS] [CrossRef] [Google Scholar]
 Rephaeli, Y., & Salpeter, E. E. 1980, ApJ, 240, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, P. M., & Taam, R. E. 2008, ApJ, 672, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, P. M., & Taam, R. E. 2012, ApJ, 746, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. A., & Spiegel, E. A. 1971, ApJ, 165, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Ruffert, M. 1994, ApJ, 427, 342 [NASA ADS] [CrossRef] [Google Scholar]
 Ruffert, M., & Arnett, D. 1994, ApJ, 427, 351 [NASA ADS] [CrossRef] [Google Scholar]
 SánchezSalcedo, F. J., & Brandenburg, A. 2001, MNRAS, 322, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Shima, E., Matsuda, T., Takeda, H., & Sawada, K. 1985, MNRAS, 217, 367 [NASA ADS] [Google Scholar]
 Slavin, J. D., Jones, A. P., & Tielens, A. G. G. M. 2004, ApJ, 614, 796 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssandier, J., Terquem, C., & Papaloizou, J. C. B. 2013, MNRAS, 428, 658 [NASA ADS] [CrossRef] [Google Scholar]
 van Dyke, M. 1982, An album of fluid motion (Stanford, CA: Parabolic Press) [Google Scholar]
 van Dyke, M. D., & Gordon, H. 1959, NASA, Technical Report R1 [Google Scholar]
 van Marle, A. J., Meliani, Z., Keppens, R., & Decin, L. 2011, ApJ, 734, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Villaver, E., & Livio, M. 2009, ApJ, 705, L81 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Hydrostatic equilibrium stratification
As mentioned in the text, behind the shock front and around the object the material is highly subsonic and can be approximated by assuming a hydrostatic equilibrium. The morphology is illustrated in more detail in Fig. A.1 where we show the stratification of different physical quantities in different directions. The data shown is extracted from a simulation with the fiducial parameters specified in Table 1.
Fig. A.1 Stratification around the moving object. Values in front of the moving object (θ = 0), perpendicular to its velocity (θ = 90), and behind the object (θ = 180) are displayed. From top to bottom the panels show gas density, pressure, Mach number, and entropy of the flow. 

Open with DEXTER 
Fig. A.2 Stratification of gas density (upper panel) and pressure (bottom panel) around the moving object. As in the previous figure, values in front of the moving object (θ = 0), perpendicular to its velocity (θ = 90), and behind the object (θ = 180) are displayed. The additional red lines denote the analytical estimates obtained from Eqs. (A.7)and (A.2), respectively. In the formula of the estimate, a value for the standoff distance of x_{SO} = 0.533 was used, in accordance with the final quasistationary state of the numerical experiment. 

Open with DEXTER 
As is clearly visible, within the standoff radius, as given by the discontinuity of the solid curve, all three cuts agree with each other and so the stratification is spherically symmetric. Because the Mach number is very close to zero, it is hydrostatic as well. Here, we derive analytical relations for this stratification that are used in the following section of the Appendix to obtain an analytical estimate of the shock’s standoff distance.
For spherical symmetry the equation of hydrostatic equilibrium reads (A.1)where φ is the gravitational potential. Now we use an adiabatic approximation for the pressure, (A.2)where K is a constant, and integrate from a reference radius r_{0} to r, yielding (A.3)Here the subscript 0 refers to the reference radius r_{0} where the values are known.
As a next step, we want to analyse the bow shock in front of the object. We can find the physical properties behind the shock from the standard jump conditions Here, the index “1” refers to the preshock values (supersonic regime) and the index “2” to the postshock (subsonic) regime. These are valid in a reference frame moving with the shock. In our case the shock is stationary in the comoving frame of the body, and we can directly apply the jump conditions if we know the preconditions. The simulations show that the preshock values of density, pressure, and temperature at the standoff radius R_{SO} are just the prescribed inflow values (at ∞; cf. Fig. A.1).
To obtain the preshock Machnumber ℳ_{1} at this radius we use the freefall condition. From the conservation of energy (kinetic and gravitational), we obtain (A.6)where the coordinate x denotes the radial distance to the shock front in units of the accretion radius, specifically x_{SO} = R_{SO}/R_{A}. Equation (A.6)matches the simulation outcome quite well (cf. Fig. A.1).
To obtain the (hydrostatic) postshock density stratification around the object we first use these normalizations to obtain the general relation (A.7)with an arbitrary point x_{0}, where ρ_{0} and c_{0} are known values at the reference radius. If we now choose x_{0} = x_{SO}, the other quantities correspond to the desired postshock quantities and can be obtained from the jump conditions using ℳ_{1} from Eq. (A.6). The required shockposition x_{SO} enters at this point as a free parameter obtained from the numerical solution.
Appendix B: Estimate of the standoff distance
Because the stratification around the object is nearly hydrostatic, some limits can be obtained on x_{SO}. Here, we use the strong shock conditions, i.e. ℳ_{1} → ∞ and obtain for the jumps in density and pressure and hence for the postshock soundspeed (in units of c_{∞}) (B.3)We now assume that for large distances, x → ∞, the density of the postshock stratification lies below the unperturbed density ρ(x → ∞) ≤ ρ_{∞} as seen in the simulations. This is shown in Fig. A.2 where the analytical postshock stratification (red solid curve) for large radii lies just below the unperturbed value ρ_{∞}. Inserting Eq. (B.3)into Eq. (A.7)yields (B.4)with ℳ_{1} from Eq. (A.6). Solving for x_{SO} finally leads to (B.5)This function is shown in Fig. 14 in comparison to data from numerical simulations. While the relation (B.5)cannot give the exact value for R_{SO} it nevertheless provides a reasonable upper limit for the standoff distance, and follows the trend of the obtained numerical results.
Appendix C: Convergence study
We ran simulations using varying grid sizes from roughly 10^{3} up to 10^{6} grid cells in total (more precisely, we varied the number of grid cells from N_{r} × N_{θ} = 50 × 23 = 1150 up to N_{r} × N_{θ} = 1720 × 790 = 1.3588 × 10^{6} grid cells). The inner radial boundary of these simulations was chosen as R = 0.01 R_{⊙}. Other simulation parameters correspond to the “fiducial setup” given in Table 1. From the final quasistationary state we determine the two most relevant physical quantities of our investigation, namely the shock’s standoff distance R_{SO} and the dynamical friction acting on the moving body.
The results of this convergence study are presented in Fig. C.1.
Fig. C.1 Convergence test of the numerical setup. Values on the horizontal axis denote the total number of grid cells N_{r} × N_{θ}. The vertical axis shows the relative deviation of the resulting standoff distance R_{SO} (solid circles) and the dynamical friction (dashed circles) from its corresponding value of the highest resolution run. 

Open with DEXTER 
Deviations are given as relative differences to the values from the highest resolution simulation using more than 10^{6} grid cells. The numerical results show a clear convergence trend. For numerical grids larger than 10^{5} grid cells, the deviations become negligibly small. To be on the very safe side, we chose N_{r} × N_{θ} = 860 × 400 = 3.44 × 10^{5} grid cells as our default resolution for the simulations performed. Simulations with a larger inner radial boundary of R = 0.1 R_{⊙} require fewer grid cells in the radial direction, especially owing to the logarithmic grid in the radial direction. Nonetheless, we also use here a comparable grid size of N_{r} × N_{θ} = 700 × 480 = 3.36 × 10^{5} grid cells as our default resolution.
All Tables
All Figures
Fig. 1 Visualization of a comparison simulation result of a sphere embedded in air flow (γ = 1.4) at a Mach number ℳ = 1.53. The homogeneous gas arrives from the top and flows in the negative zdirection. Shown is the density distribution for the final equilibrium state in the central part of the computational domain around the sphere. Upper panel: numerical results, black lines denote isodensity contours. Bottom panel: laboratory data from van Dyke (1982, Fig. 266). 

Open with DEXTER  
In the text 
Fig. 2 Quantitative comparison of the shock’s standoff distance as function of Mach number. The solid black line denotes the relationship derived from laboratory experiments as given in the main text. Small red dots represent results from numerical experiments for the same setup. Blue dots denote early numerical results by van Dyke & Gordon (1959). 

Open with DEXTER  
In the text 
Fig. 3 Force contributions in the equilibrium state according to Eq. (17) for the nongravitating, rigid sphere in air, given as a function of the Mach number. denotes the pressure force measured at the surface of the sphere, the pressure force measured at the outer boundary of the computational domain, and the momentum transport through the outer boundary. The sum of these three forces is denoted by F_{tot}. The curve describes a parabola with F ∝ ℳ^{2}. 

Open with DEXTER  
In the text 
Fig. 4 Density and velocity structure of the final quasistationary state around a gravitating object moving with supersonic motion through a gaseous homogeneous medium. The body is moving towards the positive zdirection. Density is colourcoded and black lines denote isodensity contours. Velocities are given as arrows, scaled by the speed. 

Open with DEXTER  
In the text 
Fig. 5 Fractions of the gravitational drag force f(r_{i}) to the total gravitational drag force acting on the moving object as a function of distance to the object. The vertical solid lines mark the standoff distance R_{SO} of the shock ahead of the object, and the dashed vertical lines denote the accretion radius R_{A} according to Eq. (3). Results are shown for simulations with three different Mach numbers ℳ as labelled in the bottom left corners. 

Open with DEXTER  
In the text 
Fig. 6 Force contributions for the equilibrium state according to Eq. (17) for the gravitating, rigid sphere, given as a function of the Mach number. denotes the pressure force measured at the surface of the sphere, denotes the total dynamical friction calculated according to Eq. (22), the pressure force measured at the outer boundary of the computational domain, and the momentum transport through the outer boundary. The sum of these four forces is denoted by F_{tot}. 

Open with DEXTER  
In the text 
Fig. 7 Correlation of the shock’s standoff distance R_{SO} with object velocity or Mach number ℳ (top panel) and with object mass M (bottom panel). Solid lines denote linear fits in the loglog plane. 

Open with DEXTER  
In the text 
Fig. 8 Correlation of the shock’s standoff distance R_{SO} with the nonlinearity parameter η, as given in Eq. (26), for various model sequences. In the top panel, R_{SO} is normalized by the object radius R. For η> 1, a dashed line representing the relation R_{SO}/R = η is superimposed. In the bottom panel, R_{SO} is normalized by the corresponding accretion radius. 

Open with DEXTER  
In the text 
Fig. 9 Scaling of dynamical friction with the size of the forming wake for model C in Table 1. 

Open with DEXTER  
In the text 
Fig. 10 Correlation of the ratio of dynamical friction and Coulomb logarithm with the mass of the perturber. 

Open with DEXTER  
In the text 
Fig. 11 Correlation of the ratio of dynamical friction and Coulomb logarithm with the velocity of the perturber. 

Open with DEXTER  
In the text 
Fig. 12 Correlation of the shock’s standoff distance with the Mach number for different values of the adiabatic index γ. Dots denote results from the numerical simulations. The coloured lines represent the scaling of R_{SO}/R ∝ ℳ^{2}. 

Open with DEXTER  
In the text 
Fig. 13 Correlation of the shock’s standoff distance with the nonlinearity parameter η given in Eq. (26)(top panel) and the Mach number (bottom panel). The top panel is analogous to the top panel of Fig. 8, but now with results for different values of the adiabatic index γ. Dots denote results from the numerical simulations. Coloured lines represent the relations discussed in the main text. 

Open with DEXTER  
In the text 
Fig. 14 Standoff distance in units of accretion radius R_{SO}/R_{A} as a function of the adiabatic index γ. The black line denotes the analytically estimated upper limit given in Eq. (33). 

Open with DEXTER  
In the text 
Fig. 15 Comparison of the semianalytical approximation of the function g(γ), see Eq. (35), shown as a black line with numerically obtained values. Small red dots denote the simulation series “G” with varying adiabatic index; each red dot represents a single simulation. Larger blue dots denote values derived by fitting the scaling of the standoff distance with the nonlinearity parameter for a variety of different simulation series shown in the top panel of Fig. 13; each blue dot represents such a fit value to a full simulation series. See main text for details of the derivation. 

Open with DEXTER  
In the text 
Fig. 16 Comparison of the newly derived approximate formula for dynamical friction F_{DF} with data from various numerical experiments . The corresponding ratio is given as function of the nonlinearity parameter η. The top panel shows results from a variety of simulations with adiabatic index γ = 5/3, varying the Mach number, soundspeed, and object mass. The bottom panel shows results from a variety of simulations with varying Mach number for different adiabatic indices as labelled. 

Open with DEXTER  
In the text 
Fig. 17 Normalized drag force as a function of the Mach number. The solid line resembles our newly derived formula as given in Eq. (36). The numerical data are given by the black dots, for which we we used simulation series V4. 

Open with DEXTER  
In the text 
Fig. A.1 Stratification around the moving object. Values in front of the moving object (θ = 0), perpendicular to its velocity (θ = 90), and behind the object (θ = 180) are displayed. From top to bottom the panels show gas density, pressure, Mach number, and entropy of the flow. 

Open with DEXTER  
In the text 
Fig. A.2 Stratification of gas density (upper panel) and pressure (bottom panel) around the moving object. As in the previous figure, values in front of the moving object (θ = 0), perpendicular to its velocity (θ = 90), and behind the object (θ = 180) are displayed. The additional red lines denote the analytical estimates obtained from Eqs. (A.7)and (A.2), respectively. In the formula of the estimate, a value for the standoff distance of x_{SO} = 0.533 was used, in accordance with the final quasistationary state of the numerical experiment. 

Open with DEXTER  
In the text 
Fig. C.1 Convergence test of the numerical setup. Values on the horizontal axis denote the total number of grid cells N_{r} × N_{θ}. The vertical axis shows the relative deviation of the resulting standoff distance R_{SO} (solid circles) and the dynamical friction (dashed circles) from its corresponding value of the highest resolution run. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.