Free Access
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/0004-6361/201527629
Published online 05 April 2016

© 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 higher-than-average 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 set-up of the simulations. In Sect. 4, we present comparison simulations to laboratory experiments for non-gravitating 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 FDF acting on the mass M: (1)In Eq. (1), ρ denotes the ambient constant density, G is the gravitational constant, and CA stands for the Coulomb logarithm (2)where smin and smax 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 smin and smax need to be known. Typically, for smax the maximum extension of the medium is assumed. For the minimum length, smin, 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, smin ≈ max(R,RA). Here RA stands for the accretion radius (3)Recently, Cantó et al. (2011) derived an approximation for smin 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 time-dependent 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 CA, very similar to that of Rephaeli & Salpeter (1980), which includes a time-dependent maximum radius, smax ~ Vt, 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 time-dependent phenomena. The expression for the supersonic case reads (4)where ℳ = V/c denotes the Mach number of the problem. Obviously, if we set smax = Vt, 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 time-dependent 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 smin 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 Bondi-Hoyle-Littleton 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 pre-requisite – on the determination of the minimum integration radius smin.

The validity of the formula by Ostriker has been demonstrated by Sánchez-Salcedo & Brandenburg (2001) who showed agreement for moderate Mach numbers in non-linear, isothermal simulations. However, they used a Plummer-type potential with a smoothing length much larger than RA. 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 Plummer-type 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 RA. 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 three-dimensional (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 two-dimensional (2D) flows that are non-physical, 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 CD 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 smin of the problem, and finally give a convenient expression for the dynamical friction in general.

3. Physics and numerics

3.1. Problem set-up

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 time-dependent 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 = ekin + eth the total (kinetic and thermal) energy density of the gas. The acceleration due to external forces aext 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, kB the Boltzmann constant, and mH 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 open-source code PLUTO (Mignone et al. 2007), version 4. Part of the simulations were carried out with a new in-house developed CUDA version for usage on GPUs. We use a Runge-Kutta time stepping scheme of second order with a second-order reconstruction of states in space, a van Leer limiter in the reconstruction step, and a Harten-Lax-van Leer solver for the Riemann problem.

Equations (7)to (9)are solved in the co-moving 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 z-direction. 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 non-accreting object. The computational domain extends in the radial direction from the radius R of the object up to Rdomain = 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 Rdomain, 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 zero-gradient boundary condition so that the out-flowing material can leave the computational domain without reflections. As already mentioned, we perform simulations in axial symmetry and therefore at θ = 0 (positive z-axis) and at θ = π (negative z-axis) 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. Self-gravity of the gas is not taken into account. Effects of magnetic fields are not investigated. The long-evolution problem of a time-dependent slow-down of the moving body by the acting dynamical friction is not included in this study (see Sánchez-Salcedo & Brandenburg 2001 for a discussion of this issue). The formula for dynamical friction is derived within the co-moving 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 z-direction and we have to consider the z-component of Eq. (8), which reads in equilibrium (15)where az denotes the z-component 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, FDF, 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, Rdomain. 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 non-gravitating moving bodies. To check the validity of our numerical ansatz and prove the accuracy of the code, including set-up, 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).

thumbnail 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 z-direction. 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 iso-density contours. Bottom panel: laboratory data from van Dyke (1982, Fig. 266).

Open with DEXTER

4.1. Morphology

The black-and-white photograph in van Dyke (1982), Fig. 266, shows the result of a laboratory experiment of a spherical non-gravitating 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 low-density 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 z-direction; actually, simulations are performed in the co-moving frame of the body, hence, the gas flow is initialized into the negative z-direction. 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 low-density 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 well-known analytical Rankine-Hugoniot jump conditions. In contrast to the axisymmetric and inviscid numerical result, the laboratory experiment shows that the low-density 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. Stand-off distance of the shock front

As we point out below, the shock’s stand-off distance, RSO, plays a crucial role in determining the dynamical friction on a gravitating moving body. Here, we measure RSO 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 set-up as in the previous section but with a variety of different inflow velocities V. Figure 2 shows the resulting stand-off distance of the shock, measured from the centre of the sphere, as a function of .

thumbnail Fig. 2

Quantitative comparison of the shock’s stand-off 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 set-up. Blue dots denote early numerical results by van Dyke & Gordon (1959).

Open with DEXTER

The RSO decreases with increasing and approaches asymptotically a value of RSO → 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 RSO/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 stand-off 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

thumbnail Fig. 3

Force contributions in the equilibrium state according to Eq. (17) for the non-gravitating, 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 Ftot. 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, Ftot, 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

thumbnail Fig. 4

Density and velocity structure of the final quasi-stationary state around a gravitating object moving with supersonic motion through a gaseous homogeneous medium. The body is moving towards the positive z-direction. Density is colour-coded and black lines denote iso-density 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 non-gravitating 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 non-gravitating case is that behind the moving object a wake of higher-than-average 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.

Table 1

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, FDF. 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 well-defined 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 stand-off distance as the minimum spatial scale of the forming anisotropic density structure around the moving body. We analyse and discuss this new aspect in-depth and derive a semi-analytical relation between the stand-off distance and the accretion radius from fits to the numerical data in Sect. 5.3. Finally, we combine these findings to derive a convenient semi-analytical 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 smin in Eq. (2) is closely related to the stand-off distance of the shock from the object. We then present a new formula for calculating smin 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 semi-analytical 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 smax = Vt. 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 smax.

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 smin 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 smin 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 stand-off 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 scale-freedom 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 scale-free 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 smin 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 ri with a radial shell thickness of ri + 1/2ri−1/2, (20)where the volume of shell S(ri) 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 z-component is non-zero due to the axial symmetry of the problem. The net gravitational drag will lead to a slow-down of the body, which is moving along the positive z-direction.

In Fig. 5, we show the fractions of the gravitational force f(ri) as a function of radius or distance to the moving body, respectively.

thumbnail Fig. 5

Fractions of the gravitational drag force f(ri) 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 stand-off distance RSO of the shock ahead of the object, and the dashed vertical lines denote the accretion radius RA 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 so-called stand-off distance RSO of the shock front as the searched for minimum spatial interaction scale smin: (23)The associated bow shock forming in front of the body is visible in the density morphology depicted in Fig. 4. The shock front RSO 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 RSO. 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<RSO) 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 stand-off distance has a positive sign, indicating positive acceleration of the body.

thumbnail 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 Ftot.

Open with DEXTER

Clearly, the shells inside the sphere around the object with radius RSO 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 smin is given by RSO. 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 stand-off 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 stand-off distance, which is based on the basic problem parameters only. Therefore, we make use of the apparent relationship of the stand-off distance RSO to the accretion radius RA (see Eq. (3)for the definition of the accretion radius).

Figure 5 seems to indicate that RSO is directly proportional to the accretion radius RA. Since RA scales as we check these scaling laws for the stand-off 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.

thumbnail Fig. 7

Correlation of the shock’s stand-off distance RSO with object velocity or Mach number (top panel) and with object mass M (bottom panel). Solid lines denote linear fits in the log-log plane.

Open with DEXTER

Both panels show least-square fit lines that support the fact that RSOM/ ℳ2. For larger the stand-off distance becomes smaller according to RSO ∝ ℳ-2 as shown in the top panel of Fig. 7. For fixed , lowering the object mass yields a smaller stand-off distance in agreement with RSOM 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 RSO. For high Mach shocks and for low body masses, respectively, the stand-off distance cannot decrease to smaller radii and approaches the size of the body R instead. In all other cases, RSO is directly proportional to RA: (25)An important difference between the two radii is given by the fact that the stand-off distance depends on the thermodynamics of the system. As a consequence, we study the impact of the adiabatic index γ on the stand-off distance and the final dynamical friction term in Sect. 5.7 below.

thumbnail Fig. 8

Correlation of the shock’s stand-off distance RSO with the non-linearity parameter η, as given in Eq. (26), for various model sequences. In the top panel, RSO is normalized by the object radius R. For η> 1, a dashed line representing the relation RSO/R = η is superimposed. In the bottom panel, RSO 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 stand-off distance RSO we utilize the non-linearity parameter (26)as introduced by Kim & Kim (2009). In Fig. 8 we show the stand-off distance RSO and the ratio of RSO/RA as a function of the non-linearity parameter η. As seen in the top panel, for η> 1 (larger M, lower ) we find the scaling RSO/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 long-term runs. So, we expect deviations from the following relation (27)for non-supersonic 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 stand-off distance is limited by the size of the object R, and therefore the ratio RSO/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 RSO/Rη implies (27)Relation (27)represents a suitable expression for the stand-off 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 stand-off distance and the accretion radius for different Mach numbers, object masses, and soundspeeds in the regime RSOR. We also checked that a change in the size R of the moving body does not influence the stand-off distance as long as R remains substantially smaller than RSO. 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 smax in the Coulomb logarithm CA: (28)For an infinite medium, the size of the wake smax 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 time-dependent size of the wake smax at a later time t is given by the distance travelled smax(t) = Vt.

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 slow-down of the object on a dynamical timescale of the gas-body interaction (Sánchez-Salcedo & 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 tDF = MV/FDF. 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 slow-down of the body).

Here, we model the evolution of the system on timescales much smaller than the timescale for the slow-down 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 Rdomain = 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.

thumbnail 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 stand-off 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 back-reaction, 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 M2 in the dimensional part and its further dependence within the dimensionless Coulomb logarithm CA. 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.

thumbnail 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 high-density 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, Racc, 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.

thumbnail 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 stand-off 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 iso-density 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 set-up parameters of these simulation series are given in Table 1.

The correlation of the stand-off distance with the velocity of the perturber is shown in Fig. 12, but now for a variety of different values of the adiabatic index γ.

thumbnail Fig. 12

Correlation of the shock’s stand-off 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 RSO/R ∝ ℳ-2.

Open with DEXTER

For large Mach numbers , the stand-off distance approaches the geometrical radius of the object. For lower , as indicated by the coloured lines, all simulation series confirm the scaling law of RSO/R ∝ ℳ-2 as previously determined in Fig. 7.

In Fig. 13 top panel, the correlation of the stand-off distance with the non-linearity parameter η is shown for different values of γ.

thumbnail Fig. 13

Correlation of the shock’s stand-off distance with the non-linearity 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 log-log plane of the plot. Since we found RSO/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 stand-off distance is shown in units of the accretion radius as a function of the Mach number. Here too the previously determined correlation between the stand-off 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 stand-off distance: (33)This relation is shown in comparison to a simulation series of varying adiabatic index in Fig. 14.

thumbnail Fig. 14

Stand-off distance in units of accretion radius RSO/RA 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 stand-off 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 RSO/RARSO/RA + b(34)with a and b as free fitting parameters. A least-square 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.

thumbnail Fig. 15

Comparison of the semi-analytical 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 stand-off distance with the non-linearity 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 stand-off distance, RSO, 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 stand-off distance. As shown above, the shock’s stand-off distance RSO can be approximated by (37)with g(γ) given by Eq. (35).

Strictly speaking, the formula above for the stand-off distance is only valid for a stand-off distance larger than the geometrical radius of the object. In the regime of large, low-mass bodies moving with high velocity, the stand-off distance approaches the value of the geometrical radius instead, and smin 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 non-linearity parameter η = [2 (ℳ2−1)/ℳ2R/RA] -1. The result of this comparison is shown as the ratio of the force in the experiments and the approximate formula in Fig. 16.

thumbnail Fig. 16

Comparison of the newly derived approximate formula for dynamical friction FDF with data from various numerical experiments . The corresponding ratio is given as function of the non-linearity 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 stand-off 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.

thumbnail 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, FDF, 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 CA = ln(smax/RSO) with RSO according to Eq. (32). For the maximum radius in CA we used rmax = 1000R. As can be seen, our approximative function for CA 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 FDF 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 stand-off 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 stand-off distance RSO 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 stand-off 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 stand-off distance does a wake of higher density form behind the moving body and induce a gravitational drag force.

Hence, in the well-known formulae by Ruderman & Spiegel (1971) or Ostriker (1999), it is the stand-off distance RSO of the shock that determines the minimum length scale smin required in the so-called Coulomb logarithm CA = ln(smax/smin): (39)The formation of a spherically symmetric atmosphere (with radius RSO) around a gravitating object implies that the hydrodynamical drag force vanishes in contrast to the non-gravitating case, but for very high Mach numbers the shock moves very close to the object so that RSOR. In this case, no spherically symmetric envelope around the object forms and the hydrodynamical drag becomes important again.

7.3. Relating the stand-off distance to the accretion radius

In a second step, to allow an easy computation of the stand-off distance for general dynamical friction problems, we derived the relation of the shock’s stand-off distance to the common definition of the accretion radius, (40)depending only on the mass M and velocity V of the perturber. The stand-off 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 stand-off distance (RSOR).

Via an analytical estimate of the stand-off distance and further fitting of the remaining free parameters, we gave an approximate relation of the dependence of the shock’s stand-off 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 stand-off distance to give an update of the well-known 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 smax denotes the extent of the wake and RSO can be obtained using the derived formulae above.

Acknowledgments

We thank Chris Ormel for very helpful discussions. R.K. acknowledges financial support from the Emmy-Noether-Program of the German Research Foundation (DFG) under grant No. KU 2849/3-1. 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 Baden-Wü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

  1. Armitage, P. J., & Natarajan, P. 2005, ApJ, 634, 921 [NASA ADS] [CrossRef] [Google Scholar]
  2. Billig, F. S. 1967, J. Spacecraft Rockets, 4, 822 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [NASA ADS] [CrossRef] [Google Scholar]
  4. Cantó, J., Raga, A. C., Esquivel, A., & Sánchez-Salcedo, F. J. 2011, MNRAS, 418, 1238 [NASA ADS] [CrossRef] [Google Scholar]
  5. Chandrasekhar, S. 1943, ApJ, 97, 255 [NASA ADS] [CrossRef] [Google Scholar]
  6. Chavarría, L., Mardones, D., Garay, G., et al. 2010, ApJ, 710, 583 [NASA ADS] [CrossRef] [Google Scholar]
  7. Edgar, R. 2004, New Astron. Rev., 48, 843 [NASA ADS] [CrossRef] [Google Scholar]
  8. Foglizzo, T., Galletti, P., & Ruffert, M. 2005, A&A, 435, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Grishin, E., & Perets, H. B. 2015, ApJ, 811, 54 [NASA ADS] [CrossRef] [Google Scholar]
  10. Hunt, R. 1971, MNRAS, 154, 141 [NASA ADS] [Google Scholar]
  11. Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  12. Kim, H., & Kim, W.-T. 2009, ApJ, 703, 1278 [NASA ADS] [CrossRef] [Google Scholar]
  13. Landau, L. D., & Lifshitz, E. M. 1966, Hydrodynamik (Berlin: Akademie-Verlag) [Google Scholar]
  14. Lee, A. T., & Stahler, S. W. 2014, A&A, 561, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Meyer, D. M.-A., Gvaramadze, V. V., Langer, N., et al. 2014a, MNRAS, 439, L41 [NASA ADS] [CrossRef] [Google Scholar]
  16. Meyer, D. M.-A., Mackey, J., Langer, N., et al. 2014b, MNRAS, 444, 2754 [NASA ADS] [CrossRef] [Google Scholar]
  17. Meyer, D. M.-A., Langer, N., Mackey, J., Velázquez, P. F., & Gusdorf, A. 2015, MNRAS, 450, 3080 [NASA ADS] [CrossRef] [Google Scholar]
  18. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
  19. Narayan, R. 2000, ApJ, 536, 663 [NASA ADS] [CrossRef] [Google Scholar]
  20. Ostriker, E. C. 1999, ApJ, 513, 252 [NASA ADS] [CrossRef] [Google Scholar]
  21. Rein, H. 2012, MNRAS, 422, 3611 [NASA ADS] [CrossRef] [Google Scholar]
  22. Rephaeli, Y., & Salpeter, E. E. 1980, ApJ, 240, 20 [NASA ADS] [CrossRef] [Google Scholar]
  23. Ricker, P. M., & Taam, R. E. 2008, ApJ, 672, L41 [NASA ADS] [CrossRef] [Google Scholar]
  24. Ricker, P. M., & Taam, R. E. 2012, ApJ, 746, 74 [NASA ADS] [CrossRef] [Google Scholar]
  25. Ruderman, M. A., & Spiegel, E. A. 1971, ApJ, 165, 1 [NASA ADS] [CrossRef] [Google Scholar]
  26. Ruffert, M. 1994, ApJ, 427, 342 [NASA ADS] [CrossRef] [Google Scholar]
  27. Ruffert, M., & Arnett, D. 1994, ApJ, 427, 351 [NASA ADS] [CrossRef] [Google Scholar]
  28. Sánchez-Salcedo, F. J., & Brandenburg, A. 2001, MNRAS, 322, 67 [NASA ADS] [CrossRef] [Google Scholar]
  29. Shima, E., Matsuda, T., Takeda, H., & Sawada, K. 1985, MNRAS, 217, 367 [NASA ADS] [Google Scholar]
  30. Slavin, J. D., Jones, A. P., & Tielens, A. G. G. M. 2004, ApJ, 614, 796 [NASA ADS] [CrossRef] [Google Scholar]
  31. Teyssandier, J., Terquem, C., & Papaloizou, J. C. B. 2013, MNRAS, 428, 658 [NASA ADS] [CrossRef] [Google Scholar]
  32. van Dyke, M. 1982, An album of fluid motion (Stanford, CA: Parabolic Press) [Google Scholar]
  33. van Dyke, M. D., & Gordon, H. 1959, NASA, Technical Report R-1 [Google Scholar]
  34. van Marle, A. J., Meliani, Z., Keppens, R., & Decin, L. 2011, ApJ, 734, L26 [NASA ADS] [CrossRef] [Google Scholar]
  35. 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.

thumbnail 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

thumbnail 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 stand-off distance of xSO = 0.533 was used, in accordance with the final quasi-stationary state of the numerical experiment.

Open with DEXTER

As is clearly visible, within the stand-off 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 stand-off 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 r0 to r, yielding (A.3)Here the subscript 0 refers to the reference radius r0 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 pre-shock values (supersonic regime) and the index “2” to the post-shock (subsonic) regime. These are valid in a reference frame moving with the shock. In our case the shock is stationary in the co-moving frame of the body, and we can directly apply the jump conditions if we know the pre-conditions. The simulations show that the pre-shock values of density, pressure, and temperature at the stand-off radius RSO are just the prescribed inflow values (at ; cf. Fig. A.1).

To obtain the pre-shock Mach-number 1 at this radius we use the free-fall 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 xSO = RSO/RA. Equation (A.6)matches the simulation outcome quite well (cf. Fig. A.1).

To obtain the (hydrostatic) post-shock density stratification around the object we first use these normalizations to obtain the general relation (A.7)with an arbitrary point x0, where ρ0 and c0 are known values at the reference radius. If we now choose x0 = xSO, 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 shock-position xSO enters at this point as a free parameter obtained from the numerical solution.

Appendix B: Estimate of the stand-off distance

Because the stratification around the object is nearly hydrostatic, some limits can be obtained on xSO. Here, we use the strong shock conditions, i.e. 1 → ∞ and obtain for the jumps in density and pressure and hence for the post-shock soundspeed (in units of c) (B.3)We now assume that for large distances, x → ∞, the density of the post-shock stratification lies below the unperturbed density ρ(x → ∞) ≤ ρ as seen in the simulations. This is shown in Fig. A.2 where the analytical post-shock 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 xSO 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 RSO it nevertheless provides a reasonable upper limit for the stand-off distance, and follows the trend of the obtained numerical results.

Appendix C: Convergence study

We ran simulations using varying grid sizes from roughly 103 up to 106 grid cells in total (more precisely, we varied the number of grid cells from Nr × Nθ = 50 × 23 = 1150 up to Nr × Nθ = 1720 × 790 = 1.3588 × 106 grid cells). The inner radial boundary of these simulations was chosen as R = 0.01 R. Other simulation parameters correspond to the “fiducial set-up” given in Table 1. From the final quasi-stationary state we determine the two most relevant physical quantities of our investigation, namely the shock’s stand-off distance RSO and the dynamical friction acting on the moving body.

The results of this convergence study are presented in Fig. C.1.

thumbnail Fig. C.1

Convergence test of the numerical set-up. Values on the horizontal axis denote the total number of grid cells Nr × Nθ. The vertical axis shows the relative deviation of the resulting stand-off distance RSO (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 106 grid cells. The numerical results show a clear convergence trend. For numerical grids larger than 105 grid cells, the deviations become negligibly small. To be on the very safe side, we chose Nr × Nθ = 860 × 400 = 3.44 × 105 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 Nr × Nθ = 700 × 480 = 3.36 × 105 grid cells as our default resolution.

All Tables

Table 1

Overview of the series of simulations performed.

All Figures

thumbnail 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 z-direction. 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 iso-density contours. Bottom panel: laboratory data from van Dyke (1982, Fig. 266).

Open with DEXTER
In the text
thumbnail Fig. 2

Quantitative comparison of the shock’s stand-off 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 set-up. Blue dots denote early numerical results by van Dyke & Gordon (1959).

Open with DEXTER
In the text
thumbnail Fig. 3

Force contributions in the equilibrium state according to Eq. (17) for the non-gravitating, 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 Ftot. The curve describes a parabola with F ∝ ℳ2.

Open with DEXTER
In the text
thumbnail Fig. 4

Density and velocity structure of the final quasi-stationary state around a gravitating object moving with supersonic motion through a gaseous homogeneous medium. The body is moving towards the positive z-direction. Density is colour-coded and black lines denote iso-density contours. Velocities are given as arrows, scaled by the speed.

Open with DEXTER
In the text
thumbnail Fig. 5

Fractions of the gravitational drag force f(ri) 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 stand-off distance RSO of the shock ahead of the object, and the dashed vertical lines denote the accretion radius RA 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
thumbnail 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 Ftot.

Open with DEXTER
In the text
thumbnail Fig. 7

Correlation of the shock’s stand-off distance RSO with object velocity or Mach number (top panel) and with object mass M (bottom panel). Solid lines denote linear fits in the log-log plane.

Open with DEXTER
In the text
thumbnail Fig. 8

Correlation of the shock’s stand-off distance RSO with the non-linearity parameter η, as given in Eq. (26), for various model sequences. In the top panel, RSO is normalized by the object radius R. For η> 1, a dashed line representing the relation RSO/R = η is superimposed. In the bottom panel, RSO is normalized by the corresponding accretion radius.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 10

Correlation of the ratio of dynamical friction and Coulomb logarithm with the mass of the perturber.

Open with DEXTER
In the text
thumbnail Fig. 11

Correlation of the ratio of dynamical friction and Coulomb logarithm with the velocity of the perturber.

Open with DEXTER
In the text
thumbnail Fig. 12

Correlation of the shock’s stand-off 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 RSO/R ∝ ℳ-2.

Open with DEXTER
In the text
thumbnail Fig. 13

Correlation of the shock’s stand-off distance with the non-linearity 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
thumbnail Fig. 14

Stand-off distance in units of accretion radius RSO/RA 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
thumbnail Fig. 15

Comparison of the semi-analytical 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 stand-off distance with the non-linearity 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
thumbnail Fig. 16

Comparison of the newly derived approximate formula for dynamical friction FDF with data from various numerical experiments . The corresponding ratio is given as function of the non-linearity 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
thumbnail 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
thumbnail 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
thumbnail 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 stand-off distance of xSO = 0.533 was used, in accordance with the final quasi-stationary state of the numerical experiment.

Open with DEXTER
In the text
thumbnail Fig. C.1

Convergence test of the numerical set-up. Values on the horizontal axis denote the total number of grid cells Nr × Nθ. The vertical axis shows the relative deviation of the resulting stand-off distance RSO (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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.