Bondi–Hoyle–Lyttleton accretion ﬂow revisited: Analytic solution

The time-steady equation for a 1D wind accretion ﬂow, i.e. the Bondi–Hoyle–Lyttleton (BHL) equation, is investigated analytically. The BHL equation is well known to have inﬁnitely many solutions. Traditionally, the accretion radius has been assumed to be 2 GM /v 2 ∞ , but its mathematical foundation has not been clariﬁed because of the non-uniqueness of the solution. Here, we assume that the solution curves possess physically nice characteristics, i.e. velocity and line mass-density increase monotonically with radial distance. This condition restricts the accretion radius to the range ( 0 . 71 − 1 . 0 ) × 2 GM /v 2 ∞ . Further assumptions, speciﬁcally, that the solution curves for velocity and line mass-density are convex upward, restrict the accretion radius to ( 0 . 84 − 0 . 94 ) × 2 GM /v 2 ∞ , and 0 . 90 × 2 GM /v 2 ∞ , respectively. Therefore, we conclude that the accretion radius is almost uniquely determined to be 0 . 90 × 2 GM /v 2 ∞ . .


Introduction
Accretion flow is one of the main sources of energy for a close binary system. There are two kinds of accretion flows: Roshe lobe overflow and wind accretion flow. In the present paper, we consider the latter.
If a gravitating object moves through an interstellar gas cloud or, equivalently, the gas flows past a gravitating object, the gas is attracted by and forms an accretion column behind the object. We term this phenomenon wind accretion.
When the temperature of the gas in the accretion column is zero, the accretion column becomes an infinitely thin line, i.e. an accretion line. We call this accretion flow Bondi-Hoyle-Lyttleton (BHL) accretion flow after the pioneering studies performed by Bondi, Hoyle, and Lyttleton [1][2][3].
Because the accretion rate is important in an astrophysical context, they tried to estimate it intuitively. Consider a particle flying freely in the gravitational field produced by a gravitating object. By comparing the gravitational energy and the kinetic energy of the particle, it is accreted by the object if the following condition is satisfied: where r HL is the so-called Hoyle-Lyttleton accretion radius. The present derivation is rather naive and intuitive. See Sects. 2.1 and 2.2 for a more detailed discussion.
In the above discussion, particles are assumed to be collisionless in the symmetry axis. Bondi and Hoyle [3] derived the time-steady hydrodynamic equation of cold gas (T = 0). They tried to solve the equation analytically. However, it is now known that the equation does not possess a unique solution [4,5].
The effect of gas pressure was considered by Wolfson [6,7] and Yabushita [8,9]. In this case, the pressure thickens the accretion line so that it appears as an accretion column. Yabushita [8,9] claimed that the non-uniqueness of the solution can be removed by introducing pressure. However, Horedt [10] pointed out that the argument by Yabushita [8,9] was inadequate because of the inappropriate use of the Bernoulli equation. Horedt studied the equation and claimed that there are still infinitely many solutions, and therefore the problem of non-uniqueness of the solution cannot be resolved by introducing gas pressure.
Therefore, the problem of non-uniqueness of BHL flow solutions still remains. Determining a unique solution is necessary to estimate the mass accretion rate, which is important in the astrophysical context. The purpose of this paper is to analyze the BHL equation. In the present study, we investigate the time-steady BHL equation analytically and remove the non-uniqueness of the solution by assuming physically plausible conditions. Hunt [11,12] was the first to perform numerical hydrodynamic simulations of the wind accretion flow assuming an axisymmetric geometry. It was further elaborated by Shima et al. [13] and Koide et al. [14]. Ruffert [15][16][17][18] calculated 3D accretion flows numerically. Non-uniqueness of the solution does not exist in their time-dependent numerical simulations. However, they calculated "warm" gas flows, because the conventional numerical scheme can barely treat high Mach number flows with M > 10. Note that the pressureless flow corresponds to an infinitely large Mach number.
By performing a linear stability analysis, Cowie [20] proved that the solution of the steady flow of the cold gas in the symmetry axis is unstable. This means that there is no steady solution in an exact sense. Nevertheless, we expect that the time-steady solution represents somewhat of a time-averaged flow. In this sense, the present study is motivated by mathematical rather than physical interest.
In Sect. 2, we provide an outline of steady BHL flow. We present the basic equations in Sect. 3, and discuss solutions of the equations in Sect. 4. The removal of the non-uniqueness of solutions and a comparison of our results with others given in the literature are discussed in Sect. 5.

Bondi-Hoyle-Lyttleton accretion flow
The BHL flow is basically a 3D phenomenon. If the pressure in the accretion column is sufficiently small, the accretion column becomes very thin. If the pressure is zero, the accretion column can be approximated by an infinitely thin line, and a 1D approximation can be applied. We limit our discussion to this case.  and follow a ballistic orbit, which is hyperbolic. The trajectory is expressed as [19] (see also Ref. [5]) where b is the impact parameter of the particle, and (r, θ) the polar coordinates. On the symmetry axis behind the central object (θ = 0), we have The mass of gas flowing through the annular region (b, b + db) per unit time is This means that there is a constant mass loading on the unit length of the symmetry axis. The radial component of the velocity of gas particles just before hitting the symmetry axis is v ∞ regardless of the position on the symmetry axis [19] (see also Ref. [5]). These two characteristics are very remarkable and will be used in the following discussion.

Accretion line
Gas particles are assumed to lose their angular momentum about the object in the collision on the symmetry axis. In this process, part of the kinetic energy is transformed into thermal energy, which is assumed to be radiated away immediately. In this case, the temperature and the pressure of the gas on the symmetry axis are kept zero; an accretion line with infinite density is formed. As depicted in Fig. 1, a point called the stagnation point forms where the gas velocity is zero. The stagnation point is also the bifurcation point of the streamline: The gas left of the stagnation point has negative radial velocity and eventually accretes onto the object, and the gas right of the stagnation point has positive velocity and flows away from the object. Therefore, the position of the stagnation where v e and r are the escape velocity and the radial position on the accretion line, respectively. As stated above, the radial velocity of the incoming gas particle on the accretion line is v ∞ . If v ∞ > v e , the gas particle in the accretion line would escape from the object. Based on this argument, Hoyle and Lyttleton [1,2] estimated the position of the stagnation point as where r HL is known as the Hoyle-Lyttleton radius. We define the accretion radius r a as the critical impact parameter such that material with an impact parameter smaller than this value will be accreted.
When r = r s , then b = r a . Therefore, we have the following relation between the stagnation radius r s and the accretion radius r a :

Numerical simulations of collisionless particles
We perform a numerical simulation of the motion of collisionless particles in the accretion line. The method of calculation employs a particle scheme. The computational domain is a 1D space of ε < x ≤ 100, where x is the dimensionless radial distance defined as and ε is a small positive number. A gravitating object of unit mass is placed at the origin, x = 0.
Once the position x of a particle is less than ε, the particle is removed from the computation; the particle is also removed if its position is larger than 100.
Initially there are no particles in the computational domain. At each time step t, we inject particles uniformly into this domain; the starting positions of the injected particles are determined using a random number generator.
We define the dimensionless velocity of a gas particle as The initial dimensionless velocity of the injected particle is set to unity regardless of its initial position. We follow the trajectory of all the particles in the computational domain until a steady state is reached.
The dimensionless velocity and line mass-density of the gas are calculated as follows. We divide the domain by 10 000 cells. In a cell, there are a number of particles. We calculate a mean value from the particle velocities to define the gas velocity in the cell. Because the line mass-density is associated with the number of particles in a cell, we estimate the dimensionless line mass-density by 4  counting the number of particles in the cell and scaling properly as where is the line mass-density. In Fig. 2, the dimensionless velocity y and the dimensionless line mass-density z are plotted against the dimensionless radial position x. From the figure, the velocity y is 0 at x = 1, and approaches unity for larger x, signifying that the stagnation point is x = 1 as expected. We also note that the velocity y is a monotonically increasing function of x: and that y(x) is convex upward: The line mass-density z has the same characteristics: (2.14) and We assume that these characteristics constitute the guiding principle in looking for the unique solution of the BHL equation. This principle will be discussed later.

Basic equations 3.1. Time-dependent equations
In this section, we consider the time-dependent equations of collisional gas particles, i.e. the Euler equation. The mass-conservation equation is written as where m is the line mass-density. The momentum equation is The above equations are made dimensionless using Eqs. (2.9), (2.10), and (2.11), giving (3.4)

Time-steady equation
We omit the time-dependent term in (3.3) to obtain which can be easily integrated to give Equation (3.7) is the basic equation that we try to solve in the following sections. In order to constrain the solutions further, we derive using Eqs. (3.6) and (3.7) the following auxiliary equation: Moreover, the second-order derivatives of y 0 and z 0 are  . Solutions satisfying condition I lie between these two curves.

Argument by Edgar [5]
Solving Eq. (3.7) is not easy. Apart from the fact that there is no analytic solution, a numerical treatment is also not easy because of the presence of singular points. On physical grounds, we assume that the velocity far downstream approaches unity, i.e. y 0 → 1 as x → ∞. as condition I. Hereafter, we drop the subscript 0 for the sake of brevity. Because the stagnation point, x = a, is a singular point, a solution starting from the region x < a, y < 0 has to pass the point (a, 0) smoothly to connect to a solution in the region x > a, y > 0. Edgar proved that, for the existence of solutions under condition I, the position of the stagnation point cannot be arbitrarily chosen, but must satisfy the inequality a ≥ 0.5. Following Edgar's procedure, we show two typical solution curves in Fig. 3, where we assume that a = 0.8. The upper and lower solution curves shown in Fig. 3 each pass the critical points (x, y) = (2a, y 1 ) and (x, = (2a, y 2 ), respectively, where y 1 = 1 + √ 1 − 1/2a 2 and y 2 = 1 − √ 1 − 1/2a 2. These two critical points correspond to the extrema of the curve determined by Eq. (3.14): (3.14) All possible solutions of Eq. (3.7) under condition I lie between those two curves (see Ref. [5] for a detailed discussion). The integration procedure is described in Sect. 4.     Figure 4 shows the corresponding line mass-density z for Fig. 3. The density exhibits a curious behavior in that, with increasing x, it increases, decreases, and increases again. Condition I permits solutions with this feature, but they are physically unrealistic.

Further requirements
We further assume that the line mass-density increases monotonically: We refer to this requirement as condition II. We show an example of a solution satisfying this condition in Figs. 5 and 6. In the present example, we set a = 0.8. Nevertheless, there are still infinitely many solutions even if we specify a. The method of solution will be discussed in Sect. 4. As can be seen, both y and z increase monotonically. However, these curves show wavy behavior, which is not found in Fig. 2 Fig. 6. Line mass-density corresponding to Fig. 5 showing a monotonic increase as required. Note the wavy nature of the solution near x = 1, which is not seen in Fig. 2.  so that the velocity curve is convex upward. We refer to this requirement as condition III. In Figs. 7 and 8, we present solution curves for the velocity and line mass-density solved subject to condition III. Again we set a = 0.8. We find that the velocity curve shows a smooth increase as required. However, the line mass-density curve still shows a slight wavy nature.

9/14
Downloaded from https://academic.oup.com/ptep/article-abstract/2015/11/113E01/2606831 by guest on 26 July 2018 Fig. 9. Velocity curve for a solution obeying condition IV; the solution curve for the line mass-density is convex upward. We further require the second derivative of z with respect to x to be non-positive, referring to this requirement as condition IV. An example of a solution solved under condition IV is exhibited in Figs. 9 and 10; we assume a = 0.81. This solution curve does not show any wavy nature and seems to satisfy all conditions I, II, III, and IV. If condition IV is met, all other conditions I, II, and III are automatically satisfied. This point will be discussed later. Note that the stagnation point is not x = a = 1, which occurs in collisionless flows (Fig. 2).

Determination of the stagnation point
To obtain solutions satisfying conditions I, II, III, and/or IV, we integrate the basic equation (3.7) numerically using Matlab. The procedure is as follows: 1) Firstly we specify one of conditions I, II, III, and IV.
2) We assume that the stagnation point, x = a, is arbitrarily in the range 0.5 < a < 1.1. The lower bound a = 0.5 is the minimum value with which the solution satisfying condition I can exist. The upper bound is arbitrarily chosen, but includes the value a = 1, which was anticipated by Bondi-Hoyle-Lyttleton. 3) We arbitrarily set a starting point for the integration, which is assumed to be x i = 2.5 in this study. 4) We choose the starting velocity y i at x i = 2.5 in the range 0 < y i < 1. 5) From the starting point (x i , y i ), we integrate Eq. (3.7) with decreasing x until the point x = a + ε is reached. This point is very close to the singular point. Here ε is a sufficiently small positive number. 6) At each integration step, we check if the specified condition is satisfied. 7) If the specified condition is not met, we abort the calculation. 8) From the starting point, we integrate Eq. (3.7) with increasing x up to a sufficiently large value, e.g. x = 100, in the same manner as step 6. 9) If the specified condition is not met, we abort the calculation. 10) When the calculation is successfully finished without abortion, we save the pair (a, y i ). 11) We repeat the above steps 4-10 with another y i . 12) When all values y i are exhausted, we change the stagnation point a and repeat steps 4-11. 13) We then plot the region (a, y i ) in which the solutions satisfy the specified condition.
In regard to the range 0 < x < a, we can obtain a unique solution. Taylor expansion of Eq. (3.7) around x = a gives the form of the solution that passes the singular point smoothly: where O (x − a) 2 denotes a term of second order, which we neglect in the present calculation. We start the integration at x = a − ε with the starting value y i as up to a point, which is sufficiently close to the origin. In order to integrate Eq. (3.7) numerically we use ordinary differential equation solvers, ODE45 and ODE23s, of Matlab. With ODE45 we can obtain solutions that satisfy one of the above conditions I, II, and III. ODE45 is a standard method to solve ordinary differential equations in Matlab, the method in which the solution is determined by comparing two solutions obtained by the 4th and 5th-order Runge-Kutta schemes. However, we could not obtain rigorous results under condition IV, because this condition makes Eq. (3.7) too stiff in the vicinity of the singular point to use ODE45. We use ODE23s instead, which is more suited to solving stiff equations. Figure 11 shows the zones of the starting value y i , with which the solutions satisfy conditions I, II, III, and/or IV, against the location of the stagnation point a. The abscissa represents the values for the stagnation point a, and the ordinate represents the starting velocity y i at x = 2.5. In the figure, all points in each region enclosed by curves I, II, III, or IV represent the starting values of the solutions that satisfy conditions I, II, III, or IV, respectively. An enlarged view of the region for condition IV is given in Fig. 12. The four conditions restrict the range of a as follows: 1) Condition I: 0.5 < a, 2) Condition II: 0.6627 < a < 1, 3) Condition III: 0.6975 < a < 0.8842, 4) Condition IV: 0.8086 < a < 0.8137.
The first one was already given by Edgar [5].   11. Regions in which conditions I to IV are satisfied. The abscissa represents the stagnation radius, a, and the ordinate the dimensionless velocity at the starting point of numerical integration, which is set to 2.5 in this calculation. Figures 11 and 12 clearly show that condition IV is the most restrictive. If condition IV is satisfied, all other conditions are automatically satisfied. Therefore, we can conclude that the allowed range for a is 0.8086 < a < 0.8137, or approximately a = 0.81.

Comparison with other authors
We determined the position of the stagnation point as a = 0.81 for the BHL flow by assuming that the solution should have physically reasonable characteristics, i.e. conditions I, II, III, and/or IV.
As to the value of a, Horedt [10] compiled the results of various authors. His Fig. 1 shows the distribution of a in the (a, M) map, where M is the upstream Mach number. Note that his nondimensionalization is different from ours by a factor of 2, and his a = 2 corresponds to our a = 1. Our Mach number is infinite and is not included in his figure.
The largest Mach number quoted by Horedt is M = 10. Ruffert [15][16][17][18] gives a = 0.8 − 0.9, whereas Koide et al. [14] give a = 1.2, where our normalization is used. For M = 5, Shima et al. [13] give a = 0.9, whereas Koide et al. have a = 1.1. All these results are consistent with ours.  Comparison between an analytic solution and a numerical solution. A 1D simulation of the accretion flow for a polytropic gas with a specific heat ratio, 1.2, was performed, and the time-averaged solution between t = 20 and t = 30 was calculated. The upper two curves plot velocities calculated analytically, velocity (analytic), and numerically, velocity (numerical). The lower two curves give line mass-densities calculated analytically, line mass-density (analytic), and numerically, line mass-density (numerical). The unit of numerical line mass-density in this figure is the number of particles in a cell in the numerical simulation. The analytic curve is scaled to fit the numerical curve. The stagnation radius of the numerical solution is found to be x = 0.83. We present the analytic solution in which the stagnation radius is x = 0.81, as required from condition IV. Both solutions agree fairly well.
In the above we discuss the stagnation radius, r s = a. However, in the case of cold incoming gas, we can determine the accretion radius r a by Eq. (2.8), which is written in dimensionless form as r a = √ r s . (5.1) Our stagnation radius r s = 0.81 gives r a = 0.90. Therefore, we may conclude that the dimensional accretion radius is

Comparison with numerical results
Cowie [20] pointed out that the cold BHL flow is intrinsically unstable; therefore, we could not obtain time-steady solutions. In this sense, our present consideration of time-steady solutions for BHL flows is mathematically rather than physically acceptable. To obtain a smooth and stable numerical solution, we performed a numerical simulation of a warm gas with a specific heat ratio γ = 1.2 in an accretion column with constant thickness. The numerical method used was the molecular hydrodynamic method proposed by the present authors [21,22].
The numerical solution thus obtained is not unstable, but some time variability still exists. We compute a time average during 20 ≤ t ≤ 30 to obtain a smooth solution. The stagnation point thus obtained is x = 0.83. Figure 13 compares our analytic solution, in which we assume a = 0.81, as required by condition IV; y 0i is chosen so that the analytic solution best fits the numerical solution. The upper curves show the velocities, and the lower curves the line mass-densities.
As can be seen in the figure, both solutions show strong agreement, although there is a slight discrepancy in line mass-density around x = 1.5. The reason for the agreement is probably due to 13