Radiation GRMHD Simulations of the Hard State of Black Hole X-Ray Binaries and the Collapse of a Hot Accretion Flow

Accreting black holes in X-ray binary systems (BHBs) show distinct spectral and variability states (e.g., Remillard & McClintock 2006; Done et al. 2007; Belloni & Motta 2016). The "soft" states show thermally dominated spectra peaking at a photon energy of ≃1 keV, broadly consistent with optically thick emission from a geometrically thin accretion disk (Shakura & Sunyaev 1973). The spectrum in the "hard" states is instead an inverted (hard) power law rising to a cutoff at ≃100 keV, thought to arise from inverse Compton scattering in hot, optically thin plasma (e.g., Sunyaev & Titarchuk 1980). Synchrotron emission from relativistic jet electrons may also contribute (Markoff et al. 2001).

BHB hard states are observed over a wide range of sub-Eddington luminosities, 10⁻⁴ ≲ L/L_Edd ≲ 10⁻¹ (e.g., Maccarone 2003). At the low end of this regime the infalling gas is collisionless with virial (hot) ions and colder electrons (e.g., Ichimaru 1977; Rees et al. 1982; Narayan & Yi 1994). At the high end, the radiative cooling time should become shorter than the inflow time (Rees et al. 1982). The resulting collapse of the hot accretion flow to a thin disk has long been associated with hard to soft state transitions (e.g., Esin et al. 1997). Both the radiative efficiency and maximum accretion rate of hot accretion flows remain uncertain. Their values in analytic theory depend on the use of simplified prescriptions for angular momentum transport, electron heating, and radiative cooling (e.g., Yuan & Narayan 2014).

It is now possible to address some of these shortcomings using MHD accretion theory. Radiative general relativistic MHD (GRMHD) simulations of weakly magnetized hot accretion flows with turbulent electron heating prescriptions have previously found low radiative efficiencies ≲1% and maximum luminosities L/L_Edd ≲ 10⁻⁴ (Ryan et al. 2017; Sadowski et al. 2017), far below that observed in the hard spectral state of BHBs. In the limit of strong magnetization and high black hole spin parameter, radiation GRMHD models of M87 have found higher radiative efficiency (Chael et al. 2019).

Here we perform 3D, two-temperature, radiation GRMHD simulations of strongly magnetized accretion onto a spinning stellar mass black hole with frequency-dependent Monte Carlo radiative transfer (Section 2). We show that for these parameter choices, equilibrium states reach higher luminosities than previously possible, L/L_Edd ∼ 10⁻⁴–10⁻². For L/L_Edd ≳ 10⁻³, radiative cooling reduces the scale height of the accretion flow (Section 3). The emergent spectrum forms an apparent inverted power law, typical of that seen in BHB hard states. At L/L_Edd ∼ 10⁻², our hot accretion flow solution appears to undergo a cooling runaway and collapse (Section 4). We compare our results with hard state phenomenology and analytic theory in Section 5, and briefly discuss the implications for observed state transitions in BHB systems.

We carried out radiation GRMHD simulations using the public code ebhlight ³ (Ryan et al. 2015, 2017). ebhlight uses the HARM scheme for conservative ideal GRMHD (Gammie et al. 2003; Noble et al. 2006) and a Monte Carlo treatment of the anisotropic, frequency-dependent radiation field. The electron entropy is evolved separately (Ressler et al. 2015), with heating contributions from Coulomb collisions (Sadowski et al. 2017) and grid-scale dissipation. Assuming a thermal distribution, the electrons interact with photons via synchrotron and bremsstrahlung emission and self-absorption and Compton scattering.

The simulations were initialized from a gas torus (Fishbone & Moncrief 1976) with an inner radius of r_in = 12 r_g , a pressure maximum radius of r_max = 25 r_g , and a dimensionless black hole spin parameter a = 0.9375, where r_g = GM/c² and we set M = 10M_⊙. The grid resolution of 320 × 256 × 160 cells in modified spherical Kerr–Schild coordinates was chosen to adequately resolve both the magnetorotational instability (Dexter et al. 2020a) and the time-dependent evolution of the mass accretion rate and magnetic flux accumulation onto the black hole (Dexter et al. 2020b). The torus was seeded with poloidal magnetic field, whose topology was chosen to saturate the magnetic flux on the black hole (a magnetically arrested disk, or MAD state; Bisnovatyi-Kogan & Ruzmaikin 1974; Igumenshchev et al. 2003; Tchekhovskoy et al. 2011). We adopt an electron heating prescription from particle-in-cell calculations of magnetic reconnection with a modest guide field (Werner et al. 2018). According to this prescription, electrons receive f_e = 1/4 − 1/2 of the total dissipated energy depending on the magnetization parameter σ ≡ b²/ρ, where b and ρ are the magnetic field strength and rest mass density in GM = c = 1 units.

We first ran a nonradiative GRMHD simulation for t = 9 × 10³ r_g /c. Radiation was then initialized, with a simulation mass scale M_unit chosen to approximately match target accretion rate values of $\dot{m}=\dot{M}/{\dot{M}}_{\mathrm{Edd}}={10}^{-5}$, 10⁻⁴, and 10⁻³, where $0.1{\dot{M}}_{\mathrm{Edd}}{c}^{2}={L}_{\mathrm{Edd}}$. We subsequently evolved both the radiative and nonradiative simulations for an additional (1−2) × 10³ r_g /c, sufficient to reach radiative and inflow equilibrium over the same range in radius (≲20 r_g ), where the inflow equilibrium radius at time t is defined by t = r_eq/∣v^r ∣ (e.g., Narayan et al. 2012), and v^r = u^r /u^t with u^μ the coordinate four-velocity. Superphoton packets reaching r = 40 r_g were removed from the grid and recorded, accounting for the gravitational redshift at that location (Ryan et al. 2015). No radiative processes were considered outside that radius. Table 1 lists time-averaged properties of our calculations over the final 500 r_g /c. Radially dependent quantities are averaged over the range 3−15 r_g . Angle brackets denote shell averages:

$$\begin{eqnarray}&&\langle x\rangle =\displaystyle \frac{\int \int d\theta d\phi \sqrt{-g}\rho \,x}{\int \int d\theta d\phi \sqrt{-g}\rho },\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\langle x{\rangle }_{J}=\displaystyle \frac{\int \int d\theta d\phi \sqrt{-g}J\,x}{\int \int d\theta d\phi \sqrt{-g}J},\end{eqnarray} \tag{ 2 }$$

where $\sqrt{-g}$ is the metric determinant and J is the cooling rate per volume. In practice the integration is carried out in code coordinates.

Table 1. Averaged Properties of the Final 500 r_g /c of Our ebhlight Simulations

Name	Munit (1010 g)	Nsph (109)	$\dot{m}$ (10−3)	〈H〉/r	〈β〉	〈Qem〉	〈Qsc〉	〈θe 〉J	L/LEdd (10−3)	Lsc/Lem	Γ3–20	Epeak (keV)	$L/\dot{M}{c}^{2}$
M5	1	0.50	0.02	0.20	3.44	922	342	10.58	0.02	0.7	2.10	0.01	0.08
M4	10	1.42	0.44	0.22	2.57	2318	196	4.33	0.79	5.5	1.92	57.61	0.18
M3	17	1.48	1.34	0.15	0.72	2129	159	2.84	3.22	18.4	1.76	121.75	0.24
M3f	30	1.72	1.57	0.13	0.90	2308	122	1.52	6.22	36.2	1.66	121.75	0.40

Note. The mass unit M_unit and average number of superphoton packets N_sph are independent variables, while the other quantities are simulation outcomes. We italicize the M3f results since that simulation fails to achieve a steady state.

Download table as:  ASCIITypeset image

The ebhlight radiation field is represented by a number of weighted superphoton packets, N_sph, chosen so that the average number of emission and scattering events per cooling time is large, $Q\equiv \dot{N}{u}_{e}/J\gg 1$, where $\dot{N}$ is the average rate of emission or scattering events and u_e is the electron internal energy density. The simulations presented here used a total of N_sph ≃ (5−30) × 10⁸, or ≃80–450 superphoton packets per grid cell. The quality factors for both emission (Q_em) and scattering (Q_sc) interactions are ≳100. Such large Q values should result in convergence of the radiation field and electron thermodynamics (Yao et al. 2021).

We parallelized the radiation runs using 1 MPI process per node, and OpenMP across each node. The domain is split into hemispheres, each containing the full radial grid. This results in an equal number of superphoton packets on all nodes and at all times to within ±5% when using 64 nodes. Our radiative simulations were ≃5–30× more computationally expensive than nonradiative GRMHD simulations performed with ebhlight at the same grid resolution. The cost increased with $\dot{m}$, due to both the higher number of superphoton packets used to resolve the radiation field and the higher rate of absorption and scattering interactions.

At low $\dot{m}$, a new radiative equilibrium is reached rapidly. The radiation is produced close to the black hole where timescales are short, and the cooling primarily modifies the electron temperature (e.g., Ryan et al. 2017, 2018; Chael et al. 2019). At higher $\dot{m}$, cooling changes the accretion flow structure and equilibrium is only reached after an inflow time from the outermost radius of interest. We determine the time range over which a new steady state has been reached by waiting for the flow to establish stationary radial profiles of fluid quantities, including the total cooling rate and gas and electron temperature. By these criteria, M3 appears to reach a stable equilibrium for r ≲ 20 r_g within ≲2000 r_g /c. For r > 20 r_g , the gas is still cooling and inflow equilibrium has not yet been reached. We note that this procedure is complicated by turbulent fluctuations in $\dot{M}$ on comparable timescales to our time averaging, which can produce secular changes in the normalization of J and τ_sc.

Figure 1 shows radial profiles of vertically integrated scattering optical depth and density-weighted shell averages of the electron heating ratio from Coulomb collisions and viscous dissipation (〈Q_coul〉/f_e 〈Q_visc〉), the gas and electron temperatures, and the density scale height. All runs are optically thin to electron scattering in the vertical direction, although the M3 run approaches τ_sc ∼ 1 near r ≃ 20 r_g . For that simulation, electron heating from Coulomb collisions exceeds that from grid-scale dissipation for r ≳ 3 r_g . The electron temperatures are far below those from nonradiative GRMHD: radiative cooling is important in all cases. For $\dot{m}\gtrsim {10}^{-4}$, the average electron temperature is nonrelativistic. In the M3 simulation, radiative cooling reduces the gas temperature and pressure and in turn the density scale height.

Zoom In Zoom Out Reset image size

Azimuthally averaged snapshots of fluid and radiation quantities are shown in Figure 2 for each simulation reaching a steady state. The M5 and M4 models structurally are still geometrically thick, hot accretion flows. The high magnetization chosen here results in plasma β ≡ p_gas/p_B ≳ 1 in the disk and ≪1 outside of it, where p_gas and p_B = b²/2 are the gas and magnetic pressure. For model M3, radiative cooling of the gas results in the formation of a strongly magnetized (β < 1), denser region close to the midplane for all radii r ≲ 20 r_g , with a scale height 〈H〉/r ≈ 0.15.

Zoom In Zoom Out Reset image size

In all cases, the radiative cooling by both emission (J_em) and scattering (J_sc) are concentrated toward the equatorial plane in the dense accretion flow, rather than in a surrounding corona or jet. This is partly by construction: we do not allow radiation from magnetically dominated regions where σ > 1. Compton cooling becomes increasingly important at higher $\dot{m}$, as can be seen in plots of the Compton y_C parameter (calculated using τ_sc in the radial direction, right middle panels). A region with y_C > 1 forms in the M4 and M3 simulations near the equatorial plane and extends out to large radii.

Time-averaged spectra integrated over a solid angle are shown in Figure 3. The thin lines show the separate emission and inverse Compton scattering contributions. Synchrotron radiation is the dominant emission process in all cases and peaks in the optical band (∼10¹⁵ Hz). A secondary peak due to optically thin bremsstrahlung grows in relative strength with increasing $\dot{m}$ (secondary hard X-ray emission peak near ∼10¹⁹ Hz).

Zoom In Zoom Out Reset image size

For the M5 model, inverse Compton scattering produces a subdominant contribution to the bolometric luminosity (L_sc/L_em ≲ 1 in Table 1). For the M4 and M3 models, repeated inverse Compton scattering of synchrotron seed photons is the dominant cooling mechanism. The result is a hard power-law X-ray spectrum extending to a thermal cutoff at a photon energy of ≃100 keV. The radiative efficiency also increases with $\dot{m}$. Its value of ≈24% for the M3 model is somewhat higher than that of a Novikov & Thorne (1973) thin accretion disk for this spin parameter.

For simulation M3, Coulomb coupling is a significant source of electron heating and the gas cools efficiently by inverse Compton scattering. The resulting, apparent equilibrium is poised at the edge of instability: both the scattering optical depth and Coulomb electron heating ratio are approaching unity. Once those values are exceeded, no stable two-temperature solution is expected (e.g., Rees et al. 1982; Yuan & Narayan 2014).

We demonstrate this explicitly using a simulation with a slightly higher mass unit (M3f in Table 1). This model reached an initial radiative equilibrium at a small radius r ≲ 10 r_g . The optical depth increases to τ_sc ∼ 1, and Coulomb collisions dominate the electron heating. A cooling front driven by inverse Compton scattering then propagates inward from r ≃ 20 r_g (Figure 4). The loss of pressure support is strong enough to decrease the gas density scale height, while efficient Coulomb coupling drives the ion to electron temperature ratio to unity.

Zoom In Zoom Out Reset image size

As the gas temperature decreases, the total integrated cooling rate increases. Apparently the flow is thermally unstable. We are unable to follow the further evolution past the final state shown in Figure 4. According to analytic theory, the likely result is collapse to an optically thick and geometrically thin disk.

We have presented 3D radiation GRMHD simulations of accretion onto a 10M_⊙ black hole as a function of the mass accretion rate. Our model choices of saturated magnetic flux, high spin, and efficient electron heating result in radiatively efficient accretion flows that reach luminosities of L/L_Edd ≲ 10⁻² while remaining optically thin. They provide relativistic, 3D, MHD realizations of physical regimes studied extensively using 1D accretion theory. They also exhibit radiative properties seen in a wealth of BHB observations.

In the hard state, BHBs show power-law X-ray spectra that harden from a photon index ⁴ Γ = 2 to 1.5 for increasing L/L_Edd ∼ 10⁻⁴ to 10⁻² (e.g., Remillard & McClintock 2006; Skipper & McHardy 2016). We find time-averaged photon indices of Γ = 1.93 and 1.77 from 3 to 20 keV for models M4 and M3, respectively, in good agreement with observations. The increasing spectral hardness is expected from analytic Comptonization models (e.g., Sunyaev & Titarchuk 1980; Haardt & Maraschi 1993). The spectrum is exponentially cut off at a thermal energy of ≳100 keV, again broadly in agreement with observations of the hard state (e.g., Grove et al. 1998) and without requiring any nonthermal particle acceleration.

The cooling is concentrated toward the inner radii in all cases. The average emission radius of 〈r〉_J ≲ 6 r_g is consistent with observations of a compact X-ray emitting region (Dai et al. 2010; Uttley et al. 2014). A small amount of additional luminosity would likely be produced from the neglected region r > 40 r_g . The radiated luminosity is produced by a small fraction of the total volume, and the electron temperature weighted by emissivity, 〈θ_e 〉_J , exceeds that weighted by density, 〈θ_e 〉 (Table 1 and Figure 1). For the M3 model, 〈θ_e 〉_J ≃ 3 while 〈θ_e 〉 ≃ 0.2 at the average emission radius.

The HARM algorithm requires the (artificial) injection of mass and internal energy in order to limit the fluid magnetization, set here to σ ≤ 50. We neglect electron–photon interactions whenever σ > 1 to ensure that this material does not contribute to the radiation field. The chosen cutoff value is arbitrary. At low $\dot{m}$, 〈σ〉_J ≲ 1 for the M5 model and higher σ material would likely cause order-unity changes to the emergent spectrum and bolometric luminosity (see Chael et al. 2019). As inverse Compton cooling from larger radii becomes more important at high $\dot{m}$, 〈σ〉_J decreases and the results for the M4 and M3 models should be less sensitive to the chosen cutoff value.

Axisymmetric radiative GRMHD simulations have been carried out using ebhlight for $\dot{m}={10}^{-9}\mbox{--}{10}^{-5}$ for supermassive black holes (M = 10⁸ M_⊙ or ≈(3–6) × 10⁹ M_⊙; Ryan et al. 2017, 2018). Those calculations reached L/L_Edd ≲ 10⁻⁴. Qualitatively, our results are similar. The effects of radiative cooling and Coulomb coupling become more important with increasing $\dot{m}$, and a significant fraction of the radiated luminosity originates from inverse Compton scattering at larger radii. We find radiative efficiencies a factor ≃10 higher at $\dot{m}$ = 10⁻⁵. This is most likely due to higher accretion flow temperatures in our models, resulting from the interplay between strong magnetization (lower plasma β) and our choice of a more uniform electron heating prescription. For example, the emergent radiation here is produced by the dense accretion flow rather than the jet wall. Similar radiative efficiencies of ≃10% have been found in 3D radiation GRMHD models of strongly magnetized accretion onto a rapidly spinning black hole for conditions relevant to M87 (Chael et al. 2019; Yao et al. 2021). Avara et al. (2016) further showed that the radiative efficiency of MAD accretion can significantly exceed the Novikov & Thorne (1973) value, as we find for model M3. Kinch et al. (2021) found a similar hard, power-law spectrum and high radiative efficiency in a post-processed model including a weakly magnetized, geometrically thin disk for $L/{L}_{{\rm{Edd}}}\sim {10}^{-2}$. For the accretion rate range here, our models show ≃3–5× higher radiative efficiencies than found in analytic models (Xie & Yuan 2012). Our limiting value of $\dot{m}\approx {10}^{-3}$ is also about an order of magnitude higher than Ryan et al. (2017) estimated. We see similar differences in the Coulomb to viscous heating ratio between our models and theirs, and cooling may also be somewhat more important at fixed $\dot{m}$ for higher M.

For L/L_Edd ≲ 10⁻², our hot accretion flow model M3f appears to collapse due to thermal instability in a similar fashion as long predicted by analytic accretion theory. This maximum luminosity is lower than the highest observed from hard state BHBs (L/L_Edd ≲ 10⁻¹). However, cold accretion disks appear to be present in at least some luminous BHB hard states (e.g., Miller et al. 2006; Kara et al. 2019) and the hard power law also softens for L/L_Edd ≳ 10⁻² (Skipper & McHardy 2016). Numerical methods that do not rely on sampling individual absorption and scattering events, possibly including moment closure methods (Ryan & Dolence 2020) or implicit Monte Carlo (Roth & Kasen 2015), will be needed to study the equilibrium flow structure and spectrum following collapse to an optically thick accretion disk.

Our models assume that the black hole has (i) saturated magnetic flux and (ii) high spin. Existing measurements suggest that high spin may be common in BHBs (McClintock et al. 2011). If sufficient magnetic flux is available from the binary companion, it may be efficiently advected in the quiescent and hard spectral states (e.g., McKinney et al. 2012) at least for smaller binary systems where a hot accretion flow is expected in the outer region (Begelman & Celotti 2004). Future parameter surveys can explore the dependence of the maximum $\dot{m}$ and radiative efficiency for hot accretion flows on the black hole magnetization and spin. We also neglect radiative cooling for r > 40 r_g , even though the scattering optical depth and Coulomb to viscous heating ratios are found to increase outward. The evolution of the accretion flow and the accumulation of magnetic flux could be significantly altered if the flow becomes optically thick at a larger radius.

The simulations presented here are ≈10–30× more computationally expensive than nonradiative GRMHD, and as such are limited to short durations of ≃0.1s for a 10 M_⊙ black hole. Comparisons to the rich phenomenology of timing features including low-frequency quasi-periodic oscillations and energy-dependent time lags may become possible with longer-duration simulations, especially if they can reach radiative and inflow equilibrium to larger radius.

J.D. thanks B. Ryan and G. Wong for many useful discussions about the ebhlight code. We thank B. Ryan, J. Krolik, E. Quataert, J. Tomsick, and the anonymous referee for helpful comments that improved the paper. This work was supported in part by NASA Astrophysics Theory Program grants NNX16AI40G, NNX17AK55G, and 80NSSC20K0527 and by an Alfred P. Sloan Research Fellowship (J.D.). The calculations presented here were carried out using resources supported by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.