A&A 441, 863-872 (2005)
DOI: 10.1051/0004-6361:20053091

Radiative effects in supersonic wind accretion onto gravitating objects

I. Kryukov1,2 - N. Pogorelov2,[*] U. Anzer3 - G. Börner3 - G. Bisnovatyi-Kogan4


1 - Institute for Problems in Mechanics, Russian Academy of Sciences, 101-1 Vernadskii Avenue, 119526 Moscow, Russia
2 - Institute of Geophysics and Planetary Physics, University of California, Riverside, CA 92521, USA
3 - Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany
4 - Space Research Institute, Russian Academy of Sciences, Profsoyuznaya St. 84/32, 117810 Moscow, Russia

Received 18 March 2005 / Accepted 24 June 2005

Abstract
We investigate the influence of radiative effects on supersonic wind accretion onto gravitating objects. The accreting matter is assumed to be optically thin. The physical mechanisms taken into account include cooling due to free-free and free-bound transitions, the Compton heating via X-ray scattering on electrons and the inverse Compton cooling in the regions where the temperature of the matter becomes sufficiently large to be able to transfer part of its internal energy to photons. A wide range of determining parameters was covered, including the values applicable to the Vela X-1 binary system, but our main emphasis is on the study of the effects of radiative processes on the behavior of accretion flows. It is shown that the applicability of polytropic accretion models is very limited and the actual accretion rate can be considerably lower than that provided by the Bondi-Hoyle-Lyttleton formula. The detailed consideration of the realistic radiative effects proved to be of great importance in our understanding of the accretion phenomenon, since they can substantially affect it both qualitatively and quantitatively.

Key words: stars: neutron - shock waves - methods: numerical - radiative transfer - accretion, accretion disks

1 Introduction

Accretion onto neutron stars and black holes is believed to produce the main energy supply in galactic X-ray sources. The investigation of the wind accretion onto a gravitating center, which is often called Bondi-Hoyle-Lyttleton (BHL) accretion, has long been the subject of interest in view of possible observations of black holes moving through the interstellar gas (Hoyle & Lyttleton 1939, 1940a,b,c). If one allows some oversimplification, this problem admits clear physical and mathematical formulations, which allows one to analyze in detail various specific phenomena intrinsic to the accretion process. The solutions of the wind accretion problems are expected to explain the observations of X-ray pulsars where magnetic neutron stars accrete from massive O/B companions (Benensohn et al. 1997). These sources are reported to exhibit fluctuations in their spin period (Boynton et al. 1984; Nagase et al. 1989; Chakrabarty et al. 1993), but they also show the presence of a long-term equilibrium rotational period (Börner et al. 1987; Nagase et al. 1989; Bisnovatyi-Kogan 1991; Anzer & Börner 1995).

For our investigation we use a gas dynamics approach, which was first applied to analyze the supersonic wind accretion of massive objects by Salpeter (1964), who showed the presence of standing shock fronts behind such objects, with a large pressure increase across the front. It was determined that the shock has a conical shape and the cone angle depends both on the Mach number of the flow and the polytropic index of the matter. If the gas consists of noninteracting particles (zero pressure approximation) with a number density $n_\infty$ and a velocity $V_\infty$ at infinity, then the accretion rate is given by the formula

\begin{displaymath}\dot{M}=\pi R_{\rm HL}^2 m_{\rm p} n_\infty V_\infty;\quad R_{\rm HL}=\frac{2GM}{V_\infty^2}, \end{displaymath} (1)

where M and G are the mass of the accreting star and the gravitational constant, respectively; and $m_{\rm p}$ is the proton mass (for a pure hydrogen gas). The quantity $R_{\rm HL}$ has the dimension of length and is called the Hoyle-Lyttleton radius. If the flow is axially symmetric, only particles with impact distances smaller than $R_{\rm HL}$ can be accreted. If hydrodynamic effects are taken into account, this formula is invalid, and $\dot{M}$ should depend on the unknown function of the polytropic index $\gamma$, as shown by Bisnovatyi-Kogan et al. (1979) who constructed self-similar solutions in the framework of the gas dynamical model. Due to the complexity of the problem, a series of numerical studies has been undertaken (Hunt 1971, 1979; Livio et al. 1979, 1986a,b, 1991; Okuda 1983; Shima et al. 1986; Matsuda et al. 1987, 1991, 1992). Fryxell & Taam (1988), Benensohn et al. (1997), and Pogorelov et al. (2000) studied both accretion from a uniform and nonuniform wind at infinity. In a number of these papers, the plane and three-dimensional solutions turned out to be highly unsteady, exhibiting violent oscillations of the accretion cone and its tip attached to the accreting boundary. Since the stability of the solutions depends on $\gamma$ and is basically non-axisymmetric, a detailed analysis has been performed for a number of polytropic indices, e.g., for isothermal accretion (Shima et al. 1998), and in three dimensions (Ishii et al. 1993; Ruffert 1994a,b, 1995, 1996). Boffin & Anzer (1994) applied the Smoothed Particle Hydrodynamics (SPH) method to this problem. These SPH calculations gave results which were similar to those obtained with the finite difference methods, but they produced a much lower degree of instability.

Summarizing all these results one can say that solutions of the axisymmetric problem are generally stable for most parameters of the uniform wind. Non-axisymmetric calculations can be divided into two groups. One of them simply neglects the dependence on one of the Cartesian coordinates, thus reducing the problem to the planar case (note that a gravitation source remains spherically symmetric in this case). The second, truly three-dimensional, approach requires much more computer resources and can hardly be expected to incorporate such fine grids as are possible in the planar calculations. Most of the planar calculations of supersonic wind accretion onto a spherically symmetric gravitating center exhibit oscillations. The oscillation amplitude can be extremely large. This instability is usually accompanied by a sporadic increase in the angular momentum of matter, and accretion temporarily stops. Note that the instabilities can originate spontaneously from an initially symmetric flow. Possible physical mechanisms of such a behavior have been discussed by Foglizzo et al. (2005). It should be mentioned, however, that the flip-flop instability turned out to be incapable of accounting for the periodic phenomena observed in some transient Be-type X-ray binaries, since the amplitude of fluctuations observed in the planar case is not large enough (Shima et al. 1998). Moreover, three-dimensional computations (Ruffert 1996) give amplitudes which are 3 to 10 times smaller than those observed in the planar case.

To clarify the stability issue, Foglizzo & Ruffert (1997, 1999) performed a local stability analysis of the BHL accretion and determined the growth rates of both Rayleigh-Taylor and Kelvin-Helmholtz instability modes caused by the entropy gradient generated at the shock surface. Their results convincingly show that the presence of a shock is essential for obtaining unstable flows. The overall instability is greater if this shock is detached and for stronger shocks. A WKB analysis used in those papers is only marginally applicable, since the growth time is comparable with the advection time along streamlines that cross the shock and terminate at the accretor boundary. Decreasing the accretor size, if it is considerably less than the bow shock stand-off distance, only slightly affects the stability. Foglizzo & Tagger (2000) and Foglizzo (2001, 2002) discovered an entropic-acoustic (or advective-acoustic) cycle that, in the linear approximation, feeds the above-mentioned instability. It is important that a necessary condition for the possibility of such a cycle is the acceleration of the accretion flow behind the shock to supersonic velocities. Foglizzo et al. (2005) systematized numerical results for BHL accretion onto gravitating objects and demonstrated the relevance of the advective-acoustic instabilities for cases with detached bow shocks. It was determined that the accretion pattern must be stable in the limit of a weak shock and that it must be unstable in three dimensions for polytropic indices about the adiabatic value 5/3, Mach numbers greater than 3, and an accretor size sufficiently small to allow the development of perturbations.

As shown by Kryukov et al. (2003), the use of the polytropic approximation (Kryukov et al. 2000) can considerably distort the scenario of accretion onto stellar magnetospheres. By decreasing the polytropic index below its adiabatic value 5/3, we can only approximately allow for cooling effects, provided that the flow itself is very uniform and space variations in cooling intensity are not great. In fact, by specifying the value of $\gamma<5/3$, we automatically prescribe energy losses at each point. Since accretion transforms a part of the consumed kinetic energy into radiation, one can expect related heating effects which ultimately lead to a Comptonization of the accretion flow behind the shock near the surface of the accreting boundary.

The effect of radiation pressure becomes important for star luminosities of the order of the Eddington luminosity

\begin{displaymath}L_{\rm Edd}=4\pi\frac{GMcm_{\rm p}}{\sigma_{\rm T}}, \end{displaymath} (2)

where c, $m_{\rm p}$, and $\sigma_{\rm T}$ are the speed of light, the proton mass, and the Thompson cross-section of electron scattering. The effect of the radiation drag on the Hoyle-Lyttleton accretion was analyzed by Nio et al. (1998). Blondin et al. (1990, 1991) investigated a two-dimensional gas flow in the orbital plane of a massive X-ray binary system, where the mass accretion was fuelled by a radiationally-driven wind from a companion star, and found time-dependent accretion wakes containing dense filaments of compressed gas. On the other hand, for high X-ray luminosities ( $L \geq 10^{37} ~{\rm erg}~ {\rm s}^{-1}$) and due to Comptonization effects, the accretion flow was shown to be relatively steady. Stellar wind accretion calculations taking into account radiative processes were performed by Taam et al. (1991).

In this paper, we analyze the accretion flow from a uniform supersonic wind onto a gravitating center. The physical mechanisms taken into account include cooling due to free-free and free-bound transitions, the Compton heating via X-ray scattering on electrons and the inverse Compton cooling in the regions where the temperature of the matter becomes sufficiently large to be able to transfer part of its internal energy to photons. In contrast to the parametric analysis by Taam et al. (1991), we choose the wind parameters in a range that approximately corresponds to the Vela X-1 binary, according to the observations by Kreykenbohm et al. (1999). On the basis of these data, we derive proper parameters that determine the Compton heating rate and accretion efficiency. To facilitate the analysis of the radiative effects, we solve the problem in a rather simplified statement where the flow from the supergiant star is uniform at about 20-30 $R_{\rm HL}$ from the accreting star. The assumption of a uniform flow is only approximately valid for the binary systems considered. But as long as the accretion radius is sufficiently smaller than the dimensions of the binary system this assumption is reasonable. One finds that for the Vela X-1 binary system this is the case, since the accretion radius is around $3\times 10^{11}\ {\rm cm}$ and the binary separation is $3.6\times 10^{12}\ {\rm cm}$. For the accretion flow, we take a basic velocity of $400\ {\rm km}~{\rm s}^{-1}$, which results from an orbital velocity of $280\ {\rm km}~{\rm s}^{-1}$ and a wind speed of $300\ {\rm km}~{\rm s}^{-1}$ (for the justification of reduced wind speeds, see Feldmeier et al. 1996). Then the Bondi-Hoyle-Lyttleton accretion rates were taken between 10-10 and $10^{-8}\ M_{\odot}~{\rm yr}^{-1}$. We further assume that the chosen flow rate is constant in time and ignore any possible back reaction due to X-ray heating of the stellar atmosphere. One also finds that the stellar radius in such massive binaries is much smaller than the Roche-lobe surface, therefore neglecting tidal effects may also be justifiable.

The paper is structured as follows. Section 1 gives an introduction where we summarize a number of important results related to the supersonic wind accretion onto gravitating objects. In Sect. 2, we describe the physical and mathematical formulations of the problem of the supersonic wind accretion onto gravitating objects with radiation effects taken into account in the approximation of an optically thin medium. The details of the numerical method are briefly addressed. Section 3 contains the results of our numerical modeling. Finally, we discuss the importance of the obtained results for a better understanding of the accretion process in widely separated binaries.

2 Description of the problem

We shall perform our calculations on the basis of the Euler gas dynamics equations including the source terms responsible for the gravitational attraction from the star, the radiative pressure, and the cooling and heating terms. We assume that the matter is ideal and perfect, with the caloric equation of state in the form

 \begin{displaymath}\varepsilon={p\over (\gamma-1)\rho}, \end{displaymath} (3)

where $\varepsilon$, p, $\rho$, and $\gamma=5/3$ are the specific internal energy, the pressure, the density, and the specific heat ratio, respectively. The thermal equation of state is

 \begin{displaymath}p=2\rho {kT\over m_{\rm p}}, \end{displaymath} (4)

where T and k are the temperature and the Boltzmann constant. The coefficient 2 is chosen to take into account the pressures of electrons and protons, the temperatures of which are assumed to be the same in the, on average neutral, hydrogen plasma.


  \begin{figure} \par\includegraphics[width=8cm,clip]{3091f1.eps} \end{figure} Figure 1: Computational grid around the accreting center.
Open with DEXTER

The problem is solved for axisymmetric configurations. The symmetry axis, 0z, is aligned with the velocity vector at infinity. The governing equations are solved in a Cartesian coordinate system using the nonuniform mesh which is shown in Fig. 1. The grid is constructed as follows. One of the two families of grid lines is represented by radial rays $\theta={\rm const}$. The second family is given by a super-elliptic function

 \begin{displaymath}\left({x\over a}\right)^{2n}+\left({z\over b}\right)^{2n}=1. \end{displaymath} (5)

Here we assume a=b=R. To implement a smooth transition from a circle at R=R*, where R* is the radius of the inner boundary, to a rectangle at $R=R_\infty$, we use a variable exponent in Eq. (5):

 \begin{displaymath}n(R)=1+{29\over R_\infty -R_*}(R-R_*), \end{displaymath} (6)

where R is the distance from the star to the current point of the computational region. The mesh is refined for both $R\to R_*$ and $\theta\to 0$ using the formula

\begin{eqnarray*}R_k = R_* + (R_{\infty}-R_*)\times\biggl[\frac{k - 1}{N} + \lef... ...ac{1}{N}\right) \frac{\sin(\pi (k - 1)/N)}{\sin(\pi/N)}\biggr], \end{eqnarray*}


where k is in the range between 1 and N+1, N is the number of intervals, and $\Delta R_{\rm min}$ is the minimum cell size chosen. In our calculations, the minimum radial and angular cell sizes were 0.01R* and 0.005R*, respectively. In Fig. 1, we show a sample $80\times 60$ mesh with $R_\infty=500R_*$. From now on quantities at the inner and outer boundaries will be denoted by the subscripts * and $\infty$, respectively.

The dimensional form of the system of governing equations in the coordinate system (x,z) shown in Fig. 1 is the following:

 \begin{displaymath}{\partial x\vec{U}\over\partial t} + {\partial x\vec{E}\over\... ...\vec{H}_{\rm 2}-x\vec{H}_{\rm C} +x\vec{H}_{\rm cool}=\vec{0}, \end{displaymath} (7)

where

 \begin{displaymath}\vec{U}=[\rho,~ \rho u,~ \rho w,~ e]^{\rm T}, \end{displaymath} (8)


 \begin{displaymath}\vec{E}=[\rho u,~ \rho u^2+p,~ \rho u w,~ (e+p)u]^{\rm T}, \end{displaymath} (9)


 \begin{displaymath}\vec{G}=[\rho w,~ \rho uw,~ \rho w^2+p,~ (e+p)w]^{\rm T}, \end{displaymath} (10)


 \begin{displaymath}\vec{H}_1= \left[\begin{array}{c} 0,~ -p,~ 0,~ 0 \end{array}\right]^{\rm T}, \end{displaymath} (11)


 \begin{displaymath}\vec{H}_2=\left(1-{L\over L_{\rm Edd}}\right) {\rho~ G{\rm M}... ...a\\ \sin\theta\\ u\sin\theta+w\cos\theta \end{array}\right], \end{displaymath} (12)


 \begin{displaymath}\vec{H}_{\rm C}=[0,0,0,H]^{\rm T},\quad H={k\sigma_{\rm T}\do... ...\rm p} m_{\rm e}}~ E_{\rm eff}~ {\rho (T_{\rm C}-T)\over R^2}, \end{displaymath} (13)


 \begin{displaymath}\vec{H}_{\rm cool}=[0,0,0,C]^{\rm T},\quad C=n^2\Lambda(T)={\rho^2\over m^2_{\rm p}}~ \Lambda(T). \end{displaymath} (14)

In these formulae $\vec{v}=[u,~v,~w]^{\rm T}$ is the velocity vector and $e=p/(\gamma -1) + \rho \vert\vec{v}\vert^2/2$ is the total energy per unit volume. $\vec{H}_1$ is a geometrical source term that appears since in the axisymmetric case we can exactly calculate the partial derivative $\partial/\partial y$ in the three-dimensional system of the Euler equations (the y-axis is perpendicular to the plane shown in Fig. 1). The given form of this term is valid in the half-plane $x\geq 0,\ y=0$.

The vector $\vec{H}_2$ describes the action of gravity in the presence of the radiation pressure force in the momentum equations and its work in the energy equation. Here we assumed that the radiation flux is isotropic in space. For this reason the total luminosity L is connected with the spectrum averaged luminosity  $\mathcal{L}$ per unit solid angle $\Omega$ by the formula

\begin{displaymath}L=\int_{4\pi} \mathcal{L} {\rm d}\Omega=2\pi\int^\pi_0 \mathcal{L}(\theta)\sin{\theta}{\rm d}\theta=4\pi \mathcal{L}. \end{displaymath}

We assume that

 \begin{displaymath}L=E_{\rm eff} \dot{M}c^2, \end{displaymath} (15)

where $E_{\rm eff}$ is the efficiency factor. For the case of the accretion onto black holes this factor is rather uncertain, but for the accretion onto a neutron star we have

 \begin{displaymath}E_{\rm eff}={GM\over Rc^2}\cdot \end{displaymath} (16)

For Vela X-1 van Kerkwijk et al. (1995) give $M=1.9~M_\odot$. Assuming $R=10^6\ {\rm cm}$, we obtain $E_{\rm eff}\approx 0.25$.

Formula (13) describes the heating/cooling due to the Compton effect (Levich & Sunyaev 1971). The value of the Compton temperature can be obtained (see Levich & Sunyaev 1971) from the formula

 \begin{displaymath}T_{\rm C}={1\over 4k}{\int_0^\infty \epsilon E(\epsilon){\rm d}\epsilon\over\int_0^\infty E(\epsilon){\rm d}\epsilon}, \end{displaymath} (17)

where the photon radiation energy density E depends on the photon energy $\epsilon=h\nu$. Kreykenbohm et al. (1999) analyzed observational photon spectra for the X-ray pulsar Vela X-1 obtained with the RTXE (Rossi X-ray Timing Explorer) in February, 1996. None of these spectra could be fitted with thermal bremsstrahlung or black body spectra and they used a power-law approximation modified by a so-called Fermi-Dirac cut-off:

\begin{displaymath}C(\epsilon)= \frac {\epsilon^{-\alpha}} { \exp^ { {\epsilon-\epsilon_{\rm cut}\over \epsilon_{\rm fold}} } +1 }\cdot \end{displaymath} (18)

The fitting parameters $\alpha$, $\epsilon_{\rm cut}$, and $\epsilon_{\rm fold}$ are summarized in Table 3 of their paper. As $E(\epsilon)=\epsilon C(\epsilon)$, we readily obtain $T_{\rm C}=4.6$, 3.0, 2.9, and $6.1\times 10^7~{\rm K}$ for observations 1-4. This means that the observation-fitted Compton temperature in all these cases is systematically greater than the value $\epsilon_{\rm m}/4k$, which is usually taken to characterize $T_{\rm C}$(see, e.g., Bisnovatyi-Kogan & Blinnikov 1980). The actual value of the Compton temperature is determined by the physical processes near the surface of the neutron star. In our calculations we use an effective constant Compton temperature, which is based on observations of Kreykenbohm et al. (1999) and is expected to take into account those processes. Thus, only the luminosity of the neutron star varies with the accretion rate. We calculated models for all the Compton temperatures quoted above.

The cooling term due to radiative losses is represented by Eq. (14), where $\Lambda(T)$ gives the cooling function, which is usually calculated for a low-density, optically thin gas of some chosen cosmic element abundances (see, e.g., Raymond et al. 1976). We shall give the expression approximating this function later.

Let us normalize the quantities of density, velocity, pressure and temperature to $\rho_\infty$, $V_\infty$, $\rho_\infty V_\infty^2$ and $T_\infty$, respectively. Since a wind accretion flow has its internal length scale, e.g., $R_{\rm HL}=2GM/V^2_\infty$, the time scale becomes $2GM/V^3_\infty$. Thus, the dimensionless system acquires a form similar to Eq. (7), where the structure of $\vec{U}$, $\vec{E}$, $\vec{G}$, and $\vec{H}_1$ remains unchanged. The other source terms become

 \begin{displaymath}\vec{H}_2=\left(1-{L\over L_{\rm Edd}}\right) {\rho \over 2R^... ...a\\ \sin\theta\\ u\sin\theta+w\cos\theta \end{array}\right], \end{displaymath} (19)


 \begin{displaymath}H=S_{\rm C}E_{\rm eff} \dot{\mathcal{M}}{\rho (T_{\rm C}-T)\over R^2},\quad C=S_{\rm b}\rho^2\bar{\Lambda}(T), \end{displaymath} (20)

where we used the Bondi-Hoyle-Lyttleton accretion rate (Lipunov 1992) to introduce

\begin{displaymath}\dot{\mathcal{M}}=\dot{M}/\dot{M}_{\rm BHL},\quad \dot{M}_{\rm BHL}=4\pi(GM)^2(a_\infty^2+V^2_\infty)^{-3/2}\rho_\infty. \end{displaymath} (21)

The dimensional cooling function is normalized here by some characteristic value $[\Lambda]$.

Formulae (19)-(20) involve two additional dimensionless parameters:

 \begin{displaymath}S_{\rm C}={T_\infty\over 2GMV_\infty} {k\sigma_{\rm T}\dot{M}... ...}={2GM \rho_\infty [\Lambda]\over V_\infty^5 m^2_{\rm p}}\cdot \end{displaymath} (22)

The parameter $S_{\rm C}$ measures the ratio of the Compton heating rate to the rate of potential energy gain of the accretion flow; the other parameter $S_{\rm b}$quantifies the ratio of the radiative and convective transport cooling.

The equation of state (3) preserves its form, whereas the dimensionless thermal equation of state becomes

 \begin{displaymath}p=\frac{\rho T}{\gamma M_\infty^2},\quad M^2_\infty=\frac{V_\... ...\infty^2},\quad a_\infty^2=\frac{\gamma p_\infty}{\rho_\infty} \end{displaymath} (23)

where $M_\infty$ and $a_\infty$ are the Mach number and the speed of sound in the supersonic stellar wind.

All quantities at the entrance segments of the outer boundary should be specified. In dimensionless form, we have $\rho_\infty=1$, $T_\infty=1$, $p_\infty=1/\gamma M^2_\infty$, u=0, and w=1. In fact, if we specify dimensional values of $V_\infty$ and $T_\infty$ and then find  $\rho_\infty$by choosing $\dot{M}_{\rm BHL}$, all parameters of the problem will be fixed.

Because of the symmetry of the problem, it clearly suffices to seek the solution only in the half-plane $x\geq 0$ of the Cartesian system. As conditions on the inner boundary of the computational region we use absorbing boundary conditions suggested by Pogorelov & Semenov (1997). They consist in the extrapolation of data on supersonic segments of the boundary and in applying relations in a centered rarefaction wave to accelerate the flow to a sonic velocity on subsonic segments.

The numerical solution of system (7) is performed with the help of the hydrodynamic package MHD-E-NS2D, developed for the simulation of compressible gas and plasma flows governed by ideal magnetohydrodynamic (MHD) and purely gas dynamic equations and the Navier-Stokes equations (Ivanov & Kryukov 1996; Pogorelov & Matsuda 1998; Ivanov et al. 1999). The package is characterized by a developed modular structure that allows us to obtain the required resolution in time and space and add different physical processes by combining the proper program modules. It is substantially based on numerical schemes that ensure high (second or third) order of accuracy and preserve the monotonicity of functions at discontinuities. Most of the incorporated methods can be classified as Godunov schemes and can be attributed to TVD (total variation diminishing), ENO (essentially nonoscillatory) or WENO (weighted ENO) classes, depending on the implementation of the code modules. To determine inviscid fluxes through the computational cell interfaces, either the exact or some of the approximate solutions to the Riemann problem of the disintegration of an arbitrary discontinuity are used (Kulikovskii et al. 2001). The methods are implemented on both stationary and moving structured adaptive grids. The numerical results presented below were obtained using substantially two-dimensional reconstruction procedures for determining the values of the vector $\vec{U}$ on the right- and left-hand sides of the computational cell interfaces. The time integration is performed by third- or higher-order Runge-Kutta methods, whose application allows one to improve the overall stability of the numerical method.

3 Numerical results

The code has been tested by Kryukov et al. (2003) performing a comparison with solutions of the one-dimensional spherically symmetric accretion with radiative effects (Bisnovatyi-Kogan & Blinnikov 1980) and supersonic, polytropic wind accretion onto a gravitating center (Pogorelov et al. 2000), and a good agreement has been obtained. It turned out that an implicit approximation of the source term allowed us to perform calculations with the time step determined by the stability condition of the discretized differential operator in a wide range of the dimensionless parameters $S_{\rm C}$ and $S_{\rm b}$ (in all calculations the Courant-Friedrichs-Levy stability parameter was 0.5). The choice of the initial values turned out to be of particular importance for the calculations described below; e.g., if we start with a uniform flow at large distances from the accretor, a thin low-temperature column forms in the tail region around the symmetry axis ( $T\sim 10^2\ {\rm K}$ for $T_{\rm C}=10\ {\rm keV}$). Since the Mach number in this case becomes extraordinary large (about 104), calculations with the help of an explicit scheme become impractical. This can be avoided if one starts with the solution that corresponds to accretion with only bremsstrahlung cooling. It is important to note in this connection that the restriction of the temperature values from below, as was done by Taam et al. (1991), where only values greater than $5\times 10^4$ K were allowed, can cause substantial perturbations in the flow field that result in the development of highly unsteady solutions.

Table 1: List of calculated variants.

In this paper we follow Buff & McCray (1974), who noted that if stellar radiation substantially increases the ionization degree of trace elements by the photoionization process, the radiative cooling rate can be lower than in the purely collisional case. For this reason, an appropriate cooling function (Stellingwerf & Buff 1982) has been chosen:

 
$\displaystyle \Lambda(T)=\biggl( \left[2.4\times 10^{-27} T^{1/2}+6\times 10^{-22} T^{-1/2}\right]^{-1} +2.5\times 10^{75} T^{-12} \biggr)^{-1}.$     (24)

This function obviously degenerates into the free-free, bremsstrahlung cooling if $T\to \infty$.

Prior to proceeding with a detailed description of the numerical results, we note that very long-time calculations were usually performed, with at least $6\times 10^5$ time steps. This resulted from the large difference in computational cell sizes in the vicinity of and far away from the accretor. A solution was considered as converged to a steady state if the accretion rate $\dot{M}$ had not changed more than 0.01% within 1000 consecutive time steps, which corresponds to 6 min physical time. It should be noted that a general structure of the flow at distances of about $1\ R_{\rm HL}$ turns out to be established after 10 h (compare with the Vela X-1 binary period 8.964 days given by Kreykenbohm et al. 1999). This means that the accretion time scale is considerably less than the rotation period.

We performed a rather detailed parametric investigation of the possible accretion patterns. Some of the variants calculated are summarized in Table 1, where the right two columns represent the results of our calculations. Here accretion rates are given in $M_{\odot} ~ {\rm yr}^{-1}$, velocities in ${\rm km}~{\rm s}^{-1}$, densities in ${\rm g}~{\rm cm}^{-3}$, and the star luminosity L in ${\rm erg}~{\rm s}^{-1}$. All other parameters are dimensionless. Variant 2 (the base variant) is the one which is closest to the parameters of the Vela X-1 binary. We fixed $T_\infty=10^4$ K and varied the following quantities: $V_\infty$, $T_{\rm C}$, $\rho_\infty$ (via $\dot{M}_{\rm BHL}$, since $T_\infty$ was fixed), $E_{\rm eff}$, and the form of the cooling function. The effect of different values for $V_\infty$, $T_{\rm C}$, and $\dot{M}_{\rm BHL}$ can be seen by comparing Variants 2-4; 2, 5, and 6; and 2, 7, and 8, respectively. The abbreviations BR and SB refer to the pure bremsstrahlung and the Stellingwerf-Buff cooling function. Variants 9 and 10 have the lowest and the highest calculated accreting rates. Variants 11-14 correspond to the accretion efficiency equal to 0.1, which is frequently used in calculations of the accretion onto neutron stars (Bisnovatyi-Kogan & Blinnikov 1980). The size of the computational region was chosen so that radiative effects near the left (upstream) boundary were negligible. The last two columns of the table indicate the calculated accretion rate (in terms of the Bondi-Hoyle-Lyttleton rate) and X-ray luminosity of the neutron star.


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f2.eps} \end{figure} Figure 2: Contour lines of the Mach number and the temperature logarithm. Variant 1.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f3.eps} \end{figure} Figure 3: Contour lines of the Mach number and the temperature logarithm. Variant 2.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f4.eps} \end{figure} Figure 4: Contour lines of the Mach number and the temperature logarithm. Variant 3.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f5.eps} \end{figure} Figure 5: Contour lines of the Mach number and the temperature logarithm. Variant 4.
Open with DEXTER

In Figs. 2-5 we show the distributions of the Mach number (above the symmetry axis) and temperature logarithm (below the symmetry axis) for Variants 1-4. All contour plots are shown for a fixed region of $-20\leq z\leq 20$ and $0\leq x\leq 20$. It is clearly seen that the macrostructure of the accretion flow remains similar in all four cases, only the linear scales being different. One can see an extended heated zone ahead of the star, which is caused by Compton heating and hydrodynamic focusing effects. All variants involve a bow shock and a contact discontinuity which separates cold "external'' gas from that with the higher temperature due to the Compton heating and passage through the bow shock. The tail region possess a rather complicated shock structure, typical for jet flows. This macrostructure is intrinsic for all cases given in Table 1. The exception is Variant 15, which is given for comparison purposes. In Fig. 6, we show the Mach number and temperature logarithm distributions for this variant, which corresponds to adiabatic accretion with $V_\infty$, $T_\infty$, and $\rho_\infty$ from Variant 2 without radiation. The general structure of the flow can be seen from Figs. 7 and 8 that show the distributions of streamlines and velocity vectors for Variant 2.

Unperturbed matter starts feeling the influence of the Compton heating at distances of several Hoyle-Lyttleton radii (in the investigated range of determining parameters). Radiative cooling, however, does not allow any substantial temperature growth. Starting from a certain distance, cooling can no longer compete with heating and the temperature of the gas increases and reaches 105-106 K ahead of the bow shock. The bow shock stand-off distance from the star does not change substantially from variant to variant, being approximately equal to $(0.15\div 0.25)R_{\rm HL}$. This is comparable with the value of 0.17 $R_{\rm HL}$ given by Taam et al. (1991). Behind the bow shock, the gas heating further continues as it approaches the inner boundary. The gas temperature near it is close or even higher than the Compton temperature.


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f6.eps} \end{figure} Figure 6: Contour lines of the Mach number and the temperature logarithm for the adiabatic case. Variant 15.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f7.eps} \end{figure} Figure 7: Streamlines and velocity vectors in the box 0<x<1,-1<z<1. Variant 2.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f8.eps} \end{figure} Figure 8: Streamlines in the box 0<x<20, -20<z<20. Variant 2.
Open with DEXTER

Figure 9 shows the temperature and dimensionless mass flux distributions along the inner boundary for Variant 2. It is not surprising that the main portion of the matter is accreted from the face side - this is characteristic for accretion flows with a detached bow shock. This is also due to the local temperature maximum location in the downwind region, since it is crossed by the gas that experienced longer action of the Compton heating. The space extent of the bow shock and the size of the heated gas zone strongly depend on the amount of radiative cooling. It is largest for Variant 1, where only free-free cooling is operational. In a polytropic case or if only the Compton effects are taken into account, the space extent of the bow shock can, in principle, be infinite (see, e.g., Ruffert 1994 or Pogorelov et al. 2000). In all considered cases of radiative flows, due to intensive cooling behind the bow shock, the temperature jump across it decreases rapidly with distance to the symmetry axis, which leads to the decrease in the pressure jump across it.

  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f9.eps} \end{figure} Figure 9: Angular distributions of temperature a) and dimensionless mass flux b) through the inner boundary for Variant 2.
Open with DEXTER

As seen from Table 1, radiative effects impede accretion. All calculated accretion rates are substantially less than the Bondi-Hoyle-Lyttleton rate. What is the influence of the other factors? By comparing Variants 1 and 11 with Variants 2 and 12, we see that the choice of the cooling function only slightly affects the accretion rate. On the other hand, comparison of Figs. 2 and 3 reveals that the accretion pattern can be geometrically very different. The reason for this behavior is the fact that for all temperatures which are not too large the SB cooling is much more efficient than the BR cooling. Therefore the high temperature region will be much more extended in the bremsstrahlung-only case than in the SB case; and also the shock cone is correspondingly wider. Only in those regions which have very high temperatures (i.e. close to $T_{\rm C}$) the SB formula of Eq. (24) tends towards pure bremsstrahlung cooling. Therefore in these regions the behavior for the two cases will be similar.

Comparing variants 2, 3, and 4, which all correspond to the same Bondi-Hoyle-Lyttleton accretion rate, we see that the increase in $V_\infty$ raises the accretion rate. This can be explained by the decrease of the dimensionless parameter $S_{\rm C}$, which is responsible for the Compton heating. The cooling dimensionless parameter $S_{\rm b}$also decreases. However, the Compton effects dominate in the vicinity of the accretor in the chosen range of parameters. Accretion becomes more efficient for cooler gas. We also would like to point out in this context that if $\dot{M}_{\rm BHL}$ is fixed, $\rho_\infty$ becomes larger for larger values of $V_\infty$.

The accretion rate also decreases for higher Compton temperatures (see Table 1), since this primarily increases the Compton heating intensity. Being hotter near the inner boundary, the matter expands, and the resulting decrease of density makes the mass flux onto the accretor smaller. The role of the efficiency factor $E_{\rm eff}$ is analogous.

The last column of Table 1 shows the total luminosities of the accreting star calculated from its accretion rate. It appears that the luminosities are of the order of the Vela X-1 typical luminosity $4\times 10^{36}\ {\rm erg}~ {\rm s}^{-1}$ (Kreykenbohm et al. 1999) for all variants except 7-10, which correspond to high or low densities in the wind. Note that the Eddington luminosity for Vela X-1 is $2.4\times 10^{38}\ {\rm erg}~ {\rm s}^{-1}$, that is, radiation pressure is of minor importance.


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{3091f10.eps} \end{figure} Figure 10: Dependence of the accretion rate on the efficiency of the Compton heating. The calculated rates are shown with circles, the black dot corresponds to one of the calculations of Taam et al. (1991), and the black rectangle shows the accretion rate for adiabatic accretion with parameters of Variant 2 without radiation.
Open with DEXTER

It is interesting to examine the dependence of the dimensionless accretion rate $\dot{\cal{M}}$ on the efficiency of the Compton heating shown in Fig. 10. As a measure of this efficiency we chose the dimensionless parameter $S_{\rm C}E_{\rm eff} T_{\rm C}$, which appears as factor in the corresponding source term of the energy equation. Because of Eq. (22), this parameter is proportional to the Bondi-Hoyle-Lyttleton rate. One can see that as this efficiency coefficient grows from small values to larger ones, the accretion rate first drastically decreases, and then tends to some nearly constant value independent of the choice of the cooling function. For comparison, we included a black dot in this plot. It shows the accretion rate for Variant 3 in the paper of Taam et al. (1991). It appears that our results are in good agreement with this value. For reference purposes, we also included a black rectangle marker on the ordinate axis, which corresponds to Variant 15 that preserves the attributes of Variant 2 with omitted radiation effects.

To check the applicability of the optically-thin-medium approximation, we calculated the optical depth of plasma along different rays starting at the neutron star. The result is shown in Fig. 11 for parameters of Variant 2. It appears that the maximum of the optical depth is about 0.03 for this (base) variant (integration was performed from R=0.05 to R=20). As one would expect, the maximum optical depth is in the accretion tail. The situation is very much different in Variants 8 and 9, where the largest optical depth is observed. If we integrate up to R=12.5, the optical depth in the tail region riches the value 0.65, However, there exist a farther region around the z-axis, where the density of matter increases due to decreasing cross-section of the tail jet (see the conical contact discontinuity in Fig. 4). This is caused by the intense cooling of matter at large distances ( $z>15\ R_{\rm HL}$) behind the accretor. Although matter becomes optically thick in this region, distances to the neutron star are so large that radiative effects become of little importance for accretion.


  \begin{figure} \par\includegraphics[width=6.9cm,clip]{3091f11.eps} \end{figure} Figure 11: Optical depth of the accreting matter as a function of the line-of-sight angle for Variant 2.
Open with DEXTER

It is known that numerically obtained accretion rate may depend on the size of the accretor (Matsuda et al. 1991; Foglizzo et al. 2005). We performed calculations for the base variant (Variant 2) with sizes of the inner boundary equal to 0.05, 0.025, 0.0125, and 0.00625  $R_{\rm HL}$, while preserving the grid resolution in terms of the inner radius R*. Although it was very difficult to find any difference in distributions of quantities for all these variants, it turned out that the accretion rate gradually decreases with R*, revealing the differences between consecutive variants equal to approximately 10%, 1%, and 0.1%. Since higher resolution drastically increases the calculation time, we presented all results for $R_*=0.05\ R_{\rm HL}$.

4 Conclusions

In this paper we analyzed the influence of radiative effects on the accretion pattern in a scenario where matter is accreted from a supersonic wind possessing no angular momentum. Although this approach is oversimplified, it allows us to separate the contributions of different mechanisms and make predictions that might be applicable to widely separated binary systems, like Vela X-1. We showed that the applicability of polytropic models is very limited, and the combined influence of radiative cooling and the Compton effect (both heating and cooling) is so large that it can change the topology of the accretion pattern as well as the value of the accretion rate. Moreover, it appears that the Compton effect is a dominant factor (accretion rates calculated with the bremsstrahlung or Stellingwerf-Buff cooling mechanisms differ by about 1%). It leads to flows which always have temperatures close to the Compton temperature, $T_{\rm C}$, for the entire region around the inner boundary. Such a behavior cannot be modelled with polytropic flows. In particular, choosing the polytropic index lower than 5/3 will be in conflict with such a strong heating process. We showed that steady solutions are possible in a wide range of determining parameters and the accretion rate can be considerably smaller than the Bondi-Hoyle-Lyttleton value. This reduction factor becomes particularly large for high accretion rates. For our set of variants the actual accretion rates lie between $0.6 \times 10^{-10}$ and $1 \times 10^{-9} ~{M_{\odot}} ~ {\rm yr}^{-1}$, whereas the theoretical Bondi-Hoyle-Lyttleton rates range from 10-10 to $10^{-8} ~{M_{\odot}} ~ {\rm yr}^{-1}$. This then implies that although the observed X-ray luminosities of Vela X-1 are well within the range of our present calculations, much higher X-ray luminosities measured in some other sources will be very difficult to explain by the type of models studied here.

Acknowledgements
This work was supported in part by NSF Award ATM-0428880, NASA grant NNG05GD45G and Russian Foundation for Basic Research grant 02-01-00948. N.P. and I.K. were also supported by the Visitor Program of the Max-Planck-Institut für Astrophysik, Garching bei München. G. B. thanks the University of Tokyo, RESCEU, for its hospitality.

References

 
Copyright ESO 2005