Issue 
A&A
Volume 607, November 2017



Article Number  A68  
Number of page(s)  9  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201730935  
Published online  15 November 2017 
Magnetised Polish doughnuts revisited
^{1} Departamento de Astronomía y Astrofísica, Universitat de València, Dr. Moliner 50, 46100 Burjassot (València), Spain
email: sergio.gimeno@uv.es
^{2} Observatori Astronòmic, Universitat de València, C/ Catedrático José Beltrán 2, 46980 Paterna (València), Spain
email: j.antonio.font@uv.es
Received: 4 April 2017
Accepted: 11 July 2017
We discuss a procedure to build new sequences of magnetised, equilibrium tori around Kerr black holes which combines two approaches previously considered in the literature. For simplicity we have assumed that the testfluid approximation holds, and hence we neglected the selfgravity of the fluid. The models were built assuming a particular form of the angular momentum distribution from which the location and morphology of equipotential surfaces can be computed. This ansatz includes, in particular, the constant angular momentum case originally employed in the construction of thick tori – or Polish doughnuts – and it has already been used to build equilibrium sequences of purely hydrodynamical models. We discuss the properties of the new models and their dependence on the initial parameters. These new sequences can be used as initial data for magnetohydrodynamical evolutions in general relativity.
Key words: accretion, accretion disks / magnetohydrodynamics (MHD) / black hole physics
© 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 Xray binaries and microquasars, as well as in the central engine of gammaray 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 nonmagnetised 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 powerlaw distributions of angular momentum. Moreover, the magnetised tori of Wielgus et al. (2015) were used to explore the growth of the magnetorotational instability (MRI) through timedependent numerical simulations. In particular, the longterm 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 testfluid approximation holds, thus neglecting the selfgravity 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 BoyerLindquist 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 where we have defined the following quantities, , , 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, (3)where u^{μ} is the fluid fourvelocity. The relationship between l and Ω is given by the equations (4)where g_{μν} is the metric tensor and we are assuming circular motion, that is, the fourvelocity can be written as (5)We also introduce the Keplerian angular momentum (for prograde motion) in the equatorial plane l_{K}, defined as (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 (7)where constants l_{0} and l_{ms}(r) are defined by l_{0} ≡ ηl_{K}(r_{ms}) and l_{ms}(r) ≡ l_{0} [l_{K}(r_{ms}) /l_{0}] ^{β}. 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) (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_{0} as l_{0} = 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 (9)is the energymomentum tensor of a magnetised perfect fluid, w and p being the fluid enthalpy density and fluid pressure, respectively. Moreover, is the (dual of the) Faraday tensor relative to an observer with fourvelocity u^{μ}, and b^{μ} is the magnetic field in that frame, with b^{2} = 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 energymomentum tensor with the projection tensor , we arrive at (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 (11)where . To integrate Eq. (11)we first assume a barotropic equation of state w = w(p) of the form (12)with K and κ constants. Then, we define the magnetic pressure as p_{m} = b^{2}/ 2, and introduce the definitions and , in order to write an analogue equation to Eq. (12)for (Komissarov 2006) (13)or, in terms of the magnetic pressure p_{m}(14)where K_{m} and λ are constants. This particular choice of barotropic relationships, w = w(p) and , 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 (15)On the surface of the disc, and particularly on its inner edge, the conditions are satisfied and, therefore, the integration constant is simply given by (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 (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 (18)which for a barotropic equation of state can be easily integrated to give (19)Replacing p and p_{m} by Eqs. (12)and (14), the previous equation reduces to (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)(21)and (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}(23)Following Qian et al. (2009; see also Jaroszynski et al. 1980) we can divide Eqs. (21)and (22)to obtain (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 restmass 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 restmass density at the centre of the disc and the radial profile of the logarithm of the restmass 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), (25)where w_{c} is the enthalpy density at the centre and (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, (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 restmass density distribution, which can be trivially obtained inverting the barotropic relation, p = Kρ^{κ}. We note that the fluid enthalpy density w includes the restmass 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 nonrelativistic 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 nonrelativistic 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 fourthorder RungeKutta 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.
Models for β_{mc} = 10^{3}.
Models for β_{mc} = 1.
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} = 10^{3}, 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^{3}, 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, highlymagnetised 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 restmass 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 restmass 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 restmass density, which explains the small differences between our figures and the ones presented in Komissarov (2006).
Models for β_{mc} = 10^{3}.
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 topleft 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 restmass density is higher (with respect to the value at the disc centre) as the spin increases.
Fig. 4 Radial profiles of the restmass density in the equatorial plane for β_{mc} = 10^{3} (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^{3}, 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 restmass 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 restmass density maximum increases and the maximum value of the central restmass 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 restmass 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 restmass 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 restmass density maximum does not change. However, in the interval 10^{2} ≤ β_{mc} ≤ 10^{2}, ρ_{max} changes abruptly. The larger the black hole spin the larger the drop in the maximum of the restmass 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^{3} 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^{3} 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 restmass 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 restmass 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, nonselfgravitating 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 nonconstant specific angular momentum tori. These authors limited their consideration to powerlaw 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^{3} 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 selfgravity 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 AYA201566899C21P) and the Generalitat Valenciana (PROMETEOII2014069).
References
 Abramowicz, M. A., & Fragile, P. C. 2013, Liv. Rev. Relat., 16, 1 [Google Scholar]
 Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221 [NASA ADS] [Google Scholar]
 Baiotti, L., & Rezzolla, L., 2017, Rep. Prog. Phys. 80, 096901 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Faber, J. A., & Rasio, F. A. 2012, Liv. Rev. Relat., 15, 8 [CrossRef] [Google Scholar]
 Font, J. A., & Daigne, F. 2002, MNRAS, 334, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Fragile, P. C., & Sądowski, A. 2017, MNRAS, 467, 1838 [NASA ADS] [CrossRef] [Google Scholar]
 Herdeiro, C. A. R., & Radu, E. 2014, Phys. Rev. Lett., 112, 221101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Herdeiro, C., Radu, E., & Rúnarsson, H. 2016, Class. Quant. Grav., 33, 154001 [NASA ADS] [CrossRef] [Google Scholar]
 Jaroszynski, M., Abramowicz, M. A., & Paczynski, B. 1980, Acta Astron., 30, 1 [NASA ADS] [Google Scholar]
 Komissarov, S. S. 2006, MNRAS, 368, 993 [NASA ADS] [CrossRef] [Google Scholar]
 Kozlowski, M., Jaroszynski, M., & Abramowicz, M. A. 1978, A&A, 63, 209 [NASA ADS] [Google Scholar]
 Montero, P. J., Zanotti, O., Font, J. A., & Rezzolla, L. 2007, MNRAS, 378, 1101 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, L., Abramowicz, M. A., Fragile, P. C., et al. 2009, A&A, 498, 471 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rezzolla, L., Baiotti, L., Giacomazzo, B., Link, D., & Font, J. A. 2010, Class. Quant. Grav., 27, 114105 [NASA ADS] [CrossRef] [Google Scholar]
 Sekiguchi, Y., & Shibata, M. 2011, ApJ, 737, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Shibata, M., & Taniguchi, K. 2011, Liv. Rev. Relat., 14, 6 [Google Scholar]
 Stergioulas, N. 2011, J. Phys. Conf. Ser., 283, 012036 [CrossRef] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Wielgus, M., Fragile, P. C., Wang, Z., & Wilson, J. 2015, MNRAS, 447, 3593 [NASA ADS] [CrossRef] [Google Scholar]
 Zanotti, O., & Pugliese, D. 2015, Gen. Relat. Gravit., 47, 44 [NASA ADS] [CrossRef] [Google Scholar]
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 (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 (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 (A.3)where and , and we have dropped the absolute value, as it is irrelevant to this discussion. Then, the derivative is (A.4)We use BoyerLindquist coordinates, for which the metric components read (A.5)where ρ^{2} = r^{2} + a^{2}cos^{2}θ. Therefore, we obtain ℒ = Δsin^{2}θ, where Δ = r^{2}−2Mr + a^{2}, and its derivative ∂_{θ}ℒ = 2Δsinθcosθ, which is zero for θ = π/ 2. Next, we take the derivative of (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^{2}θ or cos^{2}θ. 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 (B.1)which yields (B.2)\newpageInside the disc K_{m}ℒ^{λ−2}w^{λ} ≠ 0. Therefore, to fulfill the extremum condition (B.1)we need (B.3)To evaluate this expression, we need to compute the partial derivatives ∂_{r}ℒ and ∂_{r}w. The derivative of ℒ is straightforward, since ℒ = Δsin^{2}θ, (B.4)Let us now discuss ∂_{r}w. From Eq. (20)we have (B.5)Taking its derivative we obtain (B.6)Since at r = r_{c} the derivative of the potential is zero, ∂_{r}W(r,π/ 2)  _{r = rc} = 0, we can simplify the above expression to obtain (B.7)which can be further simplified to the following form (B.8)Next, by inserting this expression in Eq. (B.3)we obtain (B.9)Since (∂_{r}ℒ)w ≠ 0 inside the disc, the above equation leads to (B.10)Moreover, we can use Eqs. (12), (14)and (26)at the centre (remember that w_{c} = 1) to write K as K = β_{mc}K_{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 (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 restmass 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 restmass density at the centre of the disc and the radial profile of the logarithm of the restmass 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} = 10^{3}, 1, and 10^{3}. 

In the text 
Fig. 4 Radial profiles of the restmass density in the equatorial plane for β_{mc} = 10^{3} (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 restmass 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 restmass 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.