Formation of Magnetically Truncated Accretion Disks in 3D Radiation-transport Two-temperature GRMHD Simulations

Observations spanning the radio to gamma-ray frequency range show that black hole (BH) X-ray binaries (XRBs) cycle through different spectral states of accretion over the course of months to years (e.g., Esin et al. 1997; Remillard & McClintock 2006). BHXRBs spend most of their time in a dim, quiescent state. However, from time to time they go into outburst, during which they transition into orders of magnitude more luminous spectral states, some of which are associated with powerful collimated outflows, or jets. During such spectral state transitions, the spectrum of the source changes from a power-law-like "hard" spectrum to a blackbody-like "soft" spectrum. During the outbursts, the luminosity L can reach a substantial fraction of the Eddington limit, L_Edd, at which the radiation forces become comparable with gravity. In some extreme cases, such as in extragalactic ultraluminous X-ray sources (ULXs) and in the galactic XRB SS433, the inferred bolometric luminosity can even exceed the Eddington limit (e.g., Fabrika 2004; Begelman et al. 2006; Fuchs et al. 2006; Middleton et al. 2021), which can potentially be a sign of super-Eddington accretion (Sądowski & Narayan 2015). In fact, it is conceivable that most BH growth occurs during outbursts when the disk is accreting at or above the Eddington limit (e.g., Volonteri et al. 2007). In particular, the "luminous-hard state" is very interesting because it is associated with the most powerful transient jets in XRBs (e.g., Fender et al. 2004) and Fanaroff–Riley type 2 (FRII; Fanaroff & Riley 1974) jets in luminous active galactic nuclei (AGNs) (e.g., Tchekhovskoy 2015; Davis & Tchekhovskoy 2020).

Accretion disks are generally described by either a geometrically thin, optically thick disk (Novikov & Thorne 1973; Shakura & Sunyaev 1973) or a geometrically thick, optically thin radiatively inefficient accretion flow (RIAF) such as an advection-dominated accretion flow (ADAF; see Narayan & Yi 1994). In a geometrically thin disk, the density and optical depth are high. The balance between the radiative cooling and turbulent dissipation rates determines the geometric thickness, or the scale height, of the disk. Such geometrically thin disk models can explain thermal emission in the high-soft state. In an RIAF, the density is so low that the plasma decouples into a two-temperature fluid, in which the Coulomb collisions between the ions and electrons become too rare to equilibrate the temperatures of the two species. This thermal decoupling, along with the preferential heating of the ions, prevents the disk from radiating a dynamically significant amount of dissipated energy. RIAFs can explain emission in the quiescent and hard spectral states (e.g., Abramowicz & Fragile 2013 and references therein).

When a quiescent BH accretion system goes into an outburst, it transitions into a hybrid state that features both thermal "soft" and nonthermal "hard" components in its spectrum. This suggests the presence of both a standard accretion disk and an RIAF (e.g., Remillard & McClintock 2006). In the literature, "truncated" disk models form a popular paradigm describing such hybrid states (e.g., Esin et al. 1997; Ferreira et al. 2006; Done et al. 2007; Marcel et al. 2018), including their quasi-periodic variability (e.g., Ingram et al. 2009). In these truncated disk models, the outer part of the disk is geometrically thin while the inner part is geometrically thick, hot, and described by an RIAF. We refer to the inner, hot RIAF-like part of the system as the corona. The nature of the corona is actively debated, and its origin and geometry—often represented as a hot "lamppost"—are not definitively known. The corona can Comptonize its own cyclosynchrotron emission in addition to thermal emission from a thin accretion disk (e.g., Wardziński & Zdziarski 2000; Del Santo et al. 2013). It is important to note that the power-law emission attributed to the corona can make up a significant fraction of the total X-ray luminosity in BHXRBs and AGNs and hence in such cases the hot coronal gas is likely to be also dynamically important (e.g., Reynolds & Fabian 1997; Fukumura & Kazanas 2007; Wilkins & Fabian 2012, and references therein).

A major point of contention concerns the behavior of the truncation radius, ${r}_{\mathrm{tr}}$, which marks the transition from the outer geometrically thin, radiatively efficient disk to the inner geometrically thick, radiatively inefficient corona. It is not clear if the softening of the spectrum during a state transition is primarily caused by the truncation radius moving inwards (Fender et al. 2004) or the corona shrinking in size (e.g., Kara et al. 2019), and what roles Compton-upscattered X-ray emission from the base of the jet (e.g., Markoff et al. 2005) and/or relativistic winds (Beloborodov 1999) play. In addition, the presence of broadened iron lines (e.g., Reis et al. 2010), which are typically associated with a reservoir of cold optically thick gas near the innermost stable circular orbit (ISCO), in the low-hard state of several BHXRB systems is inconsistent with the standard picture of a truncated disk, which does not predict the presence of any cold gas within the truncation radius (e.g., Done & Gierliński 2006; Done et al. 2007).

Simulating truncated accretion disks remains extremely challenging because both radiation and thermal decoupling between ions and electrons is expected to play a prominent role. Current MHD simulations only partly addressed disk truncation by simplifying the (radiative) physics with an ad hoc cooling function, which led to a smooth transition (Hogg & Reynolds 2017, 2018). Here, we use the state-of-the-art radiation-transport two-temperature GRMHD simulations to demonstrate that truncated accretion disks naturally form in the presence of a large-scale net poloidal magnetic flux. In Sections 2 and 3 we describe our GRMHD code H-AMR and the initial conditions, before presenting our results and concluding in Sections 4 and 5.

We use for our simulations the graphical processing unit (GPU) accelerated GRMHD code H-AMR (Liska et al. 2018, 2019a; see Porth et al. 2019 for code comparison). H-AMR is triple-level parallelized using a hybrid CUDA-OpenMP-MPI framework that solves the equations of GRMHD in the Kerr–Schild foliation. Magnetic fields are evolved using a staggered grid approach based on Gardiner & Stone (2005). The spatial reconstruction of primitive variables from cell centers to cell faces is done using a third-order-accurate PPM method (Colella & Woodward 1984), which leads to second-order overall spatial accuracy (cross-dimensional terms are not taken into account). Integration in time is performed using the second-order IMEX2A scheme described in McKinney et al. (2013). H-AMR features a flexible adaptive mesh refinement (AMR) framework that allows for refinement both in space and time. Namely, our local adaptive time-stepping approach allows different blocks to have different time steps even at the same refinement level (Liska et al. 2019a).

We have recently implemented on-the-fly radiation transport into H-AMR using an M1 (two-moment) closure, following the implementations in Koral (Sa̧dowski et al. 2013, 2017) and HARMRAD (McKinney et al. 2013). We include the energy-averaged bound–free, free–free, and cyclosynchrotron absorption (κ_abs) and emission (κ_em) opacities in addition to the electron-scattering (κ_es) opacity as given in McKinney et al. (2017). We use gray opacities energy-averaged over a diluted blackbody spectrum (McKinney et al. 2017) to provide as accurate as possible opacities within a single-energy M1 framework. In addition, we assume that only electrons contribute and set κ_ion = 0. This is valid because, in high-density regions where bound–free opacities (κ_bf ∝ ρ) become dominant, Coulomb collisions typically equilibrate the plasma to a single temperature.

The M1 closure works well in the optically thick disk body and outside of the thin disk where the radiation moves predominantly in one direction (perpendicular to the disk surface). Because this is not always valid in the corona, we further discuss the physical implications of the M1 approximation on our results in Section 5. To approximate Compton cooling of the corona, we implemented both blackbody and photon-number-conserving Comptonization in H-AMR (Sa̧dowski & Narayan 2015). By evolving the photon number (n_rad), a more accurate radiation temperature can be calculated in the optically thin corona (T_rad ∼ (E_rad/(2.7k_b n_rad) instead of ${T}_{\mathrm{rad}}={({E}_{\mathrm{rad}}/a)}^{1/4}$). This increases the accuracy of the temperature-dependent emission, absorption, and Compton scattering coefficients in the plasma and might play an important effect when the radiation spectrum deviates from a blackbody such as in an optically thin corona.

To allow for the development of a two-phase medium, we have implemented two-temperature thermodynamics in H-AMR and combined it with our M1 radiation scheme. In H-AMR, we evolve the ion and electron entropy tracer, κ = p_e,i /ρ^Γ−1, where p_e,i is the electron and ion pressure, respectively, and ρ is the fluid-frame rest-mass density. In this work, we adopt Γ = 5/3 for both electrons and ions: We have verified that this is a safe choice by carrying out simulations with a variable polytropic index for the electrons. The total dissipation is calculated as the difference between the internal energies computed from the total energy equation and from the entropy tracers (Ressler et al. 2015). For the purposes of partitioning the dissipated energy between the ions and electrons, we assume here that the dissipation ultimately occurs via magnetic reconnection in unresolved current sheets. The dissipated energy is divided between the ions and electrons using an analytical prescription derived from particle-in-cell simulations (Rowan et al. 2017). We account for the energy exchange via Coulomb collisions between ions and electrons through an implicit source term (Sa̧dowski et al. 2017). Due to the low radiative efficiency of ions compared to electrons, when a two-temperature fluid forms, we assume that radiative cooling reduces the electron entropy and keeps the ion entropy constant.

H-AMR allows broad flexibility in the choice of a coordinate system. Here, we adopt spherical polar coordinates, r, θ, ϕ, and choose a uniform grid in the $\mathrm{log}r$, θ, and ϕ variables. We use a base resolution of 1020 × 432 × 288 cells in the three dimensions, respectively. We then use AMR to increase this resolution in the regions of interest, as we describe below. We place the inner radial boundary inside of the event horizon (with ≥5 radial cells inside of the event horizon at the base level: This ensures that the inner radial grid boundary is properly causally disconnected from the BH exterior). We place the outer radial boundary at R_out = 10⁴ r_g . To resolve the structure of the thin equatorial disk, we quadruple the resolution for 5r_g ≲ r ≲ 120r_g within 75 of the midplane to an effective resolution of 4080 × 1728 × 1152 cells by adding two layers of static mesh refinement. In addition, we use four levels of local adaptive time stepping, which adaptively sets the time step based on the local Courant condition (Liska et al. 2019a). This increases the time step in the outer disk, which speeds up and improves the accuracy of the simulations (Chatterjee et al. 2019). To prevent cell squeezing near the poles (θ = 0, π), we progressively reduce the ϕ resolution to N_ϕ = 18 within 30° from each pole (we leave the resolutions in the other two dimensions unaffected). We use outflow boundary conditions in r, transmissive boundary conditions in θ, and periodic boundary conditions in ϕ directions. The grid features in total ∼0.45 × 10⁹ cells and achieves a computational speed of effectively 1.5 × 10⁷ zone-cycles/s/GPU (taking into account the speedup from the local adaptive time stepping). This is an order of magnitude less compared to the nonradiative version of H-AMR. The grid is visualized in Figure 2.

To study the effect of large-scale magnetic flux on the structure of the inner disk, we consider two models. We initialize the model RADPOL with a purely poloidal magnetic field and the model RADTOR with a purely toroidal magnetic field. Both models feature a rapidly spinning BH of dimensionless spin a = 0.9375 and use Γ = 5/3 for ions and electrons. Because ideal GRMHD is unable to describe the physical processes responsible for injecting gas in the jet funnel, we use drift-frame density floors (Ressler et al. 2017) to ensure that ρ c² ≥ p_mag/12.5 in the jets.

Model RADPOL starts with a thick equatorial torus in hydrostatic equilibrium with the inner radius at r_in = 12.5r_g , density maximum at ${r}_{\max }=25{r}_{g}$, and outer radius r_out = 200r_g (Fishbone & Moncrief 1976). The magnetic field has the vector potential A_ϕ ∝ (ρ − 0.05)² r². We normalize this magnetic field such that $\beta =\max {p}_{\mathrm{gas}}/\max {p}_{\mathrm{mag}}=30$, where p_gas and p_mag are the gas and magnetic pressures, and "max" denotes the maximum taken over the torus. This setup is intended to represent the hard intermediate state, during which a hard-state thick torus collapses into a thin accretion disk due to radiative cooling (e.g., Esin et al. 1997). This radiative collapse is modeled by cooling the torus on the orbital timescale using a cooling function (see Section 4). This cooling process is more thoroughly described in Liska et al. (2019b), who used the same setup but ran the simulation at a lower resolution and for a shorter duration.

Model RADTOR starts with a thin equatorial accretion disk on Keplerian orbits in approximate hydrostatic equilibrium whose density distribution, $\rho (r,z)\propto {r}^{-1}\exp (-{z}^{2}/2{h}^{2})$, extends from an inner radius, r_in = 6r_g , to the outer radius, r_out = 76r_g at the constant h/r = 0.03. The magnetic field is described by a vector potential A_θ ∝ (ρ − 0.0005)r², which is normalized to give an approximately uniform β ∼ 7. This setup is intended to represent the high-soft state, where neither large-scale poloidal magnetic flux nor jets are present, presumably because all net poloidal magnetic flux has diffused out after the disk has settled into a geometrically thin disk (e.g., Begelman & Armitage 2014).

We evolve both models in two stages. During the first stage, we use neither radiation transport nor two-temperature thermodynamics. Instead, we use a cooling function (Noble et al. 2009) to cool the disk to a target thermal scale height of (h/r)_init = 0.03: This cooling process takes ∼ 10³ r_g /c for RADPOL and mimics the catastrophic cooling of a BHXRB disk undergoing an outburst as described in Liska et al. (2019b). We initialize RADTOR at the target temperature profile corresponding to the same disk scale height, (h/r)_init = 0.03. We evolve RADPOL for t ∼ 188, 840r_g /c and RADTOR for t ∼ 136, 075r_g /c, which gives both models sufficient time to reach an approximately constant accretion rate on the BH. This allows us to select a well-defined accretion rate (with respect to the Eddington rate) before we turn on radiation.

The left panel in Figure 1 shows that RADPOL accumulates large-scale poloidal magnetic flux in the inner disk: This is similar to the model described in Liska et al. (2019b), which is essentially the same as RADPOL, apart from a small disk tilt, numerical grid, and duration. Model RADTOR does not generate nor advect any significant large-scale poloidal magnetic flux as was suggested in previous work considering geometrically thick disks (Liska et al. 2020). The reason for this is unclear, but it might involve the inability of geometrically thin disks to generate and/or advect magnetic flux inwards (e.g., Lubow et al. 1994). It is also possible that the simulation was not run long enough to capture this effect for geometrically thin disks, which have a two orders of magnitude longer viscous time than the disk presented in Liska et al. (2020). The flux accumulation in RADPOL is accompanied by a decrease in density at r ≲ 20r_g and the development of a magnetically arrested disk (MAD; Tchekhovskoy et al. 2011): The accretion in the inner disk proceeds through magnetic Rayleigh–Taylor instabilities (RTI), as is expected in a MAD (e.g., Tchekhovskoy et al. 2011; McKinney et al. 2012). Neither of these two effects is observed for RADTOR (right panel of Figure 1) and highlights the importance of the net poloidal magnetic flux for driving the disk truncation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the second stage, we set M_BH = 10M_⊙, renormalize the density such that $\langle \dot{M}\rangle \sim 0.35{\dot{M}}_{\mathrm{Edd}}$, and set T_rad = T_e = T_i , where we define ${\dot{M}}_{\mathrm{Edd}}={L}_{\mathrm{Edd}}/\eta$ with η = η_NT = 0.178 and L_Edd = 4π GM_BH/κ_es. We assume solar abundances of hydrogen, helium, and metals (X = 0.7, Y = 0.28, Z = 0.02). We then run both simulations for another t ∼ 1.35 × 10⁴ r_g /c in full radiation-transport two-temperature GRMHD: for time interval t ∈ [188, 840, 203, 360]r_g /c in RADPOL and t ∈ [136, 076, 148997]r_g /c in RADTOR. We use blackbody Comptonization for both models, except during t ∈ [191, 840, 193, 650]r_g /c in RADPOL, where we use photon-number-conserving Comptonization. In addition, we use for RADPOL during t ∈ [193, 650, 196, 550]r_g /c a lower-density floor of ρ c² ≥ p_mag/25. Additionally, for the time interval of t ∈ [188, 840, 192630]r_g /c in RADPOL, we limit the ratio of the electron to ion temperatures to T_e /T_i ≳ 10⁻² instead of the usual T_e /T_i ≳ 10⁻⁴. These tests validate that the results presented in this article are not sensitive to the Comptonization routine employed or the value of the density floors as can be witnessed by the animations on our YouTube channel and in Figure 10. However, setting the temperature floor to T_e /T_i ≳ 0.01 artificially increased the temperature of the electrons in current sheets near the midplane within r ≲ 5r_g of the BH.

The top row in Figure 3 (see also our YouTube channel for animations of the figure) shows that in RADPOL the inner disk decouples into a hot, optically thin, two-temperature plasma and can thus be associated with a corona. Because the physical processes (e.g., Levinson & Cerutti 2018) responsible for mass loading of the jet, such as pair creation, cannot be modeled in GRMHD simulations, which artificially inject gas in the jet for the scheme to remain numerically stable, the jet thermodynamics fall beyond the scope of this work. Electrons in the corona reach temperatures of T_e ≳ 5 × 10⁸ K while the ions reach temperatures of T_i ≳ 10¹⁰ K irrespective of the Comptonization routine employed. Dissipation, most likely caused by magnetic reconnection in current sheets near the midplane (e.g., Ripperda et al. 2022), leads to localized heating of the plasma to temperatures of T_i ∼ 10¹² K and T_e ∼ 10⁹ K, which suggest that electrons radiate their heat locally (see also Beloborodov 2017). A significant part of this heated plasma flows out as wind along the magnetic field lines (Figure 7). On the other hand, pockets of cold gas survive and fall into the BH predominantly along the equatorial plane. As explained in Section 5 this has interesting consequences for reflection modeling.

Zoom In Zoom Out Reset image size

Figure 4 shows the cooling rate of the plasma due to emission (${\lambda }_{\mathrm{em}}=\rho {\kappa }_{\mathrm{em}}{a}_{\mathrm{rad}}{T}_{e}^{4}$) and inverse-Compton scattering (λ_sc = ρ κ_es E_rad 4k_b T_e /m_e c²), in addition to the heating rate through absorption (λ_abs = ρ κ_abs E_rad). As can be seen, the cooling is predominantly concentrated near the equatorial plane which has an optical depth of order unity (τ_es ∼ 1). Because the Compton scattering rate exceeds the emission rate, it is conceivable that most radiation actually gets Compton upscattered before leaving the corona. We estimate that the Compton-y parameter (y = τ_es 4k_b T_e /m_e c²), which measures the average fractional energy gain of photons due to Comptonization (e.g., Rybicki & Lightman 1986), can easily exceed y ≳ 1 in the T_e ≳ 5 × 10⁸ K plasma near the midplane of the corona. Pinpointing the exact value of the Compton y parameter is nontrivial because absorption will play an important role in some optically thick patches near the midplane. A more detailed calculation that takes into account emission, absorption, and scattering in the plasma along photon geodesics will be required to further constrain the Compton-y parameter. Nevertheless, this work suggests that the corona will upscatter and significantly harden thermal radiation from the disk and cyclosynchrotron emission from the corona (see also Dexter et al. 2021; Scepi et al. 2021). In contrast, the bottom row in Figure 3 shows that in RADTOR, the ions and electrons are strongly coupled in the cold disk, which remains optically thick to both scattering and absorption. Because this disk extends all the way to the event horizon, we expect a strong thermal emission component in the emergent spectrum. Compton upscattering of disk photons by the winds sandwiching the inner disk to temperatures necessary to significantly harden the spectrum is not expected because these winds are predominantly cold (T_e ≲ 5 × 10⁷ K). The hot pockets of gas with temperatures T_e ≳ 10⁸ K within r ≲ 5r_g in Figure 3 are relatively rare and seem to be more prominent at later times when the disk becomes numerically underresolved due to thermal collapse (Appendix A.2).

Zoom In Zoom Out Reset image size

We show the accretion rate $\dot{M}=-{\int }_{0}^{2\pi }{\int }_{0}^{\pi }\rho {u}^{r}\sqrt{-g}d\theta d\phi$, with the integral evaluated at r = 5r_g to avoid contamination by the density floors and dimensionless BH poloidal magnetic flux ${\phi }_{\mathrm{BH}}=0.5{\int }_{0}^{2\pi }{\int }_{0}^{\pi }| {B}^{r}| \sqrt{-g}d\theta d\phi /\langle \dot{M}{\rangle }^{1/2}$, with the integral evaluated at the event horizon, r = r_H ≃ 1.3r_g , in the upper and middle panels of Figure 5. In the RADPOL model, the accretion disk transports so much of the initial poloidal magnetic flux to the BH that it becomes as strong as gravity, obstructs the accretion, reaches a saturated value, and leads to an MAD. The normalized magnetic flux ϕ_BH ∼ 33 is about two-thirds of that in thick accretion disks (Tchekhovskoy et al. 2011) and consistent with the MAD saturation value in thin accretion disks (S. Cheng & M. Liska 2022, in preparation). The lower panel of Figure 5 shows that RADPOL launches powerful jets and winds. Namely, we show the jet (η_jet) and wind (η_wind) energy efficiencies, which we obtain by normalizing the electromagnetic plus fluid energy fluxes at r = 5r_g in the jets and winds by the mass accretion rate (see Liska et al. 2019b). We define the jet-wind boundary by the p_mag = ρ c² criterion. In RADPOL, we find η_wind ∼ 15% and η_jet ∼ 50%. In RADTOR, we find η_wind ∼ 20% and η_jet ≲ 1%. In both models, the associated winds are significantly hotter than the disk. We compute the radiative efficiency from the normalized radiation energy flux measured at r = 100r_g and find η_rad ∼ 20% in both models. This is remarkably close to the Novikov & Thorne (1973) value η_NT ∼ 17.8%.

Zoom In Zoom Out Reset image size

As Figure 3 shows, in RADPOL the excess magnetic flux remains in the inner disk, $r\lesssim {r}_{\mathrm{tr}}\approx 20{r}_{g}$. The left column in Figure 6 shows that this leads to the formation of a magnetic-pressure-supported low-density inner disk coupled to a radiation-pressure-supported high-density outer disk. The transition between these two regimes occurs at a well-defined "magnetic truncation" radius, ${r}_{\mathrm{tr}}$, which is set by the distance out to which the inner disk is flooded by the poloidal magnetic flux. Similar to the single-temperature simulation stage (Figure 1), here the accretion in the inner disk also proceeds through magnetic RTI and is the radiative analog of the MAD (see also Morales Teixeira et al. 2018). On the other hand, the disk in RADTOR remains radiation pressure dominated and does not contain any large-scale poloidal magnetic flux. We define the effective viscosity as ${\alpha }_{\mathrm{eff}}=-\langle {v}_{r}{v}_{\phi }{\rangle }_{\rho }/\langle (h/r{)}_{\mathrm{init}}^{2}{v}_{\phi }^{2}{\rangle }_{\rho }$, where 〈...〉_ρ is the density-weighted angular average, and v_r and v_ϕ are the physical radial and azimuthal velocity components, respectively. Figures 6(g), (h) shows that α_eff ∼ 10¹ in RADPOL and α_eff ∼ 10⁻² in RADTOR. This illustrates that the presence of large-scale magnetic flux, as in RADPOL, leads to more rapid accretion of gas (e.g., Ferreira et al. 2006). Please note that even when α_eff ∼ 100, the ratio between the radial and azimuthal velocities becomes v_r /v_ϕ = (h/r)² α ≲ 0.1, suggesting the infalling gas has a significant toroidal velocity component capable of producing a broadened iron line.

Zoom In Zoom Out Reset image size

In both models, the scale height of the radiation-pressure-supported part of the disk shrinks to a density-averaged scale height h/r = 〈∣θ − π/2∣〉_ρ ∼ 0.015–0.02 (Figure 6 panels g and h and YouTube channel) at which point (when the disk scale height or MRI wavelength becomes resolved by less than ∼8 cells beyond r ≳ 5r_g ) we terminate the simulation. This runaway cooling effect suggests that both disks are subject to viscous and/or thermal instabilities (Lightman & Eardley 1974; Shakura & Sunyaev 1976). However, in RADPOL, magnetic pressure stabilizes the inner disk against thermal collapse (see also Sa̧dowski 2016; Jiang et al. 2019), which settles into a scale height of h/r ∼ 0.05–0.2. Please note that the outer disk in both RADPOL and RADTOR is radiation pressure dominated and thus remains thermally unstable. Future simulations with higher resolutions and longer runtimes will investigate whether magnetic pressure or other physical processes eventually stabilize the entire accretion disk against runaway collapse (e.g., Begelman & Pringle 2007). This was supported by radiative GRMHD simulations of thin accretion disks seeded with equipartition strength magnetic fields (e.g., Sa̧dowski 2016; Lančová et al. 2019).

In this work, we have presented the first radiation-transport two-temperature GRMHD simulations of luminous sub-Eddington accretion disks that accrete close to the Eddington limit (L/L_edd ≳ 0.01). Previous work in GRMHD only addressed highly sub-Eddington disks with $\lambda =\langle \dot{M}\rangle /{\dot{M}}_{\mathrm{Edd}}\ll 0.01$ (e.g., Chael et al. 2017; Sa̧dowski et al. 2017; Chael et al. 2019; Dexter et al. 2021). More specifically, the simulations in this article accrete at the Eddington ratio $\lambda =\langle \dot{M}\rangle /{\dot{M}}_{\mathrm{Edd}}\sim 0.35$. One of the simulations, model RADTOR, had no net poloidal magnetic flux, while the other one, model RADPOL, became saturated with poloidal magnetic flux and entered an MAD state. These models were chosen to demonstrate the importance of magnetic flux saturation on the (thermo)dynamics of the plasma and do not attempt to describe a possible transition between RADPOL and RADTOR and/or to rule out other physical mechanisms that can lead to spectral hardening or disk truncation (e.g., Meyer & Meyer-Hofmeister 1994). We demonstrated that the saturation of magnetic flux in the inner disk of model RADPOL creates a sharp transition between a radiation-pressure-supported outer disk and a magnetic-pressure-supported inner corona. The sharp transition between the two occurs at a magnetic truncation radius of ${r}_{\mathrm{tr}}\sim 20{r}_{g}$. Powerful jets, winds, and radiative outflows with a combined efficiency of η_tot ≳ 90% emerge in the presence of this poloidal magnetic flux. All of these characteristics are consistent with accretion in the hard (-intermediate state) of BHXRBs.

The corona is best described by an MAD where radial magnetic pressure gradients prevent ordered quasi-axisymmetric accretion of gas (e.g., Igumenshchev et al. 2003; Narayan et al. 2003; Igumenshchev 2008; Tchekhovskoy et al. 2011; McKinney et al. 2012; Begelman et al. 2022). Accretion of gas in the corona proceeds through magnetic RTI instabilities originating at the truncation radius. While a significant fraction of the infalling gas gets heated by magnetic reconnection to T_i ≳ 10¹⁰ K and T_e ≳ 5 × 10⁸ K, we also observe cold, optically thick, clumps of gas falling into the BH (Figure 3). This could potentially explain broad iron line emission observed in the hard spectral states of XRBs and AGNs (e.g., Reis et al. 2010) without the need for the corona to recondense into a thin accretion disk near the BH (e.g., Meyer-Hofmeister & Meyer 2011).

Our model differs substantially from most truncated accretion disk models in the literature (e.g., Esin et al. 1997; Ferreira et al. 2006; Begelman & Armitage 2014). In such models, heat conducted from the corona into the upper layers of the disk leads to the disk's gradual evaporation into a corona (e.g., Meyer & Meyer-Hofmeister 1994; Liu & Osher 1998; Qian et al. 2007). However, in magnetically truncated disks, the transition from a disk into a corona occurs at a well-defined radius, which is determined by the location of the magnetospheric radius (e.g., Igumenshchev 2009). Because accretion within the magnetospheric radius proceeds through RTI modes, we expect that the power spectrum of the Comptonized emission from the corona will differ substantially from the outer disk's blackbody emission. In future work, we will address how the cold-phase gas transitions into hot-phase gas in the corona. At the moment of writing, we suspect that magnetic reconnection in equatorial current sheets heats a small fraction of the cold-phase gas to extremely high temperatures while leaving most of the cold-phase gas unaffected. The hot-phase case subsequently flows out along magnetic field lines (Appendix A.1). This is a conceivable scenario in a magnetic-pressure-supported plasma, but will require further analysis (M. Liska et al. 2022a, in preparation).

Spectral state transitions (e.g., Fender et al. 2004) can potentially be explained by the outward diffusion of magnetic flux (e.g., Lubow et al. 1994; Begelman & Armitage 2014). From an observational perspective, the truncation radius needs to move inwards to explain the softening of the spectrum during a hard-to-soft state transition (e.g., Esin et al. 1997) while the magnetic flux needs to diffuse out to explain the lack of any jets in the high-soft state (e.g., Fender et al. 2004). This has never been observed in GRMHD simulations because magnetic flux diffusion is expected to occur on much longer timescales than can be captured by state-of-the-art GRMHD models. Nevertheless, it is conceivable that outward diffusion of the magnetic flux in our model will progressively reduce the magnetic truncation radius, which will cause the spectrum to become softer. However, even when the entire disk leaves the MAD state, residual large-scale poloidal magnetic flux might still be present and launch powerful winds from the inner disk. Such winds could be significantly hotter compared to the winds in RADTOR and could potentially explain the hard emission tail observed in the soft (-intermediate) states of BHXRBs. In addition, other processes such as shocks (Musoke et al. 2022; M. Liska et al. 2022b, in preparation) and emission from jets (e.g., Markoff et al. 2005) can lead to significant hardening of the spectrum.

Interestingly, Kinch et al. (2021) used an iterative postprocessing approach (Kinch et al. 2019) to generate spectra of thin accretion disks, which were simulated in GRMHD using a cooling function (Noble et al. 2009, 2011). In these simulations, the complex radiative processes responsible for cooling the disk and corona (defined as gas outside of the τ_es ≲ 1 surface) were approximated with cooling functions (Noble et al. 2009; Kinch et al. 2020). The temperature of the radiation-emitting electrons in the corona was determined by matching the cooling rate in the corona to the Compton emission rate of the electrons. Making a direct comparison between Kinch et al. (2021) and RADPOL/RADTOR simulations is difficult because both the initial configuration and numerical methods differ substantially. The thin disk in Kinch et al. (2021) is, for example, much thicker (h/r ∼ 0.05–0.06) and contains a significant amount of poloidal magnetic flux, but stays below the MAD saturation value, suggesting neither RADPOL or RADTOR can describe it. In addition, our radiative GRMHD simulations model the radiative processes in the disk and corona such as emission, absorption, and scattering self-consistently. We also account for thermal decoupling between ions and electrons in the corona and include feedback between radiation and gas. For example, in our radiative models, the local dissipation rate is not guaranteed to match the local emission rate, which can lead to runaway thermal collapse of the accretion disk (Section 4).

By ray-tracing our simulation data (e.g., Krawczynski 2021) we will be able to place unique constraints on the relative contributions of the cold-phase and hot-phase gas. However, generating self-consistent spectra will face several challenges. For example, magnetic reconnection in collisionless plasma leads to the formation of nonthermal particle populations (e.g., Sironi & Spitkovsky 2014; Beloborodov 2017; Sironi & Beloborodov 2020; Sridhar et al. 2021). Because ideal GRMHD is unable to model the full energy distribution of such particles, deviation from the thermal spectrum is a free parameter that should be quantified by future particle-in-cell (PIC) simulations. In addition, it is likely that there are some inaccuracies in the electron temperatures presented in this work. These stem from the intrinsic challenges associated with quantifying the dissipation rate in a highly magnetized plasma, and from the inherent limitations related to our two-moment (M1) radiation scheme, which treats radiation as a highly collisional fluid. This works best in regions of high optical depth where emission and scattering of photons are localized. However, in the optically thin corona, there are multiple radiation fields that M1 cannot capture. These include thermal radiation from the cold-phase gas and disk, Compton-scattered radiation from the hot-phase gas in the corona, and reflected coronal radiation from the cold-phase gas and disk. This can lead to irradiation of the dense gas clumps in the corona by each other and the hot-phase gas, which can potentially disperse the plasma to a less clumpy state than predicted by M1. In other words, future work will be necessary to determine if the clumpy structure in the corona survives when properly accounting for all radiation fields and if the irradiated clumps will indeed be able to produce a relativistic broadened iron line (which will require an elaborate calculation to determine the ionization state of the plasma).

We thank Michiel van der Klis, Adam Ingram, Sera Markoff, Ramesh Narayan, and Eliot Quataert for insightful discussions. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program under award PHY129. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC05-00OR22725. This research also used HPC and visualization resources provided by the Texas Advanced Computing Center (TACC) at The University of Texas at Austin, which contributed to our results via the LRAC allocation AST20011 and through Texascale Days (http://www.tacc.utexas.edu). We acknowledge PRACE for awarding us access to JUWELS Booster at GCS@JSC, Germany. M.L. was supported by the John Harvard Distinguished Science Fellowship, G.M. is supported by a Netherlands Research School for Astronomy (NOVA), Virtual Institute of Accretion (VIA) postdoctoral fellowship, and A.T. by the National Science Foundation grants AST-2009884, AST-1815304, and AST-1911080.

A.1. Formation of Outflows

Both RADPOL and RADTOR form radiative and mechanically driven outflows. In Figure 7 we use streamlines to illustrate the direction of the radiation velocity (red) and gas velocity (black). Both of these outflows move out radially or vertically from the disk's midplane. In future work, we plan to address the mass and energy ejection rates of these outflows and if these outflows indeed become gravitational unbound. Note that radiation is treated as a single fluid in the M1 approximation and thus only the net radiative energy flow is captured. In reality, the radiation field would consist of a blackbody component originating in the disk and Compton-scattered radiation originating from the corona, which can be reflected by the disk.

Zoom In Zoom Out Reset image size

A.2. Convergence

To quantify the convergence of our models, we plot in Figure 8 the radial dependence of the mass accretion rate and MRI quality factors. In model RADPOL, the mass accretion rate reaches a steady value within r ≲ 25r_g . This suggests RADPOL achieved inflow equilibrium in the MAD region of the disk and corona. In RADTOR, no inflow equilibrium can be achieved because the disk collapses into an infinitely thin slab due to radiative cooling. The MRI quality factors (Q_r , Q_θ , Q_ϕ ) are defined as the number of cells per MRI wavelength (λ_MRI) weighted by b² ρ (${Q}_{r,\theta .\phi }=| {\lambda }_{\mathrm{MRI}}{| }_{{b}^{2}\rho }/{N}_{r,\theta ,\phi }$). Resolving the MRI wavelength by at least 10 cells, and ideally by more than 20 cells (e.g., Shiokawa et al. 2012; Porth et al. 2019), is necessary to capture the growth of the longest-wavelength MRI modes in the θ and ϕ directions. Both RADPOL and RADTOR have Q_ϕ ≳ 20 throughout the domain. However, in RADTOR, Q_r and Q_θ reach a marginal Q_r ∼ Q_θ ∼ 6–10 for r ≲ 10. Because the initial setup only contains a toroidal magnetic field, it is not unexpected that the poloidal components of the magnetic field become challenging to resolve, especially after the inner disk has undergone thermal collapse. We suspect that this can potentially heat up the plasma above the photosphere in RADTOR (Figure 9).

Zoom In Zoom Out Reset image size

Download figure:

Video Standard image High-resolution image

A.3. Effect of Comptonization and Density Floors