# MRI turbulence and thermal instability in accretion disks

###### Abstract

A long-standing puzzle in the study of black-hole accretion concerns the presence or not of thermal instability. Classical theory predicts the encircling accretion disk is unstable, as do self-consistent MHD simulations of the flow. Yet observations of strongly accreting sources generally fail to exhibit cyclic or unstable dynamics on the expected timescales. This paper checks whether turbulent fluctuations impede thermal instability. It also asks if it makes sense to conduct linear stability analyses on a turbulent background. These issues are explored with a set of MRI simulations in thermally unstable local boxes in combination with stochastic equations that approximate the disk energetics. These models show that the disk’s thermal behaviour deviates significantly from laminar theory, though ultimately a thermal runaway does occur. We find that the disk temperature evolves as a biased random walk, rather than increasing exponentially, and thus generates a broad spread of outcomes, with instability often delayed for several thermal timescales. We construct a statistical theory that describes some of this behaviour, emphasising the importance of the ‘escape time’ and its associated probability distribution. In conclusion, turbulent fluctuations on there own cannot stabilise a disk, but they can weaken and delay thermal instability.

###### keywords:

accretion, accretion disks — instabilities – MHD — turbulence – X-rays: binaries## 1 Introduction

The assumption that the turbulent stress is proportional to pressure in accretion disks (the -model) is a fundamental ingredient of classic accretion disk theory (Shakura & Sunyaev 1973). It has been partially justified by successful application to quasi-steady systems, such as the thermal spectrum of dwarf novae and X-ray binaries (e.g. Warner 1995, Gierliński & Done 2004); but when applied to more delicate time-dependent dynamics, such as instabilities, the model has encountered difficulties. For instance, the modelling of dwarf novae outbursts requires different alphas for the high and low states (Smak 1984). Another example involves radiation-pressure dominated accretion flows which the alpha model predicts are subject to thermal and viscous instability (Shakura & Sunyaev 1976, Lightman & Eardley 1974). X-ray observations, however, fail to find variability on the timescales expected (Gierliński & Done 2004), with only the exceptional luminous source GRS 1915+105 and the intermediate black hole HLX-1 exhibiting anything like cyclic behaviour driven by thermal instability (Belloni et al. 1997, Done et al. 2004, Sun et al. 2016, Wu et al. 2016). While it is possible the disks are stabilised by an additional but unknown cooling mechanism, it may be that the heating depends on temperature in a weaker way than the alpha model assumes. For example, the turbulent stresses may be proportional to gas pressure rather than total pressure, or not on pressure at all (Gierlinski & Done 2004).

Of course, one can bypass the alpha model (and its assumptions) and simulate the MHD turbulence in these disks self-consistently, the turbulence then supplied directly by the magnetorotational instability (MRI, Balbus & Hawley 1991). And in fact early work indicated that radiation-pressure dominated flows were thermally stable (Hirose et al. 2009), in agreement with most observations. However, recent local and global simulations do exhibit thermal runaways (Jiang et al. 2013, Sadowski 2016, Mishra et al. 2016), though these expose additional complications that may suppress instability, such as numerical effects (especially box size) and the impact of a mean magnetic flux. Obviously, both observations and simulations indicate that the onset and development of thermal instability is far less straightforward than predicted by the classical laminar theory.

One very clear complication is the fact that the stress-pressure relationship can change, for both numerical and physical reasons, yielding instability or stability depending on conditions in the disk and in the simulation. A quite separate issue is the assumption that thermal instability can be appropriately defined at all, at least when dealing with quasi-steady turbulent states. The alpha-theory treats the turbulence as a static eddy viscosity, and hence the equilibrium state as laminar. However, if the state hosts vigorous fluctuations it may not be well-defined, or even make sense, to add a small linear perturbation on top the stochastic background field and subsequently calculate a growth rate. One envisages that, at the very least, non-model, non-exponential behaviour ensues. Indeed, Jiang et al. (2013) report delayed runaway and algebraic growth rather than exponential growth in their MHD simulations, while Janiuk & Mishra (2012) show via a stochastic 1D model that fluctuations induce random luminosity variations rather than the regular outbursts expected. It is to this aspect of the problem that this paper is devoted, focussing on the constructive and destructive interference of turbulence on thermal instability.

In order to isolate the essence of the problem we employ an idealised model of MRI turbulence and of thermal instability. Unstratified shearing box simulations are performed using the code RAMSES on a state that is MRI turbulent and thermally unstable (at least according to the laminar alpha theory). Note that radiation pressure is omitted and the gas cools due to a simple cooling function. We find that the turbulent fluctuations induce thermal behaviour substantially different to that expected from laminar theory. In particular, the evolution of the temperature resembles more a biased random walk than an exponential runaway, with a wide range of trajectories possible: the temperature in some simulations departs from the laminar equilibrium relatively rapidly, whereas in others it can ‘hang around’ for several thermal times.

This motivates a probabilistic interpretation of instability, and we develop a simple statistical framework based on the model of geometric Brownian motion. A key idea is that of the ‘escape time’ (which replaces the e-folding time). It describes how long it takes for the system to deviate significantly from the equilibrium. Reduced models involving both white noise and the power spectrum of the MRI show that the probability distribution of possesses a long tail. Thus there is a reasonable chance in any given simulation that thermal runaway is delayed. It should be stressed that ultimately realistic models of disks still undergo thermal runaways: turbulent fluctuations can impede instability but it cannot destroy it. The stabilisation witnessed in Janiuk & Misra (2012) we attribute to the peculiarities of their stochastic model and a very large noise amplitude.

Another feature of our MHD simulations is thermal fragmentation when the cooling rate is too small, and hence the laminar thermal instability timescale too short. The disk can then break up into hot and cold clouds. This occurs when the thermal mixing (by turbulence or radiative diffusion) is inefficient compared to thermal instability. On a sufficiently long lengthscale this is always presumably the case, but how this relates to the onset of instability in hot accretion flows is unclear. Estimates of both radiative diffusion and turbulent mixing suitable for the inner regions of X-ray binaries indicate that fragmentation is a marginal possibility.

The paper is organised as follows. In Section 2 we discuss a number of issues pertinent to thermal instability, stochastic fluctuations, and the limits of the alpha theory (not all of which we take up in the paper). The third and fourth sections contain our numerical MHD model and the corresponding results, respectively. We then explore stochastic models in Section 5 and construct a statistical theory to help explain the MHD simulations. Our conclusions are then presented in Section 6.

## 2 Theoretical issues

### 2.1 Stress-pressure relationship

The exact dependence of the turbulent stress, , on pressure is key to the onset of thermal stability in radiation pressure dominated flows. The instability occurs when , but does not occur when (Piran 1978). Here and denote gas and radiation pressure, respectively. The failure to observe signatures of thermal instability in most X-ray observations have been attributed to the stress depending on gas pressure alone, or possibly the geometric mean of gas and radiation pressure (Gierlinski & Done 2004). It has also been speculated that only exceptionally luminous flows could lead to a situation where , explaining the outbursts of GRS 1915+105 (though what exactly causes this shift in the stress’s behaviour is unclear) (Gierlinski & Done 2004, but see also King and Ritter 1998). Numerical simulations of the MRI have since complicated this picture, as they indeed exhibit instability (Jiang et al. 2013), and we are left with the task of numerically tracing out the non-straightforward behaviour of in different conditions.

While earlier unstratified shearing box simulations without radiation pressure found only a weak dependence of stress on pressure (Sano et al. 2004), recent work has shown that when the following conditions are satisfied: the computational domain is sufficiently large, explicit dissipation is included, and there is no net magnetic field (Ross et al. 2016, hereafter RLG16). Small boxes restrict the size of the turbulent eddies and prevent them from fully responding to an increase in pressure (be it gas or radiation). This restriction no doubt played a role in the failure of early radiation-pressure dominated simulations to show thermal instability (Hirose et al. 2009): the radial domain was too small, the eddies unnaturally confined, and as a result the stress unable to respond to changes in the total pressure. When larger radial boxes are used, as in recent work, the instability does in fact materialise (Jiang et al. 2013).

Another numerical effect uncovered by RLG16 was a sensitivity to the grid in simulations with no net magnetic field and no explicit diffusion. As in isothermal runs (Fromang & Papaloizou 2007), the stress is proportional to grid size and this leads to a significantly weaker stress-pressure relationship. Local boxes are not the only domains that exhibit the effect; recent vertically stratified simulations also show that the stress depends on the grid length (Ryan et al. 2017). The weaker dependence, both from numerical dissipation and from the box size, leads to artificially more stable systems (as can be shown from dimensional analysis). It is therefore necessary for both of these numerical complications to be considered when simulating thermal instability involving gas pressure. It is likely that global disk simulations of the MRI also suffer from strong numerical effects, though these have yet to be fully explored.

Another intriguing result from RLG16 is that the stress-pressure relationship depends on the existence and strength of any imposed magnetic flux. The stronger a mean toroidal flux, the weaker the relationship. If this behaviour generalises to other field configurations and to the radiation-dominated regime, then one might speculate that highly magnetised flows are less prone to thermal instability. Indeed global simulations of the MRI suggest that strong magnetic fields impede thermal instability (Sadowski 2016), possibly because they weaken the connection between the stress and pressure. This raises the interesting prospect that to assess susceptibility to thermal instability we must also account for the build up and evacuation of large-scale magnetic flux, in addition to the turbulent dynamics (e.g. Guilet & Ogilvie 2012a, 2012b).

### 2.2 Time lags

The alpha theory is a turbulence closure model, supplying a simple eddy viscosity in place of the complicated and chaotic time-dependent dynamics of the flow. On timescales and lengthscales much longer than the characteristic scales of the turbulence, this approximation provides an adequate description, but its performance worsens the shorter the scales of interest. Since thermal instability can possess growth rates of tens of orbits (or less), then the detailed turbulent dynamics could potentially interfere with its onset.

Simulations, as expected, show that on shorter timescales the stress-pressure relationship is more complicated than a mere proportionality. One interesting feature is a time lag of a few orbits between the stress and pressure. Moreover, it is the pressure that follows the stress on shorter times, rather than the other way around. And so the dependency is opposite to that assumed by the alpha theory: bursts in stress are followed by jumps in pressure (Hirose et al. 2009). What is happening here is that the bursts in stress drive fluctuations in the heating rate (once their associated energy has reached the dissipation scale) and hence cause bursts in pressure a short time later.

The effect of this time delay on a stable thermal equilibrium was first explored by Hirose et al. (2009), who argued that if the direction of causation was from the stress to the pressure then thermal instability may not work. The argument fails, however, to acknowledge that the stress-pressure dynamics exhibit different timescales, with longer timescales ( orbits) characterised by a pressure-dependent stress (as in the classical theory), while the short timescales show the lags described above (Latter & Papaloizou 2012, RLG16).

Follow up work by Lin et al. (2011) and Ciesielki et al. (2012) present linear instability analyses of an alpha disk model with a time-delayed alpha. These show that the simulated delay of 1-10 orbits is insufficiently long to stabilise a thermally unstable state. They also reveal potential inconsistencies in such simple models: for example, infinite growth rates are possible for certain time-lags, a physical impossibility. Of course, these models are also incomplete, as they include only the either the short term or long term dependencies, but not both concurrently.

### 2.3 Temporal fluctuations

MRI turbulence has strong variation over a range of timescales, from tens of orbits to a few shear times (Sano et al. 2004, Lesur & Ogilvie 2008). In what sense does a thermal equilibrium exist in such a system? Time varying perturbations are constantly emerging which lead to a shifting balance of heating and cooling that the system continually responds to. On the other hand, if we assume that there is well-defined mean equilibrium, then it is awash in finite amplitude fluctuations. How can one then undertake a linear instability analysis? Is it meaningful to add a tiny perturbation ontop of a sea of finite amplitude perturbations and check if it grows or not? On long length and time scales this might work, but certainly not on shorter scales.

Putting aside the difficulty of interpreting linear stability analyses, a stochastic system exhibits a range of complicated and sometimes unexpected behaviour. A classic example is the destabilisation of fixed points deemed stable by laminar theory. Originally studied in biological population dynamics (e.g. Levins 1969, May 1973), this feature of noisy systems appears in numerous applications, such as atmospheric modelling (e.g. De Swart & Grasman 1987), where the unresolved short time and lengths dynamics are represented by stochastic terms (see Majda et al. 1999, 2003). On the other hand, the influence of stochasticity on an otherwise unstable fixed point has been studied in financial mathematics, where geometric Brownian motion can be used to model volatile stock prices in a rising market. Despite the mean trend of increasing prices, stochasticity can depress the price of some stock dramatically, if not stabilising the fixed point then delaying a runaway in price for some period of time.

The very last example is perhaps the most relevant for our study of thermal instability, as it possesses the key ingredients of (a) an unstable fixed point (according to a deterministic or ‘laminar’ theory) and (b) stochastic fluctuations. The competition between them gives rise to behaviour one might liken to a biased random walk. In between the kicks delivered by the turbulence, the system drifts according to the deterministic unstable dynamics. One can then imagine certain limits: when the characteristic frequency of the turbulence is much greater than the thermal instability growth rate then we may expect an unbiased random walk, and the system will only weakly sense the underlying thermal physics. In the opposite limit, when the growth rate is much greater than the turbulent frequency, the deterministic laminar dynamics should be reproduced. It is in the intermediate regime, explored in this paper, that interesting nontrivial behaviour manifests. There are also other key ratios, such as the size of the kicks relative to the magnitude of the fixed point or the initial condition. If these are too small, then we return to the laminar case. But for intermediate values, as exhibited by our MRI simulations, system trajectories can deviate markedly from both the laminar behaviour and a simple unbiased random walk.

### 2.4 Spatial fluctuations

In the previous subsection we considered only temporal fluctuations on the system variables, implicitly regarding them as ‘box averaged’ or mean quantities. Indeed, the model assumes a homogenised temperature over , the disk scale height. Turbulent heating, however, is spatially inhomogeneous with strong dissipation occurring in current and vorticity sheets and minimal dissipation in the surrounding regions. These spatial fluctuations also complicate the picture of thermal instability, especially when the instability growth rate is large.

Some form of thermal mixing is necessary to homogenise the temperature of the fluid. This can take multiple forms, such as turbulent advection, radiative diffusion, or thermal conduction. When thermal instability is present, the assumption of a uniform temperature is reasonable as long as the instability time scale is longer than the mixing timescale, . For less efficient mixing or stronger instability, there exists a maximal thermal coupling length scale, . Regions separated by more than this only weakly interact thermally during an efolding time. This implies that regions of a disk separated by more than can undergo thermal runaways independently, and the disk fragments into cold and hot clouds. Regions of strong kinetic and magnetic dissipation are likely to heat catastrophically, while those with weak dissipation will cool catastrophically. How relevant this scenario is in realistic disks is unclear, though perhaps marginally possible in X-ray binaries. It is certainly possible in numerical simulations as we show later.

## 3 Numerical tools and setup

### 3.1 Formulation

We wish to explore the essential features of thermal stability and its onset in MRI driven turbulence and so we choose an idealised set-up to isolate it. We adopt the local shearing box model (Goldreich & Lynden-Bell 1965). To prevent complications such as buoyancy, mass loss and disk expansion we consider the unstratified case. With this model, MRI driven turbulence can be obtained when a Keplerian flow profile is assumed (Hawley et al. 1995). As is conventional, are the radial, azimuthal and vertical spatial variables and , , are the corresponding unit vectors. This frame of reference co-rotates with the disk at some radius with angular frequency . The ideal compressible MHD equations are hence

(1) | |||

(2) | |||

(3) | |||

(4) |

where is the mass density, is the velocity, is the gas pressure, is the magnetic field and the internal energy is denoted by . Heating is represented by and cooling by . This set of equations is then closed by relating the internal energy to the pressure by assuming an ideal gas so that

(5) |

where is the adiabatic index, taken to be . The sound speed is then given by and the pressure scale height by .

Ideally we would include viscosity and Ohmic diffusion and so would be given by the sum of the physical dissipative processes; however, to resolve the diffusion length scales requires a higher resolution than is practical for this study. Instead, we rely on numerical dissipation for heating. By solving Equations (1) - (4) in conservative form, the kinetic and magnetic energy dissipated by the grid is converted to internal energy. Energy that is extracted from the background shear is converted to internal energy and ultimately removed via the cooling function .

Finally, for our cooling function, we take a power law of pressure,

(6) |

where and are both constants. Though this choice is mainly for convenience, it might crudely approximate an optically thin medium.

### 3.2 Numerical methods

All of the simulations that we perform are carried out using RAMSES, a finite-volume Godunov code based on the MUSCL-Hancock algorithm (Teyssier 2002; Fromang et al. 2006). The HLLD Riemann solver (Miyoshi & Kusano 2005), and the multidimensional slope limiter described in Suresh (2000) are used in all the simulations presented in this paper.

Rather than solving for the total -momentum, we evolve the equivalent conservation law for the angular momentum fluctuation , with the Keplerian velocity. An upwind solver is used for solving the azimuthal advection arising from . The tidal and Coriolis forces are treated as source terms and implemented following the Crank-Nicholson algorithm described in Stone & Gardiner (2010).

The algorithm solves for the total fluctuation energy and its conservation law is written as

(7) |

where is the total pressure tensor

(8) |

The left hand side of Equation (3.2) comprises the usual energy conservation law, which we solve using the MUSCL-Hancock algorithm. The treatment of the two terms on the right hand have been modified: the azimuthal advection of energy is solved with an upwind solver, and the second term involving the Maxwell and Reynolds stresses is added as a source term.

For the set of simulations shown in this paper, we used a box size of with a resolution of . is a reference scale height which is close to but not exactly the same as the scale height at the start of a simulation. It will be defined in detail later. The grid scale is defined to be . We set and in code units, where is the initial sound speed. The resolution is low in comparison to other MHD shearing box simulations, however it is sufficient to capture the basic properties of the MRI, in particular the turbulent fluctuations. Importantly, it is computationally inexpensive allowing for the simulations to run for orbits.

### 3.3 Thermal equilibrium

In our shearing box model, the energy is injected into the computational domain by the second term on the right side of Equation (3.2) which represents the liberation of shear energy by the total stress:

(9) |

Simulations have found that the box-averaged stress is roughly proportional to the box-averaged gas pressure to a given power (Sano et al. 2004, RLG16), though be aware of the caveats given in Section 2. Therefore, in numerical simulations, the averaged heating rate may be approximated by

(10) |

where and may be calculated from the simulation. Note that the former is not the same as the parameter. The exponent depends on the field geometry as well as the numerical parameters (RLG16). For our set-up, , as long as the pressure is sufficiently small (see later).

Combining this approximation with our cooling prescription, the evolution of the volume averaged pressure, , is determined by

(11) | ||||

(12) |

where

(13) |

In the formulation of Equation (11) we have made the approximation that , which is reasonable as long as the variance of within the box at fixed time is not too large. This system has two fixed points

(14) | ||||

(15) |

In shearing box simulations the thermodynamics have an additional complication: stress is independent of pressure once the scale height is larger than the radial and vertical box sizes (RLG16). The exponent may then be viewed as a function of the disk temperature, equal to for cool disks, while asymptoting to 0 as the gas heats and the scale height equals the box size . This behaviour introduces the possibility of an additional stable fixed point if . In Fig. 1 we plot (a) the heating rate as a function of pressure from a simulation (appearing in RLG16), (b) a smooth fit to this curve, and (c) overlay a power law cooling. Where the latter two curves intersect give three thermal equilibria. These are the and fixed points, described above, in addition to a third equilibrium, , which arises from the finite size of the box. Note that a similar fixed point will also appear in vertically stratified simulations.

The above volume averaged analysis is deterministic, but, turbulent flows are not. The stress fluctuates around its mean value as a result of the formation and break up of coherent structures within the flow. This means that can no longer be expressed in as simple a form as Equation (10). Instead, the stress is determined by the sum of a deterministic term and a fluctuating term. By using Equation (11) we can obtain an estimate of the equilibrium pressure that the system feels at any given time

(16) |

### 3.4 Thermal instability

According to a linear analysis of Equation (4), the thermal instability criterion is

(17) |

which, given our power law expressions for and , leads to the simple condition

for , and the opposite stability for and .

We calculate the growth rate as

(18) |

which can be used for comparison with the MHD simulations. The associated instability timescale we define to be .

If the system was truly laminar the instability criterion would be fully determined by and and the instability timescale from (18). Fluctuating systems, however, exhibit non-exponential and indeed non-monotonic behaviour which is poorly approximated by the laminar model. In particular, may be an unsuitable measure for the instability timescale. In its place we introduce the ‘escape time’ of a simulation, which is defined to be the last instance that the system lies within a pressure interval containing the fixed point. Mathematically it may be defined via

(19) |

where is one of the three equilibria introduced earlier, and is the interval size, outside of which we consider the system to have unequivocally departed from the fixed point. We are free to specify the size of , and it might reflect the particular problem of interest. For disk transitions between states differing by many orders of magnitude in temperature we might be generous with , permitting it to be up to 10 times the fixed point pressure. For smaller transitions, then must be smaller. Note that the system may stochastically dip in and out of this interval, but will capture the time when the system finally leaves it forever. If the disk is laminar, then

(20) |

where is the initial pressure of the system.

### 3.5 Initial conditions

If we want a thermal equilibrium to be achievable within our box and for the stress to depend appreciably on pressure, the initial conditions must be chosen carefully. An initial state that is too hot means the box size will unduly influence the stress and weaken ’s dependence on . As a consequence, thermal instability will fail to occur.

Stress is only observed to be a strong function of pressure () when (Sano et al. 2004, RLG16). Therefore, must be sufficiently low so that the stress is increasing with pressure, but sufficiently high so that the characteristic lengthscale of the turbulence is not on (or below) the grid. That way we can be assured that the thermal instability is able to appear in our numerical set-up. We choose

(21) |

where is the pressure at which the scale height equals the box size and is the scale height when . From this we can obtain the required value of for a given choice of :

(22) |

To achieve a turbulent state we first initialize a zero net-flux simulation with an initial field of and no cooling. In code units which we set with . Once a turbulent state is reached we switch on cooling with and . This choice of is to obtain a which we expect to be stable, as discussed in Section 3.4. These parameters lead to where is the initial pressure, having used that . During this steady state we calculate . This fully turbulent state in thermal equilibrium is used as our initial condition. The parameters and are then changed as appropriate.

## 4 Results

### 4.1 Stable system

To fix ideas we first consider the stable thermal equilibrium associated with the initial condition, , (simulation R1). We run the simulation for over orbits and the resulting stresses and pressures are shown in Fig. 2. As predicted from the laminar linear analysis, the equilibrium is ‘stable’ in the sense that both the stress and pressure fluctuate within some interval enclosing . Though the stress shows substantial variation during the simulation (with a maximum value of times its mean), there is no runaway or mean drift during the orbits. The variability in the pressure is less extreme, with a maximum value of times its mean value. The turbulence administers random ‘kicks’ to the system, but the deterministic physics always draws it back to the vicinity of . The stochastic fluctuations in stress are a result of dynamo cycles and the formation and breakup of coherent structures. The energy in the flow then cascades down to the dissipation scale where the magnetic and kinetic energies are converted to thermal energy leading to changes in pressure.

A short time delay of a few orbits () between the stress and pressure is clearly visible in Fig. 3, a result of the finite time taken for kinetic and magnetic energies to reach the dissipation scale from the injection scales (Hirose et al. 2009). The variations in pressure are smoother and longer than those in the stress. However, this short time dependency of pressure on stress is distinct to the longer time dependence of stress on pressure that drives instability/stability (Latter & Papaloizou 2012).

Though the pressure seems to be drawn back to the mean of , in actual fact the equilibrium balance that the system feels at any given instance is changing with time. A crude approximation to this is the instantaneous fixed point , calculated using Equation (16). But because the thermal timescale is longer than the turbulent timescale, a better approximation is the average of over the last thermal time. Using Equation (18) we estimate this time to be orbits. The resulting smoothed and time-dependent equlibrium is plotted in Fig. 4 in blue, alongside the actual pressure of the system in red. It is clear that follows the equilibrium with a time lag of orbits. Despite the fluctuations the system senses the stable fixed point and is attracted towards it, though because of those same fluctuations it can never come to rest upon it.

### 4.2 Unstable systems

Having explored the stable case we now proceed to an unstable thermal equilibrium. For this set of simulations we choose which we expect to be unstable based on the argument in Section 3.4.

We present 8 simulations (R2a-R2h) that have the same fully turbulent initial state. In order to vary the initial condition between runs, these simulations each have slightly different and so have slightly different . Each of these is within a few percent of the initial pressure. In practice this means that though each simulation starts from the same initial condition, each corresponds to a slightly different perturbation from equilibrium.

In Fig. 5 we show the evolution of the box-averaged pressures alongside the trajectory derived from the laminar linear theory. Unsurprisingly the turbulent fluctuations lead to a diversity of outcomes but do not indefinitely prevent thermal runaway. For example, a large kick can cause the system to escape from the fixed point on a shorter timescale than the laminar timescale, or the system may remain close to the fixed point for extended periods of time, longer than . Our most ‘stable’ simulation remains close to equilibrium for . None of the simulations can be well modelled by the laminar theory. For , the behaviour is closer to algebraic growth than exponential runaway. In fact, the behaviour of the system is strongly influenced by the fluctuations that occur on timescales comparable to the instability timescale, for these parameters. A more apt description of the system could be a biased random walk, when strong fluctuations in stress repeatedly perturb the system, while in between kicks the system drifts according to the deterministic physics.

If we consider some characteristic interval around the fixed point then we can find the very last time, , the system was within this band, Equation (19). If we define the interval to be rather narrow, , then this ‘escape time’ can be compared to . We find a wide range of escape times in our MHD simulations, from to . Ideally, we would calculate a probability distribution function for the escape time but this would require substantially more simulations, which, at this time is impractical. The choice of band is arbitrary, but it should be chosen to suit the system.

Once the simulations become cooler than the laminar model prediction. This can be attributed to the box size beginning to influence the evolution. At this point, the dependence of stress on pressure decreases (see Fig. 1) and hence the laminar model is an over estimate. This numerical effect introduces a third equilibrium point, , as discussed in Section 3.3. We do not observe a plateau in pressure associated with the system being attracted to , but, we expect that if run for a sufficient duration then a plateau would appear.

That a runaway heating is limited by the box size is an important problem, both for these simulations and stratified radiation MHD simulations. To explore this further, we rescale the simulations by choosing a larger value of , and initialize a simulation with the same initial turbulent state as in R2a-R2h. However, now the initial condition is much further away from the putative equilibrium, by a factor of some 3, and yet box effects remain negligible, because Equation (21) is still satisfied. In these runs the thermal runaway was somewhat faster than witnessed in Fig 5 for the hotter systems. These few simulations illustrate the numerical limitations inherent in any simulation of catastrophic heating undertaken in a finite domain.

### 4.3 Thermal fragmentation

In a turbulent system, heat is not deposited uniformly throughout the box, but rather is localised in coherent structures such as current and vorticity sheets. If the instability timescale, , is large compared to the mixing timescale, , then the heating inhomogeneity will have little effect as temperature fluctuations will be smoothed out. Conversely, if then regions of fluid can evolve independently of each resulting in localised runaway. In Fig. 6 we plot the pressure probability distribution function, , for two R2 simulations at . Inefficient turbulent mixing would result in the spreading of during runaway. The figure exhibits minimal evidence of spreading in the heating run, and the pressure ensemble evolves with a well defined and relatively narrow shape. The pressure in the cooling simulation is attracted to and in fact we see a further narrowing of the distribution as the equilibrium is approached.

To show localised runaway and fragmentation, we consider a case which we expect to be very unstable, choosing and (simulation R3). In Fig. 7 we plot at instances of time. The width of the distribution quickly increases, indicating localised thermal runaway: initially Var(), but after and orbits this grows to and respectively. During this time itself varies little, which emphasises that the box-averaged properties no longer give a satisfactory description of the state of the system. Soon after the final snapshot very small pressures occur resulting in the termination of the simulation, preventing further exploration. In Fig. 8 we show - slices in pressure before and during thermal fragmentation. Prior to thermal fragmentation there are strong acoustic waves propagating in the radial direction. These appear to break up into patches that undergo rapid thermal runaway independent of each other. Because this behaviour appears to be significantly nonlinear and disordered, we do not attribute it to the action of a linear instability mode with non-zero and .

Though it is straightforward enough to achieve fragmentation in an unstratified local box, how likely is this in the inner regions of an X-ray binary? In this context the two major contributors to mixing are radiative diffusion and turbulent advection. We first consider radiative diffusion. In the hot dense gas, the opacity is largely dominated by Thomson scattering with opacity . The radiative diffusion time across a length of is then where is the speed of light and is the mean free path. The instability timescale we estimate to be of order, but bounded below, by the thermal timescale, so that . If , we have the following condition on

(23) | ||||

(24) |

where is surface density. Regions separated by more than will be unable to mix sufficiently well on the instability timescale. Typical values for an X-ray binary indicate that regions apart my in fact thermally fragment.

What about turbulent mixing? We assume that the turbulent transport of heat by the MRI is similar in efficiency to its transport of angular momentum (though this is a point that has not been studied in detail). If we are permitted this assumption then the turbulent diffusion time scale is . If we next assume that the relevant eddies are of size (plausible if some form of MHD convection is operative) then on these outer scales. These very rough scalings indicate that turbulent heat transport is somewhat more efficient than radiative transport, and moreover that it may be sufficient to preclude fragmentation — though a simple order of magnitude treatment is unable to determine precisely when this might occur. Only realistic simulations themselves can decide on this issue, and in fact Jiang et al. (2013) do not find fragmentation. Our simulations are marginally susceptible, in agreement with the above argument, but they omit important physical effects such as buoyancy, which may be crucial here, and enhanced compressibility effects in radiation-dominated flow.

## 5 Reduced stochastic models

Because it is impossible to run a sufficient number of simulations to build reliable statistics, especially regarding the distribution of , we turn to simpler approximate models that illustrate more fully the effects of the fluctuations on stability and which also permit analytical results.

We work primarily with the averaged energy equation (12) but model the fluctuating turbulent stress via a random function . Our model is related to the logistic equation, and may be written as

(25) |

This can be derived from (12) by an appropriate rescaling, with and representing pressure/temperature and time, respectively. Constant parameters are , , and . To simplify the analysis while not losing much generality, we set in much of what follows.

Equation (25) admits the trivial steady state and the more interesting equilibrium . In the ‘laminar’ case of this equilibrium is unstable when with modes possessing the growth rate .

For a random but continuous function, (25) is the Bernoulli equation with analytic solution {dmath} xxeq=ϕ(τ)[(x0xeq)^1-m-(1-m)∫^τ_0ϕ(s)^m-1ds]^1/(1-m) where

(26) |

and is the initial value of . As it stands, the analytic solution is too unwieldy to be useful, even for basic prescriptions for , but it does illustrate clearly the competition between instability and stochasticity. These manifest as the two terms in the exponent of . The first describes the deterministic exponential runaway, while the second stochastic term potentially impedes this tendency. In fact, if behaves like a random walk, then its standard deviation will be proportional to , and so on short to intermediate times the second stochastic term can outcompete the first instability term. On longer times, however, will always defeat and the system will approach the unstable laminar solution, given by

(27) |

### 5.1 Geometric Brownian motion (GBM)

Before presenting an analysis of the full equation (25), it is worthwhile examining the simpler case of and . The resulting system isolates cleanly all the main characteristics of more realistic systems — an unstable fixed point ( now) and stochastic noise — while being analytically tractable. Equation (25) becomes

(28) |

with initial condition . For smooth random , the solution is

(29) |

When studying stochastic dynamical systems, white noise is a convenient choice for modelling the variability. For white noise to be a good approximation, the fluctuation timescales of the system should be much less than the characteristic time of interest. For MRI-driven turbulence, this choice is not ideal given that the spectrum of has preferential frequencies. However, we use this as our starting point as it makes a number of results especially clear. With this choice of , Equation (28) must be written in differential form

(30) |

where is the volatility coefficient (or noise amplitude) and is white noise. Here, for simplicity, we have interpreted the calculus in the Ito sense. Equation (30) actually describes geometric Brownian motion and is frequently used in financial modelling. Its solution is

(31) |

Given that , the solution remains strictly positive as a result of the multiplicative form of the noise. We plot sample trajectories in Fig. 9 along with the th and th-percentiles. When calculating these trajectories we use the Euler-Maruyama method (Kloeden & Planten 1992). A feature of this collection of sample paths is the wide variation between them, an attribute that is shared with the simulations shown in Fig. 5. They are also non-monotonic; fluctuations ‘kick’ the system towards or away from the fixed point. The light blue curve is particularly striking, exhibiting a trajectory that remain close to equilibrium, , up to time . Note that at a purely deterministic model would have predicted to be , an order of magnitude greater. The stochastic term in Equation (31) has ‘balanced out’ the deterministic drift, at least on these shorter times.

The reader may note that when , the stability of switches. It becomes an attractor, and the system is stabilised. We stress, however, this effect is an artefact of multiplicative white noise in combination with the Ito calculus, and is not to be expected in real turbulent systems. For instance, the stabilisation vanishes in the Stratonovich calculus and/or with noise models with memory and which are not multiplicative. We certainly do not expect MHD turbulence to exhibit the combination of special features that leads to this stabilisation. Indeed, it makes little physical sense that ‘shaking’ an unstable system more vigorously ultimately leads to zero fluctuations. Moreover, the required amplitude of the fluctuations must be extremely large, in our case this would require negative which is impossible. It is worth pointing out that the auto-regressive stochastic model employed by Janiuk & Misra (2012) shares the same stabilising property as white noise in the Ito calculus; consequently, we view their stabilisation of thermal instability as an artefact of their model, and not representative of a real fluctuating disk system.

The probability distribution of the solution trajectories may be obtained by solving the associated Fokker-Planck equation. If is the probability of finding a path lying between and at time , then

(32) |

The expectation, , and variance, Var, can also be calculated

(33) | ||||

(34) |

The mean is independent of the volatility and hence agrees with the laminar model. However, the variance contains a factor which grows exponentially! This means that as time progresses the expectation value grows less and less meaningful because the distribution becomes increasingly wide and flat. Overall, trajectories move away from the unstable equilibrium point, but individual trajectories can deviate from the laminar model dramatically.

Next we turn to the statistics of the escape time (or ‘last hitting time’), which can be defined as follows:

(35) |

where . This gives us the last time that a sample path is within the interval . Thus provides a measure for the effective instability timescale, more accurate than the inverse of the laminar growth rate. If the system was laminar (), however, the escape time would be simply .

For geometric Brownian motion, it is possible to derive the probability distribution of analytically (Kennedy 2010, Profeta 2010). Because this is not a standard calculation we go through its details in the Appendix. Denoted by where , it is a modified form of the Rayleigh distribution:

(36) |

where the parameters , , and have conveniently combined into . In Fig. 10 we plot for a range of noise amplitudes.

The mean and variance of can be calculated analytically

(37) | ||||

(38) |

It is clear that the mean of the distribution does not depart greatly from the laminar escape time. Though the variance does increase with the volatality of the noise, the spread is not especially dramatic. Perhaps a more illuminating quantity is the kurtosis:

(39) |

which varies from 3, when , to 9, when . A Gaussian possesses a kurtosis of 3, so our escape time distribution can in fact be exceptionally ‘fat-tailed’, or ‘leptokurtic’, meaning it generates a significant number of outliers. Fig 10 illustrates this point, with larger giving rise to very skewed and broad distributions. Inserted in the upper right of the figure we also plot the cumulative distribution function for various values of . For , in particular, we see that only 30% of systems possess an escape time equal to or less than the laminar escape time (), while over 10% of systems possess an escape time of four times or more the laminar escape time. In fact, the probability that a trajectory possesses a greater than some value may be approximated by the expression

for large . When and , the probability is only about 1%. However, this rises to 10% when , indicating the strongly nonlinear dependence of the statistics on noise amplitude. In summary, through the statistics of the escape time, one can observe stochasticity impeding instability over some initial period of the evolutio.

### 5.2 Random logistic equation

We return to the more general problem by reintroducing the cooling term with

(40) |

Again assuming white noise for and setting , we rewrite the ODE in differential form

(41) |

The equation was solved numerically, and sample trajectories are shown in Fig. 11. With a fixed initial condition, the system can either undergo runaway heating or cooling. Note that the initial perturbation does not give a good indication of which direction the sample paths eventually go. The trajectories in Fig. 11 qualitatively resemble those in Fig. 5 from the MHD simulations. Both show trajectories that remain close to their equilibrium values for multiple thermal timescales along with trajectories that escape from the equilibrium faster than the laminar model. One substantial difference exists however: Fig. 11 shows no indication that the trajectories will be in general slower than the laminar model as is the case in Fig. 5.

Unlike geometric Brownian motion, it is difficult to analytically calculate the probability density function, and instead this is accomplished by solving the Fokker-Planck equation numerically

(42) |

In order to do this, an approximation for the initial distribution must be made. Rather than using a Dirac- function as the initial condition it is necessary to use a somewhat smoother function, a Gaussian distribution with variance . An example solution for the probability distribution is plotted in Fig. 12. As earlier, the distribution becomes increasingly wide and flat as time progresses, indicating the broad range in possible evolutionary paths.

With the introduction of cooling comes the chance of a sample path going to zero, hence, our escape time definition must be modified to include a lower boundary

(43) |

where . For a given sample path, the escape time is the last instance when the solution is within the interval , containing . In Fig. 13 we show estimated probability density functions of for a range of . These are calculated from sample paths for each . As expected, for small we approach the laminar escape time, while as increases, undergoes asymmetric broadening, shifting its maximum to lower and developing a tail at long . Qualitatively this behaviour mirrors that shown for geometric Brownian motion, Fig 10. And so our conclusions for GBM carry over to the more realistic logistic case.

#### 5.2.1 Variable and

In Equation (25), and are free parameters which can be chosen to fit the system of interest. This flexibility allows it to approximately describe a range of different scenarios. For instance, the classical Shakura & Sunyaev disk (Shakura & Sunyaev 1973) can be represented with and . Alternatively, and could be chosen so as to model the vertically stratified, radiation-pressure dominated simulations of Jiang et al. (2013, 2016). In this case a large surface density disk yields and , while a less dense disk gives and . Finally, when the temperature of the gas ensures the opacity is influenced by the ‘iron bump’, the scaling of cooling with central pressure is found to greatly exceed that of the classical Shukura & Sunyaev model and (Jiang et al. 2016). In Fig. 14 we plot the probability distribution function of the escape time for various choices of and . As approaches from below, the escape time distribution becomes increasingly elongated to large . The deterministic component weakens and the stochastic drift becomes more important, dominating the results.

To better quantify this effect we compute the probability that a given trajectory possesses a greater than various multiples of and plot these probabilities as a function of for fixed . These results appear in Fig. 15. The curves show that as approaches the escape time can be significantly enhanced. For there is approximately chance that it is five times the laminar prediction, and chance that is ten times greater. It should be noted that as approaches the laminar instability timescale itself can be significantly longer than the thermal time , which further separates the expected turbulent from the thermal time. The stability uncovered in Jiang et al. (2016), on the timescales of their simulations, can be easily explained by the enhanced delay witnessed in such marginally unstable systems.

### 5.3 MHD power spectral density

The next step is to replace the white noise with the power spectral density (PSD) of calculated from the MHD simulations of Section 4. In Fig. 16 we show the PSD of from simulation R1. We then represent the fluctuating term in Equation (25) by

(44) | ||||

(45) |

where are the frequencies, logarithmically spaced between , and are phase shifts, chosen randomly from a uniform distribution on the interval . Finally, the constant amplitudes are calculated from a two slope power law that fits our MHD simulations (constant for and otherwise), which is then multiplied by a random number generated from a uniform distribution on the interval .

To improve the comparison with the simulations R2a-R2h, we set and and thus our model equation is

(46) |

where is defined in Equation . This equation is evolved forward in time times for each choice of in order to derive adequate statistics, especially for the escape time.

We show the resulting escape time distribution in Fig. 17. We find similar behaviour to both the GBM and the random logistic equation. As increases, the distribution broadens and its tail at large increases. We estimate that gives approximately the correct fluctuation amplitude when compared to the stress in Fig. 2. For this choice, Fig. 17 shows a clear tail reaching .

We conclude that fat-tailed distributions are a generic feature of the escape time in turbulent but thermally unstable systems. Such systems produce a broad range of outcomes, and instability can be delayed for several instability timescales. Being fat-tailed they also exhibit significant outliers — systems that ‘hang around’ for surprisingly long times before wandering away.

## 6 Discussion and Conclusions

We have performed a set of idealised shearing box simulations of the MRI in order to explore the effects of turbulent variability on thermal instability. Our main aim was to check if turbulence interferes with the thermal runaway predicted by laminar theory. In our simulations heating comes from numerical dissipation while a power law cooling imitates optically thin radiative cooling. Relatively large computational domains are used in order to ensure a strong dependence of stress on pressure.

Simulations with an expected stable thermal equilibrium are found to fluctuate around their fixed points. When the cooling power law exponent is decreased, and the laminar analysis predicts instability, our simulations indeed show thermal runaways. However, the system trajectories deviate from the laminar theory. The turbulent fluctuations can stall the onset of instability in a large number of cases, for multiple thermal times. To better account for the instability timescale, we introduce the concept of the escape time which we define to be the last time the system leaves an interval encapsulating the equilibrium in phase space. Our simulations show a large range of escape times ranging from to times the laminar thermal timescale, for relatively narrow intervals.

Further reducing the cooling power law exponent results in disk fragmentation. This is a due to localised imbalances between heating and cooling; the instability timescale is shorter than the mixing timescale and hence distant pockets of fluid evolve independently. Very rough estimates indicate that thermal fragmentation on length scales larger than are at best marginally possible in the inner regions of X-ray binaries.

To better understand our results we construct a probabilistic theory centred on simple stochastic equations that approximate the box-averaged thermal equation. These present us with a much larger sample of possible system outcomes to analyse than the MHD simulations can afford. First we analyse geometric Brownian motion (GBM), which possesses the main ingredients of interest (an unstable fixed point and stochastic variations), while remaining analytically tractable. The distribution function of exhibits a variance proportional to the square of the noise amplitude and considerable kurtosis. In general the distribution is ‘fat-tailed’, permitting many instances of delayed thermal instability, and outliers for whom the escape time can be thermal times. Our second model introduces a logistic term to incorporate a power law cooling, and a third model replaces the white noise with the power spectral density of the stress from a thermally stable shearing box simulation. In both cases we obtain qualitatively similar behaviour to before, which instils confidence that its behaviour is generic to systems sustaining noise and instability. GBM may be thought of as the model equation for such systems.

In our GBM model, very large amplitude noise can stabilize the thermal instability, but we believe this is physically implausible. In fact, this behaviour arises from the special combination of Ito calculus and multiplicative white noise. It is not generic. While GBM is a convenient model for ‘unresolved’ turbulence, this dynamical peculiarity must be discarded when applying results to real systems. The stochastic process employed by Janiuk & Misra (2012) shares the same properties and thus suffers the same dynamical artefact. The stabilisation witnessed in their simulations we hence view as unphysical.

The bulk of X-ray binary observations show no indication of limit cycles that could correspond to radiation-pressure induced thermal instability (Gierliński & Done 2004). Notable exception are GRS 1915+105 (Done et al. 2004) and HLX-1 (Sun et al. 2016). We show that it is unlikely that turbulent fluctuations alone stabilise these disks. However, turbulence can delay and weaken instability. This weakening might be significant in combination with additional stabilising mechanisms, for example: energy lost in the disk corona or by outflows, magnetic buoyancy effects, alteration of the pressure-stress relationship by strong magnetic fields, or opacity shifts near the iron bump (Svensson & Zdziarski 1994, Sadowski 2016, Jiang and Stone 2016). Finally, separate and important physical effects may arise from the global nature of the flow, especially from accretion. The viscous timescale is probably longer than a typical , as measured in this paper. But if further delayed by additional physics, could in some circumstances approach the accretion time. If this occurs then the classic limit cycle behaviour expected from thermal instability could well be impeded, and/or pushed to smaller radii.

## Acknowledgements

The authors thank the anonymous reviewer for a set of helpful comments. They also thank Yanfei Jiang and Jim Stone for discussions that helped sharpen the paper. This work was partially funded by STFC grants ST/L000636/1 and ST/K501906/1. Some of the simulations were run on the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E- Infrastructure capital grant ST/K000373/1 and STFC DiRAC Op- erations grant ST/K0003259/1. DiRAC is part of the UK National E-Infrastructure.

## References

- [1] Balbus S. A. and Hawley J. F., 1991, ApJ, 376, 214-233.
- [2] Belloni T., Méndez M., King A. R., van der Klis M., van Paradijs J., 1997, ApJL, 488, L109-L112.
- [3] Ciesielski A., Wielgus M., Kluźniak W., Sa̧dowski A., Abramowicz M., Lasota J.-P. and Rebusco P., 2012, A& A, 538, A148.
- [4] De Swart H. E., Grasman J., 1987, Tellus, 39A, 10-24.
- [5] Done C., Wardziński G., Gierliński M. 2004, MNRAS, 349,393-403.
- [6] Fromang S., Hennebelle P., Teyssier R., A&A, 2006, 457, 371-384.
- [7] Fromang S., Papaloizou J. C. B., A&A, 2007, 476, 1113-1122.
- [8] Gierliński M., Done C., MNRAS, 2004, 347, 885-894.
- [9] Goldreich P., Lynden-Bell D, 1965, MNRAS, 130, 125.
- [10] Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440.
- [11] Hirose S., Krolik J. H., Blaes O, 2009, ApJ, 691, 16-31.
- [12] Janiuk A., Misra R., 2012, A&A, 540, A114
- [13] Jiang Y. F., Stone S. M., Davis. S. W., 2013, ApJ, 778, 65.
- [14] Jiang Y. F., Stone S. M., Davis. S. W., 2016, ApJ, 778, 10.
- [15] Kennedy D., 2010, Stochastic Financial Models, Chapman & Hall/CRC.
- [16] King A. R., Ritter H. R. 1998, MNRAS, 293, L42-L48.
- [17] Kloeden P. E., Platen E., 1992, Numerical Solution of Stochastic Differential Equations, Applications of Mathematics, Vol 23, Springer Berlin Heidelberg.
- [18] Latter H. N., Papaloizou J. C. B., 2012, MNRAS, 426, 1107-1120
- [19] Lesur G., Ogilvie G. I., 2009, A&A, 488, 451-461.
- [20] Levins R., 1969, PNAS, 62, 1061-1065.
- [21] Lightman A. P, Eardley D. M., 1974, ApJ, 187, 13-66.
- [22] Lin D.-B., Gu W.-M., Lu J.-F., 2011 ,MNRAS ,415, 2319-2322.
- [23] Majda A. J., Kramer P. R., 1999, Phys. Reports, 314, 237-574.
- [24] Majda A. J., Timofeyev I., Vanden-Eijinden E., 2003, J. Atmos. Sci., 60, 1705-1722.
- [25] May R. M., 1973, Stability and complexity in model ecosystems, Princeton Univ. Press, Princeton, NJ.
- [26] Mishra B., Fragile P. C., Johnson L. C. and Kluźniak W., 2016, MNRAS, 463, 3437-3448.
- [27] Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315.
- [28] Nowak M. A. et al., 1999, ApJ, 510, 874.
- [29] Piran T., 1978, ApJ, 221, 652-660.
- [30] Profeta C., Roynette B., and Yor M., Option Prices as Probabilities: A New Look at Generalized Black-Scholes Formulae, Springer, 2010.
- [31] Ross J., Latter H. N., Guilet J., 2016, MNRAS, 455, 1, 526-539.
- [32] Sadowski A., 2016, MNRAS, 459, 4397-4407.
- [33] Sano T., Inutsuka S., Turner N. J., Stone J. M., 2004. ApJ, 605, 321.
- [34] Shakura N. I., Sunyaev R. A., 1973, A&A., 24:337.
- [35] Shakura N. I., Sunyaev R. A., 1976, MNRAS, 175, 613-632.
- [36] , Smak J., 1984, PASP, 96, 5-18.
- [37] Stone J., Gardiner T., ApJSS, 2010, 189, 142-155.
- [38] Sun M., Gu W.-M., Yan Z., Wu Q., Liu T., 2016, MNRAS, 463, L99-L102.
- [39] Suresh A., 2000, SIAM J. Sci. Comp., 22, 1184.
- [40] Svensson R., Zdziarski A. A., 1994, ApJ, 436,599-606.
- [41] Teyssier R., 2002, A&A, 364, 337-364.
- [42] Warner B., 1995, Catacylsmic variables, Cambridge University Press, Cambridge.
- [43] Wu, Q., Czerny, B., Grzedzielski, M., Janiuk, A., Gu, W.-M., Dong, A.-J., Cao, X.-F., You, B., Yan, Z. and Sun, M.-Y., 2016, ApJ, 833, 79.

## Appendix A Escape time distribution for GBM

Following Kennedy (2010), we derive the probability density function (PDF) for the escape time, or last hitting time, for geometric Brownian motion. To achieve this, we calculate the PDF of the last hitting time of Brownian motion with drift and then perform a change of variables to obtain the last hitting time for GBM.

In the appendix will refer to a standard Brownian motion, the amplitude of the fluctuations can be represented by a pre-factor , which we refer to as the ‘volatility’. To begin, consider the first hitting time of a random function , denoted by and defined by

(47) |

(Here the infimum can be thought of as the minimum.) Physically this represents the first time the function reaches the value . In the special case of corresponding to standard Brownian motion with drift, , the probability density function, , of is

(48) |

See Kennedy (2010) for a derivation. With this result, we can then determine the last hitting time of Geometric Brownian motion with drift.

Let denote the last time that the standard Brownian motion hits the line , for some constant . Then the probability density function of is given by

(49) |

To see why this is true consider the probability that , for some :

In the second equality we have used the fact that is also a Brownian motion. The probability distribution function of is given by the derivative of the above expression. By considering the probability distribution function of where the drift is and the threshold , we conclude that

(50) |

from which the result follows. Notice that is identical to the probability density function of the last time Brownian motion with drift hits the line .

Let us now finally define the last hitting time, or escape time, for a random function :

(51) |

for some threshold . If corresponds to Brownian motion with drift , for constant , then the probability distribution of is

(52) |

If however, is a GBM

(53) |

and a constant, then the distribution of , for some constant threshold , is easy to obtain from (52). Consider the probability for GBM that , for some . This corresponds to

which is precisely the probability for Brownian motion with drift, but with and . This supplies us with an expression for the distribution of for GBM. A rescaling of and a renormalisation obtains Equation (36).