The stress–pressure lag in MRI turbulence and its implications for thermal instability in accretion discs

ABSTRACT

1 INTRODUCTION

The alpha-disc model has underpinned the study of accretion discs for several decades (Shakura & Sunyaev 1973). It has permitted researchers to develop semi-analytic theories with which to interpret observations and also direct numerical simulations. Despite its successes, this essentially mean field model makes strong assumptions that problematise its application to phenomena on lengthscales of order the disc scale height or shorter, and on time-scales shorter than the viscous time (Balbus, Gammie & Hawley 1994; Ogilvie 2003). Instabilities, in particular, are potentially misrepresented, with the most famous example being thermal instability in radiation-pressure dominated accretion flows: while alpha-disc models predict instability and subsequent non-linear oscillations (Lightman & Eardley 1974; Shakura & Sunyaev 1976; Honma, Matsumoto & Kato 1991; Szuszkiewicz & Miller 2001), X-ray observations generally fail to find variability on the time-scales expected (e.g. Gierliński & Done 2004).

Several solutions have been proffered to explain the thermal instability conundrum (e.g. see discussion in Ross, Latter & Tehranchi 2017), but in this paper, we focus on the time lag between variations in pressure and variations in the turbulent stress, an effect that can weaken the instability (Lin, Gu & Lu 2011; Ciesielski et al. 2012). The alpha disc prescription assumes that the stress responds instantly to the pressure, but in reality it will take time for variations in one field to be communicated to the other through the complicated tangle of turbulent motions and magnetic fields. This lag is connected (and probably equal to) the stress’s relaxation time [e.g. Crow (1968)], and if this is similar to the time-scales of interest it may be necessary to solve a separate evolution equation for the stress itself, rather than adopt the oversimple alpha prescription (Kato & Yoshizawa 1993; Ogilvie 2003; Pessah, Chan & Psaltis 2006).

In this paper, we measure the relaxation time of the turbulent stress in local simulations of the zero-net-flux magnetorotational instability (MRI), with the finite volume Godunov code pluto (Mignone et al. 2007). We construct numerical experiments in which we artificially drive oscillatory variations in pressure on various time-scales, and then see how quickly or slowly the turbulent MRI stress responds. The choice to drive such cycles is not motivated by any particular astrophysical phenomena, but because such cycles make the calculation of the intrinsic lag easier. Our results show that driving on intermediate to long time-scales (of order or greater than the thermal time), the turbulent stress lags behind the pressure by only about 5 orbits. If these results carry over to more realistic disc environments (especially those incorporating radiation pressure) then the lag is too short to impact on thermal instability unless the turbulent alpha is ≳ 0.1 (i.e. the dynamical and thermal times are similar), which can be the case in certain substates of X-ray binaries in outburst (e.g. Done, Gierliński & Kubota 2007). The lag then weakens thermal instability, and in combination with other processes (magnetic fields, the stress’s inherent stochasticity, etc.) may potentially suppress it (Begelman & Pringle 2007; Oda et al. 2009; Ross et al. 2017). Finally, though only tangentially applicable to viscous overstability (Kato 1978; Blumenthal, Lin & Yang 1984), our simulations suggest that our measured lags may be more than sufficient to complicate its onset in MRI-turbulent gaseous discs (Ogilvie 2001).

The structure of the paper is as follows. In Section 2, we review the theoretical context and discuss thermal instability and viscous overstability. In Section 3, we discuss our set up, initial conditions, and diagnostics. Section 4 presents our results. First, we confirm the stress–pressure relationship reported in Ross, Latter & Guilet (2016) in simulations that employ either no explicit cooling (heating runs), or otherwise constant time cooling (cooling runs). Next, we show simulations with driven thermal oscillations in the mean pressure and study the resultant stress response, and any time lags that appear. Finally, we discuss the results and conclude in Section 5.

2 BACKGROUND

This section sketches out the context and motivation for this work, and explains some theoretical expectations and applications.

2.1 Pressure and stress variations

In local simulations of the MRI, the turbulent stress and the pressure can engage in a two-way interaction, but with each reacting upon the other on slightly separated time-scales. On shorter times, of order an orbit or so, the MRI’s turbulent stress undergoes stochastic fluctuations of non-negligible amplitude, and in particular in strong bursts. These bursts can be communicated to the pressure via heating, in which case one observes the stress acting on the pressure (Hirose, Krolik & Blaes 2009; Latter & Papaloizou 2012). But it is important to recognize that the stress fluctuations are essentially random, and cumulatively don’t impart any mean trend if the gas is in thermal equilibrium. A slight lag exists between bursts in stress and bursts in pressure, because of the finite time it takes for the excess energy to tumble down the turbulent cascade and to thermalize.

If, however, on longer times there is a mean change in the average pressure, caused by an external mechanism or by instability for example, then the MRI’s turbulent stress follows the pressure. This has been demonstrated in local simulations (of sufficiently large radial size) describing the pure heating or excessive cooling of the turbulent gas (Ross et al. 2016) and of thermal instability (Jiang, Stone & Davis 2013; Ross et al. 2017). The link between pressure and stress these simulations show is bound up in the MRI saturation mechanism, though its exact details are yet to be fully understood. At the very least, it appears that compressibility limits the size of the largest MRI turbulent eddies and hence the magnitude of the subsequent stresses, though there remains evidence that the smallest scales are also involved (Fromang & Papaloizou 2007; Fromang et al. 2007; Ross et al. 2016; Ryan et al. 2017; Ross & Latter 2018).

It has yet to be established numerically that there exists a time lag in the second slower interaction described, i.e. between the response of the stress to the pressure. One might speculate that such a time lag should always be of order or somewhat longer than the dynamical time-scale (Ogilvie 2003). As one expects the lag to be set by the MRI saturation mechanism, at the very least it should be bounded below by the longest eddy turnover time ∼Ω⁻¹ and bounded above by the thermal time ∼(αΩ)⁻¹, where Ω is the local orbital frequency and α is the alpha parameter.

2.2 Thermal instability

We now show how a time lag in the stress impacts on thermal instability in discs (Lightman & Eardley 1974; Shakura & Sunyaev 1976), using a simple but transparent 1D toy model, based on analyses by Lin et al. (2011) and Ciesielski et al. (2012). A more complete account might include two coupled (possibly stochastic) equations for the stress and the energy (cf. Hirose et al. 2009).

Consider the energy equation for a specified blob of gas in the disc: de/dt = Γ − Λ, where e is the blob’s thermal energy, Γ is its heating rate (proportional to the turbulent viscous stress), and Λ is the cooling rate. Suppose the viscous stress depends on pressure (and thus e) – but on its value a constant time τ in the past (as discussed in Section 2.1). It follows that the heating rate will also depend on the energy a time τ ago, and so we write Γ = Γ [e(t − τ)]. We also take Γ to be a monotonically increasing function of e. It is assumed, however, that the cooling depends on the instantaneous thermal energy, and thus Λ = Λ(e).

If there is no lag (τ = 0) in equation (1), then we have the usual result: s = Γ′ − Λ′, with instability occurring when Γ′ > Λ′. If satisfied, a small increase in energy will be reinforced because the change in heating outstrips the change in cooling; similarly, a small decrease in energy will be exacerbated. We hence denote s₀ = Γ′ − Λ′.

Figure 1.

Representative growth rate of thermal instability s as a function of the time lag τ. Here s₀ is the growth rate achieved with no lag, and Γ′ represents the inverse thermal time-scale. We have chosen Λ′/Γ′ = 1/2 in (2).

Open in new tabDownload slide

Given the above analysis, we argue that the lag time τ poses a problem for thermal instability when it is of order or greater the thermal time 1/Γ′, i.e. if Γ′τ ≳ 1. We next add to this the assumption τ ≳ Ω⁻¹ (cf. Section 2.1), i.e. the lag is naturally of order or greater than the dynamical time. Combining the two estimates, yields a relatively straightforward (and rather strict) condition for the importance of the lag: the thermal time must be sufficiently close to the dynamical time. In a modified alpha-disc prescription with a retarded pressure dependence, we have Γ′ = αΩ, and thus the condition for the lag’s importance is roughly α > 0.1.

Note that in previous local MRI simulations of thermal instability (e.g. Jiang et al. 2013; Ross et al. 2017) the measured alpha was ∼0.01, and thus too small for the lag effect to be noticeable.¹ On the other hand, in order to match observed light curves of X-ray binaries we must have α ≳ 0.1 in the outbursting state (Lasota 2001), which includes the radiation-dominated substates that should be thermally unstable (e.g. Done et al. 2007). It is in these contexts that we might expect the lag to reduce the instability growth rate and lengthen its characteristic time-scale, though certainly not by more than an order of magnitude.

2.3 Viscous overstability

Though this paper restricts itself to measuring the delay between pressure variations and stress variations, a lag in those quantities should be indicative of lags in strain and density as well, and more generally the stress’s finite relaxation time. Thus a second application of our investigations is to the viscous overstability (Kato 1978; Blumenthal et al. 1984), which is sensitive to a non-negligible stress relaxation. The viscous overstability is essentially a growing density wave (f-mode; Ogilvie 1998): the wave produces an oscillatory perturbation in the stress, which extracts energy from the background orbital shear and directs it back into the wave. The instability can lead to eccentricity growth in both narrow rings and extended gaseous discs (Borderies, Goldreich & Tremaine 1985; Papaloizou & Lin 1988; Lyubarskij, Postnov & Prokhorov 1994; Ogilvie 2001), global oscillations in X-ray binaries (Chan 2009; Miranda, Horák & Lai 2015), and fine-scale structure in Saturn’s rings (Schmit & Tscharnuter 1995; Latter & Ogilvie 2009).

In order for there to be an appreciable energy flow from the background shear to the wave, the turbulent stress oscillation and the dynamical oscillation need to be sufficiently in phase: a significant lag between them can quench instability (Ogilvie 2001; Latter & Ogilvie 2006). In dense planetary rings, the two quantities respond almost instantaneously, as shown by N-body simulations and kinetic theory (Salo 2001; Salo, Schmidt & Spahn 2001; Latter & Ogilvie 2008), and hence instability can proceed unproblematically. But, as discussed in Section 2.1, MRI turbulence will probably exhibit a relaxation/lag time of order or longer than the underlying inertial-acoustic oscillation, and this is likely to be sufficient to complicate the onset of overstability in turbulent gaseous discs (Ogilvie 2001).

3 METHODS

3.1 Governing equations

3.1.1 Cooling Prescription

Using this simple time-dependent cooling prescription, we can control the volume-averaged pressure on time-scales comparable to τ_H and τ_C, and in fact force it to oscillate between a minimum and maximum value of 〈P〉₊ and 〈P〉₋, respectively. Thus cooling times are free parameters which we vary in order to control the period of the pressure oscillations: these can be significantly longer than the orbital period or pushed to a time-scale near an orbit. In practice, different combinations of τ_C and τ_H yield similar thermal oscillations. The choices that appear later in the paper were arrived upon by trial and error.

3.2 Numerical set-up

3.2.1 Code

We use the conservative, finite-volume code pluto (Mignone et al. 2007) with the HLLD Riemann solver, 2nd-order-in-space linear interpolation, and the 2nd-order-in-time Runge–Kutta algorithm. In order to enforce the condition that |$\nabla \cdot \mathbf {B}=0$|⁠, we employ constrained transport (CT), and the UCT-contact algorithm to calculate the EMF at cell edges. To allow for longer time steps, we take advantage of the fargo scheme (Mignone et al. 2012). When explicit viscosity and resistivity are included, we further reduce the computational time via super-time-stepping (STS) (Alexiades, Amiez & Gremaud 1996). Ghost zones are used to implement the boundary conditions.

pluto solves the governing equations in conservative form, not the primitive form given by equations (4)–(7). Moreover, it evolves the total energy equation rather than the thermal energy equation. Note that due to the conservative implementation of the equations, kinetic and magnetic energy is not lost to the grid but converted directly into thermal energy. In pluto, the cooling term Λ is not implemented directly in the total energy equation. Instead, it is included on the right-hand side of the thermal energy equation, which is then integrated (in time) analytically.

3.2.2 Initial conditions

All our simulations are initialized from an equilibrium exhibiting spatially uniform density and pressure profiles, ρ₀ and P₀. The background velocity is |$\mathbf {u} = -(3/2) \Omega _0 x \, \mathbf {e}_y$|⁠. At initialization, we perturb all the velocity components with random noise exhibiting a flat power spectrum. The perturbations |$\delta \mathbf {u}$| have maximum amplitude of about |$0.1\, c_{s0}$|⁠, unless stated otherwise. Here |$c_{s0}=\sqrt{\gamma P_0/\rho _0 }$| is the sound speed at initialization.

All simulations are initialized with a zero-net-flux magnetic field configuration: |$\mathbf {B}_0 = B_0\sin {(k_x x)}\mathbf {\hat{e}}_z$|⁠. We take the radial wavenumber to be k_x = 4(2π/L_x), where L_x is the radial size of the numerical domain (see Section 3.2.4). The field strength B₀ is set by the ratio of gas pressure to magnetic pressure at the mid plane, |$\beta _0 \equiv P / (B_0^2 /2)$|⁠. In our simulations we set β₀ ≡ 1000.

3.2.3 Units

All quantities in the following sections are given in terms of dimensionless code units. The unit of time is selected so that Ω₀ = 1, where from now on, we drop the subscript on the angular frequency. The length unit is chosen so that the initial mid plane sound speed c_s0 = 1, which in turn defines a reference scale-height H₀ ≡ c_s0/Ω₀ = 1. Note, however, that the sound speed (and the scale height) is generally a function of both space and time. Finally the mass unit is set by the initial density so that ρ₀ = 1.

3.2.4 Diffusivities

We employ both an explicit magnetic diffusivity of η = 2 × 10⁻⁴ and an explicit viscosity of ν = 8 × 10⁻⁴. These values correspond to a Reynolds number of Re = c_s0H₀/ν = 1250, a magnetic Reynolds number of Rm = c_s0H₀/η = 5000, and a magnetic Prandtl number of Pm = ν/η = 4. The magnetic Prandtl number was chosen to be sufficiently high so that the MRI does not switch off (Fromang & Papaloizou 2007). These values are the same as those employed in the 256³ simulations of Ross et al. (2016) with explicit diffusion coefficients.

3.2.5 Box size, resolution, and boundary conditions

Our boxes take a size of L_x × L_y × L_z, with L_x = L_y = L_z = 4H₀, with a resolution of 256³ (i.e. 64 cells per scale height). We use shear-periodic boundary conditions in the x-direction (see Hawley et al. 1995) and periodic boundary conditions in the y- and z-directions. In principle, mass should be conserved in this set up, however we discovered that periodic boundary conditions in pluto paired with the fargo algorithm does lead to some very small mass loss (⁠|$\lt 0.3{{\ \rm per\ cent}}$| in 500 orbits).

3.3 Diagnostics

To calculate the cross correlation between the volume-averaged stress and pressure, we first smooth the stress (which is quite bursty on time-scales ∼1 orbit) using a lag-free rolling window average (see Pandas Development Team 2020). We then measure the cross correlation (Pearson r) between the smoothed stress and the pressure by shifting the pressure time series and repeatedly measuring the correlation between pressure and stress (see Cheung 2019).

4 RESULTS

We first reproduce previous results by (Ross et al. 2016), to make sure that our MRI turbulent stresses follow the pressure in pure heating and pure cooling runs. The rest of the section describes thermally driven runs, using the prescription in equation (9), in which the pressure (and hence the stress) exhibits long time oscillations. Analysing these time series, we can characterize the lag between stress and pressure.

4.1 Pure heating and cooling runs

We first carry out several runs with heating only (i.e. Λ = 0). In these simulations, the pressure increases monotonically (due to MRI turbulent dissipation). Our aim is to check that the stress also increases and to measure its power law dependence on pressure, i.e. Π_xy ∝ P^q.

Our fiducial heating simulation adopts the parameters presented in the last section. We find that its volume-averaged pressure increases by a factor of around 7 over 40 orbits, reaching 〈P〉 ∼ 4.8 just before the box-dominated regime begins (i.e. when the sound speed is of order or greater than L_xΩ). The stress exhibits more complicated behaviour, it fluctuates over short time-scales (≲ 1 orbit), nevertheless it too clearly increases after the MRI’s initial breakdown into turbulence. Once the simulation reaches the box-dominated regime, the stress plateaus and remains quasi-steady in time, though still undergoing fairly large fluctuations of around 0.01 in magnitude. In summary, there is a correlation between the stress and pressure, as in Ross et al. (2016).

To determine the stress-pressure scaling, we plot the stress and pressure against each other in Fig. 2 (the red curve). The curve, though complicated, can be fitted reasonably well with a power law of exponent q ≈ 1.4, which is somewhat steeper than the slope of q ≈ 1 measured in the non-ideal simulation of Ross et al. (2016). The exact relation is difficult to obtain, not least because it only appears in the transient phase of these heating simulations, and is thus subject to the initial condition and the simulation’s numerical particulars. Moreover, the natural stochastic and bursty variation in the stress can partially destroy the correlation during the transient phase if one is unlucky. To check on this, we carried out four additional heating runs with identical set ups but different random initial conditions and calculated qs ranging from 0.7 to 1.8 (full details not shown). The average of the set is ≈1.1.

Figure 2.

Log–log plot of total stress against pressure from a pure heating (red) and a pure cooling (blue) run. The slopes (dashed lines) are q ∼1.40 (heating) and q ∼1.50 (cooling).

Open in new tabDownload slide

Next, we performed simulations that cooled monotonically in order to check that the stress follows the pressure not only when the latter increases but also when it decreases. All simulations were started from the box-dominated regime of our fiducial heating run (around orbit 50) and employed the non-linear cooling prescription given by equation (8), with the cooling parameters set to m = 2 and θ = 7.66 × 10⁻³, respectively. An example cooling track is plotted in blue in Fig. 2. For this example, we find a power law exponent q closer to 1.5, before the box settles on a stable fixed point in the thermal dynamics around 〈P〉 = 2.51. This demonstrates that the stress depends on the pressure in a similar way both when the gas heats up and when it cools down.

4.2 Runs with oscillatory cooling

In this section, we investigate how the stress responds to oscillatory variations in the mean pressure, and determine any time lags that appear. To do so, we employ the explicit cooling prescription described in Section 3.1.1. By adjusting the two cooling time-scales τ_H and τ_C, we can control the behaviour of the mean pressure, allowing it to increase when the cooling time-scale τ_H is sufficiently long, so that heating is greater than cooling (the ‘heating phase’), and to decrease when τ_C is sufficiently short, so that cooling overwhelms heating (the ‘cooling phase’). We present three different cases: a simulation in which the mean thermal (i.e. pressure) oscillation is ‘intermediate’, i.e. exhibiting a period of ∼30 – 40 orbits, and thus lying between the dynamical time (∼1/Ω) and thermal time [∼1/(αΩ) ∼ 100/Ω]; ‘fast’ oscillations, with period ∼10 orbits, such that the stress cannot keep up with the changes in the pressure; and ‘slow’ oscillations, with a period of order or greater than the thermal time.

4.2.1 Intermediate thermal oscillations

In Fig. 3, we show the time evolution of stress and pressure in our intermediate case. The cooling in this run was turned on shortly after the MRI’s initial breakdown (orbit 5). The long and short cooling time-scales were set to τ_H = 30 orbits and τ_C = 5 orbits, and the mean pressure was permitted to vary between 〈 P〉 ₋ = 1 and 〈 P〉 ₊ = 2.5. It is plotted in blue. The black curve shows the stress normalized by P₀, and multiplied by 100 to make comparison with the pressure easier. To better track its longer time-scale oscillations, we smooth out the shorter stochastic variation with a lag-free rolling mean (red curve).²

Figure 3.

Time series of total stress and pressure from a simulation exhibiting ‘intermediate’ thermal oscillations. The black curve shows the volume-averaged stress (multiplied by 100), the red curve shows the the lag-less rolling average of the stress, and the blue curve shows the volume-averaged pressure. The dashed vertical line shows the time at which cooling was turned on. The inset shows a time interval between orbit 150 and orbit 200: the lag can clearly be seen by eye.

Open in new tabDownload slide

We begin by analyzing the time series of the mean pressure. After the cooling function is turned on, the pressure continues to increase, though at a slower rate than before. At around orbit 25, the pressure reaches its maximum permitted value at which point the cooling time-scale is switched to the shorter time-scale τ_C. The pressure consequently drops until 〈P〉₋, at which point the cycle repeats. On the other hand, after an initial transient (up to orbit 75), the stress tracks the pressure oscillation relatively closely, but with a noticeable lag of several orbits. This means during the start of each cooling phase the stress hits its peak, and because of the extra heating this provides, and the stress’s inherent stochasticity, the cooling phases lengthen irregularly.

It is important to check whether our runs are in the box-dominated regime or not. In our heating-only run, we found that the simulation entered the box-dominated regime once the stress reached around 0.06–0.07 (cf. the right-hand panel of Fig. 2) after which the stress plateaued. In Fig. 3, Π_xy reaches ∼0.07 only at the very tips of some of the outbursts. Thus we are confident that our oscillatory results are not contaminated by the box size.

To quantify the lag between the stress and the pressure, we calculated the cross correlation between the two, shown in Fig. 4. The correlation was calculated using data between orbits 66 and 350 in order to omit residual and transient behaviour following initialization. We measure a peak correlation coefficient of r = 0.47. This peak occurs when the stress lags behind the pressure by around 4.54 orbits.

Figure 4.

Cross-correlation coefficient between stress and pressure from the run exhibiting ‘intermediate’ thermal oscillations. To the right of the black vertical dashed line is the correlation for a lag in which the stress follows the pressure, while to the left of the line the correlation is for the stress leading the pressure. The peak correlation (red vertical dashed line) occurs when the stress lags behind the pressure by about  + 4.54 orbits.

Open in new tabDownload slide

4.2.2 Fast thermal oscillations

Next, we investigate how the stress tracks the pressure when the latter changes over time-scales that are much shorter than those appearing in the ‘intermediate’ case. To instigate such rapid variations, we require that the pressure increases quickly during the heating-dominated phase (thus raising τ_H), and that the pressure decreases quickly during the cooling-dominated phase (lowering τ_C). Thus we set the cooling time-scales to τ_H = 30 orbits and τ_C = 3 orbits. This simulation is restarted from orbit 80 of our ‘intermediate’ simulation, in order to avoid the initial transient phase, and was run until orbit 200.

In Fig. 5, we show the time-evolution of the total stress multiplied by 100, the rolling-averaged total stress, and the pressure. After the restart position (orbit 80) the mean pressure varies rapidly due to the shorter cooling time-scale τ_C. The stress struggles to track the pressure in this run, and after a few cycles, the stress and pressure appear to be almost completely uncorrelated (e.g. from around orbit 120). This is confirmed by the cross correlation (not shown), which shows that the Pearson r has no clear peak, and instead ‘wiggles’ between small positive and negative values between −15 orbits to +15 orbits.

Figure 5.

Time series of total stress and pressure from simulation exhibiting ‘fast’ thermal oscillations. The black curve shows the volume-averaged stress (multiplied by 100), the red curve shows the lag-less rolling average of the stress, and the blue curve shows the volume-averaged pressure. Note that after the second restart (vertical dashed red line), the stress is unable to keep up with the rapid variation in the mean pressure.

Open in new tabDownload slide

4.2.3 Slow thermal oscillations

Finally, in our ‘slow’ simulation we look at the other extreme: when the thermal driving occurs on a time-scale of order the thermal time, and thus much longer than the dynamical time-scale. This was achieved by lowering the ‘long’ cooling time-scale to τ_H = 16 orbits (resulting in longer rise times of 〈 P〉 during the heating-dominated phase) and increasing the ‘short’ cooling time-scale to τ _C = 5.75 orbits (thus prolonging the decrease of 〈 P〉 during the cooling-dominated phase). The simulation was run for 500 orbits.

Our results appear in Fig. 6. The pressure undergoes an irregular, somewhat bursty, cycle of period ∼170 orbits. The stress exhibits clear features of this thermal oscillation, and thus can be deemed to follow the pressure on this long time-scale. As in previous simulations, the stress undergoes a great deal of large-amplitude stochastic variation on a range of shorter times. In fact, during the long cooling periods between heating bursts, we observe that the pressure experiences associated variations, because it is buffeted around by the stochastic heating. On these shorter times, the pressure possesses a random component that follows the stress, and in fact lags an orbit or so behind.

Figure 6.

Time series of total stress and pressure from a simulation exhibiting ‘slow’ thermal oscillations. The black curve shows the volume-averaged stress (multiplied by 100), the red curve shows the lag-less rolling average of the stress, and the blue curve shows the volume-averaged pressure.

Open in new tabDownload slide

To illustrate these two competing effects, we plot the cross correlation of the two time series in Fig. 7. We find significant correlation in a region encompassing a lag of zero, but the profile is rather blurred. We have a small peak at −1.11 orbits, which we associate with the component of the pressure that follows the stress’s short-time stochastic variability during cooling phases. But this peak joins up to a plateau extending to a lag of about  + 6 orbits, a feature that we associate with the stress following the pressure on the long time-scale of the driven thermal cycle. This ‘blurring’ was not apparent in the intermediate case due to the dominance of the thermal oscillation, and it is likely that if we ran our ‘slow’ cycle run for more periods, we would observe a cleaner peak at a time lag of ≈+5.

Figure 7.

Same as Fig. 4 but the cross correlation is taken from a simulation with ‘slow’ thermal oscillations.

Open in new tabDownload slide

5 CONCLUSION

We have carried out unstratified shearing box zero-net-flux MRI simulations in pluto in order to characterize the temporal relationship between turbulent stress and pressure, with an eye towards the existence, size, and sign of a time lag in the variation of these two fields. There are several astrophysical applications of this work, but we have focused on its relevance for thermal instability and viscous overstability.

Our main set of simulations were run with a time-dependent cooling prescription which could drive long thermal oscillations in the gas, the period of which we could control. These would force oscillations in the pressure, which potentially could force oscillations in the turbulent stress, but with a lag, we could measure. Simulations exhibiting thermal cycles on an intermediate time-scale, with a period ∼50 orbits, show the stress clearly lagging the pressure with a peak cross correlation at about 4.5 orbits. Runs with a faster thermal cycle, periods ∼10 orbits, possess a turbulent stress poorly correlated with the pressure: the stress struggles to ‘keep up’ with the thermal oscillations because the measured time lag is too close to the oscillation period. Finally, when we attempted to drive long cycles the relationship between the stress and pressure is more complicated: while on the longer period of the cycle, we find the stress following the pressure with a lag of ∼5 orbits, as before, on shorter times the stress exhibits considerable stochastic variability, which the stress follows with a lag of around 1 orbit.

The main application of these results is to thermal instability, particularly in low-mass X-ray binaries. Linear theory shows that its growth rate is significantly decreased if the time lag τ is of the order or greater than the thermal time-scale, i.e. if τ ≳ (αΩ)⁻¹ (cf. Section 2.2; Lin et al. 2011; Ciesielski et al. 2012). If our numerical results generalize, and τ ≈ 5 orbits, then thermal instability is only impeded appreciably for large α ≳ 0.1, which should indeed be the case in radiation-pressure dominated states in outburst (Lasota 2001). A precise quantitative estimate for the reduction in the thermal instability’s growth rate is not possible to derive from our simulations. It is unlikely, however, that thermal instability is completely quenched by the lag; yet it could combine with other physics, such as the equilibrium’s inherent stochasticity and magnetic fields (Begelman & Pringle 2007; Oda et al. 2009; Ross et al. 2017), to delay the onset of thermal instability to the point that it is no longer dynamically relevant. An important caveat when applying our results is that our simulations were run with gas pressure only and, being zero-net flux, supported relatively low-level MRI turbulence, with α ∼ 0.01. Future work should attempt to generalize our calculations to conditions more representative of the inner radii of low-mass X-ray binaries in outburst.

Finally, we note that the response of the stress to the pressure is only one aspect of the turbulent dynamics. More generally, the turbulent stress will depend on strain and density as well as pressure, and may relax towards a Navier–Stokes description on some characteristic time-scale (e.g. Ogilvie 2001, 2003). It is likely that the time lags measured in our simulations correspond to this more general relaxation time. If such an identification can be made, then our results also bear on the onset of viscous overstability, and in fact indicate that the time-lags we find are sufficient to stabilise the disc.

ACKNOWLEDGEMENTS

This work was funded by a Science and Technologies Facilities Council (STFC) studentship. Simulations were run on the Fornax cluster at the Department of Applied Mathematics and Theoretical Physics (DAMTP) in Cambridge, and on the Sakura and Cobra clusters at the Max Planck Computing and Data Facility (MPCDF) in Garching. We would like to thank the referee for comments which helped us clarify certain concepts in the paper.

DATA AVAILABILITY

The data underlying this article will be shared on a reasonable request to the corresponding author.

Footnotes

REFERENCES