Relativistic wind accretion on to a Schwarzschild black hole

TL;DR

This paper introduces a new relativistic wind accretion model onto a Schwarzschild black hole, extending classical theories with analytic solutions, and confirms its accuracy through numerical simulations, showing increased accretion rates at high wind speeds.

Contribution

It provides the first analytic relativistic extension of the classical wind accretion model, including explicit expressions and numerical schemes, validated against hydrodynamic simulations.

Findings

Relativistic accretion rates are significantly higher than Newtonian predictions at high wind speeds.

Analytic streamlines and accretion rates match well with numerical simulations within 10%.

The model is valid for steady, axisymmetric, ballistic flows around Schwarzschild black holes.

Abstract

We present a novel analytic model of relativistic wind accretion on to a Schwarzschild black hole. This model constitutes a general relativistic extension of the classical model of wind accretion by Bondi, Hoyle, and Lyttleton (BHL). As in BHL, this model is based on the assumptions of steady state, axisymmetry, and ballistic motion. Analytic expressions are provided for the wind streamlines while simple numerical schemes are presented for calculating the corresponding accretion rate and density field. The resulting accretion rate is greater in the relativistic model when compared to the Newtonian BHL one. Indeed, it is two times greater for asymptotic wind speeds $v_\infty \ge 0.4\,\mathrm{c}$ and more than an order of magnitude greater for $v_\infty \ge 0.8\,\mathrm{c}$. We have compared this full relativistic model versus numerical simulations performed with the free GNU Public Licensed hydrodynamics code AZTEKAS and found a good agreement for the streamlines in the upstream region of the flow and also, to within 10 per cent, for the corresponding accretion rates.

Figures

Equations

$$r$$

$$\dot{r}$$

$$\dot{\theta }$$

$$\rho$$

$$E^{\prime }=\frac{1}{2}\dot{r}(\pi {)}^2-\frac{\mathrm{G}M}{r(\pi )}=\frac{1}{2}{v}_{\infty }^2-\frac{2(\mathrm{G}M{)}^2}{b^2{v}_{\infty }^2}.$$

$$b_c=\frac{2\mathrm{G}M}{v_{\infty }^2},$$

$${\dot{M}}_{\mathrm{BHL}}=\pi b_c^2{\rho }_{\infty }{v}_{\infty }=4\pi {\rho }_{\infty }\frac{(\mathrm{G}M{)}^2}{v_{\infty }^3}.$$

$$\frac{\mathrm{d}t}{\mathrm{d}\tau }=\mathcal{E}{\left(1-\frac{2M}{r}\right)}^{-1},$$

$$\frac{\mathrm{d}r}{\mathrm{d}\tau }=-{\left[{\mathcal{E}}^2-\left(1-\frac{2M}{r}\right)\left(1+\frac{h^2}{r^2}\right)\right]}^{1/2},$$

$$\frac{\mathrm{d}\theta }{\mathrm{d}\tau }=\frac{h}{r^2},$$

$$\mathcal{E}=\frac{1}{\sqrt{1-v_{\infty }^2}}={\gamma }_{\infty },$$

$$h=\frac{b{v}_{\infty }}{\sqrt{1-v_{\infty }^2}}=b{V}_{\infty },$$

$$\frac{\mathrm{d}r}{\mathrm{d}\theta }=\frac{\sqrt{\mathcal{R}(r)}}{h},$$

$$\mathcal{R}(r)=r\left[V_{\infty }^2{r}^3+2M{r}^2-b^2{V}_{\infty }^2(r-2M)\right].$$

$$b_0=\frac{M\sqrt{27V_{\infty }^4+18V_{\infty }^2+{\gamma }_{\infty }{\left(1+9V_{\infty }^2\right)}^{3/2}-1}}{\sqrt{2}{V}_{\infty }^2}.$$

$$r=\frac{r_1(r_3-r_2)-r_2(r_3-r_1){\mathrm{cn}}^2\left(\xi ,k\right)}{r_3-r_2-(r_3-r_1){\mathrm{cn}}^2\left(\xi ,k\right)},$$

$$\xi =\frac{\sqrt{r_3(r_3-r_2)}}{2b}(\theta +{\theta }_{\infty }),$$

$$k=\sqrt{\frac{r_2(r_3-r_1)}{r_3(r_2-r_1)}},$$

$${\theta }_{\infty }=\frac{2b}{\sqrt{r_3(r_3-r_2)}}{\mathrm{cn}}^{-1}\sqrt{\frac{r_3-r_2}{r_3-r_1}},$$

$$r=\frac{\delta r_1[1-\mathrm{cn}(\tilde{\xi },\tilde{k})]}{\delta -\eta -(\eta +\delta )\mathrm{cn}(\tilde{\xi },\tilde{k})},$$

$$\tilde{\xi }=\frac{\sqrt{\eta \delta }}{b}(\theta +{\tilde{\theta }}_{\infty }),$$

$$\eta =\sqrt{(r_2-r_1)(r_3-r_1)},$$

$$\delta =\sqrt{r_2{r}_3},$$

$${\tilde{k}}^2=\frac{(\eta +\delta {)}^2-r_1^2}{4\eta \delta },$$

$${\tilde{\theta }}_{\infty }=\frac{b}{\sqrt{\eta \delta }}{\mathrm{cn}}^{-1}\left[\frac{\delta -\eta }{\eta +\delta },\tilde{k}\right].$$

$${\nabla }_{\mu }(\rho U^{\mu })=0,$$

$$\frac{\partial }{\partial x^i}(r^2\sin \theta \rho U^i)=0.$$

$$\rho r{U}^{\theta }{S}_{\theta }={\rho }_{\infty }{\gamma }_{\infty }{v}_{\infty }{S}_{\infty },$$

$$\rho r^2{U}^{\theta }\sin \theta \mathrm{d}r={\rho }_{\infty }{\gamma }_{\infty }{v}_{\infty }b\mathrm{d}b.$$

$$\rho =\frac{{\rho }_{\infty }}{\sin \theta }{\left(\frac{\partial r}{\partial b}\right)}^{-1}.$$

$$r(\pi {)}^3+[2M-r(\pi )]b_c^2=0.$$

$$b_c^{\mathrm{fit}}=\frac{2M}{v_{\infty }^2}\left(1+c_1{v}_{\infty }+c_2{v}_{\infty }^2+c_3{v}_{\infty }^3\right),$$

$$\dot{M}=\pi b_c^2{\rho }_{\infty }{v}_{\infty }{\gamma }_{\infty }.$$

$${\nabla }_{\mu }{T}^{\mu \nu }=0,$$

$$T^{\mu \nu }=\rho h{U}^{\mu }{U}^{\nu }+p{g}^{\mu \nu },$$

Full text

††pagerange: Relativistic wind accretion onto a Schwarzschild black hole–A.2

1 Introduction

Accretion physics has become a basic tenet of astrophysics ever since it was recognized that the process of accretion, especially when compact objects are involved, is one of the most efficient mechanisms for converting rest mass energy into luminosity at work in our universe (Frank et al., 2002). Indeed, it has been established that a thin accretion disc around a non-rotating black hole can reprocess as much as 10 per cent of the accreted gas rest mass into electromagnetic radiation, while for a maximally rotating black hole this figure can reach up to 46 per cent (see e.g. Longair, 2011).

The simplest accretion scenario consists of the stationary, spherically-symmetric solution first discussed by Bondi (1952), where he considered an infinitely large homogeneous gas cloud steadily accreting onto a central gravitational object. The general relativistic extension of this model was developed by Michel (1972) who took a Schwarzschild black hole as the central accretor.

In the so-called wind accretion scenario, the spherical symmetry approximation is relaxed by considering a non-zero relative velocity between the central object and the accreted medium (cf. Edgar, 2004; Romero & Vila, 2014). As it turns out, even after assuming a steady state and axial symmetry, the problem becomes sufficiently complex as not to admit a full analytic solution in general. In the pioneering work of Hoyle & Lyttleton (1939) and Bondi & Hoyle (1944) (BHL hereafter), the authors provided an analytic model for supersonic wind accretion by adopting the so-called ballistic approximation, i.e. by neglecting pressure gradients and assuming that the fluid’s dynamics is solely dictated by the central object’s gravitational field. This approximation is well suited to describe highly supersonic flows given that, within this regime, incoming fluid elements cannot oppose pressure gradients readily, effectively following nearly free-fall trajectories.

Further analytic solutions to wind accretion problems have been found for a perfect fluid with a stiff equation of state in general relativity (Petrich et al., 1988) and for the corresponding non-relativistic case of an incompressible fluid (Tejeda, 2018).

The problem of wind accretion has also been the focus of various numerical studies, from both a Newtonian perspective (Hunt, 1971; Shima et al., 1985; Ruffert & Arnett, 1994; El Mellah & Casse, 2015; El Mellah et al., 2018) and in general relativity (Petrich et al., 1989; Font & Ibáñez, 1998; Zanotti et al., 2011; Lora-Clavijo & Guzmán, 2013; Gracia-Linares & Guzmán, 2015; Cruz-Osorio & Lora-Clavijo, 2016; Cruz-Osorio et al., 2017).

In this article we introduce a simple, analytic model for a supersonic wind accreting onto a non-rotating black hole (Schwarzschild spacetime). This model is a general relativistic extension of the BHL model and follows closely the methodology outlined in Mendoza et al. (2009), Tejeda et al. (2012, 2013). In that series of works, we presented an analytic model of the accretion flow of a rotating dust cloud infalling towards a central object, first in a Newtonian regime (Mendoza et al., 2009) and then in general relativity for a Schwarzschild black hole (Tejeda et al., 2012) and for a Kerr black hole (Tejeda et al., 2013). See Schroven et al. (2017) for a recent extension of this model considering the case of charged dust particles accreting onto a Kerr-Newman black hole.

The article is organised as follows. In Section 2 we give a brief review of the BHL model. In Section 3 we introduce the new relativistic wind model, and in Section 4 we present the comparison against relativistic hydrodynamic numerical simulations performed with the free GPL hydrodynamics code aztekas (Olvera & Mendoza, 2008; Aguayo-Ortiz et al., 2018).111aztekas.org ©2008 Sergio Mendoza & Daniel Olvera and ©2018 Alejandro Aguayo-Ortiz & Sergio Mendoza. Finally, in Section 5 we summarize the work and present our conclusions. The Appendix A presents the results of the benchmark test of aztekas against Michel’s analytic model, as well as self-convergence and various resolution tests relevant to the present study.

2 Bondi-Hoyle-Lyttleton model

In this Section we give a brief overview of the BHL accretion model as its general relativistic extension constitutes the scope of the present work. The BHL model deals with a steady, supersonic wind accreting onto a gravitational object of mass $M$. Considering that the accretor is held fixed at the centre of coordinates, infinitely far away from the central object, the wind is homogeneous and characterized by a uniform density $\rho_{\infty}$ and wind speed $v_{\infty}$.

Under the ballistic approximation, expressions for the streamlines, velocity field, and density for this model are given by (Bisnovatyi-Kogan et al., 1979)

[TABLE]

where $b$ is the impact parameter of a given incoming fluid element. Note that we are taking here a reference frame in spherical coordinates such that the polar axis is aligned with the incoming wind direction. The wind comes asymptotically from the direction of the $\theta=0$ axis. Figure 1 shows a schematic representation of the model setup and streamlines.

Given the uniform wind condition at infinity, all of the incoming trajectories have initially a common specific total energy $E=v^{2}_{\infty}/2$. In other words, the wind fluid elements follow energetically unbound (hyperbolic) trajectories. From Eqs. (2.1) and (2.2) we see that the incoming trajectories reach the downstream axis at $r(\pi)=b^{2}v^{2}_{\infty}/(2\mathrm{G}M)$ with radial velocity $\dot{r}=v_{\infty}$. It is expected that streamlines coming from mirror-reflected points with respect to the symmetry axis will collide with one another along this axis (see Figure 1). Following this collision, the BHL model envisages that the component of the velocity perpendicular to the axis is instantly transformed into thermal energy and, eventually, lost as radiation.222More precisely, in the Hoyle & Lyttleton (1939) model the fluid streamlines are assumed to focus downstream onto the symmetry axis leading to an infinite density there. The accretion then proceeds along this line. Instead, Bondi & Hoyle (1944) envisioned an accretion flow that spreads out onto a finite density region that they referred to as accretion column. In this work we are adopting the former description. After this loss of kinetic energy, each fluid element is left with a new specific total energy

[TABLE]

By equating this energy to zero, the following critical value for the impact parameter is found

[TABLE]

such that any fluid element following a streamline with $b<b_{c}$ is energetically bound to the central object and, hence, eventually accreted. On the other hand, any fluid element following a streamline with $b>b_{c}$ is energetically unbound to the central object and, therefore, ultimately escapes to infinity. Following this argument, and accounting for all of the material within the cylinder of radius $b_{c}$, the total accretion rate onto the central object in the BHL model is given by

[TABLE]

Subsequent numerical studies of this problem have found stationary accretion rates that agree remarkably well with those predicted by the simple BHL model (see e.g. Hunt, 1971). The resulting flows obtained in full-hydrodynamic simulations of supersonic wind accretion show the development a bow shock around the central accretor at which the incoming wind streamlines abruptly decelerate and become subsonic. Clearly, the simple analytic description of the streamlines provided by the BHL model is no longer valid inside the bow shock, nevertheless, it provides a qualitatively good description of the streamlines in the supersonic, upwind region.

3 Relativistic wind model

In this section we present the extension of the BHL wind model to the case in which the central accretor is a non-rotating black hole of mass $M$. We shall assume that the mass-energy content of the accreting gas is negligible as compared to the mass of the central black hole and, thus, that the overall spacetime metric corresponds to the Schwarzschild solution. For constructing this model we are closely following the methodology described in Tejeda et al. (2012). For the remainder of this work we adopt a geometrised system of units for which $\mathrm{G}=\mathrm{c}=1$.

3.1 Velocity field

Adopting the ballistic approximation for test particles in general relativity amounts to describe the incoming streamlines as geodesic trajectories. Thanks to the symmetries of Schwarzschild spacetime, the trajectory of a test particle is restricted to a plane and governed by the equations of motion (Frolov & Novikov, 1998)

[TABLE]

where $\mathcal{E}$ is the conserved specific (relativistic) total energy and $h$ is the conserved specific angular momentum. Assuming a uniform wind velocity at infinity $v_{\infty}$, these conserved quantities for a given streamline with impact parameter $b$ are given by

[TABLE]

where $\gamma_{\infty}$ is the wind’s Lorentz factor at infinity and where we have introduced the shorthand notation $V_{\infty}=\gamma_{\infty}\,v_{\infty}$.

In order for the velocity field in Eqs. (3.1)-(3.3) to be useful in practice, we need to provide an expression for the streamlines of the form $r(\theta,b)$ as we discuss in the next subsection.

3.2 Streamlines

An expression for the streamlines can be obtained by combining Eqs. (3.2) and (3.3) as

[TABLE]

where

[TABLE]

As discussed in detail in Tejeda et al. (2012), Eq. (3.6) can be solved in terms of elliptic integrals. For the problem at hand, we need to distinguish between two types of trajectories: 1) unbound trajectories that reach a minimum distance in their descent towards the central object before turning back to infinity and 2) trapped trajectories that plunge onto the black hole’s event horizon (located at $r=2M$). Specifically, streamlines with $b\geq b_{0}$ belong to type 1) while those with $b<b_{0}$ belong to type 2) where

[TABLE]

In case 1), the polynomial in Eq. (3.7) has three non-trivial real roots. For the particular boundary condition that we have adopted here, it can be proved that one of these roots is negative and the other two are positive. We call them $r_{1}$, $r_{2}$ and $r_{3}$, such that $r_{1}<0<r_{2}<r_{3}$. These roots can be explicitly given in terms of the constants of motion, see e.g. Tejeda et al. (2012). In terms of these roots, the equation for the streamlines is given by

[TABLE]

where $\mathrm{cn}\left(\xi,\,k\right)$ is a Jacobi elliptic function with modulus (Lawden, 1989)

[TABLE]

and

[TABLE]

is the polar phase setting $\theta=0$ as the incoming wind direction, i.e. $r(\theta=0)\rightarrow\infty$.

In case 2) the polynomial in Eq. (3.7) has as roots a negative real number $r_{1}$ and a complex conjugate pair $r_{2}=r_{3}^{*}$. In this case the equation for the streamlines can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

Equations (3.9) and (3.13) constitute the required expressions for the streamlines. As demonstrated by Tejeda et al. (2012), these expressions for the streamlines reduce to the usual Newtonian conic-sections in the non-relativistic limit, i.e. for $v_{\infty}\ll 1$ Eq. (3.9) reduces to Eq. (2.1), while the velocity field in Eqs. (3.2) and (3.3) reduce to Eqs. (2.2) and (2.3), respectively.

3.3 Density field

For calculating the density field we follow here a similar strategy as in Tejeda et al. (2012, 2013). We start from the continuity equation333Here and in what follows Greek indices run over spacetime components, Latin indices run over spatial components only, and Einstein’s summation convention over repeated indices is adopted.

[TABLE]

where $U^{\mu}=\mathrm{d}x^{\mu}/\mathrm{d}\tau$ is the four-velocity, $\rho$ is the rest mass density and $\nabla_{\mu}$ stands for the covariant derivative.

Writing Eq. (3.19) for a Schwarzschild spacetime and adopting the stationary condition results in

[TABLE]

Let us now integrate Eq. (3.20) over the spatial volume element delimited by a sufficiently small set of streamlines that start from an area element $S_{\infty}=2\pi\,b\,\mathrm{d}b$ located at infinity and that end at the intersection with the conic surface defined by $\theta=\mathrm{const.}$, i.e.  $S_{\theta}=2\pi\,r\,\sin\theta\,\mathrm{d}r$, as shown in Figure 2. By construction, fluid elements enter the integration volume only across $S_{\infty}$ and leave across $S_{\theta}$. Therefore, by means of Gauss’s theorem, it follows that

[TABLE]

from where

[TABLE]

Finally, substituting $U^{\theta}$ from Eq. (3.3) into Eq. (3.22) results in

[TABLE]

Calculating explicitly the partial derivative in Eq. (3.23) is something trivial to do in the Newtonian case. Indeed, this calculation is the step leading from Eq. (2.1) to Eq. (2.4). However, this same calculation in the relativistic case represents a complex procedure involving the derivative of an elliptic function with respect to its argument and modulus. Following Tejeda et al. (2012), we do not attempt here to calculate this derivative analytically but rather use a finite difference scheme to compute it numerically. A suitable grid for performing this calculation can be constructed as follows: Start from a collection of streamlines separated by uniform intervals of $\Delta b$ at infinity. Follow these streamlines from $\theta=0$ to $\theta=\pi$ storing the different values of $r$ at uniform steps of $\Delta\theta$. Use these grid values for estimating $\partial r/\partial b$ as a finite difference of the radial coordinate between neighboring streamlines. Such a grid is schematically represented in Figure 2.

In Figure 3 we show an example of a wind accreting at $v_{\infty}=0.5$ onto a Schwarzschild black hole. The figure shows the flow streamlines as expressed analytically by Eqs. (3.9) and (3.13) together with the density field as calculated numerically from Eq. (3.23).

In Figure 4 we compare the resulting streamlines in Schwarzschild spacetime with those coming from the Newtonian BHL model for different values of the asymptotic wind speed. From this figure it is clear that the accretion cylinder, and thus the total accretion rate, in the relativistic model is greater than the corresponding Newtonian value. In the next subsection we discuss this in further detail.

3.4 Accretion rate

Just as for the density field, the procedure described in Section 2 to compute the total accretion rate in the Newtonian case (BHL model) is not as simple to implement for the relativistic problem. However, the logic behind this calculation remains the same: when the collision of mirror-reflected streamlines along the symmetry axis $(\theta=\pi)$ occurs, the flow’s kinetic energy associated to the normal component of the velocity ($U^{\theta}$ in this case) is lost from the system as thermal energy and/or radiation. The remaining relativistic energy $\mathcal{E}^{\prime}$ of a given streamline with impact parameter $b$ will determine whether the gas traveling along the streamline is accreted ($\mathcal{E}^{\prime}<1$) onto the central black hole or lost to infinity ($\mathcal{E}^{\prime}>1$). The streamline that is left marginally bound with $\mathcal{E}^{\prime}=1$ is characterized by the critical impact parameter $b_{c}$. Using the normalization condition $U_{\mu}U^{\mu}=-1$, together with $U^{\theta}=0$ and Eq. (3.2), it is simple to show that the condition $\mathcal{E}^{\prime}=1$ is equivalent to

[TABLE]

Unfortunately, after substituting $r(\pi)$ from either Eqs. (3.9) or (3.13), the resulting equation is highly non-linear in $b_{c}$ and it does not seem possible to solve it for $b_{c}$ explicitly. Nonetheless, it is straightforward to solve Eq. (3.24) numerically using a root-finding algorithm. We have done this using the bisection method for 1000 values of $v_{\infty}$ uniformly distributed between 0.001 and 0.999 to a precision of 10*-8*. Now, based on these numerically calculated values, we found the following fit that approximates $b_{c}$ to an accuracy better than $3\%$ for $v_{\infty}<0.999$ and to within $0.5\%$ for $v_{\infty}<0.98$:

[TABLE]

with $c_{1}=0.081135$, $c_{2}=3.452826$ and $c_{3}=-1.758438$.

All of the streamlines within the cylinder $b\leq b_{c}$ contribute to the total accretion rate. According to the right-hand side of Eq. (3.22), we can therefore express $\dot{M}$ as

[TABLE]

In Figure 5 we show the resulting accretion rate as a function of the wind speed at infinity and compare it to the Newtonian value found in the BHL model (Eq. 2.7). As can be seen from this figure, the relativistic accretion rate is greater than the corresponding BHL one for all values of $v_{\infty}$, with the difference being more than double for $v_{\infty}\gtrsim 0.4$ and more than ten times larger for $v_{\infty}\gtrsim 0.8$. As expected, the relativistic accretion rate converges to the Newtonian value in the non-relativistic limit, i.e for $v_{\infty}\ll 1$. Contrary to the Newtonian case in which the accretion rate was a monotonically decreasing function of $v_{\infty}$, the relativistic result has an inflection point at around $v_{\infty}\simeq 0.8$ and becomes arbitrarily large as $v_{\infty}\rightarrow 1$. The culprit behind this behaviour is the special relativistic effect of Lorentz contraction that compresses the fluid volume elements along the direction of the wind and is represented by the Lorentz factor in Eq. (3.26).

4 Comparison with numerical simulations

In this section we compare the analytic wind accretion model presented in the previous section against numerical hydrodynamic simulations performed with aztekas, a free GNU Public Licensed (GPL) code for solving any conservative set of equations, in particular, relativistic (with a non-trivial fixed metric) and non-relativistic hydrodynamic equations. Some aspects of aztekas, together with numerical and convergence tests relevant to the present work are discussed in the Appendix A. Further details and tests of the aztekas code will be presented elsewhere.

For all of the simulations discussed in this work, we considered a perfect fluid evolving on top a fixed Schwarzschild background metric described in terms of the horizon-penetrating, Kerr-Schild coordinates. Since the problem under study presents axial symmetry, all of the simulations were performed in 2D using polar coordinates $r$ and $\theta$. All of the results presented below were obtained after the numerical simulations had evolved in time from a uniform initial state condition until a relaxed stationary state was reached.

The relativistic hydrodynamic equations consist on the continuity equation (3.19) together with the local conservation of energy-momentum

[TABLE]

where we take the energy-momentum tensor $T^{\mu\nu}$ corresponding to a perfect fluid (Landau & Lifshitz, 1975)

[TABLE]

with rest mass density $\rho$, pressure $p$, specific internal energy $\epsilon$ and specific enthalpy $h=1+\epsilon+p/\rho$. Moreover, we adopt a Bondi-Wheeler equation of state (Tooper, 1965) of the form

[TABLE]

where $\Gamma$ is the polytropic index.

In order to integrate numerically Eqs. (3.19) and (4.1) with aztekas, we recast them in a conservative form based on the 3+1 formalism (see e.g. Font, 2000; Alcubierre, 2008), as follows

[TABLE]

where $\gamma$ and $g$ are the determinants of the spatial 3-metric $\gamma_{ij}$ and the spacetime 4-metric $g_{\mu\nu}$, respectively. $\mathbf{Q}$ is the conservative variable vector, $\mathbf{F}^{i}=\left\{\mathbf{F}^{r},\mathbf{F}^{\theta}\right\}$ are the fluxes along the $r$ and $\theta$ coordinates and $\mathbf{S}$ is the source vector. All of these quantities depend on the primitive variables $\mathbf{U}$. The functional form of these vectors are presented below:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\Gamma_{\mu\nu}^{\lambda}$ are the usual Christoffel symbols, $v^{i}$ is the three-velocity as measured by local Eulerian observers, $W=\left(1-\gamma_{ij}v^{i}v^{j}\right)^{-1/2}$ is the associated Lorentz factor, and

[TABLE]

[TABLE]

are the lapse function and the shift vector, respectively, and correspond to the 3+1 decomposition of Schwarzschild spacetime in Kerr-Schild coordinates.

For the numerical flux calculation, we use in aztekas a high resolution shock capturing method (HRSC) with an HLLE approximate Riemann solver (Harten et al., 1983), combined with a monotonically centered (MC) second order reconstructor at cell interfaces. For the time integration, we use a second order Runge-Kutta method of lines in the total variation diminishing (TVD) version (Shu & Osher, 1988). Finally, we adopt a constant time-step defined through the CFL condition $\Delta t=C\,\min\left(\Delta r,r\,\Delta\theta\right)$, with $C=0.1$.

For the relativistic wind simulations, we ran two sets of simulations exploring asymptotic wind speeds from $v_{\infty}=0.1$ to $0.9$ for two different polytropic indices: $\Gamma=4/3$ and $5/3$. We took as numerical domain $[r_{\mathrm{min}},\ r_{\mathrm{max}}]\times[0,\pi]$, with $r_{\mathrm{min}}=0.5\,r_{\mathrm{acc}}$, $r_{\mathrm{max}}=10\,r_{\mathrm{acc}}$ and

[TABLE]

where $a_{\infty}$ is the asymptotic speed of sound. The radius $r_{\mathrm{acc}}$ is commonly used in the literature (cf. Font et al., 1999; Cruz-Osorio et al., 2012), as a good estimate for the extension of the numerical domain necessary for numerical convergence.444Note that with our choice of $r_{\mathrm{min}}$, simulations with $v_{\infty}\geq 0.5$ are such that $r_{\mathrm{min}}<2M$, which aztekas can handle without any problem thanks to our choice of horizon penetrating coordinates. On the other hand, for those simulations with $r_{\mathrm{min}}>2M$, we performed trial tests to make sure that the resulting steady-state accretion rate was unchanged independently of whether the black hole event horizon ($2M$) was part of the numerical domain or not. Our choice of this numerical domain is based on the convergence and resolution tests discussed in detail in the Appendix A.

For most of these simulations we used a fixed Mach number $\mathcal{M}=v_{\infty}/a_{\infty}=5$ in order to ensure that we were considering the same supersonic conditions in all cases. As a sanity check, we also considered $\mathcal{M}=10$ for a reduced number of cases. We used a uniform grid of $400\times 400$ cells except for the cases $v_{\infty}=0.1$, $0.8$ and $0.9$, where we had to use a larger grid of $1000\times 1000$ in order to find converged solutions.

As external boundary condition at the sphere $r=r_{\mathrm{max}}$, we enforced a constant, uniform inflow from the northern hemisphere $\theta\in[0,\pi/2)$ and free outflow from the southern one $\theta\in[\pi/2,\pi]$. As internal boundary condition, we set free outflow across the sphere at $r=r_{\mathrm{min}}$. The incoming wind at the external boundary has uniform thermodynamical variables $\rho_{\infty}$ and $p_{\infty}$. We set $\rho_{\infty}=10^{-10}$ in arbitrary units and take $p_{\infty}$ consistent with the chosen Mach number $\mathcal{M}$ and the equation of state, as in Cruz-Osorio et al. (2012):

[TABLE]

For the velocity field, we set an incoming wind with a constant velocity $\sqrt{v_{i}\,v^{i}}=v_{\infty}$, and components:

[TABLE]

[TABLE]

where $g_{rr}$ and $g_{\theta\theta}$ are the radial and polar components of the metric respectively. As for the initial conditions of the simulation, we set them equal to the constant boundary values over all the numerical domain.

The simulations were left to run until a stationary state was reached. This was monitored by keeping track of the mass accretion rate, which was computed on the fly by integrating the relativistic radial mass flux across a control sphere of radius $r$ according to (Petrich et al., 1989)

[TABLE]

It is important to remark that, once the stationary state is reached across all the numerical domain, the mass accretion rate as calculated from Eq. (4.15) has to have a constant value (to within numerical precision) independently of where the control radius $r$ is located (see Figure 12).

The typical simulation time at which steady state was reached depends on the wind velocity at infinity and the domain extension roughly as $t_{c}\approx r_{\mathrm{acc}}/v_{\infty}$. As an example, in Figure 6 we show the mass accretion rate as a function of time for $v_{\infty}=0.5$ and $\Gamma=4/3,\ 5/3$, measuring it at the event horizon $r=2M$.

In Figure 7, we show the steady-state of the density field and streamlines of a wind accretion flow with $v_{\infty}=0.5$ and polytropic index $\Gamma=4/3$ on the left side panel and $\Gamma=5/3$ on the right side panel. As expected for a supersonic flow, a bow shock is formed downstream around the accretor, with a smaller shock cone for $\Gamma=4/3$ than for $\Gamma=5/3$.

In Figure 8, we compare the streamlines as obtained from the two simulations discussed in the previous paragraph against the ones from the analytic model. As can be seen from this figure, the ballistic approximation provides a qualitatively good description of the resulting streamlines upstream of the flow and in the region outside the bow shock, while inside the shock cone the streamline behavior is notably different.

We can also notice from both Figure 6 and Figure 8, that there is a better agreement between the analytic model and the $\Gamma=4/3$ case than for the $\Gamma=5/3$ one. This is related to the fact that the relative degree of incompressibility or stiffness of a polytropic fluid is directly proportional to the adiabatic index and, thus, a $\Gamma=5/3$ fluid resists more effectively the compression due to the gravitational field of the central object (geodesic focusing) than a fluid with $\Gamma=4/3$.

In order to compare the numerical mass accretion rate with the ballistic model (Eq. 3.24), we computed the mean value of Eq. (4.15) across all the radial domain. In Figure 9 we show the comparison between the analytic model presented in the previous section and the simulations for both $\Gamma=4/3$ and $5/3$.

Moreover, we also explored the dependence of the mass accretion rate on the Mach number by performing 10 additional simulations with $\mathcal{M}=10$. In Figure 9, the square marks show the mass accretion rate for these simulations. Note that these results overlay with the $\mathcal{M}=5$ cases (the relative difference between both results is less than $0.7\%$). The resulting bow shocks are also very similar in both cases, with slightly smaller cones for the $\mathcal{M}=10$ case.

As can be seen in Figure 9, the resulting accretion rate from the numerical simulations is consistently slightly larger for the $\Gamma=4/3$ simulations than for the $\Gamma=5/3$ ones. Moreover, a good agreement is found between the accretion rate as predicted by the relativistic model and the one obtained from the aztekas simulations, with a relative error less than 10% for all of the 28 simulations presented in this work.

5 Summary

We have presented a full relativistic, analytic model of a supersonic wind accreting onto a Schwarzschild black hole. In addition to the assumptions of stationarity and axisymmetry, the model is based on the ballistic approximation in which the streamlines of the accretion flow correspond to geodesic trajectories of a Schwarzschild spacetime. Following the methodology presented in Tejeda et al. (2012), the streamlines of the model were described analytically in terms of Jacobi elliptic functions. The density field of the resulting accretion flow and the corresponding accretion rate were calculated using simple numerical schemes.

The model presented in this paper constitutes the relativistic generalisation of the Newtonian wind accretion model by Bondi-Hoyle-Lyttleton (BHL). Naturally, the relativistic model recovers the BHL model in the non-relativistic limit $v_{\infty}\ll\mathrm{c}$. The enhanced gravitational field of the accreting object in general relativity, together with the special relativistic effect of Lorentz contraction, contribute to a larger accretion rate as compared to that of the BHL model (see Figure 5). This difference becomes substantial (by a factor of 10-100) for asymptotic wind speeds $v_{\infty}$ close to the speed of light c. Although these large velocities are not expected to be common in astrophysical settings, they can appear in extreme cases such as the velocity kick imparted onto a newborn black hole following an asymmetrical supernova explosion (Janka, 2013) or after the merger of two rotating black holes (Gerosa & Moore, 2016).

We have compared the new relativistic model against numerical simulations performed with the aztekas code. This code solves numerically the full hydrodynamic evolution of a perfect fluid in Schwarzschild spacetime starting off from a uniform condition until a stationary state is reached. We have used two different polytropic equations of state ($\Gamma=4/3$ and $\Gamma=5/3$) and asymptotic wind speeds from $v_{\infty}=0.1\,\mathrm{c}$ to $v_{\infty}=0.9\,\mathrm{c}$. We have considered two different values for the asymptotic Mach number, $\mathcal{M}=5$ and 10, finding virtually no difference between these two values. We have found a good agreement (to within 10%) between the accretion rate predicted by the relativistic model and the aztekas simulations (see Figure 9). As expected for these supersonic flows, the ballistic streamlines of the analytic model agree quite well with the resulting streamlines of the numerical simulations in the upwind region outside the bow shock.

6 Acknowledgements

We thank Sergio Mendoza, Alejandro Cruz-Osorio, Olivier Sarbach and Francisco S. Guzmán for useful discussions and comments on the manuscript. We also thank the anonymous referee for helpful suggestions and remarks. This work was supported by DGAPA-UNAM (IN112616 and IN112019) and CONACyT (CB-2014-01 No. 240512; No. 290941; No. 291113) grants.

Appendix A Validation of the numerical hydrodynamic code aztekas

In this Appendix we present several numerical tests intended to validate our use of aztekas in this work.

A.1 Spherical accretion

In order to test aztekas in the general relativistic hydrodynamic regime, it is important to compare with a benchmark solution. For this, we employ the analytic model of spherical accretion developed by Michel (1972).

In Figure 10, we compare the outcome of numerical simulations performed with aztekas against Michel’s analytic solution for the same values of polytropic index used in the wind accretion problem ($\Gamma=4/3,5/3$) and for an asymptotic sound speed of $a_{\infty}=0.01$.

For these simulations, we took a 2D spherical axisymmetric grid of $400\times 400$ uniformly distributed radial bins $r\in\left[0.5M,20M\right]$ and polar bins $\theta\in\left[0,\pi/2\right]$. We set the boundary conditions at $r=20M$ by imposing Michel’s solution there. As initial conditions we populate the entire numerical domain with the boundary constant value and let the system evolve until a stationary regime is reached for $t\gtrsim 100M$. As can be seen from Figure 10, an excellent agreement is found between Michel’s analytic solution and the aztekas simulation results. We compute the mass accretion rate and the relative error between both solutions is below 0.2%.

Also using this benchmark solution, we looked at the convergence rate of aztekas by computing the $L^{1}$ norm of the density error for different numerical resolutions $N_{r}$,555Since no angular dependence is found for this spherically symmetric test, we only varied the radial resolution while keeping a constant $N_{\theta}=400$. i.e.

[TABLE]

where $\rho_{\mathrm{num}}$ and $\rho_{\mathrm{exact}}$ are the numerical and exact values of the density, respectively. In Figure 11, we show the result of this test from where we obtain a convergence rate consistent with a second order, which is to be expected for smooth solutions and for our use of the MC reconstructor.

A.2 Dependence on domain extension and resolution

We also want to make sure that the numerical solutions obtained with aztekas have converged to physical values and that the results are independent from both numerical resolution and domain extension.

Due to the finite extension of the numerical domain $r_{\mathrm{max}}$, the total accretion rate calculated from the numerical simulations shows a small dependence on $r_{\mathrm{max}}$ that gets weaker as larger domain extensions are considered. For example, in Figure 12 we show the resulting mass accretion rate for simulations with $v_{\infty}=0.5$, $\Gamma=5/3$ and seven different domain extensions, at $t=1000M$. As can be seen from this figure, our choice of a domain extension of $r_{\mathrm{max}}=10\,r_{\mathrm{acc}}$ leads to an overestimation of the mass accretion rate of the order of $\sim 5\%$, which is an acceptable margin of error for us in this work. We expect the rest of the simulations to behave in a similar way.

Similarly, we also looked at this same domain extension dependence but now by monitoring the mass accretion rate across a fixed radius (in this case across the event horizon at $r=2M$) as a function of time. As can be seen from Figure 13, the difference between a domain extension of $10\,r_{\mathrm{acc}}$ and of $25\,r_{\mathrm{acc}}$ maintains the same margin error of $\sim 5\%$ along time.

In what regards the numerical resolution, we found that the grid size does not affect as directly the value of the resulting accretion rate. However, it does contribute to the smoothness of the solution, with larger resolutions leading to less numerical noise. In Figure 14 we show an example of this again for the case $v_{\infty}=0.5$, $\Gamma=5/3$ and five different numerical resolutions: $\mathcal{R}_{1}=200\times 200$, $\mathcal{R}_{2}=300\times 300$, $\mathcal{R}_{3}=400\times 400$, $\mathcal{R}_{4}=500\times 500$ and $\mathcal{R}_{5}=600\times 600$. As discussed in the Section 4, for most of the simulations presented in this work we settled for $\mathcal{R}_{3}$ as a good compromise between accuracy and performance.

Finally, in order to estimate the convergence rate of the simulations, and given that there is no analytical solution for the full hydrodynamic wind accretion, we follow Lora-Clavijo & Guzmán (2013) and measure the self-convergence rate using different resolutions. To this purpose, we ran two sets of three simulations each, with resolutions

[TABLE]

and

[TABLE]

Note that, for each set separately, there is a factor of $k=1.5$ and $k=2$ between successive resolutions. With the simulation results at hand, we compute the self-convergence factor $Q$ according to

[TABLE]

with $L^{1}(\rho_{i},\rho_{j})$ as defined in Eq. (A.1) and where $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$ are the densities along the accretion axis $(\theta=\pi)$ for each of the employed resolutions. The results for this test are shown in Figure 15. As can be seen from this figure, for both set of resolutions we find an order of convergence between 1.5 and 2 in the steady state region, which is to be expected for this kind of algorithms due to the presence of shocks.