Revisited equilibrium solution of Fishbone and Moncrief torus for extended GRMHD simulations

General Relativity (GR) has been widely accepted due to its profound predictions and observational verification (e.g., Will, 2014). More recently, the detection of the first gravitational waves by the Laser Interferometer Gravitational-Wave Observatory (LIGO) from the binary black hole (BH) merger (Abbott et al., 2016) and the images of the supermassive BH by the Event Horizon Telescope (EHT) at the center of the radio galaxy M 87 (Akiyama et al., 2019a) and at the center of the Milky Way Sgr A^∗ (Akiyama et al., 2022a) established GR to be one of the well-defined theories of gravity. However, the GR itself has seen some fundamental issues in the observations and theoretical perspectives, such as the existence of the singularity at the center of the BH spacetime (Hawking & Ellis, 1973), the accelerated expansion of the Universe (Riess et al., 2001), and the renormalization of the theory of GR (’t Hooft & Veltman, 1974; Goroff & Sagnotti, 1986). Therefore, testing the theory of gravity around compact objects is very important to clarify it.

The general relativistic magnetohydrodynamic (GRMHD) simulation has been an essential tool in the study of the dynamics of astrophysical phenomena near the extreme gravitational field of the compact object in relativistic regimes (e.g., Banyuls et al., 1997; Koide et al., 1999; Gammie et al., 2003; De Villiers et al., 2003). It has been used to study numerous astrophysical phenomena, including the formation of relativistic jets, winds, hot spots, the structure of the hot accretion flow around the BH, and neutron star systems (e.g., Narayan et al., 2012; Liska et al., 2018; Parfrey & Tchekhovskoy, 2017; Nathanail et al., 2019; Davis & Tchekhovskoy, 2020; Porth et al., 2021; Dihingia et al., 2021; Mizuno et al., 2021; Begelman et al., 2022; Dihingia et al., 2022; Chatterjee & Narayan, 2022; Mizuno, 2022; Das & Porth, 2024; Dihingia & Fendt, 2024). The EHT has also provided a large library of 3D GRMHD simulations of the magnetized hot accretion flows in the stationary, axisymmetric spinning vacuum solution of the Einstein field equation: Kerr spacetime (Akiyama et al., 2019b, 2022b). We can test GR by comparing the synthetic images produced by these libraries with the EHT observed images (Akiyama et al., 2019c, 2022c). A diverse work on the semi-analytical model has also been studied to test the GR (e.g., Narayan et al., 2019; Psaltis et al., 2020; Kocherlakota et al., 2021; Bauer et al., 2022; Ozel et al., 2022; Ayzenberg, 2022; Younsi et al., 2023; Vagnozzi et al., 2023; Uniyal et al., 2023, 2024; Jiang et al., 2024b). However, further testing of the theory of gravity is needed to go beyond Kerr spacetime. The simulation community as well as EHT have already made some efforts in this direction (e.g., Mizuno et al., 2018; Olivares et al., 2020; Fromm et al., 2021; Akiyama et al., 2022c; Nampalliwar et al., 2022; Röder et al., 2023). These simulations and calculations show that images have high degeneracy when it comes to the deformation in Kerr spacetime. Therefore, further studies of horizon-scale dynamics and imaging are inevitable (e.g., Chatterjee et al., 2023b, a; Jiang et al., 2024a).

In order to perform the GRMHD simulations in non-Kerr spacetime, we first need to build initial conditions that are consistent with them. It depends on spacetime, and therefore constructing the initial setup, such as the hydrodynamical equilibrium torus in a general framework is an essential requirement. So far, in the GRMHD simulations of magnetized accretion flow around Kerr black holes, most of the work has been done by considering the Fishbone and Moncrief (FM) torus solution (Fishbone & Moncrief, 1976) which is applicable not only to Kerr spacetime but also to any stationary and axisymmetric spacetime in GR. In this paper, we apply the idea of the FM torus to any general stationary and axisymmetric spacetime and show the stability of the torus. Previously, initial conditions of equilibrium torus for non-Kerr spacetimes were provided by Font & Daigne (2002) which is commonly known as FD torus. Although the FD torus has already been used for several GRMHD simulations of non-Kerr spacetime (Mizuno et al., 2018; Gimeno-Soler et al., 2019; Olivares et al., 2020; Teodoro et al., 2021; Gimeno-Soler et al., 2021; Röder et al., 2023; Cassing & Rezzolla, 2023), the FM torus is widely used to perform GRMHD simulations in Kerr spacetime. It is important to note that other stable tori solutions also exist in the literature such as magnetized tori with a toroidal magnetic field (Komissarov, 2006), tori with non-constant angular momentum (De Villiers et al., 2003; Gimeno-Soler & Font, 2017; Lei et al., 2009) and tori with the alternative initial conditions (Penna et al., 2013). Therefore, building stationary FM tori in generic stationary and axisymmetric backgrounds could provide an alternative framework to study accretion flow in non-Kerr spacetimes.

To test the generic FM torus solution in non-Kerr spacetime, we perform general relativistic hydrodynamic (GRHD) simulations in the Johannsen-Psaltis (JP) parametrized metric (Johannsen & Psaltis, 2011; Johannsen, 2013a, b) which deforms the Kerr metric by the deformation parameters. We show that the general FM torus is stable in long-term evolution, which is a desired characteristic of initial conditions.

This paper is structured in the following ways: We start with the mathematical formulation of the general FM torus in Section 2. We then show the stability of the general FM torus by GRHD simulations in Section 3 and conclude our discussion in Section 4. We perform the simulations in geometrized units: the speed of light, $c=1$, the gravitational constant, $G=1$, and the mass of the central black hole, $M=1$ using the GRMHD code BHAC (Porth et al., 2017; Olivares et al., 2019). The length and time are expressed in units of $r_{g}=GM/c^{2}$ and $t_{g}=r_{g}/c$, respectively.

We consider a general form of the stationary, axisymmetric, asymptotically vacuum or non-vacuum flat spacetime in spherical polar coordinates $(t,r,\theta,\phi)$ as

where, $\nu$, $\psi$, $\omega$, $\mu_{1}$, and $\mu_{2}$ are the function of the coordinates. In this case, all the coefficients can be written in terms of the metric components in the Boyer-Lindquist coordinates as,

where all the metric components depend on the radial coordinate $r$ and polar coordinate $\theta$. To construct the hydrostatic equilibrium torus, we would like to have the conserved angular momentum per unit inertial mass ($l=u_{\phi}u^{t}$), pressure ($p$), and the rest-mass density of the matter ($\rho$). We start with the integral form of the relativistic Euler equation for isentropic fluid $dS=0$ (Fishbone & Moncrief, 1976),

where $u^{\rho}=(u^{t},0,0,u^{\phi})$ is the four-velocity of the fluid and $h=(e+p)/\rho$ is the specific enthalpy with total energy density, $e$. In the Newtonian limit, $W$ is the effective potential with contributions from the centrifugal and gravitational forces. Therefore, in the hydrodynamic equilibrium, $W$ can be defined as the equation of the balance by the pressure gradient in the following way,

here $w=e+p$. The subscript $"in"$ hereon will represent the value of the physical quantity in the inner edge ($r_{in}$) of the torus at the equatorial plane ($\theta=\pi/2$) where the pressure vanishes. Using the polytropic equation of state $p=K\rho^{\gamma}$, where $\gamma$ and $K$, are the adiabatic index and polytropic constant respectively, we can integrate the above equation which gives,

The distribution of density and pressure of the torus are expressed as,

where

and

Note that the equation of state here considered is the polytropic equation of state which satisfied $\displaystyle{\lim_{P\to 0}h=\displaystyle{\lim_{P\to 0}\frac{\rho+P}{n}}=1}$, $n$ is the baryon number density. Hence, the boundary of the disk will be the surface where $\ln(h)=0$ and therefore here $\ln(h_{in})=0$ (Fishbone & Moncrief, 1976). This further helps in calculating $W_{in}$ which will be the first three terms in equation (8) at $r=r_{in}$ and $\theta=\pi/2$. The angular momentum per unit inertial mass $l$ can then be expressed as the solution of the pressure gradient equation ($[\ln(h)]_{,r}=0$ at $\theta=\pi/2$) following (Fishbone & Moncrief, 1976),

Here $,r$ denotes the partial derivative with respect to radial coordinate ($r$) at $\theta=\pi/2$. The equation gives two solutions for the angular momentum per unit inertial mass, one representing the counter-rotating disk and the other for the co-rotating disk. The co-rotating disk solution leads to the following expression:

where

The above equation is the generalized version of Eq.$(3.8)$ of Fishbone and Moncrief for a generic stationary and axisymmetric spacetime (Fishbone & Moncrief, 1976). The four-velocity components are also needed to construct the complete structure of the torus. First, we set the radius of the maximum density of the torus ($r_{\rm max}$). By supplying $r_{\rm max}$ to Eq. (10), we obtain the constant angular momentum for the torus. It is important to note that since the pressure is required to increase outwards from $r_{\rm in}$, i.e., $[\ln(h)]_{,r}>0$ at $r=r_{\rm in}$, we can not fix the angular momentum at $r=r_{in}$. This further restricts the inner edge of the disk to be $r_{\rm in}>r_{\rm ph_{+}}$, where $r_{\rm ph_{+}}$ is the retrograde photon radius of the circular photon orbit. Note that for $r<r_{\rm ph_{+}}$, the formation of the circular orbits needs acceleration due to the pressure gradient outwards. However, at $r=r_{\rm in}$, it is directed inwards. Such criteria also prevent the disk from having any cusp at the inner edge (Kozlowski et al., 1978; Font & Daigne, 2002). In the case of the massive particle, as it approaches the photon sphere, it moves more and more rapidly, which results in the divergence of angular momentum and energy per unit mass. Therefore, the positive root of $B=0$ will provide the photon sphere radius. A general analytical expression is difficult to find, thus we provide the prograde and retrograde photon sphere radius for the lowest order parameter variation of the JP metric in Tabel 1.

Further, we use this angular momentum to calculate the initial conditions of the four-velocity components. We can calculate the azimuthal component of the contravariant four-velocity by the following transformations:

where $u_{(\phi)}$ is the azimuthal velocity component in the locally non-rotating frame (LNRF) and the basis vectors in this frame can be transformed as (Bardeen et al., 1972),

Therefore, with the help of the above transformations, we can write the four-velocity transformations as follows (Bardeen et al., 1972),

which reduce $u^{\phi}$ to be,

Following Fishbone & Moncrief (1976), $u_{(\phi)}$ can be expressed as,

The time-component of the contravariant four-velocity, $u^{t}$ can be calculated using the normalization condition $g_{\alpha\beta}u^{\alpha}u^{\beta}=-1$ where $u^{r}=u^{\theta}=0$ while constructing the torus. Therefore, the azimuthal component of the four-velocity of the torus can be constructed using Eq. (12) and the above-constructed relations. The pressure can be calculated within the torus by using Eq. (6). Therefore, the structure of the torus has been formulated in this section. We will test the stability of the constructed generic FM torus in the next section (Sec. 3).

In the previous section, we described the procedure to construct general FM torus in non-Kerr spacetime. In this section, we test the stability of the newly constructed general FM torus through numerical simulations. To represent the non-Kerr metric, we use the Johannsen-Psaltis (JP) parametrized metric in the Boyer-Lindquist (BL) coordinates, which is given as (Johannsen, 2013a):

the non-zero metric components are given by,

where $A_{1}$, $A_{2}$, and $A_{5}$ are the function of radial coordinate ($r$) only and,

The coefficients $\epsilon_{n},\alpha_{1n},\alpha_{2n}$, and $\alpha_{5n}$ are the deviation parameters. For our purposes, we stick to the first-order contribution of the deformation parameters and vary one of these parameters in simulations with $\epsilon_{3}\equiv\epsilon$. It is important to note that, as the deformation parameter vanishes, it returns to the Kerr metric.

We perform 2D GRHD simulations of general FM torus in Kerr and non-Kerr spacetimes on a spherical polar grid where the grid spacing is logarithmic in the radial direction and linear in the polar direction. We implemented the horizon penetrating form of the JP metric in the code BHAC (Porth et al., 2017) using Eq. (111) of (Johannsen, 2013a) to represent the non-Kerr spacetime.

As an initial condition, we apply the general FM torus, which is in hydrostatic equilibrium. We use the Kerr metric and the metric with deformation from Kerr spacetime. In order to make a fair comparison, we fix the total mass content within the torus. Here we choose the inner radius of the torus at $r_{in}=6\,r_{g}$ and the density maximum radius at $r_{\rm max}=12\,r_{g}$ which provides the angular momentum. In that case, the mass content inside the torus remains the same. To check the stability of the torus, we add a random perturbation in the pressure up to $4\%$. In all the simulation runs, we use the same random seed number so that we can have the same initial setup. The simulation here is performed for black hole spin $a=0.9375$, adiabatic index $\gamma=5/3$, and polytropic constant $K=0.001$.

We treat the outside of the torus as an atmosphere that is filled with very low-density matter. We fix the velocity in the Eulerian frame to be zero. Pressure and density set the very low values of $\rho=\rho_{\rm min}r^{-1.5}$ and $p=p_{\rm min}r^{-2.5}$, where $\rho_{\rm min}=10^{-5}$ and $p_{\rm min}=10^{-7}$ respectively, to avoid the unphysically small or negative values during the simulation run.

The effective grid resolution of the simulations is $1024\times 512$ grids with two levels of static mesh refinement. At the inner and outer radial boundaries, we apply standard inflow/outflow boundary conditions by copying the physical quantities. At the polar boundaries, reflecting boundary conditions are applied. For azimuthal direction, periodic boundary conditions are employed.

Fig. 1 shows the distribution of the density for the Kerr spacetime (panels 1(a) and 1(b)) and non-Kerr spacetime with the deformation $\alpha_{22}=0.5$ (panels 1(c) and 1(d)) at $t=0$ and $t=10,000\,t_{g}$. We note that for the non-Kerr spacetime case, we set zero for the other deformation parameters. From the $2D$ density distributions, it is clearly shown that the torus remains stable, but the matter outside of the torus falls onto a BH for long-term evolution. We expect that some small changes may happen due to the given initial perturbation within the torus. Nevertheless, the results remain similar in both cases. We found similar results for all non-Kerr spacetime cases by adding different deformation parameters. We do not display them explicitly here to avoid repetition.

We further check the stability of the torus structure quantitatively by the vertically integrated radial density and pressure profiles, which are shown in Fig. 2. We varied the one deformation parameter ($a_{13}$, $a_{22}$, $a_{52}$, or $\epsilon$) for each non-Kerr case and compared their radial density and pressure profiles at $t=0\,t_{g}$ (solid black line) and $t=10,000\,t_{g}$ (solid blue line). The values of the non-zero Kerr parameters are written on each panel. We observed that the main torus profiles remain stable in the long-term evolution, i.e., unchanged in their initial structure even at time $t=10,000\,t_{g}$. However, due to the small perturbation in the torus, we see some small changes at the surface of the torus in all the cases, including the case for Kerr spacetime (Fig. 2a).

The final confirmation of the stability check comes from the oscillation of maximum density $(\rho_{\rm max})$ within the torus caused by the initial perturbation, as shown in Fig. 3. In panel (a), we present the evolution of $\rho_{\rm max}$ with time, where we observe the random small oscillations clearly that decay (damping) in long-term evolution for all cases. Even up to $t=10,000\,t_{g}$, we see such random oscillations. The value of $\rho_{\rm max}$ increases with time since all the perturbed high-pressure regions started to move to the location of the density maximum to make the torus equilibrium. The general trend has not changed in the different spacetimes. To further characterize these oscillations, we show the corresponding power spectral density (PSD) for different spacetime cases in panel (b). From this figure, we do not observe any dominant peaks in the PSD due to the provided random perturbation. Dominate peaks in PSD usually correspond to quasi-periodic oscillations (QPOs) in the emission profile (Zanotti et al., 2005; de Avellar et al., 2018), However, such studies need to be done in the MHD, where the turbulent magnetic field can give rise to more oscillation in the evolution of the initially constructed torus. These results will be reported soon in the upcoming work. The important point is that for all the cases, the PSD is mostly consistent, showing the robustness of the initial torus. Overall, the long-term GRHD simulations show that generic FM torus remains stable and can be used for further studies of magnetized accretion flows in the GRMHD simulations. To make our claim more robust, we study the stability of the generalized FM torus by changing values of $r_{\rm in}$ and $r_{\rm max}$ and also in one more considered metric (Kerr-Sen metric). The results of these studies are shown in Appendix A and B, respectively.

For the sake of completeness, we compare the initial structures ($t=0$) of two torus solutions, the FD torus and the FM torus for non-Kerr spacetime in Appendix C. Both torus have the same inner and outer edge radius. However, in the interior of the torus, there are some small deviations. We will further expand the study to investigate the influence of these differences by using GRMHD simulations, which will be reported soon.

Testing the theory of gravity has become more important recently. To understand the dynamics of matter near the central compact objects, GRMHD simulation is one of the key tools. High-resolution $3D$ GRMHD simulations in Kerr spacetime have been used to study the magnetized accretion flows near the supermassive BH in Sgr A^∗ and M 87^∗ (Akiyama et al., 2019b, 2022b). However, the testing theory of gravity will need to perform GRMHD simulations in non-Kerr spacetimes. To do it, the initial conditions need to be formed carefully, which is consistent with the spacetime characteristics. In Kerr spacetime, we have widely used one of the hydrostatic equilibrium torus solutions, the so-called FM torus to perform GRMHD simulations. In this paper, we present more general accretion tori solutions using the FM approach in stationary, axisymmetric, non-Kerr spacetimes.

To investigate the stability of our torus solution, we performed 2D GRHD simulations of FM tori in non-Kerr spacetimes. We fixed the inner radius as well as the mass content within the torus for all cases and introduced a $4\%$ of random perturbation to the gas pressure in the torus. The JP parametrized metric is used as an example for the deformation from Kerr spacetime. By checking the density and pressure distribution, we show that the torus remains stable for long-term evolution, even when the deformation parameters are non-zero. Furthermore, we looked into the oscillation of the maximum density within the torus and its corresponding PSD to study the effect of the perturbation on the torus structure. We found that the FM torus in non-Kerr spacetime behaves similarly to the Kerr spacetime. Therefore, we conclude that generic FM tori are stable and ready to apply for the GRMHD simulations of magnetized accretion flows in non-Kerr spacetime. It is important to note that the simulation performed in the paper is only for Kerr-like spacetime however the generality of the procedure is valid in all axisymmetric, stationary, and vacuum spacetimes. Therefore we can apply the same FM tori solution in any non-Kerr spacetime.

In the past, the simulation library of magnetized accretion flows onto a black hole in Kerr spacetime has been well-developed and used to understand the near-horizon physics from the black hole shadow images observed by EHT (Akiyama et al., 2019b, 2022b). However, the simulation library of magnetized accretion flows onto a black hole beyond GR spacetime has not been much developed yet. Although it is essential for testing the theory of gravity via black hole shadow images. Therefore, the development of a new simulation library using non-Kerr spacetime is an urgent task. It is important to note that a horizon-penetrating form of the most general axisymmetric, stationary, asymptotically flat spacetime to use generic FM torus is needed for the long-term evolution of GRMHD simulations. Recently, several new horizon-penetrating forms of the most general axisymmetric, stationary, asymptotically flat spacetime have been proposed (e.g., Johannsen, 2013a; Konoplya et al., 2021; Kocherlakota et al., 2023; Ma & Rezzolla, 2024). Here, we computed generic initial conditions for the FM torus solution for any non-Kerr stationary, axisymmetric, asymptotically vacuum, or non-vacuum flat spacetime. Further application work using the MHD simulations will be reported soon.