© The Authors 2024

1. Introduction

Binary supermassive black holes (SMBHs) are expected to be formed in galaxy mergers. Such binary systems are usually embedded in a gaseous environment at the center of active galactic nuclei (AGNs). Stellar dynamical friction and torques from gas are expected to bring the binary to the subparsec scale (Milosavljević & Merritt 2003). In a binary separated by approximately a parsec, the system emits low-frequency gravitational waves (GWs; on the order of a millihertz) through the inspiral and merger phases, which may be detectable by future observatories such as the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2023). So far, several candidate binary SMBHs have been identified in AGNs and quasars using different observational methods, including NuSTAR and Chandra X-ray observations (Saade et al. 2020, 2024) and periodic features in AGN light curves (Graham et al. 2015). The current strongest candidates are OJ 287 and PG 1302−102, which both display periodicity in their light curves (Bogdanović et al. 2022).

The accretion of gas by the binary can result in appreciable electromagnetic (EM) radiation. On the other hand, the orbital frequency of the binary can be affected by the torque imposed by the gas environment, which may cause a phase shift in the GW signal. The EM emissions and the GW signal together may provide us with a useful probe of the gas in galaxy cores and a diagnostic of the physics of black hole (BH) accretion. In recent years, several analytical and numerical studies have been carried out to investigate this EM emission, mostly in binary systems with equal masses in general relativistic (GR) hydrodynamics (Farris et al. 2011; Shapiro 2013) and GR magnetohydrodynamics (MHD; Noble et al. 2012; D’Ascoli et al. 2018; Bowen et al. 2018; Combi et al. 2022; Avara et al. 2023).

Several groups have studied the interaction of the gaseous circumbinary disk with the binary system via 2D and 3D nonrelativistic hydrodynamic simulations. They include the viscous gas described by the Navier-Stokes fluid equations interacting with an equal-mass circular binary, assuming the binary orbit remains fixed (Moody et al. 2019; Tiede et al. 2020; Mahesh et al. 2023) or evolves with time (Franchini et al. 2022, 2023a). Depending on the assumptions for the live binary orbit, grid resolution, and circumbinary disk parameters, they find that the exchange of angular momentum between the disk and the binary may lead to binary expansion or shrinkage.

Although equal-mass binary studies are more common in the literature, the GW phase shift due to the environmental effect is more significant for the SMBHs with extreme mass ratio inspirals (EMRIs; q ∼ 10⁻⁴ − 10⁻³) and more likely to be detected with the future laser interferometric observatories (Barausse et al. 2014). Based on semi-analytical studies by Yunes et al. (2011) and Kocsis et al. (2011), the accretion disk environment imposes measurable imprints on the GW signal in the EMRIs. These studies show that, depending on the disk parameters, the perturbation of the GW phase is between 10 and 1000 radians per year, which will be detectable by LISA. Barausse & Rezzolla (2008) explored the effects of the BH spin (including prograde and retrograde orbits) and the disk’s mass and density on the orbital evolution of EMRIs. Derdzinski et al. (2021) performed 2D nonrelativistic numerical simulations, inspired by the work of Duffell et al. (2014), to measure the torques on a Jupiter-like planet embedded in a protoplanetary disk. They applied the same code to the AGN context and found that if the disk’s surface density is high enough, the phase shift would exceed a few radians during a 5-year LISA observation. Most recently, Garg et al. (2022) carried out a parameterized study on the gas torque measurements using an analytical approach for intermediate-mass BH binaries embedded in α disks. They quantified the torque over a wide range of Shakura–Sunyaev disk characteristics, such as the surface density and Mach number, as well as the primary BH’s mass and the binary’s mass ratio. Cardoso & Maselli (2020) and Speri et al. (2023) determined how well the environmental effects can be measured using GW observations from LISA, and a new fully relativistic formalism has been introduced by Cardoso et al. (2022) to study GW emission by EMRIs in non-vacuum curved space-times.

In this paper we focus on the case of extreme mass ratio binaries. We simulated a geometrically thin, Keplerian disk orbiting in the plane of a spinning BH. The low-mass companion BH is explicitly not included in our numerical simulation; instead, we used an analytic estimation to measure the viscous torque on the secondary BH. The radiative cooling is neglected in our simulations, and therefore our models are applicable to radiatively inefficient accretion flows (Narayan & Yi 1994), for example those of low-luminous AGNs such as NGC 5548 (Bentz et al. 2007; Crenshaw et al. 2009). Interestingly, in this source the observed profile variability of broad emission lines (Shapovalova et al. 2004) suggests an inspiral effect and that formation of a hot spot due to a collision of the disk with a passing star or a binary BH system (Bon et al. 2016).

Previous works that addressed the environmental effects of the disk on a GW signal assumed the artificial α prescription as the mechanism for the angular momentum transport (Shakura & Sunyaev 1976). In this approach, the α viscosity is assumed to be constant (typically ∼0.01−0.1) for the entire disk and for the entire time of the inspiral phase. However, in a more realistic approach, one needs to include the magnetic field evolution to provide the physical mechanism for the angular momentum transport caused by the magnetorotational instability (MRI; Balbus & Hawley 1991). To include this realistic assumption, we performed full GR MHD simulations. We seeded the initial magnetic field to have a weakly magnetized plasma and evolved the disk in time. We quantified the effective value of α directly from the Reynolds and Maxwell contribution to the stress-energy tensor, as computed for different parts of the disk and over time. Most our simulations are axisymmetric, but to check the importance of non-axisymmetric effects, we also performed a 3D simulation. As the measured torque varies over time and radius, our calculations impose new constraints on the GW phase shift estimation derived with the α disk approach. We evolved the system with GR MHD equations, and therefore our study includes curved space-time and BH spin effects. In addition to the spin, we studied the effects of other physical parameters, such as the magnetic field strength and its configuration.

This paper is organized as follows. In Sect. 2 we describe the initial configuration of our simulations and the unit conversions. The numerical results are presented in Sect. 3 along with a detailed discussion on MRI analysis and α-parameter computation. Section 4 is devoted to our estimations of the viscous torque and its fluctuations. In Sect. 5 we discuss the importance of our results for future GW detections by LISA, with a rough estimation of possible de-phasing due to the gaseous environment, and we give a brief comparison of our α values with the previous studies. Finally, a summary and our conclusions are given in Sect. 6.

2. Methods and setups

2.1. Numerical methods and initial configuration

We used our version of the GR MHD code HARM described in Janiuk (2017) and Sapountzis & Janiuk (2019), which uses numerical algorithms developed initially by Gammie et al. (2003) and Noble et al. (2006). HARM uses a conservative shock-capturing scheme, with fluxes calculated using classical Harten–Lax–van Leer method. The constrained transport maintains divergence free magnetic field. The background space-time is fixed by the Kerr metric of the BH with constant mass and angular momentum. The hydro equations are evolved in the modified spherical Kerr–Schild (KS) coordinates that are nonsingular on the horizon. The following radial and angular maps from KS coordinates to modified Kerr–Schild (MKS) coordinates are used, which increase the resolution close to the BH and the equatorial plane, respectively, to resolve the thin disk accurately: r_KS = exp(r_MKS), ${\theta }_{\mathrm{KS}}=\pi {\theta }_{\mathrm{MKS}}+\frac{(1-h)}{2}sin(2\pi {\theta }_{\mathrm{MKS}})$. The coordinate parameter h is set to 0.3 for all models in this paper.

The gas pressure is calculated using a polytropic equation of state P = κρ^Γ, with κ = 0.1, and Γ = 4/3. The initial density configuration is based on Dihingia et al. (2021) prescription for thin disks. The distribution of density on the equator is defined as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{e}}={\left(\frac{{\mathrm{\Theta }}_0}{\kappa }\right)}^{\frac{1}{(\mathrm{\Gamma }-1)}}{\left(\frac{f(x)}{x^2}\right)}^{\frac{1}{4(\mathrm{\Gamma }-1)}}·\end{array}$$(1)

The disk is truncated at the innermost stable circular orbit radius, r_ISCO, and the density on the entire grid is derived from

$$\begin{array}{c}\hfill \rho (r,\theta )={\rho }_{\mathrm{e}}\exp (-\frac{{\alpha }_{\mathrm{disk}}^2{z}^2}{{\mathcal{H}}^2});z=rcos(\theta ).\end{array}$$(2)

The Θ₀ parameter is the dimensionless temperature set to Θ₀ = 0.001, the parameter α_disk = 2 is chosen to keep the disk thin (H/R ≈ 0.05), and $x=\sqrt{r}$. We refer the readers to Eqs. (4)–(13) from Dihingia et al. (2021) for definitions of f(x) and ℋ.

The initial poloidal magnetic field configuration is based on Zanni et al. (2007) with the nonzero azimuthal component of the magnetic vector potential defined as

$$\begin{array}{c}\hfill A_{\varphi }=r^{3/4}\frac{m^{5/4}}{{(m^2+{cos}^2\theta )}^{5/8}}·\end{array}$$(3)

Here m is a constant and defines the inclination angle of the initial magnetic field. For this study, we chose m = 0.1 (low inclination angle) and m = 0.5 (high inclination angle) for different cases as shown in Fig. 1.

Fig. 1. 2D profile of density and magnetic field lines at the initial time, for β50-m0.5-a0.94 (top) and β50-m0.1-a0.94 cases (bottom).

We normalized the initial field strength with a given value of the ratio of the maximum gas pressure to the maximum magnetic pressure, $\beta =P_{\mathrm{gas}}^{\max }/P_{\mathrm{mag}}^{\max }$. The radial profile of P_g/P_B ratio on the equator at the initial time is shown in Fig. 2.

Fig. 2. Radial profile of the β parameter on the equator at the initial time for all cases.

We summarize the initial parameters of our simulations in Table 1. The main simulations are performed in 2D with axisymmetric assumption and the grid resolution is 1056*528 points in radial and polar directions, respectively. We present a couple of test cases to study the higher resolution and the 3D model in Sect. 3.3. The outer boundary is located at r = 1000 r_g. All cases are evolved for about t ∼ 60 000 t_g. The geometric units r_g and t_g and their relations with physical units are explained in Sect. 2.2.

2.2. Physical scales and physical units

In HARM, we use geometric units G = c = M = 1. To convert quantities to physical units, we follow Janiuk (2019), where the spacial and time units are scaled with the mass of the primary BH as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& L_{\mathrm{unit}}=\frac{\mathit{GM}}{c^2}=1.48\times {10}^5\frac{M}{M_{\odot }}\mathrm{cm},\hfill \\ \hfill & T_{\mathrm{unit}}=\frac{r_{\mathrm{g}}}{c}=4.9\times {10}^{-6}\frac{M}{M_{\odot }}\mathrm{s}.\hfill \end{array}\end{array}$$(4)

The density scale is related to the length unit by ${\rho }_{\text{unit}}=M_{\text{scale}}/L_{\text{unit}}^3$, and the mass scale is adopted to M_scale = 1 × 10⁻⁵ M_⊙ for β50-m0.1-a0.94 case, and M_scale = 2 × 10⁻⁶ M_⊙ for the rest of the models. The density scales are chosen to create disks with surface density Σ ∼ 10³ g cm⁻² around the radius of r ∼ 100 r_g as suggested by Derdzinski et al. (2021). These scaling factors give the measured accretion rate ∼ 0.01 − 1 Ṁ_Edd for our simulations (see Sect. 3 for detailed discussion on the accretion rates).

If we assume the primary BH has a mass of M = 10⁶ M_⊙, the outer boundary is located at r = 1000 r_g ∼ 1.5 × 10¹² cm, and the evolution time is about t ∼ 60 000 t_g ∼ 3.5 days based on Eq. (4). According to these scales, we can claim that our model covers only the inner part of the AGN disks (which extends to 10¹⁴ − 10¹⁶ cm according to Frank et al. 2002), and the evolution time covers only a fraction of the LISA observational time for SMBH inspiral (∼several years; Derdzinski et al. 2021). This evolution time is insufficient to make the entire disk turbulent for the selected resolution. Therefore, we considered only the region inside r < 200 r_g for our torque measurements (assuming the hypothetical secondary BH is orbiting over this range of binary separations). The fluid completes more than 20 orbits at this radius based on the Keplerian orbital frequency.

According to observations, the Seyfert 2 galaxy GSN 069 at a redshift of z = 0.018 with nine-hour X-ray quasi-periodic eruptions is a candidate for extreme mass ratio binary BHs. The primary BH is a low mass SMBH of a few times 10⁵ M_⊙ with a relatively high Eddington ratio of about 0.5 (Miniutti et al. 2019). The variability observed in GSN 069 may be explained by the interaction between an existing accretion disk and an orbiting secondary body according to Linial & Metzger (2023), Franchini et al. (2023b), Tagawa & Haiman (2023), and Miniutti et al. (2023). However, another possible explanation has been suggested based on self-gravitational-lensing in SMBHs by Ingram et al. (2021) for low-mass AGNs such as GSN 069 and RX J1301.9+2747. Our moderately high Eddington ratio models can resemble this type of object; therefore, we scaled our results for such a low-mass primary SMBH (M_p ∼ 10⁵ − 10⁶ M_⊙) with an intermediate-mass companion BH.

3. Numerical results

3.1. Disk evolution and MRI analysis

We evolved our models for about t ≈ 60 000 t_g. This time is long enough to let the MRI be active to form the turbulent structure at the inner part of the disk (r < 150 r_g), and short enough to avoid the magnetic field dissipation due to the anti-dynamo effect in a 2D simulation (the 3D simulation of β50-m0.1-a0.94 case supports this claim; see Sect. 3.3 for more details). In Fig. 3 we show the final configuration of the density and magnetic field lines at the end of the simulations for different models.

Fig. 3. 2D profiles of density and field lines at the final time for different models. The dense and strongly magnetized inner torus is visible in cases β50-m0.1-a0.94 and β10-m0.1-a0.7 separated from the turbulent thin disk. The β50-m0.1-a0.7 has multiple plasmoid structures in the inner part. The model β1-m0.5-a0.7 enters the MAD state at an early time in the evolution, with the magnetic field preserved in a vertical configuration for most parts of the disk.

At the earlier evolution time, for cases β10-m0.1-a0.7, β50-m0.1-a0.7, and β50-m0.1-a0.94, the magnetic field is amplified exponentially due to the MRI and magnetic winding, which creates turbulent magnetized fluid with quite high accretion rates. As a result, the disk expands vertically, launches the outflows and becomes slightly thicker geometrically compared to the initial configuration. However, in these cases, the thin structure of the disk is preserved during evolution, creating a stream that flows radially toward the BH on the equator. This observation suggests that the MRI channel solution is developed in our simulations (channel solutions represent a specific form of poloidal MRI, which is characterized by prominent radial extended features; Balbus & Hawley 1991; Stone & Norman 1994; Dibi et al. 2012).

At later times, we observe that the disks are divided into two distinguishable parts (for these three cases). The inner region is denser, geometrically thicker and more magnetized, while the outer part is geometrically thinner, less dense and less magnetized. The formation of the inner “mini-torus” has been observed before in magnetized thick disks in 2D simulations (De Villiers et al. 2003). In order to investigate the MRI effect and the formation of the inner torus, we compared the radial profiles of the MRI fastest growing mode wavelength, λ_MRI, to the scale height of the disk in Fig. 4. The instability is suppressed when λ_MRI exceeds the scale height (McKinney et al. 2012; White et al. 2019), which happens for r < 20 r_g in β50-m0.1-a0.94 case, where the mini-torus is formed. On the other hand, the visualized data show that the magnetic field lines loop inside the inner torus and form a magnetic barrier at its boundaries, and therefore, create a plasmoid structure. A detailed study on the formation of plasmoids due to magnetic reconnection and their observational effects should be done with high-resolution 3D simulations (Ripperda et al. 2022). With our current 2D moderate-resolution simulation we can explain the existence of the inner torus as a result of two physical processes intensifying each other: (i) the effective MRI at the outer radii makes the fluid lose its angular momentum and be dragged inward, while the less effective MRI at the inner radii makes the hot magnetized fluid slow down and pile up over time, and (ii) at the same time the magnetic field is amplified and creates loops in the inner region causing the plasma trapped and disconnected from the rest of the disk. Therefore, the inner torus is formed and becomes stable till the end of the simulation. For β10-m0.1-a0.7 and β50-m0.1-a0.7 cases, multiple loop structures are visible. Figure 4 shows the comparison of the MRI fastest growing mode’s wavelength λ_MRI, with the disk scale height H at the final time. The scale height is defined as H = c_s/Ω_rot, and the MRI fastest growing wavelength is λ_MRI ≈ 2π|v_θ, A|/|Ω_rot|, where c_s is the adiabatic sound speed, Ω_rot is the fluid rotational frequency, and v_θ, A is the θ-directed Alfvén velocity (cf. McKinney et al. 2012). We observe that for all models, the MRI is suppressed locally at the smaller radii, where the λ_MRI > H. The bottom panel of the same figure shows the surface density profile over the same radial region at the same time. The dense inner torus is distinguishable for β10-m0.1-a0.7, β50-m0.1-a0.7, and β50-m0.1-a0.94 models.

Fig. 4. Comparison of the fastest-growing MRI mode wavelength with the disk scale height at the final time (top) and the surface density at the same time for all the models (bottom). The high-density plasmoid structures are formed at the inner radii, where the MRI is suppressed. The surface density is scaled for the central BH mass, 106 M⊙.

Completing our MRI analysis, we investigated the possibility of a transition to a magnetically arrested disk (MAD), where the disk becomes magnetically dominated and the MRI is suppressed. The MAD status is considered a probable scenario for the disks at the center of a galaxy. The 3D numerical simulations done by Liska et al. (2020) suggested that for a long enough evolution, the disk eventually enters the MAD state, and in this state, the final disk’s characteristics are not sensitive to the exact initial conditions it started with. The recent observations by the Event Horizon Telescope confirm that the MAD state is more favorable for observed sources such as M87 (Event Horizon Telescope Collaboration 2021). To investigate this, we measured the ratio of the magnetic flux to the square root of the mass flux at the BH horizon ${\mathrm{\Phi }}_B/\sqrt{\dot{M}}$ for different cases. Based on the literature, the MAD state happens when this ratio is high enough, ∼15 (Tchekhovskoy et al. 2011). Figure 5 shows that β50-m0.1-a0.94 case, for instance, becomes magnetically arrested for a part of the evolution. At this period of time, the MRI does not act as an effective process for the angular momentum transport, and the accretion rate drops significantly as illustrated in Fig. 6. However, the accretion rate is maintained at around ∼10⁻² M_⊙ yr⁻¹ for β10-m0.1-a0.7 and β50-m0.1-a0.94 cases at the end of the simulations, where the MAD condition is marginally satisfied or below the threshold (i.e., ${\mathrm{\Phi }}_B/\sqrt{\dot{M}}\le 15$).

Fig. 5. Ratio of the magnetic flux to the square root of the mass flux computed at the BH horizon for all the cases. The disks enters the MAD state when this ratio goes above 15. The horizontal line shows the MAD criterion, i.e., the ratio is 15.

At this point, we would like to highlight our β1-m0.5-a0.7 case, which started with a higher inclination angle for the magnetic field initial configuration. The changes are quite dramatic for this case, and it enters the MAD state at an earlier time and in an episodic way (see Figs. 5 and 6). Such episodic accretion rates are commonly observed in MADs (Igumenshchev 2008; Dihingia et al. 2021). More specifically, at the inner part, the magnetic winding causes high magnetic pressure and creates a magnetic barrier that reduces the accretion rate to below 0.01 Ṁ_Edd, while the other cases have accretion rate higher than 0.1 Ṁ_Edd for a long period of evolution. At the further radii, as illustrated in Fig. 4, in a tiny fraction of the disk we have a turbulent structure with the λ_MRI standing below the disk scale height and yet high enough to be resolved with current resolution. However, the vertical configuration of the magnetic field is preserved for the most part of the disk (r > 50 r_g) till the end of the simulation, and makes the disk expand vertically and become less dense compared to the other cases (see Fig. 3). To investigate this further, we performed a test with a similar setup but a weaker initial magnetic field (i.e., β = 50; the result of this test is not presented in our figures). The result shows that even with a weaker initial magnetic field, the MRI cannot be triggered at larger radii and the GR MHD evolution keeps the magnetic field in vertical geometry in most parts of the disk. Therefore, the vertical configuration most likely results in either a weak MRI or a MAD state. Similar works in the literature confirm that a modest change in the initial field configuration may signify MAD structure (Narayan et al. 2003; McKinney et al. 2012; Dihingia et al. 2021).

Fig. 6. Accretion rate for different cases compared with the Eddington accretion limit. The results are scaled for a central BH with a mass of 106 M⊙.

Overall, the results presented in Figs. 4 and 5 determine that the MRI and its effects are not only suppressed locally at the inner part of the disk but are also suppressed globally during the evolution when the disk enters the MAD state.

3.2. α measurement

The MRI analysis presented in Sect. 3.1 shows that this instability does not remain active everywhere in the disk and for the entire time of the evolution. Therefore, the viscous effects driven by MRI vary with radius and time as well. In this section, we compute the equivalent α viscosity caused by the MRI in the turbulent fluid and compare it with the constant values used in the literature.

Following the prescription given by McKinney et al. (2012) for turbulent relativistic fluid, we calculated the equivalent α viscosity by considering the dominant Reynolds and Maxwell terms in the stress-energy tensor as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \alpha & ={\alpha }_{\mathrm{R}}+{\alpha }_{\mathrm{M}},\hfill \\ \hfill {\alpha }_{\mathrm{R}}& \approx \frac{{\rho }_0\delta u_{\mathrm{r}}\delta u_{\varphi }\sqrt{g^{\varphi \varphi }}}{P_{\mathrm{tot}}},\hfill \\ \hfill {\alpha }_{\mathrm{M}}& \approx -\frac{b_{\mathrm{r}}{b}_{\varphi }\sqrt{g^{\varphi \varphi }}}{P_{\mathrm{tot}}}·\hfill \end{array}\end{array}$$(5)

In these equations P_tot = P_gas + P_mag is the total pressure, and b^μ is the magnetic field four-vector (see Eq. (8) from Gammie et al. 2003 for the definition). The four-velocity fluctuations are defined by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \delta u_{\mathrm{r}}(r,\theta ,t)& =u_{\mathrm{r}}(r,\theta ,t)-\overline{u_{\mathrm{r}}(r,t)},\hfill \\ \hfill \delta u_{\varphi }(r,\theta ,t)& =u_{\varphi }(r,\theta ,t)-\overline{u_{\varphi }(r,t)}.\hfill \end{array}\end{array}$$(6)

The average of quantity Q (i.e., velocity and viscosity components) is taken vertically and weighted by density as

$$\begin{array}{c}\hfill \overline{Q(r,t)}=\frac{{\int }_{-H}^H\sqrt{-g}\rho Q\mathrm{d}z}{{\int }_{-H}^H\sqrt{-g}\rho \mathrm{d}z},\end{array}$$(7)

where H is the disk’s scale height. The comparison between the Maxwell and Reynolds contributions to the total volume averaged α over time is shown in Fig. 7 for case β10-m0.1-a0.7. The α values are vertically averaged according to Eq. (7). Figure 8 shows the radial profiles of α_R at three time snapshots, t = 4 × 10⁴, 4.5 × 10⁴, and 5 × 10⁴. These figures show that α_M has a bigger contribution to the averaged α (more than 90%), while α_R fluctuates significantly due to the turbulence. At some radii, α_R may vary in the range of [−0.08, 0.08]. The torque’s fluctuation is discussed in Sect. 5 in detail.

Fig. 7. Volume-averaged αM and αR versus time for case β10-m0.1-a0.7. The αR contribution to the total α is less than 10%.

Fig. 8. Fluctuations in αR at time snapshots t = 40 000, 45 000, and 50 000 computed for β10-m0.1-a0.7. We zoom in on the data for radii [80, 200].

In Fig. 9 we show the volume averaged of α versus time for all the models except β1-m0.5-a0.7, which enters the MAD state periodically at the early time of the evolution. The volume average is taken over the turbulent part of the disk (inner radius of the grid to r = 150 r_g). The comparison between the models for time averaged α is illustrated in Fig. 10. The time average is taken over almost the second half of the evolution (25 000 < t < 60 000) to make sure that the MRI is being triggered and making the disk turbulent up to r ∼ 150 r_g. These results show that the average value of α may vary by factor of 2 for highly magnetized cases such as β10-m0.1-a0.7 and β50-m0.1-a0.94. It reaches 0.3 at the highest and finally converges to α ≈ 0.12 − 0.15 for all the cases including β50-m0.1-a0.7 at the late evolution. The average α drops (with a delay) in all these three cases after entering the MAD status. On the other hand, the radial profile in Fig. 10 demonstrates that α changes significantly over radius: it becomes very small at the inner radii where the MRI is suppressed and gradually increases at the larger radii, reaching α ≈ 0.2 at r ∼ 150 r_g for β50-m0.1-a0.94 and β10-m0.1-a0.7 cases, while the less magnetized case β50-m0.1-a0.7 settles down to α ≈ 0.03 for 100 r_g < r < 200 r_g.

Fig. 9. Volume-averaged α weighted by density versus time for the β10-m0.1-a0.7, β50-m0.1-a0.7, and β50-m0.1-a0.94 cases computed inside the turbulent disk r < 150 and −H < z < H.

Fig. 10. Time-averaged radial profile of total α for the β10-m0.1-a0.7, β50-m0.1-a0.7, and β50-m0.1-a0.94 cases. The time average is taken from the second half of the evolution.

3.3. 3D evolution and resolution effects

The MRI growth rate, saturation and unstable modes are highly dependent on the grid resolution and accurate evolution of the magnetic field in 3D (Shi et al. 2016; Oishi et al. 2020). Based on the Cowling anti-dynamo theorem, the axisymmetric magnetic fields cannot be maintained by axisymmetric motions of a conducting fluid by dynamo action (James et al. 1980; Nunez 1996). On the other hand, it is well known that the disk’s evolution under the MAD condition in 3D is different from 2D (James et al. 2022; White et al. 2019; Tchekhovskoy et al. 2011). To investigate the resolution and the third dimension effects, we performed a high-resolution 2D simulation with 2016*1056 grid points (labeled 2D-High), and a 3D simulation with a moderate resolution of 288*256*96 (labeled 3D) and compare them with the standard 2D resolution (labeled 2D-Std) of the β50-m0.1-a0.94 case.

The first numerical sanity check is whether the MRI turbulence is well captured in our simulations with the proper scaling of the cell sizes. For this purpose we computed the MRI resolution in the polar direction as Res_MRI ≡ λ_MRI, θ/Δθ. The meridional slice of this quantity at t = 50 000 is shown in Fig. 11 for the standard resolution of β50-m0.1-a0.94 case and compared with 3D and 2D-High test cases. We find that our grid always provides at least 10 cells per MRI fastest growing mode’s wavelength in the equatorial region, where the thin disk exists, providing evidence that the MRI turbulence is well resolved within our standard 2D and 3D resolutions.

Fig. 11. MRI resolution defined as ResMRI ≡ λMRI/Δθ in logarithmic scale measured at t = 50 000, for 2D-Std (the standard 2D resolution grid), 2D-High (the high-resolution 2D grid) and 3D (3D grid) models with the same initial setup as β50-m0.1-a0.94. Almost the entire simulation domain has roughly 10 or more grid points per MRI’s fastest growing wavelength, indicating MRI turbulence should be well captured.

Figure 12 shows the density and field line configuration for the 2D-High case, and the meridional and equatorial slices of the same quantities from the 3D simulation at the final time (the color scale is different in the figures). Though in the 3D case, the disk remains geometrically thin and the inner torus with a looping magnetic field is not formed. The other significant difference between the 2D and the 3D simulation is the accretion rate, which stays close to 0.5 − 1 Ṁ_Edd for the entire simulation (see Fig. 13). It is also obvious that there is no sudden transition to the MAD state at the early evolution compared to the 2D simulation. The MAD state develops later, where the ${\mathrm{\Phi }}_B/\sqrt{\dot{M}}$ ratio gradually increases (Fig. 13), though the simulation is not long enough to observe the magnetic barrier close to the horizon and therefore causes no significant drop in the accretion rate. The 2D-High case, on the other hand, resulted in a similar qualitative evolution as the standard resolution: a highly magnetized, plasmoid structure is created at the inner radii, distinguishable from the turbulent thin disk with an MRI channel solution at the outer radii. It also shows frequent transitions to the MAD state and significant drops in the accretion rate.

Fig. 12. Density profile and magnetic field lines for the 2D-High case (left), the meridional slice at ϕ = 0 (middle), and the equatorial slice of the same quantities for the 3D simulation (right) at the final snapshot.

The volume averaged α viscosity measured for each case is shown in the bottom panel of Fig. 13¹. Both the 3D and high-resolution 2D cases demonstrate deviations from the standard 2D case at the earlier time of the evolution, which is expected for MRI. However, these values start converging in the middle of the simulation. The 3D evolution of α follows the standard 2D case closely with some fluctuations, while the average α for 2D-High goes through minor changes over time and it settles down to α ≈ 0.12 during the second half of the evolution. These tests show that, despite the different early evolution, our ⟨α⟩ measurements for 2D simulations are reliable for the second half of the evolution and are not affected by the anti-dynamo theorem.

Fig. 13. Ratio of the magnetic flux to the square root of the mass flux on the horizon (top), the accretion rate (middle), and the volume averaged α viscosity measured for the 2D-Std (the standard 2D resolution grid), 2D-High (the high-resolution 2D grid), and 3D (3D grid) models with the same initial setup as β50-m0.1-a0.94 (bottom).

4. Analytical results: Viscous torque measurements

In this section, we apply the numerical results from our GR MHD simulations, reported in Sect. 3, to estimate the viscous torque in a 1D GR approach. Since the secondary BH is not explicitly included in our simulation, we can make different assumptions about its mass and different scenarios it may go through while inspiraling in the turbulent gaseous environment. The results in this section are scaled for a primary BH mass of 10⁶ M_⊙ and mass ratio of 10⁻³ orbiting within the disk’s inner turbulent part, at r < 200 r_g.

In a realistic numerical model, where the secondary BH is included, one needs to calculate all the different components of the total torque exerted on the secondary – gravitational torque, accretion on the secondary, and viscous torques – to measure how much the secondary’s orbit is influenced by the environment (Derdzinski et al. 2021). However, there are analytical approximations in the literature to estimate the torque. Historically, these approximations were applied to explain the planet radial migration in the protoplanetary disks: According to Tanaka et al. (2002) and Lyra et al. (2010) the linear torque can be estimated for two types of migration: (1) type-I migration for a very low mass companion with mass ratio (q < 10⁻⁴), and (2) type-II migration for a higher mass ratio such as q ∼ 10⁻³. In type-II migration the companion is massive enough to carve a low-density region (gap) inside the disk and therefore, experiences a lower torque caused by the viscosity (Lin & Papaloizou 1986). A similar 1D type-II migration model was adopted by Armitage & Natarajan (2002) for SMBHs merging in accretion disks, which was further advanced by Shapiro (2013). In the latter, the evolution equation of the surface density was derived in the curved space-time, considering the gas effect on the secondary BH (viscous torque) and the secondary BH’s gravitational tidal effects on the disk (tidal torque). Our assumption for the hypothetical mass ratio fits the type-II migration; therefore, our torque estimations are limited to the viscous torque computed relativistically using the Shapiro 1D approach.

In Sect. 4.1 we present the viscous torque measurements from the time-averaged values, and in Sect. 4.2 we discuss the viscous torque fluctuations and their effects on the orbital evolution.

4.1. Viscous torque: 1D GR-hybrid thin disk model

The Newtonian 1D thin disk prescription used by Garg et al. (2022) and Derdzinski et al. (2021) to compute the viscous torque for type-II migration followed

$$\begin{array}{c}\hfill T_{\nu ,\mathrm{Newt}}=-3\pi r^2{\mathrm{\Omega }}_2\nu \mathrm{\Sigma },\end{array}$$(8)

where ${\mathrm{\Omega }}_2\approx \sqrt{M_{\mathrm{p}}/r^3}$ is the orbital frequency of the secondary BH assuming circular orbit for the extreme mass ratio case, M_p is the mass of the primary BH, $\mathrm{\Sigma }={\displaystyle {\int }_{-H}^H}\rho \mathrm{d}z$ is the disk surface density and ν = αc_sH is the kinematic viscosity.

In a relativistic formalism, we followed the simple 1D GR-hybrid model given by Shapiro (2013). In this approach, the rest mass accretion rate due to the viscous torque is

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{GR}}=2\pi \left[\frac{\mathcal{G}}{\mathcal{Q}}3r^{1/2}\frac{\partial }{\partial r}\left(r^{1/2}\nu {\mathrm{\Sigma }}_{\mathrm{GR}}\frac{{\mathcal{D}}^2}{\mathcal{C}}\right)\right],\end{array}$$(9)

where 𝒢, 𝒬, 𝒞, 𝒟 and relativistic surface density, Σ_GR, are defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{G}& =\mathcal{B}{\mathcal{C}}^{-1/2},\hfill \\ \hfill L^{+}& =M_{\mathrm{p}}x\mathcal{C}(1-2a{x}^{-3}+a^2{x}^{-4}),\hfill \\ \hfill \mathcal{Q}& =2x^{1/2}\partial L^{+}/\partial r,\hfill \\ \hfill \mathcal{B}& =1+a{x}^{-3},\hfill \\ \hfill \mathcal{C}& =1-3x^{-2}+2a{x}^{-3},\hfill \\ \hfill \mathcal{D}& =1-2x^{-2}+a^2{x}^{-4},\hfill \\ \hfill {\mathrm{\Sigma }}_{\mathrm{GR}}& ={\displaystyle {\int }_{-H}^H}\rho u^t\sqrt{-g}\mathrm{d}z,\hfill \\ \hfill \end{array}\end{array}$$(10)

for the Kerr metric in the Boyer–Linquist coordinates, where M_p and a are the mass and spin of the primary BH and $x=\sqrt{r}$. The relativistic viscous torque will be computed as

$$\begin{array}{c}\hfill T_{\nu ,\mathrm{GR}}=-{\dot{M}}_{\mathrm{GR}}{r}^2{\mathrm{\Omega }}_2.\end{array}$$(11)

We get nearly identical results as from Eq. (8) in the weak field approximation, where 𝒢, 𝒬, 𝒞, and 𝒟 approach unity. Our results show that the Newtonian approach (Eq. (8)) overestimates the viscous torque, especially at the inner regions, for instance, the relativistic viscous torque is ∼30% lower at r ∼ 100 r_g.

On the other hand, the binary is inspiraling due to losing energy by emitting GWs. So, one can define the effective GW torque to be compared with the viscous torque:

$$\begin{array}{c}\hfill T_{\mathrm{GW}}=\frac{1}{2}q{M}_{\mathrm{p}}r{\dot{r}}_{\mathrm{GW}}{\mathrm{\Omega }}_2,\end{array}$$(12)

where q is the mass ratio and ṙ_GW is derived by the quadruple approximation for the evolution of the orbital separation as follows (Peters 1964):

$$\begin{array}{c}\hfill {\dot{r}}_{\mathrm{GW}}=-\frac{64}{5}\frac{{(GM)}^3}{c^5}\frac{1}{1+q^{-1}}\frac{1}{1+q}\frac{1}{r^3}·\end{array}$$(13)

Figure 14 shows the ratio T_ν, GR/T_GW computed with q = 0.001 for different models in our simulations. The values for Σ, α and c_s used in Eq. (9) are taken from our simulations’ time-averaged numerical measurements. This result shows that the average viscous torque may reach up to a few per cent of the gravitational torque and produce a measurable phase shift in the GW signal (for the selected mass ratio and M_p). Moreover, considering the results from β50-m0.1-a0.7 and β1-m0.5-a0.7 cases, we can claim that there are particular magnetic field strengths and configurations, which result in less effective MRI turbulence, and hence, lower viscous torques.

Fig. 14. Ratio of the viscous torque to the effective GW torque for all models scaled for a fixed primary BH mass of Mp = 106 M⊙ and a mass ratio of q = 0.001.

In comparison with results from an analytical (and Newtonian) study done by Garg et al. (2022), which considered a constant value α = 0.01, we conclude that our numerical results provide higher values for α. Therefore, the secondary BH experiences a higher viscous torque at binary separations around r ∼ 100 r_g. Obviously, T_ν is scaled by other quantities such as the surface density, sound speed and scale height. Hence, the GW phase shift can still be used to probe the density and temperature of AGN disks for α ∼ 0.1 at this binary separation.

4.2. Torque fluctuations

What we estimate as a torque in the previous section is the time-averaged value of viscous torque, or equivalently we can call it the linear torque. However, in a realistic scenario, we do not have a 1D laminar gas flow. Instead, we have to deal with nonlinear turbulent fluid that continuously interacts with the secondary orbiting BH. These nonlinear hydrodynamic interactions lead to rapid changes in density, velocity and magnetic fields, which enhance or suppress the torque value over time. The deviation from the linear torque can affect the orbital evolution of the binary and introduce an additional phase shift in GWs. Zwick et al. (2022) presented a detailed study on the stochastic torque or time-variable torque estimations and their measurable effects in the GW signals. As they suggested, there are two important sources for these fluctuations: one is the disk-driven fluctuations experienced by the low-mass orbiting object from the turbulent fluid, and the second, is the perturber-driven fluctuations occurring due to asymmetries in the gas flow near a sufficiently massive secondary BH.

The perturber-driven fluctuations may depend on many physical processes including the stochastic accretion and tidal effects of the secondary BH, as well as gas friction and small-scale gas dynamics. On the other hand, for fluctuation studies, one needs to measure the torque directly from a relatively long simulation in the presence of the perturber and output it frequently for the Fourier analysis. The Newtonian torque measurements and its Fourier analysis have been done by Nelson (2005) for turbulent, magnetized protoplanetary disks. They found the torque fluctuations contain high-amplitude, low-frequency components, and the stochastic migration dominates over type-I migration for their models evolved for 100−150 planet orbits. For a relativistic approach, one can follow the direct torque computations from Farris et al. (2011, see their Appendix B for details of the relativistic derivation of torque components).

In our study, the torque fluctuations do not appear in the average values we presented in Sect. 4.1, and overall we observe that the average viscous torque is negative (by its definition) and therefore, it facilitates the shrinkage of the binary orbit. However, it is still interesting to investigate the fluctuations of the viscous torque and have a qualitative discussion on this topic. (Similar to in Sect. 4.1, we assume that the dominant torque component is the viscous part for our chosen mass ratio, and we limited our study of the fluctuations to the viscous torque.) Taking a close look at the case β50-m0.1-a0.94, for instance, shows that the viscous torque may deviate from the time-averaged torque by up to factor of 5. Figure 15 shows the ratio of the viscous torque over the time-averaged torque versus time for two selected radii (r = 50 r_g and r = 75 r_g). As we observe this ratio can change dramatically over time and it even may have a different sign during evolution. Measuring a positive viscous torque may seem incorrect. However, we should remind the readers that when α is computed from the Reynolds and Maxwell contributions, it may take negative values by the definition in a turbulent fluid (see α_R in Fig. 8 for instance). We defer the direct computations of the torque and its fluctuations to future studies in which the orbiting perturber will be included in our simulations.

Fig. 15. Ratio of the viscous torque to the time-averaged viscous torque versus time at radii r = 50 rg and r = 75 rg for the case β50-m0.1-a0.94. The fluctuations in the measured torque affect the orbital evolution of the binary system.

5. Discussion: Phase shift due to gas, α viscosity, and other important physical scenarios

Using α and the linear torque computed in Sects. 3.2 and 4, we can estimate the orbital evolution and GW de-phasing due to the gas environment. Here we follow the analysis given by Zwick et al. (2022) derived based on Tanaka et al. (2002) and Ward (1997), so the flux of angular momentum induced by the viscous torque is estimated from local α, density and sound speed as

$$\begin{array}{c}\hfill {\dot{L}}_{\mathrm{T}}=-\alpha \frac{6\pi r^{7/2}{c}_{\mathrm{s}}{(r)}^3\rho (r)}{\sqrt{\mathit{GM}}},\end{array}$$(14)

and we can derive the variation in the binary separation from

$$\begin{array}{c}\hfill \dot{r}={\dot{r}}_{\mathrm{GW}}+2\frac{{\dot{L}}_{\mathrm{T}}}{\mathit{Mq}}\sqrt{\frac{r}{\mathit{GM}}}\equiv {\dot{r}}_{\mathrm{GW}}+{\dot{r}}_{\mathrm{T}}.\end{array}$$(15)

Finally, the phase shift is approximated by

$$\begin{array}{c}\hfill \delta \varphi ={\varphi }_{\mathrm{vac}}-\varphi \approx 2\pi {\displaystyle \int f_{\mathrm{GW}}}(r)\frac{{\dot{r}}_{\mathrm{gas}}}{{\dot{r}}_{\mathrm{GW}}^2}\mathrm{d}r.\end{array}$$(16)

Approximating the GW frequency for a binary with mass ratio 0.001 and M_p = 10⁶ M_⊙ from Derdzinski et al. (2021), we predict the de-phasing would be roughly around ∼10 radians for about 10⁵ inspiral orbits. In comparison with the analytical work by Garg et al. (2022) and the numerical study by Derdzinski et al. (2021), which used constant α ∼ 0.01, our computed viscous torque is slightly larger because our directly measured α is larger by more than one order of magnitude. Generally, for LISA, there are sources with an S/N from the order of unity up to a few hundred, and the phase of the GW signal can be reconstructed within the accuracy of 1 over the S/N (Thorpe et al. 2019). Therefore, the predicted phase shift must be detectable by LISA for a few years of observational time (Kocsis et al. 2011; Garg et al. 2022; Derdzinski et al. 2021).

At this point, it is worth having a brief discussion over the α viscosity we computed and comparing it with similar works from the literature. Our numerical measurement for average α viscosity as ⟨α⟩=⟨α_M⟩+⟨α_R⟩ from the mid-evolution time estimates this value around ∼0.1 − 0.25 for different cases, which varies from almost zero close to the horizon and reach up to ∼0.2 at radii around 150 r_g after taking the time average. However, in the 1D analytical approach and viscous hydrodynamic simulations, it is very common to assume smaller constant values on the order of 0.001 − 0.02 (Shibata et al. 2017) for α. This assumption is based on high-resolution shearing box simulations presented in the literature such as Shi et al. (2016) and Salvesen et al. (2016). In 2007 King et al. (2007) estimated the α viscosity as α ≈ 0.1 − 0.4 for different astrophysical systems based on the observation. This study explained a factor of 10 discrepancy between observational and theoretical estimates of α as a result of the restrictions in the local models, and suggested undertaking fully 3D, global simulations for a realistic measurement. Since then, several convergence test studies have been done for global accretion disks including Hawley et al. (2013) and Suzuki & Inutsuka (2014). The nonrelativistic Maxwell stress (−⟨B_rB_ϕ⟩/4πP_gas) reported in Hawley et al. (2013) varies from 0.17 to 0.46 for different models, and what they defined as the Shakura–Sunyaev viscosity, or α_SS, is the density-weighted α_M and reported as α_SS ≈ 0.03 − 0.04 for their highest resolution. They also claimed that the resolution also affects the value of Maxwell stress or α_mag, though not dramatically. Our 2D models also show slight changes in volume-integrated α for a different resolution; however, the relativistic definition of α is different in our measurements, and it includes the Reynolds stress as well. Suzuki & Inutsuka (2014) has reported the density-weighted α_mag, with the same nonrelativistic definition, within the same order of magnitude (∼0.01) for 3D global disk simulation with initial vertical field configuration. Moreover, the most recent studies by Mishra et al. (2020, 2022) reported the value of α_mag in the range of ∼0.2 − 0.4 for a magnetized thin disk model with H/R ≈ 0.05.

We should emphasize that all the results presented in this work are order-of-magnitude estimations for torque measurements and their effects on the GW detections. This study is limited in many ways, most importantly for a realistic approach, one needs to include the radiative cooling for optically thin AGN disks. In such a model, radiative cooling competes with viscous heating to reach thermal equilibrium, which overall affects the thermal evolution and the geometry of the disk significantly (Narayan et al. 1998). However, our study can still be considered as a good approach for AGNs with lower radiation efficiency. So far, several GR MHD groups have included radiative transfer processes to study the emission spectrum from the accretion disks in their simulations including Noble et al. (2011), McKinney et al. (2014), Dexter (2016), Younsi et al. (2016), Bronzwaer et al. (2018, 2020), Mościbrodzka (2020), Chatterjee et al. (2020), and Dihingia et al. (2023). The low-mass orbiting object is another missing part of our simulation. based on the discussion given by Zwick et al. (2022) the perturber-driven torque fluctuations can be affected by the disk parameters and become noticeably important in GW detections. Moreover, the presence of the low-mass secondary BH provides an opportunity to study Lindblad or orbital resonance torques in numerical simulations. Armitage & Natarajan (2002) have observed the formation of spiral waves that mediate angular momentum exchange between the secondary and the disk. It is worth mentioning that besides gas-induced torques such as the MRI-driven viscosity and Lindblad resonance, the tidally distorted primary BH’s horizon can gravitationally couple to the orbiting low-mass secondary, transferring energy and angular momentum from the BH to the orbit and causing an additional phase shift in the GW signal (O’Sullivan & Hughes 2014).

Evolving a supermassive binary BH in magnetized fluid for hundreds of orbits (or longer) in a high-resolution grid and with all the relevant physics included (such as radiation) requires extremely expensive GR MHD simulations in dynamic space-time. However, for future studies, it is possible to simplify the problem in the case of extreme mass ratios by evolving only MHD equations in a fixed space-time and adding the secondary BH as a perturber. Such approximations are taken by Combi et al. (2021) and Suková et al. (2021) for adjusting metric and hydrodynamic equations, respectively. The secondary’s accretion can also be added as an additional sink term in the source part of the hydro equations (see Eq. (6) from Derdzinski et al. 2021 as an example).

6. Summary and conclusions

This study was one step toward a realistic estimation of a disk’s environmental effects on possible future GW detections. We evolved several magnetized thin disk models to quantify the viscous α parameter in a turbulent fluid developed by MRI. We used the results of these simulations to estimate the viscous torque experienced by the hypothetical low-mass secondary BH inspiraling the primary BH inside an AGN disk. Finally, we estimated the phase shift in the GW signal caused by the disk’s environmental effect based on the torque magnitude.

We observe that the disks with well-resolved MRI have an average α viscosity that varies around 0.1 − 0.25 during the second part of the evolution. The MRI is suppressed at the inner part of the disk, close to the primary BH, so the value of the α viscosity is negligible in this region. However, the time-averaged α reaches ≈0.1 − 0.2 at larger radii, where the fluid is turbulent and the MRI’s fastest-growing mode is resolvable.

Altering the initial conditions, we found that an initial magnetic field with a lower inclination angle plays an important role in triggering and sustaining MRI, while a field configuration with a higher inclination angle brings the disk into the MAD state, with episodic high magnetic flux at the horizon. In the initially weakly magnetized case, the MRI saturates at a later time and overall has a smaller α and, therefore, a lower viscous torque compared to the other cases. Moreover, we find that the BH spin does not change the results significantly.

Since 3D evolution is essential for capturing all MRI features, we ran one 3D simulation. The results of this simulation show a different dynamical evolution of the disk, especially regarding the development of the MAD state. Nevertheless, the average value of the α viscosity agrees well with the 2D simulation for a long enough evolution.

We applied the numerical results from the GR MHD simulations to estimate the viscous torque using the GR-hybrid approach for the GR 1D thin disk. We find that the time-averaged viscous torque can be as high as ∼1% of the GW torque for a mass ratio of q = 10⁻³ at radii around r ∼ 100 r_g, where α is maximal. This extra torque from the environment appears as a faster shrinkage of the binary’s orbit and a phase shift in the GW signal. We also observe that the Newtonian-calculated torque deviates from the relativist one by up to 30% (higher) at radii around r ∼ 100 r_g.

Monitoring the viscous torque at different radii shows that the fluctuations in the torque values can change dramatically, and sometimes the torque’s sign even changes or it deviates from the average value by a factor of 5. The study of torque fluctuations is essential to understanding a binary’s orbital evolution and should be considered in future numerical studies where the perturbation from a secondary BH is explicitly included.

Acknowledgments

The authors thank Andrea Derdzinski, Bożena Czerny, Hector Olivares, Scott Noble, Milton Ruiz and Jose Font for helpful discussions and advice throughout this project. This work was supported by grant No. 2019/35/B/ST9/04000 from the Polish National Science Center, Poland. We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2023/016071. We also acknowledge support from the Interdisciplinary Center for Mathematical Modeling of the Warsaw University.

References

All Tables

All Figures

	Fig. 1. 2D profile of density and magnetic field lines at the initial time, for β50-m0.5-a0.94 (top) and β50-m0.1-a0.94 cases (bottom).
In the text

	Fig. 2. Radial profile of the β parameter on the equator at the initial time for all cases.
In the text

Fig. 3. 2D profiles of density and field lines at the final time for different models. The dense and strongly magnetized inner torus is visible in cases β50-m0.1-a0.94 and β10-m0.1-a0.7 separated from the turbulent thin disk. The β50-m0.1-a0.7 has multiple plasmoid structures in the inner part. The model β1-m0.5-a0.7 enters the MAD state at an early time in the evolution, with the magnetic field preserved in a vertical configuration for most parts of the disk.

In the text

	Fig. 4. Comparison of the fastest-growing MRI mode wavelength with the disk scale height at the final time (top) and the surface density at the same time for all the models (bottom). The high-density plasmoid structures are formed at the inner radii, where the MRI is suppressed. The surface density is scaled for the central BH mass, 106 M⊙.
In the text

	Fig. 5. Ratio of the magnetic flux to the square root of the mass flux computed at the BH horizon for all the cases. The disks enters the MAD state when this ratio goes above 15. The horizontal line shows the MAD criterion, i.e., the ratio is 15.
In the text

	Fig. 6. Accretion rate for different cases compared with the Eddington accretion limit. The results are scaled for a central BH with a mass of 106 M⊙.
In the text

	Fig. 7. Volume-averaged αM and αR versus time for case β10-m0.1-a0.7. The αR contribution to the total α is less than 10%.
In the text

	Fig. 8. Fluctuations in αR at time snapshots t = 40 000, 45 000, and 50 000 computed for β10-m0.1-a0.7. We zoom in on the data for radii [80, 200].
In the text

	Fig. 9. Volume-averaged α weighted by density versus time for the β10-m0.1-a0.7, β50-m0.1-a0.7, and β50-m0.1-a0.94 cases computed inside the turbulent disk r < 150 and −H < z < H.
In the text

	Fig. 10. Time-averaged radial profile of total α for the β10-m0.1-a0.7, β50-m0.1-a0.7, and β50-m0.1-a0.94 cases. The time average is taken from the second half of the evolution.
In the text

Fig. 11. MRI resolution defined as ResMRI ≡ λMRI/Δθ in logarithmic scale measured at t = 50 000, for 2D-Std (the standard 2D resolution grid), 2D-High (the high-resolution 2D grid) and 3D (3D grid) models with the same initial setup as β50-m0.1-a0.94. Almost the entire simulation domain has roughly 10 or more grid points per MRI’s fastest growing wavelength, indicating MRI turbulence should be well captured.

In the text

	Fig. 12. Density profile and magnetic field lines for the 2D-High case (left), the meridional slice at ϕ = 0 (middle), and the equatorial slice of the same quantities for the 3D simulation (right) at the final snapshot.
In the text

	Fig. 13. Ratio of the magnetic flux to the square root of the mass flux on the horizon (top), the accretion rate (middle), and the volume averaged α viscosity measured for the 2D-Std (the standard 2D resolution grid), 2D-High (the high-resolution 2D grid), and 3D (3D grid) models with the same initial setup as β50-m0.1-a0.94 (bottom).
In the text

	Fig. 14. Ratio of the viscous torque to the effective GW torque for all models scaled for a fixed primary BH mass of Mp = 106 M⊙ and a mass ratio of q = 0.001.
In the text

	Fig. 15. Ratio of the viscous torque to the time-averaged viscous torque versus time at radii r = 50 rg and r = 75 rg for the case β50-m0.1-a0.94. The fluctuations in the measured torque affect the orbital evolution of the binary system.
In the text