Self-gravitating magnetised tori around black holes in general relativity

TL;DR

This paper models stationary, self-gravitating, magnetised tori around black holes in general relativity, exploring how magnetic fields influence the structure and properties of these accretion disks.

Contribution

It introduces the first numerical models of magnetised, self-gravitating tori around black holes, including magnetic fields in the Einstein and MHD equations.

Findings

Magnetic pressure shifts density maxima inward.

Magnetic fields influence total mass and angular momentum.

Models cover a wide range of magnetisation levels.

Abstract

We investigate stationary, self-gravitating, magnetised disks (or tori) around black holes. The models are obtained by numerically solving the coupled system of the Einstein equations and the equations of ideal general-relativistic magnetohydrodynamics. The mathematical formulation and numerical aspects of our approach are similar to those reported in previous works modeling stationary self-gravitating perfect-fluid tori, but the inclusion of magnetic fields represents a new ingredient. Following previous studies of purely hydrodynamical configurations, we construct our models assuming Keplerian rotation in the disks and both spinning and spinless black holes. We focus on the case of a toroidal distribution of the magnetic field and build a large set of models corresponding to a wide range of values of the magnetisation parameter, starting with weakly magnetised disks and ending at configurations in which the magnetic pressure dominates over the thermal one. In all our models, the magnetic field affects the equilibrium structure of the torus mainly due to the magnetic pressure. In particular, an increasing contribution of the magnetic field shifts the location of the maximum of the rest-mass density towards inner regions of the disk. The total mass of the system and the angular momentum are affected by the magnetic field in a complex way, that depends on the black hole spin and the location of the inner radius of the disk. The non-linear dynamical stability of the solutions presented in this paper will be reported elsewhere.

Tables

Table 1. Table 1: Parameters of the numerical solutions. In all cases we assumed the polytropic exponent γ = 4 / 3 𝛾 4 3 \gamma=4/3 , the magnetisation law parameter n = 1 𝑛 1 n=1 , and the black hole mass parameter m = 1 𝑚 1 m=1 . From left to right the columns report: the black hole spin parameter a 𝑎 a , the coordinate inner radius of the disk r 1 subscript 𝑟 1 r_{1} , the circumferential inner radius of the disk r C , 1 subscript 𝑟 C 1 r_{\mathrm{C,1}} , the coordinate outer radius of the disk r 2 subscript 𝑟 2 r_{2} , the circumferential outer radius of the disk r C , 2 subscript 𝑟 C 2 r_{\mathrm{C,2}} , the maximum rest-mass density within the disk ρ max subscript 𝜌 max \rho_{\mathrm{max}} , the parameter C 1 subscript 𝐶 1 C_{1} appearing in the magnetisation law, the total ADM mass m ADM subscript 𝑚 ADM m_{\mathrm{ADM}} , the mass of the black hole m BH subscript 𝑚 BH m_{\mathrm{BH}} , the angular momentum of the disk J 1 subscript 𝐽 1 J_{1} , and the magnetisation parameter β mag subscript 𝛽 mag \beta_{\mathrm{mag}} .

No.	$a$	$r_1$	$r_{C,1}$	$r_2$	$r_{C,2}$	${\rho }_{\max }$	$C_1$	$m_{\mathrm{ADM}}$	$m_{\mathrm{BH}}$	$J_1$	${\beta }_{\mathrm{mag}}$
1a	$-0.5$	8.0	9.2	35.3	36.7	$5\times {10}^{-5}$	0	1.33	1.01	1.7	$\infty$
1b	$-0.5$	8.0	9.2	35.3	36.7	$5\times {10}^{-5}$	0.01	1.34	1.01	1.75	30.5
1c	$-0.5$	8.0	9.3	35.3	36.8	$5\times {10}^{-5}$	0.1	1.40	1.01	2.08	3.49
1d	$-0.5$	8.0	9.4	35.3	36.9	$5\times {10}^{-5}$	1	1.51	1.02	2.65	0.21
1e	$-0.5$	8.0	9.4	35.3	36.9	$5\times {10}^{-5}$	1.3	1.50	1.02	2.58	$5.3\times {10}^{-2}$
1f	$-0.5$	8.0	9.3	35.3	36.9	$5\times {10}^{-5}$	1.42	1.49	1.02	2.54	$1.3\times {10}^{-3}$
2a	0	8.1	9.3	35.1	36.5	$5\times {10}^{-5}$	0	1.33	1.02	1.64	$\infty$
2b	0	8.1	9.3	35.1	36.5	$5\times {10}^{-5}$	0.01	1.34	1.02	1.69	29.4
2c	0	8.1	9.3	35.1	36.5	$5\times {10}^{-5}$	0.1	1.40	1.02	2.02	3.37
2d	0	8.1	9.4	35.1	36.7	$5\times {10}^{-5}$	1	1.52	1.03	2.61	0.19
2e	0	8.1	9.4	35.1	36.7	$5\times {10}^{-5}$	1.3	1.51	1.03	2.55	$3.0\times {10}^{-2}$
2f	0	8.1	9.4	35.1	36.7	$5\times {10}^{-5}$	1.37	1.50	1.03	2.52	$5.8\times {10}^{-4}$
3a	0.9	3.0	4.4	20.0	21.7	$3.5\times {10}^{-4}$	0	1.52	1.00	2.04	$\infty$
3b	0.9	3.0	4.4	20.0	21.7	$3.5\times {10}^{-4}$	0.01	1.52	1.00	2.05	$75.8$
3c	0.9	3.0	4.4	20.0	21.7	$3.5\times {10}^{-4}$	0.1	1.55	1.00	2.17	$8.38$
3d	0.9	3.0	4.4	20.0	21.7	$3.5\times {10}^{-4}$	1	1.57	1.01	2.23	0.96
3e	0.9	3.0	4.4	20.0	21.6	$3.5\times {10}^{-4}$	2	1.47	1.00	1.79	0.26
3f	0.9	3.0	4.4	20.0	21.5	$3.5\times {10}^{-4}$	2.74	1.39	1.00	1.45	$5.88\times {10}^{-4}$
4a	0.99	0.8	2.41	20.1	21.9	$1.5\times {10}^{-3}$	0	1.70	1.00	2.31	$\infty$
4b	0.99	0.8	2.41	20.1	21.9	$1.5\times {10}^{-3}$	0.01	1.70	1.00	2.30	805.5
4c	0.99	0.8	2.40	20.1	21.9	$1.5\times {10}^{-3}$	0.1	1.68	1.00	2.24	80.3
4d	0.99	0.8	2.38	20.1	21.7	$1.5\times {10}^{-3}$	1	1.51	1.00	1.64	7.72
4e	0.99	0.8	2.35	20.1	21.5	$1.5\times {10}^{-3}$	2	1.32	1.00	1.01	3.07
4f	0.99	0.8	2.33	20.1	21.3	$1.5\times {10}^{-3}$	3	1.17	1.00	0.52	1.31
4g	0.99	0.8	2.32	20.1	21.3	$1.5\times {10}^{-3}$	4	1.08	1.00	0.24	0.39
4h	0.99	0.8	2.32	20.1	21.2	$1.5\times {10}^{-3}$	4.5	1.05	1.00	0.17	0.11
4i	0.99	0.8	2.32	20.1	21.2	$1.5\times {10}^{-3}$	4.7	1.05	1.00	0.15	2.28 $\times {10}^{-2}$

Equations

$$T_{\mu \nu }=(\rho h+b^2)u_{\mu }{u}_{\nu }+\left(p+\frac{1}{2}{b}^2\right)g_{\mu \nu }-b_{\mu }{b}_{\nu },$$

$$b_{\mu }{u}^{\mu }=0.$$

$$g=g_{tt}d{t}^2+2g_{t\phi }dtd\phi +g_{rr}d{r}^2+g_{\theta \theta }d{\theta }^2+g_{\phi \phi }d{\phi }^2,$$

$$b_t=-\frac{u^{\phi }}{u^t}{b}_{\phi }=-\Omega b_{\phi },$$

$$g_{tt}+2g_{t\phi }\Omega +g_{\phi \phi }{\Omega }^2=-\frac{1}{(u^t{)}^2}.$$

$$b_{\phi }^2=-(u^t{)}^2\mathcal{L}{b}^2,$$

$${\nabla }_{\mu }T\text{indices}{}_{\nu }^{\mu }=\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left(\sqrt{-g}T\text{indices}{}_{\nu }^{\mu }\right)-\frac{1}{2}({\partial }_{\nu }{g}_{\alpha \beta })T^{\alpha \beta }=0$$

$${\partial }_{\nu }\left(p+\frac{1}{2}{b}^2\right)-\frac{1}{2}({\partial }_{\nu }{g}_{\alpha \beta })\left[(\rho h+b^2)u^{\alpha }{u}^{\beta }-b^{\alpha }{b}^{\beta }\right]=0.$$

$$\frac{1}{2}{u}^{\alpha }{u}^{\beta }{\partial }_{\nu }{g}_{\alpha \beta }=\frac{{\partial }_{\nu }{u}^t}{u^t}-u^t{u}_{\phi }{\partial }_{\nu }\Omega$$

$$\frac{1}{2}{b}^{\alpha }{b}^{\beta }{\partial }_{\nu }{g}_{\alpha \beta }=b^2\left(\frac{{\partial }_{\nu }{u}^t}{u^t}+\frac{{\partial }_{\nu }\mathcal{L}}{2\mathcal{L}}-u^t{u}_{\phi }{\partial }_{\nu }\Omega \right).$$

$$\rho h\left(u^t{u}_{\phi }{\partial }_{\nu }\Omega -\frac{{\partial }_{\nu }{u}^t}{u^t}\right)+{\partial }_{\nu }p+\frac{{\partial }_{\nu }(b^2\mathcal{L})}{2\mathcal{L}}=0,$$

$$u^t{u}_{\phi }{\partial }_{\nu }\Omega -\frac{{\partial }_{\nu }{u}^t}{u^t}+\frac{{\partial }_{\nu }h}{h}+\frac{{\partial }_{\nu }(b^2\mathcal{L})}{2\rho h\mathcal{L}}=0,$$

$$\int u^t{u}_{\phi }d\Omega +\ln \left(\frac{h}{u^t}\right)+\int \frac{d(b^2\mathcal{L})}{2\rho h\mathcal{L}}=C,$$

$$\int j(\Omega )d\Omega +\ln \left(\frac{h}{u^t}\right)+\int \frac{d(b^2\mathcal{L})}{2\rho h\mathcal{L}}=C.$$

$$g=(-{\alpha }^2+{\beta }_i{\beta }^i)d{t}^2+2{\beta }_{i}d{x}^{i}dt+{\gamma }_{ij}d{x}^{i}d{x}^j,$$

$$n^{\mu }=\frac{1}{\alpha }(1,-{\beta }^r,-{\beta }^{\theta },-{\beta }^{\phi }),n_{\mu }=(-\alpha ,0,0,0).$$

$$D_j(K^{ij}-{\gamma }^{ij}K)=8\pi j^i$$

$$\frac{1}{2}\left(R+K^2-K_{ij}{K}^{ij}\right)=8\pi {\rho }_{\mathrm{H}},$$

$$K_{ij}=-\frac{1}{2\alpha }\left({\partial }_t{\gamma }_{ij}-{\mathcal{L}}_{\beta }{\gamma }_{ij}\right),$$

$${\mathcal{L}}_{\beta }{\gamma }_{ij}={\beta }^k{\partial }_k{\gamma }_{ij}+{\gamma }_{ik}{\partial }_j{\beta }^k+{\gamma }_{kj}{\partial }_i{\beta }^k.$$

$${\rho }_{\mathrm{H}}=n^{\mu }{n}^{\nu }{T}_{\mu \nu },j_i=-P\text{indices}{}_i^{\alpha }{n}^{\beta }{T}_{\alpha \beta },$$

$${\partial }_t{K}_{ij}-{\mathcal{L}}_{\beta }{K}_{ij}$$

$$S_{ij}=P\text{indices}{}_i^{\mu }P\text{indices}{}_j^{\nu }{T}_{\mu \nu },$$

$${\rho }_{\mathrm{H}}$$

$$S_{ij}$$

$$S$$

$$j_{\phi }=\alpha \rho h{u}^t{u}_{\phi }.$$

$$g$$

$$K_{ij}=\frac{1}{2\alpha }\left({\gamma }_{ik}{\partial }_j{\beta }^k+{\gamma }_{kj}{\partial }_i{\beta }^k\right),$$

$$D_{j}K\text{indices}{}_l^j=8\pi j_l.$$

$$K_{\phi r}$$

Full text

I Introduction

Compact accretion disks (or tori) around black holes are astrophysical transient systems that can form in a number of situations. Examples include the core-collapse of massive stars woosley:2006 , the merger of compact binaries consisting of either two neutron stars or a black hole and a neutron star (see e.g. baiotti2017 and references therein), and the gravitational collapse of a supermassive star rees1984 ; SMS . Observations of the formation and evolution of black hole–torus systems are challenging, either using neutrino, electromagnetic or gravitational-wave approaches. Since the recent breakthrough observation of gravitational waves from a binary neutron star (BNS) merger by Advanced LIGO and Virgo GW170817 ; GW-GRB one may hope that this cosmic messenger may offer the best possibility of observing black hole–torus systems in the near future.

Numerical relativity is the best approach to study the dynamical formation of black hole–torus systems from ab-initio simulations. Long-term simulations of BNS mergers that include the late inspiral of the two neutron stars and account for the relevant physics (i.e. relativistic gravity, general-relativistic magnetohydrodynamics (GRMHD), and neutrino transport) are in general fairly expensive. Therefore, building equilibrium initial data of black hole–torus systems is highly motivated, as it allows to carry out follow-up studies of the last stages of the merger in a less expensive way and in a more controlled environment, sidestepping the computation of the late inspiral and early merger phase. Equilibrium models must therefore be as faithful as possible to the end-products of the numerical evolutions, increasing their realism as new physical ingredients are incorporated (see abramowicz:2013 and references therein). Numerical works have shown that the mass of the tori may be large enough to render necessary to account for the disk self-gravity in order to properly describe its dynamics. This is particularly true for the case of unequal-mass BNS mergers baiotti2017 ; rezzolla10 . Motivated by these results, we present in this paper new families of self-gravitating disks around black holes.

A few authors have previously investigated this issue nishida ; ansorg ; shibata ; stergioulas . In their seminal work, Nishida and Eriguchi nishida computed self-gravitating toroids around stars and black holes using Komatsu-Eriguchi-Hachisu’s (KEH) self-consistent-field method komatsu . Elliptic-type field equations were converted into integral equations using Green’s functions. Later on, Ansorg and Petroff ansorg built solutions of black holes surrounded by uniformly rotating rings of constant density using the same approach as nishida , but solving the equations with a highly accurate multi-domain, pseudo-spectral method. A similar strategy was followed by Stergioulas stergioulas to construct general-relativistic models of self-gravitating, constant angular momentum tori around black holes with KEH’s self-consistent-field method. An important ingredient of this approach was the use of a compactified radial coordinate, which improved the enforcement of the boundary conditions asymptotically. Among all previous studies, the most relevant one for our work is that of Shibata shibata , since we follow very closely his procedure. Shibata’s work departs from the other three approaches in that it builds self-gravitating tori around rotating black holes adopting the so-called puncture framework to describe the spacetime of a rotating black hole, and hence avoiding potential numerical issues when dealing with the curvature singularity at the origin. The models reported by shibata are purely hydrodynamical (i.e. with no magnetic field), and they are characterized by the constant angular momentum. (Non-constant angular momentum tori were considered in kiuchi2011 albeit around non-rotating black holes.) On the contrary, the models presented in this paper incorporate a toroidal distribution of the magnetic field, a physically motivated Keplerian rotation law kkmmop ; kkmmop2 and rotating black holes.

In addition to self-gravity, our disks also incorporate magnetic fields, within the ideal GRMHD approach. To the best of our knowledge, equilibrium sequences of self-gravitating and magnetised disks around black holes in general relativity have not yet been reported in the literature. (Notice, however, that Appendix A of shibata outlines the procedure to build a self-gravitating magnetised disk with a toroidal magnetic field, but no examples are provided.) There exist a number of previous works where equilibrium solutions of magnetised disks around black holes have been built komissarov:2006 ; montero:2007 ; gf but, to the best of our knowledge, all of them are restricted to the test-fluid approximation (i.e. neglecting self-gravity). Komissarov komissarov:2006 first presented a general procedure to build magnetised “Polish doughnuts” (constant angular momentum tori) using a barotropic equation of state and the assumption that the specific enthalpy of the fluid is close to unity. This restrictive condition on the thermodynamics was relaxed in the work of Montero et al. montero:2007 , who also performed dynamical evolutions of those tori. More recently, Gimeno-Soler and Font gf built new sequences of equilibrium magnetised tori around Kerr black holes assuming a form of the angular momentum distribution proposed in qian:2009 that departs from the constant case of komissarov:2006 and from which the equipotential surfaces can be easily computed.

The study of the stability of equilibrium solutions of accretion tori under perturbations has received considerable numerical attention (see abramowicz:2013 for a review). In particular, constant angular momentum disks have been found to be generically unstable. On the other hand, in most BNS merger simulations the final black hole–torus system does not manifest signs of dynamical instabilities on short dynamical timescales (see (baiotti2017, ) and references therein). Specifically, the simulations of rezzolla10 indicate that the angular velocity $\Omega$ of tori formed from unequal-mass BNS mergers follows Keplerian profiles, $\Omega\propto r^{-3/2}$, where $r$ denotes the distance from the rotation axis, which explains the scaling of the specific angular momentum as $r^{1/2}$. This provides firm evidence that tori produced self-consistently are dynamically stable. However, despite their non-constant angular momentum profiles make them stable against the development of the so-called runaway instability (abramowicz, ; Font2002, ; montero, ), on longer timescales non-axisymmetric instabilities (e.g. the Papaloizou-Pringle instability (PPI) (PPI, )) set in (korobkin2011, ; kiuchi2011, ; mewes16a, ; mewes16b, ). Recently Bugli et al. bugli studied the development of the PPI in tori threaded by weak toroidal magnetic fields and how this instability may be affected by the concurrent development of the magnetorotational instability (MRI). Their simulations, within the test-fluid limit, showed that the magnetic fields provide local viscous stresses through turbulence and global angular momentum transport, leading to the suppression of large-scale PPI modes. The self-gravitating, magnetised tori we built in the present work may thus be used in the future to investigate the generality of the findings of bugli beyond the test-fluid limit.

The paper is organised as follows: in Sec. II we discuss the mathematical aspects of our procedure, presenting the Euler-Bernoulli equations, the Einstein equations, and the Keplerian rotation law. The masses and angular momentum of the black hole–torus spacetime are discussed in Sec. III. Section IV briefly describes our numerical method and the results are discussed in Sec. V. Finally, Sec. VI gives a summary of this work. In the Appendix we provide expressions for the Kerr metric in quasi-isotropic coordinates. We use geometric units with $G=c=1$, where $G$ is Newton’s constant and $c$ is the speed of light, and assume the signature of the metric $(-,+,+,+)$. Spacetime dimensions are labeled with Greek indices, $\mu=0,1,2,3$, while Latin indices are used to denote spatial dimensions, $i=1,2,3$.

II Equations

We start by deriving the equations describing stationary, axially symmetric, self-gravitating, magnetised toroids rotating around black holes. The black hole (which can be spinning) is included in the system using the puncture method. The torus is described in terms of ideal GRMHD; we restrict ourselves to toroidal magnetic fields and barotropic equations of state. In specific numerical examples discussed in Sec. V we assume polytropic fluids and a Keplerian rotation law introduced recently in kkmmop ; kkmmop2 .

As mentioned in the introduction the formulation presented in this paper is based on the approach to modeling stationary perfect-fluid disks around black holes derived originally in shibata for the case with no magnetic fields. In particular, in the derivation of the equations we follow closely the steps taken in shibata . Because the terms connected with the magnetic field appear irregularly in many of the equations, we repeat the corresponding calculations also in this paper. The magnetic field enters the description of stationary disks in two places: in the stationary Euler equation and in the “source terms” of the Einstein equations.

Within the framework of ideal GRMHD the energy-momentum tensor has the form

[TABLE]

where $\rho$ is the baryonic density, $h$ is the specific enthalpy, $p$ is the (thermal) pressure, $u^{\mu}$ denotes components of the four-velocity of the fluid, and $b^{\mu}$ is the four-vector of the magnetic field. We denote $b^{2}=b_{\mu}b^{\mu}$. Note that the quantity $p_{\mathrm{mag}}=\frac{1}{2}b^{2}$ plays the role of a magnetic pressure. It is assumed that

[TABLE]

In this case the dual of the Faraday tensor relative to an observer with four-velocity $u^{\mu}$, ${}^{\ast}F^{\mu\nu}=b^{\mu}u^{\nu}-b^{\nu}u^{\mu}$, satisfies $\nabla_{\mu}{}^{\ast}F^{\mu\nu}=0$.

We will work in spherical coordinates $(t,r,\theta,\varphi)$. It is convenient to start with a general, stationary and axially symmetric metric of the form

[TABLE]

where the metric potentials $g_{tt}$, $g_{t\varphi}$, $g_{rr}$, $g_{\theta\theta}$, $g_{\varphi\varphi}$ depend on $r$ and $\theta$ only. We consider an axially symmetric, stationary configuration with $u^{r}=u^{\theta}=b^{r}=b^{\theta}=0$. It follows from Eq. (2) that

[TABLE]

where $\Omega=u^{\varphi}/u^{t}$. Note that the normalization of the four-velocity $u_{\mu}u^{\mu}=-1$ yields

[TABLE]

It can be easily shown that

[TABLE]

where $\mathcal{L}=g_{\varphi\varphi}g_{tt}-g_{t\varphi}^{2}$.

II.1 Euler-Bernoulli equation

The way of deriving the Euler-Bernoulli equation (or the first integral of the Euler equations) for the ideal GRMHD energy-momentum tensor is described in Appendix A of shibata . The computation of the four-divergence

[TABLE]

yields

[TABLE]

The above equation is trivially satisfied for $\nu=t$ and $\nu=\varphi$. Nontrivial information is contained in Eq. (8) for $\nu=r,\theta$. Following shibata , one can show that

[TABLE]

and

[TABLE]

Combining the above expressions one gets

[TABLE]

or, dividing by $\rho h$,

[TABLE]

where we have used the fact that $dh=dp/\rho$. Therefore, it is possible to search for a solution in the form

[TABLE]

which is Eq. (A11) of shibata (note a misprint in the last term of the equation given in shibata ). The above equation is also equivalent to Eq. (11) of gf . We define the angular momentum per unit inertial mass $\rho h$ as $j=u^{t}u_{\varphi}$ and write

[TABLE]

This stays in agreement with the purely hydrodynamical case, where a functional relation $j=j(\Omega)$ (the rotation law) is an integrability condition of the Euler equations. Here we assume that $j=j(\Omega)$. It follows that also $b^{2}\mathcal{L}$ must be a function of $\rho h\mathcal{L}$. Further details of the Euler-Bernoulli equation, as well as the specific choices regarding the equation of state, the rotation law, and the prescription of the magnetic field will be discussed in Sec. II.3.

II.2 Einstein equations

Following shibata we derive the set of equations corresponding to a stationary black hole–torus spacetime from the standard $3+1$ formulation of the Einstein equations. The 3+1 metric reads

[TABLE]

where $\alpha$ is the lapse function, $\beta_{i}$ is the shift vector, and $\gamma_{ij}$ denote the components of the spatial metric. The vector normal to a surface of constant time $\Sigma_{t}$ is given by

[TABLE]

The Einstein constraint equations read

[TABLE]

and

[TABLE]

where $D_{i}$ and $R$ denote, respectively, the covariant derivative and the scalar curvature with respect to the metric $\gamma_{ij}$, induced on the slices $\Sigma_{t}$. The extrinsic curvature $K_{ij}$ of a slice $\Sigma_{t}$ is defined by

[TABLE]

where the Lie derivative of the three-metric is given by

[TABLE]

We denote $K=\gamma^{ij}K_{ij}$ and use a standard convention that spatial indices are raised and lowered using the induced metric $\gamma_{ij}$. The source terms $\rho_{\mathrm{H}}$ and $j^{i}$ are defined as

[TABLE]

where $P\indices{{}^{\mu}_{\nu}}=\delta\indices{{}^{\mu}_{\nu}}+n^{\mu}n_{\nu}$ is the spatial projection operator. The evolution equation for the extrinsic curvature $K_{ij}$ is

[TABLE]

Here $R_{ij}$ is the Ricci tensor with respect to the metric $\gamma_{ij}$. The tensor $S_{ij}$ is defined as

[TABLE]

and $S=\gamma^{ij}S_{ij}$.

We start by computing the source terms $\rho_{\mathrm{H}}$, $j_{i}$, $S_{ij}$ and $S$. Assuming the energy momentum-tensor (1), one gets

[TABLE]

Finally, the only nonvanishing component of $j_{i}$ is

[TABLE]

Note that there is no explicit magnetic contribution to $j_{\varphi}$. All these formulas can be obtained quite generally, assuming the metric of the form (3) and the conditions $u^{r}=u^{\theta}=b^{r}=b^{\theta}=0$.

As in shibata we assume from now on a metric in quasi-isotropic form:

[TABLE]

Thus, we need to provide equations for the four metric potentials appearing in Eq. (28), $\alpha$, $\beta$, $\psi$, and $q$, or, as we shall see, suitable combinations of these quantities (as in shibata ). In Eq. (28) $\beta_{\varphi}=\psi^{4}r^{2}\beta\sin^{2}\theta$, $\beta\equiv\beta^{\varphi}$, and $\alpha$ denotes the lapse function [as in Eq. (15)]. Since $\partial_{t}\gamma_{ij}=0$ and $\beta^{k}\partial_{k}\gamma_{ij}=0$, we get

[TABLE]

and $K=\gamma^{ik}K_{ik}=0$, i.e., we are in fact working in a maximal slicing. Therefore, the momentum constraint (17) can be written as

[TABLE]

The only nonvanishing components of $K_{ij}$ read

[TABLE]

To compute the momentum constraint we use a standard formula for symmetric tensors $K_{ij}$:

[TABLE]

where $\gamma=\det({\gamma_{ij}})$. Assuming the metric of the form (28) one obtains $\sqrt{\gamma}=\psi^{6}r^{2}\sin\theta e^{2q}$. The only nontrivial component of Eqs. (30) is obtained for $l=\varphi$,

[TABLE]

Following shibata the shift vector $\beta$ can be now split into two parts, $\beta=\beta_{\mathrm{K}}+\beta_{\mathrm{T}}$, where subindex K indicates the Kerr metric and subindex T the torus. We require that

[TABLE]

[TABLE]

and assume $H_{\mathrm{E}}$ and $H_{\mathrm{F}}$ corresponding to the Kerr metric (see Appendix A). More precisely, we choose $H_{\mathrm{E}}$ and $H_{\mathrm{F}}$ so that for the Kerr metric written in the form (28) one has

[TABLE]

These functions satisfy the momentum constraint of the form

[TABLE]

If a self-gravitating torus is present, we compute $\beta_{\mathrm{K}}$ from the relation

[TABLE]

which does not yield the Kerr form, as the conformal factor $\psi$ contains a contribution from the torus.

Inserting the expressions for $K_{r\varphi}$ and $K_{\theta\varphi}$ into Eq. (34) we obtain, after some algebra, an elliptic-type equation for $\beta_{\mathrm{T}}$

[TABLE]

where $\Delta$ denotes the flat Laplacian operator in spherical coordinates. Again, as in shibata we replace the lapse function $\alpha$ by the combination $\Phi=\alpha\psi$ and rewrite the previous equation as

[TABLE]

which is the same as Eq. (20) of shibata .

The equation for the conformal factor follows from the Hamiltonian constraint, Eq. (18), which for metric (28) reads

[TABLE]

It can be easily shown that

[TABLE]

where we use the short-hand notation

[TABLE]

The Ricci scalar can be computed as

[TABLE]

where

[TABLE]

This allows us to write the Hamiltonian constraint (18) in the form (Eq. (19) of shibata )

[TABLE]

The next equation follows from the general evolution equation for $K$. It can be obtained by computing the trace of Eq. (22):

[TABLE]

Using the Hamiltonian constraint, this equation can also be written as

[TABLE]

Since

[TABLE]

we obtain an elliptic-type equation for $\Phi$

[TABLE]

or, in terms of $\tilde{R}$,

[TABLE]

which is Eq. (18) of shibata .

The last equation (for the potential $q$) is obtained from Eq. (22). We define

[TABLE]

Consider the equation

[TABLE]

It yields, in particular, the term

[TABLE]

and the equation

[TABLE]

This allows us to compute $\tilde{R}$. The result can be combined with Eq. (48), yielding a new form of the elliptic equation for the conformal factor $\psi$

[TABLE]

which corresponds to Eq. (30) of shibata . A direct calculation then gives

[TABLE]

which generalizes Eq. (31) of shibata .

From the technical point of view, the black hole is introduced by specifying suitable boundary conditions. This can be done in an elegant manner by adapting the above equation to the “puncture” form (see brandt:1997 ; krivan:1998 ). Assuming that the puncture is located at $r=0$, we define

[TABLE]

where $r_{\mathrm{s}}=\frac{1}{2}\sqrt{m^{2}-a^{2}}$ is a radius of a coordinate sphere corresponding to the black hole horizon. Parameters $m$ and $a$ are chosen in such a way that in the Kerr spacetime with the asymptotic mass $m$ and the spin parameter $a$ the event horizon would be located at $r=r_{\mathrm{s}}$.

The above definitions lead to the following equations for the functions $\phi$, $B$, $\beta_{\mathrm{T}}$, and $q$:

[TABLE]

where

[TABLE]

which replace Eqs. (44–47) of shibata when a toroidal magnetic field is present in the disk. In the above formulas we denoted $B=e^{\tilde{b}}$. Equation (40) can be written as

[TABLE]

Notice that it does not yield the Kerr form, as there are contributions from the torus in both $B$ and $\phi$.

In our numerical approach we assume equatorial symmetry and solve equations (61), (62) and (63) in the domain defined by $r\in(r_{\mathrm{s}},r_{\infty})$, $\theta\in(0,\pi/2)$. Here $r_{\infty}$ is large, but finite.

The boundary conditions at $r=r_{\mathrm{s}}$ read

[TABLE]

It can be shown that Eq. (61d) requires a more stringent condition, which we set as $\beta_{\mathrm{T}}=O[(r-r_{\mathrm{s}})^{4}]$, or equivalently, $\beta_{\mathrm{T}}=\partial_{r}\beta_{\mathrm{T}}=\partial_{rr}\beta_{\mathrm{T}}=\partial_{rrr}\beta_{\mathrm{T}}=0$ at $r=r_{\mathrm{s}}$. In this choice, reflecting a freedom of fixing the splitting $\beta=\beta_{\mathrm{T}}+\beta_{\mathrm{K}}$, we follow shibata ; this choice has consequences in the definition of the angular momentum of the black hole (cf. Sec. III).

With the above boundary conditions, the two-surface $r=r_{\mathrm{s}}$ embedded in a hypersurface of constant time $\Sigma_{t}$ becomes a Marginally Outer Trapped Surface (MOTS) or the so-called apparent horizon. This can be easily demonstrated as follows. A MOTS is defined as a two-surface $\mathcal{S}$ embedded in $\Sigma_{t}$ on which the scalar expansion of the outgoing null geodesics

[TABLE]

vanishes. Here $H=D_{i}m^{i}$ denotes the mean curvature of the surface $\mathcal{S}$, and $m^{i}$ is a unit vector tangent to $\Sigma_{t}$ and normal to $\mathcal{S}$. For the two-surface $r=r_{\mathrm{s}}$, the components of the three-vector $m^{i}$ are given by $m^{i}=(m^{r},m^{\theta},m^{\varphi})=(\psi^{-2}e^{-q},0,0)$. Consequently, at $r=r_{\mathrm{s}}$,

[TABLE]

since both terms $K_{ij}m^{i}m^{j}$ and $K$ vanish. Using Eq. (60), one can show that

[TABLE]

at $r=r_{\mathrm{s}}$. It is now clear that the boundary conditions assumed at $r=r_{\mathrm{s}}$ imply that $\theta_{+}=H=0$. Note that the surface $r=r_{\mathrm{s}}$ is not only an apparent horizon, but it is also a minimal surface.

At the axis $\theta=0$ we assume regularity conditions $\partial_{\theta}\phi=\partial_{\theta}B=\partial_{\theta}\beta_{\mathrm{T}}=0$. Local flatness implies that $q=0$ at $\theta=0$. At the equator, we require symmetry conditions $\partial_{\theta}q=\partial_{\theta}\phi=\partial_{\theta}B=\partial_{\theta}\beta_{\mathrm{T}}=0$.

The asymptotic expansions of $q$, $\phi$, $B$, and $\beta_{\mathrm{T}}$ are discussed in Sec. III. They are used to impose boundary conditions at $r=r_{\infty}$. Further details on the numerical implementation of the boundary and asymptotic conditions can be found in kkmmop2 .

II.3 Details of the Euler-Bernoulli equation

We next discuss details of the Euler-Bernoulli equation, Eq. (14). The following three components have to be specified in order to obtain a solution: the equation of state, the rotation law $j=j(\Omega)$, and a prescription of the distribution of the magnetic field.

We assume the Keplerian rotation law derived in kkmmop ; kkmmop2 , i.e.,

[TABLE]

This is an exact formula that characterizes the motion of circular geodesics at the equatorial plane of the Kerr spacetime, in which case $w^{2}=m$, where $m$ is the Kerr mass. For self-gravitating tori $w^{2}\neq m$, in general. In the Newtonian limit rotation law (68) yields the standard Keplerian prescription of the angular velocity $\Omega=w/(r\sin\theta)^{\frac{3}{2}}$. It also agrees (for $a=0$) with the post-Newtonian Keplerian prescription proposed in mm2015 . The advantage of using this rotation law is that it allows one to obtain numerical solutions in a wide range of the parameters describing the torus kkmmop ; kkmmop2 .

The angular velocity $\Omega$ can be obtained by solving the relation $j(\Omega)=u^{t}u_{\varphi}$ for $\Omega$. In more explicit terms this relation reads

[TABLE]

or

[TABLE]

where $j(\Omega)$ is given by Eq. (68). We assume a convention with $\Omega>0$. This can correspond both to a torus corotating with the black hole, if $a>0$, or counterrotating, for $a<0$.

The specification of the magnetic term

[TABLE]

with $\mathcal{L}=-\alpha^{2}\psi^{4}r^{2}\sin^{2}\theta$, is somewhat more arbitrary, in the sense that there seem to be no physical “hints” concerning its prescription. Assuming a functional relation of the form $b^{2}|\mathcal{L}|=f(x)$, where $x=\rho h|\mathcal{L}|$ (note that this functional relation fulfills the general relativistic version of the von Zeipel condition for a purely toroidal magnetic field Zanotti15 ), we obtain

[TABLE]

Suppose that we would like to get

[TABLE]

where $C_{1}$ and $n$ are constants. This yields

[TABLE]

and a solution of the form

[TABLE]

For $x=0$ we get $f(x=0)=C_{2}$. Consequently, we set $C_{2}=0$. This ensures that the magnetic field vanishes for vanishing $\rho$. We have, finally

[TABLE]

We assume the above prescription of the magnetic field in this paper.

Equation (5) with the metric terms of Eq. (28) yields

[TABLE]

The Euler-Bernoulli Eq. (14) can be now written in the form

[TABLE]

In this paper we work with the polytropic equation of state of the form $p=K\rho^{\gamma}$. This yields the expression for the specific enthalpy

[TABLE]

Note that the magnetic distribution was chosen in such a way that in the limit $\rho\to 0$, the Euler-Bernoulli equation has the same form as in the absence of the magnetic field, i.e.,

[TABLE]

In our numerical procedure this form is used to establish the constants $w$ and $C^{\prime}$, assuming that the torus is characterized by some fixed equatorial coordinate radii $r_{1}$ and $r_{2}$.

An important numerical aspect concerns the specification of the polytropic constant $K$ of the equation of state. It is adjusted during the numerical iterative procedure so that the maximum value of the density $\rho$ within the torus is fixed at an a priori prescribed value (see Sec. IV for an exact discussion of this point).

We note that another possibility of setting up the details of the Euler-Bernoulli equation is to assume the rotation law of the form $\tilde{j}=hu_{\varphi}=\tilde{j}(\Omega)$. This is, for instance, the choice used in shibata . The Euler-Bernoulli equation is then given by Eq. (A10) of shibata . Such a formulation would suggest a different profile of $b^{2}$, given for instance by Eq. (A14) of shibata .

III Masses and angular momenta

The asymptotic Arnowitt-Deser-Misner mass can be computed as shibata

[TABLE]

where

[TABLE]

and $m$ is the Kerr mass. Defining the mass of the black hole is less straightforward. The central black hole is surrounded by a minimal two-surface located at $r=r_{\mathrm{s}}$ in the puncture method, on a fixed hypersurface of constant time. There is a collection of quantities that can be used to characterize the geometry of the horizon $r=r_{\mathrm{s}}$. The area of the horizon is given by

[TABLE]

where the integral is evaluated at $r=r_{\mathrm{s}}$. One can also define (at $r=r_{\mathrm{s}}$) $\Omega_{\mathrm{H}}=-\beta=-\beta_{\mathrm{K}}=\mathrm{const}$, and the surface gravity $\kappa=\partial_{r}\alpha\psi^{-2}e^{-q}=\mathrm{const}$. It can be easily shown that $\kappa=Be^{-4\phi-q}/8r_{\mathrm{s}}$. The angular momentum of the black hole can be defined as (see below)

[TABLE]

A mass defined at $r=r_{\mathrm{s}}$ as

[TABLE]

obeys the Smarr formula

[TABLE]

The mass of the black hole used in this paper is defined differently. Following shibata we adopt Christodoulou’s formula christodoulou . We define first the so-called irreducible mass

[TABLE]

Then the mass of the black hole is defined as

[TABLE]

From the ADM mass and the black hole mass we can define the torus mass as $M_{\mathrm{T}}=M_{\mathrm{ADM}}-M_{\mathrm{BH}}$, as was done in kkmmop ; kkmmop2 . There is, however, another possibility for the mass measure of the torus:

[TABLE]

It satisfies the relation

[TABLE]

We use this relation as a test of the accuracy of our numerical solutions. A direct computation yields

[TABLE]

In the numerical code, we compute the above quantity as

[TABLE]

Correspondingly, the angular momentum of the torus is defined as

[TABLE]

This is a standard definition corresponding to the Killing vector $\eta^{\mu}=(0,0,0,1)$, and the conservation law $\eta^{\nu}\nabla_{\mu}T\indices{{}^{\mu}_{\nu}}=\nabla_{\mu}(T\indices{{}^{\mu}_{\nu}}\eta^{\nu})=0$ leshouches . Note that $T\indices{{}^{t}_{\varphi}}=(\rho h+b^{2})u^{t}u_{\varphi}-b^{t}b_{\varphi}=\rho hu^{t}u_{\varphi}$, i.e., the contributions from the magnetic terms cancel. The total, asymptotic angular momentum reads $J=J_{\mathrm{H}}+J_{1}$.

The value of the angular momentum $J_{\mathrm{H}}$ depends on the assumed boundary conditions for $\beta_{\mathrm{T}}$. In our case $\beta_{\mathrm{T}}=\partial_{r}\beta_{\mathrm{T}}=\partial_{rr}\beta_{\mathrm{T}}=\partial_{rrr}\beta_{\mathrm{T}}=0$ at $r=r_{\mathrm{s}}$, and consequently $J_{\mathrm{H}}=am$.

The mass $M_{1}$ and the angular momentum $J_{1}$ are related to the asymptotic behaviour of the metric functions $\phi$ and $\beta_{\mathrm{T}}$, namely,

[TABLE]

as $r\to\infty$. The asymptotic behaviour of the two remaining functions $B$ and $q$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

We use the above asymptotic expansions to set the boundary conditions at the outer boundary of the numerical grid, i.e., at $r=r_{\infty}$.

IV Numerical method

To construct our models of self-gravitating, magnetised tori around rotating black holes we need to solve numerically the equations derived in Sec. II. The metric functions are described by Eq. (63) and Eqs. (61) with the source terms given by Eqs. (62). The angular velocity $\Omega$ must satisfy Eq. (70) with $j(\Omega)$ given by Eq. (68). The distribution of the enthalpy $h$, rest-mass density, and the pressure $p$ are obtained from Eq. (78) and from the polytropic relation (79). The quantity $A^{2}$ appearing in expressions (62) is defined by (45) where $K_{r\varphi}$ and $K_{\theta\varphi}$ should be computed according to formulas (35) and (36).

The numerical code used to obtain the solutions presented in this paper is a modification of the code described and tested in kkmmop2 to which the interested reader is addressed for details. It is an iterative method, where in each Newton-Raphson iteration one solves Eq. (70) for the angular velocity $\Omega$, Eq. (78) for the density $\rho$ (or the specific enthalpy $h$), and then Eqs. (61) for the metric functions. The latter are solved with 2nd-order finite differences. We take advantage of the banded matrix structure of the resulting linear equations and use LAPACK lapack . The changes introduced with respect to the version of the code described in kkmmop2 are only related to the presence of the magnetic field. While the inclusion of the magnetic terms in Eqs. (61) is straightforward, solving Eq. (78) with the magnetic terms is more troublesome. To describe it we need to discuss details connected with the treatment of Eqs. (70) and (78).

Each iteration is started with a Newton-Raphson procedure that gives the values of constants $w$ and $C^{\prime}$, assuming that the inner and outer equatorial radii of the disk ($r_{1}$ and $r_{2}$, respectively) are fixed. This procedure solves Eqs. (70) and (80) at points $(r,\theta)=(r_{1},\pi/2),(r_{2},\pi/2)$. These are four equations for the four unknowns $w$, $C^{\prime}$, $\Omega_{1}=\Omega(r_{1},\pi/2)$, $\Omega_{2}=\Omega(r_{2},\pi/2)$; in this step we assume that the metric functions are known from the previous iteration, or from the initial guess. In the next step we compute the values of $\Omega$ in a region which is large enough to contain the disk, but smaller than the domain covered by the numerical grid. In this way we can avoid problems with finding solutions to Eq. (70) in the vicinity of the symmetry axis $\theta=0$. Equation (70) is also solved with a Newton-Raphson scheme. The next stage consists in solving Eq. (78) for the specific enthalpy $h$, also by a Newton-Raphson procedure. The problem that one encounters here (which is absent in the purely hydrodynamical case) is that Eq. (78) contains a density term $\rho$, and in order to obtain a solution for $h$ one has to specify the value of the polytropic constant $K$. On the other hand, in kkmmop2 we found that the possibility of obtaining a convergent solution increases considerably if the solution is parameterized by a maximum value of the rest-mass density $\rho_{\mathrm{max}}$ within the disk. In the purely hydrodynamical case the value of the polytropic constant is then adjusted at each iteration so that the maximum value of the specific enthalpy $h$ obtained from Eq. (78) (with no magnetic terms) corresponds to the maximum of $\rho$ equal to an a priori prescribed value $\rho_{\mathrm{max}}$. This approach is not straightforward in the present GRMHD case. Therefore, we instead take the value of the polytropic constant $K$ inherited from the previous iteration, solve Eq. (78) for $h$, and then assume a value of $K$ so that the maximum in the specific enthalpy $h$ corresponds to the maximum in $\rho$ equal to $\rho_{\mathrm{max}}$. This approach leads to convergent solutions.

All stationary solutions of self-gravitating, magnetised disks obtained in this work have been computed on a numerical grid with approximately 800 nodes in the radial direction and 200 nodes in the angular direction. Specifically, the nodes in the grid are distributed according to

[TABLE]

in the radial direction, and

[TABLE]

where $\Delta\mu=1/(N_{\theta}-2)$, in the angular direction. We choose, in particular, $\Delta r=r_{\mathrm{s}}/50$, $f=1.01$, $N_{r}=800$, $N_{\theta}=200$. The above grid specification is similar to the one used in shibata .

The number of iterations required to obtain a solution depends mainly on the resolution of the grid, but also on the parameters of the solution kkmmop2 . Obtaining the solutions collected in Table 1 required typically $\sim 10^{4}$ to $\sim 2\times 10^{4}$ interations. Highly magnetised disks denoted in Table 1 as 4e–4i are exceptional, and required up to $\sim 10^{5}$ iterations.

V Results

The numerical solutions are specified by the following set of parameters: the black hole mass and angular momentum parameters, $m$ and $a$, the inner and outer radii of the disk $r_{1}$, $r_{2}$, the polytropic exponent of the equation of state $\gamma$, the maximum rest-mass density within the disk $\rho_{\mathrm{max}}$, and the constants $C_{1}$ and $n$ that characterize the prescription of the magnetic field. We note that this parameterization does not specify solutions uniquely. In fact, even in the case with no magnetic field one can observe a bifurcation: two solutions corresponding to different asymptotic masses can be obtained for fixed $m$, $a$, $r_{1}$, $r_{2}$, $\gamma$, and $\rho_{\mathrm{max}}$. Usually, one of these solutions corresponds to a case with the mass of the torus much larger than the mass of the central black hole kulczycki . This effect is interesting per se, and it will be the subject of a separate study.

The values $r_{1}$ and $r_{2}$ refer to coordinate radii of the torus. The simplest way of obtaining a geometrical size measure would be to use the circumferential radius $r_{\mathrm{C}}$ related to $r$ at the equatorial plane by the coordinate transformation $r_{\mathrm{C}}=\psi^{2}r$. In the following, by $r_{\mathrm{C,1}}$ and $r_{\mathrm{C,2}}$ we will denote the circumferential radii corresponding to coordinate radii $r_{1}$ and $r_{2}$, respectively. It should be kept in mind that in the strong gravitational-field regime the relation between $r$ and $r_{\mathrm{C}}$ may be not monotonic. In fact, numerical solutions representing self-gravitating tori with a maximum of $r_{\mathrm{C}}$ occurring within the torus, and not at its outer edge, were computed in labranche .

We measure the impact of the magnetic field by computing a magnetisation parameter $\beta_{\mathrm{mag}}$ defined as

[TABLE]

and evaluated at the maximum of the thermal pressure $p$.

We have computed a sample of numerical solutions. Their parameters are reported in Table 1. This table also provides the values of a handful of quantities that can be used to characterize the solutions: the total ADM mass $m_{\mathrm{ADM}}$, the mass of the black hole $m_{\mathrm{BH}}$, the angular momentum of the torus $J_{1}$, the circumferential inner and outer radii of the torus $r_{\mathrm{C,1}}$, $r_{\mathrm{C,2}}$, and the magnetisation parameter $\beta_{\mathrm{mag}}$. In all our numerical examples we set $m=1$. This can be viewed as fixing the system of units; in practice, setting $m=1$ yields $m_{\mathrm{BH}}\approx 1$. Table 1 is divided into parts which group models with the same values of $a$, $r_{1}$, $r_{2}$, $\rho_{\mathrm{max}}$ but different values of $C_{1}$. In particular, each group corresponds to a different value of $a$; we chose specifically $a=-0.5$, [math], $0.9$, and $0.99$. A negative value of $a$ means counterrotation, i.e. we adhere to a convention with $\Omega>0$. For simplicity, we fix in all our solutions the value of the polytropic exponent to $\gamma=4/3$ and the parameter $n$ of the magnetisation law to $n=1$. Except for the models with the fastest spinning black hole ($a=0.99$), in all investigated cases we were able to obtain solutions with the magnetisation parameter $\beta_{\mathrm{mag}}$ ranging from $\infty$ (no magnetic field) to the level of the order of $10^{-3}$ to $10^{-4}$ (i.e. highly magnetised models). The case with ($a=0.99$) listed in Table 1 is exceptional: the tori characterized by small values of $\beta_{\mathrm{mag}}$ are fairly light. Moreover, a large number of numerical iterations ($\sim 10^{5}$) is required in order to converge to a solution.

It is not unreasonable to assume that stable solutions should have $r_{1}$ larger than the location of the Innermost Stable Circular Orbit (ISCO). Because of the self-gravity of the torus, the location of the ISCO deviates from the value characteristic for the Kerr spacetime with a given mass $m$ and spin parameter $a$. Nevertheless, the Kerr values can be still used to get a rough estimate of the location of the ISCO. For $m=1$ and $a=-0.5$, the circumferential radius of the ISCO in the Kerr spacetime is $r_{\mathrm{C,ISCO}}=7.57$. For $m=1$ and $a=0,0.9$ and $0.99$, we obtain, respectively, $r_{\mathrm{C,ISCO}}=6,2.63$, and $2.11$. Of course, given a numerical solution, the true location of the ISCO can also be computed, for instance using the formalism described in sasaki . We have actually checked that the solutions listed in Table 1 satisfy the condition $r_{\mathrm{C,1}}>r_{\mathrm{C,ISCO}}$. A detailed analysis of the influence of the self-gravitating torus on the location of the ISCO will be given elsewhere.

From the inspection of all solutions listed in Table 1 we conclude that the configurations with smaller values of $\beta_{\mathrm{mag}}$ (relatively stronger magnetic fields) tend to have the maxima of the density shifted towards smaller radii. This behavior is illustrated in Fig. 1, which depicts the morphology of a subset of models by plotting the logarithm of their rest-mass density $\rho$ in the meridional half-plane. Figure 2, in which we plot radial profiles of the density at the equatorial plane, shows this trend in a more clear way. For clarity, Fig. 2 displays the profiles both in linear and logarithmic scales (the latter in the insets). We note that the shift of the maximum of the density towards smaller radii in magnetised disks has already been observed in the test-fluid models built in gf .

On the other hand, the analysis of the solutions with different values of the constant $C_{1}$ appearing in Eq. (78) shows that the larger $C_{1}$ the smaller the thermal pressure $p$, even in cases in which the maximum of the baryonic density $\rho$ is fixed. In general we also observe an increase of the absolute values of the magnetic pressure $p_{\mathrm{mag}}=\frac{1}{2}b^{2}$. Both factors lead to a rapid decrease of the magnetisation parameter $\beta_{\mathrm{mag}}$. In some sense, with an increase of $C_{1}$, the role of the (gradient of the) thermal pressure in counter-balancing gravity is taken over by the magnetic pressure. This effect is illustrated in Fig. 3 for models 2c, 2e, 3c and 3e. Note that in both cases the maximum of the largest of the two pressures (thermal or magnetic) is about two orders of magnitude smaller than the maximum of the rest-mass density.

The presence of a magnetic field affects the total ADM mass of the system (mainly by influencing the mass of the torus) in a nontrivial way. Although a direct contribution of the terms related with the magnetic field to the mass, as computed e.g. from Eq. (89), is small, the magnetic field changes the total mass of the system by affecting the distribution of the rest-mass density within the disk. The nature of the changes of the ADM mass with an increasing magnetic field contribution depends on the spin of the black hole and on the distance between the black hole and the torus (mainly on the location of the inner edge of the torus $r_{1}$). Disentangling these two factors is difficult, since the location of the ISCO depends predominantly on the spin of the black hole, and we want our models to satisfy the condition $r_{\mathrm{C,ISCO}}<r_{\mathrm{C,1}}$. The dependence of the ADM mass with the black hole spin in the examples collected in Table 1 is as follows: for $a=-0.5$ the ADM mass grows with increasing $C_{1}$, for $a=0$ and $a=0.9$ the behaviour of the ADM mass is not monotonic with $C_{1}$, and for $a=0.99$ the ADM mass decreases with increasing $C_{1}$. In the latter case this effect is strong. The ADM mass drops from $m_{\mathrm{ADM}}=1.70$ for $C_{1}=0$ (no magnetic field) to $m_{\mathrm{ADM}}=1.05$ for $C_{1}=4.7$ (magnetisation parameter $\beta_{\mathrm{mag}}=2.28\times 10^{-2}$).

The importance of the effects of the disk self-gravity can be estimated by plotting the deviation of the lapse function at the equatorial plane with respect to its value for an isolated Kerr black hole, $(\alpha_{\mathrm{K}}-\alpha)/\alpha$ (here $\alpha_{\mathrm{K}}$ is computed using Eq. (107)). These radial profiles are plotted in Fig. 4 for the same subset of models depicted in Fig. 1. We observe a correlation of the maximum deviations with the rest-mass density $\rho_{\mathrm{max}}$ of the models, the deviations being larger the larger the value of $\rho_{\mathrm{max}}$, namely $~{}20\%$ for model 4a, $~{}8.5\%$ for model 3d and $~{}3.5\%$ for models 1d and 2d. Furthermore, it can be seen that, as expected, the deviation grows if the fraction of the mass stored at the torus is greater. This fraction can be easily inferred from Table 1.

A close inspection of plots in Figs. 2 and 4 reveals a somewhat unexpected similarity of certain features of models 1a–1f and 2a–2f. These two families of models are characterized by the same values of coordinate radii $r_{1}$ and $r_{2}$, and the same maximal densities $\rho_{\mathrm{max}}$, but they differ in the assumed values of the black hole spin parameter ($a=-0.5$ and $a=0$, respectively). The plots of the rest-mass density shown in Fig. 2 and corresponding to models belonging to both families are nearly indistinguishable from one another. One can also hardly spot any difference between models 1a–1f and 2a–2f in the plots of the differences between the lapse functions (the actual lapse and the lapse corresponding to the Kerr metric) shown in Fig. 4. However, both classes of models are actually different. The differences stand out clearly in the plots of the shift vector $\beta$ shown in Fig. 5. Of course, the absolute values of the lapse function characterizing the models belonging to both classes are also different.

VI Summary

We have presented general-relativistic models of stationary, axisymmetric, self-gravitating, magnetised disks (or tori) rotating around spinning black holes. They have been obtained by solving numerically the coupled system of the Einstein equations and the equations of ideal GRMHD. The mathematical formulation of our approach has closely followed the work of Shibata shibata , who built purely hydrodynamical self-gravitating, constant angular momentum tori around black holes in the puncture framework. The inclusion of magnetic fields represents the first new ingredient of our approach. On the other hand, building on previous studies of configurations with no magnetic field kkmmop ; kkmmop2 ; kmm , we have constructed our magnetised models assuming a Keplerian rotation law in the disks, which departs from the constant angular momentum disks reported by shibata . The use of the Keplerian rotation law is the second new ingredient of our procedure. We have focused on toroidal distributions of the magnetic field and presented a large set of models corresponding to a wide range of values of the magnetisation parameter, starting with weakly magnetised disks and ending at configurations in which the magnetic pressure dominates over the thermal one.

The impact of the magnetic field on the disk structure is mainly related to the magnetic pressure. In all investigated models, we have observed a shift of the location of the maximum of the rest-mass density towards the central black hole. The impact of the magnetisation on the total mass of the system (or the mass of the disk) depends on the black hole spin and on the geometry of the disk. It is possible to obtain classes of solutions in which the mass of the torus decreases with the decreasing magnetisation parameter, but the converse can also be true.

All our solutions have been obtained for the polytropic equation of state with the polytropic exponent $\gamma=4/3$, and for a specific choice of the magnetisation law. These choices can, of course, have an impact on the obtained results. Furthermore, the values of the ratio of the black-hole mass to the total mass of the system reported in this paper were kept within a reasonable range: the obtained disks are massive enough so that the effects connected with the self-gravity become important. On the other hand, in this work we have not considered disks with masses exceeding the mass of the central black hole. It is known that allowing for sufficiently large disk masses can lead to the occurrence of several general-relativistic effects, characteristic of the strong gravitational-field regime. In ansorg:2006 Ansorg and Petroff showed that a perfect fluid torus rotating (rigidly) around a black hole can create its own ergoregion, disconnected from the ergoregion of the black hole. (This effect is also present in more exotic objects, for instance the scalar hairy black holes described in Herdeiro15 , for which the scalar field has a toroidal distribution.) In labranche Labranche, Petroff, and Ansorg gave examples of perfect fluid tori (with no central object) in which the circumferential radius attains its maximum inside the torus, and not at its outer edge. We expect these effects to be present also within our formulation.

The results presented in this paper can be extended in several directions. On the one hand we plan to investigate the influence of the self-gravity of the torus on the location of the ISCO of a rotating black hole. In addition, we will also study the non-linear stability properties of our solutions through numerical-relativity simulations in a dynamical spacetime setup. There are a number of instabilities that may affect the disks, such as the runaway, the Papaloizou-Pringle and the magneto-rotational instabilities. In particular, the development of the PPI in tori threaded by toroidal magnetic fields may be significantly affected by the concurrent development of the MRI, as shown recently by bugli for non-self-gravitating disks. The self-gravitating, magnetised tori we have built in this work can be used to investigate the generality of those findings beyond the test-fluid limit.

The presented solutions, and models of self-gravitating magnetised disks in general, should be also relevant to the ongoing attempts to estimate the amount of angular momentum within a given volume. This is an interesting area of research within mathematical relativity dain ; khuri . Recent works focusing on such estimates for stationary Keplerian self-gravitating disks around black holes include kmm ; kmmpx .

Acknowledgements.

PM was partially supported by the Polish National Science Centre grant No. 2017/26/A/ST2/00530. JAF acknowledges support from the Spanish MINECO (grant AYA2015-66899-C2-1-P) and from the European Union’s Horizon 2020 RISE programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. MP was supported in part by “PHAROS”, COST Action CA16214; LOEWE-Program in HIC for FAIR.

Appendix A Kerr metric in quasi-isotropic coordinates

For completeness, we express the Kerr metric in the quasi-isotropic coordinates of the form (28) shibata ; brandtseidel .

Define

[TABLE]

The Kerr metric can be expressed as

[TABLE]

where

[TABLE]

The functions $H_{\mathrm{E}}$ and $H_{\mathrm{F}}$ corresponding to the Kerr metric read

[TABLE]