# Numerical simulations of relativistic wind accretion on to black holes using Godunov-type methods

## 1 Introduction

The term “wind” or hydrodynamic accretion refers to the capture of matter by a moving object under the effect of the underlying gravitational field. The canonical astrophysical scenario in which matter is accreted in such a non-spherical way was suggested originally by Bondi and Hoyle , who studied, using Newtonian gravity, the accretion on to a gravitating point mass moving with constant velocity through a non-relativistic gas of uniform density. Such process applies to describe mass transfer and accretion in compact X-ray binaries, in particular in the case in which the donor (giant) star lies inside its Roche lobe and loses mass via a stellar wind. This wind impacts on the orbiting compact star forming a bow-shaped shock front around it.

The problem was first numerically investigated in the early 70’s. Since then, contributions of a large number of authors using highly developed Godunov-type methods extended the simplified analytic models (see, e.g.,  and references there in). These Newtonian investigations helped develop a thorough understanding of the hydrodynamic accretion scenario, in its fully three-dimensional character, revealing the formation of accretion disks and the appearance of non-trivial phenomena such as shock waves or flip-flop instabilities.

We have recently considered hydrodynamic accretion on to a moving black hole using relativistic gravity and the “test fluid” approximation . We present here a brief summary of the methodology and results of such simulations. We integrate the general relativistic hydrodynamic equations in the fixed background of the Kerr spacetime (including its non-rotating Schwarzschild limit) and neglect the self-gravity of the fluid as well as non-adiabatic processes such as viscosity or radiative transfer. In the black hole case the matter flows ultimately across the event horizon and becomes causally disconnected of distant observers . Near that region the problem is intrinsically relativistic and the gravitational accelerations significantly deviate from the Newtonian values.

## 2 Equations

The general relativistic hydrodynamic equations can be cast as a first-order flux-conservative system describing the conservation of mass, momentum and energy. Formulations of this sort are given, e.g. in . In this work we follow the approach laid out in  for a perfect fluid stress-energy tensor $`T^{\mu \nu }`$. The system of equation then reads:

$$\frac{1}{\sqrt{g}}\left(\frac{\sqrt{\gamma }𝐮}{x^0}+\frac{\sqrt{g}𝐟^i}{x^i}\right)=𝐬$$
(1)   
($`x^0=t`$; $`x^i`$ spatial coordinates, $`i=1,2,3`$) where $`𝐮𝐮(𝐰)`$ are the evolved quantities, $`𝐮=(D,S_j,\tau )`$ and $`𝐟^i`$ are the fluxes

$$𝐟^i=(D\left(v^i\frac{\beta ^i}{\alpha }\right),S_j\left(v^i\frac{\beta ^i}{\alpha }\right)+p\delta _j^i,\tau \left(v^i\frac{\beta ^i}{\alpha }\right)+pv^i),$$
(2)   
$`v^i`$ being the 3-velocity and $`p`$ the pressure. The corresponding sources $`𝐬`$ are given by

$`𝐬=(0,T^{\mu \nu }\left({\displaystyle \frac{g_{\nu j}}{x^\mu }}\mathrm{\Gamma }_{\nu \mu }^\delta g_{\delta j}\right),\alpha \left(T^{\mu 0}{\displaystyle \frac{\mathrm{ln}\alpha }{x^\mu }}T^{\mu \nu }\mathrm{\Gamma }_{\nu \mu }^0\right)).`$ (3)   
We note the presence of geometric terms in the fluxes and sources which appear as the local conservation laws of the density current and stress-energy are expressed in terms of partial derivatives. These terms are the lapse function $`\alpha `$, the shift vector $`\beta ^i`$ and the connection coefficients $`\mathrm{\Gamma }_{\nu \mu }^\delta `$ of the 3+1 spacetime metric

$$ds^2g_{\mu \nu }dx^\mu dx^\nu =(\alpha ^2\beta _i\beta ^i)dt^2+2\beta _idx^idt+\gamma _{ij}dx^idx^j$$
(4)   
Additionally $`gdet(g_{\mu \nu })`$ is such that $`\sqrt{g}=\alpha \sqrt{\gamma }`$ and $`\gamma det(\gamma _{ij})`$.

The vector $`𝐰`$, representing the primitive variables, is given by $`𝐰=(\rho ,v_i,\epsilon )`$ where $`\rho `$ is the density and $`\epsilon `$ the specific internal energy. The evolved quantities are defined in terms of the primitive variables as $`D=\rho W`$, $`S_j=\rho hW^2v_j`$ and $`\tau =\rho hW^2pD`$, $`W`$ being the Lorentz factor $`W=(1v^2)^{1/2}`$, with $`v^2=\gamma _{ij}v^iv^j`$, and $`h`$ the specific enthalpy, $`h=1+\epsilon +p/\rho `$. A perfect fluid equation of state $`p=(\mathrm{\Gamma }1)\rho \epsilon `$, $`\mathrm{\Gamma }`$ being the constant adiabatic index, closes the system.

In our computations we specialize the above expressions to the Kerr line element which describes the exterior geometry of a rotating black hole. We use the Kerr-Schild form of the Kerr metric, which is free of coordinate singularities at the black hole horizon. Computations using the more standard Boyer-Lindquist (singular) form of the metric are presented in . Pertinent technical details concerning the specific form of these metrics are given in .

## 3 Numerical scheme

Our hydrodynamical code performs the numerical integration of system (1) using a Godunov-type method. The time update from $`t^n`$ to $`t^{n+1}`$ proceeds according to the following algorithm in conservation form:

$`𝐮_{i,j}^{n+1}=𝐮_{i,j}^n{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x^k}}(\widehat{𝐟}_{i+1/2,j}\widehat{𝐟}_{i1/2,j})+\mathrm{\Delta }t𝐬_{i,j},`$ (5)   
improved with the use of (second-order) conservative Runge-Kutta sub-steps to gain accuracy in time . The numerical fluxes are computed by means of Marquina’s flux-formula . After the update of the conserved quantities the primitive variables are computed via a root-finding procedure.

The flux-formula makes use of the complete characteristic information of system (1), eigenvalues (characteristic speeds) and right and left eigenvectors. Generic expressions are collected in .

The state variables, $`𝐮`$, must be computed (reconstructed) at the left and right sides of a given interface, out of the cell-centered quantities, prior to compute the numerical fluxes. In relativistic hydrodynamics one has the freedom to reconstruct either $`𝐰`$ (primitive variables) or $`𝐮`$ (evolved variables). For efficiency and accuracy considerations we reconstruct the first set, from which the remaining variables are obtained algebraically. The code uses slope-limiter methods to construct second-order TVD schemes by means of monotonic piecewise linear reconstructions of the cell-centered quantities. We use the standard minmod slope which provides the desired second-order accuracy for smooth solutions, while still satisfying the TVD property.

## 4 Results

The classical solution for an asymptotically uniform wind of presureless gas past a compact source (modeled analytically by a point mass) was obtained by . In this solution the material is focussed at the rear part of the object as a result of the gravitational pull. For a pressureless gas, the density at this symmetry line could reach an infinite value and matter would flow on to the hole along this accretion line. However, when pressure is included in the model, a cylindrical shock forms around this line and the accretion proceeds along an accretion column of high density and pressure shocked material. The predicted final accretion pattern consists of a stationary conical shock with the material inside the accretion radius being captured by the central object. An schematic representation of this solution is depicted in Fig. 1.

A numerical evolution of relativistic wind accretion past a rapidly-rotating Kerr black hole ($`a=0.999M`$, $`a`$ specific angular momentum, $`M`$ black hole mass) is depicted in Fig. 2 (left panel). This simulation shows the steady-state pattern in the equatorial plane of the black hole. The tail shock appears stable to tangential oscillations, in contrast to Newtonian simulations with tiny accretors (see, e.g.,  and references there in; see  for a related discussion). The accretion rates of mass and linear and angular momentum also show a stationary behavior (see, e.g., ). As opposed to the non-rotating black hole, in the rotating case the shock becomes wrapped around the central accretor, the effect being more pronounced as the black hole angular momentum $`a`$ increases. The inner boundary of the domain is located at $`r=1.0M`$ (inside the event horizon which, for this model, is at $`1.04M`$) which is only possible with the adopted regular coordinate system. The flow morphology shows smooth behavior when crossing the horizon, all matter fields being regular there.

The enhancement of the pressure in the post-shock zone is responsible for the “drag” force experienced by the accretor. The rotating black hole redistributes the high pressure area, with non-trivial effects on the nature of the drag force. The pressure enhancement is predominantly on the counter-rotating side. We observe a pressure difference of almost two orders of magnitude, along the axis normal to the asymptotic flow direction. The implication of this asymmetry is that a rotating hole moving across the interstellar medium (or accreting from a wind), will experience, on top of the drag force, a “lift” force, normal to its direction of motion (to the wind direction). Although different in origin this feature bears a superficial resemblance with the Magnus effect of classical fluid dynamics.

The right panel of Fig. 2 shows how the accretion pattern would look like were the computation performed using the more common (though singular) Boyer-Lindquist coordinates. The transformation induces a noticeable wrapping of the shock around the central hole. The shock would wrap infinitely many times before reaching the horizon. As a result, the computation in these coordinates would be much more challenging than in Kerr-Schild coordinates, particularly near the horizon. Since the last stable orbit approaches closely the horizon in the case of maximal rotation, the interesting scenario of co-rotating extreme Kerr accretion would be severely affected by the strong gradients which develop in the strong-field region. This will most certainly affect the accuracy and, potentially, also the stability of numerical codes.

### Acknowledgments:

J.A.F. acknowledges financial support from a TMR fellowship of the European Union (contract nr. ERBFMBICT971902).