© 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 X-ray 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 so-called thin disk model (Shakura & Sunyaev 1973; Novikov & Thorne 1973; Abramowicz & Fragile 2013). However, when the accretion rate drops well below the so-called Eddington accretion rate Ṁ_Edd = 1.39 × 10¹⁸(M∕M_⊙)g s⁻¹, 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 magneto-hydrodynamic (MHD) fluid without radiation back-reaction (see e.g., Dexter et al. 2009). A more accurate model, however, requires the inclusion of radiative cooling and heating, and the partially two-temperature nature of the fluid in this accretion mode, see for example Ryan et al. (2017) and references therein. Consequently, a number of global general-relativistic MHD studies were conducted showing that the MRI-driven turbulence (Balbus & Hawley 1991, 1998) provides the angular momentum transport in the disk, and many features predicted by analytical or semi-analytical models emerge with only a small amount of additional ad-hoc 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¹ 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 space-time 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 time-scales. This instability was later shown to disappear when full self-gravity 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 space-times (e.g., in the Schwarzschild space-time) 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 space-time metric. A comma before an index defines a partial derivative with respect to the appropriate coordinate, whereas a preceding semi-colon 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, self-gravitation, 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 four-velocity, μ its rest-mass 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 magneto-fluid into an axially symmetric and stationary space-time with the metric (6)

Some useful symbols we use are the rotation frequency of zero-angular-momentum 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 radius-like quantity .

Now let us assume that all the properties of the magneto-fluid 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 particle-conservation 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 space-times

In static space-times (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² 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 right-hand side of Eq. (10) is a coordinate gradient, and the same must be true also for the left-hand 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 magneto-thermodynamic effective potentials respectively, Φ_c has clearly the meaning of a centrifugal potential, and Φ_(Z) is, in the currently discussed case of static space-times, 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.

  Analytically integrating the potentials from Eq. (12)–(14);

• 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 ideal-gas w = μ + 5μk_BT∕(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 space-times.

Nonetheless, we find it important to state these results explicitly because there is a number of studies that construct tori in static space-times and use somewhat complicated numerical methods tailored for stationary space-times, even though the above-stated general closed solution is available (e.g., Daigne & Font 2004; Narayan et al. 2012; Penna et al. 2013).

2.3 Stationary space-times

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 space-times. Furthermore, Fishbone & Moncrief (1976) used this form even in the case of stationary space-times to obtain equilibria with constant ℓ*. However, we found it much more fruitful to use Eq. (9) in the general g_tφ ≠ 0 stationary space-times.

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 magneto-thermodynamic potentials are specified by this procedure not only as functions of coordinates but also of an as-of-yet 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), Gimeno-Soler & 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⁽ⁿ⁻¹⁾(Ω)∕P⁽ⁿ⁻¹⁾(Ω), 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 ℓ₀, κ, λ are some constants.

We chose this particular form of the parametrization of ℓ(Ω) because the case κ = 0, λ = 0 corresponds to the well-known ℓ = const. = ℓ₀ 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. = ℓ₀ 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 space-time the choice λ = 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.

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 space-time 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 set-up”, 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² = 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), long-term 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 rest-mass 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 angular-momentum relation Eq. (20) in dependence on the various values of the parameter κ. We set λ = 0 and determine ℓ₀ 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 ideal-gas 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 von-Zeipel radii by Abramowicz et al. (1978); Kozłowski et al. (1978) that the ℓ = const. toroids in the Kerr space-time 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 ℓ = ℓ₀ = const. tori, we can then assume that for moderate magnitudes of κ our solutions will also have an angular velocity fall-off, 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 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 fall-off. 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 long-term 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 quasi-steady 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 quasi-stationary 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 angular-velocity 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 fastest-growing MRI mode (maximum of the MRI growth rate in the perturbation wave-vector 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 space-times see e.g., Semerák 1998).

The resulting MRI growth-rate estimates are given in Fig. 6, where we see that by variation of κ we get an MRI e-folding 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 e-folds 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 e-folding 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 quasi-stationary “asymptotic³” states wheretheir relation to the properties of initial conditions has mostly been washed out, perhaps up to the energy and angular-momentum 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 fastest-growing 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 mass-averaged 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 magnetic-field 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 high-energy fluid elements into the black hole. This is demonstrated in Fig. 9, where wesee that undergoes erratic “build-ups” 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 build-ups 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 “super-energetic” 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 axi-symmetric 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 super-energetic 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 closed-form 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 semi-Keplerian 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 fastest-growing 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 disk-corona-funnel 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 then obtain (A.4)

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 φ, t-restriction 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)

We now compute (A.7)

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 stress-energy 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 Fishbone-Moncrief form of acceleration

The Fishbone-Moncrief 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 four-velocity components of the circular flow are given as

These expressions are then directly substituted into the four-acceleration to obtain (A.12)

A.2.2 “Polish-donut” form of acceleration

The Polish-donut form of acceleration is derived by realizing that the acceleration can be rewritten as (A.13)

We now use the four-velocity normalization to express u_t = −(1 + u^φu_φ)∕u^t and obtain (A.14)

The last step is to use the four-velocity 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 ℓ = (ℓ₀ + λΩ)∕(1 + κℓ₀Ω), we obtain the angular rotation frequency for κ≠1 as (B.1)

When we have ℓ = ℓ₀ (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 coordinate-dependent 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⁴ 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 ℓ₀ using either well-known 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

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 fastest-growing 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