ACCRETION OF SUPERSONIC WINDS ONTO BLACK HOLES IN 3D: STABILITY OF THE SHOCK CONE

The equations of general relativistic hydrodynamics are written as a flux conservative system, which assumes a standard 3 + 1 decomposition of the spacetime, known as the Valencia formulation (Banyuls et al. 1997):

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{u}}+{\partial }_{{x}^{i}}{{\boldsymbol{F}}}^{(i)}({\boldsymbol{u}})={\boldsymbol{S}}({\boldsymbol{u}}),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{u}}=\{D,{S}_{i},\tau \}$ is the set of conserved variables, D is the generalized rest-mass density of the fluid, S_i are the generalized momenta in each direction, τ is the internal energy, ${{\boldsymbol{F}}}^{(i)}({\boldsymbol{u}})$ is the fluxes, and ${\boldsymbol{S}}({\boldsymbol{u}})$ is the source. In terms of the primitive variables of the fluid elements, ρ the rest-mass density, p the pressure, vⁱ the fluid 3-velocity, and the specific internal energy, the conserved variables are defined by

$$\begin{eqnarray}D & = & \sqrt{\gamma }W\rho \\ {S}_{i} & = & \sqrt{\gamma }\rho {{hW}}^{2}{v}_{i}\\ \tau & = & \sqrt{\gamma }\left(\rho {{hW}}^{2}-p\right)-D,\end{eqnarray} \tag{ 2 }$$

where γ is the determinant of the spatial 3-metric ${\gamma }_{{ij}}$ of the 3D spatial slices that the spacetime is foliated with, $h\equiv 1+\epsilon +p/\rho$ is the specific enthalpy and $W={(1-{\gamma }_{{ij}}{v}^{i}{v}^{j})}^{-1/2}$ is the Lorentz factor. As usual, the system of equations is closed with an equation of state. In terms of primitive and conservative variables, the fluxes and sources are:

$$\begin{eqnarray*}\ {{\boldsymbol{F}}}^{i}({\boldsymbol{u}}) & = & \left(\begin{array}{c}(\alpha {v}^{i}-{\beta }^{i})D\\ (\alpha {v}^{i}-{\beta }^{i}){S}_{j}+p{\delta }_{j}^{i}\\ (\alpha {v}^{i}-{\beta }^{i})\tau +{{pv}}^{i}\end{array}\right),\\ \ {\boldsymbol{S}}({\boldsymbol{u}}) & = & \left(\begin{array}{c}0\\ {T}^{\mu \nu }\left({\partial }_{\mu }{g}_{\nu j}+{{\rm{\Gamma }}}_{\mu \nu }^{\delta }{g}_{\delta j}\right)\\ \alpha \left({T}^{\mu 0}{\partial }_{\mu }\mathrm{ln}\alpha -{T}^{\mu \nu }{{\rm{\Gamma }}}_{\nu \mu }^{0}\right)\end{array}\right),\end{eqnarray*}$$

where ${g}_{\mu \nu }$ are the components of the 4-metric and ${{\rm{\Gamma }}}^{\delta }{}_{\mu \nu }$ are the Christoffel symbols of the spacetime, α is the lapse function, and ${\beta }^{i}$ are the components of the shift vector associated to the 3 + 1 decomposition of the spacetime.

We solve the system of Equations (1) using the Cactus Einstein Toolkit (ETK; Einstein Toolkit 2005; Löffler et al. 2012), which provides the necessary computational tools to evolve this relativistic fluid in full 3D. In particular we use the GRHydro Thorn (Baiotti et al. 2005), which contains high-resolution shock capturing methods to solve these equations. For these methods we specifically use the HLLE and Marquina numerical flux formulae and minmod and MC linear reconstructors. For the integration in time we use a fourth order Runge–Kutta method. As part of the diagnostics tools, we use the Outflow Thorn that allows us to measure the flux of rest-mass across a given spherical surface and monitor the behavior of the accretion rate.

We describe the spacetime of the rotating black hole using Kerr–Schild (KS) horizon-penetrating coordinates, with the axis of rotation parallel to $\hat{z}$ as follows:

$$\begin{eqnarray}{{ds}}^{2} & = & \left({\eta }_{\mu \nu }+\displaystyle \frac{2{{Mr}}^{3}}{{r}^{4}+{a}^{2}{z}^{2}}{l}_{\mu }{l}_{\nu }\right){{dx}}^{\mu }{{dx}}^{\nu },\\ {l}_{\mu } & = & \left(1,\displaystyle \frac{{rx}+{ay}}{{r}^{2}+{a}^{2}},\displaystyle \frac{{ry}-{ax}}{{r}^{2}+{a}^{2}},\displaystyle \frac{z}{r}\right),\\ r & = & \sqrt{\displaystyle \frac{{r}_{*}^{2}-{a}^{2}+\sqrt{{({r}_{*}^{2}-{a}^{2})}^{2}+4{a}^{2}{z}^{2}}}{2}},\\ {r}_{*} & = & \sqrt{{x}^{2}+{y}^{2}+{z}^{2}},\end{eqnarray} \tag{ 3 }$$

where M is the mass of the black hole and we have assumed geometric units $G=c=1$ to hold.

The fluid is set initially as a spatially constant rest-mass density ideal gas, moving toward the black hole and with a given asymptotic velocity ${v}_{\infty }^{2}={v}^{i}{v}_{i}.$ We assume the fluid obeys a gamma-law equation of state $p=({\rm{\Gamma }}-1)\rho \epsilon$ and thus we introduce the asymptotic speed of sound ${c}_{s\infty }.$ Once we define the value of c_s∞ and assume the density to initially be a constant $\rho ={\rho }_{\mathrm{ini}},$ the pressure can be written as ${p}_{\mathrm{ini}}={c}_{{\rm{s}}\infty }^{2}{\rho }_{\mathrm{ini}}/({\rm{\Gamma }}-{c}_{{\rm{s}}\infty }^{2}{{\rm{\Gamma }}}_{1}),$ where ${{\rm{\Gamma }}}_{1}={\rm{\Gamma }}/({\rm{\Gamma }}-1).$ In order to avoid negative and zero values of the pressure, the condition ${c}_{{\rm{s}}\infty }\lt \sqrt{{\rm{\Gamma }}-1}$ has to be satisfied. Finally, with this value for ${p}_{\mathrm{ini}},$ the initial internal specific energy is reconstructed using the equation of state.

In order to parameterize the initial data, we define the relativistic Mach number at infinity ${{\mathcal{M}}}_{\infty }^{R}={{Wv}}_{\infty }/{W}_{s}{c}_{{\rm{s}}\infty },$ where W is the Lorentz factor of the gas and W_s is the Lorentz factor calculated with the speed of sound. When ${{\mathcal{M}}}_{\infty }^{R}$ is bigger than one, the wind is said to be supersonic and otherwise subsonic.

According to Foglizzo et al. (2005), the two parameters deciding on the stability of the process are the Mach number of the gas ${{\mathcal{M}}}_{\infty }^{R}$ and the relative size of the accretor with respect to the Bondi accretion radius r_acc. In our case the compact object is a black hole and we define an approximate black hole radius to be ${r}_{\mathrm{hor}}=2\;M,$ even though our hole is rotating. Thus the relative size of the accretor is given by ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}={v}_{\infty }^{2}+{c}_{{\rm{s}}\infty }^{2}.$

This formula indicates that a relatively small black hole is equivalent to having a case with the wind moving with a small velocity v_∞. This explains why the cases where ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}$ is small are difficult to track, and even more difficult in 3D. The reason is that the numerical domain is set such that it contains at least a sphere of radius r_acc in order to see the accretion process. Slower winds imply a bigger r_acc, therefore demanding a bigger numerical domain and consequently more computational resources.

In order to observe the behavior of the process in a general scenario, we choose five representative different directions of the wind, which we characterize in terms of its orientation with respect to $\hat{z}$, which is the axis of rotation of the black hole denoted by $\Uparrow .$ These five different orientations of the wind are represented by $\uparrow \;\downarrow \;\to \;\swarrow \;\nwarrow ,$ and correspond to directions parallel to $\hat{z},\;-\hat{z},\;\hat{x},\;-\hat{x}-\hat{z}$, and $-\hat{x}+\hat{z}$ respectively.

The ETK allows the solution of the equations in Cartesian coordinates, and thus the domain is a cubic box of a given size. In the case of the evolution of a wind, two boundary conditions are required (see, e.g., Dönmez et al. 2011; Cruz-Osorio et al. 2012), namely an inflow boundary condition in the upstream region of the domain, where the gas is being pumped into the box, and outflow boundary conditions in the piece of the boundary where the wind is expected to exit the numerical domain.

The ETK does not provide the tools to automatically incorporate different boundary conditions on different faces of the cubic domain and we programmed a module that allows one to do it.

There is another very important boundary to consider. Since we are using horizon-penetrating coordinates we are allowed to perform an excision inside the black hole horizon as performed in Cruz-Osorio et al. (2012) and Lora-Clavijo & Guzmán (2013). The excision method consists of the removal of a chunk of the numerical domain inside the black hole horizon. The removed chunk is usually assumed to be a lego sphere that defines a lego-shaped boundary inside the horizon. Considering the fact that the light cones in KS coordinates remain open and point toward the singularity inside the black hole horizon, no boundary conditions are needed, and extrapolation of the fluid variables suffices to guarantee that the fluid is not reflected back toward the exterior of the black hole. This method is already implemented in the ETK, as described in Hawke et al. (2005).

In the first set of experiments we use the following parameters: ${\rho }_{0}={10}^{-6},$ ${\rm{\Gamma }}=4/3,$ c_s∞ = 0.1. The velocities are supersonic with ${{\mathcal{M}}}_{\infty }=3,\;5$ and we use two values of the wind velocity ${v}_{\infty }=0.3,\;0.5,$ which correspond to the sizes of the accretor ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.1$ and 0.26, respectively. These parameters correspond to mildly small black holes with sizes bigger than those expected for objects showing FF instability according to Foglizzo et al. (2005). We also carried out these simulations for two values of the black hole angular momentum $a=0.5,\;0.95.$

The evolution was performed using resolution ${\rm{\Delta }}{xyz}=0.25\;M,$ on a domain ${[-40\;M,40\;M]}^{3},$ where M is the mass of the black hole. This is a modest resolution compared to that used in the analysis using s-lab symmetry (Font & Ibañez 1998; Dönmez et al. 2011; Cruz-Osorio et al. 2012); however, this should trigger an instability faster, if it triggers any at all. In all of the cases, that is, for the two velocities, the two values of the black hole angular momentum and the five orientations of the wind, we did not find the FF instability, as expected for these parameters. In Table 1 we list the set of configurations for which we track the evolution of the wind.

Table 1.  Configurations Analyzed for ${\rho }_{0}={10}^{-6},$ ${\mathcal{M}}=3,\;5,$ ${\rm{\Gamma }}=4/3,$ c_s∞ = 0.1, and ${v}_{\infty }=0.3$

Config.	${\mathcal{M}}=3$			${\mathcal{M}}=5$
	${\dot{M}}_{3}({10}^{-4})$	${\dot{M}}_{3\mathrm{BHL}}({10}^{-4})$		${\dot{M}}_{5}({10}^{-4})$	${\dot{M}}_{5\mathrm{BHL}}({10}^{-4})$
a = 0.5
$\uparrow \Uparrow$	8.29	⋯		3.81	⋯
$\downarrow \Uparrow$	8.24	⋯		3.80	⋯
$\to \Uparrow$	8.21	4.00		3.63	1.00
$\swarrow \Uparrow$	8.27	⋯		3.60	⋯
$\nwarrow \Uparrow$	8.32	⋯		3.61	⋯
a = 0.95
$\uparrow \Uparrow$	8.23	⋯		3.64	⋯
$\downarrow \Uparrow$	8.22	⋯		3.64	⋯
$\to \Uparrow$	8.82	4.00		3.63	1.00
$\swarrow \Uparrow$	8.57	⋯		3.64	⋯
$\nwarrow \Uparrow$	8.50	⋯		3.58	⋯

Note. The value of $\dot{M}$ corresponds to the value measured on a spherical surface located at $r=3\;M$ after the process of shock cone formation has been settled down. No instability was found in all of the presented cases. In order to compare how well the BHL formula works, we also show the accretion rate for each case.

Download table as:  ASCIITypeset image

In these simulations there are already some previously unknown results, for example, the case $\to \Uparrow ,$ which is the generalization of the s-lab studies shown in, e.g., Font & Ibañez (1998); Dönmez et al. (2011); and Cruz-Osorio et al. (2012). In Dönmez et al. (2011) it was shown that there was a FF instability, whereas in Cruz-Osorio et al. (2012) the use of horizon-penetrating coordinates to describe the black hole indicated that there is no such instability. However, both analyses were performed using s-lab symmetry and the result in 3D was still uncertain. We show, for the first time, in Figure 1, a snapshot of the cone in 3D for the case with ${v}_{\infty }=0.5,$ ${c}_{s\infty }=0.1$ and black hole angular momentum a = 0.95. No signs of instability were found up to t = 10,000 M, which is a comparable time with that in Font & Ibañez (1998; Dönmez et al. (2011); and Cruz-Osorio et al. (2012).

Zoom In Zoom Out Reset image size

Aside from the morphology of the shock cone that looks very stable, we monitor the stability of the accretion process by measuring the accretion rate $\dot{M}$ at a spherical surface located at $r=3\;M.$ In all of our experiments we find this rate to stabilize after an initial transient corresponding to the shock cone formation lapse. We measure this rate for the different cases and check that the slower the wind, the bigger the accretion rate, as expected from the Bondi formula. We also show the accretion rate in Table 1 for the two cases ${{\mathcal{M}}}_{\infty }=3$ and ${{\mathcal{M}}}_{\infty }=5$, and for comparison we additionally show the values obtained using the BHL accretion rate formula (Edgar 2004)

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{BHL}}=\displaystyle \frac{4\pi \rho }{{({v}_{\infty }^{2}+{c}_{s\infty }^{2})}^{3/2}};\end{eqnarray} \tag{ 4 }$$

we label the results with ${\dot{M}}_{3\mathrm{NHL}}$ and ${\dot{M}}_{5\mathrm{NHL}}$ for the two wind velocities considered. We point out that the accretion rate, as measured for a black hole spacetime in our simulations for this size of black hole, is considerably bigger than the accretion rate calculated with the BHL formula.

3.1.1. Particular Cases

Some special configurations may be interesting. For instance one may wonder whether there is a whirlpool behind the black hole when the wind moves parallel or antiparallel to the axis of rotation. In Figure 2 we show slices of the velocity field at three planes of z = constant, and a lateral view of the process at $t=8000\;M.$ In this case the black hole size is ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.26$ and no dragging effects are noticeable in the velocity field due to the rotation of the black hole. What can actually be seen is the convergence of the stream lines at the shock cone boundary.

Zoom In Zoom Out Reset image size

Another case that has never been shown in 3D is the accretion of a diagonal wind. We show a snapshot of the shock cone density at $t=5000\;M$ in Figure 3, for the wind orientation $\swarrow \Uparrow .$

Zoom In Zoom Out Reset image size

A second set of tests consists of the evolution of the wind for a smaller black hole, particularly three cases: ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.011,0.012,0.025.$ Because these numbers are smaller than the threshold of 0.05, the FF instability is expected to happen as suggested in Foglizzo et al. (2005). According to the analysis in Foglizzo et al. (2005), even axisymmetric simulations should become unstable when considering ${\rm{\Gamma }}=5/3,4/3$ and ${{\mathcal{M}}}_{\infty }=2.4,\;3.0.$ We evolved these particular cases in 3D for all the wind orientations described before. These cases are presented in Table 2. In this case we only perform the simulations for the bigger black hole angular momentum a = 0.95. We use the same resolution ${\rm{\Delta }}{xyz}=0.25\;M,$ however, because a domain containing at least a sphere of radius r_acc is required, we used the domain ${[-100\;M,100\;M]}^{3}.$ Once again we did not find the FF instability after 10,000 M. We point out that unlike the mildly small accretor size, the accretion rate for this size of accretor is considerably smaller than the accretion calculated with the BHL formula.

Table 2.  Configurations Analyzed for ${\rho }_{0}={10}^{-6}$ and a = 0.95

${\rm{\Gamma }}=5/3$	${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.011$	${\mathcal{M}}=2.4$
	${\dot{M}}_{2.4}(\times {10}^{-2})$	${\dot{M}}_{\mathrm{BHL}}({10}^{-2})$
$\uparrow \Uparrow$	⋯	⋯
$\downarrow \Uparrow$	⋯	⋯
$\to \Uparrow$	0.14	1.00
$\swarrow \Uparrow$	⋯	⋯
$\nwarrow \Uparrow$	⋯	⋯
${\rm{\Gamma }}=4/3$	${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.012$	${\mathcal{M}}=2.5$
	${\dot{M}}_{2.5}(\times {10}^{-2})$	${\dot{M}}_{\mathrm{BHL}}({10}^{-2})$
$\uparrow \Uparrow$	⋯	⋯
$\downarrow \Uparrow$	⋯	⋯
$\to \Uparrow$	0.48	1.00
$\swarrow \Uparrow$	⋯	⋯
$\nwarrow \Uparrow$	⋯	⋯
${\rm{\Gamma }}=4/3$	${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.025$	${\mathcal{M}}=3$
	${\dot{M}}_{3.0}(\times {10}^{-2})$	${\dot{M}}_{\mathrm{BHL}}({10}^{-3})$
$\uparrow \Uparrow$	⋯	⋯
$\downarrow \Uparrow$	⋯	⋯
$\to \Uparrow$	0.34	3.72
$\swarrow \Uparrow$	⋯	⋯
$\nwarrow \Uparrow$	⋯	⋯

Note. The value of $\dot{M}$ corresponds to the value measured on a spherical surface located at $r=3\;M$ after the process of shock cone formation has been settled down. None of the cases show instability. In order to compare how well the BHL formula works, we also show the accretion rate for each case. Unlike in the case of mildly small black holes, in this case the numerical accretion rate is about the same, independent of the orientation of the wind.

Download table as:  ASCIITypeset image

In order to illustrate how the dragging effects due to the rotation of a small black hole on the fluid are more significant compared to the bigger black hole case, in Figure 4 we show a snapshot of the density for the case with orientation $\to \Uparrow ,$ ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.011,$ ${\rm{\Gamma }}=5/3$, and ${\mathcal{M}}=2.4$ from Table 2. In Figure 5 we also show a snapshot on a vertical plane of the diagonal case, which shows a more dragged fluid than in the case of Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we did not find any unstable case, in addition to the analysis of simply tracking the evolution of the gas and expecting the numerical error to trigger the FF instability, we designed a numerical experiment in which we evolve the gas until the shock cone forms and then apply a perturbation to the velocity field of the form

$$\begin{eqnarray}&&{v}^{i}={v}^{i}+k\delta (t),\end{eqnarray} \tag{ 5 }$$

where $k\geqslant 0$ is a random number between 0 and A. We applied this random perturbation perpendicular to the direction of the shock cone and expected for it to trigger the instability. We performed this experiment using the two extreme relative accretor sizes in this paper, ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.011$ and ${r}_{\mathrm{hor}}/{r}_{\mathrm{acc}}=0.26,$ and all the orientations we have considered here. In all of the cases no instability was triggered for all the directions of the wind.

In Figure 6 we show the accretion rate for the perturbed case, where a feature can be seen at the moment at which we apply the perturbation. We used $A=0.2,$ which means that the velocity perpendicular to the shock cone has been added randomly, even one-fifth of the speed of light, and nevertheless the flux restores and becomes stable. In Figures 7 and 8 we show snapshots of the density before, during, and after the perturbation is applied, and show how the morphology of the cone restores.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We wanted to keep the physical system as free as possible of numerical sophistication in order to avoid undesirable effects. Thus, in the production runs we used a single refinement level because we wanted to avoid the numerical dissipation required to damp instabilities at refinement level interfaces, which could hide the potential unstable modes of the shock cone. This has been discussed as a potential damping of the instability in FF-type phenomena (Foglizzo et al. 2005).

In order to show the effects due to the use of fixed mesh refinement (FMR) in Figure 9 we show the accretion rate when using FMR with two refinement levels and when using a single refinement level. The deviation is of the order of 3% in the mass accretion rate. Thus, even though FMR allows  more efficient execution of long runs, the accuracy is not a good as with unigrid, which is important in long runs.

Zoom In Zoom Out Reset image size

We also verified our results on the stability of the shock cone in all our runs using various combinations of approximate Riemann solvers and reconstruction methods. We evolved all the configurations in the tables using four combinations, namely: (1) Marquina with MC and minmod and (2) HLLE with MC and minmod. In all of the cases we found the shock cone to be stable in all of the described orientations and velocities of the wind. In Figure 9 we compare the accretion rate estimated using the four combinations of fluxes and reconstructors for a particular case. The difference among the accretion rate lies within 1%, which is smaller than the error introduced when using FMR.