The dynamics of accretion flows near to the innermost stable circular orbit

Accretion flows are fundamentally turbulent systems, yet are classically modelled with viscous theories only valid on length scales significantly greater than the typical size of turbulent eddies in the flow. We demonstrate that, while this will be a reasonable bulk description of the flow at large radii, this must break down as the flow approaches absorbing boundaries, such as the innermost stable circular orbit (ISCO) of a black hole disc. This is because in a turbulent flow large velocity fluctuations can carry a fluid element over the ISCO from a finite distance away, from which it will not return, a process without analogy in conventional models. This introduces a non-zero directional bias into the velocity fluctuations in the near-ISCO disc. By studying reduced random walk problems, we derive a number of implications of the presence of an absorbing boundary in an accretion context. In particular, we show that the average velocity with which a typical fluid element crosses the ISCO is much larger than is assumed in traditional theories. This enhanced velocity modifies the thermodynamic properties of black hole accretion flows on both sides of the ISCO. In particular, thermodynamic quantities for larger ISCO stresses no longer display pronounced cusps at the ISCO in this new formalism, a result with relevance for a number of observational probes of the intra-ISCO region. Finally, we demonstrate that these extended models reproduce the trans-ISCO behaviour observed in GRMHD simulations of thin discs.


INTRODUCTION
The accretion of material onto astrophysical black holes liberates vast amounts of energy and is the process through which some of the brightest objects in the Universe are powered.The modelling of these astrophysical accretion flows represents one of the original probes of the strong field regime of gravity, through which the properties of numerous black holes have been constrained (e.g., Reynolds 2013;McClintock et al. 2014).One such prediction of general relativity, the existence of an innermost stable circular orbit (hereafter ISCO; Bardeen et al. 1972), profoundly modifies the dynamic and thermodynamic properties of accreting fluids at short distances from a black hole.Within this ISCO radius circular motion is unstable to inwards perturbations, and test particles plunge towards the singularity at r = 0.As a fundamental prediction of general relativity, the development of theoretical descriptions of the observational characteristics of this strong field regime may be leveraged in the future to derive tighter constraints on black hole properties (see e.g., Reynolds & Begelman 1997;Wilkins et al. 2020, for an example of this philosophy applied to iron line fitting).For this to be a viable approach we must be sure that our theoretical descriptions of accretion accurately reflect the physical conditions of astrophysical sources.
Accretion flows are fundamentally turbulent systems, a result of the magneto-rotational instability (MRI; Balbus & Hawley 1991).However, all conventional analytical modelling of accretion flows assume that they can be described by an effective "viscous" redistribution of angular momentum (e.g., Shakura & Sunyaev 1973;Novikov & Thorne 1973).These classical models of accretion flows ⋆ E-mail: andrew.mummery@physics.ox.ac.uk must break down on length scales shorter than the typical size of a disc turbulent eddy, which for an accretion flow is macroscopic -of order the disc's scale height H.When there is no radial length scale which probes the disc fluctuation scale (i.e, far out in the main body of the disc where observational diagnostics can be safely averaged over many disc scale heights), the effects of this simplification are likely inconsequential (which is why conventional thin disc theory works so well in this limits).However, as we argue here, once there is a relevant radial length scale with which to contrast with the eddy scale (such as the length scale associated with the inner disc edge), this classical prescription should break down, and new descriptions should be developed which may well result in different, and observationally relevant, predictions.
In an attempt to move away from these classical descriptions, in this paper we examine the properties of a series of reduced random walk models.Systems undergoing a random walk are more mathematically flexible than purely viscous systems, and are well suited to modelling physical systems with both a global diffusive character (like an accretion flow on large scales), but also with large amplitude velocity fluctuations (like an accretion flow at turbulent eddy scales).These random walk calculations highlight that the typical velocity with which a fluid element crosses the ISCO may be orders of magnitude higher than predicted from classical thin disc models with a finite ISCO stress (e.g., Agol & Krolik 2000).It was recently demonstrated (Mummery & Balbus 2023) that the trans-ISCO velocity plays a key role in the thermodynamic evolution of fluid flows inside of the ISCO, and this velocity amplification therefore produces important modifications to these profiles in the near-ISCO region.
While we motivate this calculation on purely physical grounds, it could be equally well motivated from a purely model comparison perspective.Models of accretion with a large finite stress (e.g., Agol & Krolik 2000;Mummery & Balbus 2023) show pronounced cusps at the ISCO.These cusps, or to be more precise discontinuities in the gradients of thermodynamic properties, are almost certainly unphysical, and a result purely of the governing assumptions of classical accretion theory.Indeed, GRMHD simulations of accretion flows display a smooth evolution across the ISCO, with no indication of cusps (e.g., Shafee et al. 2008;Noble et al. 2010;Zhu et al. 2012;Schnittman et al. 2016;Liska et al. 2022), even for those simulations which display large ISCO stresses.It seems likely that ISCO-cusps in analytical models of accretion could produce systematic effects when such models are fit to data (e.g., Reynolds & Begelman 1997;Wilkins et al. 2020), and therefore it is of general modelling interest to examine the physical causes of such behaviour, and improve on the underlying modelling assumptions in such instances.
The layout of this paper is as follows.In section 2 we review the classical physics of accretion relevant for this near-ISCO study.In section 3 we introduce the mathematical random walk models we shall consider in this paper, before solving some explicit reduced problems in section 4. In section 5 we put these results in an astrophysical context, and derive modified global thin disc solutions which take into account the insight gained from the random walk models.We compare these models to some GRMHD simulations of thin discs in section 6, finding good agreement, before concluding in section 8.

NEAR-ISCO ACCRETION FLOWS
In this section we construct an argument based on classical disc theory which highlights that in the presence of absorbing boundaries the classical description of the mean fluid flow of an accretion disc must break down.
The classical theory of extra-ISCO relativistic accretion proceeds by first defining a stress energy tensor which describes the accretion flow and then by constructing mass, energy and momentum conservation equations.The classical accretion disc stress energy tensor is the following (e.g., Novikov & Thorne 1973;Balbus 2017) where ρ is the rest mass density, e the internal energy density and P the pressure of the fluid.The 4-velocity of the flow is u µ , while uµ is its covariant counterpart.The final pair of terms represent the energy-momentum flux carried out of the system by photons, where q µ is the photon flux 4-vector.With this stress energy tensor defined, one solves the equations of mass, angular momentum and energy conservation where in these expressions ∇µ is a covariant derivative with respect to Kerr metric coordinate x µ .These three constraints are sufficient to determine three key quantities: the governing equation for the evolution of the disc surface density, the radial velocity of the flow, and the energy flux out of the upper and lower disc surfaces (see e.g., Novikov & Thorne 1973;Balbus 2017, for discussions and various derivations).
The principal theoretical simplification employed in deriving the thin disc solutions of these coupled equations pertains to a series of approximations regarding the properties of the disc fluid's velocity.In particular, the solutions to these three equations are derived by making the following important assumption: the total disc 4-velocity u µ (as well as uµ) may be decomposed into a mean component U µ and vanishing-mean fluctuating component δU µ : which satisfy asymptotic scalings (Balbus & Papaloizou 1999;Balbus 2017) While the fluctuations are an asymptotic scale larger than the mean radial flow of the disc, they are assumed to vanish on average where ∆t is a time long compared to the timescale upon which turbulent fluctuations are induced in the flow, but much shorter than the timescale upon which the mean disc quantities evolve.Physically, this timescale should be thought of as a few times the orbital period at a radius r.While the fluctuations themselves are assumed to vanish on average, their average correlations are in general non-zero.
In particular, accretion is ultimately driven by the non-zero correlation of the components of the turbulent velocity fluctuations, which produce a turbulent stress tensor W µν : where the angled brackets denote the same averaging procedure introduced above.As the first order fluctuations in the disc velocity vanish on average, and the second order drift velocity is assumed to be extremely small, the zeroth order motion of the disc is well approximated by that of precisely circular motion, i.e., U 0 , U0, U ϕ and U ϕ are equal to the test particle circular motion solutions of the Kerr metric.By making certain assumption about the local properties of W µν , classical relativistic accretion models (e.g., the Novikov & Thorne 1973, Page & Thorne 1974, and Shakura & Sunyaev 1973 solutions) allow the mean second-order radial accretion velocity of the flow at a given radius to be determined as a function of the physical parameters of the system (e.g., the black hole mass and spin, the mass accretion rate, and the disc stress α-parameter).However, it is important to recall that accretion flows are turbulent, and that the typical scale of the turbulent fluctuations are of a different (larger) asymptotic scale than the mean drift velocity U r ≪ δU r . (7) In the main bulk of the disc (far from the ISCO radius rI ), these velocity fluctuations vanish on average In effect this statement follows from the fact that there is no preferred direction in which the turbulent fluctuations occur, and fluctuations (e.g.) outwards in the disc are all compensated by turbulent fluctuations inwards.However, close to the ISCO itself, these velocity fluctuations will develop a non-zero directional bias, with fluctuations across the ISCO in effect absorbed into the black hole due to the lack of rotational support within the ISCO.This means that fluctuations across the ISCO from the main body of the disc are no longer compensated for by fluctuations back from the intra-ISCO region.This favouring of fluctuations in a specific direction will mean that the average defined above will no longer vanish in the near-ISCO region, and instead Therefore, careful attention must be paid to the precise value of the the trans-ISCO velocity used in computing the evolution of the intra-ISCO thermodynamic quantities.
Of course, the ISCO does not represent a truly perfect absorbing boundary in a black hole accretion flow (only the event horizon is truly such a boundary).It is in principle possible for a fluid element to cross the ISCO and then return to the main body of the disc, and in fact to some small degree this must happen in a real accretion flow (the ISCO stresses measured in GRMHD simulations are an angular momentum flux sourced from within the ISCO after all).However, the crux of the argument put forward in this paper rests on the assumption that the ISCO acts sufficiently strongly like a oneway gate in the flow that the fluctuation-averaging integrals do not cancel to zero (eq.5).Given the asymptotic scale separation between radial fluctuations and mean drift, it is only necessary to perturb the precise cancellation in the fluctuation velocity integral to a relatively minor degree for the effects of the boundary to become apparent.
In purely gravitational dynamics the ISCO is indeed a one-way gate.If a test particle is on a circular orbit at the ISCO, and is perturbed infinitesimally in the radial direction, then it's subsequent radial velocity is given by the solution of the relativistic energy equation gµν U µ U ν = −c 2 , or explicitly Mummery & Balbus (2022) Here notation is standard, rg = GM•/c 2 .Clearly there are only real solutions of this constraint for r ≤ rI (i.e., inwards perturbations), and the velocity increases rapidly inwards.For a Schwarzschild black hole U r is already ∼ 0.01c at r = 0.9rI , which is a plausible trans-ISCO turbulent perturbation scale in a relatively thin H/r ∼ 0.1 disc.This radial inflow velocity is already much larger than the typical sound speed in the disc (as will be demonstrated in later sections), and it seems unlikely therefore that the fluid element has a significant probability of returning to the main body of the disc.Of course, if significant non-gravitational forces are present in the flow (for example for the extreme magnetic fields produced in a magnetically arrested disc) this argument may well break down.While bearing in mind the inherent simplifications employed, we will for the remainder of the paper treat the ISCO as a perfect absorber.
The trans-ISCO velocity plays a key role in the thermodynamic evolution of disc quantities in the region surrounding the ISCO (Mummery & Balbus 2023), and an incorrect value for this velocity can lead to unphysical discontinuous behaviour at the ISCO itself.It turns out that this effect is particularly relevant in the physical regime corresponding to larger ISCO stresses.This is because the mean Novikov & Thorne (1973) flow velocity decreases with ISCO stress WI (proof in Appendix A of Mummery & Balbus 2023) An increased ISCO stress on the other hand causes various disc quantities (for example the central temperature) to increase.As the acceleration within the ISCO is driven almost entirely by gravity (and therefore is to leading order independent of the local thermodynamics), classical models of finite ISCO stress discs show cuspy behaviour at the ISCO, which is almost certainly unphysical.See 1 Ṁedd , α = 0.1, δ J = 0.1; see section 5 for a description of the physical meaning of these free parameters).In the upper panel we display various normalised (by their ISCO values) disc quantities on both sides of the ISCO (ISCO displayed by vertical black line).The evolution across the ISCO is not smooth in these models, a result of the very low trans-ISCO radial velocity of the disc (lower panel).We zoom into the inner 2rg either side of the ISCO (black dashed curve, r I = 6rg).
Figure 1 for an example of this cuspy behaviour.Indeed, GRMHD simulations generally show a smooth evolution of disc quantities over the ISCO (Shafee et al. 2008;Noble et al. 2010;Zhu et al. 2012;Schnittman et al. 2016;Liska et al. 2022).
The speed of sound, however, increases with a larger ISCO stress a result of the higher inner disc temperatures.The speed of sound is likely a good measure of the typical turbulent velocity fluctuation scale (in fact this scaling is assumed in classic α-models), and therefore while the typical mean drift velocity will decrease, the typical trans-ISCO velocity may well increase as a function of ISCO stress.An increased trans-ISCO velocity in high ISCO stress discs would remove the cuspy behaviour displayed in Fig. 1.The relevant scale of the trans-ISCO velocity is the focus of the analysis in this paper.
In effect, the argument we are putting forward here regards the breakdown of classical "viscous" descriptions of turbulent fluids when they are within a few turbulent eddy scale lengths of absorbing boundaries.These classical descriptions, independent of the magnitude of the assumed viscosity, fundamentally cannot capture behaviour on scales below the eddy length.When other physics (in this case rapid gravitational acceleration) imprints characteristic length scales on the problem, this finite length scale of fluctuations will be imprinted in the flow dynamics, which may well be relevant for observational modelling of black hole accretion flows.

RANDOM WALK MODELS: THEORETICAL SETUP AND GENERAL CONSIDERATIONS
Having identified that the presence of an ISCO (or more generally any absorbing boundary) will act to modify the typical velocity scale with which a fluid element is observed to cross radii close to that boundary, we move to a toy model framework which allows us to , where we plot the position of the particle X k after k jumps.The dashed line shows the final jump of the particle, and the dot shows the position of the particle after this last jump.The arrow shows the direction of the imposed drift (the average of the jump distribution, ϵ), ensuring that all particles eventually cross X = 0. Clearly random walk models can incorporate the physics of large amplitude velocity fluctuations.Lower: The evolving (numerical, averaged over 10 4 trajectories) probability density function of a group of evolving random walk particles at different time-steps k.When averaged over large numbers of particles random walk calculations reproduce properties of diffusive evolution (in the limit of large k, with X 0 / √ k finite).
probe directly these effects in a controlled manner.The toy models we turn to are random walks.Random walk systems are well suited to modelling systems with diffusive properties, but also systems with large scale velocity fluctuations, as we now discuss.Consider a one dimensional discrete time random walk.These systems consist of a single particle moving in one dimension, in steps of finite time duration ∆t.The position of a particle currently at position X k after the next time step is given by where the velocity of each jump is sampled from some distribution An example of 10 particle trajectories, which all start at X0 = 15, are shown in the upper panel of Fig. 2. For this figure we used the Laplace jump distribution described in later sections (eq.32), which enforces a small drift velocity downwards towards X = 0.For each trajectory the final particle jump is displayed by a dashed line, highlighting the ability of these models to capture large fluctuation dynamics.
In the lower panel of Figure 2 we highlight the ability of random walk models to capture the diffusive dynamics of systems of many particles.We plot p(X|k) the probability density of finding a hopping particle at a given X coordinate after k steps for an ensemble of N = 10 4 hopping particles, which all start at X = 30 and evolve with the same Laplacian jump distribution as the upper figure.This figure highlights the well known result that in the limit ∆t → 0, random walks asymptote to Brownian motion, and are described by a diffusion equation.For finite ∆t the average evolution of many particles still captures many of the properties of diffusive dynamics.This is an important property of these systems, as we know that the average time-dependent evolution of accreting systems is diffusive (at radii large compared to the ISCO), and therefore random walk models will capture the gross large-radii properties of accreting systems, with the addition of better modelling large amplitude fluctuations in the inner disc.
As we have seen, for a given choice of f (v), and a choice of boundary conditions, a system of many random walk particles can be fully described by some particle density function p(x, t) which determines the number density of hopping particles in the region x → x + dx at time t.Keeping in mind the astrophysical setting we are examining, in this notation the coordinate x should be thought of as the distance from the ISCO (i.e., x = 0 corresponds to the ISCO).
We wish to answer the question "what is the velocity distribution of a particle given that it is crossing a boundary at x = x ′ ".We denote this quantity f (v | x ′ ).This can be rather generally written (the notation p(A|B) denotes the probability of A given B, and we suppress any explicit time dependence) where we split the contributions from positive and negative velocities via the Heaviside theta function Θ(z < 0) = 0, Θ(z > 0) = 1.In deriving this expression we have made the assumption that each of the particles evolves independently of all other particles.In this expression the term on the left of the numerator describes all trajectories moving outwards in the disc (x ′ > x), while the term on the right describes all trajectories moving inwards (x ′ < x).We define the crossing velocity through a surface at radius x = x ′ as Our definition of the crossing velocity is closely related to the "leapover distance", which is a well studied quantity in the random walk literature (e.g., Koren et al. 2007).
If we have an absorbing boundary at x ′ = 0, for example an ISCO whereafter a fluid element is extremely unlikely to fluctuate back into the stable part of the disc, then p(x < 0) = 01 , and and therefore the boundary crossing velocity is It is important to note that only negative velocities contribute to this result (only negative velocities are able to cross the boundary).In addition, higher speeds contribute to a more significant degree than smaller speeds, as they have an increased available crossing distance.
We will demonstrate in the following section that explicit models (of f (v)) offer insight into the physics we are looking to describe.Finally, it is important to highlight a key assumption we shall be making when solving explicit random walk models.We shall assume that each of the fluid elements in the disc can be treated as undergoing its own random walk, which does not interact with its neighbouring fluid elements.Of course, the dynamics of a real fluid is characterised by complex non-linear interactions, and each fluid element in the disc naturally interacts with its neighbours.The random walk framework is of course an approximation, allowing us to insert "by hand" the turbulent fluctuations into the dynamics in a controllable fashion.To leading order, we expect the interactions in a real fluid to make the jump distribution f also a function of the density of particles in the local vicinity and position f (x, v, p).Such a modification would render the random walk problems unsolvable in a closed form, but fortunately we find that our results are insensitive to the precise functional form of f , providing confidence in this general approach.We reiterate that the random walk framework is intended to probe an argument based on physical insight (see section 2), and is not intended to be a quantitative derivation of fluid properties.

EXPLICIT MODELS
As can be seen in equation 16, the mean crossing velocity at a location x ′ depends on the probability of finding a particle at all locations x within the disc.The density of particles p(x, t) at position x and discrete time t = k∆t evolves according to Formally the velocity integration limits extend to ±∞, but constraints can be placed on f (v) to ensure adherence with causality (or relativity, etc.).In this paper we are interested in steady state models of accreting flows.In particular, when t → ∞, one finds that the steady-state density p(x) satisfies the integral equation (20) where we have defined x ′ ≡ x − v∆t, and used p(x ′ < 0) = 0 to determine the integration limits.This steady state expressions can be interpreted as either the large time behaviour of an initially uniform density of particles p(x, t = 0) = p0 = const, or the large time behaviour of a system which is continuously fed a stream of particles at a large distance from the inner boundary (i.e., each new particle starting in the limit x → ∞).Equivalently, it can be interpreted as the relative fraction of time a hopping particle starting at large radii x → ∞ spends at each radius on its trajectory to x = 0.This integral equation cannot be solved in general, but we now consider two cases where exact solutions can be found.

Simple toy model
Consider the following model for the distribution of velocities with which a particle can move This simple model represents a flow with small mean −ϵ and large fluctuations V which can occur with random (equally probable) direction.This is our simplest first approximation to an accretion flow.Substituting the above jump distribution into the governing integral equation leaves This equation has analytical solutions split across two radial regimes, the first where x < β: with solution p1(x) = b1 exp(x/D1), and the second where with solution p2(x) = 1 − exp(−x/D2).The coefficients are given by the solutions of the following equations (determined from the governing equation and simple matching conditions at x = β) Where the upper implicit equation must be solved numerically (Appendix B).At large distances from the absorbing boundary, the probability of finding a fluctuating particle at a given radius will tend to unity p(x ≫ D2) → 1, therefore whereas at x ′ = 0, we have Both of these results are readily interpretable.Firstly, at large radii, the typical surface crossing velocity is twice the mean drift velocity of the particle.The fact that the crossing velocity is order the mean drift velocity is an entirely expected result.The factor two simply results from the smaller distance on the positive jump side from which particles can cross a given boundary (as the mean and fluctuating velocities are working against each other) compared to the negative jump direction (where they work together).On the other hand, the inner boundary crossing velocity is given simply, and intuitively, by the velocity of the negative jump.In Fig. 3 we display the crossing velocity as a function of location x, showing the clear increase in ⟨vcross⟩ as the absorbing boundary is approached.

An exponential random walk
The delta-function jump distribution considered previously naturally represents an unphysical simplification.In this sub-section we demonstrate that the key results are unchanged by considering a more complex jump distribution.Consider the case of the exponential random walk, with jump velocity distribution given by the Laplace distribution (e.g., Comtet & Majumdar 2005;Majumdar et al. 2006) where ϵ > 0 is the magnitude of the drift (in the direction of lowering x) , √ 2V is the standard deviation of the noise and |z| represents the absolute value of a variable z.This specific jump distribution f (v) has the following property which will be useful to compute the steady state properties.We can make progress by taking two derivatives with respect to x on both sides of Eq. ( 20), yielding Using the property in Eq. ( 33), we find In other words, starting from an integral equation, we have derived a nonlocal differential equation for p(x).It will again be of use to define the two length scales χ ≡ ϵ∆t and β ≡ V ∆t.We then make the ansatz (Majumdar et al. 2012) which upon substitution leads to the following implicit equation for the length scale which has a unique solution D > 0 for every χ, β > 0 (Appendix B).Note that in the limit of small drift velocities ϵ → 0 (the relevant limit for accretion), we have which can be thought of as the length scale at which the flow starts to "learn" about the boundary.Inserting this expression for p(x) back into the integral equation ( 20), we find which gives a condition to determine b.This determination is most easily performed at Each integral is now elementary (see Appendix A), and we have found an exact solution for the steady state probability density function of a discrete time random walk.We plot the solution p(x) for different values of ϵ in Fig. 4 (including both the astrophysically relevant ϵ ≪ V limit, but also larger ϵ).The smaller the drift velocity the earlier the flow "learns" about the absorbing boundary.This can be seen in Fig. 4, where each solution differs from the constant p(x) ≃ 1 at ϵx ≃ 3, meaning that the length scale at which the flow begins to "learn" about the boundary scales as δ ∼ 1/ϵ, as predicted from eq. ( 38).
With the probability density p(x) determined, the crossing velocity may be found straightforwardly from equation ( 16).We display the crossing velocity as a function of x for the exponential jump distribution in Fig. ( 5), for a number of different drift velocities ϵ (including both the astrophysically relevant ϵ ≪ V limit, but also larger ϵ).As expected, far from the inner boundary the crossing velocity is of order ϵ, and independent of the typical velocity fluctuation scale V .However, as the inner boundary is approached, the crossing velocity increases substantially and becomes effectively independent of the drift velocity ϵ at the location of the boundary edge.
The crossing velocity across the absorbing boundary can be computed exactly, and is given by the following expression which is also an elementary integral (presented in full in Appendix A).The relevant result in the accretion context is the ϵ → 0 limit,  where we find We find that the inner boundary crossing velocity is always significantly larger than ϵ when ϵ ≪ V .Figure 6 shows the inner boundary crossing velocity as a function of drift velocity2 , including results from numerical simulations.The potentially counter intuitive result that the boundary crossing velocity is more than twice the typical fluctuation scale can be understood by noting that particles travelling towards the boundary with higher speeds contribute to a more significant degree than those crossing with smaller speeds, as they have an increased available crossing distance and are therefore preferentially selected.

Crossing velocity far away from the boundary
It is important to verify that the crossing velocity as defined in this work has the "proper" behaviour at large radii, i.e., that it tends to a constant value corresponding to the mean drift of the flow, independent of the properties of the (unknown) jump distribution, f (v).
Far away from the absorbing boundary the stationary density p(x) is approximately constant (as there is no length scale in the problem).As a consequence, from Eq. ( 16) we find We perform the change of variable v → v ′ = v + ϵ, where −ϵ is the mean of the jump distribution f (v), here assumed small compared to the fluctuation scale.This yields where the notation ⟨.⟩ v ′ denotes an average over v ′ .Note that ⟨v ′ ⟩ v ′ = 0 by construction.We also assume that the distribution of v is symmetric around v = −ϵ, i.e., that f (v In the limit where the variance of f (v) is much larger than ϵ, the contribution to these integrals will come from velocities |v ′ | ≫ ϵ, and therefore we may expand where we have used the fact that ⟨|v ′ |v ′ ⟩ v ′ = 0 and ⟨sgn(v ′ )⟩ v ′ = 0, as a consequence of the symmetry of the distribution of v ′ .Note that for distributions with ⟨v the result above would not be valid in general.

IMPLICATIONS OF THESE RESULTS IN AN ACCRETION CONTEXT
Both the toy random walk calculations, and the general insight of section 2, suggest that the typical velocity with which a fluid element crosses the ISCO is much better approximated by the turbulent velocity scale than the mean drift velocity of the flow.In this section we highlight the effects this insight has on the thermodynamic disc quantities on either side of the ISCO.

Typical scales of relevant parameters
In this section we estimate the relevant scales of the fluctuation and drift velocities in a standard accretion flow.We perform this analysis in the Newtonian limit, as the idea is only to understand the scales involved, not perform a full rigorous analysis.We follow the notation of section 4, where we denote by ϵ the (small) radial drift velocity of the flow, and V the typical turbulent velocity fluctuation scale.Using standard α-type scaling arguments, we note that the typical velocity fluctuation scale is assumed to be This follows from the definition of the usual Shakura & Sunyaev (1973) alpha prescription and by assuming δv r and δv ϕ have similar magnitudes.The typical drift velocity is (e.g., Pringle 1981) and the fluctuation timescale is of order the orbital timescale (i.e., turbulent fluctuations are excited over the shortest timescale in the problem) In these expressions we have used the approximate (Newtonian) solution of hydrostatic equilibrium to relate the sound speed to the orbital speed In the accretion context we are therefore well into the ϵ/V → 0 limit: for typical (thin) disc parameters.
The toy exponential jump distribution model of section 4 demonstrates that there is a radial scale at which the crossing velocity begins to deviate from the drift velocity, or in effect there is a radial scale at which the flow starts to "learn" of the absorbing boundary.Using the scaling highlighted by eq. ( 38), we find i.e., the flow starts to learn of the ISCO roughly ∼ one ISCO radius away from the ISCO.Within an innermost layer of size (see fig. 5) the crossing velocity is given by its x → 0 asymptotic value of ∼ α 1/2 cs.

Energy conservation
We have argued in this paper that the classical calculation of the radial velocity of an accretion flow must be modified near to the ISCO radius, as the turbulent velocity fluctuations, assumed to vanish in the main body of the disc on average, develop a non-zero directional bias as a result of the absorbing boundary in the flow.To characterise the thermodynamic properties of the flow one begins by solving the constraints of energy conservation.It is interesting to highlight how a non-zero directional bias in δU r , provided it remains smaller than the orbital velocity scales, does not modify the dominant energy balance equation, and therefore the classical (Novikov & Thorne 1973;Page & Thorne 1974) temperature profile of the disc is unchanged.The energy balance in the main body of the disc can be determined from the conservation of the stress-energy tensor of the disc T µν .We present a full relativistic calculation of this energy balance in Appendix C, while here simply quoting the key result, namely: where ρ is the density of the disc, W µν is the turbulent stress tensor, and q z is the heat radiated out of the disc surfaces by photons.We see that there is a term proportional to U r + δU r in this energy balance equation, which describes the effects of advection, and so at first it may appear that the work in this paper modifies the temperature profiles of standard theory.However, the prefactor of the advection term is itself of order the radial velocity scale as it vanishes identically in the Kerr midplane for circular orbits (see Appendix 3 of Mummery & Balbus 2019, for a formal proof).
As such, modifications to U r + δU r which remain sub-orbital U r + δU r ≪ rU ϕ do not substantially modify the energy balance constraint, as the advection term enters at order O U r /rU ϕ .The above identity will cease to be true for extreme values of U r , when the energy balance becomes advection dominated (and a "slim disc" regime is entered Abramowicz et al. 1988).However, we shall show that the typical α 1/2 cs values of thin discs at moderate accretion rates are of order ∼ 10 −3 c, while the rotational velocities at the ISCO are of order rU ϕ ∼ c, meaning this reasoning is robust.Within the ISCO, of course, this argument completely breaks down and a new formalism must be employed (Mummery & Balbus 2023).
The final steps of the derivation of the radiative temperature of these disc solutions are therefore unchanged from the classic (Novikov & Thorne 1973;Page & Thorne 1974) calculation, and we do not repeat them here.The final result, in the steady state, is that the radiative temperature TR depends on the remaining parameters through where M•, a•, Ṁ are the black hole mass, dimensionless spin parameter and disc accretion rate respectively, η is the square root of the radius normalised to the gravitational radius η ≡ rc 2 /GM•, and R(η) is the relativistic correction function due to Page & Thorne (1974) where and The parameter δJ expresses the inner boundary condition of the ISCO stress, where δJ = 0 represents a vanishing ISCO stress.
The parameter δJ corresponds physically to the fraction of its ISCO angular momentum an accreting fluid element is able to pass back to the main body of the disc over its plunge.

Solving for the thermodynamic profiles
As we have just demonstrated, upon specifying the free parameters of the disc theory ( Ṁ , M•, a•, α, δJ ) we have a predetermined radiative temperature profile TR(r).We now derive the full solutions of the thermodynamic disc properties with a modified radial velocity profile.
From mass conservation the radial velocity is related to the surface density through Vertical hydrostatic equilibrium gives the scale height of the disc (Abramowicz et al. 1997) which relates the discs surface density Σ to the disc density The pressure of the disc is given by the sum of the gas and radiation pressures Finally, the central and radiative temperatures are related through the optical depth and we shall assume that electron scattering opacity dominates within the flow κ ≃ κes.All of the above expressions are completely standard.The new addition, using the insight gained from the random walk calculations, is that the trans-ISCO velocity of the flow will be where cs,I is the ISCO speed of sound.It turns out that this is sufficient to close the full set of disc equations at the ISCO.Remembering that TR(rI ) ≡ TR,I is known, the above equations can be manipulated into an algebraic equation for cs,I in terms of TR,I , Ṁ , M, a and α.Explicitly, the ISCO speed of sound satisfies (see Appendix D for a derivation) (66) which is trivial to solve numerically.With this boundary condition determined, the remaining thermodynamic disc profiles can be computed once a specification of the disc's radial velocity is given.

Radial velocity profile
We define a new velocity profile of the accretion flow solutions where U r NPT denotes the classical solution (Novikov & Thorne 1973;Page & Thorne 1974), we define U r < 0 for accretion (U r NPT is negative), and I(r) is an interpolation function, satisfying the constraints The rational here of course is that the drift velocity is the asymptotically correct value for large radii, but that the non-zero directional bias in the fluctuations modifies the disc flow near to the ISCO.There are of course any number of functions which satisfy these asymptotic constraints.There is an additional requirement on the interpolation function however, which helps narrow down the choice somewhat.We require that I(r) transitions from zero to one over a well defined scale length (which in the disc will be set by the turbulent eddy scale ∼ H).We therefore turn to exponential functions, which have controllable length scales.
We experimented numerically with simple exponential decays All of which allow the radial scale over which the flow deviates from the classical description to be controlled through ∆r (the exponent 4 here is chosen so that I(r) goes to zero sufficiently quickly so as to not modify the disc thermodynamics on large scales).We found no real qualitative difference between different choices of interpolation functions which shared a common tuneable length scale ∆r.This is not particularly surprising, as each disc model has the same behaviour at the ISCO (as they share the same radial velocity), and the same behaviour at large radii (given by the Novikov & Thorne 1973, disc solution).As such, we simply display solutions with an interpolation function given by eq. ( 69), while noting that the properties of I(r) will be best constrained through future comparison to numerical GRMHD experimentation.
With I(r) specified, we can solve fully for the disc thermodynamic properties.One specifies the disc and black hole parameters Ṁ , M•, a•, α and the ISCO stress δJ , which specifies the radiative temperature profile TR(r).One then solves the boundary condition constraint for cs,I , which gives the radial velocity profile through the above parameterisation.Mass conservation then gives the disc surface density, from which the central temperature, pressure, density and scale height of the disc can be calculated.
Within the ISCO we employ the formalism of Mummery & Balbus (2023), using the values of the various thermodynamic quantities at the ISCO as a boundary condition.

Example solutions
In this section we display the properties of some example solutions of the disc equations in this new framework.We take a moderately high ISCO stress (although we choose a value in the middle of the range found from GRMHD simulations), and choose other parameters suitable for a comparison to a typical X-ray binary.Explicitly we take M• = 10M⊙, a• = 0, α = 0.1, and Ṁ = 0.1 Ṁedd .We take a dimensionless ISCO stress parameter of δJ = 0.03.Thermodynamic disc properties for this solution, and a comparison to conventional models, are presented in Fig. 7.
We display with dashed curves the classical finite ISCO stress solutions, which display pronounced cusps at the ISCO (cf.Fig. 1), and by dotted curves the corresponding vanishing ISCO stress solution, which show an unphysical radial velocity dispersion U r → ∞ at the ISCO.Finally, by solid curves we display the new solutions de- rived in this work, with a trans-ISCO velocity set by α 1/2 cs, we see that the cuspy nature of the finite ISCO stress solutions has been removed, and the transition across the ISCO is significantly smoother.Each panel displays a different thermodynamic quantity, and the vertical axes of these plots are normalised by the value of the thermodynamic quantity at 10rg, except for the radial velocity (lowest panel) which is plotted in units of the speed of light.The new formalism put forward in this work nicely bisects the two traditional approaches.
In Figure 8 we plot various different thermodynamic properties of a disc evolving about a Schwarzschild black hole with differing values of the dimensionless ISCO stress parameter δJ .The ISCO radius is denoted by the vertical black dashed line.For reference, a vanishing ISCO stress solution is displayed by grey dots.The values of δJ used are log 10 δJ = −5, −4, −3, −2.3, −2, −1.3, −1, where higher ISCO stresses can be identified by larger ISCO values of the radiative temperatures (centre left panel).The other parame- ters used in constructing this solution are typical for galactic X-ray binaries M• = 10M⊙, Ṁ = 0.1 Ṁedd , α = 0.1.Interestingly, in this new formalism, certain disc quantities show a much reduced dependence on the ISCO stress.This is particularly true for the disc surface density Σ, which shows barely any dependence on δJ despite it being varied by 4 orders of magnitude.This is because the surface density is set entirely by the radial velocity (through mass conservation Ṁ = 2πrΣU r ), and the trans-ISCO fluctuation velocity turns out to be only weakly dependent on the ISCO stress in this new formalism.
Other disc parameters remain more sensitively dependent on the local physics of the ISCO stress.This is most notable for the radiative temperature of the flow TR, which shows a dependence on ISCO stress within r ≲ 10rg, and is extremely sensitive to the ISCO stress at radii within the ISCO.Similarly, the increased temperatures of these solutions leads to greater pressure support and notably different scale heights of each solution.This filters through to a much reduced intra-ISCO density (ρ) for larger ISCO stresses.This is of potential observational interest, as the disc density ρ determines the ionisation fraction (ξ) of the flow if the flow is illuminated by an incident X-ray flux FX , ξ ∝ FX /ρ.Lower densities from larger ISCO stresses will filter through to higher ionisation fractions, and correspondingly reduced iron line fluorescence.
It is clear from Fig. 8 that some disc quantities still display a slight cusp at the ISCO, even within this new framework.This will remain an unavoidable effect of analytical models of trans-ISCO flows which involve the piecewise joining of intra-and extra-ISCO flows.The present work minimises the presence of these kinks to as much of a degree as possible, and we do not believe that any remaining cusps will dramatically influence inferences from the fitting of observational data.

COMPARISON TO GRMHD SIMULATIONS
In this section we take our extended global thin disc solutions and compare their thermodynamic profiles to those extracted from GRMHD simulations.A full comparison to dedicated numerical simulations is postponed to a future work, and we for now concentrate on potentially observable profiles extracted from previously published experiments.
The two main potentially observable properties of a thin disc are the radial dependence of the "effective" (radiative) temperature profile, and the density of the flow.The radiative temperature is the key parameter of interest for so-called continuum fitting modelling of Galactic X-ray binaries (see e.g., McClintock et al. 2014), as it directly determines the locally liberated flux in the fluids rest frame.As we discussed earlier, the density profile of a flow sets the ionisation fraction of the material when illuminated by an external X-ray flux FX , with ionisation fraction ξ ∝ FX /ρ, and is therefore of direct interest to iron line studies (see e.g., Reynolds 2013).
We first compare the radiative temperature profile of our model to those published in Zhu et al. (2012).Zhu et al. (2012) extracted a radiative temperature from the local cooling rate computed in the GRMHD simulations run by Penna et al. (2010), and used them to examine some effects of the (neglect of the) plunging region on continuum fitting spin measurements.The radiative temperature profiles of Zhu et al. (2012) are displayed in Figure 9 by black dots, for three different black hole spins a• = 0 (lower panel), a• = 0.7 (middle panel) and a• = 0.9 (upper panel).Zhu et al. (2012) model an M• = 10M⊙ black hole accreting at roughly Ṁ ∼ 0.1 Ṁedd .We take these parameters as input to our analytical model.Various "effective" α parameter was reported by Zhu et al. (2012) and we take their values (α = 0.1 for a• = 0.7 and a• = 0.9, α = 0.01 for a• = 0) for simplicity (Zhu et al. 2012, their figure 6).We then only have the ISCO parameters to fit, namely δJ .
We overplot in Figure 9 the vanishing ISCO stress radiative temperature curve (red dashed curves), which are forced to zero at the ISCO contrary to the simulation results, and in blue (solid curves) the model developed in this paper.We determine the appropriate value of the ISCO stress by minimising the loss function L = i (Tsim(ri) − T disc (ri; ∆r, δJ )) 2 /Tsim(ri) 2 , where Tsim and T disc are the simulation and analytical effective temperatures respectively.The ISCO stress parameters are δJ = 0.0055 for a• = 0.9, δJ = 0.0085 for a• = 0.7 and δJ = 0.007 for a• = 0. We found no sensitivity to the fluctuation length scale parameter ∆r, which we set to equal to the scale height of the disc in the solutions.We see that we recover the global properties of the GRMHD simulations rather well (note the inflection points in the radiation temperature profiles around the ISCO).This is an important result and motivates future development of extended continuum fitting models which include intra-ISCO emission.
Simulating the density profiles of thin GRMHD accretion flows with radiative transport effects included has only recently become computationally feasible (e.g., Liska et al. 2022;White et al. 2023).We extract the density profile from a SANE (i.e., "standard and normal evolution"; the low magnetic field limit relevant for comparing to thin discs) simulation run by Liska et al. (2022), which was run for a spin a• = 0.9375 black hole, with mass M• = 10M⊙ (Liska et al. 2022, refer to this simulation as RADTOR in their paper).The accretion rate in this simulation was set to be Ṁ ∼ 0.35 Ṁedd , and Liska et al. (2022) extract an "effective" α parameter which depended on radius but was roughly equal to α ∼ 0.03 a value we take in this work.The density profile extracted from Liska et al. ( 2022) is plotted in Figure 10.
Again, as the physical parameters of the main body of the disc are prespecified, we only fit the intra-ISCO parameters of our model to the data.The fit is performed by minimising the loss function L = i (ρsim(ri) − ρ disc (ri; ∆r, δJ )) 2 /ρsim(ri) 2 , where ρsim and ρ disc are the simulation and analytical density respectively.This loss function is dependent on the parameters ∆r and δJ only.The fit was performed only for r ≤ 50rg.We find that the parameters δJ = 0.01, ∆r/rI = 0.03 fits the global density profile well.The deviation at large radii (r ≳ 75rg) is likely a result of the finite mass reservoir in the GRMHD simulation.
It is too early to say whether the fact that δJ ∼ 10 −2 best describes two different simulations which measure different disc quantities represents an interesting result, or is simply coincidence.
As we move into the era where radiative GRMHD simulations of thin black hole discs become computationally feasible, we expect to perform many future tests of the models developed here.

DISCUSSION
In this paper we have identified a potential issue with classical modelling of black hole accretion flows -by neglecting the finite turbulent eddy size of discs, classical models do not capture a directional bias in the turbulent velocity field induced by the ISCO, which acts much like a one-way gate for fluid elements.Ultimately, this bias has the dynamical effect that accretion flows typically cross the ISCO at significantly enhanced velocities (compared to the predictions of classical models).
In aiming to elucidate this effect in a controlled manner, we have examined the properties of a set of random walk models.Random walk models are flexible, in that they allow both finite perturbations and global diffusive evolution to be examined simultaneously.They remain, of course, a simplification of the true dynamics of a turbulent MHD flow.In particular, the approximation that each fluid element behaves independently of all the others is an over simplification of the complex non-linear physics of a fluid.As these non-linear effects most likely mean that any one choice for a random walk velocity distribution (e.g., eq.32) will be inaccurate, it is important that we have found that the qualitative results of the random walk models are insensitive to this assumption.We also stress that the gross behaviour we are discussing can be intuited on purely physical grounds (see section 2).
As the velocity profile of a black hole accretion flow will deviate from the classical α-model predictions at short length scales from the ISCO (e.g.Novikov & Thorne 1973), and the trans-ISCO velocity will be of order ∼ α 1/2 cs, we have introduced an interpolation function I(r) into our thin disc model radial velocity.This interpolation function must satisfy I(r → ∞) → 0, and I(r → rI ) → 1, with the transition expected to occur over a length scale of order one ISCO radius out from the ISCO, but is otherwise not currently well understood.Considering a class of exponential interpolation functions with tuneable transition scales, we did not find much sensitivity of the disc thermodynamic properties to the precise functional form chosen.
Fortunately, when confronted with numerical data (e.g., Figure 10), the tuneable parameters in the interpolation functions can be constrained (and the new model fits well).As we move into the era where radiative GRMHD simulations of thin black hole discs become computationally feasible, we anticipate that the properties of the radial velocity transition will be well constrained by simulations, which are better suited to this task than analytical derivation.
The fact that these extended models reproduce properties of full GRMHD simulations (Figs. 9, 10) is extremely promising.We stress that the modifications and extensions presented in this work are exclusive to the inner disc regime, and we expect classical α-modelling to remain a good description of accretion flows at large radii.This can most clearly be seen in the reproduction of the large radius GRMHD simulation density profile in Figure 10.The good model fit extends well into the classical regime, where the modifications put forward in this paper have no effect.
Where classical α-modelling is clearly breaking down however is at small radial scales, typically at ∼ 1 − 2 ISCO radii from the black hole.This can most clearly be seen in the cusps produced by classical models when ISCO stress parameters typically found in numerical simulations are substituted (Fig. 1, 7); these cusps are never themselves reproduced in simulations.While the new model derived in this paper smooths out these cusps to a significant degree (e.g., Figs. 7,8), there are still some small scale cusps which remain at precisely the ISCO.This is an artifact which results from the piecewise joining of two distinct accretion solutions together, and we do not expect this behaviour to be physical.It seems likely to us that in a real fluid these cusps would be smoothed out, as is seen in GRMHD simulations (Fig. 9).
In addition, the classical "zero stress" boundary condition imposed on α-models clearly does not reproduce the properties of simulated accretion flows at small radii (Fig. 9).The ISCO stress parameters we infer from fitting to GRMHD data are at the scale δJ ∼ 0.01, which is at the smaller end of the range estimated from previous simulations (for example Noble et al. 2010, found δJ ∼ 0.1, while Penna et al. 2010 found δJ ∼ 0.02).This is, however, certainly sufficiently large to show substantial deviations from disc profiles which assume zero stress (Fig. 9).
The stress at the ISCO is fundamentally magnetic in origin (e.g., Gammie 1999;Agol & Krolik 2000), and likely increases sharply with magnetisation (e.g., the Gammie 1999, model of the ISCO stress).As these new solutions reproduce properties of GRMHD simulations at a fraction of the computational cost, it is likely that we are heading into a future where observational constraints (like those inferred from continuum fitting and iron line modelling) will be able to place constraints on δJ , the value of which has been a long controversial theoretical question.

CONCLUSIONS
In this paper we have examined the near-ISCO behaviour of thin black hole accretion flows, with a particular focus on the appropriateness of the classical treatment of the mean fluid flow in this limit.We have argued that because a turbulent flow (like MRI driven accretion) will display macroscopic (of order the disc scale height) perturbations in its velocity field, the classical description (which does not distinguish scales above or below the turbulent eddy scale) becomes an increasingly poor model as absorbing boundaries, such as the ISCO, are approached.Physically, this argument stems from the fact that in a turbulent flow large velocity fluctuations can carry a fluid element over the ISCO from a finite distance away, from which it will not return, a process without analogy in a classical disc model.This introduces a non-zero directional bias into the velocity fluctuations in the near-ISCO disc, a property which is ordinarily ignored.
To examine the effects of this non-zero directional bias in the velocity fluctuations of accreting flows, we have examined the properties of some random walk models.Random walk models are more mathematically flexible than purely viscous systems, and are well suited to modelling physical systems with both a global diffusive character (like an accretion flow at large radii), but also with large amplitude velocity fluctuations (relevant for an accretion flow at small radii).It transpires that in the astrophysically relevant limit where the mean drift velocity of a fluid element ϵ is much smaller than the typical turbulent velocity fluctuations V (in thin discs this ratio is of order ϵ/V ∼ α 1/2 (H/R) ≪ 1) the typical velocity with which a fluid element crosses an absorbing boundary is of order the fluctuation scale ∼ O(V ) ≫ ϵ.
This increased radial velocity modifies the local thermodynamic quantities of this disc on either side of the ISCO.A practical appli-cation of this work is that it removes (as far as possible) cusps at the ISCO which are present in previous models of finite ISCO stress discs (e.g.,Figs. 1,7).This will be of practical importance when it comes to fitting analytical intra-ISCO models to data (e.g., Reynolds & Begelman 1997;Zhu et al. 2012;Wilkins et al. 2020) as such discontinuities might drive the overall fitting procedure.
In this framework thin disc accretion around black holes is comprised of three fundamental regimes.At large radii r/rI ≫ 1 accretion is diffusive (or "viscous") dominated, and standard models work well.At an order unity distance from the ISCO r/rI ∼ 2 the flow transitions to a fluctuation dominated state, and the fluid begins to "learn" about the ISCO.Within the ISCO r/rI < 1 the flow transitions to near geodesic motion, and is gravitationally dominated.The analytical models we describe in this paper smoothly transition between the three regions.
Finally, we have demonstrate that these new models are in good accord with the outputs of GRMHD simulations of thin discs.In Figure 9 we demonstrate that the radiative temperature of these new solutions is in good accord with the results of Zhu et al. (2012), who computed the locally liberated flux from the GRMHD simulations of Penna et al. (2010).This local radiative temperature is the chief physical parameter which determines the thermal X-ray emission observed in X-ray binary soft states, and can therefore be directly probed with observations (e.g.McClintock et al. 2014).
Similarly, we reproduce (Fig. 10) the density profile of the thin disc weak magnetic field simulation of Liska et al. (2022).Density profiles of discs can in principle be probed by the iron line fitting technique (Reynolds 2013), to which the plunging region provides a non-negligible contribution (Reynolds & Begelman 1997;Wilkins et al. 2020).It is our intention to use the models developed here and in Mummery & Balbus (2023) to extend "continuum fitting" and iron-line analysis procedures, to include emission and material inside the ISCO.
For metrics which do not depend explicitly on coordinate x γ (the Kerr metric does not depend on t and ϕ) ∂γgµν = 0, and so Expanding, and using ∇µg µν ≡ 0, and assuming that ρc 2 ≫ P + e, we are left with √ gρ (U µ U0 + W µ 0 ) + √ g(q µ U0 + U µ q0) = 0, (C8) which upon expanding (and neglecting the asymptotically small q0 term) is These are, in order, the definition of the speed of sound, the solution of vertical hydrostatic equilibrium, the definition of the disc density and pressure, the approximate solution of radiative transfer in the disc atmosphere, and the conservation of mass in the disc.Note that we have used the exact result U 2 ϕ,I − a 2 c 2 (1 − U 2 0,I ) = 2GM•rI in the equation of hydrostatic equilibrium (proof in Mummery & Balbus 2023).
Substituting D5 into D4, before substituting D4 and D3 into D1 leaves to a simple expression for cs,I in terms of known quantities (α, Ṁ , rI ) and the quantities HI and ΣI .This final expression can then be expressed entirely in terms of cs,I and known quantities by the substitution of D2 and D6.The final result is as shown in the main body of the paper, namely (D7)

Figure 1 .
Figure1.An example of "cuspy" behaviour in classical finite ISCO stress models of black hole discs (parameters: M• = 10M ⊙ , a• = 0, Ṁ = 0.1 Ṁedd , α = 0.1, δ J = 0.1; see section 5 for a description of the physical meaning of these free parameters).In the upper panel we display various normalised (by their ISCO values) disc quantities on both sides of the ISCO (ISCO displayed by vertical black line).The evolution across the ISCO is not smooth in these models, a result of the very low trans-ISCO radial velocity of the disc (lower panel).We zoom into the inner 2rg either side of the ISCO (black dashed curve, r I = 6rg).

Figure 2 .
Figure2.Upper panel: an example of individual random walk particle trajectories (solid lines), where we plot the position of the particle X k after k jumps.The dashed line shows the final jump of the particle, and the dot shows the position of the particle after this last jump.The arrow shows the direction of the imposed drift (the average of the jump distribution, ϵ), ensuring that all particles eventually cross X = 0. Clearly random walk models can incorporate the physics of large amplitude velocity fluctuations.Lower: The evolving (numerical, averaged over 10 4 trajectories) probability density function of a group of evolving random walk particles at different time-steps k.When averaged over large numbers of particles random walk calculations reproduce properties of diffusive evolution (in the limit of large k, with X 0 /

Figure 3 .
Figure 3.The absolute value of the crossing velocity as a function of x.The two horizontal black lines are the asymptotic results | ⟨vcross(x = 0)⟩ | = V + ϵ and | ⟨vcross(x ≫ 0)⟩ | = 2ϵ.For this figure we take V = 1, ϵ = 0.02, and time step ∆t = 0.1.The discontinuous gradient of ⟨vcross⟩ near x ≃ 0.1 is due to the discrete nature of the jump distribution f (v).The absorbing boundary is placed at x = 0.

Figure 4 .Figure 5 .
Figure 4.The probability density function p(x) as a function of ϵx for different values of the drift velocity ϵ.For this Figure the fluctuation velocity scale was normalised to V = 1, and the time step was ∆t = 0.1.The absorbing boundary is placed at x = 0.

Boundary crossing velocity v b Figure 6 .
Figure 6.The inner boundary crossing velocity, ⟨v b ⟩, as a function of the mean drift velocity ϵ.For this Figure the fluctuation velocity scale was normalised to V = 1.The black dot-dashed curve is the analytical result derived in this paper, whereas the red crosses are numerical averages of N = 10 5 individual particle trajectories.In the small drift velocity regime relevant for accretion the boundary crossing velocity is always ⟨v b ⟩ ≫ ϵ.

Figure 7 .
Figure7.Various disc thermodynamic quantities, normalised by their values at 10rg, for differing models of the inner disc flow.We display with dashed curves the classical finite ISCO stress solutions, which display pronounced cusps at the ISCO, and by dotted curves the corresponding vanishing ISCO stress solution, which show an unphysical radial velocity dispersion U r → ∞ at the ISCO.Finally, by solid curves we display the new solutions derived in this work, with a trans-ISCO velocity set by α 1/2 cs, we see that the cuspy nature of the finite ISCO stress solutions has been removed, and the transition across the ISCO is significantly smoother.In the lowest panel we show the radial velocities of the three models, normalised by the speed of light.The ISCO radius is denoted by a vertical black dot-dashed line.

Figure 8 .
Figure 8. Various different thermodynamic properties of a disc evolving about a Schwarzschild black hole with differing values of the dimensionless ISCO stress parameter δ J .The ISCO radius is denoted by the vertical black dashed line.For reference, a vanishing ISCO stress solution is displayed by grey dots.The values of δ J used are log 10 δ J = −5, −4, −3, −2.3, −2, −1.3, −1, where higher ISCO stresses can be identified by larger ISCO values of the radiative temperatures.The other parameters used in constructing this solution are M• = 10M ⊙ , Ṁ = 0.1 Ṁedd , α = 0.1.

Figure 9 .
Figure9.The radiative temperature profiles ofZhu et al. 2012 (black points; see text), compared to the models developed in this paper (blue solid curves) for three values of the black hole spin (displayed on plot).Also shown are vanishing ISCO stress models for comparison (red dashed curves).The simple analytical models developed in this paper reproduce the gross behaviour of full GRMHD profiles.

Figure 10 .
Figure10.The density profile of the "RADTOR" simulation presented inLiska et al. 2022 (purple points), overplotted with the analytical model developed in this paper (blue curve).We find good agreement across the ISCO plunge, and also at large radii.
∂αgµγ) S µα ≡ 0, (C6) vanishes for any symmetric tensor S µα , since the metric derivative are anti-symmetric in µ and α, while S µα is symmetric in these indices.As a result of this identity, the conservation of disc angular momentum U0∂zq z = 0. (C9)The first term in square brackets in the above expression is just mass conservation within the disc, and is zero.Thus, energy conservation leads to∇µ (T µ 0 ) = ρU r ∂rU0 + 1 √ g ∂µ ( √ gρW µ 0 ) + U0∂zq z = 0. (C10)Identical reasoning as the above produces a symmetric (with 0 replaced by ϕ) equation of conservation of angular momentum∇µ T µ ϕ = ρU r ∂rU ϕ + 1 √ g ∂µ √ gρW µ ϕ +U ϕ ∂zq z = 0. (C11)To derive the energy equation of the flow, take U 0 times the energy conservation equation and add it to U ϕ times the angular momentum conservation equation.This procedure leavesρU r U 0 ∂rU0 + U ϕ ∂rU ϕ + U 0 = −(U 0 U0 + U ϕ U ϕ )∂zq z , (C12)which is the result used in the paper.APPENDIX D: BOUNDARY CONDITION FOR CS,ITo derive the governing boundary condition expression for the ISCO speed of sound, one must solve the coupled equations