Free Access
Issue
A&A
Volume 559, November 2013
Article Number A116
Number of page(s) 8
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201219666
Published online 25 November 2013

© ESO, 2013

1. Introduction

Accretion flows onto black holes are typically magnetized and turbulent, so general relativistic magnetohydrodynamic (GRMHD) simulations have played an influential role in model building. The initial conditions for many of these simulations are the same: a rotating torus of fluid held together by gravity, pressure gradients, and centrifugal forces (Fishbone & Moncrief 1976; Fishbone 1977; Kozlowski et al. 1978; Abramowicz et al. 1978; Chakrabarti 1985). The entropy and angular momentum distribution of the torus are chosen arbitrarily and then hydrostatic equilibrium and an equation of state fix the fluid’s density, velocity, and pressure. At the start of a simulation, the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998) develops and the torus becomes turbulent. Turbulence transports angular momentum outward, allowing the fluid to accrete inwards, and the inner edge of the torus becomes an accretion flow. The torus typically persists throughout the simulation and provides a reservoir of fluid feeding the outer edge of the accretion flow.

Simulations based on the equilibrium torus solutions of Fishbone & Moncrief (1976) and Chakrabarti (1985) have found many applications: thin disk models (Shafee et al. 2008; Noble et al. 2009; Penna et al. 2010) black hole spin evolution (Gammie et al. 2004), radio emission from Sgr A* (Noble et al. 2007; Dibi et al. 2012; Shcherbakov et al. 2012), black hole jets (McKinney 2006, 2005; Nagataki 2009; Tchekhovskoy et al. 2011; Tchekhovskoy & McKinney 2012), computations of spectra (Hilburn et al. 2010), magnetized accretion with neutrino losses (Barkov & Baushev 2011; Shibata et al. 2007; Barkov 2008), pair production in low luminosity galactic nuclei (Mościbrodzka et al. 2011), numerical convergence studies (Shiokawa et al. 2012), tilted disk evolution (Fragile et al. 2007; Fragile & Blaes 2008; Henisey et al. 2012), binary black hole mergers (Farris et al. 2011, 2012), and magnetically arrested disks (De Villiers et al. 2003; McKinney et al. 2012).

The equilibrium tori of Fishbone & Moncrief (1976) and Chakrabarti (1985) are designed to be simple. They assume unphysical, power law angular momentum distributions in order to keep the solutions analytical. When they are used as the initial condition for GRMHD simulations, one hopes the turbulent accretion flow “forgets” unrealistic features of the initial torus. However, this does not always seem to be the case. For example, the Bernoulli parameter of the initial torus appears to persist through to the final accretion flow (Fig. 1).

thumbnail Fig. 1

Initial (dotted) and final (solid) midplane Bernoulli parameters, Be, of two GRMHD accretion simulations. The Bernoulli parameter is measured in units of rest mass energy. The energetically unbound flow (blue) is the A0.0BtN10 model of McKinney et al. (2012). The energetically bound flow (black) is the ADAF/SANE model of Narayan et al. (2012). The unbound simulation reached t = 96   796M and the bound simulation reached t = 200   000M. In both cases, the final Be appears to depend on the initial Be.

Open with DEXTER

The Bernoulli parameter is the sum of the kinetic energy, potential energy, and enthalpy of the gas (at least in Newtonian dynamics, where this splitting can be made precise. See Sect. A.2 for a discussion of the GRMHD Bernoulli parameter.) At large distances from the black hole, the potential energy vanishes. Since the other two terms are positive, gas at infinity must have Be ≥ 0. Furthermore, in steady state and in the absence of viscosity, Be is conserved along streamlines. Hence any parcel of gas that flows out with a positive value of Be can potentially reach infinity. A flow with positive Be is called unbound and a flow with negative Be is called bound. Unbound flows are more likely to generate outflows.

As shown in Fig. 1, the Bernoulli parameter of the accretion flows in some GRMHD simulations appears to be set by the initial torus. This is a concern, as it suggests the strength of simulated outflows, such as winds, might be sensitive to arbitrary choices in the initial conditions. GRMHD simulators should choose the initial Be with some care. For example, the true Bernoulli parameter of an accretion flow onto a supermassive black hole is probably related to its value at the outer edge of the accretion flow, where Be ~ 0. So it would be reasonable to simulate this flow with an initial torus with Be ~ 0. However, in the earlier solutions of Fishbone & Moncrief (1976) and Chakrabarti (1985), the Bernoulli parameter is tied to other parameters of the torus, such as its thickness, which one would like to vary independently. The solution in this paper is more flexible and allows for varying the Bernoulli parameter and torus thickness independently.

The Bernoulli parameter is not the only arbitrary feature of GRMHD initial conditions that may affect the final results. It is known that GRMHD simulation results depend on the initial magnetic field strength and topology (Beckwith et al. 2008; Penna et al. 2010; Tchekhovskoy et al. 2011; McKinney et al. 2012). The rotation rate of the initial torus may also be important. For example, simulations which start from slowly rotating tori will tend to be more convectively unstable than simulations which start from rapidly rotating tori, as rotation tends to stabilize accretion flows against convection.

Our paper is organized as follows. In Sect. 2, we construct our equilibrium torus solution. In Sect. 3, we obtain approximate analytical formulae for the radius of the outer edge, the radius of the pressure maximum, the Bernoulli parameter, and the geometrical thickness of the torus. In Sect. 4, we summarize our results. In Appendix A, we describe a magnetic field configuration for the torus and discuss the Bernoulli parameter of the magnetized torus. The magnetic field consists of multiple poloidal loops and is constructed so that the magnetic flux and gas-to-magnetic pressure ratio are the same in each loop. This setup could be useful for GRMHD simulations.

2. New torus solution

We consider a fluid torus in hydrostatic equilibrium around a Kerr black hole. We assume the flow is axisymmetric and stationary. We use an ideal gas equation of state and assume the fluid has a constant entropy. Self-gravity is ignored.

We work in Boyer-Lindquist coordinates (t,r,θ,φ) and units for which G = c = 1. The solution is a function of r and θ only (i.e., it is stationary and axisymmetric). There are six adjustable parameters. They are defined below but may be summarized here: rin is the radius of the inner edge of the torus, r1 and r2 control the shape of the angular momentum distribution, ξ controls the rotation rate, κ sets the entropy, and Γ is the adiabatic index.

The components of the Kerr metric we need are where M and a are the mass and spin of the black hole and Σ = r2 + a2cos2θ.

Let uμ be the fluid four-velocity. We assume the fluid’s angular momentum density,  ≡ uφ/ | ut |, is constant on von Zeipel cylinders (Abramowicz 1971). The radius, λ, of a von Zeipel cylinder with angular momentum is (Chakrabarti 1985): (3)In Newtonian gravity, λ = rsinθ, the usual cylindrical radius. In the Schwarzschild metric, . In the Kerr metric, one must solve Eq. (3)numerically.

The angular momentum density of circular, equatorial orbits in the Kerr metric is (Novikov & Thorne 1973) (4)where, We choose the angular momentum distribution of the fluid torus to be: (5)There are three regions. In the inner and outer regions of the torus, the angular momentum density is a constant independent of radius. The size of these regions is set by λ1 and λ2. At intermediate radii, the angular momentum density is a fraction ξ of the Keplerian distribution (4). This is chosen to be close to the expected, sub-Keplerian angular momentum distribution of a real accretion flow. It is more physical than the power law angular momentum distributions of the Fishbone & Moncrief (1976) and Chakrabarti (1985) solutions.

The angular velocity of the torus is (6)We have assumed ur = uθ = 0, so this fully determines the velocity.

Given the velocity of the torus, we can determine its density and pressure from the Euler equation. Let (7)which sets the gravitational force felt by the fluid. The Euler equation (Abramowicz et al. 1978) is then (8)Our notation is standard: ρ0, U, and p are the mass density, internal energy, and gas pressure of the fluid in its rest frame. Euler’s equation describes the balance between pressure gradients (LHS) and gravitational and centrifugal forces (RHS) required for hydrostatic equilibrium.

To solve the Euler equation, it is helpful to introduce the effective potential (Kozlowski et al. 1978) (9)where, (10)The lower limit of the integral, rin, is the radius of the inner edge of the torus. The upper limit, rλm(r,θ), is the equatorial radius of the Von Zeipel cylinder containing (r,θ). For example, rλm(r,π/2) = r.

In terms of the effective potential, the Euler equation is (11)We can compute the effective potential because it depends only on . So the RHS is known. The boundary of the torus is the isopotential surface W(r,θ) = Win ≡ W(rin/2).

The specific enthalpy is w = 1 + ϵ + p/ρ0, where ϵ = U/ρ0 is the specific internal energy1. For an isentropic torus, the Euler equation can be integrated to obtain (Kozlowski et al. 1978) (12)We assume the equation of state p = ρ0   ϵ   (Γ − 1), so the specific internal energy is (13)The rest mass density and pressure are The entropy, κ, is a free parameter. The torus is now fully determined.

To summarize, we first choose the angular momentum distribution (5). The angular momentum distribution determines the angular velocity and effective potential. The effective potential determines the enthalpy of the torus through Euler’s equation. Fixing an ideal gas equation of state and assuming an isentropic torus, gives the density, pressure, and internal energy. There are six free parameters: the radius of the inner edge of the torus, rin, the break radii in the angular momentum distribution, r1 and r2, the normalization of the angular momentum, ξ, the entropy, κ, and the adiabatic index, Γ. The entropy simply sets the density scale (Eq. (15)) and, as there is no self-gravity, it has no effect on the dynamics.

3. Approximate analytical formulae

The solution of Sect. 2 can be implemented in GRMHD codes numerically. However, for physical understanding, it is useful to have approximate formulae that describe the torus analytically. In this section, we obtain approximate analytical formulae for the outer edge, pressure maximum, geometrical thickness, and Bernoulli parameter of the torus.

thumbnail Fig. 2

Counterclockwise from the upper left panel, we show the radius of the outer edge of the torus, the radius of the pressure maximum, the Bernoulli parameter at the pressure maximum, and the geometrical thickness at the pressure maximum as a function of the rotation parameter, ξ, for several choices of rin/M: 8 (blue), 9.993 (green), 10 (red), 10.01 (gray), and 12 (yellow). The pressure maximum position, rmax/M, is independent of rin. The horizontal scale of the upper left panel is different from the other panels because rout is only finite when Be < 0.

Open with DEXTER

The most complicated feature of the exact solution is the integral in Eq. (10). To obtain approximate analytical formulae for the torus, we need to simplify this integral. Let us first rewrite it as (16)where we have used dΩ = d(Ω) − Ωd.

In the inner region of the torus (λ < λ1), the angular momentum is constant and the integrand in Eq. (16)vanishes. So (17)In the outer region of the torus (λ > λ2), we may approximate the integral by plugging the Newtonian formula into Eq. (5), and using the Newtonian angular velocity Ω(λ) = (λ)/λ2. Now integrating from λ1 to λ2, we obtain (18)where, (19)Equation (18)is approximate because we ignored special relativistic contributions to and Ω. But these are small in the outer regions of the torus. We can use our analytical approximation of F to obtain simple formulae describing the torus.

3.1. Radius of the outer edge

The boundary of the torus is the isopotential surface W(r,θ) = Win = W(rin/2). We can simplify W at the outer edge by using a Newtonian description there. The equation for the outer radius, rout, becomes (20)The solution for the outer radius is (21)where (22)Figure 2 (top left panel) shows the variation of rout with ξ for several choices of rin and fixed r1 = 42M and r2 = 1000M. The outer radius increases as the inner radius decreases because the boundary of the torus is moving to larger isopotential surfaces. In the region of parameter space shown in Fig. 2, the Bernoulli parameter is small and negative, as we will see below. When the Bernoulli parameter is small and negative, the outer radius is very sensitive to ξ.

Solutions with positive Bernoulli parameter have rout at infinity. They are unbound and not very physical. In fact, they are not even tori, they are infinite cylinders kept confined by a nonzero pressure at infinity. Unphysical, infinite cylinders based on the Fishbone & Moncrief (1976) and Chakrabarti (1985) solutions have at times been used as the initial condition for GRMHD simulations.

3.2. Pressure maximum

The pressure gradient in the Euler Eq. (8)is zero at the pressure maximum. So the fluid must move along geodesics there. In other words, the pressure maximum is where (23)where K is the angular momentum of circular, equatorial geodesics (Eq. (4)) and rmax is the radius of the pressure maximum. For sub-Keplerian flow (ξ < 1), the pressure maximum must be in the inner region of the torus (λ < λ1), because the angular momentum is strictly sub-Keplerian for λ > λ1.

thumbnail Fig. 3

Left column: log   ( − Be) and | h/r | for the new equilibrium torus solution described in the text. The Bernoulli parameter is shown as a function of (r,θ) and as a function of r at the midplane. The torus is energetically bound, i.e., Be < 0. Right column: log   (Be) and | h/r | for a solution of Chakrabarti (1985) which is energetically unbound, i.e., Be > 0. Note that the unbound torus is actually thinner than the bound torus in this example.

Open with DEXTER

Equation (23)gives (24)where − 4M is the leading order relativistic correction. Figure 2 (lower left panel) shows the dependence of rmax on ξ for λ1 = 42M. In the Keplerian limit ξ → 1, the pressure maximum approaches λ1. Lowering ξ moves the pressure maximum toward the inner edge of the torus.

3.3. Bernoulli parameter

The relativistic Bernoulli parameter is (Novikov & Thorne 1973) (25)Tori with Be > 0 are energetically unbound and tori with Be < 0 are energetically bound. In the inner region of the torus (λ < λ1), the Bernoulli parameter is: (26)Figure 2 (lower right panel) shows how Be(rmax) depends on ξ for several choices of rin and fixed break radii r1 = 42M and r2 = 1000M. The outer edge of the torus goes to infinity as Be → 0 (from below).

3.4. Thickness

The boundary of the torus is the isopotential surface W(r,θ) = Win. In the inner region of the torus (λ < λ1), the surface of the torus is (27)θ is measured from the polar axis, so the scale height of the torus is Θ = π/2 − θ. Figure 2 (upper right panel) shows Θ(rmax) as a function of ξ for several choices of rin.

3.5. Bernoulli parameter and thickness

Our solutions tend to be more bound than earlier solutions. We show this with an example. Figure 3 (upper left panel) shows the Bernoulli parameter of one of our solutions in the (r,θ) plane2. The torus is energetically bound: Be < 0. In the same figure (upper right panel), we show the Bernoulli parameter of a solution of Chakrabarti (1985)3. This torus is energetically unbound.

In the bottom panels of Fig. 3, we show the density scale heights, (28)of both solutions. Thinner tori tend to be more bound than thicker tori, but in this example the unbound torus is actually thinner than the bound torus. In other words, our solutions tend to be more bound than earlier solutions. This is an advantage, because unbound solutions are unphysical as initial conditions for GRMHD simulations. They are infinitely extended cylinders supported by pressure at infinity.

A second lesson to draw from this example is that the torus thickness, | h/r |, does not fix even the sign of the Bernoulli parameter (much less its value). GRMHD simulators typically focus on the thickness of the initial torus but this may not be enough if, for example, the Bernoulli parameter influences the strength of accretion disk winds. In other words, it may be important to know the thickness and the Bernoulli parameter of the initial torus independently.

4. Conclusions

We have constructed a new equilibrium torus solution. It provides an alternative initial condition for GRMHD disk simulations to the earlier solutions of Fishbone & Moncrief (1976) and Chakrabarti (1985). The angular momentum density is constant in the inner and outer regions of the torus and follows a sub-Keplerian distribution in between. The entropy is constant everywhere.

The Bernoulli parameter, rotation rate, and geometrical thickness of the solution can be varied independently. The torus tends to be more energetically bound than earlier solutions, for the same thickness | h/r |. In fact, we have shown that it is possible to generate equilibrium tori with the same | h/r |, but for which one is bound and the other unbound. So as GRMHD initial conditions, the new solutions might lead to weaker (or non-existent) disk winds. Future GRMHD simulations will need to explore how the simulation results depend on the Bernoulli parameter, rotation rate, and geometrical thickness of the initial equilibrium torus.

In real accretion flows, the Bernoulli parameter is probably set at the outer edge of the flow. If this is far from the black hole (for example, at the Bondi radius), then one expects Be ~ 0. It is not possible to obtain converged GRMHD simulations out to the Bondi radius. One would need to run the simulations for a timescale of order the Bondi radius viscous time, which is impractical. Computational resources limit the duration of even the longest run simulations to t ~ 200   000M, which corresponds to the viscous time at r ~ 200M. Probably the best one can hope to do is choose an initial condition with a realistic Bernoulli parameter (Be ~ 0) at the outset. The solution in this paper is flexible enough to allow independent control over the thickness and Bernoulli parameter of the torus, so it is an ideal initial condition for GRMHD simulations. In fact, while this paper was in preparation, there appeared several simulations based on our solution (Narayan et al. 2012; Sa¸dowski et al. 2013; Penna et al. 2013a,b). We expect more soon.


1

We caution that ϵ is used for two different but related concepts in the literature. In older papers, such as Kozlowski et al. (1978), ϵ is the total energy density, ρ0 + U. In more recent literature, such as De Villiers et al. (2003), ϵ is the specific internal energy, U/ρ0. We follow the latter convention.

2

The free parameters are rin = 10M, r1 = 42M, r2 = 1000M, κ = 0.00766, and ξ = 0.708.

3

The free parameters are rin = 10M, κ = 0.00136M, dlog /dlog λ = 0.4, and pressure maximum at rmax = 40.

Acknowledgments

R.F.P was supported in part by a Pappalardo Fellowship in Physics at MIT.

References

Appendix A: Adding a magnetic field

Implementing the new equilibrium torus as a GRMHD initial condition requires adding a magnetic field to the torus. Here we record one possible magnetic field, a series of poloidal magnetic loops. We discuss the Bernoulli parameter of the magnetized torus in Sect. A.2.

Appendix A.1: Magnetic field solution

We construct the magnetic field so that each loop carries the same magnetic flux and β = pgas/pmag is roughly constant. Simulations often require initial conditions that minimize secular variability during the run, so these features can be useful.

Three free parameters appear in the solution: rstart, rend, and λB. The first two set the inner and outer boundaries of the magnetized region and the third controls the size of the poloidal loops.

We define the field through the vector potential, Aμ. The magnetic loops are purely poloidal, so only Aφ is nonzero. To keep β close to a constant, the magnetic field strength tracks the fluid’s internal energy density. Define (A.1)where, The function q is defined to give q = 1 at the midplane and q → 0 away from the midplane. The factor sin3θ smooths the vector potential as it approaches the edges of the torus. Dropping this factor leads to a torus with highly magnetized edges.

Further define (A.4)The vector potential is then The sinusoidal factor in Aφ breaks the poloidal field into a series of loops. The function f(r) gives each loop the same magnetic flux. The number of loops is controlled by λB. The overall normalization of Aφ has not been specified so it can be tuned to give any field strength.

We give an example in Fig. A.1. The equilibrium torus is as in Fig. 3 and the field parameters are rstart = 25, rend = 550, and λB = 15/4. In this example there are eight magnetic loops. The magnetic flux, Aφ, peaks at the center of each loop, measures the flux carried by the loop, and is the same across the torus. The magnetization β = pgas/pmag peaks at loop edges and drops at loop centers, but is roughly constant across the torus. In this example β ~ 100.

thumbnail Fig. A.1

Top panels: mid-plane density and the enclosed magnetic flux as a function of radius. Lower panels: density and the magnetization parameter β of the torus in the poloidal plane. Each loop carries the same magnetic flux. The magnetization is roughly constant throughout the torus.

Open with DEXTER

Appendix A.2: GRMHD Bernoulli parameter

The stress energy tensor of a magnetized fluid is where bμ is the magnetic field in the fluid’s rest frame and hμν = gμν + uμuν is the projection tensor. The magnetic field contributes b2/2 to the total internal energy and b2/2 to the total pressure, and introduces a stress term, − bμbν.

The Euler equations, h·(∇·T) = 0, become: where a = ∇uu is the fluid’s acceleration. We have used ∇·b = 0 to simplify the last term on the RHS.

Assume the flow is stationary and adiabatic, project the Euler equations along ξ = t, and combine terms using the first law of thermodynamics. This leads to: For the field configuration of Sect. A.1, b is purely poloidal and ξ   h is purely toroidal. So the RHS is zero.

We thus obtain a straightforward generalization of Eq. (25): (A.7)The unmagnetized torus has Be ~  − wgas. Adding magnetic fields to the torus, with gas-to-magnetic pressure ratio β = pgas/pmag, changes the Bernoulli parameter by terms of order 1/β. Simulations typically have initial β ~ 100, in which case the magnetic contribution to the Bernoulli parameter is of order 1%.

All Figures

thumbnail Fig. 1

Initial (dotted) and final (solid) midplane Bernoulli parameters, Be, of two GRMHD accretion simulations. The Bernoulli parameter is measured in units of rest mass energy. The energetically unbound flow (blue) is the A0.0BtN10 model of McKinney et al. (2012). The energetically bound flow (black) is the ADAF/SANE model of Narayan et al. (2012). The unbound simulation reached t = 96   796M and the bound simulation reached t = 200   000M. In both cases, the final Be appears to depend on the initial Be.

Open with DEXTER
In the text
thumbnail Fig. 2

Counterclockwise from the upper left panel, we show the radius of the outer edge of the torus, the radius of the pressure maximum, the Bernoulli parameter at the pressure maximum, and the geometrical thickness at the pressure maximum as a function of the rotation parameter, ξ, for several choices of rin/M: 8 (blue), 9.993 (green), 10 (red), 10.01 (gray), and 12 (yellow). The pressure maximum position, rmax/M, is independent of rin. The horizontal scale of the upper left panel is different from the other panels because rout is only finite when Be < 0.

Open with DEXTER
In the text
thumbnail Fig. 3

Left column: log   ( − Be) and | h/r | for the new equilibrium torus solution described in the text. The Bernoulli parameter is shown as a function of (r,θ) and as a function of r at the midplane. The torus is energetically bound, i.e., Be < 0. Right column: log   (Be) and | h/r | for a solution of Chakrabarti (1985) which is energetically unbound, i.e., Be > 0. Note that the unbound torus is actually thinner than the bound torus in this example.

Open with DEXTER
In the text
thumbnail Fig. A.1

Top panels: mid-plane density and the enclosed magnetic flux as a function of radius. Lower panels: density and the magnetization parameter β of the torus in the poloidal plane. Each loop carries the same magnetic flux. The magnetization is roughly constant throughout the torus.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.