© ESO, 2017

1. Introduction

Matter accretion on to black holes is the most efficient form of energy production known in nature. The conversion of the gravitational energy of the infalling matter into heat and radiation may reach efficiencies of about 43% in the case of maximally rotating (Kerr) black holes. For this reason, systems formed by a black hole surrounded by an accretion disc are deemed responsible for many of the most energetic astronomical phenomena observed in the cosmos. In particular, (geometrically) thick accretion disks (or tori) are believed to be present in quasars and other active galactic nuclei, some X-ray binaries and microquasars, as well as in the central engine of gamma-ray bursts. The latter are alluded to in connection with mergers of neutron star binaries and black hole neutron star binaries, as well as with the rotational collapse that ensues at the end of the life of some massive stars. As numerical simulations show, such events often result in a black hole surrounded by a torus (see e.g. Rezzolla et al. 2010; Sekiguchi & Shibata 2011; Faber & Rasio 2012; Shibata & Taniguchi 2011; Baiotti & Rezzolla 2017).

Investigation of such systems, either by analytical or numerical means, may rely on the ability to construct suitable representations based on physical assumptions. The construction of equilibrium models of stationary disks around black holes has indeed a long tradition (see Abramowicz & Fragile 2013, and references therein). In particular, the “Polish doughnuts” (Abramowicz et al. 1978; Kozlowski et al. 1978) provide a very general method to build equilibrium configurations of perfect fluid matter orbiting around a Kerr black hole. This fully relativistic model assumes that the disc is non-magnetised and that the matter obeys a constant specific angular momentum distribution. This method was later extended by Komissarov (2006) by adding a purely azimuthal magnetic field to build magnetised tori around rotating black holes. Dynamical evolutions of magnetised tori built with the Komissarov solution were first reported by Montero et al. (2007). On the other hand, assuming different distributions of angular momentum in the disks, Qian et al. (2009) presented a method to build sequences of black hole thick accretion disks in dynamical equilibria, restricted however to the purely hydrodynamical case.

In this paper we have combined the two approaches considered in Komissarov (2006) and Qian et al. (2009) to build new sequences of equilibrium tori around Kerr black holes. Building on these works, we present here the extension of the models of Qian et al. (2009) to account for disks endowed with a purely toroidal magnetic field. In our procedure we hence assumed a form of the angular momentum distribution that departs from the constant case considered by Komissarov (2006) and from which the location and morphology of the equipotential surfaces can be numerically computed. As we shall show below, for the particular case of constant angular momentum distributions, our method is in good agreement with the results of Komissarov (2006). We also note the recent work of Wielgus et al. (2015) where Komissarov’s solution was extended for the particular case of power-law distributions of angular momentum. Moreover, the magnetised tori of Wielgus et al. (2015) were used to explore the growth of the magneto-rotational instability (MRI) through time-dependent numerical simulations. In particular, the long-term evolution of those tori has been recently investigated by Fragile & Sa¸dowski (2017), who paid special attention to the decay of their magnetisation.

The organization of the paper is as follows: Sect. 2 presents the analytic framework to build the disks while Sect. 3 explains the corresponding numerical procedure. The sequences of models are discussed in Sect. 4. Finally, the conclusions are summarised in Sect. 5, where we also briefly indicate potential extensions of this work we plan to inspect. In our work we assume that the test-fluid approximation holds, thus neglecting the self-gravity of the fluid, and further assume that the spacetime is described by the Kerr metric. In mathematical expressions below Greek indices are spacetime indices running from 0 to 3 and Latin indices are only spatial. We use geometrised units where G = c = 1.

2. Framework

Equilibrium tori around Kerr black holes were built assuming that the spacetime gravitational potentials and the fluid fields are stationary and axisymmetric. In all the derivations presented below, the Kerr metric is implicitly written using standard Boyer-Lindquist coordinates. It is convenient to introduce a number of relevant characteristic radii, such as the radius of the marginally stable circular orbit, r_ms, and the radius of the marginally bound circular orbit, r_mb, given by $\begin{array}{ccc}{\mathit{r}}_{\mathrm{ms}}& \mathrm{=}& \mathit{M}\hspace{0.17em}{\left(\mathrm{3}\mathrm{+}{\mathit{Z}}_{\mathrm{2}}\mathrm{-}\mathrm{\left[}\mathrm{\right(}\mathrm{3}\mathrm{-}{\mathit{Z}}_{\mathrm{1}}\mathrm{\left)}\mathrm{\right(}\mathrm{3}\mathrm{+}{\mathit{Z}}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{Z}}_{\mathrm{2}}\mathrm{)}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\mathit{,}\\ {\mathit{r}}_{\mathrm{mb}}& \mathrm{=}& \mathrm{2}\mathit{M}\hspace{0.17em}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{\ast }}}{\mathrm{2}}\mathrm{+}\sqrt{\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathrm{\ast }}}\right)\hspace{0.17em}\mathit{,}\end{array}$where we have defined the following quantities, ${\mathit{Z}}_{\mathrm{1}}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathrm{\ast }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{[}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{a}}_{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{]}$, ${\mathit{Z}}_{\mathrm{2}}\mathrm{=}\mathrm{(}\mathrm{3}{\mathit{a}}_{\mathrm{\ast }}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, and a_∗ = a/M, with a and M being the spin Kerr parameter and the black hole mass, respectively.

2.1. Distribution of angular momentum

We introduce the specific angular momentum l and the angular velocity Ω employing the standard definitions, $\mathit{l}\mathrm{=}\mathrm{-}\frac{{\mathit{u}}_{\mathit{\phi }}}{{\mathit{u}}_{\mathit{t}}}\mathit{,}\mathrm{\Omega }\mathrm{=}\frac{{\mathit{u}}^{\mathit{\phi }}}{{\mathit{u}}^{\mathit{t}}}\mathit{,}$(3)where u^μ is the fluid four-velocity. The relationship between l and Ω is given by the equations $\mathit{l}\mathrm{=}\mathrm{-}\frac{\mathrm{\Omega }{\mathit{g}}_{\mathit{\phi \phi }}\mathrm{+}{\mathit{g}}_{\mathit{t\phi }}}{\mathrm{\Omega }{\mathit{g}}_{\mathit{t\phi }}\mathrm{+}{\mathit{g}}_{\mathit{tt}}}\mathit{,}\mathrm{\Omega }\mathrm{=}\mathrm{-}\frac{\mathit{l}{\mathit{g}}_{\mathit{tt}}\mathrm{+}{\mathit{g}}_{\mathit{t\phi }}}{\mathit{l}{\mathit{g}}_{\mathit{t\phi }}\mathrm{+}{\mathit{g}}_{\mathit{\phi \phi }}}\mathit{,}$(4)where g_μν is the metric tensor and we are assuming circular motion, that is, the four-velocity can be written as ${\mathit{u}}^{\mathit{\mu }}\mathrm{=}\mathrm{\left(}{\mathit{u}}^{\mathit{t}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}{\mathit{u}}^{\mathit{\phi }}\mathrm{\right)}\hspace{0.17em}\mathit{.}$(5)We also introduce the Keplerian angular momentum (for prograde motion) in the equatorial plane l_K, defined as ${\mathit{l}}_{\mathrm{K}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{a}{\mathit{M}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}{{\mathit{r}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{M}{\mathit{r}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{+}\mathit{a}{\mathit{M}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}$(6)Jaroszynski et al. (1980) argued that the slope of the specific angular momentum should range between two limiting cases, namely l = const. and Ω = const. Following Qian et al. (2009) we assume an angular momentum distribution ansatz given by $\mathit{l}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ {\mathit{l}}_{\mathrm{0}}{\mathrm{(}{\frac{{\mathit{l}}_{\mathrm{K}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{l}}_{\mathrm{0}}}}^{\mathrm{)}}}^{\mathit{\beta }}{\sin }^{\mathrm{2}\mathit{\gamma }}\mathit{\theta }& \mathrm{for}\mathit{r}\mathrm{\ge }{\mathit{r}}_{\mathrm{ms}}\\ {\mathit{l}}_{\mathrm{ms}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\sin }^{\mathrm{2}\mathit{\gamma }}\mathit{\theta }& \mathrm{for}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{ms}}\mathit{,}\end{array}$(7)where constants l₀ and l_ms(r) are defined by l₀ ≡ ηl_K(r_ms) and l_ms(r) ≡ l₀ [l_K(r_ms) /l₀] ^β. Therefore, the model for the distribution of angular momentum has three free parameters, β, γ and η, whose range of variation is given by (Qian et al. 2009) $\mathrm{0}\mathrm{\le }\mathit{\beta }\mathrm{\le }\mathrm{1}\mathit{,} \mathrm{-}\mathrm{1}\mathrm{\le }\mathit{\gamma }\mathrm{\le }\mathrm{1}\mathit{,} \mathrm{1}\mathrm{\le }\mathit{\eta }\mathrm{\le }{\mathit{\eta }}_{\max }\mathit{,}$(8)with η_max = l_K(r_mb) /l_K(r_ms). In this paper, and as it is done for hydrodynamical disks in Qian et al. (2009), we choose η = η_max, and then we can write l₀ as l₀ = l_K(r_mb). For this choice of η, we can find the location of the cusp of the disc within the range r_mb ≤ r_cusp ≤ r_ms (for 0 ≤ β ≤ 1). The cusp is defined as the circle in the equatorial plane on which the pressure gradient vanishes and the angular momentum of the disc equals the Keplerian angular momentum. Moreover, this value of η guarantees that constant angular momentum disks (β = γ = 0) with their inner edge located at the cusp (r_in = r_cusp) have no outer boundary, i.e they are infinite disks (see also Font & Daigne (2002)).

2.2. Magnetised disks

The equations of ideal general relativistic MHD are the following conservation laws, ∇_μT^μν = 0, ∇_μ^∗F^μν = 0, and ∇_μ(ρu^μ) = 0, where ∇_μ is the covariant derivative and ${\mathit{T}}^{\mathit{\mu \nu }}\mathrm{=}\mathrm{(}\mathit{w}\mathrm{+}{\mathit{b}}^{\mathrm{2}}\mathrm{)}{\mathit{u}}^{\mathit{\mu }}{\mathit{u}}^{\mathit{\nu }}\mathrm{+}\left(\mathit{p}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{b}}^{\mathrm{2}}\right){\mathit{g}}^{\mathit{\mu \nu }}\mathrm{-}{\mathit{b}}^{\mathit{\mu }}{\mathit{b}}^{\mathit{\nu }}\mathit{,}$(9)is the energy-momentum tensor of a magnetised perfect fluid, w and p being the fluid enthalpy density and fluid pressure, respectively. Moreover, ${}^{\mathrm{\ast }}\mathit{F}^{\mathit{\mu \nu }}\mathrm{=}{\mathit{b}}^{\mathit{\mu }}{\mathit{u}}^{\mathit{\nu }}\mathrm{-}{\mathit{b}}^{\mathit{\nu }}{\mathit{u}}^{\mathit{\mu }}$ is the (dual of the) Faraday tensor relative to an observer with four-velocity u^μ, and b^μ is the magnetic field in that frame, with b² = b^μb_μ. Assuming the magnetic field is purely azimuthal, that is b^r = b^θ = 0, and taking into account that the flow is stationary and axisymmetric, the conservation of the current density and of the Faraday tensor follow. Contracting the conservation law for the energy-momentum tensor with the projection tensor ${\mathit{h}}_{\hspace{0.17em}\hspace{0.17em}\mathit{\beta }}^{\mathit{\alpha }}\mathrm{=}{\mathit{\delta }}_{\hspace{0.17em}\hspace{0.17em}\mathit{\beta }}^{\mathit{\alpha }}\mathrm{+}{\mathit{u}}^{\mathit{\alpha }}{\mathit{u}}_{\mathit{\beta }}$, we arrive at $\mathrm{(}\mathit{w}\mathrm{+}{\mathit{b}}^{\mathrm{2}}\mathrm{)}{\mathit{u}}_{\mathit{\nu }}{\mathit{\partial }}_{\mathit{i}}{\mathit{u}}^{\mathit{\nu }}\mathrm{+}{\mathit{\partial }}_{\mathit{i}}\left(\mathit{p}\mathrm{+}\frac{{\mathit{b}}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{-}{\mathit{b}}_{\mathit{\nu }}{\mathit{\partial }}_{\mathit{i}}{\mathit{b}}^{\mathit{\nu }}\mathrm{=}\mathrm{0}\mathit{,}$(10)where i = r,θ. Following Komissarov (2006) we rewrite this equation in terms of the specific angular momentum l and of the angular velocity Ω, to obtain ${\mathit{\partial }}_{\mathit{i}}\mathrm{\left(}\ln {\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{)}\mathrm{-}\frac{\mathrm{\Omega }{\mathit{\partial }}_{\mathit{i}}\mathit{l}}{\mathrm{1}\mathrm{-}\mathit{l}\mathrm{\Omega }}\mathrm{+}\frac{{\mathit{\partial }}_{\mathit{i}}\mathit{p}}{\mathit{w}}\mathrm{+}\frac{{\mathit{\partial }}_{\mathit{i}}\mathrm{\left(}\mathrm{ℒ}{\mathit{b}}^{\mathrm{2}}\mathrm{\right)}}{\mathrm{2}\mathrm{ℒ}\mathit{w}}\mathrm{=}\mathrm{0}\mathit{,}$(11)where $\mathrm{ℒ}\mathrm{=}{\mathit{g}}_{\mathit{t\phi }}^{\mathrm{2}}\mathrm{-}{\mathit{g}}_{\mathit{tt}}{\mathit{g}}_{\mathit{\phi \phi }}$. To integrate Eq. (11)we first assume a barotropic equation of state w = w(p) of the form $\mathit{p}\mathrm{=}\mathit{K}{\mathit{w}}^{\mathit{\kappa }}\mathit{,}$(12)with K and κ constants. Then, we define the magnetic pressure as p_m = b²/ 2, and introduce the definitions $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{p}}_{\mathrm{m}}\end{array}\mathrm{=}\mathrm{ℒ}{\mathit{p}}_{\mathrm{m}}$ and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{w}\end{array}\mathrm{=}\mathrm{ℒ}\mathit{w}$, in order to write an analogue equation to Eq. (12)for $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{p}}_{\mathrm{m}}\end{array}$ (Komissarov 2006) $\begin{array}{c}\mathrm{˜}\\ {\mathit{p}}_{\mathrm{m}}\end{array}\mathrm{=}{\mathit{K}}_{\mathrm{m}}\stackrel{}{{{\mathit{w}}_{\mathrm{m}}^{\mathit{\lambda }}}_{\mathrm{˜}}}\mathit{,}$(13)or, in terms of the magnetic pressure p_m${\mathit{p}}_{\mathrm{m}}\mathrm{=}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{w}}^{\mathit{\lambda }}\mathit{,}$(14)where K_m and λ are constants. This particular choice of barotropic relationships, w = w(p) and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{w}\end{array}\mathrm{=}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{w}\end{array}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{p}}_{\mathrm{m}}\end{array}\mathrm{\right)}$, fulfill the general relativistic version of the von Zeipel theorem for a toroidal magnetic field (von Zeipel 1924; Zanotti & Pugliese 2015), that is, the surfaces of constant Ω and constant l coincide.

We can now integrate Eq. (11)to obtain $\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{l}}\frac{\mathrm{\Omega d}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{p}}\frac{\mathrm{d}\mathit{p}}{\mathit{w}}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\begin{array}{c}\mathrm{˜}\\ {\mathit{p}}_{\mathrm{m}}\end{array}}\frac{\mathrm{d}\begin{array}{c}\mathrm{˜}\\ {\mathit{p}}_{\mathrm{m}}\end{array}}{\begin{array}{c}\mathrm{˜}\\ \mathit{w}\end{array}}\mathrm{=}\mathrm{const}\mathit{.}$(15)On the surface of the disc, and particularly on its inner edge, the conditions $\mathit{p}\mathrm{=}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{p}}_{\mathrm{m}}\end{array}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{u}}_{\mathit{t}}\mathrm{=}{\mathit{u}}_{\mathit{t,}\mathrm{in}}\mathit{,l}\mathrm{=}{\mathit{l}}_{\mathrm{in}}$ are satisfied and, therefore, the integration constant is simply given by $\mathrm{const}\mathit{.}\mathrm{=}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}{\mathrm{\int }}_{{\mathit{l}}_{\mathrm{in}}}^{\mathit{l}}\frac{\mathrm{\Omega d}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\mathrm{·}$(16)We can also introduce the total (gravitational plus centrifugal) potential W (Abramowicz et al. 1978) and write the integral form of the equation of motion (the relativistic Euler equation) as $\mathit{W}\mathrm{-}{\mathit{W}}_{\mathrm{in}}\mathrm{=}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t,}\mathrm{in}}\mathrm{\right|}\mathrm{-}{\mathrm{\int }}_{{\mathit{l}}_{\mathrm{in}}}^{\mathit{l}}\frac{\mathrm{\Omega d}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\mathit{,}$(17)where W_in is the potential at the inner edge of the disc (in the equatorial plane). With this definition, we can write Eq. (15)as $\mathit{W}\mathrm{-}{\mathit{W}}_{\mathrm{in}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{p}}\frac{\mathrm{d}\mathit{p}}{\mathit{w}}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\begin{array}{c}\mathrm{˜}\\ {\mathit{p}}_{\mathrm{m}}\end{array}}\frac{\mathrm{d}\begin{array}{c}\mathrm{˜}\\ {\mathit{p}}_{\mathrm{m}}\end{array}}{\begin{array}{c}\mathrm{˜}\\ \mathit{w}\end{array}}\mathit{,}$(18)which for a barotropic equation of state can be easily integrated to give $\mathit{W}\mathrm{-}{\mathit{W}}_{\mathrm{in}}\mathrm{+}\frac{\mathit{\kappa }}{\mathit{\kappa }\mathrm{-}\mathrm{1}}\frac{\mathit{p}}{\mathit{w}}\mathrm{+}\frac{\mathit{\lambda }}{\mathit{\lambda }\mathrm{-}\mathrm{1}}\frac{{\mathit{p}}_{\mathrm{m}}}{\mathit{w}}\mathrm{=}\mathrm{0.}$(19)Replacing p and p_m by Eqs. (12)and (14), the previous equation reduces to $\mathit{W}\mathrm{-}{\mathit{W}}_{\mathrm{in}}\mathrm{+}\frac{\mathit{\kappa }}{\mathit{\kappa }\mathrm{-}\mathrm{1}}\mathit{K}{\mathit{w}}^{\mathit{\kappa }\mathrm{-}\mathrm{1}}\mathrm{+}\frac{\mathit{\lambda }}{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{K}}_{\mathrm{m}}\mathrm{(}\mathrm{ℒ}\mathit{w}{\mathrm{)}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{0}\mathit{,}$(20)which relates the distribution of the potential with the distribution of the enthalpy density.

3. Method

To construct our models of magnetised disks we have followed the approach described in Qian et al. (2009). First, we find the radial location of the cusp and of the centre of the disc in the equatorial plane, r_cusp and r_c, defined as the solutions to the equation l(r)−l_K = 0. Next, we compute the partial derivatives of the potential, Eq. (17)${\mathit{\partial }}_{\mathit{r}}\mathit{W}\mathrm{=}{\mathit{\partial }}_{\mathit{r}}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}\frac{\mathrm{\Omega }{\mathit{\partial }}_{\mathit{r}}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\mathit{,}$(21)and ${\mathit{\partial }}_{\mathit{\theta }}\mathit{W}\mathrm{=}{\mathit{\partial }}_{\mathit{\theta }}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}\frac{\mathrm{\Omega }{\mathit{\partial }}_{\mathit{\theta }}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\mathrm{·}$(22)Then, we integrate the radial partial derivative of the potential along the segment [r_cusp,r_c] (assuming W_cusp = 0) at the equatorial plane, thus obtaining the equatorial distribution of the potential between r_cusp and r_c${\mathit{W}}_{\mathrm{eq}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{cusp}}}^{{\mathit{r}}_{\mathrm{c}}}\left({\mathit{\partial }}_{\mathit{r}}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}\frac{\mathrm{\Omega }{\mathit{\partial }}_{\mathit{r}}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\right)\mathrm{·}$(23)Following Qian et al. (2009; see also Jaroszynski et al. 1980) we can divide Eqs. (21)and (22)to obtain $\mathit{F}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{{\mathit{\partial }}_{\mathit{r}}\mathit{W}}{{\mathit{\partial }}_{\mathit{\theta }}\mathit{W}}\mathrm{=}\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}\mathit{r}}\mathit{,}$(24)where function F(r,θ) is known in closed form for the Kerr metric once an assumption for the angular momentum distribution has been made (cf. Eq. (7)). For this reason, Eq. (24) also takes the form of an ordinary differential equation for the surfaces of constant potential, θ = θ(r), which, upon integration, yields the location of those surfaces. We note that in Qian et al. (2009) these are surfaces of constant fluid pressure instead, since their disks are purely hydrodynamical.

Fig. 1 Comparison with Komissarov’s solution for a constant angular momentum model. The left panel shows the rest-mass density distribution in logarithm scale for all radii and all angles while the middle and right panels show, respectively, the angular profile of the logarithm of the rest-mass density at the centre of the disc and the radial profile of the logarithm of the rest-mass density at the equatorial plane.

Next, we choose all the initial radial values for the integration of Eq. (24)to lie between r_cusp and r_c (θ = π/ 2). Since we are only interested in the equipotential surfaces inside the Roche lobe of the disc, our choice of initial values provides us a mapping of the equipotential surfaces of the torus. Given that we have already obtained both the equipotential surfaces θ(r) which cross the segment [r_cusp,r_c] at the equatorial plane and the values of the potential in that segment, we can obtain the complete potential distribution for the torus (outside of that segment). Once we have the potential distribution, we can find the fluid pressure at the centre of the disc from Eq. (20), ${\mathit{p}}_{\mathrm{c}}\mathrm{=}{\mathit{w}}_{\mathrm{c}}\mathrm{(}{\mathit{W}}_{\mathrm{in}}\mathrm{-}{\mathit{W}}_{\mathrm{c}}\mathrm{)}{\left(\frac{\mathit{\kappa }}{\mathit{\kappa }\mathrm{-}\mathrm{1}}\mathrm{+}\frac{\mathit{\lambda }}{\mathit{\lambda }\mathrm{-}\mathrm{1}}\frac{\mathrm{1}}{{\mathit{\beta }}_{{\mathrm{m}}_{\mathrm{c}}}}\right)}^{-1}\mathit{,}$(25)where w_c is the enthalpy density at the centre and ${\mathit{\beta }}_{\mathrm{m}}\mathrm{=}\mathit{p}\mathit{/}{\mathit{p}}_{\mathrm{m}}\mathit{,}$(26)is the magnetisation parameter (β_mc being the magnetisation parameter at the centre of the disc). Using this definition, we can find the magnetic pressure at the centre, ${\mathit{p}}_{{\mathrm{m}}_{\mathrm{c}}}\mathrm{=}{\mathit{p}}_{\mathrm{c}}\mathit{/}{\mathit{\beta }}_{{\mathrm{m}}_{\mathrm{c}}}\mathit{.}$(27)With both pressures known at the centre, fluid and magnetic, we can now find the constants K and K_m using Eqs. (12)and (14). Therefore, for a given inner radius of the disc r_in we can obtain the potential W_in. We now finally have all the ingredients required to find the enthalpy density distribution (Eq. (20)), the fluid pressure and magnetic pressure distributions (Eqs. (12)and (14)), and the rest-mass density distribution, which can be trivially obtained inverting the barotropic relation, p = Kρ^κ. We note that the fluid enthalpy density w includes the rest-mass density ρ, and can thus be defined as w = ρh where h is the specific enthalpy. The relativistic definition of this quantity is h = 1 + ε + p/ρ, where ε is the specific internal energy. From a thermodynamical point of view a non-relativistic fluid satisfies h = 1. Therefore, since we use polytropic equations of state (relating p with either w or ρ in the same functional form) we are implicitly assuming that h = 1, that is, the disks we build in our procedure are non-relativistic from a thermodynamical point of view.

For the integration of Eq. (23)we have used the composite Simpson’s rule. It is important to use a very small integration step because the slope is very steep. In this work, we have used a step Δr = 10^-6. Using the analytic, constant angular momentum case for comparison, we tested that a larger value of Δr gives unacceptable accuracy losses. On the other hand, to integrate the ordinary differential Eq. (24)we have employed a fourth-order Runge-Kutta method. Again, it is also important here to choose a suitable step of integration, especially at the outer end of the disc, because Eq. (24)diverges at the equatorial plane (the equipotential surfaces cross the equatorial plane perpendicularly). We show the proof of this statement in Appendix A.

Fig. 2 Isodensity distributions for all models of Table 3. From left to right: the columns correspond to increasing values of the black hole spin: a = 0.5, 0.9 and 0.99. From top to bottom: the rows correspond to the following parameter combinations: i) γ = β = 0; ii) γ = β = 0.5; iii) γ = β = 0.9; iv) γ = 0.9, β = 0; v) γ = 0, β = 0.9. Note that the spatial scale of the plots differs as it has been chosen to facilitate the visualization of the disks.

Fig. 3 Effects of the magnetisation in the structure of the disc. From left to right: the values are βmc = 103, 1, and 10-3.

4. Results

To reduce the number of free parameters to build the initial models we have done as in Komissarov (2006) and have fixed the values of the equation of state exponents, κ and λ, and of the enthalpy density at the disc centre, w_c. More precisely, we have chosen κ = λ = 4/3 and w_c = 1. This still leaves us with five parameters to control the size, shape, thickness, and magnetisation of the disc, namely the magnetisation parameter β_mc, the parameters of the angular momentum distribution β and γ, the black hole spin parameter a, and the inner radius of the disc r_in. We have built a series of 45 models, whose main features are summarised in Table 1 (for β_mc = 10³, that is, models where the effects of the magnetisation are unimportant), Table 2 (β_mc = 1, that is, equipartition models), and Table 3 (β_mc = 10^-3, that is, highly-magnetised models). We note that we can build models with very small values of β_mc, such as for example β_mc = 10^-20 (extremely high magnetisation) without encountering numerical difficulties. No qualitative differences in the structure of the disks are found once β_mc becomes smaller than a sufficiently small value, about β_mc = 10^-3. This is the reason why we have chosen that particular value as our lower limit in the figures and tables of the manuscript.

We note in particular that the radii of the maximum fluid pressure and magnetic pressure, r_max and r_m,max, respectively, for the same value of β_mc, never coincide. This also reflects the fact that constant fluid pressure surfaces do not coincide with constant magnetic pressure surfaces.

We start by assessing our procedure by first building an extra disc model that can be directly compared with one of the two models in Komissarov (2006). This is shown in Fig. 1, which corresponds to the same constant specific angular momentum model A presented by Komissarov (2006). The parameters of this model are a = 0.9, β_mc = 0.1, and l = 2.8. The visual comparison shows that our approach can reproduce those previous results with good agreement. The rest-mass density distribution in the (rsinθ,rcosθ) plane shown in the left panel of Fig. 1 is remarkably similar to that shown in the left panel of Fig. 2 in Komissarov (2006). Not only the morphology of both models is nearly identical but also the range of variation of the rest-mass density and the location of the disc centre agree well in both cases. This can be most clearly seen in the middle and right panels of Fig. 1 which show, respectively, the angular profile at r_c and the radial profile at θ = π/ 2. The middle figure can be directly compared with Fig. 3 of Komissarov (2006). It is relevant to mention that very small changes in the location of the inner radius r_in have a significant effect on the maximum value of the rest-mass density, which explains the small differences between our figures and the ones presented in Komissarov (2006).

A representative sample of the isodensity surfaces for some of our models appears in Fig. 2. We plot, in particular, the models of Table 3, for which the magnetisation is highest (β_mc = 10^-3). From left to right the columns of this figure correspond to increasing values of the black hole spin, namely, a = 0.5, 0.9 and 0.99, while from top to bottom the rows correspond to different combinations of the γ and β parameters that characterise the ansatz for the angular momentum distribution of Qian et al. (2009, the particular values are indicated in the caption and legends of the figure). A rapid inspection shows that the structure of the disks noticeable changes when the parameters change. Notice that the spatial scale in all of the plots of this figure has been chosen so as to facilitate the visualization of the disks (which can be fairly small in some cases) and, as such, is different in each plot. A typical Polish doughnut is represented by the model in the top-left panel. As the black hole spin increases, the disks are closer to the black hole and its relative size with respect to the black hole is smaller for all values of γ and β. For γ = 0 the disks are infinite, irrespective of the value of a and β. As γ and β increase from 0 to 0.9 (compare the three top rows) the disks become significantly thinner and radially smaller (see also Fig. 5). It is also interesting to note that the maximum value of the rest-mass density is higher (with respect to the value at the disc centre) as the spin increases.

Fig. 4 Radial profiles of the rest-mass density in the equatorial plane for βmc = 103 (black curve), 1 (blue curve), and 10-3 (green curve).

Fig. 5 Radial profiles at the equatorial plane for models with a = 0.9, βmc = 10-3, and same combination of the γ and β parameters as in the rows of Fig. 2 (the specific values are indicated in the legend).

Figure 3 shows the effects of changing the parameter β_mc (the magnetisation) in the structure of the disks. From left to right the values plotted in each panel are β_mc = 10³, 1, and 10^-3. The model chosen corresponds to a = 0.9 and γ = β = 0.5. The larger the value of β_mc the less important the effects of the magnetisation in the disc structure. Figure 3 shows that, at least for the kind of magnetic field distribution we are considering, the effects of the magnetisation are minor. The disc structure remains fairly similar for all values of β_mc and the only quantitative differences are found in the location of the centre of the disc (which moves inward with decreasing β_mc) and of the range of variation of the isodensity contours (the maximum being slightly larger with decreasing β_mc). This can be more clearly visualised in Fig. 4 which displays the radial profile at the equatorial plane of the rest-mass density for the same three cases plotted in Fig. 3.

Next, we show in Fig. 5 the corresponding radial profiles at the equatorial plane for the models with black hole spin a = 0.9. We consider the case β_mc = 10^-3 (highest magnetisation) and the same combination of the γ and β parameters that we employed in Fig. 2. Therefore, this plot depicts the radial profiles of the models occupying the central column of Fig. 2. These profiles allow for a clearer quantification of the radial extent of the disks with γ and β. As mentioned in the description of Fig. 2 as γ and β increase from 0 to 0.9 (compare black, blue, and green curves) the disks become gradually smaller. At the same time, the radial location of the rest-mass density maximum increases and the maximum value of the central rest-mass density (and pressure) decreases. It is worth noticing that the radial profiles of some of the models overlap below a certain radius. More precisely, the models with equal β overlap below r ~ r_c (which corresponds to ρ ~ 1). In addition, this figure clearly shows that, besides modifying the thickness of the disc, β is the parameter which determines the value of the rest-mass density (and pressure) maximum. On the other hand, the parameter γ is only relevant for controlling the radial size of the disks.

In order to provide a quantitative comparison of the structural differences that appear along the sequence of equilibrium models, we plot in Fig. 6 the variation of the maximum of the rest-mass density (left panel), the radial location of the maximum of the fluid pressure (middle panel), and the radial location of the maximum of the magnetic pressure (right panel) as a function of the magnetisation parameter β_mc. All quantities plotted are measured at the equatorial plane. Let us first consider the left panel. All curves show the same monotonically decreasing trend with β_mc, yet for small enough and large enough values of β_mc the value of the rest-mass density maximum does not change. However, in the interval 10^-2 ≤ β_mc ≤ 10², ρ_max changes abruptly. The larger the black hole spin the larger the drop in the maximum of the rest-mass density. For sufficiently large values of β_mc (small magnetisation) the value of ρ_max stays constant to the same value irrespective of the black hole rotation.

In the middle panel of Fig. 6 we show the radial location of the maximum of the fluid pressure as a function of β_mc. The black hole spin is a = 0.99 and the values of β and γ are indicated in the legend. For all values of β and γ, the maximum of the location of the disc fluid pressure decreases with decreasing β_mc. Below β_mc ~ 10^-3 the location of the maximum does no longer change. The panel also shows that above β_mc ~ 10³ the constant value of r_max is the same irrespective of β and γ, as expected, because for purely hydrodynamic disks, r_max = r_c.

The dependence of the location of the maximum of the magnetic pressure with β_mc is shown in the right panel of Fig. 6. While this location also decreases with decreasing β_mc, the constant value achieved for values above β_mc ~ 10³ does depend on the specific values of β and γ, contrary to what happens with the fluid pressure. It must be noted that the location of the maximum of the magnetic pressure is identical for all models considered when β_mc ≡ 1/(λ−1) = 3, as we show in Appendix B. At this value of β_mc all models cross at r_mmax = r_c, as it is more clearly shown in the inset of the right panel of Fig. 6.

As a final remark we note that in the left and middle panels of Fig. 6 we do not show the data for models with β ≠ γ as they coincide with the respective models with the same value of β. Moreover, the middle and right panels only show the data for a = 0.99 because changing the value of the spin parameter of the black hole only yields a change of scale. Therefore, if we chose the range of the graph accordingly, the plots would be identical (same relative differences between the curves) irrespective of the value of a.

Fig. 6 Effects of βmc. Left panel: variation of the value of the maximum rest-mass density with respect to the logarithm of the magnetisation parameter. The solid, dashed, and dotted lines refer to a = 0.99, 0.9 and 0.5, respectively, and the colour code refers to the models shown in the legend. Middle panel: variation of the location of the rest-mass density maximum (and fluid pressure) with respect to the logarithm of the magnetisation parameter for the three models shown in the legend with a = 0.99. Right panel: variation (also for a = 0.99) of the location of the maximum of the magnetic pressure with respect to the logarithm of βmc.

5. Conclusions

In this paper we have presented a procedure to build equilibrium sequences of magnetised, non-self-gravitating disks around Kerr black holes which combines the two existing approaches of Komissarov (2006) and Qian et al. (2009). On the one hand we have followed Qian et al. (2009) and have assumed a form of the angular momentum distribution in the disc from which the location and morphology of surfaces of equal potential can be computed. As a limiting case, this ansatz includes the constant angular momentum case originally employed in the construction of thick tori – or Polish doughnuts (Abramowicz et al. 1978; Kozlowski et al. 1978) – and was already used by Qian et al. (2009) to build equilibrium sequences of purely hydrodynamical models. On the other hand, our disks are endowed with a purely toroidal magnetic field, as in the work of Komissarov (2006), which provides the methodology we have followed to handle the magnetic terms. On a similar note, we cite the work of Wielgus et al. (2015) who have recently extended the solution of Komissarov (2006) to include non-constant specific angular momentum tori. These authors limited their consideration to power-law distributions and were particularly focused on the stability of such tori to the MRI. The approach discussed in our work differs from that of Wielgus et al. (2015) and, moreover, we have explored a much wider parameter space.

We have discussed the properties of the new models and their dependence on the initial parameters, namely the magnetisation parameter β_mc, the parameters β and γ describing the angular momentum distribution, the black hole spin parameter a, and the inner radius of the disc r_in. We have shown the effects of changing β and γ beyond the purely hydrodynamical case considered in Qian et al. (2009). The morphology of the solutions we have built no longer changes for magnetisation values above β_mc ~ 10³ and below β_mc ~ 10^-3. These cases can thus be considered as the hydrodynamical and MHD limiting cases, respectively. The new sequences of magnetised disks around black holes presented in this work can be used as initial data for magnetohydrodynamical evolutions in general relativity. In the near future, we plan to extend this work along two main directions, namely (i) including the self-gravity of the disks, following the approach laid out in Stergioulas (2011); and (ii) constructing accretion disks around hairy black holes, both, with scalar and Proca hair (Herdeiro & Radu 2014; Herdeiro et al. 2016).

Acknowledgments

It is a pleasure to thank P. C. Fragile and N. Stergioulas for useful discussions and comments. We thank the referee for his constructive report which improved the manuscript. Work supported by the Spanish MINECO (grant AYA2015-66899-C2-1-P) and the Generalitat Valenciana (PROMETEOII-2014-069).

References

Appendix A: Divergence of Eq. (24) at θ = π/ 2

To prove that Eq. (24)diverges at θ = π/ 2 we need to show that Eq. (22)vanishes at the equator, that is $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{\theta }}\mathit{W}\mathrm{=}{\mathit{\partial }}_{\mathit{\theta }}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{-}\frac{\mathrm{\Omega }{\mathit{\partial }}_{\mathit{\theta }}\mathit{l}}{\mathrm{1}\mathrm{-}\mathrm{\Omega }\mathit{l}}\mathrm{=}\mathrm{0.}& & \end{array}$(A.1)Then, we have to prove that either the two terms of this equation are equal or they are identically zero at θ = π/ 2. We shall see that the latter is true. First, we take the derivative of the angular momentum distribution ∂_θl. From Eq. (7) we can write the angular momentum distribution as l(r,θ) = l(r)sin^2γθ, so that the derivative reads $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{\theta }}\mathit{l}\mathrm{=}\mathrm{2}\mathit{\gamma l}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\sin }^{\mathrm{2}\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{\theta }\cos \mathit{\theta }\mathit{.}& & \end{array}$(A.2)Taking θ = π/ 2 leads to ∂_θl = 0, so the second term equals to zero. To show that the first term is also zero, we write it as $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{\theta }}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{=}{\mathit{\partial }}_{\mathit{\theta }}\ln {\left(\frac{\mathrm{ℒ}}{\mathrm{𝒜}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}& & \end{array}$(A.3)where $\mathrm{ℒ}\mathrm{=}{\mathit{g}}_{\mathit{t\phi }}^{\mathrm{2}}\mathrm{-}{\mathit{g}}_{\mathit{tt}}{\mathit{g}}_{\mathit{\phi \phi }}$ and , and we have dropped the absolute value, as it is irrelevant to this discussion. Then, the derivative is $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{\theta }}\ln \mathrm{\left|}{\mathit{u}}_{\mathit{t}}\mathrm{\right|}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{𝒜}}{\mathrm{ℒ}}\frac{\mathrm{𝒜}{\mathit{\partial }}_{\mathit{\theta }}\mathrm{ℒ}\mathrm{-}\mathrm{ℒ}{\mathit{\partial }}_{\mathit{\theta }}\mathrm{𝒜}}{{\mathrm{𝒜}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{𝒜}{\mathit{\partial }}_{\mathit{\theta }}\mathrm{ℒ}\mathrm{-}\mathrm{ℒ}{\mathit{\partial }}_{\mathit{\theta }}\mathrm{𝒜}}{\mathrm{𝒜ℒ}}\mathrm{·}& & \end{array}$(A.4)We use Boyer-Lindquist coordinates, for which the metric components read $\begin{array}{ccc}& & {\mathit{g}}_{\mathit{tt}}\mathrm{=}\mathrm{-}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathit{Mr}}{{\mathit{\rho }}^{\mathrm{2}}}\right)\mathit{,}{\mathit{g}}_{\mathit{t\phi }}\mathrm{=}\mathrm{-}\frac{\mathrm{2}\mathit{Mar}{\sin }^{\mathrm{2}}\mathit{\theta }}{{\mathit{\rho }}^{\mathrm{2}}}\mathit{,}\\ & & {\mathit{g}}_{\mathit{\phi \phi }}\mathrm{=}\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{2}\mathit{M}{\mathit{a}}^{\mathrm{2}}\mathit{r}{\sin }^{\mathrm{2}}\mathit{\theta }}{{\mathit{\rho }}^{\mathrm{2}}}\right){\sin }^{\mathrm{2}}\mathit{\theta ,}\end{array}$(A.5)where ρ² = r² + a²cos²θ. Therefore, we obtain ℒ = Δsin²θ, where Δ = r²−2Mr + a², and its derivative ∂_θℒ = 2Δsinθcosθ, which is zero for θ = π/ 2. Next, we take the derivative of $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{\theta }}\mathrm{𝒜}& \mathrm{=}& {\mathit{\partial }}_{\mathit{\theta }}{\mathit{g}}_{\mathit{\phi \phi }}\mathrm{+}\mathrm{2}\mathrm{\left(}{\mathit{\partial }}_{\mathit{\theta }}\mathit{l}\mathrm{\right)}{\mathit{g}}_{\mathit{t\phi }}\mathrm{+}\mathrm{2}\mathit{l}{\mathit{\partial }}_{\mathit{\theta }}{\mathit{g}}_{\mathit{t\phi }}\mathrm{+}\mathrm{2}\mathit{l}\mathrm{\left(}{\mathit{\partial }}_{\mathit{\theta }}\mathit{l}\mathrm{\right)}{\mathit{g}}_{\mathit{tt}}\mathrm{+}{\mathit{l}}^{\mathrm{2}}{\mathit{\partial }}_{\mathit{\theta }}{\mathit{g}}_{\mathit{tt}}\\ & \mathrm{=}& {\mathit{\partial }}_{\mathit{\theta }}{\mathit{g}}_{\mathit{\phi \phi }}\mathrm{+}\mathrm{2}\mathit{l}{\mathit{\partial }}_{\mathit{\theta }}{\mathit{g}}_{\mathit{t\phi }}\mathrm{+}{\mathit{l}}^{\mathrm{2}}{\mathit{\partial }}_{\mathit{\theta }}{\mathit{g}}_{\mathit{tt}}\mathit{,}\end{array}$(A.6)where we have used the previous result that ∂_θl = 0 at θ = π/ 2. By inspecting the metric components, it is easy to see that all terms depending on θ are functions of sin²θ or cos²θ. This means their θ derivatives will have at least a cosθ multiplying. Then at θ = π/ 2. Therefore, ∂_θln | u_t | is also zero at θ = π/ 2.

Appendix B: Value of the magnetisation parameter for r_mmax = r_c

In this appendix we derive the condition β_mc = 1/(λ−1) for r_mmax = r_c. First, we can use Eq. (14)to write $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{p}}_{\mathrm{m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{\partial r}}\mathrm{=}{\mathit{\partial }}_{\mathit{r}}\mathrm{\left(}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{w}}^{\mathit{\lambda }}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(B.1)which yields $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{r}}{\mathit{p}}_{\mathrm{m}}\mathrm{=}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{2}}{\mathit{w}}^{\mathit{\lambda }}\mathrm{\left[}\mathrm{\right(}\mathit{\lambda }\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right(}{\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}\mathrm{)}\mathit{w}\mathrm{+}\mathit{\lambda }\mathrm{ℒ}\mathrm{(}{\mathit{\partial }}_{\mathit{r}}\mathit{w}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\mathrm{0.}& & \end{array}$(B.2)\newpageInside the disc K_mℒ^λ−2w^λ ≠ 0. Therefore, to fulfill the extremum condition (B.1)we need $\mathrm{(}\mathit{\lambda }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{\left(}{\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}\mathrm{\right)}\mathit{w}\mathrm{+}\mathit{\lambda }\mathrm{ℒ}\mathrm{\left(}{\mathit{\partial }}_{\mathit{r}}\mathit{w}\mathrm{\right)}\mathrm{=}\mathrm{0.}$(B.3)To evaluate this expression, we need to compute the partial derivatives ∂_rℒ and ∂_rw. The derivative of ℒ is straightforward, since ℒ = Δsin²θ, ${\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{r}\mathrm{-}\mathit{M}\mathrm{)}{\sin }^{\mathrm{2}}\mathit{\theta }\mathit{.}$(B.4)Let us now discuss ∂_rw. From Eq. (20)we have $\mathit{w}\mathrm{=}{\left[\mathrm{-}\mathrm{\Delta }\mathit{W}\left(\frac{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{\lambda }}\right)\frac{\mathrm{1}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right]}^{\frac{\mathrm{1}}{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\mathrm{·}$(B.5)Taking its derivative we obtain $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{r}}\mathit{w}& \mathrm{=}& \frac{\mathit{w}}{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\left[\mathrm{-}\mathrm{\Delta }\mathit{W}\left(\frac{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{\lambda }}\right)\frac{\mathrm{1}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right]}^{-1}\\ & & \mathrm{\times }[\left(\mathrm{-}{\mathit{\partial }}_{\mathit{r}}\mathit{W}\left(\frac{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{\lambda }}\right)\frac{\mathrm{1}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right)\\ & & \mathrm{+}\left(\mathrm{-}\mathrm{\Delta }\mathit{W}\left(\frac{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{\lambda }}\right)\frac{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{2}}{\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}}{\mathrm{(}\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}}\right)]\mathrm{·}\end{array}$(B.6)Since at r = r_c the derivative of the potential is zero, ∂_rW(r,π/ 2) | _{r = rc} = 0, we can simplify the above expression to obtain $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{r}}\mathit{w}& \mathrm{=}& \frac{\mathit{w}}{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\left[\mathrm{-}\mathrm{\Delta }\mathit{W}\left(\frac{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{\lambda }}\right)\frac{\mathrm{1}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right]}^{-1}\\ & & \mathrm{\times }\left(\mathrm{-}\mathrm{\Delta }\mathit{W}\left(\frac{\mathit{\lambda }\mathrm{-}\mathrm{1}}{\mathit{\lambda }}\right)\frac{\mathrm{1}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right)\\ & & \mathrm{\times }\left[\frac{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{2}}{\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right]\mathit{,}\end{array}$(B.7)which can be further simplified to the following form ${\mathit{\partial }}_{\mathit{r}}\mathit{w}\mathrm{=}\mathrm{-}\mathit{w}\left[\frac{{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{2}}{\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}}{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}}\right]\mathrm{·}$(B.8)Next, by inserting this expression in Eq. (B.3)we obtain $\mathrm{\left(}{\mathit{\partial }}_{\mathit{r}}\mathrm{ℒ}\mathrm{\right)}\mathit{w}\mathrm{\left[}\mathrm{\right(}\mathit{\lambda }\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right(}\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}\mathrm{)}\mathrm{-}\mathit{\lambda }{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}\mathrm{]}\mathrm{=}\mathrm{0.}$(B.9)Since (∂_rℒ)w ≠ 0 inside the disc, the above equation leads to $\mathit{K}\mathrm{(}\mathit{\lambda }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{=}{\mathit{K}}_{\mathrm{m}}{\mathrm{ℒ}}^{\mathit{\lambda }\mathrm{-}\mathrm{1}}\mathit{.}$(B.10)Moreover, we can use Eqs. (12), (14)and (26)at the centre (remember that w_c = 1) to write K as K = β_mcK_mℒ^λ−1. Therefore, this allows us to obtain the following simple expression for the value of the magnetisation parameter at the centre of the disc ${\mathit{\beta }}_{{\mathrm{m}}_{\mathrm{c}}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\lambda }\mathrm{-}\mathrm{1}}\mathrm{·}$(B.11)It is relevant to note that the only reference to the explicit form of the metric is done to show that ℒ ≠ 0 and ∂_rℒ ≠ 0 inside the disc. This means that our result holds for any stationary and axisymmetric metric if the aforementioned inequalities hold. Moreover, the result is also true for any angular momentum distribution that allows for the existence of a cusp and a centre (and r_c ≠ r_cusp).

All Tables

All Figures

Fig. 1 Comparison with Komissarov’s solution for a constant angular momentum model. The left panel shows the rest-mass density distribution in logarithm scale for all radii and all angles while the middle and right panels show, respectively, the angular profile of the logarithm of the rest-mass density at the centre of the disc and the radial profile of the logarithm of the rest-mass density at the equatorial plane.

In the text

Fig. 2 Isodensity distributions for all models of Table 3. From left to right: the columns correspond to increasing values of the black hole spin: a = 0.5, 0.9 and 0.99. From top to bottom: the rows correspond to the following parameter combinations: i) γ = β = 0; ii) γ = β = 0.5; iii) γ = β = 0.9; iv) γ = 0.9, β = 0; v) γ = 0, β = 0.9. Note that the spatial scale of the plots differs as it has been chosen to facilitate the visualization of the disks.

In the text

	Fig. 3 Effects of the magnetisation in the structure of the disc. From left to right: the values are βmc = 103, 1, and 10-3.
In the text

	Fig. 4 Radial profiles of the rest-mass density in the equatorial plane for βmc = 103 (black curve), 1 (blue curve), and 10-3 (green curve).
In the text

	Fig. 5 Radial profiles at the equatorial plane for models with a = 0.9, βmc = 10-3, and same combination of the γ and β parameters as in the rows of Fig. 2 (the specific values are indicated in the legend).
In the text

Fig. 6 Effects of βmc. Left panel: variation of the value of the maximum rest-mass density with respect to the logarithm of the magnetisation parameter. The solid, dashed, and dotted lines refer to a = 0.99, 0.9 and 0.5, respectively, and the colour code refers to the models shown in the legend. Middle panel: variation of the location of the rest-mass density maximum (and fluid pressure) with respect to the logarithm of the magnetisation parameter for the three models shown in the legend with a = 0.99. Right panel: variation (also for a = 0.99) of the location of the maximum of the magnetic pressure with respect to the logarithm of βmc.

In the text