A full relativistic thin disc -- the physics of the plunging region and the value of the stress at the ISCO

The widely used Novikov-Thorne relativistic thin disc equations are only valid down to the radius of the innermost stable circular orbit (ISCO). This leads to an undetermined boundary condition at the ISCO, known as the inner stress of the disc, which sets the luminosity of the disc at the ISCO and introduces considerable ambiguity in accurately determining the mass, spin and accretion rate of black holes from observed spectra. We resolve this ambiguity by self-consistently extending the relativistic disc solution through the ISCO to the black hole horizon by calculating the inspiral of an average disc particle subject to turbulent disc forces, using a new particle-in-disc technique. Traditionally it has been assumed that the stress at the ISCO is zero, with material plunging approximately radially into the black hole at close to the speed of light. We demonstrate that in fact the inspiral is less severe, with several (~4-17) orbits completed before the horizon. This leads to a small non-zero stress and luminosity at and inside the ISCO, with a local surface temperature at the ISCO between ~0.15 and 0.3 times the maximum surface temperature of the disc, in the case where no dynamically important net magnetic field is present. For a range of disc parameters we calculate the value of the inner stress/surface temperature, which is required when fitting relativistic thin disc models to observations. We resolve a problem in relativistic slim disc models in which turbulent heating becomes inaccurate and falls to zero inside the plunging region.


INTRODUCTION
Black hole accretion discs offer an exceptional insight into the relativistic physics governing material being dragged into black holes in strong gravity. Accretion discs are incredibly luminous with emitted spectra that are sensitive to the mass, spin and accretion rate on to the black hole McClintock 2006 andDone et al. 2007). This makes accretion disc observations particularly useful in determining the properties of the population of black holes in our Universe and the impact they have on their environments. However, this makes it all the more important for our disc models to be as accurate as possible.
The concept and equations governing accretion discs were first outlined in a series of seminal papers by Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974), and extended to relativistic discs by Novikov & Thorne (1973) and Page & Thorne (1974). These disc models still form the basis of the majority of spectral models used today. One shortcoming with Novikov-Thorne (NT) relativistic disc models is that they implicitly assume that at each radius the orbiting material is on an approximately circular orbit and so the equations become invalid inside the radius of the innermost stable circular orbit (ISCO), where circular orbits cease to exist. This E-mail:willpotter37@gmail.com leaves two problems: the equations are not valid in the plunging region between the ISCO and event horizon and so this section of the disc is not modelled and is assumed to have no contribution to luminosity; we are left with an undetermined inner boundary condition at the ISCO that leaves the stress (or torque), luminosity and other disc properties undetermined at the ISCO. This ambiguity leads to a variety of different choices of inner boundary conditions in widely used spectral models from: the most commonly assumed zero-stress inner boundary condition, e.g. ezdiskbb , leaving the stress at the ISCO as a user-defined parameter, e.g. KERRBB (Li et al. 2005) and BHSPEC (Davis et al. 2005), or choosing a relatively large inner stress, e.g. diskbb (an XSPEC model, Arnaud 1996, based upon Mitsuda et al. 1984and Makishima et al. 1986). This uncertainty is particularly problematic because the region around the ISCO dominates the disc luminosity and spectrum. Thus the inferred black hole properties and calculated spectrum are very sensitive to the particular choice of inner boundary condition .
In Novikov & Thorne (1973), it was assumed that the stress at the ISCO would vanish because material inside the ISCO would plunge rapidly, approximately radially into the black hole, at close to the speed of light. This would mean the material inside the ISCO would be very diffuse, with minimal shear, so would not provide any significant stress/torque on material at the ISCO. The case for a negligible inner stress in thin discs with H/r 1 has also been argued by Paczyński (2000). The work by Agol & Krolik (2000) showed explicitly what the algebraic form of the undetermined inner stress would be if it is chosen to be non-zero for a Novikov-Thorne solution, motivated by the importance of strong net magnetic fields that could be advected through the disc and enhanced by flux freezing and differential rotation in the plunging region (Gammie 1999). They pointed out that whilst the disc gas might not exert a significant torque on the inner disc, strong magnetic fields may do so.
The advent of general relativistic magnetohydrodynamic (GRMHD) simulations (Gammie et al. 2003 andDe Villiers &Hawley 2003), has allowed huge progress in our understanding of the plunging region. They have shown that large net vertical magnetic fields (the direction perpendicular to disc plane, or zdirection) can indeed extract rotational energy from a spinning black hole via the Blandford-Znajek effect (Blandford & Znajek 1977) and due to flux freezing poloidal fields can become dynamically important in the plunging region producing jets and disc stress (e.g. McKinney 2006, Hawley & Krolik 2006, McKinney & Blandford 2009. In extreme cases, merely the advection and compression of large net poloidal fields can create large enough magnetic pressure gradients to stall the disc in magnetically arrested discs, e.g. Narayan et al. (2003), Tchekhovskoy et al. (2011) and McKinney et al. (2012). Unfortunately, however, such simulations do not offer a definitive resolution to the inner stress problem because the undetermined inner stress is traded for an undetermined initial magnetic field in the gas. Another important problem for 3D magnetohydrodynamic (MHD) simulations of accretion discs in general is that it is too computationally intensive to run simulations of truly thin discs (H/R 0.1, the cutting-edge simulation by Liska et al. 2019 had H/R ∼ 0.03) or to run simulations for long enough to reach a true equilibrium solution. This is because thinner discs require finer spatial resolution to resolve their vertical structure and have longer viscous evolution time-scales than thick discs, and additionally the time taken to reach a steady state increases rapidly with radius outwards through the disc. Also, the strength of MHD turbulence in the disc, usually parameterised by the famous α-parameter (Shakura & Sunyaev 1973), does not numerically converge with increasing resolution in general Papaloizou 2007 andRyan et al. 2017) and depends upon the initial net magnetic field (Hawley et al. 1995), resistive and viscous dissipation coefficients (Fromang et al. 2007 andSimon &, and vertical size of the simulation (Ross et al. 2016 andShi et al. 2016). This means that currently GRMHD simulations are not able to uniquely determine the properties of steady-state relativistic thin disc solutions. They are also too computationally expensive and complex to run to be practical for them to be the basis of widely used spectral fitting models.
Resolving the question of the stress present at the ISCO and inside the plunging region is an important question and has been the subject of many previous studies. Slim disc models include additional physics relevant when terms of order (H/R) 2 are not neglected and discs are radiatively inefficient (Paczynski & Bisnovatyi-Kogan 1981, Muchotrzeb & Paczynski 1982and Abramowicz et al. 1988). These additional effects include radial pressure gradients, radial momentum conservation and heat advection, and crucially, allow for angular velocities that differ from circular orbits. Slim disc solutions are calculated by conserving the mass, momentum and energy flux through the disc on a spatial grid and searching for solutions that satisfy these conservation laws and smoothly pass through sonic points, matching on to an outer bound-ary. A general result from slim disc studies is that the stress at the ISCO decreases as the disc becomes thinner (see e.g. , in agreement with Paczyński (2000).
General relativistic slim disc equations were introduced by Lasota (1994), Abramowicz et al. (1996), Peitz & Appl (1997), Beloborodov et al. (1997), Igumenshchev et al. (1998) and , which can be solved using a method of relaxation on a numerical grid (e.g. Narayan 1991 and. Relativistic slim disc models allow solutions to extend smoothly through the ISCO down to the black hole horizon thereby allowing the stress at the ISCO and physics of the plunging region to be calculated. However, because slim disc models were primarily introduced to investigate radiatively inefficient thick discs, it was reasonably assumed that in this regime physical effects relevant to radiatively efficient discs could be neglected or approximated (Abramowicz et al. 1996). In particular, the angular momentum transported away by radiation for an optically thick disc is generally neglected and various approximations to the full velocity shear tensor are used that break down in the plunging region, resulting in inaccurate turbulent heating rates. Whilst these effects may not be significant in radiatively inefficient discs, they become important when considering radiatively efficient thin discs. Investigations using slim disc models, understandably, do not focus on a detailed resolution of the value of inner stress for thin discs, usually stating that such stress is negligible, e.g. Abramowicz et al. (2010)), consistent with the arguments of Paczyński (2000). This means that a comprehensive investigation of the inner stress and plunging region of relativistic thin discs is still required to offer an accurate quantitative answer to the question of how large the ISCO stress is for thin discs and to what extent Novikov-Thorne solutions are accurate.
Aside from slim disc models other attempts to resolve the problem include extending the Novikov-Thorne solution through the ISCO by assuming the gas is in free-fall in the plunging region and thereby calculating the shear and heating rate in the plunging zone, e.g. Penna et al. (2012). This assumption is sensible when the stress and heating rate are negligible in the plunging region, however, any significant stress would lead to accelerations violating the free-fall assumption and so the method does not in itself allow a self-consistent calculation of the stress at the ISCO. Extensive progress has been made in 3D GRMHD simulations of thin discs that generally find increased values of stress and heating at and inside the ISCO and differences from the Novikov-Thorne solutions, due primarily to additional magnetic stresses (e.g. Beckwith et al. 2008b, Noble et al. 2009, Noble et al. 2010and Penna et al. 2010, Noble et al. 2011). The problem is that these stresses depend sensitively on the details of the simulation, such as initial magnetic field strength and geometry, and numerical methods (Penna et al. 2010). This leads to some disagreement about what is responsible for the significant differences between simulations with similar initial conditions. Some simulations have found that as thinner discs are simulated, the luminosity and stress in the plunging zone seem to decrease, with differences to the Novikov-Thorne solution decreasing Penna et al. 2010), whilst other simulations find no systematic dependence on disc thickness (Noble et al. 2010). It was suggested by Penna et al. (2010) that the differences between these simulations could be attributed to different initial disc magnetic field and definitions of disc scale height over which averages were taken, whereas Noble et al. (2011) suggested differences might be due to different numerical grid resolutions, amongst other possibilities. It therefore remains unclear from simulations what the value of the stress and luminosity in the plunging zone for a realistically thin steady-state disc should be.
The purpose of this paper is to resolve the following important questions for steady-state relativistic thin discs: What is the value of stress at the ISCO and in the plunging region, and how does it depend on disc properties? What are the physical properties of the gas inside the ISCO? How rapid is the plunge from the ISCO into the black hole? To answer these questions we develop a new particle-indisc technique that calculates the full general relativistic equation of motion of a gas particle experiencing turbulent stresses as it spirals into the black hole whilst conserving mass, momentum and energy in the disc. This produces stable steady-state solutions, extending self-consistently through the ISCO to the black hole horizon. This eliminates the undetermined boundary condition of the inner stress at the ISCO and allows the physics of the plunging region and value of stress at the ISCO to be calculated.

RELATIVISTIC DISC EQUATIONS
The conservation equations for relativistic thin accretion discs were first derived by Novikov & Thorne (1973) and Page & Thorne (1974), where the disc structure was solved by using the conservation equations for mass, angular momentum and energy. The equations for conservation of rest mass (or more accurately particle number) and relativistic energy-momentum are where the rest mass density is ρ = nm, n is the number density of particles measured in the fluid rest frame and m is the average rest mass of the particles (with m ∼ 0.615mp for solar chemical abundances, Frank et al. 2002, and mp the proton mass), U µ the fluid four-velocity and T µν the stress-energy tensor. The spacetime metric for a rotating black hole in a vacuum is given by the Kerr metric. We make the standard assumption that the mass in the accretion disc is not sufficient to perturb the metric significantly. The Kerr metric is given by +ρ 2 dθ 2 + r 2 + a 2 r 2 g + 2r 3 g ra 2 sin 2 θ ρ 2 sin 2 θdφ 2 , where M is the mass of the black hole, c the speed of light in vacuum, rg = GM/c 2 is the gravitational radius, and the dimensionless spin a is related to the angular momentum J of the black hole via J = M argc, with −1 < a < 1. Following Novikov and Thorne we make the assumption that the rotation axis of the accretion disc is parallel to the rotation axis of the black hole, i.e. θ ≈ π/2. We also make the assumption that our disc is relatively thin so that its height is much smaller than the radius from the black hole (H/r 1), which allows us to simplify the metric by introducing rdθ ≈ dz. The Kerr metric thereby simplifies to Hereafter we set c = 1, to simplify our equations.

Stress-energy tensor
The total stress-energy tensor T µν of the disc can be broken down into a linear sum of contributions from: the thermal part of the gas motion T µν g , the turbulent component T µν t , heat transfer via radiation T µν h , and the electromagnetic stress-energy T µν EM (Novikov & Thorne 1973). For the sake of brevity a more detailed discussion and derivation of these components is given in Appendix A.
where e is the internal energy density of the plasma, p is the isotropic pressure, q ν = (0, 0, 0, q z ) is the radiative heat flux (Eckart 1940), ν is a measure of the average correlation in the turbulent velocity, and magnetic field induced by the velocity shear of the fluid σ µν , where ν = αcsH can also be interpreted as an effective turbulent kinematic viscosity, with the gas sound speed c 2 s = p/ρ and disc scale height H. The internal energy and pressure associated with turbulence and radiation are small compared to the thermal gas contribution for the disc parameters studied in this paper because the turbulent component is of order α smaller than the gas pressure, and for our range of accretion rates (0.001 − 0.1Ṁ Edd ) the radiation pressure and energy density are also smaller than those of the thermal gas (except in Solution 5, with a = 0.99, where the radiation pressure slightly exceeds the gas pressure around the temperature maximum). Thus for an ionised non-relativistic gas we have the standard result e = 3p/2, p = nkBTc, with kB the Boltzmann constant and Tc the central temperature of the disc.

Mass conservation
The equation for rest mass conservation follows from the conservation of particle number (1) ∇µ(ρU µ ) = 0.
The equation can be simplified using the identity (Misner et al. 1973) valid for an arbitrary four-vector, X µ , and where g is the determinant of the metric ( √ −g = r in our case, eq. 3), Integrating over four-dimensional volume, dtdrdφdz, of a thin annulus of width dr and using Gauss's theorem in four dimensions, the equation for the conservation of mass in the disc becomes Evaluating the first term explicitly, where we have used the definition of the partial derivative and integrated the density and velocity over the disc scale height. The surface density of the disc is defined by Σ = ∞ 0 ρdz ≈ H 0 ρdz ≈ ρH, and all quantities will be vertically averaged, except where explicitly stated, defined by (11) For simplicity of notation we drop the angle brackets from vertically averaged quantities in the rest of the paper because it will be made clear from the context and derivation whether particular equations and quantities are vertically averaged or not. The third term is zero for an axisymmetric disc and the second and fourth terms in (8) can be integrated and evaluated similarly to the first term to give (13) The fourth term above, representing mass loss due to a vertical wind in the disc, has a slightly different form since we are directly interested in the mass flux leaving the disc at its upper and lower surfaces. We require the bounds on the z integration to be constant and sufficiently large to enclose the entire vertical extent of the disc with substantial density at all radii, however, since we know that in practice the disc density falls off exponentially with scale height H, the density falls off rapidly for z > H(r). This means that the bounds on our integral z = ±∞ can be sensibly replaced by z ≈ ±H(r). In terms of vertical fluxes leaving the disc (e.g. winds and radiation), we are interested in the vertical flux leaving the disc surface at z ≈ ±H that reaches infinity, or does not return to the disc at another radius. It is sensible to also replace the bounds of these integrals to z ≈ ±H(r) to represent the vertical flux leaving the disc surface instead of the flux reaching z = ±∞ at that radius, which would lead to unnecessary confusion about travel times and the original radius of the disc the flux was emitted from, since realistic winds and emitted fluxes will not travel purely vertically. Substituting the three terms back into equation (8) we find our equation for mass conservation We can simplify the equation by assuming that the vertical mass flux, ρU z , through the upper and lower surfaces of the disc is symmetric about the disc plane, i.e. equal and opposite (ρ In the case where we look for a steady-state solution with no mass loss from a disc wind this simplifies to the standard mass conservation equation (Novikov & Thorne 1973) − 4πrΣU r =Ṁ , where √ −g = r for the metric (3) and the constant,Ṁ , is the steady-state mass accretion rate.

Geodesic equation and four-acceleration
The general relativistic equation of motion for a particle with fourvelocity dX µ /dτ = U µ in the presence of a four-force exerting a four-acceleration A µ is given by where Γ a bc are the Christoffel symbols of the metric, D/Dτ is the absolute derivative following the particle and d/dτ = U µ ∂µ. In the absence of an explicit four-force and four-acceleration, the equation reduces to the geodesic equation. In the rest frame of a gas particle the four-acceleration is always perpendicular to the velocity (i.e. the rest mass of the particle remains unchanged by the acceleration) so that A µ Uµ = 0.

Relativistic fluid equations
Our goal is to calculate the acceleration acting on an average gas particle in the disc as it accretes and spirals into the black hole through the disc. Starting from the equation of energy-momentum conservation (1), Rearranging the first term on the right, where we have used conservation of mass (1) to simplify. Substituting this into (18) and rearranging to find the fluid four-acceleration using A ν = U µ ∇µU ν (17) This is the relativistic version of the Euler equation in which, analogously to the Newtonian equations, forces appear as the divergence of stress tensors. The terms on the right-hand side of the equation are the accelerations due to: rate of change of momentum of the internal energy and pressure, pressure gradient force, forces due to turbulent shear stresses, and forces from the emission of radiation. Before proceeding further and evaluating the turbulent and radiative stress terms, it is necessary to define the standard decomposition of the covariant derivative of the velocity field into the rotational, compressional, shearing and accelerating components (Eckart 1940 andMisner et al. 1973).

Decomposition of velocity derivative
The Newtonian origin of the shear tensor is found by decomposing the strain-rate tensor ∂jvi, into (1) -an antisymmetric tensor ωij representing the rotation of the velocity field; (2) -θ the scalar trace of the tensor representing uniform expansion or compression; and (3) -σij the traceless symmetric part of the tensor representing the pure shear rate. The generalisations to the fully relativistic expressions are (Misner et al. 1973) where P αβ is the projection tensor that projects a vector on to a 3D surface orthogonal to U µ (i.e. 3D spatial coordinates in a local rest frame, Misner et al. 1973). The relativistic expression also contains an additional term AαU β due to the four-acceleration of the fluid (17), which encapsulates information contained in the time derivative of the velocity field and which does not appear in the 3D spatial Newtonian strain-rate tensor. Using eq. 25 for the velocity shear tensor we can calculate the turbulent viscous prescription of the internal stress tensor in eq. 4. It is important to mention that a Newtonian viscosity, in which the stress-energy is directly proportional to the velocity shear tensor, is problematic in relativistic fluids since Newtonian viscosity, described by a standard diffusion equation, allows instantaneous information communication across the fluid (i.e. the wave speed associated with the viscosity exceeds the speed of light). The correct relativistic generalisation of heat flow and viscosity is an active area of research, with second-order, stable, causal theories being difficult to implement in practice due to the requisite detailed knowledge of the fluid behaviour and their inherent complexity, Israel & Stewart (1979) and Carter (1989) (see e.g. Lopez-Monsalvo 2011 for a general overview, and Szuszkiewicz 1994 and for a discussion of causal viscosity in accretion discs). Since in this work we are interested in a steady state equilibrium, it is assumed that sufficient time has passed for information of the turbulent stresses to have spread through the disc many times over, such that the disc has had time to settle into a steady-state. In this case the issue of the speed of propagation of information is not important to our results and we assume that in the local rest frame of the fluid the correct general expression for viscosity reduces approximately to the classical Newtonian viscosity in our steady-state limit. Any change in viscous stresses resulting from finite propagation speed effects will also be compounded by turbulent time delays (the time taken for turbulent stresses and dissipation to adjust to changes in the velocity shear and changes in fluid properties). This means that a more complex and realistic model of turbulence than the α model, where such adjustments are instantaneous, would be required to properly deal with these effects. This is beyond the scope of this work, but may be worth pursuing in the future. As explained in Section 4, by deliberate construction, our numerical method integrates the equations inwards, from large to small radii, propagating information inwards only (i.e. derivatives are calculated only using values at the current radius and larger radii). This choice eliminates the problems with causality associated with information propagation outwards through the sonic point or event horizon from using a diffusive viscosity with formally infinite wave speed. However, because our solutions are steady state and the inner boundary of the solution is contained within the event horizon (and because we do not consider the accumulation of net electric or magnetic fields on to the black hole), our solutions only depend on information propagating inwards from the outer boundary and so our numerical choice to eliminate outwards propagating information does not affect the accuracy or validity of our steady-state solutions.

Vertical velocity and compression
There is some subtlety in choosing how to deal with the vertical velocity and vertical compression in the disc, since the detailed vertical structure is assumed to always be close to hydrostatic equilibrium and the detailed dynamics are not solved. Nevertheless, it is clear that if the disc height H changes this must be due to a redistribution of vertical density resulting from a non-zero vertical velocity. The simplest way in which to define the average vertical velocity of the disc is to calculate the vertical velocity required for the disc height to change correctly as plasma flows through the disc, i.e.
This is the vertical velocity of the disc surface required to keep the disc surface at z = H, as H changes with radius. This allows us to calculate the vertical component of the average disc compression θz This is the average volumetric vertical compression rate for the disc, useful in determining work done by pressure forces for example. This allows the total compression rate of the disc to be calculated by equation 24.
where we have used equation 6. The simplest approximation for the density averaged vertical velocity, consistent with the thin disc assumptions, is to assume that in the disc mid-plane the compression is uniform. This allows the additional components of the velocity shear tensor and turbulent stress tensor to be estimated via Since the vertical disc structure is not solved in detail here, the primary purpose of this approximate expression is to estimate the relative importance of changes to the disc height on disc heating and momentum transport. Using (29) we can estimate the zz, zr and rz components of the velocity shear tensor. It is found that the vertical compression has a substantial contribution to disc heating once inside the ISCO (see Figs. 7 and 5b), becoming the dominant heating source close to the event horizon. This is because the disc rapidly decreases in height due to increased gravitational forces close to the black hole, with the gravitational force doing work against pressure by vertically compressing the disc. However, including the vertical velocity and vertical compression in our calculation of the turbulent stresses via the velocity shear-tensor (25) has a negligible effect on the four-acceleration around and outside the ISCO, making a small contribution to the total four-acceleration, A φ , in the vicinity of the event horizon (∼ 4%). This is important because it means that a more detailed treatment of the vertical structure and vertical velocity is not required to calculate accurate solutions (only the rate of change in total disc height H is required to calculate disc heating), and the usual thin disc assumptions remain valid outside of the ISCO. Based on this information we only include the effect of vertical compression when calculating the compressional disc heating (47) and the four-momentum carried by this heat (46), where the effect is important. The vertical velocity and vertical compression are neglected when calculating the turbulent stresses (30) and turbulent heating rate where its effects are much smaller.
It is worth briefly commenting that the density averaged vertical velocity (26) does not tell us about any mass flux launched from the disc surface in the form of a wind, since the average vertical velocity simply tells us the rate of change of the scale height of the disc and the average velocity of the majority of mass. It is possible to have a disc that has a reducing scale height and is on average undergoing compression, but which is still launching mass in a wind from the disc surface, i.e. Fm ≈ (ρU z )(z ≈ H) = (ρ U z )(z = H). This means that mass flow in the form of a wind and the vertical velocity of the wind require separate treatment.
The consideration of vertical velocity and compression are usually neglected in previous work on relativistic thin and slim discs (e.g. Abramowicz et al. 1996). This is a reasonable assumption outside of the plunging region because these effects are generally small outside of the ISCO (see Fig. 7). However, as we shall later demonstrate, this situation does not remain valid in the plunging region where both the radial and vertical disc velocities increase rapidly, leading to substantial cooling via radial expansion and heating via vertical compression (see Figs. 5b and 7).

Turbulent and radiative stresses
Evaluating the divergence of the turbulent stress is relatively straightforward The divergence of the radiative heat stress-energy is In the fluid rest frame the average radiative energy flux is assumed to be travelling in the vertical direction. This means that performing local Lorentz transformations in the r and φ directions to convert from the rest frame to the lab frame (where the fluid possesses only average radial and azimuthal velocities) does not change the direction or magnitude of q µ so that q µ = q µ = (0, 0, 0, q z ). Simplifying equation 31, The first, third and fourth terms are only relevant when considering the vertical structure and vertical acceleration of the disc. Only the second term is significant for our vertically-averaged disc model and is retained. Substituting equations 30 and 32 into 21, where we have used the definition of a covariant derivative on a scalar function, ∇µp = ∂µp. Finally, the equation needs to be vertically integrated and averaged. The vertical averaging procedure is slightly complicated by the changing scale height of the disc with radius and so some care is needed when evaluating vertically integrated derivatives. A vertically integrated partial derivative of an arbitrary function f , where µ = z, will have a form which takes account of the different scale heights evaluated at different radii. For partial derivatives in z we find following from the definition of integration as the inverse of the derivative. The integral of the z-derivative therefore gives the value of the function at large distances from the disc (in the case of conserved fluxes this is equivalent to the flux leaving the disc surface that will reach large distances), with f = ρU z this is the mass flux leaving the surface as a wind, Fm, or with f = q z gives the radiative heat loss from the surface, F . This means that in the case of turbulent fields and velocities, which are confined within the disc, there is no transfer of momentum by the turbulence out through the disc surface in the z direction to large distances, i.e.
Vertically integrating eq. 33 using equations 34, 35 and 36, where i only runs over t, r and φ, we have used eq. 6 to simplify the first term in brackets on the right, and it is assumed that there is no mass or enthalpy flux through the disc surface or turbulent or magnetic stresses acting outwards through the disc surface so that the integral of the µ = z terms vanish with the exception of the radiative flux F = q z (z = H) = σT 4 S , where TS is the surface temperature of the disc and σ is the Stefan-Boltzmann constant. Because of the assumed symmetry of the disc about its mid-plane the vertically averaged A z must necessarily be zero and it is not necessary to calculate A t since we can use U µ Uµ = −1 in its place (although A t can be trivially calculated by swapping φ for t in eq. 37). It remains to calculate the energy flux term using the energy equation.

Energy equation
It is necessary to specify an explicit equation for the conservation of energy to find solutions that balance the local rate of turbulent energy dissipation (heating), radiative cooling, compressive heating and the rate of change of thermal energy in the fluid. Similarly to the non-relativistic fluid equations, the energy equation is found by contracting the equation of four-momentum conservation with the fluid four-velocity (see for example Novikov & Thorne 1973) Uν (∇µT µν ) = 0.
Rearranging using the product rule for derivatives, Contracting equation (4) with Uν we find where the pressure terms cancel, we have used q µ Uµ = 0 and, Uν σ µν = 0, which follows from the definition of the shear tensor (25). Taking the covariant derivative of the expression above we calculate using the Leibniz rule. From conservation of particle number (1) this can be simplified to Contracting the energy momentum tensor with the relativistic strain-rate tensor (22) and using σ µν Finally, substituting in equations 43 and 44 back into 40 and simplifying by swapping dummy indices and using the symmetry of σµν we find The right-hand side of the equation contains the terms responsible for generating heat that from left to right are heating due to turbulent shear stresses, compressive heating (work done against pressure) and a relativistic term given by the perceived inertia or change in the rest energy of the heat flux q µ when the rest frame is accelerating (responsible for the temperature decrease with height in an atmosphere subject to gravity, Tolman 1930 andTolman &Ehrenfest 1930). The left-hand side tells us what happens to this dissipated energy, from left to right these terms are the change in internal energy/temperature of the plasma and the heat transported out of the volume. Integrating and averaging in the vertical direction and using equations 4, 6 and 35 this becomes where we have assumed no mass flux in the form of a disc wind is present and used eq. 6. The term associated with the redshifting of the heat flux Azq z dz is neglected since the vertical acceleration is small in a thin disc (Novikov & Thorne 1973), the radiative flux is F = q z (z = H) = σT 4 S and the average compression θ is defined in eq. 28. This term is substituted into equations 37 and 38 to solve for the four-acceleration of an average disc particle. In order to calculate the change in temperature of the plasma due to heating and cooling, we can use conservation of mass (16) to simplify the vertically integrated term (1/r)∂r(rHeU r ) = Σ(dΠ/dτ ), by introducing the internal energy per unit rest-mass density, Π = e/ρ, In the standard thin disc model it is assumed that the internal energy of the gas is negligible and so all locally dissipated energy is immediately radiated away (e ∼ Π ∼ 0) and that heat generated by shear stresses dominates over compressive heating effects (Novikov & Thorne 1973). Whilst these assumptions remain valid away from the ISCO, it is not clear whether they will hold inside the plunging region where the radial velocity becomes large and the infall timescale becomes small. These conditions are likely to lead to the increased importance of heat advection, vertical compressive heating and radial cooling via expansion on the gas. For these reasons we keep all terms in equations 37, 38 and 47, including terms usually associated with radiatively inefficient thick/slim discs.

Vertical structure and opacity
In order to close the set of equations derived so far the vertical structure and opacity of the disc need to be specified. For optically thick, geometrically thin discs the standard procedure is to calculate the vertical radiative transfer assuming a plane-parallel atmosphere approximation with a frequency-averaged grey opacity. The standard result under these assumptions, in the rest frame of the gas is (Novikov & Thorne 1973) The disc opacity κ is usually approximated by the sum of the electron scattering opacity and Kramers opacity (Novikov & Thorne 1973), following NT, κ k is determined by the free-free opacity, which is expected to exceed the bound-free opacity for disc temperatures, Tc 8 × 10 5 K e.g. Novikov & Thorne (1973) and Ginzburg (1967). The analytic Novikov-Thorne solutions are calculated in different regimes in which either gas or radiation pressure is dominant and either electron scattering or Kramers opacity is dominant (Novikov & Thorne 1973). To facilitate a fair comparison between our solutions and the NT solutions, our solutions use only the opacity that is dominant in the region with maximum luminosity, just outside the ISCO, to compare to the relevant NT solution.

Disc scale height
The characteristic exponential scale height of the disc H is the approximate vertical distance over which the density at the mid-plane of the disc varies. In this work we use the relativistic expression derived by Abramowicz et al. (1997), which improved upon the original expression for scale height in Novikov & Thorne (1973), by remaining accurate and non-singular up to and including the event horizon,

SOLVING THE EQUATIONS
The geodesic equations in the Kerr metric are sufficiently complex that they do not possess general analytic solutions (even with no additional four-force present). Therefore it is necessary to resort to numerical methods to find solutions. We require an unconventional numerical method to solve this system of equations and have developed a test particle approach to solve for the motion of a gas particle in the disc. The particle four-acceleration and four-velocity are determined by the forces acting upon it (i.e. turbulent internal torques, pressure gradient forces, radiative forces, etc.), whilst the particle four-velocity and its gradient simultaneously determine the disc properties (surface density, temperature, etc.) as it spirals in through the disc. The set of simultaneous equations we need to solve is: the general relativistic four-acceleration of the fluid eqs. 37 and 38, conservation of mass (16), the energy equation balancing heating and cooling (46) and (47), and the vertical structure of the disc under radiative cooling and hydrostatic equilibrium (48).
To demonstrate the equivalence and validity of the particle-in-disc method developed in this paper to the standard thin disc equations, we present a detailed derivation of the Newtonian thin disc equations using the particle-in-disc method in Appendix B.

Numerical methods
The system of equations describing the inspiralling motion of a gas particle in a relativistic disc characteristically produce quasielliptical oscillations, with a period approximately that of the orbital period, when a force acts on the orbiting particle. For simplicity, we shall refer to these oscillating perturbations from a smooth inspiral as elliptical oscillations, even though relativistic effects mean that the oscillations are not truly elliptical. These oscillations are problematic when attempting to solve our governing equations simultaneously, since both the turbulent four-acceleration and turbulent heating depend sensitively on the velocity shear (and its derivative in the case of the four-acceleration). The averaged velocity shear of the disc gas cannot be sensibly calculated from a particle undergoing radial oscillations superimposed upon a slower smooth inspiral because the radial elliptical oscillations become a dominant component of the calculated velocity shear and its derivative. Furthermore, since the turbulent forces depend on the derivative of the velocity shear (30), without modification this leads to an unstable feedback loop whereby the turbulent forces induce elliptical oscillations that then lead to larger turbulent forces causing larger oscillations etc. Fortunately, in real gas discs, any such bulk elliptical oscillations of gas would be quickly damped by pressure forces and converted into excess heat (which is then radiated away) and would not occur in steady state. Another reason to remove these elliptical oscillations is that whilst we are solving the equation of motion for a fluid element, the full turbulent disc contains a vast number of individual gas particles, all with different randomised phases of elliptical oscillations and so when we calculate the net radial velocity oscillation of all particles in an annulus of the disc, this would average to approximately zero.
To address the issue of elliptical oscillations and mimic the physics of their damping in a real disc, we have introduced an artificial radial damping force into our radial four-acceleration equation, A r damp , whose purpose is to dampen out radial oscillations, converting the extra energy into heat, whilst not changing the angular velocity of the particle. Inside the ISCO stable orbits are no longer possible and so the damping term is chosen to exponentially decay to a negligible value. The additional radial damping term is chosen to be where b is a constant that sets the magnitude of the damping and our solutions are insensitive to the precise form of the exponential decay. A legitimate concern is that this damping may be too strong and have a dominant effect on the radial velocity of the infalling particle. To confirm that this is not the case in our solutions, we have checked that increasing or decreasing the damping coefficient, b, by an order of magnitude has no noticeable effect on our results (Fig. 1d), except in the decay of initial transients close to the initial radius of the particle at large distances. Even with an implicit integration method incorporating damping of elliptical oscillations, however, solving the equation of motion to determine the particle velocity simultaneously with the disc properties is numerically unstable. This is because the fouracceleration from turbulent stresses depends on the gradient of the velocity shear. This causes any undamped elliptical oscillations or numerical errors in the four-velocity to grow via a feedback loop (as outlined in the first paragraph of this section). It is inevitable that some numerical errors and (as yet) undamped elliptical oscillations will exist in the four-velocity resulting in this numerical instability.
Our solution to this problem is to decouple the calculation of the motion of the gas particle from the calculation of the accretion disc properties. This prevents the unstable feedback loop between elliptical oscillations in the gas motion and increased disc stresses occurring. This decoupling is achieved by choosing a smooth analytic function for the gas four-acceleration A µ (r), under which the particle's motion and four-velocity can be calculated stably. The calculation is now stable because we have chosen a smooth four-acceleration function in advance, thereby removing its dependence on the particle motion and elliptical oscillations. For the chosen analytic four-acceleration function the particle orbit and four-velocity are calculated by numerically integrating equation 17. The particle four-velocity and its derivatives are then used to compute the accretion disc properties corresponding to this choice of four-acceleration. The actual, self-consistent four-acceleration experienced by gas in a disc with these properties is then calculated from equations 37 and 38. These two values of the fouracceleration can then be compared and iterated, changing the analytic fit to the four-acceleration until it closely matches the selfconsistent four-acceleration calculated directly from the relativistic disc equations. When the fitted four-acceleration function matches the self-consistently calculated four-acceleration this corresponds to a solution to the relativistic disc equations, i.e. the chosen analytic four-acceleration is equal to the actual gas four-acceleration in the steady-state disc solution.
We have tested that the disc solutions calculated using this method do not depend on the numerical integration scheme used (whether implicit or explicit, Bulirsch-Stoer or Runge-Kutta see e.g. Press et al. 1992). The solutions presented in this paper are calculated using an implicit Bulirsch-Stoer integration method with adaptive time steps. To solve the equations implicitly we use a 2D Newton-Raphson method to calculate the four-velocities U r and U φ , from the acceleration equations 37 and 38. U t is calculated analytically from U µ Uµ = −1, except close to the ergosphere at r = 2rg where gtt = 0 and a bisection method is used instead. It is clear that for our method to work we require that the calculated disc properties are not sensitive to small changes or errors in the analytic fit to the gas four-acceleration. It is demonstrated later that this is not a problem.
One advantage of using this dynamic method of calculating the inspiral of a fluid element through the disc, over methods traditionally used to solve slim disc equations, is that the solution smoothly passes through any sonic points and does not form shocks. This is because, in the rest frame of the test particle, the fluid equations do not change when a sonic point is passed through. In fact, in the rest frame, the fluid still remains in local causal contact with the surrounding fluid and there is nothing particularly special about the sonic point as viewed locally. This is because the location of a sonic point is somewhat arbitrary as it can be changed by a Lorentz transformation, but this does not change the local physics, i.e. the sonic point for an observer moving relative to the black hole is different from the sonic point for an observer stationary to the hole. Of course, this does not mean that information can pass outwards through a sonic point because, although the fluid is in local causal contact with surrounding fluid, the speed of information propagating between the fluid elements in the rest frame (the sound speed) is by definition slower than the speed with which the frame is falling into the black hole inside the sonic point. There are no shocks present in our solutions. For shocks to form in our solutions there would either have to be: a natural discontinuity present in the fluid equations, i.e. the inspiralling gas would have to spontaneously form a shock with itself mid-stream, for which there is no physical mechanism in the case we are studying; or else we must assume pre-existing gas at a certain radius with discontinuous properties to the inspiralling gas, which the inspiralling gas could shock against when it is encountered (this is clearly not a sensible method of finding steady-state solutions).
We have chosen to use Boyer-Lindquist coordinates to facilitate comparison with Novikov-Thorne thin disc and slim disc models. However, in Boyer-Lindquist coordinates there is a coordinate singularity at the event horizon and the metric coefficients become ill-defined. To avoid this singularity we stop integrating our solutions just outside the event horizon at r = 1.02revent. In order to extend solutions through the horizon an infalling system of coordinates would be required instead, e.g. Gammie et al. (2003). However, since gas this close to the horizon will be so severely redshifted that it is unlikely to be observable, we have chosen to use Boyer-Lindquist coordinates for this study.

Numerical tests
Since we are not aware of this particle-in-disc method being previously used to solve the accretion disc equations, it is necessary to confirm that it produces accurate solutions. Since the analytic fit to the four-acceleration will not be exact and the calculation of the four-acceleration is subject to numerical errors, for our solutions to be accurate we require that: (a) the solutions do not depend on the precise initial conditions of the gas particle at large distances (convergence); (b) the accuracy of the solutions does not depend on the size of the damping coefficient used to dampen elliptical oscillations in the disc; (c) the disc solutions are not sensitive to small differences between the analytic fit to the four-acceleration and the self-consistent four-acceleration; and (d) Our disc solutions accurately reproduce the Shakura-Sunyaev and Novikov-Thorne disc solutions at radii where their assumptions are valid. (Our chosen analytic fitting function for the four-acceleration consists of a sum of different power laws, with a break mediated by an exponential decay around the ISCO.) In Figure 1 we demonstrate that the numerical method and disc solutions are accurate and meet requirements (a) through (d) set out above. (a) We have confirmed that the numerical solutions are not sensitive to the initial conditions of the gas particle by checking that solutions with different initial particle radii and fourvelocities quickly converge. This is because initial differences from the equilibrium four-velocity at large distances rapidly decay and are damped, as shown in Figs. 1d and 1e. This damping is physical in origin since it occurs both because of the physically motivated damping of elliptical oscillations and because of the thermalviscous stability of the accretion disc equations themselves (in the case of gas-pressure-dominated discs with Kramers/electron scattering opacities), which causes perturbations to the steady-state solution to decay with time. (b) Fig. 1d shows calculated disc solutions that differ only by the magnitude of the damping coefficient, b (see eq. 51). Whilst the damping coefficient obviously affects the time taken for initial transients to decay, it has a negligible impact on the converged solution away from the initial transient. (c) Figs. 1b and 1c display three disc solutions in which only the analytic fit to the four-acceleration, A φ (r), differs. Analytic fits that bound the self-consistent four-acceleration from above and below were chosen, confirming that the disc solution and self-consistent four-acceleration are not particularly sensitive, at any radius, to small inaccuracies in our analytic fit function. The accuracy of our numerical fit to the four-acceleration is demonstrated in Fig.  1a, through a comparison to the self-consistent four-acceleration and the four-acceleration calculated from the Novikov-Thorne so-lution. (d) Our numerical solutions accurately reproduce the standard analytic thin Newtonian disc solutions at large distances and the Novikov-Thorne disc solution outside the ISCO. This is demonstrated in Figs. 2a-2d (see figure caption for details).
Because of the small size of the radial pressure gradient in thin discs, there is no noticeable change in our solutions whether we accurately fit A r , or set it to zero (excluding the damping term, eq. 51). The Newtonian thin disc model we use is the standard Shakura-Sunyaev α-disc model (Shakura & Sunyaev 1973, Frank et al. 1992, with the zero-stress inner boundary chosen to be at the location of the ISCO radius of the disc.

Dynamics of the plunging region
The most common assumption in solving the relativistic disc equations is to assume that inside the ISCO, gas rapidly plunges into the black hole on approximately radial orbits and so contributes a negligible stress or torque, on outer parts of the disc. This rapid radial plunge is typically assumed to consist of gas free-falling inwards from the disc edge at the ISCO with large radial velocities, such that the inspiral into the black hole is very rapid and predominantly radial, e.g. Novikov & Thorne (1973) (but not literally a radial trajectory). Our particle-in-disc method allows us to accurately calculate the steady-state thin disc properties in the plunging region. The inspiral from the ISCO into the black hole horizon is shown for all our solutions in Figure 3. Importantly, we find that the disc plasma does not plunge approximately radially into the horizon at close to the speed of light from the ISCO, but instead executes a steady inspiral consisting of between ∼4 and 17 complete orbits, depending primarily of the spin of the black hole. This is in contrast to the usual assumptions of thin disc models and explains why we should not expect zero stress at the ISCO, because the material at and inside the ISCO is not particularly diffuse (the radial velocity is still subrelativistic) and a large velocity shear is present at and inside the ISCO, see Fig. 5f, as opposed to the shear calculated between relativistic circular orbits that vanish at the ISCO.
The radial velocity and dynamics of the plunging region depend on several factors, including the spin, which determines the location of the ISCO, the accretion rate and black hole mass. The ISCO shrinks as the spin increases from negative to positive values and this influences the radius at which the radial gravitational pull begins to overwhelm the centrifugal barrier, causing rapid radial acceleration in the plunging zone. The acceleration of the radial velocity around the ISCO behaves similarly, independent of the spin as seen in Fig. 2f. The accretion rate and black hole mass mainly affect the radial velocity of the gas as it nears the ISCO. The radial velocity remains decidedly subrelativistic at the ISCO for all the values of spin, black hole mass and accretion rate we tested, with the velocity increasing steadily towards the horizon as expected.
The dynamics of counter-rotating discs (a < 0) are particularly interesting because the angular velocity of the disc material must swap direction to corotate before it enters the ergosphere, as shown in Fig. 2e (by definition all material inside the ergosphere must rotate in the same direction as the black hole). Within the ergosphere the turbulent stresses continue to act changing the angular velocity of the disc plasma, however, this force is always much weaker than the gravitational force and so once inside the ISCO it does not change the path of the infalling material significantly as it falls into the black hole (as shown by the relatively constant value of U φ and Ut for all solutions within the ISCO, shown in Fig. 6).  Figure 1. Numerical tests to confirm the accuracy and convergence of our numerical method and disc solutions. (a)-Our numerical fit to the self-consistent fluid four-acceleration is good, demonstrating that a satisfactory self-consistent disc solution has been found. Furthermore, our solution agrees closely with the four-acceleration calculated from the Novikov-Thorne solution (except close to the ISCO where the Novikov-Thorne solution becomes invalid) which confirms the accuracy of our solution outside of the ISCO. (b)-The effect of multiplying our numerical fit to the four-acceleration component, A φ , by 1.1 or 0.9 is shown on the self-consistent acceleration (37). This demonstrates that the self-consistent acceleration is not overly sensitive to small changes in our fitting function and validates our method of solution. (c)-The surface temperature of the disc is not overly sensitive to changes in the numerical fit to the four-acceleration component, A φ , which demonstrates that small errors in our numerical fit to the four-acceleration (±10%) remain only small errors in the disc parameters (± ∼ 4%T S , or a ± ∼ 16-20% change in luminosity). (d)-Changing the size of the radial damping coefficient b (51) only affects the damping of initial transients in the solution close to the starting radius of the particle and so has no effect on the accuracy of our solutions at radii of interest (which is why the initial radius was chosen to be so far outside the ISCO at r = 10 4 rg). (e)-Outside of the ISCO, where U r is small, we expect Ut and U φ to have values very close to those of circular orbits as used in Novikov-Thorne (U t (NT) and U φ (NT) ). This is seen to be the case, once initial transients have damped away at large radii (r ∼ 1000 − 10 4 rg). (f)-The aspect ratio of the disc H/r is shown for our full solution using both the Abramowicz et al. (1997) prescription for H (50), alongside H/r ratios for the Newtonian and Novikov-Thorne solutions. There is good agreement at radii where Newtonian and Novikov-Thorne solutions are valid, with the disc becoming thinner close to the horizon as gravitational forces increase and Solution 1 cools (Fig. 2a). © RAS, MNRAS 000, 1-27  The angular four-velocity of disc gas for a variety of black hole parameters. The angular velocity steadily accelerates as the particle spirals inwards at large radii, rapidly accelerating as it approaches the horizon. The parameter with the most significant effect on the angular velocity is the spin. It is interesting to see the change in direction from counter to corotation for our negative spin solution (Solution 6), a necessity as it approaches the ergosphere at r = 2rg. (f) The radial four-velocity component for a variety of disc parameters. The radial velocity remains subrelativistic around the ISCO but rapidly accelerates to relativistic velocities once it enters the plunging region.  . The inspiral of an average disc particle from the ISCO to 1.02revent. The disc gas completes many full orbits as it inspirals from the ISCO into the black hole horizon, as opposed to the rapid, approximately radial plunge that is often assumed. This slower plunge combined with the large increasing shear (Fig. 5f) results in a non-negligible amount of heat and stress generated at and within the ISCO. For the range of parameters calculated here the number of full orbits completed ranges between ∼ 4 and 17. Slightly counter-intuitively, a larger ISCO radius and plunging region actually result in fewer completed orbits before the horizon (i.e. higher spin results in more orbits). The solution with negative spin shows the disc gas reverse its direction of orbit before it is forced to corotate within the ergosphere.  . The surface temperature of the disc is shown for a variety of black hole parameters, alongside the Newtonian and Novikov-Thorne solutions for comparison. Away from the ISCO, the Novikov-Thorne solution closely agrees with our full solution as we would expect. This is because the Novikov-Thorne assumptions are valid until close to the ISCO, where the angular velocity starts to differ from a pure circular orbit (Fig. 1e). The surface temperature drops from its peak as the particle nears the ISCO and plunging zone, where the radial velocity increases, the gas becomes less dense and the extracted angular momentum reduces. Since the velocity shear continues to increase through the plunging zone the surface temperature remains non-negligible in this region, increasing or decreasing as the particle approaches the event horizon depending primarily on black hole spin. However, a combination of increasing gravitational redshift and decreasing surface area at smaller radii means that we expect the contribution of luminosity to the overall observed spectrum to be small from the plunging region.  (f) The velocity shear squared for our different solutions Figure 5. (a) A comparison of our full solutions to solutions calculated using thin disc approximations (no heat advection, effects of internal energy, pressure forces, or compression) reveals that the thin disc approximation holds accurately until the plunging region. Inside the plunging region the full solutions maintain a higher surface temperature because they include the reservoir of internal energy of the gas that causes heat to be advected into the plunging region and radiated. (b) Outside the ISCO surface temperature is determined by a close balance between turbulent heating and cooling, whilst inside the ISCO heat advection and vertical compressive heating provide the dominant source of thermal energy whilst cooling via radial expansion (due to radial acceleration) also becomes important. (c) The sonic point, where the thermal sound speed is equal to the radial velocity, occurs very close to the ISCO radius for all of our disc solutions. This is in agreement with relativistic slim disc models. (d) The four-acceleration A φ due to turbulent stresses, radiative angular momentum loss and changes in the inertia/internal energy of the gas. Turbulent stresses dominate the four-acceleration outside the ISCO, with inertial effects dominating at the ISCO and radiative deceleration becoming steadily more important at smaller radii as the emitted radiative flux increases. (e) The radiative contribution to the total fouracceleration in φ, due to angular momentum carried away by radiation, is significant both outside the ISCO, ∼ 5 − 10% at r ∼ 20rg, and in the plunging region even exceeding ∼ 100%, depending on disc properties. This means that the radiative angular momentum loss should not be neglected when dealing with radiatively efficient discs. (f) -The full velocity shear squared σ µν σµν calculated for each of our solutions, which is closely related to the turbulent heating rate (47). The shear squared smoothly rises as radius decreases. Small differences appear between solutions as the ISCO is approached that are primarily due to different black hole spins.  Table 1. Accretion disc parameters for the solutions presented in this paper where a is the dimensionless black hole spin (−1 < a < 1), M is black hole mass, M mass accretion rate of the black hole, the Shakura-Sunyaev α-parameter is defined in equation 30, r ISCO is the radius of the innermost stable circular orbit (ISCO), and revent is the radius of the event horizon of the black hole. M is the solar mass andṀ Edd is an estimate of the Eddington accretion rate chosen to have a value ofṀ Edd /M = 7.02 × 10 −16 s −1 . Dominant opacity refers to the opacity dominant in the inner disc close to the ISCO, either Kramers or electron scattering. T S ISCO and T S max are the surface temperatures at the ISCO and maximum surface temperatures as measured in the rest frame of the disc gas respectively. (b) Figure 6. The specific angular momentum U φ and specific energy Ut of gas particles are calculated for the different solutions. Inside the plunging region these quantities are almost constant meaning that the gas is close to free-fall with gravity dominating over the acceleration caused by disc stresses (since U φ and Ut are conserved along geodesics).

Full solutions to the relativistic disc equations
Now that we have verified the accuracy of our method and solutions and established that the inspiral into the black hole from the ISCO is not the dramatic radial plunge usually assumed but a more gradual inspiral, it remains to calculate just how much radiation is emitted in the plunging region and what the disc properties are.
The central temperature, Tc, surface density, Σ, and surface temperature, TS, for Solution 1 are calculated in Figures 2a, 2b and 2d. The Newtonian and Novikov-Thorne solutions are also calculated to show where differences appear between the full solution and the two most commonly used thin disc solutions, illustrating where the assumptions of non-relativistic effects and relativistic circular orbits break down. To be clear, our purpose in choosing to include the standard Newtonian disc is not to criticise its accuracy in comparison to more complex relativistic models, especially since its elegant simplicity and accuracy make it well suited for fitting observations of many astrophysical accretion disc systems, but simply to demonstrate how and where our solutions differ from the existing most widely used and best understood disc models. Fig. 2c demonstrates that relativistic effects start to become important for radii r 1000rg, resulting in differences appearing between the Newtonian disc solution and relativistic models. The Novikov-Thorne solution remains accurate outside of the region immediately around the ISCO, r 1.1rISCO, see Fig. 2d. This is expected because the assumption of circular orbits used in the Novikov-Thorne solution only breaks down close to the ISCO, as shown in Fig. 1e. In the Novikov-Thorne solution all of the disc properties either diverge to ±∞ or become 0 at the radius of the ISCO and the solutions cease inside the ISCO. Our full solutions remain finite, smooth and continuous through the ISCO down to the event horizon. The most important result is that the temperature and stress at the radius of the ISCO are finite and non-negligible.
Whilst the disc luminosity and stress do fall off around the ISCO, the surface temperature at the ISCO is still approximately one-quarter of the maximum surface temperature of the disc, see Fig. 4 and Table 1. This intermediate surface temperature in the plunging region is significant if the plunging region were to be resolved and observed on horizon scales (Event Horizon Telescope Collaboration 2019) but is unlikely to produce a significant contribution to observable flux (F = σT 4 S ) if observations are blended with hotter material outside the ISCO. Because of the sudden acceleration of the radial velocity and subsequent rarefaction of the disc gas around the ISCO, the surface temperature, central temperature and surface density all initially decrease when entering the plunging region. However, away from the ISCO region the behaviour of the surface temperature shows surprising sensitivity to the spin of the black hole inside the plunging region with the surface tempera- (1/r).∂ r (rU r ) -(1/r).∂ r (rU r ) -θ z r ISCO Figure 7. A breakdown of the radial and vertical components of the dimensionless scalar compression, θ = ∇µU µ , for Solution 0, into components associated with radial expansion, ∇tU t + ∇rU r + ∇ φ U φ = (1/r)∂r(rU r ), and vertical compression, θz = ∇zU z , see eq. 28. The radial component is weakly compressing at large distances away from the ISCO since the radial velocity increases less rapidly than r −1 (Fig. 2f). The radial velocity accelerates rapidly as it approaches the ISCO, resulting in rapid radial expansion and leading to substantial cooling by pressure forces doing work on the expanding gas (Fig. 5b). The vertical component is always negative (compressive) since the disc height H continuously decreases towards the horizon resulting in vertical compression and compressive heating (Fig. 1f). Vertical compression is a small effect at large radii but increases rapidly towards the horizon as the radial velocity and vertical gravitational force increase. Deeper in the plunging region radial expansion is offset by increasing vertical compression, causing θ to change sign, resulting in significant compressive heating (see also Fig. 5b).  Figure 8. A comparison of popular approximations to the full velocity shear tensor contraction used by slim disc models (see eqs. 53 and 55). The approximations all assume that only the local rφ components are important and so become inaccurate inside the ISCO, diverging substantially from the true value, σ µν σµν , in the plunging region by either passing through zero or becoming infinite. This is important because the volumetric heating rate by turbulence in the disc is 2ρνσ µν σµν and an inaccurate heating rate will lead to inaccuracies in the calculated disc properties and emitted spectrum.
ture rising towards the horizon for more negative spins and falling for higher positive spins (Fig. 4). The rate of heating and cooling becomes complicated in the plunging region where heating via vertical compression, the advection of heat, and cooling via radial expansion also become important in addition to the usual turbulent heating and radiative cooling processes, see Fig. 5b. Although the turbulent heating remains similar in all of our solutions (Fig. 5f), the radial extent of the plunging region and the radial and vertical velocity profiles depend on spin, giving rise to the different surface temperature behaviour in the plunging region.

The value of the inner stress
In our solutions we calculate the relativistic equation of motion of an average fluid parcel subject to turbulent stresses and so we are able to extend the equilibrium relativistic thin disc solutions through the ISCO to the event horizon. One of the main purposes of this work is to demonstrate that for a relativistic thin α-disc, with no net magnetic field, the stress at the ISCO is not an undetermined parameter but has a calculable and definitive value that can be used when fitting spectral models to observations. The surface temperature of the disc is calculated in Fig. 4 for a variety of accretion disc parameters and the surface temperature at the ISCO for each solution is shown in Table 1. The free parameter in the Novikov-Thorne solutions set at the inner ISCO boundary can be expressed in terms of a stress, a surface temperature or any other disc parameter, since these are all related by the relativistic equations. For the purpose of convenience and clarity we choose to express this free parameter/inner boundary condition in terms of the ratio of the surface temperature at the ISCO to the maximum surface temperature of the disc. This is because in radiatively efficient thin discs the disc luminosity and surface temperature are independent of the detailed disc physics (such as which opacity or pressure prescriptions are used), being determined solely by conservation of angular momentum and energy. The surface temperature is also the most useful quantity when calculating the disc spectrum and so a ratio of surface temperature at the ISCO to the maximum surface temperature is a convenient and intuitive way to quantify the size and relevance of the stress at the ISCO. In Appendix C the form of the inner boundary condition and the relation between its expression in terms of the stress and surface temperature for an NT disc is calculated, for those wanting to switch between surface temperature and inner stress.
In Fig. 2d the surface temperature is displayed in a magnified region around the ISCO so that the surface temperature at the ISCO can be easily determined. The significant result is that for our chosen range of disc parameters the surface temperature at the ISCO is between 0.15 and 0.32 times the maximum surface temperature of the disc, see Fig. 4. This means that the luminosity from the ISCO region and plunging region is small compared with the flux (F = σT 4 S ) emitted by the most luminous disc regions but not negligible, especially considering that horizon scale imaging is now becoming possible with projects such as the Event Horizon Telescope (Event Horizon Telescope Collaboration 2019). In Table  1, the surface temperature at the ISCO and its ratio to the maximum surface temperature are shown for the different accretion disc solutions. It is somewhat surprising that the ratio of surface temperature at the ISCO, to the maximum surface temperature of the disc, varies only by a factor of 2 between solutions spanning a wide range of spins, black hole masses and accretion rates. In fact, the ratio TS ISCO/TS max seems only very weakly dependent on the accretion disc parameters: it increases very slightly from 0.281 to 0.312 as the spin goes from -0.9 to 0.99; it decreases from 0.298 to 0.218 when the accretion rate drops by a factor of 100; and it decreases from 0.298 to 0.147 when black hole mass is increased by a factor of 10 8 .
The relatively weak dependence of the surface temperature at the ISCO on accretion disc parameters means that this ratio is not particularly sensitive to accretion disc parameters and so a reasonable estimate of this ratio can be made based on interpolating our limited results (this can be refined in future when a wider variety of parameters have been investigated). This estimate, or better, a full relativistic solution with the correct parameters, can be used to set the inner stress parameter when fitting a relativistic thin disc model to spectra. This will improve the accuracy of black hole parameters estimated by spectral fitting compared to allowing an undetermined arbitrary value of the inner stress, although it is worth emphasizing that our calculation is based on the assumption of a disc with no physically significant large net magnetic fields or disc winds. We leave a detailed calculation of the observed spectrum to a future work since this involves relativistic ray-tracing, however, it seems reasonable to comment that spectral models based on a relativistic NT model with zero inner stress are likely to be a good approximation to our solutions (excluding the plunging region), whereas models with a significant inner stress/temperature do not accurately represent thin discs with moderate accretion rates and without significant large-scale net magnetic fields.

Do the thin disc assumptions hold in the plunging region?
Analytic thin disc solutions have been remarkably successful and influential in our understanding and modelling of black hole accretion discs. It is thought that for observable discs with moderate accretion rates (∼ 0.001 − 0.1Ṁ Edd ) thin disc assumptions hold reasonably well, with advection becoming important at higher accretion rates (e.g. slim discs, Sądowski 2009), and inefficient electron-ion thermal coupling leading to radiative inefficiency and advection-dominated accretion flows (ADAFs) at lower accretion rates,Ṁcrit ADAF ∼ 0.2α 2 (Yuan & Narayan 2014). In this paper, we are interested in extending thin disc solutions down to the event horizon and so we have naturally restricted the accretion rates of our solutions to the regime in which we expect thin discs to be a good assumption. However, since previous thin disc solutions have generally ignored the plunging zone, it is sensible to check when and where these assumptions hold and to what extent the assumptions underlying our equations and solutions remain accurate.
Let us outline the relevant thin disc assumptions: (1) the disc is geometrically thin (H/r 1), this implicitly underpins many of the other assumptions; (2) the heating/cooling time-scale is much smaller than the infall time-scale and efficiently radiates all heat away as it is generated (radiatively efficient, no heat advection or significant internal energy); (3) the disc material is optically thick so that the heat is radiated as an approximate blackbody; (4) pressure gradient forces can be neglected since the gas sound speed is much smaller than the orbital speed; and (5) heating/cooling by compression/expansion is negligible.
(1) -In Fig. 1f we show the H/r ratio as a function of radius throughout the disc. It is clear that the disc remains thin throughout, even in the plunging zone (this is true for all of our solutions). Since the ratio of the gas sound speed to orbital speed is given approximately by the ratio H/r, this means that the gas sound speed is always much smaller than the orbital speed.
(2) -Heat advection is unimportant outside of the plunging zone, however, inside the plunging zone the increased radial velocity means that a substantial amount of heat is advected. This is shown by the rate of change of internal energy of the gas ΣdΠ/dτ in Fig. 5b.
(3) -For all of our solutions the optical depth of the disc is larger than 1, so that it is optically thick, even in the plunging region. In order to make accurate comparisons to the Shakura-Sunyaev and Novikov-Thorne solutions, we have similarly chosen only one opacity to be dominant for the entirety of the disc, see Table 1.
(4) -Radial pressure gradient forces have no significant effect on the solutions, as expected for thin discs, since the pressure gradient forces are always much smaller than the gravitational forces and have a negligible effect on our solutions. This is closely related to the reason why the discs are thin-the sound speed is much smaller than the orbital speed. The disc pressure has a maximum outside the ISCO, roughly when the temperature is a maximum, so that if the pressure gradient were to be important, we would expect it to influence the disc solution around the pressure maximum rather than inside the plunging zone.
(5) -Outside of the ISCO where the radial velocity is subsonic, compressive effects are negligible as usually assumed in thin discs. However, inside the plunging region, where the radial velocity rapidly accelerates, both cooling via radial expansion and heating via vertical compression of the disc become very important, as displayed in Figs. 5b and 7.
To conclude, it is found that many of the thin disc assumptions hold surprisingly well throughout the entire disc for moderate accretion rates, with the two exceptions being that heat advection and heating/cooling by compression/expansion become important within the plunging region. The similarity outside the ISCO and differences within the plunging region between our full solutions and our solutions retaining only thin disc physics are demonstrated in Fig. 5a.

Comparison to relativistic slim disc models
The relativistic conservation equations for rest mass and stressenergy are in principle the same as those derived and used in relativistic slim discs. Our work differs primarily in: (1) our method of solution; (2) our retention of the momentum transported by optically thick quasi-blackbody radiation and gas internal energy; (3) our inclusion of vertical compression; and (4) our evaluation of the full velocity shear tensor to calculate the turbulent disc stress and heating rate. It is therefore worth calculating the significance of these differences, to find out to what extent our solutions and results differ from those of previous relativistic slim disc investigations.
(1) -For a given mass accretion rate slim disc solutions are specified by the location of any sonic points (and its regularity condition), the angular momentum flux passing through the horizon, and appropriately joining on to an analytic disc solution at the outer boundary. Solutions are found by requiring constant mass, momentum and energy fluxes through the disc, often using a method of relaxation on a spatial grid (see  for a clear explanation of the equations and method). Our particle-in-disc method solves the full equation of motion (four-acceleration) of a gas particle given an accretion rate, black hole mass and spin. Our method has the advantage that it converges to a unique stable steady-state solution because it is dynamic and the accretion disc equations are stable (see discussion in Section 4). This means that our solutions have no additional free parameters and naturally pass regularly through any sonic points and converge to physically sensible NT disc solutions at large radii. The downside is that our numerical method is relatively complicated due to the inherent numerical instability of the equation of motion due to elliptical oscillations (see Section 4.1) and so our solutions require iteratively fitting a smooth four-acceleration function to the selfconsistent four-acceleration. Despite these differences both methods should agree when the same conservation equations with the same physical terms are solved, provided the outer boundary condition is appropriately specified for the slim disc.
To illustrate these core similarities, in our solutions the sonic point of the disc occurs at or just inside the ISCO, see Fig. 5c. This agrees with the results of previous relativistic slim disc models Sądowski (2009) and Abramowicz et al. (2010), which found that for thin disc solutions the sonic point occurs very close to the ISCO due to the rapid radial acceleration that occurs here.
(2) -The momentum transported by radiation is neglected in most slim disc solutions Lasota (1994), Abramowicz et al. (1996), ), Igumenshchev et al. (1998), Sądowski (2009), Sądowski et al. (2011, etc. (we are not aware of a case where angular momentum transport by radiation from an optically thick disc is included). This is generally a sensible assumption, considering that slim disc solutions were created to deal with radiatively inefficient thick discs and it also simplifies the equations, however, in radiatively efficient thin discs this term is important and retained (Novikov & Thorne 1973). As the amount of heat generated per unit rest mass in the gas increases at smaller radii and the temperature and radiative flux increase, the effect of angular momentum loss by radiation in thin discs becomes more important. Fig. 5e demonstrates that this effect becomes important even at comparatively large radii ∼ 20rg affecting the A φ component of the four-acceleration of the plasma by ∼ 5 − 10%. Within the plunging region the importance of this effect generally increases affecting A φ by up to ∼ 100%, depending on the disc parameters of our solution. Neglecting radiative angular momentum loss in relativistic slim disc models will therefore introduce inaccuracies when calculating solutions for radiatively efficient thin discs. Figure 5d also demonstrates that the four-acceleration due to the rate of change of internal energy of the gas (change in gas inertia) becomes important in the plunging region. This is caused by the non-negligible momentum associated with the internal energy of the gas and the resulting four-acceleration that acts in response to changes to the internal energy in order to conserve the total fourmomentum of the gas and radiation.
(3) -Vertical disc velocities are usually neglected in relativistic thin and slim disc solutions because the detailed vertical structure of the disc is not solved dynamically and the vertical structure is assumed to be in hydrostatic and thermal equilibrium (e.g. Abramowicz et al. 1996). The effect of vertical velocities on turbulent stresses and turbulent heating is relatively small (see Section 3.2.4), however, heating caused by the rapid vertical compression of the disc in the plunging region becomes the dominant heating source close to the horizon, see Figs. 5b and 7. Vertical compression is a small effect outside of the plunging region but increases rapidly towards the event horizon, as vertical gravity increases, compressing the disc height more rapidly, i.e. Fig. 1f. The heating from vertical compression starts to dominate over cooling from radial expansion for Solution 0 when r 3rg and this leads to the sign change in the total compressive heating, see Fig. 5b. This shows that vertical compression becomes an important effect when calculating heating and luminosity in the plunging region but can be safely neglected outside of the ISCO.
(4) -In the seminal paper by NT, the relativistic form of the velocity shear tensor and relativistic extension of the alpha model for disc stress were laid out. At large radii the velocity shear is dominated by the radial derivative of the orbital velocity and so in NT only the local rφ components of the stress and shear tensors were retained. The local rφ velocity shear σ (r)(φ) (where the brackets indicate the local spatial coordinates are not cylindrical but Cartesian and are aligned with the lab frame φ and approximate r direc-tions) is obtained by Lorentz boosting from the local Minkowski rest frame to the lab frame. This process can be formalised using tetrads, i.e. σ (r)(φ) = σµν e µ (r) e ν (φ) , see e.g.  for more details. The velocity shear square scalar found by contracting the local shear tensor, in the case where only the rφ and φr components are non-zero, is where η µν is the Minkowski metric. This scalar is the simplest and most relevant quantity to allow comparison between this approach and the full velocity shear tensor that we calculate. NT calculated approximate formulae for σ (r)(φ) that are valid outside the ISCO for thin discs where the total Lorentz factor is γ = U t (DNT/ANT) 1/2 (Novikov & Thorne 1973) and Another approximation set out in NT is (56) These NT assumptions were imported into slim disc models with only the local rφ components of the velocity shear assumed to be significant and the approximate NT formulae in eq. 53 or 55 used to calculate them. Equation 53 is probably the most widely used, as set out in the original relativistic slim disc papers by Lasota (1994) and Abramowicz et al. (1996), Equation 55 is used more recently in the work by e.g. Sądowski et al. (2011) (more precisely, here the turbulent heating rate was assumed to be ∝ σ (r)(φ) αP ). In the comprehensive work by , a full calculation of σ (r)(φ) was carried out, although other components of the tensor were still neglected. In this work, we have calculated the full velocity shear tensor and made no assumptions about which components may be neglected. A comparison of the contraction of the shear tensor with itself, which is proportional to the turbulent heating rate, 2νΣσ µν σµν , (47) is shown for these various approaches in Fig. 8. To facilitate comparison we have also calculated the local rφ components of the velocity shear tensor corresponding to our full lab frame shear tensor σ (r)(φ) using tetrads. The comparison of our σ (r)(φ) with σ NT1 demonstrates that it is indeed a good approximation to the local rφ shear tensor, however, clearly neglecting all other components except rφ leads to large errors within the plunging region, as shown by the difference with the full velocity shear square. Figure 8 demonstrates that compared to the full tensor contraction used in this work, σ µν σµν , the various approximations used in relativistic slim disc models all start to diverge away from the full shear square at approximately the radius of the ISCO, becoming very inaccurate and unreliable within the plunging region. The main problems caused by these approximations are that the heating rate is inaccurate and either incorrectly decreases to zero in the plunging region, or diverges to infinity, and that the turbulent stresses (angular momentum transport) also become inaccurate by neglecting all components except the local rest frame rφ components. The approximation σ NT2 is clearly less accurate than σ NT1 and becomes infinite when CNT → 0, causing the heating rate and stress to become infinite at some radius within the plunging zone (e.g. CNT(r = 3rg, a = 0) = 0). This means that it is an acutely poor approximation to use in the plunging region. Whilst it is argued that the value of the turbulent stress and heating rate in the plunging region has little effect on thick discs, where heat transport is dominated by advection ), this will not remain true for radiatively efficient discs. This is because radiatively efficient discs quickly radiate away the heat generated by disc stresses, so that the disc temperature more closely follows the heating rate. This means that if the heating rate becomes inaccurate, or drops to zero, this will lead to inaccurate thin disc solutions.
In this work, we have calculated the full velocity shear tensor and made no assumptions about which components may be neglected. This resolves the problem of the zero or infinite turbulent heating rates found when using only the local rφ components of the shear or approximations (Fig. 8). In contrast, our turbulent heating rate is always positive and well behaved, as shown in Figure 5f. What causes the large differences in Fig. 8? Contrary to the usual assumption that only the rφ components of the shear are significant, we find that most components of the velocity shear tensor, excluding z components, become important inside the plunging region where relativistic effects such as time dilation, length contraction, frame dragging, etc. kick in to spoil our Newtonian intuition. This means that in the plunging region the direction of the velocity shear measured in the rest frame is rotated and is no longer aligned to the lab frame rφ direction. Therefore σ (r)(φ) gives only a component of the rotated total shear and so it can pass through zero and become negative depending on the size of the rotation. The fact that σ (r)(φ) becomes negative in the plunging region for thin discs has been known for some time (e.g. Anderson & Lemos 1988, Abramowicz & Prasanna 1990, however, the importance of other components of the shear tensor in the plunging region has not been previously demonstrated in slim or thin disc models (to our knowledge). Therefore an accurate calculation of the four-acceleration and turbulent heating rate requires the full velocity shear tensor to be evaluated. This problem means that previous relativistic slim disc calculations of the stress and luminosity in the radiatively efficient thin disc regime became inaccurate in the plunging region, with the degree of inaccuracy depending upon the particular form of σ (r)(φ) used.

DISCUSSION
In this paper, we have calculated the properties of a fully relativistic steady-state thin accretion disc down to the event horizon of the black hole. This has allowed us to determine the value of the turbulent stress and luminosity produced at and inside the ISCO. Previous authors have reasoned that for a thin disc the inner stress at the ISCO will be negligible because material at the ISCO lacks centrifugal support and so will plunge rapidly, almost radially into the black hole, e.g. Novikov & Thorne (1973), Abramowicz et al. (2010). In this scenario, the very rarefied gas present at the ISCO would exert minimal outwards stress. We have shown that this picture is well reasoned but slightly exaggerated. The stress at the ISCO is indeed rapidly dropping compared to larger radii due to the rapid radial acceleration of gas near the ISCO, however, the radial velocity is still very much subrelativistic at the ISCO and the plunge from the ISCO into the event horizon takes many orbits (∼ 4 − 17), instead of being close to radial. This means that the disc stress and luminosity are small but not vanishing around and inside the ISCO, with a small but non-negligible amount of radiation emitted from inside the ISCO, within the plunging region. For a relativistic thin disc with no significant net magnetic fields, the inner stress is not arbitrary and should be set to a value such that the surface temperature at the ISCO is approximately 15 − 30% of the maximum surface temperature, depending primarily on the mass accretion rate and black hole mass (see Table 1). The calculation of the value of inner stress for a relativistic thin disc is an important result because the chosen value of inner stress has a significant impact on the calculated spectrum of the disc and the black hole properties inferred from observations .
To be able to extend the relativistic thin disc equations through the ISCO requires dropping the assumption of circular orbits and calculating the equation of motion of an average gas particle in the disc subject to turbulent stresses/torques. To do this we solved the relativistic equation of motion for a gas parcel, whilst simultaneously conserving rest mass, momentum and energy, as the gas spirals inwards through the disc into the black hole. This new particlein-disc method allows the disc properties to be calculated inside the ISCO up to the horizon. Using our full solutions we confirmed that the Novikov-Thorne solutions remain accurate down to within ∼ 5% of the ISCO radius, where the assumption of circular orbits starts to break down. We found that the Novikov-Thorne solution with an inner surface temperature ≈ 0.15 − 0.3TS max is a good approximation to the full solution outside of the ISCO, though the solution inside the plunging region is entirely omitted in the Novikov-Thorne solution. In a future paper, we will calculate the observed luminosity of the plunging region for a variety of disc parameters.
Our relativistic particle-in-disc method has advantages and disadvantages when compared to 3D GRMHD simulations. Thin discs are much more problematic and computationally expensive to simulate in 3D than thick discs. This is because the spatial resolution required to resolve the vertical structure of a thin disc is much higher and the viscous evolution time-scale longer. This means that truly thin discs (H/r 1) are currently not simulated in global 3D MHD simulations. Similarly, it is too computationally expensive to run a large radial disc for a sufficient number of orbits to reach a true steady-state equilibrium. Disc properties are sensitive to the initial condition of the mass reservoir, usually a torus, and any net magnetic field threading this torus. Although the value of net magnetic field is not known, it has a dominant effect on both disc turbulence and jet launching, e.g. Beckwith et al. (2008a). This means that 3D GRMHD simulations are, as yet, not well suited to simulate equilibrium thin discs. In contrast, our method is able to quickly and accurately calculate the full relativistic steady-state solution for a thin disc. The code takes ∼20 seconds to run on a single CPU core (∼ 0.006CPU hours) compared to many millions of CPU hours required for a high-resolution 3D GRMHD simulation. However, the price we pay is that of vastly simplified dynamics, 1D versus 3D, with the properties of MHD turbulence boiled down to a single parameter, α. Since MHD turbulence and α are not numerically converged in 3D simulations as spatial resolution is increased (Fromang & Papaloizou 2007) and both 3D turbulence and the dynamics of the inner magnetosphere depend sensitively on the unknown initial condition of net vertical field in the disc, the most problematic aspects of 1D models are not solved by moving to 3D and remain problematic issues in large 3D simulations as well. This suggests that there remains a useful and complementary role for relativistic 1D disc models, particularly in fitting to observations.
In this work, we retained all of the slim disc effects such as heat advection, pressure gradients, etc. since it was not clear initially whether these effects would become important in the plung-ing region. For a range of moderate accretion rates, ∼ 0.1Ṁ Edd − 0.001Ṁ Edd , we have shown that the optically thick, geometrically thin assumptions remain accurate for the entirety of our disc solutions. However, inside the ISCO, in the plunging region, heat advection and compressive heating/cooling become important. Our work makes several improvements over previous relativistic slim disc models: by calculating the full velocity shear tensor, instead of only retaining the local σ (r)(φ) component, we resolve a problem in relativistic slim disc models whereby the turbulent heating rate becomes inaccurate and falls to zero in the plunging region; the angular momentum carried away by radiation from an optically thick disc is included, which is important for radiatively efficient thin discs; the effects of four-momentum carried by internal energy of the gas and vertical compression are also included and shown to be important in the plunging region. In agreement with relativistic slim disc investigations we find that the disc sonic point, in the case of thin discs, occurs very close to the ISCO.
Finally, an important limitation of α disc models is that the standard thin disc equations describe a turbulent zero-net flux magnetic field. If a significant net poloidal magnetic field is present in the disc, it can dramatically affect the strength of turbulence, mass and angular momentum loss via magnetocentrifugally launched winds and can be advected through the disc, potentially contributing dynamically significant magnetic pressure and stresses in the most extreme case of a magnetically arrested disc (Narayan et al. 2003). Needless to say that these complex 3D effects are not captured in standard 1D relativistic thin disc models and are best calculated by 3D GRMHD simulations. From simulations it seems that for thin discs with weak or zero net poloidal magnetic fields and modest accretion rates, MHD winds are weak, with turbulent magnetic fields acting similarly to analytic relativistic thin disc assumptions Penna et al. 2010). This suggests that for such discs our solutions should be accurate. It is not clear whether discs with small net magnetic fields are common in the universe and so investigating solutions with dynamically important net magnetic fields and winds is an important next step.

CONCLUSIONS
The highly influential Novikov-Thorne relativistic thin disc solutions are only valid up to the ISCO, where the assumption of circular orbits breaks down. This leads to an undetermined inner boundary condition at the ISCO corresponding to an unknown inner temperature and stress of the disc. The chosen value of this inner stress has a significant impact on the calculated spectrum of the disc and the inferred black hole properties . By developing a new particle-in-disc method that calculates the relativistic equation of motion of an average gas parcel as it spirals inwards through the disc, we have calculated full relativistic disc solutions, extending through the ISCO up to the black hole horizon. This eliminates the arbitrary boundary condition at the ISCO and allows us to calculate the correct value of luminosity and stress at the ISCO and throughout the plunging region for a steady-state thin disc.
It is found that the stress and surface temperature at and inside the ISCO are small but non-negligible, with a surface temperature at the ISCO roughly ∼ 0.15 − 0.3 times the maximum surface temperature of the disc, depending primarily on accretion rate and black hole mass. This demonstrates that for a thin disc model the inner stress has a definitive, calculable value, which should be used when fitting spectral models to observations. We calculate the inner stress/inner surface temperature for a range of accretion disc parameters with modest accretion rates and zero net magnetic flux. Using our full solutions we confirm that the Novikov-Thorne solutions remain accurate down to within ∼ 5% of the ISCO radius, where the assumption of circular orbits starts to break down. The Novikov-Thorne solution, with the addition of an inner stress/surface temperature at the ISCO chosen to match the value calculated in our full solutions, is found to be a good approximation to our solutions down to the ISCO radius, though the plunging region inside the ISCO is entirely omitted in the Novikov-Thorne solution.
It is often assumed that the disc material makes a rapid radial plunge from the ISCO into the black hole, which would lead to no significant stress at the ISCO and justifies the choice of zero stress at the ISCO. In contrast, we calculate that disc gas completes many (∼ 4 − 17) full orbits inside the ISCO before falling into the black hole, with a subrelativistic radial velocity at the ISCO. It is because of this gradual inspiral that the stress and luminosity of the disc remain non-negligible at the ISCO and within the plunging region. Our method allows rapid, accurate calculations of the dynamics and properties of a full steady-state relativistic disc, including the plunging zone, without the need for expensive 3D GRMHD simulations. Our model improves upon relativistic slim disc studies by resolving a problem with inaccurate turbulent heating and stresses in the plunging region and by including the angular momentum carried away by optically thick radiation, both of which are important in the case of radiatively efficient thin discs. This method can be used as the basis for accurate dynamic and spectral relativistic disc models.

ACKNOWLEDGEMENTS
WJP thanks the referee, Ramesh Narayan, for his many helpful comments that improved the clarity of the paper. WJP acknowledges support from the University of Oxford.

DATA AVAILABILITY
Data available on request. by.
where the number density has been modified from eq. A1 because the number density of particles differs between that measured in the rest frame of the gas and that measured in our own rest frame because of length contraction of the gas. Again it is easy to verify that in our rest frame the stress-energy has the desired properties, with the correct energy density and momentum density. For particles with rest mass m, the energy density is T tt = nγm (the number density of particles multiplied by the individual particle energy), and the momentum density is T ti = nγmu i (the number of particles multiplied by the individual particle momentum in i, where the lower case u refers to the spatial 3-velocity). The stress-energy tensor for a collection of ideal gas particles is the sum of individual particle stress-energy densities (A2) where the sum in parentheses (i) is over N particles in a infinitesimal spatial 3-volume element dV measured in the local rest frame.

A1 Thermal gas (non-relativistic particles)
An isotropic ideal gas is composed of an isotropic distribution of particle velocities, with a thermal distribution of energy/velocities. If the gas particles have rest mass m (i) , with particle 4-velocity and the total stress-energy can be expressed in terms of the averaged properties (averaged over solid angle and velocity distribution, indicated by angled brackets) In the case where particle rest masses are the same m (i) = m, or can be sensibly replaced by an averaged rest mass m and average rest mass density ρ = nm, and using the Lorentz factor of particle i, measured in the fluid rest frame γ (i) = U µ V (i)µ , this simplifies to If the velocity distribution is isotropic then there is no correlation between the three orthogonal spatial components of the 4-velocity so that V i V j γ = 0, where i and j only run over the 3 spatial indices. In this case the stress-energy tensor is easily calculated where eg is the internal energy of the gas, the spatial 3-velocity v i is given by V i = γiv i , which in the isotropic case gives v 2 = v 2 x + v 2 y + v 2 z = 3v 2 x and pg the isotropic pressure. A representation of this tensor in terms of the fluid 4-velocity U µ and metric g µν is which can be verified in the rest frame where the local metric is Minkowski g µν = η µν . In the non-relativistic limit we can expand γ ≈ 1 + v 2 /(2c 2 ) to show that eg ≈ ρ v 2 /2 and pg ≈ ρ v 2 /3, which agrees with the standard result from kinetic theory. Finally, we could allow for f extra degrees of internal freedom in the gas particle (make it a molecule which can oscillate for example). From kinetic theory we know that on average an energy kBT /2 is distributed to every degree of freedom in the gas in thermal equilibrium. This additional energy will appear in the rest mass of the gas molecule such that, γn(m + δm)c 2 = ρc 2 + (f + 3)pg/2. From this we find that the internal energy has the familiar relation to the pressure eg = (f + 3)pg/2 = (f + 3)nkBT /2.

A2 Thermal gas (photons or relativistic particles)
If the gas is a photon gas (or the particles are relativistic) then equation A3 can instead be simplified to where the photon energy is E = P t . For an isotropic ideal relativistic gas (P x ) 2 +(P y ) 2 +(P z ) 2 = 3(P i ) 2 = E 2 , since P µ Pµ = 0, so that with all other components zero. From which it follows that er = 3pr. In the case of a thermal photon gas the radiation pressure has the well-known form pr = 4σT 4 /(3c). The stress-energy of a noninteracting isotropic relativistic gas has the same perfect-fluid form as the non-relativistic case, with these different values of internal energy and pressure.

A3 Transport of 4-momentum
A gradient in the average 4-velocity will cause 4-momentum to be transferred through a gas by the action of collisions as gas molecules from parts of the fluid with different average velocities collide with each other. In the simplest prescription this process causes the gas velocities to become correlated in proportion to the velocity gradient ∇µUν and corresponds to the action of heat conduction and viscosity (see Section 3.2.3 for the breakdown of the spatial velocity gradient into compressional, rotational and shearing components). The momentum transport due to compression is proportional to the bulk viscosity (not important in very sub-sonic incompressible gases but important in shocks), the momentum transport due to rotation is proportional to rotational viscosity (always assumed to be negligible since the moment of inertia and angular momentum of individual gas particles is tiny compared to the fluid angular momentum, hence the stress-energy tensor is assumed to be symmetric, Misner et al. 1973) and the momentum transport due to velocity shear σ µν (25) is proportional to the shear viscosity of the gas µ. In this work we neglect the effects of heat transport by collisions between gas particles with different energies caused by spatial gradients in the average particle energy U t because heat transport is dominated by radiation (or convection) in radiatively efficient thin discs. The stress-energy contribution of viscous shear stresses can be included (Misner et al. 1973) For accretion disc plasmas the molecular gas viscosity and radiative viscosity are insignificant Shakura & Sunyaev (1973) and so will be neglected in this paper µg ∼ µr ∼ 0 (along with the bulk viscosity of the gas since we assume shocks are not present in our steady-state solutions). From eq. A6 we can see that the velocity correlations brought about by molecular viscosity and velocity shear are The photon viscosity µr (also negligible compared to the anisotropic turbulent stresses Shakura &Sunyaev 1973 andPotter &Balbus 2014) can be defined via eq. A12 − 2µrσ ij = n P i P j E . (A17)

A4 Electromagnetic stress-energy
The electromagnetic stress-energy is given by Misner et al. (1973) 4πT µν where the electromagnetic field tensor (Faraday tensor) is F µν and its dual (Maxwell tensor) F * µν = (1/2) µνλσ F λσ can be conveniently defined in terms of the electric field and magnetic field 4-vectors e µ and b µ F µν = U µ e ν − U ν e µ + µνλσ U λ bσ, where the Levi-Civita tensor, µνλσ , is related to the Levi-Civita symbol, ε(µνλσ), by µνλσ = √ −gε(µνλσ), µνλσ = −(−g) −1/2 ε(µνλσ) because it is a tensor density (Misner et al. 1973), where ε(txyz) = 1. These equations can be verified in the fluid rest frame where e µ = (0, E) and b µ = (0, B) where E and B are the usual electric and magnetic field 3-vectors. The stressenergy corresponding to the fields is given by where e 2 = e µ eµ and b 2 = b µ bµ. The field evolution equations are given by the usual electromagnetic equations ∇µF νµ = 4πJ ν , ∇µF * µν = 0, where J µ is the current density 4-vector. In thin accretion discs with moderate accretion rates, the ideal MHD limit is a good approximation see e.g. Wardle (1999). In this case the abundance of free charge in the fluid rest frame means that any electric fields are quickly screened by rearrangements of charged particles so that in the fluid rest frame e µ ≈ 0. In the ideal MHD case the stressenergy of magnetic fields are therefore given by e.g. Komissarov (1999) and Gammie et al. (2003) 4πT µν = b 2 U µ U ν + 1 2 b 2 g µν − b µ b ν .
and therefore require either direct numerical 3D MHD calculations, or a simplifying model such as the α model (Shakura & Sunyaev 1973) to obtain approximate solutions. Combining the stress-energy from the turbulent gas velocity and turbulent magnetic fields, the total turbulent stress-energy tensor is (A30) The standard assumption in accretion disc modelling is that the strength of turbulence is linearly proportional to the velocity shear. It seems likely that this linear proportionality is an oversimplification with dedicated simulations suggesting that α ∝ q n , with the shear parameter q = −d ln Ω/d ln r, and n ∼ 6, see e.g. Pessah et al. (2008) and Penna et al. (2013). In order for our model to be easily comparable to the vast body of existing work and because this point is not the main purpose of this investigation, we will stick to the standard assumption of linear proportionality. (In fact this may not be a poor assumption since for a thin Kerr disc the calculated linear growth rate of the MRI appears to be largely unaffected by relativistic effects and remains close to Keplerian value, with the maximum linear growth rate proportional to the velocity shear square, σ µν σµν , Araya-Góchez 2002, Gammie 2004 and Fig. 8).
That the turbulent stresses should increase with the rate of shear at all is perhaps more obvious for magnetic fields, since they are approximately frozen-in, so that a larger shear will result in more rapid stretching and a larger induced anisotropic correlation in the direction of the velocity shear. For turbulent velocities a similar process acts as fluid is moved across the average velocity shear with typical lengthscales of order the disc scale height, resulting in anisotropic velocity correlations proportional to the velocity shear (analogous to how regular viscosity works and produces anisotropic velocity correlations). It is standard in work on thin discs for the strength of the anisotropic stresses and correlations to be set proportional to α times a maximum effective turbulent viscosity (Shakura & Sunyaev 1973). The maximum effective turbulent viscosity would have a mean free path given by the size of turbulent eddies (whose maximum size is limited by the disc scale height, H) and a maximum speed of these turbulent eddies being approximately the sound speed cs (since supersonic turbulence dissipates efficiently through shocks). This means that we anticipate α < 1 and its value encapsulates our ignorance of turbulence (Shakura & Sunyaev 1973). Under these common assumptions the turbulent stresses are with σ µν the velocity shear tensor and ν the effective kinematic viscosity corresponding to the turbulent stresses.

A6 Stress-energy of heat
The correct relativistic form of the heat and viscosity are contentious because the Newtonian descriptions in general have infinite wave speeds and the most common relativistic prescriptions are the covariant generalisations of the Newtonian expressions Eckart (1940) and Misner et al. (1973). Addressing these problems is beyond the scope of the current work and so the usual relativistic expressions are adopted (for a discussion on the issue of faster than light viscous propagation in discs see e.g. Szuszkiewicz 1994 and). We do not expect these issues to be particularly problematic in our solutions because firstly we are looking at steady-state solutions in which sufficient time has passed for information to have propagated multiple times through the entire extent of the disc and secondly because our numerical method integrates the equations from large to small r so that information only propagates inwards (i.e. no information is passing outwards through the sonic points or event horizon). The standard stress-energy expression for heat was proposed by Eckart (1940) and takes the form where q µ is the heat energy flux 4-vector measured in the fluid rest frame. Following Novikov & Thorne (1973), the heat transfer in a optically thick thin disc is assumed to be primarily due to radiation and the heat primarily escapes the disc vertically where it is effectively radiated away from the last scattering surface of the disc where the vertical optical depth τ = Σκ ∼ 1. Under these assumptions in the fluid rest frame the heat 4-vector has the form q µ (z = ±H) = (0, 0, 0, σT 4 S ). In thin disc models the detailed vertical structure is not solved dynamically and so only the value of the radiation heat flux at the disc surface is required to obtain solutions.

A7 Final results
The total stress-energy tensor of the fluid, in the zero net magnetic flux case, is the sum of the thermal, turbulent and heat stress-energy tensors found by combining equations A11, A14, A30 and A32.

APPENDIX B: PARTICLE-IN-DISC EQUATIONS IN THE NEWTONIAN REGIME
Since our approach to solving the thin disc equations differs through solving the 4-acceleration of an average fluid parcel, it is a useful exercise to explicitly demonstrate that in the Newtonian regime we recover the standard thin disc equations. In the Newtonian nonrelativistic limit the linearised Schwarzschild metric is ds 2 ≈ −c 2 1 − 2rg r dt 2 + 1 + 2rg r dr 2 + dz 2 + r 2 dφ 2 , (B1) when r 2rg. When gravitational forces are calculated the first order terms rg/r are retained, however, in the Newtonian regime deviations from flat space time are minimal, rg/r 1 and so timedilation and length-contraction aren't important when calculating quantities such as the velocity shear tensor so that the cylindrical Minkowski metric is then appropriate ds 2 ≈ −c 2 dt 2 + dr 2 + dz 2 + r 2 dφ 2 .
The only non-zero Christoffel symbols in cylindrical coordinates are Solving the relativistic geodesic equation of motion for orbiting disc gas in the Newtonian thin disc regime, rg/r 1, U r rU φ , necessarily reduces to the equations of circular Keplerian orbits (which are accurate at large distances). The Keplerian angular velocities are as standard U φ = Ω = GM r 3 1/2 , U φ = g φφ U φ = r 2 U φ = (GM r) 1/2 .
calculateU φ Noting that where we have again used the relation E = ΩL . Expanding the derivative in (C14) and substituting in the above result f (E − ΩL) 2 Ω = L (E − ΩL)dr + C, where C is an arbitrary constant determining the stress and radiative flux at the inner boundary of the disc The standard Novikov-Thorne solution corresponds to the solution with a vanishing inner stress which we label wNT (see PT for analytic expression), the inner stress w0 can be changed to an arbitrary value by changing the value of C.
where a subscript zero indicates the value of a quantity evaluated at the inner boundary of the disc (the radius of the ISCO). For a given inner stress, w0, the radiative flux emitted by the disc is From the relation of flux to surface temperature, F = σT 4 S , and using FNT(r = rISCO) = 0 the surface temperature at the ISCO is given by This allows us to set the inner boundary condition via either the surface temperature or stress at the ISCO and to calculate one in terms of the other. It is convenient for our purposes to be able to directly set the surface temperature at the ISCO as opposed to the stress.
In the NT and PT derivation it is assumed that w0 = 0 (i.e. the inner stress is zero) and this remains the most commonly used assumption in modelling disc spectra. In this paper we present calculations extending the solution through the ISCO up to the event horizon for a steady-state relativistic thin disc. This means we are able to calculate the actual value of the disc stress at the ISCO and inside the plunging region. It is therefore useful to be able to compare how well a Novikov-Thorne disc with non-zero stress/surface temperature at the ISCO corresponds to our full solution in the range r ≥ rISCO, and what value of TS(r = rISCO) or w0 is most appropriate to use when fitting relativistic thin disc models to data.