# The general relativistic thin disc evolution equation

###### Abstract

In the classical theory of thin disc accretion discs, the constraints of mass and angular momentum conservation lead to a diffusion-like equation for the turbulent evolution of the surface density. Here, we revisit this problem, extending the Newtonian analysis to the regime of Kerr geometry relevant to black holes. A diffusion-like equation once again emerges, but now with a singularity at the radius at which the effective angular momentum gradient passes through zero. The equation may be analysed using a combination of WKB, local techniques, and matched asymptotic expansions. It is shown that imposing the boundary condition of a vanishing stress tensor (more precisely the radial-azimuthal component thereof) allows smooth stable modes to exist external to the angular momentum singularity, the innermost stable circular orbit, while smoothly vanishing inside this location. The extension of the disc diffusion equation to the domain of general relativity introduces a new tool for numerical and phenomenolgical studies of accretion discs, and may prove to be a useful technique for understanding black hole X-ray transients.

###### keywords:

accretion, accretion discs — black hole physics — turbulence^{†}

^{†}pagerange: The general relativistic thin disc evolution equation–The general relativistic thin disc evolution equation

^{†}

^{†}pubyear: 2017

## 1 Introduction

In thin disc accretion theory, the constraints of angular momentum and mass conservation may be combined into a single evolutionary equation for the disc surface density, a classic result first emphasised and discussed at length by Lynden-Bell & Pringle (1974) (see also the review of Pringle 1981). In its original implementation, the equation takes the form of a diffusion equation, with a diffusion
coefficient proportional to an ad hoc turbulent viscosity. Balbus & Papaloizou (1999) later showed
how the same evolution equation emerges without the need to introduce an explicit viscosity. By writing the velocity field as a sum of a mean plus a fluctuation (with vanishing mean), an effective diffusion coefficient emerges which is proportional to the correlation in the radial and azimuthal velocity fluctuations^{1}^{1}1The presence of a magnetic field
can be incorporated into this formalism, with the product of the radial and azimuthal Alfven velocities subtracted from the kinetic velocity fluctuations in the diffusion coefficient (resulting in an additional positive stress)..

The evolutionary equation has heretofore been used in the regime of Newtonian gravity (e.g. Pringle 1981). Solutions of the equation show that matter in accretion discs drifts inward, while angular momentum is transported outward, sustained by a vanishingly small mass fraction of the disc. The extension of the evolutionary equation to include general relativistic gravity has not yet been done, and it is not without interest. It is the purpose of this paper to derive and analyse the general relativisitic version of the thin disc evolutionary equation. We present a very general global asymptotic analysis (assuming small stress and/or rapid modal time scales), which can be equally well applied in the Newtonian limit.

The inner regions of neutron star and black hole discs are dynamically complex. It would be naive to apply simple thin disc dynamics uncritically. The formal thin disc problem is nevertheless quite interesting, first as an illustration of how the diffusion dynamics breaks down at the innermost stable circular orbit (ISCO) of the disc, second of how the diffusion equation extends to ISCO-free Kerr orbits in general relativity, and third as a useful analytical tool for understanding numerical simulations. It is especially noteworthy that while the effective diffusion coefficient of the disc equation becomes singular at the ISCO, the solution is nevertheless mathematically well-behaved. The global normal modes include exponentially growing modes confined to the zone within the ISCO, which completely disrupt the interior disc structure, leaving the outer disc intact. This is in accord with numerical simulations.

The plan of this paper is as follows. In §2 we first derive the form of the disc evolution equation that follows from the conservation of particle number and the azimuthal component of the stress energy tensor. This reduces to the Lynden-Bell—Pringle (1974) equation in the Newtonian limit. A solution of the general equation is presented in §3 for modes with exponential time dependence, using WKB, local analysis, and matched asymptotic expansions. The cases of both finite and vanishing stress at the location of the ISCO are presented, and we argue that thin discs will evolve to a state at which the vanishing stress boundary condition is achieved. We use the modal solutions to construct a general Green’s function solution. Finally, in §4 we summarise the presentation. This scope of this paper is to present a mathematical treatment of the equation. Astrophysical applications will be explored in a separate study.

We observe the following conventions. The speed of light is set to unity throughout this work. Greek indices generally denote spacetime coodinates. The exception is , which is always the azimuthal angular coordinate. The time coordinate is labelled . The metric in local inertial coordinates is Other notation is standard: is the gravitational constant, the central black hole mass, and the gravitational radius.

## 2 Fundamental equations

### 2.1 Conserved fluxes

The two conserved quantities of interest are the particle number current , where is the rest frame number density and the contravariant 4-velocity, and the azimuthal component of the stress energy tensor . Following Page & Thorne (1974), we will work in “cylindrical Boyer-Lundquist” coordinates in the Kerr metric, ultimately using the equations in their height-integrated form. This involves ignoring the higher order curvature terms of order near the equatorial plane. We assume that the disc is axisymmetric and thin. The conservation equation of particle number before height integration is simply

(1) |

where the semi-colon denotes a covariant derivative. If we now integrate over , and assume that the velocities are independent of height, the remaining 4-velocity components are , and . With rest mass per particle , and the integrated column density

the particle conservation equation for an axisymmetric disk is

(2) |

where , the absolute value of the metric tensor determinant.

The equation for conservation of the stress-energy tensor is , with taking the form of an ideal fluid plus a contribution from the radiation field, denoted :

(3) |

The radiation stress is given by

(4) |

where is the radiative energy flow vector, which satisifies (Page & Thorne 1974). We ignore (for the moment) the contribution of any additional stress tensor that may be present. The angular momentum carried by radiated photons is not negligible when rotational velocities are of order the speed of light (Novikov & Thorne 1973). Then, assuming axisymmetry , and defining

(5) |

the conservation equation becomes

(6) |

Here is the rest energy density (including in principle a thermal contribution, which is however ignored in the thin disk limit), is the thermal pressure (which shall likewise be ignored), and is the affine connection. For axisymmetric metrics this is:

(7) |

Therefore, for any symmetric tensor , the combination is

(8) |

since the metric derivatives are antisymmetric in and while is symmetric^{2}^{2}2Knowledgeable readers will recognise in equation (8) a Killing vector calculation. I thank C. Gammie for drawing my attention to this point..
By contrast, is not independent of , so that

(9) |

a result we use below. Equation (6) now reduces to:

(10) |

The disc turbulence is represented by writing the 4-velocity as a mean flow plus a fluctuation with vanishing mean, . In particular,

(11) |

The asymptotic scalings of the fluctuations satisfy:

(12) |

i.e., the orbital velocity and angular momentum are much larger than their associated fluctuations, and the inward mean radial drift velocity is yet an asymptotic order smaller than the fluctuations in either the radial velocity or orbital velocity. The two fluctuations are assumed to be of comparable order, suitably dimensionalised. In common with Newtonian theory (Balbus & Papaloizou 1999), we expect (the product of a zeroth order rotational velocity and a second order radial drift) to be of the same asymptotic order as (the product of two first order fluctuations). As always, it is important to distinguish contravariant (angular 4-velocity) from covariant (angular 4-momentum):

where we have ignored as negligibly small.

### 2.2 “Stress by strain” and radiation

#### 2.2.1 Equilibrium models

For the equilibrium models under consideration, Page & Thorne (1974) present a relationship between the disc shear, a tensor coupling like viscosity, and the energy radiated from its surface. In our notation, this relation reads:

(13) |

where

(14) |

is the angular velocity measured by an observer at infinity, and is the radiated energy flux in the local rest frame. In essence, this states that the energy extracted from differential rotation and put into turbulent fluctuations is locally radiated away at the same rate. We will make use of this relation in §2.3 below, which also holds in our case because of the assumption that the thermal timescale is more rapid than the evolutionary timescale. It is of some technical interest to revisit this important relationship in more detail in an out-of-equilibrium context, which we do in the following section (see also Balbus & Hawley 1998, Balbus & Papaloizou 1999). The reader willing to adopt equation (13) directly may wish to skip directly to §2.3 below, without loss of continuity.

#### 2.2.2 Free energy from shear

The radial conservation equation is given, with the help of equation (9) and particle number conservation, by

(15) |

We have multiplied by with the aim of assembling a fluctuation energy equation; the radial component of the radiation stress is small, but retained here to maintain a covariant formulation.

Next, we write the velocities as a mean plus fluctuating . The largest contributions from the final term of equation (15) comprise the equilibrium solution and are not of interest; they cancel out. Retaining the next largest group of terms, our equation becomes

(16) |

where, in the final term, we have used symmetry. It is convenient for now to retain in the divergence term without separating its mean and fluctuation. (The radiation stress will likewise contain a fluctuating component, which is not shown explicitly.) The term involving is an asymptotic order smaller than the others, and may be dropped. We arrive at:

(17) |

Next, the equation for (angular momentum conservation) is

(18) |

Using particle number conservation and remembering that depends only upon , this becomes

(19) |

Following the thin disc Novikov and Thorne (1973) models, the radiation flux is assumed to be dominated by its vertical component, and in particular by the term. We rewrite the last term to obtain:

(20) |

The equation is handled similarly:

(21) |

which expands to

or

(22) |

The final equation is simple and straightforward, as there is by assumption no mean flow:

(23) |

We now sum over equations (16), (20), (22) and (23) to obtain, after some algebra and index shifting:

(24) |

The final step is to use , which gives

Using this in (24), averaging to form and collecting terms, we obtain

(25) |

where

(26) |

is just the relativistic analogue of the shear gradient (e.g. Page and Thorne 1974).

Equation (25) is a relationship for the rate at which stress extracts energy from the shear, involving via the first-order correlated velocities that are residual fluctuations from circular motion. As written, however, there appears at first to be a gross mismatch: the left side of the equation is smaller than the right by a factor of order the ratio of the drift velocity to the rotation velocity. But this assumes that the length scales associated with the gradients on either side of the equation are comparable. Because the extracted free energy is in fact locally dissipated, and dissipation is dominated by the smallest scales, the gradient length scale on the left side provides the balance from the input of the right side by being, in effect, tiny. The analysis of nonrelativistic discs presented in Balbus & Hawley (1998) shows that when explicit dissipation terms are included in the energy fluctuation equation from the start, the balance struck is between the right side of equation (25) and explicit viscous (or resistive) dissipation. These energy loss terms are unimportant for large scale transport (and can be ignored for this purpose), but they represent the “thermal processor” between the extracted large-scale mechanical free energy and the disc’s radiative energy losses. As we have already noted, the thermal timescale over which this occurs is assumed to be rapid compared with the evolutionary time scale of the disc. This implies that the height-integrated, volume specific source term on the right side of (25) satisfies

(27) |

i.e. the total energy extracted over the local disk thickness is equal to the energy radiated through the upper and lower surfaces.

### 2.3 Large scale evolution

Henceforth, we drop the bars on the 4-velocities, and take these non- quantities to be understood as time-averaged means. If we integrate over height, assume axisymmetry, and ignore the pressure contributions, the equation of angular momentum conservation (10) now expands to:

(28) |

In the first term of (28), we have assumed that and are prescribed functions of only. The final term is obtained by integrating over height, which is now the angular momentum radiated from each side of the disc (Page & Thorne 1974). Using equation (2) for and simplifying, we obtain:

(29) |

where . Using now (29) back in equation (2), we find:

(30) |

The final step is to use equation (13) for . With

(31) |

this brings us to our governing equation:

(32) |

This is the equation we have been seeking. The first term on the right in square brackets is a straight translation from the Newtonian equation, while the second term is a relativistic correction stemming from the photon angular momentum.

A final point. We have been assuming that is a specified function of . If, however, has functional dependence upon , then would be implicitly time-dependent. In that case, equation (32) should be modified to:

(33) |

a form that holds more generally.

## 3 Solution of the evolutionary equation

### 3.1 Preliminaries

Let us introduce a more compact formulation. Define by

(34) |

Equation (32) becomes

(35) |

Next, with

(36) |

our governing equation takes the form of a pure diffusion equation

(37) |

This has a steady-state solution of . Equations (29) and (27) together imply

(38) |

where is the time-steady accretion rate and contains an additive constant boundary condition embodying the vanishing stess location (conventionally the ISCO radius). This is, in fact, the Novikov & Thorne (1973) solution in its entirety! This reader may wish to verify this for the relatively simple Schwarzschild limit with and

### 3.2 Modal solution

#### 3.2.1 Global WKB

We seek time-dependent solutions of the form . Then equation (37) becomes

(39) |

When is small (a not unphysical choice) or sufficiently large, equation (39) has the formal (unnormalised) WKB solution (Bender & Orszag 1978):

(40) |

When we should of course interpret this in terms of trigonometric functions. Returning to in preference to , and in preference to , we obtain

(41) |

Consider first an unstable mode, . At the ISCO location , our WKB solution formally breaks down, but as we shall see, it is still valid rather close to it. Let , so that positive and negative define regions of stable and unstable circular orbits, with and respectively. On physical grounds we certainly should not expect much of a disc-like structure to prevail for , but it is of interest to see how the equation discovers this on its own.

For , we have and the solution that is well-behaved as takes the form

(42) |

(We have chosen for later convenience a lower limit of integration to be . For a convergent integral, this simply amounts to setting the normalisation factor.) In the WKB limit, this is a sharply cut-off function in the bulk of the disc . For , and we may write down a formal solution:

(43) |

Here, the amplitude and phase are determined by the requirement that the solution join smoothly onto the exponentially cut-off solution for . This is already enough to see that unstable modes have significant amplitudes only inside of the ISCO, a physically very sensible result.

For stable () solutions, it is clear from our general solution (41) that the exponential cut-off behaviour now occurs inside the ISCO, , while for , the bulk of the disk hosts a spectrum of spatially-oscillatory, temporally-decaying modes.

#### 3.2.2 Local ISCO structure for nanvanishing

The ISCO is an apparent singularlity of our equation, which can, however, be treated rigorously. In the local vicinity of the ISCO equation (35) may be written

(44) |

where , , and the notation means that , , and are all evaluated at the ISCO . Note that the term is subdominant and has actually disappeared from the local ISCO-centred equation. (It will reappear as part of a locally determined normalisation factor.) We shall first assume that does not vanish, with the ultimate intent of showing the opposite: on physical grounds, it must vanish in a thin disc. With finite , the (unnormalised) solution to this equation for is

(45) |

where is the derivative of the Airy function. As in our WKB solution (42), positive values of the argument correspond to exponentially cut-off behaviour (the solution not chosen, Bi, rises exponentially), whereas negative values correspond to oscillatory behaviour. (See figure [1].) The “dispersion relation” we have found, the -definition of equation (45), may be written

(46) |

and exhibits violent instabilites on the smallest scales. This is a compelling reason to seek physically viable solutions with the ISCO boundary condition .

Before we do, however, we note a point of some mathematical consequence. The WKB solution (40) depends upon large for its validity, whereas the local solution merely requires . These are not mutually exclusive restrictions. There is no reason why they both cannot be valid in an overlapping domain. In this shared asymptotic regime, the two solutions must take one and the same form. To verify this is indeed so, note that the large argument expansion of the Ai function is (up to an overall normalisation):

(47) |

which is exactly the same form as equation (42) in the limit , (once again up to an overall normalisation):

(48) |

Equations (47) and (48) have exactly the same functional form, as was sought.

#### 3.2.3 A uniformly valid solution

The agreement between the two solutions in an overlapping asymptotic zone suggests the possibility that there may be a single analytic formula that is valid everywhere. Such a solution is known to exist for a certain class of “one-turning-point problems,” in quantum mechanical solutions of the Schrödinger equation (Bender & Orszag 1978). Rather than derive this function, it is simplest just to write it down, and then verify that it reduces to each of our asymptotic forms in the appropriate limits.

Define by

(49) |

Then, our (unnormalised) uniformly valid solution is

(50) |

To verify this, we assume first that . In the limit , Ai has the asymptotic form

(51) |

so the factors cancel in (50), and we are led directly to equation (42) for . Next, when and , we expand , and becomes

and then reduces to equation (45): , since is locally equal to the constant . Finally, when away from the ISCO, then . Multiply by unity, written as . Then,

Switching the limits and , this is the same as

i.e, is a purely real negative quantity, despite all of the complex-valued exponents and nested fractional powers. Then, use the standard large negative argument for Ai (Bender & Orszag 1978):

(52) |

It is easy to see that with equation (50), by adjusting and the above leads to a precise match with equation (43). Equation (50) is therefore a uniformly valid solution to equation (32). Figure (1) shows an explicit solution for Schwarzschild geometry with , chosen for ease of analytics. (Recall that correlates angular momentum and radial velocity fluctuations, so there is some sense for it to increase slowly with .) In this case, the needed integral over is (see the end of §3.1):

where is .

### 3.3 Modal solutions for vanishing

#### 3.3.1 Exterior region,

If there is any finite stress at the the ISCO, then there are extremely unstable modes present on small scales. This is a compelling argument in favour of the customary boundary condition of setting for . Let us see how this removes the unstable behaviour.

We shall assume that vanishes. In equilibrium, . From the definitions of and in (36), it follows that . The question then is what are the solutions of (32) near the ISCO with this behaviour for ?

Set the local stress . Then, the local ISCO equation is

(54) |

or

(55) |

By equation (31), itself must now vanish at . If , then there are two formal solutions to this equation in the region :

(56) |

But the solution is not well-behaved for large , and the solution does not vanish at , so in fact there are no solutions compatible with the boundary conditions. In other words, there are no unstable solutions for .

Consider next . Then, the solution satisfying the vanishing boundary condition at is

(57) |

where is the Bessel function of order . The corresponding solution with the Bessel function does not vanish at . Hence, there is a well-determined stable set of modal solutions with vanishing at the ISCO.

Once again, there is an overlap zone near the ISCO in which the WKB solution is valid together with the small local solution. The large argument expansion of (57) is

(58) |

The WKB solution follows from (43), now outside the ISCO:

(59) |

where and are once again arbitrary. In the limit , we require and . The integral in (59) is then precisely . With the proper choice of and , there is a complete agreement of functional form between (59) and (58).

Finally, there is once again a simple, uniformly valid solution. With now defined by

(60) |

the desired solution is

(61) |

To verify this, simply expand the above: first for large (recovering [59]), then for small (recovering [57]), and then simultaneously for large and small (recovering [58]). It is readily seen that the function (61) reduces to all proper asymptotic forms.

#### 3.3.2 Interior region,

For , there are no stable solutions that are well-behaved. The unstable , but spatially well-behaved, interior solution is now easy to construct, since it has precisely the same mathematical form as the exterior solution. Moreover, vanishing at together with its first derivative, this solution seems to join smoothly onto the stable exterior solution. The smoothness is maintained even though the growth rate is different on either side of ! How does it make sense to have a global “mode” with two different growth rates, one with positive , the other with negative , in different regions?

Of course a single mode cannot have different growth rates in different disk regions. What we have been discussing is in reality a superposition of two modes. This points to the resolution of our problem. The location is a branching singularity of the governing equation, and lacks a unique prescription for traversing it. It is an “improper node.” All modal solutions have vanishing and at . In particular, a smooth solution to our problem is that equation (61) holds for , and for , a stable mode that lives entirely in the exterior bulk of the disk. Similar considerations hold for its “dual,” unstable solution, in this case with vanishing for . Thus the answer to the question posed at the end of the previous paragraph is that the peculiar global solution described is not, in fact, one mode, but a superposition of two. What is perhaps unusual is that each mode vanishes identically where the other is present! The stable exterior mode is unique, and for astrophysical purposes, the mode of interest.

### 3.4 Green’s function solution

With given by equation (61), we may construct more general solutions by superposing modes. Formally, we may write

(62) |

where is whatever appropriate function we choose. (For convenience, we have replaced by in this role as a dummy variable, and the explicit dependence is exhibited in .) Consider next the integral (Gradshteyn & Ryzhik 2014):

(63) |

where is the Bessel function of order and the corresponding modified Bessel function. In the limit , this integral represents a delta function type of concentration at which then spreads as increases. This behaviour, together with symmetry, is what we seek for a Green’s function response, initially concentrated at . We are assuming that our global WKB solution holds over all in the integral, an assumption that must break down at . But in the limit of small , this will affect only the detailed behaviour at very late times; the small contribution to the integral is otherwise negligible.

Combining the results of equations (61), (62), and (63) allows us to write down the (unnormalised) Green’s function solution to our equation. With and

(64) |

At early times , the asymptotic behaviour of the terms on the final line of (64) simplifes to:

(65) |

This takes on the classic form for the diffusion of an initially very localised concentration. The local-frame surface emissivity is then given directly by (13):

(66) |

## 4 Discussion

In this paper, we have derived a form of the thin disc diffusion equation that is valid for general relativisitc spacetimes. We assume only that the metric tensor is axisymmetric, so that our equation is suitable for both the Schwarzschild and Kerr geometries. Remarkably, the metric itself enters into the calculation only in the form of a determinant, which is then absorbed as a multiplicative factor of our surface density variable. It then disappears entirely from the calculation.

The physics of an evolving thin relativistic disc compels unstable modes trapped inside to rapidly destroy their host in this “Rayleigh-unstable” zone. Quasi-stable equilibrium circular orbits are mathematical fantasies here: without a retaining potential, the orbits simply plunge. The exterior modal solutions, by contrast, are always stable. On the other hand, by imposing the boundary condition of a vansihing stress tensor at (), stable modes exist exclusively in the stable region, vanishing together with their first derivatives as from positive values. Modes with finite at must penetrate, at least exponentially, into the plunging region, and in thin disc models would not be supported.

By using a combination of WKB techniques, local analysis and matched asymptotic expansions, it is possible to solve very generally the disc diffusion equation in terms of quadratures. This is the scope of the current paper. Using these methods, we have been able to calculate the Green’s function solution, which is expected to be valid up until very late times when a quasi-equilibrium is reached. These findings may be useful as numerical diagnostics, but only if the thin disc condition is well-satisfied, a limit that has yet to be convincingly simulated. The most interesting astrophysical application of this work is likely to be to black hole transients. These include dramatic state changes in which the inner regions of the disc are thought to disappear and then reform, and tidal disruption events in which a disc forms and subesquently accretes from the debris of a mangled star. In principle, these events may be modelled by the one-dimensional evolution equation (32) with appropriate boundary conditions, and the time-dependent surface emission calculated in the observer’s reference frame. These interesting possibilities will be the subject of future investigations.

## Acknowledgements

It is a pleasure to acknowledge important discussions with W. Potter during the early formative stages of this work, and very constructive comments from C. McKee, M. Rees and C. Gammie on an earlier manuscipt. I am grateful for support from the Royal Society in the form of a Wolfson Research Merit Award, and from STFC (grant number ST/N000919/1).

## References

- [1] Balbus, S. A., & Hawley, J. F. 1998, RMP, 70, 1
- [2] Balbus, S. A. & Papaloizou, J. C. B. 1999, ApJ, 521, 650
- [3] Bender, C., & Orszag, S. 1978, Advanced Mathematical Methods for Scientists and Engineers, (New York: McGraw-Hill)
- [4] Gradshteyn, I. S., & Ryzhik, M. 2014, Table of Integrals, Series, and Products, (New York: Academic Press) [eq. 6.633 (2.)]
- [5] Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
- [6] Novikov, I. D., & Thorne, K. S. 1973, Black Holes—Les Astres Occlus, ed. C. De Witt, (New York: Gordon and Breach), p. 346
- [7] Page, D. N., & Thorne, K. S. 1974, ApJ, 191, 499
- [8] Pringle, J. E. 1981, ARAA, 19, 137