Issue 
A&A
Volume 614, June 2018



Article Number  A75  
Number of page(s)  11  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201732361  
Published online  15 June 2018 
New closed analytical solutions for geometrically thick fluid tori around black holes
Numerical evolution and the onset of the magnetorotational instability
Center of Applied Space Technology and Microgravity (ZARM), Universität Bremen,
Am Falturm 2,
28359
Bremen, Germany
email: vojtech.witzany@zarm.unibremen.de; paul.jefremow@zarm.unibremen.de
Received:
25
November
2017
Accepted:
24
February
2018
Context. When a black hole is accreting well below the Eddington rate, a geometrically thick, radiatively inefficient state of the accretion disk is established. There is a limited number of closedform physical solutions for geometrically thick (nonselfgravitating) toroidal equilibria of perfect fluids orbiting a spinning black hole, and these are predominantly used as initial conditions for simulations of accretion in the aforementioned mode. However, different initial configurations might lead to different results and thus observational predictions drawn from such simulations.
Aims. We aim to expand the known equilibria by a number of closed multiparametric solutions with various possibilities of rotation curves and geometric shapes. Then, we ask whether choosing these as initial conditions influences the onset of accretion and the asymptotic state of the disk.
Methods. We have investigated a set of examples from the derived solutions in detail; we analytically estimate the growth of the magnetorotational instability (MRI) from their rotation curves and evolve the analytically obtained tori using the 2D magnetohydrodynamical code HARM. Properties of the evolutions are then studied through the mass, energy, and angularmomentum accretion rates.
Results. The rotation curve has a decisive role in the numerical onset of accretion in accordance with our analytical MRI estimates: in the first few orbital periods, the average accretion rate is linearly proportional to the initial MRI rate in the toroids. The final state obtained from any initial condition within the studied class after an evolution of ten or more orbital periods is mostly qualitatively identical and the quantitative properties vary within a single order of magnitude. The average values of the energy of the accreted fluid have an irregular dependency on initial data, and in some cases fluid with energies many times its rest mass is systematically accreted.
Key words: accretion, accretion disks / black hole physics – magnetohydrodynamics (MHD) / methods: analytical / methods: numerical
© ESO 2018
1 Introduction
Accretion of matter onto black holes can lead to the release of large amounts of radiation and is responsible for some of the most energetic phenomena in the Universe such as Xray binaries or active galactic nuclei. When the accretion rate on a black hole is sufficiently high while radial convection of internal energy sufficiently low, the accretion disk surrounding the black hole reaches a state which is well described by the socalled thin disk model (Shakura & Sunyaev 1973; Novikov & Thorne 1973; Abramowicz & Fragile 2013). However, when the accretion rate drops well below the socalled Eddington accretion rate Ṁ_{Edd} = 1.39 × 10^{18}(M∕M_{⊙})g s^{−1}, as is believed to be, for example, the case for the black hole in our Galactic center, the accretion disk enters a rather different mode that is characterized as geometrically thick, optically thin, “hot” in the sense that thermal energies are nonnegligible as compared to gravitational binding energies, and radiatively inefficient. This also implies that the energy losses from the disk are dominated by either direct advection through the horizon or outflows (see Yuan & Narayan 2014, and references therein).
Amongst other things, this means that the radiatively inefficient disk can be, to a first approximation, modeled as a magnetohydrodynamic (MHD) fluid without radiation backreaction (see e.g., Dexter et al. 2009). A more accurate model, however, requires the inclusion of radiative cooling and heating, and the partially twotemperature nature of the fluid in this accretion mode, see for example Ryan et al. (2017) and references therein. Consequently, a number of global generalrelativistic MHD studies were conducted showing that the MRIdriven turbulence (Balbus & Hawley 1991, 1998) provides the angular momentum transport in the disk, and many features predicted by analytical or semianalytical models emerge with only a small amount of additional adhoc input (see e.g., De Villiers et al. 2003; McKinney & Gammie 2004; Narayan et al. 2012).
For simulations of this type of accretion disks one usually uses initial conditions which are in a smooth equilibrium state. If such equilibria are to be given completely in closed form, then they are predominantly either the fluid tori with constant specific angular momentum ℓ ≡−u_{φ}∕u_{t} of Kozłowski et al. (1978); Abramowicz et al. (1978), or with constant ℓ* ≡ u_{φ}u^{t} of Fishbone & Moncrief (1976).
However, the angular momentum or angular velocity distributions play a decisive role in the onset of various instabilities. The MRI growth rate is proportional to the angular velocity gradient, and another notable instability, the centrifugal or Rayleigh–Taylor^{1} instability, has a growth rate proportional to the gradient of angular momentum.
Another case where the angular momentum distribution was formerly believed to make a crucial difference is the gravitational runaway instability as studied for example by Font & Daigne (2002); Daigne & Font (2004). In the particular model chosen by the authors of the aforementioned study, the mass parameter of the black hole spacetime was allowed to dynamicallygrow by the amount of accreted matter, possibly leading to an accretion runaway. Then, they found that if specific angular momentum grew with radius, the runaway instability was suppressed, whereas for the ℓ = const. tori it caused the accretion disk to be fully accreted on orbital timescales. This instability was later shown to disappear when full selfgravity of the torus is accounted for (Montero et al. 2010; Rezzolla et al. 2010; Mewes et al. 2016). However, the study of Font & Daigne (2002); Daigne & Font (2004) demonstrates that the precise shape of the angular momentum distribution can often have a critical influence on the evolution obtained within a given model.
Is it possible that we are missing new dynamical effects or different states of the accretion disks in our simulations because we are restricting ourselves only to a limited number of initial conditions? To answer this question, a variety of closed form solutions corresponding to different gradients of either ℓ or Ω in the initial fluid configurations would be useful.
We present such solutions in this paper in Sect. 2, and demonstrate how these lead, at least in principle, to differences in numerical simulation results in Sect. 3 by implementing and evolving our solutions in the HARM 2D code (Gammie et al. 2003; Noble et al. 2006). A reader looking for an explicit recipe for the constructionof the tori will find all the needed formulas in Appendix B. Of interest is also the essentially unknown fact that the fully general solutions for the toroidal equilibria in static spacetimes (e.g., in the Schwarzschild spacetime) are expressible in closed form, which we briefly describe in Sect. 2.2.
2 Analytical solutions for fluid tori near black holes
We use the G = c = 1 geometrized units and the −+++ signature of the spacetime metric. A comma before an index defines a partial derivative with respect to the appropriate coordinate, whereas a preceding semicolon a covariant derivative with respect to the coordinate.
2.1 Basic equations
The plasma orbiting the black hole is modeled by ideal MHD where molecular dissipation, electric resistivity, selfgravitation, and radiative or nonequilibrium effects are neglected in the dynamics. As a result, the system is governed by the set of equations (cf. Anile 1989)
where u^{μ} is the fluid fourvelocity, μ its restmass density, w its enthalpy density, and P its pressure. The vector b^{μ} is the magnetic field vector in the rest frame of the fluid.
We place the magnetofluid into an axially symmetric and stationary spacetime with the metric (6)
Some useful symbols we use are the rotation frequency of zeroangularmomentum observers (ZAMOs), ω ≡ g^{tφ}∕g^{tt} = −g_{tφ}∕g_{φφ}; the potential whose gradient generates the acceleration of ZAMOs, Φ_{(Z)} ≡−ln(−g^{tt})∕2; minus the determinant of the tφ part of the metric, ; and a radiuslike quantity .
Now let us assume that all the properties of the magnetofluid are axisymmetric and stationary, the flow is purely circular, u^{r} = u^{ϑ} = 0, and that the magnetic field is purely toroidal, b^{r} = b^{ϑ} = 0. Then the Faraday Eq. (3) and the particleconservation Eq. (2) are trivially fulfilled and we need to solve only the momentum balance equation (relativistic Euler Eq. (1)).
The magnetic part of this equation can under the given assumptions be simplified as (7)
This expression was first derived by Komissarov (2006), but we provide, in our view, a more direct and clear rederivation in Appendix A.
There are many different ways how to rewrite the fluid part of the Euler equation, one of them is the form which can be attributed to Fishbone & Moncrief (1976) (8)
where is the linear velocity of the fluid as measured in the ZAMO frame.
However, one of the best known forms of the fluid part of the Euler equation for circular equilibria is due to Kozłowski et al. (1978); Abramowicz et al. (1978) and reads (9)
where ℓ ≡−u_{φ}∕u_{t} and Ω = u^{φ}∕u^{t}. In all of the Eqs. (7), (8) and (9) the nontrivial components are of course only those with ν = r, ϑ. We provide brief descriptions of the derivations of the expressions Eq. (8) and (9) in Appendix A.
2.2 Static spacetimes
In static spacetimes (g_{tφ} = 0) the full Euler equation using Eq. (7) and (8) simplifies as (10)
where we have introduced the notation .
Let us now further assume that the configuration of the fluid fulfills the barotropic^{2} condition of coincident surfaces of constant pressure and enthalpy P = P(w) and the “magnetotropic” condition of coincident constant and , (see Komissarov 2006). Then the righthand side of Eq. (10) is a coordinate gradient, and the same must be true also for the lefthand side. Thus we conclude that under the given conditions we necessarily have and the general solution of the Euler equation reads. (11)
where we have defined the effective potentials
and where the potentials are determined only up to integration constants. The functions W and can be understood as thermodynamic and magnetothermodynamic effective potentials respectively, Φ_{c} has clearly the meaning of a centrifugal potential, and Φ_{(Z)} is, in the currently discussed case of static spacetimes, the gravitational potential experienced by static observers.
General solutions for the toroidal equilibria can then be constructed by
 1.
Prescribing arbitrary distributions (e.g., as an “enthalpic polytrope” , see Komissarov 2006; Wielgus et al. 2015) and an equally arbitrary rotation profile ;
 2.
 3.
Solving Eq. (11) for w(r, ϑ);
 4.
Using and the relation b^{μ}u_{μ} = 0 to determine b^{t}, b^{φ};
 5.
Postulating an equation of state such as the idealgas w = μ + 5μk_{B}T∕(2m), where T is the temperature and m the particle mass, to derive the density and entropy field from fundamental thermodynamic relations and the function P(w).
The results presented in this section are probably not original, even though it is hard to find their explicit statement in the literature. They can be, however, seen as easily obtainable consequences of the results of either Abramowicz (1971) or Fishbone & Moncrief (1976) when restricted to static spacetimes.
Nonetheless, we find it important to state these results explicitly because there is a number of studies that construct tori in static spacetimes and use somewhat complicated numerical methods tailored for stationary spacetimes, even though the abovestated general closed solution is available (e.g., Daigne & Font 2004; Narayan et al. 2012; Penna et al. 2013).
2.3 Stationary spacetimes
We see above that the form of the Euler equation coming from Eq. (8) can be used to construct the entirely general barotropic and magnetotropic solutions in static spacetimes. Furthermore, Fishbone & Moncrief (1976) used this form even in the case of stationary spacetimes to obtain equilibria with constant ℓ*. However, we found it much more fruitful to use Eq. (9) in the general g_{tφ} ≠ 0 stationary spacetimes.
We obtain the Euler equation under the assumption of barotropicity and magnetotropicity as (15)
Once again, this equation quickly leads to an integrability requirement that either Ω_{,μ} = 0 or ℓ = ℓ(Ω). This result is known as the relativistic von Zeipel theorem (Abramowicz 1971).
However, finding closed analytical solutions is more complicated than in the static case; if we postulate the function ℓ(Ω), we obtain (16)
where we have used the fact that and we define L(Ω) analogously to the potentials Eq. (12)–(14), (17)
In other words, the thermodynamic and magnetothermodynamic potentials are specified by this procedure not only as functions of coordinates but also of an asofyet unspecified function Ω = Ω(r, ϑ).
To obtain Ω(r, ϑ) from the postulated ℓ(Ω) and thus find the full explicit solution, one has to use the relation (18)
which leads to the equation for Ω (19)
In general, this equation has to be solved numerically, see examples for nonmagnetized tori in Chakrabarti (1985), Daigne & Font (2004), Qian et al. (2009), Penna et al. (2013), and also more recent examples for magnetized tori in Wielgus et al. (2015), GimenoSoler & Font (2017). Nevertheless, if we instead want to obtain Ω(r, ϑ) as a closed analytical expression, we must impose some constraints on the form of ℓ(Ω). For instance, if Ω is to be a root of a polynomial equation of order n, it is easy to see that necessarily ℓ = P^{(n−1)}(Ω)∕P^{(n−1)}(Ω), where P^{(k)}(x) are some polynomials of order k.
We have chosen here to focus on the general form of ℓ(Ω) which leads to the expression for Ω as a root of a quadratic equation (n = 2) (20)
where ℓ_{0}, κ, λ are some constants.
We chose this particular form of the parametrization of ℓ(Ω) because the case κ = 0, λ = 0 corresponds to the wellknown ℓ = const. = ℓ_{0} Polish donuts of Kozłowski et al. (1978); Abramowicz et al. (1978), and, on the other hand, κ = 1, λ = 0 corresponds to ℓ* ≡ u_{φ}u^{t} = const. = ℓ_{0} of Fishbone & Moncrief (1976). The other choices of κ, λ, however, represent a continuous class connecting and extending these two particular solution families.
The expressions for Ω and the thermodynamic potential as functions of coordinates corresponding to Eq. (20) are given in Appendix B; examples of the obtained rotation curves and density profiles are discussed in the next section. In particular, it is shown that in Kerr spacetime the choice λ = 0, ℓ_{0} > 0 and a small κ > 0 generically leads to ℓ growing with radius, whereas κ < 0 to an angular momentum falling of with radius.
Even more specifically, we will see in the next section that κ < 0 leads to very thick tori with strongly nonKeplerian rotation profiles. However, since the solutions themselves do not have any direct physical meaning, and since they are used only as initial conditions for an MHD evolution that completely changes their character, we do not exclude any choice of parameters, save perhaps for those leading to outright pathological tori. In return, this allows us to potentially discover new states and behaviors of accretion disks.
Simulation parameters.
Initial condition parameters.
3 Numerical simulations
Because the solutions presented in the previous section are in closed form, it is easy to set them up as initial conditions for the evolution in Kerr spacetime in existing numerical codes and see whether they bring any variations in the properties of the accretion disks that evolve from them. Hence, we implemented our initial conditions in the 2D HARM code of Gammie et al. (2003); Noble et al. (2006) and conducted the numerical study in a “standard HARM setup”, as described for example in McKinney & Gammie (2004).
We placed the tori with various choices of the solution parameters near a Kerr black hole with spin a = 0.9375M, perturbed them with a small poloidal magnetic field with P∕b^{2} = 100 at the pressure maximum of the torus, and evolved the tori for a time period of 3000 M with a 512 × 512 resolution. As indicated already in Balbus & Hawley (1991, 1998), longterm MHD turbulence requires at least some magnetic field threading the disk vertically and this is the reason why the weak poloidal field is included in the initial conditions. Toroidal fields also trigger the MRI, but on their own they eventually decay (see also Wielgus et al. 2015; Fragile & Sądowski 2017). We have thus decided not to include a toroidal component of the magnetic field in the initial conditions () and leave the study of its effects to other works.
What follows is a description of some of the properties of the chosen initial conditions and a discussion of the simulation results. We used the default HARM outputs of energy, angular momentum, and restmass accretion rates to diagnose the evolutions. The parameters of the simulation are summarized in Table 1.
3.1 Properties of initial conditions
Our main target was to investigate the properties of the tori with the angularmomentum relation Eq. (20) in dependence on the various values of the parameter κ. We set λ = 0 and determine ℓ_{0} by requiring that the toroids have their inner and outer edge at r = 5 M and r = 12 M in Boyer–Lindquist coordinates respectively.
The thermodynamic relations were then determined by setting up the gas as initially isentropic, nonmagnetized (), and having a polytropic equation of state P = Kμ^{γ}, where μ is mass density and we set γ = 13∕9 to roughly mimic an idealgas mixture of relativistic electrons and mostly nonrelativistic ions. Furthermore, we used the HARM convention , where μ_{c} is the density of the toroid at the pressure maximum (see Appendix B).
We constructed 24 such initial conditions with κ in the interval [−1.11, 1.31]; a summary of the resulting parameters is given in Table 2, and the properties of the initial conditions can be inspected in Figs. 1–4.
We know from the analysis of vonZeipel radii by Abramowicz et al. (1978); Kozłowski et al. (1978) that the ℓ = const. toroids in the Kerr spacetime will have an angular velocity along the equator falling off with radius, at least as long as we are above the photon orbit. Since κ = 0, λ = 0 corresponds to the ℓ = ℓ_{0} = const. tori, we can then assume that for moderate magnitudes of κ our solutions will also have an angular velocity falloff, at least if we are far away from the photon sphere. This expectation is fulfilled in the case of our tori, as can be seen on a few examples in Fig. 1.
Furthermore, by examining the angular momentum relation Eq. (20) at λ = 0, ℓ_{0} > 0, Ω > 0 and by assuming that angular velocity is falling off with radius, we conclude that κ > 0 means angular momentum growing with radius, and a small κ < 0 (such that the expression Eq. (20) does not become singular) an angular momentum falloff. This is demonstrated in Fig. 2.
As for the dependence of the vertical structure of the toroid on the parameter κ, that seems to be impossible to characterize by any simple analytical argument. We can only offer the observation that the angular momentum distribution is such that for κ ≳ 1 the rotation curve in the equatorial plane is closer to Keplerian values and thus the disk should be closer to the canonical “thin disk” also inthe vertical structure. This can be seen to hold in Figs. 3 and 4. The κ ≳ 1 (λ = 0) disks are thus the ones closest to the state assumed to arise from longterm MHD simulations and, as such, could be preferred as initial conditions in astrophysical studies.
Last but not least, the HARM convention leads to the specific enthalpy h ≡ w∕μ to be equal to h = (2γ − 1)∕(γ − 1) = 4.25 at the pressure maximum of every torus. This fact leads to the relativistic Bernoulli parameter being much larger than one at the pressure maxima of the tori, a property common to hot accretion flows (Narayan & Yi 1994, 1995). The values of the Bernoulli parameter at the pressure maxima of the tori as a function of κ are plotted in Fig. 5.
Fig. 1 Profile of the angular velocity Ω inside the tori along the equatorial plane for different values of the parameter κ. In this and all the following graphs we assume λ = 0 and ℓ_{0} is determined by the constraint that the torus has its inner and outer edge at r = 5 M and r = 12 M respectively. 
Fig. 2 Profile of specific angular momentum ℓ inside the tori along the equatorial plane for different values of the parameter κ. 
Fig. 3 Vertical extent z = rcosϑ of the tori as a function of the parameter κ. 
3.2 Simulation results
The general scenario of the evolution of any of our tori is up to some intermittent episodes (see Sect. 3.3) qualitatively equivalent to the default HARM simulation with the ℓ* = const. tori as described for example in McKinney & Gammie (2004). We see the onset of the MRI and a rapid transition to turbulence throughout the disk starting from the inner edge. In parallel, the inner edge extends all the way toward the black hole, forming a plunging region, and the accretion disk reaches a quasisteady state with a magnetically dominated corona layer above it, and a nearly evacuated magnetized funnel around the poles.
The simulation time t = 3000 M is about ten orbital periods of the pressure maxima and it is sufficient to see the establishment of these quasistationary structures. However, it is not enough to see the ultimate fate of the disk in terms of the accumulation of magnetic flux around the black hole and a possible transition to a magnetically arrested disk, or the spreading of the disk to large radii through angular momentum exchange (for a long simulation see Narayan et al. 2012).
The mathematical analysis of Balbus & Hawley (1991) showed that the MRI grows with a rate which is dependent both on the wave vector of the disturbance and the local angularvelocity gradients in the fluid. To obtain a single numerical estimate of some kind of “typical strength” of the MRI for each torus, we have chosen to use the growth rate of the fastestgrowing MRI mode (maximum of the MRI growth rate in the perturbation wavevector space) Im (ω_{max}) evaluated at the location of the pressure maxima of the tori. Gammie (2004) showed that the rate of the fastest growing mode can be computed as half the shear rate in the case of relativistic Keplerian disks, and, even though this result has not been rigorously proven to apply to general relativistic flows, we have used it here as our estimate (for definition and computation of shear rates in Kerr spacetimes see e.g., Semerák 1998).
The resulting MRI growthrate estimates are given in Fig. 6, where we see that by variation of κ we get an MRI efolding time from 10 M to about 20 M. As can be seen on a few selected examples of accretion rates in Fig. 7, this leads to the first wave of matter arriving at the black hole at around 100–200 M. This certainly does not mean that the MRI needs approximately ten efolds to take effect because this time also includes the inspiral time required for the matter released by the instability to arrive to the black hole horizon. In fact, the waves arrive within a few tens of M away from each other, which suggests that matter is released from the tori within approximately five efolding times of the MRI.
This observation is further documented in Fig. 8, where we see that the average accretion rates over the first 500 M of evolution are correlated with the initial MRI estimates. This may be explained by noticing from Fig. 7 that the “first accretion waves” arrive earlier from the tori with higher shear (smaller κ), but also that these waves bring larger amounts of matter with them. However, we also see that the correlation starts to be broken for the mostshearing tori (with largest MRI growth rates), which is probably because they are already entering the nonlinear mode of evolution in the chosen averaging interval.
Furthermore, if we try to see whether the accretion rates in the last 1000 M of our evolution are correlated to the initial Im(ω_{max}), we see in the second part of Fig. 8 that there is no tangible correlation. This is because by that time the tori have evolved to the quasistationary “asymptotic^{3}” states wheretheir relation to the properties of initial conditions has mostly been washed out, perhaps up to the energy and angularmomentum content. For further details see the next subsection.
Fig. 4 Meridional sections of the density of the tori for different values of κ. The progression goes from geometrically very thick tori with large pressure gradients at their inner edge up to thinner toroids witha cusp (vanishing pressure gradient) gradually forming at their inner radius. 
Fig. 5 Relativistic Bernoulli parameter at the pressure maxima of the tori ( corresponds to matter with enough energy to escape to infinity) as a function of the parameter κ. In our case the thermodynamic properties are chosen such that h = 4.25 at the pressure maxima of the toroids, so the behavior of is due to the fact that the maxima are getting farther from the black hole with growing κ. 
Fig. 6 Growth rates of the fastestgrowing MRI mode Im(ω_{max}) evaluated at the position of the pressure maxima of the tori plotted as a function of the parameter κ. 
Fig. 7 Some examples of the numerically obtained mass accretion rate as a function of time. The height of the first peak of the accretion rate is correlated with κ, but the late evolutions are harder to relate to the parameters of the initial conditions. 
3.3 Advected specific energy and nominal efficiency
We also computed the “nominal efficiency” as was done in McKinney & Gammie (2004) but we found negative efficiencies for all of our tori at late times. This is not a mistake but, as will be seen from the following, a robust property of this type of simulations with an analytical underpinning. Hence, we must warn against the usage of as any sort of efficiency of the accretion process. Let us now discuss the arguments which explain this phenomenon.
The total energy in the disk is defined as (21)
where V is the spatial volume, and is the energy in the magnetic fields (initially almost zero). The total rest mass is, on the other hand, defined as (22)
The expression is just the coordinate mass density, so for instance is simply the massaveraged Bernoulli parameter.
Thus, we can see that the ratio corresponds to the average Bernoulli parameter with contributions to energy from the internal stress in the gas, and the magneticfield energy. A similar argument leads us to the realization that is the average Bernoulli parameter of the accreted fluid elements plus similar terms. However, since we start with a toroid with most of its elements having (see Fig. 5), we should not strictly expect the advection dominated flow to accrete exclusively the portion of the matter with and spit the rest out in an outflow. Thus, even when ignoring the pressure and magnetic contribution to energy, we should not be surprised by .
The simulation of McKinney & Gammie (2004) did not find the average to be larger than 1 but we believe that this is only thanks to the fact that their simulation halted at t = 2000M. In contrast, we found that for times t > 2000M the irregular advection mechanism is able to push even the highenergy fluid elements into the black hole. This is demonstrated in Fig. 9, where wesee that undergoes erratic “buildups” whose origin seems to be linked with the behavior of the magnetic fields and the coronal layer (see Fig. 10).
Specifically, visual inspection of the plots of the strength of magnetic fields reveals that the buildups in are associated with the coronal layer slowly raising all the way to the axis and filling the polar funnel; the abrupt drop in accreted specific energy is then associated with a sudden evacuation of the axis region and the reemergence of the funnel. This of course suggests that what is accreted in the “superenergetic” periods is mostly the magnetically dominated matter from the coronal region around the poles. On the other hand, it also suggests that this phenomenon might be associated with the strictly axisymmetric structure of the simulation and that it would not be reproduced in a 3D simulation. Alternatively, it could be caused by the small topological defect which is introduced in HARM to regularize the coordinate singularity at the axis (see Gammie et al. 2003).
One thing which is clear to us is the fact that this superenergetic accretion cannot be easily related to any parameter of the initial conditions, as we show in Fig. 11. All the properties of the initial conditions vary as smooth, typically monotonous functions of the parameter κ, so if the late is to depend on any simple property of the initial conditions, it should be obvious from its dependence on κ. Even though the late average can be recognized in Fig. 11 as a continuous function of κ in certain regions, the sampling is not sufficient to understand the structure any further. The only other observation, inferred from Figs. 2 and 11, is that the extremely high values of specific energy seem to be associated with initial rotation curves which are very far away from a Keplerian profile.
We believe that this dependence of on initial conditions should be studied in more detail in future works. Additionally, a more careful study should separate from the actually accreted Bernoulli parameter because the former involves a significant, or even dominant contribution of the energy of the magnetic field at late evolution times. However, implementing and running diagnostics that would allow us to make such a distinction and shed more light on this numerically observed phenomenon is beyond the scope of the current paper.
Fig. 8 Numerically obtained accretion rates of rest mass , energy and angular momentum averaged over the evolution time intervals [0, 1000 M] and [2000 M, 3000 M] as a functionof the initial MRI growth rate estimate Im(ω_{max}). We note that a lower negative value of corresponds to a higher loss rate of the quantity from the torus and thus a higher accretion rate. The displayed values correspond to the parameter κ in the interval [−1, 1.31]. 
Fig. 9 Accreted specific energy as a functionof time for chosen values of the parameter κ. 
Fig. 10 Snapshots of the structure of the magnetic fields for different tori evolutions plotted in Boyer–Lindquist coordinates. The color indicates , where is the maximal strength of the magnetic field in the respective image, the black circle is the interior of the black hole, and the right equatorial edge of the image is at r = 12M. The first twoimages on the left are snapshots for the κ = −0.48 torus before and after the first peak, at t = 1200 M and t =1700 M respectively (c.f. Fig. 9). The other two images on the right are snapshots for the “quiescent” κ = 1.2 torus at t = 1200 M and t =1700 M respectively. It appears that the behavior of is linked with the evolution of the coronal layer. 
Fig. 11 Accreted specific energy averaged over the interval t ∈ [2000 M, 3000 M] as a functionof κ. The vertical axis has its origin at 1, which means that all our flows return a negative nominal efficiency at late simulation times. 
4 Conclusions
We presented generalizations of the known closedform analytical solutions of geometrically thick fluid tori near black holes; our solutions are easy to construct and implement in numerical codes. Unlike the previously known solution families, our family provides a rich spectrum of possibilities for the rotation curve and geometrical shapes even if parameters such as the inner and outer radii are constrained. In particular, our κ ≳ 1 disks have semiKeplerian rotation profiles and could thus be preferred for initial conditions of astrophysical simulations rather than the ℓ = const. (Abramowicz et al. 1978; Kozłowski et al. 1978) or ℓ* = const. (Fishbone & Moncrief 1976) tori.
The MRI growth rates for different tori from our new class have different magnitudes and the tori exhibit different numerical evolutions accordingly. Namely, the average normalized accretion rate over the first 500 M of the simulation (three or more orbital periods) turns out to be linearly proportional to the growth rate of the fastestgrowing MRI mode at the pressure maximum of the torus.
On the other hand, the asymptotic states of the accretion disk which are achieved for an evolution time ≳ 2000 M do not exhibit any correlation of accretion rates with the initial MRI growth rates. Additionally, all of the asymptotic accretion disks tend to have the same qualitative diskcoronafunnel structure as described by for example McKinney & Gammie (2004) and the normalized accretion rates of their angular momentum, energy and mass vary within a single order of magnitude.
Nevertheless, we see that the disks vary quite wildly in the specific energy of the fluid they accrete; for particular choices of initial conditions and at the most extreme periods, the energy of the accreted elements approaches ten times their rest mass. This seems to be connected with the behavior of the strongly magnetized coronal layer because the extreme accretion episodes are associated with the corona pushing into the funnel, sometimes even up to the point of the short disappearance of the evacuated region. This behavior, however, appeared only for initial conditions far away from a Keplerian rotation profile. Further insights into this question should be obtained by future studies which should also involve more realistic matter models for the disk and better diagnostics ran during the simulation.
Acknowledgements
We would like to thank Claus Lämmerzahl and Volker Perlick for supervision and useful discussions on the topic. We are grateful for the support from a Ph.D. grant of the German Research Foundation (DFG) within Research Training Group 1620 “Models of Gravity”. Additionally, PJ kindly acknowledges the support from the Erasmus Mundus Joint Doctorate IRAP program.
Appendix A Deriving various forms of Euler equation
A.1 Simplifying
Omitting terms which are zero due to the stationarity and axisymmetry of all quantities, the circularity of velocity field, and toroidality of magnetic field, we obtain (A.1)
Let us now rewrite the a_{ν} and b_{ν;μ}b^{μ} terms as
We now realize that u_{μ} and are two normalized orthogonal vectors exclusively in the t − φ direction and we thus have (A.5)
where g_{αβ}_{(φt)} is the φ, trestriction of the metric (thatis, g_{rr}_{(φt)} = g_{θθ}_{(φt)} = 0 but otherwise the same as the metric). We then see that the magnetic part of the structural equations reads (A.6)
where in the last step we have used the formula for the derivative of a determinant of a matrix, and we denote . Using (A.7) we then obtain that the magnetic part of the stressenergy tensor can be written as (A.8)
A.2 Simplifying
Under the symmetry assumptions and circularity of the flow we obtain (A.9)
Let us now briefly derive the two possible forms of a_{ν} useful for the construction of analytical solutions.
A.2.1 FishboneMoncrief form of acceleration
The FishboneMoncrief form is obtained by expressing the acceleration in terms of the velocity in the ZAMO frame. The linear velocity of the flow in the ZAMO frame is , and the fourvelocity components of the circular flow are given as
These expressions are then directly substituted into the fouracceleration to obtain (A.12)
A.2.2 “Polishdonut” form of acceleration
The Polishdonut form of acceleration is derived by realizing that the acceleration can be rewritten as (A.13)
We now use the fourvelocity normalization to express u_{t} = −(1 + u^{φ}u_{φ})∕u^{t} and obtain (A.14)
The last step is to use the fourvelocity normalization to obtain the identity u^{t} u_{t} = −1∕(1 − Ωℓ) and thus u_{φ}u^{t} = ℓ∕(1 − Ωℓ) which leads to Eq. (9).
Appendix B Construction of tori
When we postulate ℓ = (ℓ_{0} + λΩ)∕(1 + κℓ_{0}Ω), we obtain the angular rotation frequency for κ≠1 as (B.1)
When we have ℓ = ℓ_{0} (such as when κ = λ = 0), the formerly quadratic equation for Ω in Eq. (19) becomes linear and the quadratic root (B.1) has a formal singularity. In that case we can either take an appropriate limit or directly substitute into Eq. (18) to obtain (B.2)
The function L(Ω) from Eq. (17) is easily integrated as (for a general choice of parameters) (B.3)
There are degeneracies in the parametrization when we allow for λ≠ 0 and we thus recommend to set λ = 0 for most practical purposes. Then the function L simplifies as (B.4)
The final expression for the thermodynamic potentials (for λ = 0) then reads (B.5)
where we of course need to substitute the λ = 0 version of (B.1) to obtain a purely coordinatedependent expression.
The expressions for mass and enthalpy density of an isentropic, nonmagnetized, and polytropic fluid P = Kμ^{γ} read (cf. Rezzolla & Zanotti 2013)
where W_{s} is the value of W at the surface of the torus (conventionally the inner edge).
The dynamical properties of the polytropic gas are invariant with respect to K (see e.g., Rezzolla & Zanotti 2013), so the choice of K is somewhat arbitrary. We choose K in the same way as the original HARM code and that is so that (B.8)
where μ_{c} is the maximum density^{4} in the torus (i.e., at the pressure maximum). This leads immediately to .
However, as was mentioned in the main text and can be seen from (B.6) and (B.7), the choice of K does influence the absolute values of specific enthalpy h = w∕μ throughout the fluid. Hence, some of the important physical characteristics of the toroids such as their total specific energy or total specific angular momentum depend on the value of K even at fixed spatial extent and rotation curves!
If we then fix κ (and λ = 0), and we want to confine the torus between some r_{in} and r_{out}, we must numerically solve the following equation for ℓ_{0} using either wellknown algorithms or software such as Mathematica or Matlab
The parameters which we obtained in the case of r_{in} = 5M, r_{out} = 12M are given in Table 2.
References
 Abramowicz, M. A. 1971, Acta Astron., 21, 81 [NASA ADS] [Google Scholar]
 Abramowicz, M. A., & Fragile, P. C. 2013, Living Rev. Relativ., 16, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Abramowicz, M., Jaroszyński, M., & Sikora, M. 1978, A&A, 63, 221 [NASA ADS] [Google Scholar]
 Anile, A. M. 1989, Relativistic Fluids and MagnetoFluids: With Applications in Astrophysics and Plasma Physics (Cambridge University Press) [Google Scholar]
 Balbus, S. A.,& Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A.,& Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
 Chakrabarti, S. K. 1985, ApJ, 288, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Daigne, F., & Font, J. A. 2004, MNRAS, 349, 841 [NASA ADS] [CrossRef] [Google Scholar]
 De Villiers, J.P., Hawley, J. F., & Krolik, J. H. 2003, ApJ, 599, 1238 [NASA ADS] [CrossRef] [Google Scholar]
 Dexter, J., Agol, E., & Fragile, P. C. 2009, ApJ, 703, L142 [NASA ADS] [CrossRef] [Google Scholar]
 Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962 [NASA ADS] [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]
 Gammie, C. F. 2004, ApJ, 614, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444 [Google Scholar]
 GimenoSoler, S., & Font, J. A. 2017, A&A, 607, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gourgouliatos,K. N., & Komissarov, S. S. 2018, MNRAS, 475, L125 [NASA ADS] [CrossRef] [Google Scholar]
 Komissarov, S. S. 2006, MNRAS, 368, 993 [NASA ADS] [CrossRef] [Google Scholar]
 Kozłowski, M., Jaroszyński, M., & Abramowicz, M. 1978, A&A, 63, 209 [NASA ADS] [Google Scholar]
 McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977 [NASA ADS] [CrossRef] [Google Scholar]
 Mewes, V., Font, J. A., Galeazzi, F., Montero, P. J., & Stergioulas, N. 2016, Phys. Rev. D, 93, 064055 [NASA ADS] [CrossRef] [Google Scholar]
 Montero, P. J., Font, J. A., & Shibata, M. 2010, Phys. Rev. Lett., 104, 191101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Narayan, R., & Yi, I. 1994, ApJ, 428, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., & Yi, I. 1995, ApJ, 444, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., Sądowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241 [NASA ADS] [CrossRef] [Google Scholar]
 Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna L. 2006, ApJ, 641, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Novikov, I., & Thorne, K. S. 1973, Black holes, 6, 343 [Google Scholar]
 Penna, R. F., Kulkarni, A., & Narayan, R. 2013, A&A, 559, A116 [NASA ADS] [CrossRef] [EDP Sciences] [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., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford University Press) [Google Scholar]
 Rezzolla, L., Baiotti, L., Giacomazzo, B., Link, D., & Font, J. A. 2010, Class. Quant. Grav., 27, 114105 [NASA ADS] [CrossRef] [Google Scholar]
 Ryan, B. R., Ressler, S. M., Dolence, J. C., et al. 2017, ApJ, 844, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Semerák, O. 1998, Gen. Rel. Grav., 30, 1203 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Tooper, R. F. 1965, ApJ, 142, 1541 [NASA ADS] [CrossRef] [Google Scholar]
 Wielgus, M., Fragile, P. C., Wang, Z., & Wilson, J. 2015, MNRAS, 447, 3593 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529 [NASA ADS] [CrossRef] [Google Scholar]
In many contexts the “Rayleigh–Taylor instability” is used exclusively for the instability of interfaces of fluids with different densities. Here, we have used the word as also applying to the similar effect in continuous media (see e.g., Gourgouliatos & Komissarov 2018).
The fluid fulfilling the condition P = P(w) is called barotropic by Kozłowski et al. (1978); Abramowicz et al. (1978), while other authors would rather call a fluid where the conditions P = P(μ) is fulfilled barotropic (see e.g., Tooper 1965). Both conditions coincide in flows with constant entropy per particle and in the Newtonian limit.
All Tables
All Figures
Fig. 1 Profile of the angular velocity Ω inside the tori along the equatorial plane for different values of the parameter κ. In this and all the following graphs we assume λ = 0 and ℓ_{0} is determined by the constraint that the torus has its inner and outer edge at r = 5 M and r = 12 M respectively. 

In the text 
Fig. 2 Profile of specific angular momentum ℓ inside the tori along the equatorial plane for different values of the parameter κ. 

In the text 
Fig. 3 Vertical extent z = rcosϑ of the tori as a function of the parameter κ. 

In the text 
Fig. 4 Meridional sections of the density of the tori for different values of κ. The progression goes from geometrically very thick tori with large pressure gradients at their inner edge up to thinner toroids witha cusp (vanishing pressure gradient) gradually forming at their inner radius. 

In the text 
Fig. 5 Relativistic Bernoulli parameter at the pressure maxima of the tori ( corresponds to matter with enough energy to escape to infinity) as a function of the parameter κ. In our case the thermodynamic properties are chosen such that h = 4.25 at the pressure maxima of the toroids, so the behavior of is due to the fact that the maxima are getting farther from the black hole with growing κ. 

In the text 
Fig. 6 Growth rates of the fastestgrowing MRI mode Im(ω_{max}) evaluated at the position of the pressure maxima of the tori plotted as a function of the parameter κ. 

In the text 
Fig. 7 Some examples of the numerically obtained mass accretion rate as a function of time. The height of the first peak of the accretion rate is correlated with κ, but the late evolutions are harder to relate to the parameters of the initial conditions. 

In the text 
Fig. 8 Numerically obtained accretion rates of rest mass , energy and angular momentum averaged over the evolution time intervals [0, 1000 M] and [2000 M, 3000 M] as a functionof the initial MRI growth rate estimate Im(ω_{max}). We note that a lower negative value of corresponds to a higher loss rate of the quantity from the torus and thus a higher accretion rate. The displayed values correspond to the parameter κ in the interval [−1, 1.31]. 

In the text 
Fig. 9 Accreted specific energy as a functionof time for chosen values of the parameter κ. 

In the text 
Fig. 10 Snapshots of the structure of the magnetic fields for different tori evolutions plotted in Boyer–Lindquist coordinates. The color indicates , where is the maximal strength of the magnetic field in the respective image, the black circle is the interior of the black hole, and the right equatorial edge of the image is at r = 12M. The first twoimages on the left are snapshots for the κ = −0.48 torus before and after the first peak, at t = 1200 M and t =1700 M respectively (c.f. Fig. 9). The other two images on the right are snapshots for the “quiescent” κ = 1.2 torus at t = 1200 M and t =1700 M respectively. It appears that the behavior of is linked with the evolution of the coronal layer. 

In the text 
Fig. 11 Accreted specific energy averaged over the interval t ∈ [2000 M, 3000 M] as a functionof κ. The vertical axis has its origin at 1, which means that all our flows return a negative nominal efficiency at late simulation times. 

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.