Radiation RMHD accretion flows around spinning AGNs: a comparative study of MAD and SANE state

Here, we present ideal RMHD governing equations for the interaction between matter and electromagnetic (EM) fields in cylindrical coordinates $(r,\phi,z)$. The governing equations are as follows (Melon Fuksman & Mignone, 2019)

where $\rho$, $\bm{v}$, $\bm{m}$, $h$, $\gamma$ and $\bm{B}$ are the mass density, fluid velocity, momentum density, specific enthalpy, Lorentz factor and mean magnetic field, respectively. $\bm{E}=-\bm{v}\times\bm{B}$ is the electric field. Moreover, $E_{\rm rad}$, $\bm{F}_{\rm rad}$ and $P_{\rm rad}$ are the radiation energy density, the radiation flux, and the pressure tensor as moments of the radiation field, respectively. Here, the total pressure, momentum density, and energy density of matter are given by

Also, here the radiation-matter interaction terms ($G^{0},\bm{G}$) is given by

where all the fields are measured in the fluid’s comoving frame. Here $\kappa$ and $\sigma$ are the frequency-averaged absorption and scattering opacity, respectively. In this work, we adopt free-free absorption opacity as $\kappa=1.7\times 10^{-25}m_{p}^{-2}\rho T^{-7/2}~{}{\rm cm}^{2}~{}{\rm g}^{-1}$ and scattering opacity as $\sigma=0.4~{}{\rm cm}^{2}~{}{\rm g}^{-1}$, respectively (Igarashi et al., 2020). Here, $T$ is the temperature of the fluid and $a_{R}$ is the radiation density constant. In this work, the radiation transfer equations are solved based on the gray approximation and imposing the M1 closure (Melon Fuksman & Mignone, 2019).

In order to close the system of equations (1-7), it is required additional sets of relations. One of them is the equation of state (EoS) which provides closure between thermodynamical quantities. In this work, we consider constant-$\Gamma$ EoS, and is given as

where, $\Theta=\frac{P_{\rm gas}}{\rho}$ and $\Gamma$ is the adiabatic index. Further, another closure relation is needed for the radiation fields which relates to pressure tensor $P_{\rm rad}^{ij}$, $E_{\rm rad}$ and $\bm{F}_{\rm rad}$. The closure relation is as follows (Melon Fuksman & Mignone, 2019)

where, $\bm{n}=\bm{F}_{\rm rad}/|{\bm{F}_{\rm rad}}|$ and $f=|{\bm{F}_{\rm rad}}|/E_{\rm rad}$. Here $\delta^{ij}$ is the Kronecker delta.

In this work, the gravitational effect around a spinning black hole is modeled using effective Kerr potential proposed by Dihingia et al. (2018). The effective Kerr potential is as follows

where $R=\sqrt{r^{2}+z^{2}}$ = spherical radial distance, $\Delta=a_{k}^{2}+R^{2}-2R$, $\Sigma=\frac{a_{k}^{2}z^{2}}{R^{2}}+R^{2}$ and $A=\left(a_{k}^{2}+R^{2}\right)^{2}-\frac{a_{k}^{2}r^{2}\Delta}{R^{2}}$. Here, $\lambda$ and $a_{k}$ are specific flow angular momentum and black hole spin, respectively. The Keplerian angular momentum can be obtained from equation (17) in the equatorial plane $(z\rightarrow 0)$ as $\lambda_{K}=\sqrt{r^{3}\frac{\partial\Phi^{\rm eff}}{\partial r}|_{\lambda\rightarrow 0}}$. The angular frequency is obtained as $\Omega=\lambda/r^{2}$. The effective Kerr potential is implemented in PLUTO code in a similar way as Aktar et al. (2024). It is to be mentioned that Dihingia et al. (2018) examined various comparative studies between GR and effective Kerr potential based on analytical studies. They showed an excellent close agreement of the transonic properties adopting effective Kerr potential in the semi-relativistic regime and GR around the Kerr black hole. Therefore, our numerical approach qualitatively retains the essential space-time features around black holes by implementing effective Kerr potential. Moreover, our simulation model could achieve higher spatial resolutions than GRMHD simulations for a given computing resource.

The initial equilibrium torus is constructed by adopting the Newtonian analog of relativistic tori as prescribed by Abramowicz et al. (1978). In this work, we adopt the same formalism to set up the initial equilibrium torus as described by Aktar et al. (2024). The density distribution of the torus is obtained by considering constant angular momentum flow (Matsumoto et al., 1996; Hawley, 2000; Kuwabara et al., 2005; Aktar et al., 2024)

where ‘$\mathcal{C}$’ is the constant and it is calculated considering zero-gas pressure $(P_{\rm gas}\rightarrow 0)$ at $r=r_{\rm min}$ at the equatorial plane. Here, $r_{\rm min}$ is the inner edge of the torus. The density distribution inside the torus is obtained using the adiabatic equation of state $P_{\rm gas}=K\rho^{\Gamma}$ as

where $K$ is determined considering the density maximum condition $(\rho_{\rm max})$ at $r=r_{\rm max}$ at the equatorial plane and is given by

On the other hand, the density distribution outside the torus can be obtained by assuming hydrostatic equilibrium condition (Matsumoto et al., 1996; Kuwabara et al., 2005; Aktar et al., 2024). Therefore, the density outside the torus is assumed to be the isothermal, nonrotating, high-temperature halo, and is given by (Manmoto et al., 1997; Kuwabara et al., 2005; Igarashi et al., 2020)

where, $\cal{H}$ = $1/c_{s}^{2}$ and $c_{s}=\Gamma P_{\rm gas}/\rho$ is the sound speed. Here, we consider $\eta=10^{-4}$ and $\cal{H}$=2 for the purpose of representation (Aktar et al., 2024).

In this paper, we initialize a poloidal magnetic field by applying a purely toroidal component of vector potential (Hawley & Krolik, 2002). The expression of the toroidal component of vector potential is as follows (Hawley & Krolik, 2002; Aktar et al., 2024)

where, $B_{0}$ is the normalized initial magnetic field strength and $\rho_{\rm min}$ is the minimum density in the torus. Here the initial magnetic field strength is parameterized by the ratio of the gas pressure to the magnetic pressure known as plasma-$\beta$ parameter $\beta_{0}=\frac{2P_{\rm gas}}{B_{0}^{2}}$. Also, we define the magnetization parameter as $\sigma_{\rm M}=\frac{B^{2}}{\rho}$ (Dihingia et al., 2021, 2022; Dhang et al., 2023; Curd & Narayan, 2023; Aktar et al., 2024). The magnetization parameter essentially represents the ratio of magnetic energy to the rest mass energy in the flow.

In this work, to simulate radiation relativistic MHD accretion flows, we adopt PLUTO simulation code (Melon Fuksman & Mignone, 2019). Here, we consider the computational domain for our simulation model in the radial and vertical direction as $1.5r_{g}\leq r\leq 200r_{g}$ and $-100r_{g}\leq z\leq 100r_{g}$, respectively. We adopt a comparatively high resolution of grid sizes as $(n_{r},n_{z})=(920,920)$ and uniform grid spacing in both radial and azimuthal directions. Here, we use the HLL Riemann solver, second-order-in-space linear interpolation, and the second-order-in-time Runge-Kutta algorithm. We impose the hyperbolic divergence cleaning method for solving the induction equation $\nabla\cdot\bm{B}=0$ (Mignone et al., 2007). The inner boundary is set to be absorbing boundary conditions around black hole (Okuda et al., 2019, 2022, 2023; Aktar et al., 2024). Here, we employ the inner boundary at $R_{\rm in}=2.5r_{g}$ which acts as the event horizon location for our model. Further, we set axisymmetric boundary conditions at the origin, and all the rest are set to outflow boundary conditions (Aktar et al., 2024).

To set up the initial equilibrium torus, we consider the inner edge of the torus at $r_{\rm min}=40r_{g}$ and the maximum pressure surface at $r_{\rm max}=60r_{g}$. We set the minimum density of the torus $\rho_{\rm min}=0.5\rho_{\rm max}$, where $\rho_{\rm max}$ is the maximum density of the torus. We choose $\rho_{\rm max}=\rho_{0}$ for our simulation, where $\rho_{0}$ is the reference density. It is widely accepted in the literature that there are consistent findings regarding the sizes of the torus for the MAD and SANE states in GRMHD simulations (Wong et al., 2021; Fromm et al., 2022). For the MAD state, the inner edge ($r_{\rm min}$) and maximum pressure location ($r_{\rm max}$) are typically situated at approximately $r_{\rm min}\approx 20r_{g}$ and $r_{\rm max}\approx 40r_{g}$ respectively. Conversely, for the SANE state, $r_{\rm min}\approx 6r_{g}$ and $r_{\rm max}\approx 12r_{g}$. Additionally, the initial magnetic field is set differently for the MAD and SANE states (Wong et al., 2021; Fromm et al., 2022; Dihingia et al., 2023). However, the primary focus of this study is to analyze the impact of initial magnetic fields on accretion states (MAD and SANE) when considering the same initial magnetic field configuration (see sub-section 2.5), to compare luminosity, and to investigate the spectral evolution of black holes without altering the geometrical configuration of the torus. For solving radiation fields, we need to supply additionally the radiation energy density $E_{\rm rad}$. In this work, we fix the reference density $\rho_{0}=10^{-10}$ g cm^-3 and the radiation energy density $E_{\rm rad}=5\times 10^{-10}$ erg cm^-3. Also, the mass of the black hole is chosen as $10^{8}M_{\odot}$, and the constant adiabatic index is $\Gamma=4/3$. The conversion relation between the code unit and C.G.S. unit system is summarized in Aktar et al. (2024).

To avoid the situation of encountering negative density and pressure in supersonic and highly magnetized flows, we fix the floor values of density and pressure as $\rho_{\rm floor}=10^{-6}$ and $P_{\rm floor}=10^{-8}$, respectively. Along with that, we also turn on two flags (i) SHOCK FLATTENING to MULTID and (ii) FAILSAFE to YES in PLUTO module (Aktar et al., 2024). We also impose the magnetization condition as $\sigma_{\rm M}<4\pi v_{p}^{2}$ for all the radii for the entire computational domain (Proga & Begelman, 2003; Aktar et al., 2024). Here, $v_{p}$ is the poloidal fluid velocity. Therefore, the floor value of magnetization is $\sigma_{\rm M}^{\rm floor}<4\pi v_{p}^{2}$ in our simulation model.

To investigate the magnetic state of our model, we investigate magnetic flux accumulated at the horizon. In general, we define dimensionless normalized magnetic flux threading to the black hole horizon $(\dot{\phi}_{\rm acc})$ known as the MAD parameter. Here, the MAD parameter is defined as Tchekhovskoy et al. (2011); Narayan et al. (2012); Dihingia et al. (2021); Dhang et al. (2023); Aktar et al. (2024)

where $R_{\rm in}$ is the inner boundary representing the horizon for our model. $B_{r}$ is the radial component of the magnetic field. In our model, we define normalized magnetic flux in Gaussian units by multiplying the factor of $\sqrt{4\pi}$ (Narayan et al., 2022). Here, $\dot{M}_{\rm acc}$ is the mass accretion rate and given as

where the negative sign refers to the inward direction of the accretion flow. In Fig. 2a, we represent the temporal variation of mass accretion rate in Eddington units. Every panel represents different initial plasma-$\beta$ parameters such as $\beta_{0}=10,25,50$, and 100, respectively (see Table 1). Here, we fix the spin of the black hole at $a_{k}=0.98$. We observe that the mass accretion rate saturates close to Eddington mass accretion rate ($\dot{M}_{\rm acc}\sim\dot{M}_{\rm Edd}$) for $\beta_{0}=10$ and 25. However, the mass accretion rate remains in the sub-Eddington limit ($\dot{M}_{\rm acc}<<\dot{M}_{\rm Edd}$) for $\beta_{0}=50$ and 100. One of the good indicators for examining the magnetic state is to investigate the value of the normalized magnetic flux entering through the horizon. Therefore, we also present the corresponding normalized magnetic flux $(\dot{\phi}_{\rm acc})$ in Fig. 2b. The accretion state goes into the MAD state when the saturated value of $(\dot{\phi}_{\rm acc})$ reaches a critical value. The general consensus is that the MAD state is achieved when the critical value of $\dot{\phi}_{\rm acc}\gtrsim 50$ in Gaussian units (Tchekhovskoy et al., 2011; Narayan et al., 2012; Dihingia et al., 2021; Jiang et al., 2023; Aktar et al., 2024). Here, we observe that the $\beta_{0}=10$ model enters into MAD state at $t\gtrsim 9500t_{g}$. Also, the model $\beta_{0}=25$ goes into the MAD state at $t\gtrsim 21000t_{g}$. On the other hand, the model $\beta_{0}=50$ remains in the high magnetized SANE state, and the model $\beta_{0}=100$ represents the low magnetized SANE state $(\dot{\phi}_{\rm acc}<50)$. We also investigate the variability of mass accretion rate and normalized magnetic flux following Narayan et al. (2022). The mean ($\mu$) and its variance around the mean ($\sigma^{2}$) can be defined as $\mu=\frac{1}{N}\sum_{i=1}^{N}x_{i}$ and $\sigma^{2}=\frac{1}{N-1}\sum_{i=1}^{N}(x_{i}-\mu)^{2}$. Here $N$ is number of data points. We find that the variability for $\dot{M}_{\rm acc}$ and $\dot{\phi}_{\rm acc}$ are 2.0483 and 1.1010 for MAD state ($\beta_{0}=10$), respectively. We observe that the variability is much higher for our model compared to the GRMHD simulation in Narayan et al. (2022) because our simulation model has a much higher resolution.

Now, we investigate the behavior of the MAD state and the transition from SANE to MAD state. In Fig. 3, we present the various temporal snapshots for the MAD state of $\beta_{0}=10$. The upper row and lower row in Fig. 3 represent the distribution of gas density $(\rho)$ and plasma-$\beta$ ($\beta$) parameter, respectively. The velocity field lines are indicated by grey lines. Here, we consider the time evolution of torus at $t=7500,9500,25000$, and 50000 $t_{g}$, respectively. At the early simulation time, the accretion rate tends to increase and the flow velocity streamlines remain smooth in nature towards the horizon, as depicted by velocity streamlines in the first column of Fig. 3 at time $t=7500t_{g}$. The accretion flow remains in the SANE state at time $t\leq 9000t_{g}$ as normalized magnetic flux is also below the threshold value ($\dot{\phi}_{\rm acc}<50$) of MAD state, shown in Fig. 2b for $\beta_{0}=10$. Also, the magnetic pressure ($\beta$) remains very low at the horizon at time $t=7500t_{g}$, as depicted in the lower row, the first column of Fig. 3. We also represent the transition time, when the magnetic flux crosses the critical value to become a MAD state, shown in second column at time $t=9500t_{g}$. It is found that the velocity field streamlines start to deviate from a smooth connection to the horizon and the disk region as magnetic pressure increases near the horizon. With time, enough magnetic flux accumulated at the horizon to resist further increase of accretion rate, and the accretion flow becomes MAD in nature, as shown by plasma-$\beta$ distribution in Fig. 3 at times $t=25000t_{g}$ and $t=50000t_{g}$. The increase of magnetic pressure $(\beta<<1)$ near the horizon is also observed, as shown in Fig. 3 at time $t=25000t_{g}$ and $t=50000t_{g}$. Moreover, it is observed that the velocity field lines of the inflowing streamlines turn outward or vertically into a prominent outflow component in the MAD state (Aktar et al., 2024).

Now, we compare the 2D distribution of the other flow variables for MAD and SANE state in Fig. 4. The upper and lower rows represent the MAD ($\beta_{0}=10$) and SANE state ($\beta_{0}=100$), respectively. Here, we compare the distribution of temperature $(T)$, magnetization parameter $(\sigma_{\rm M})$, radiation energy density $(E_{\rm rad})$ and magnitude of vertical velocity ($|v_{z}|$) in the first, second, third and fourth columns for MAD and SANE state, respectively at time $t=50000t_{g}$. The initial dense torus with a temperature of approximately $\sim 10^{9}$ K is surrounded by hot rarefied gas with a temperature of around $\sim 10^{13}$ K, as illustrated in Figure 1b. As time progresses, MRI activity causes an increase in magnetic pressure and synchrotron emission above and below the equator near the horizon. As time further increases, the magnetic field is significantly amplified by MRI in the MAD state compared to the SANE state. As the MAD state approaches its final stage, the torus expands increasingly and a high-velocity jet forms in the funnel region. Consequently, no remnant hot gas is left in the MAD state. However, for the SANE state, the magnetic pressure remains relatively low due to the initially weak magnetic field. The overall flow is smooth, and the initial hot rarefied gas still exists outside the torus for SANE state. We observe that the magnetization parameter is very high ($\sigma_{\rm M}\gtrsim 10$) near the funnel region for the MAD state. Conversely, we find that the magnetization parameter is much lower ($\sigma_{\rm M}<<1$) for SANE. The radiation energy density is much higher ($\gtrsim 10^{3}$) in the MAD state compared to the SANE state. In general, a high magnetization parameter ($\sigma_{\rm M}\gtrsim 1$) indicates the ejection of relativistic jets (Dihingia et al., 2021, 2022; Narayan et al., 2022; Jiang et al., 2023). To investigate further, we present the distribution of vertical velocity magnitude $(|v_{z}|)$ in the fourth column in Fig. 4. Interestingly, we find that the vertical velocity is very high $(|v_{z}|\lesssim c)$ in the funnel region and extends towards the vertical outer boundary for the MAD state. However, the vertical velocity remains at low values (${|v_{z}|}\lesssim 0.01c$) for the SANE state. This indicates that the MAD state is capable of producing highly relativistic magnetized jets or mass outflows (Tchekhovskoy et al., 2011; Narayan et al., 2012; Dihingia et al., 2021; Jiang et al., 2023; Aktar et al., 2024). On the other hand, the SANE state fails to produce high relativistic jets and only produces non-relativistic magnetized mass outflows, as shown in Fig. 4 (Aktar et al., 2024).

In this section, we investigate the effect of black hole spin on the highly magnetized accretion flow. For completeness, we examine the effect of black hole spin on the MAD state. Here we vary only the black hole spin ($a_{k}$) for fixed initial magnetic field $\beta_{0}=10$. We vary the black hole spin as $a_{k}=0.99,0.50,0.0$ and -0.99. In Fig. 11a, we show the variation of mass accretion rate $(\dot{M}_{\rm acc})$ with time by varying the black hole spin $a_{k}$. The corresponding magnetic flux is shown in Fig. 11b. We observe that the accretion rate saturates near the Eddington rate ($\dot{M}_{\rm acc}\sim\dot{M}_{\rm Edd}$) in the MAD state for all the spin cases. Interestingly, we find the prograde flow reaches the MAD state much faster than the retrograde flow, as depicted in Fig. 11b. The MAD state is achieved at $t\sim 9500,11000,13000$ and $15000t_{g}$ for $a_{k}=0.99,0.50,0.0$ and $-0.99$, respectively . Further, we present the distribution of density, temperature, plasma-$\beta$, and magnetization parameters in the first, second, third, and fourth row of Fig. 12, respectively. Here, we consider the simulation time at $t=40000t_{g}$. It is observed that the black hole spin has minimal effect on the accretion process in the torus (Jiang et al., 2023; Aktar et al., 2024). Only the initial torus size slightly reduces for the lower spin values (Aktar et al., 2024). We also observe that MAD state is possible for prograde as well retrograde flow and the accretion properties are almost similar for all the spin values as depicted in Fig. 12. We find high magnetic pressure ($\beta<<1$) near the horizon for all spin values for the MAD state, shown in the third column of Fig. 12. Further, we observe high magnetization parameter $(\sigma_{\rm M}\sim 10)$ near the horizon for all the spin caese which is very common for the MAD state. The high magnetization region is the region for generating highly relativistic jets (Dihingia et al., 2021, 2022; Jiang et al., 2023; Narayan et al., 2022).