Demonstration of a magnetic Prandtl number disc instability from first principles

TL;DR

This paper demonstrates from first principles that the magnetic Prandtl number can induce a thermal-viscous instability in accretion discs, potentially explaining state transitions in black hole binaries.

Contribution

First self-consistent MHD simulation showing that the magnetic Prandtl number dependence causes a thermal-viscous instability in accretion discs.

Findings

The $oldsymbol{ m extit{ extbf{α}}}$-Pm dependence persists with increased resolution.

The Pm dependence can lead to a local thermal-viscous instability.

The instability manifests as an unstable limit cycle in the disc.

Abstract

Understanding what determines the strength of MHD turbulence in accretion discs is a question of fundamental theoretical and observational importance. In this work we investigate whether the dependence of the turbulent accretion disc stress ($\alpha$) on the magnetic Prandtl number (Pm) is sufficiently sensitive to induce thermal-viscous instability using 3D MHD simulations. We first investigate whether the $\alpha$-Pm dependence, found by many previous authors, has a physical or numerical origin by conducting a suite of local shearing-box simulations. We find that a definite $\alpha$-Pm dependence persists when simultaneously increasing numerical resolution and decreasing the absolute values of both the viscous and resistive dissipation coefficients. This points to a physical origin of the $\alpha$-Pm dependence. Using a further set of simulations which include realistic turbulent heating and radiative cooling, and by giving Pm a realistic physical dependence on the plasma temperature and density, we demonstrate that the $\alpha$-Pm dependence is sufficiently strong to lead to a local instability. We confirm that the instability manifests itself as an unstable limit cycle by mapping the local thermal-equilibrium curve of the disc. This is the first self-consistent MHD simulation demonstrating the Pm instability from first principles. This result is important because a physical Pm instability would lead to the global propagation of heating and cooling fronts and a transition between disc states on timescales compatible with the observed hard/soft state transitions in black hole binaries.

Tables

Table 1. Table 1: Table showing the simulation parameters used in this work, where β x = ( c s / v A ​ x ) 2 subscript 𝛽 𝑥 superscript subscript 𝑐 𝑠 subscript 𝑣 𝐴 𝑥 2 \beta_{x}=(c_{s}/v_{A\,x})^{2} is the plasma beta of the initial net magnetic field in the x 𝑥 x -direction and κ 𝜅 \kappa is the electron scattering opacity used in the cooling function in cgs units.

Simulation set	${\beta }_{\varphi }$	${\beta }_z$	$\mathrm{Pm}$	Rm	Resolution (x,y,z)	E.O.S	$\rho /{\rho }_0$	$\kappa$	Pminitial
IsoBphi1	$100$	n.a.	1/2, 1, 2, 4	variable	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsoBphi2	$100$	n.a.	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsoBphi3	1000	n.a.	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsoBphi4	10000	n.a.	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsoBphi5	10000	n.a.	1/16, 1/4, 1, 4, 16	12800	$128\times 200\times 128$	isothermal	n.a.	n.a.	n.a.
IsoBphi6	10000	n.a.	1/16, 1/4, 1, 4, 16	25600	$128\times 200\times 128$	isothermal	n.a.	n.a.	n.a.
IsoBphi7	10000	n.a.	1/16, 1/4, 1, 4, 16	51200	$256\times 400\times 256$	isothermal	n.a.	n.a.	n.a.
IsoBz1	n.a.	$100$	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsoBz2	n.a.	1000	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsoBz3	n.a.	10000	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IsononetB	n.a.	n.a.	1/16, 1/4, 1, 4, 16	12800	$64\times 100\times 64$	isothermal	n.a.	n.a.	n.a.
IdealBphi1	10000	n.a.	variable (1/16-16)	25600	$128\times 200\times 128$	ideal	0.6$\dot{7}$, 0.7$\dot{3}$, 0.8, 0.9$\dot{3}$, 1.0, 1.0$\dot{6}$	0.4	16
IdealBphi2	10000	n.a.	variable (1/16-16)	25600	$128\times 200\times 128$	ideal	0.6$\dot{7}$, 0.7$\dot{3}$, 0.8, 0.9$\dot{3}$, 1.0, 1.0$\dot{6}$	0.4	1/16

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla \cdot (\rho \mathbf{vv}-\mathbf{BB})+\nabla P_{\mathrm{tot}}=\rho {\mathbf{g}}_s-2{\Omega }_0\hat{\mathbf{z}}\times \rho \mathbf{v}+\nabla \cdot \Pi ,$$

$$\frac{\partial E}{\partial t}+\nabla \cdot [(E+P_{\mathrm{tot}})\mathbf{v}-(\mathbf{v}\cdot \mathbf{B})\mathbf{B}]$$

$$\frac{\partial \mathbf{B}}{\partial t}-\nabla \times (\mathbf{v}\times \mathbf{B})=-\nabla \times (\eta \cdot \mathbf{J}),$$

$$\mathbf{J}=\nabla \times \mathbf{B},P_{\mathrm{tot}}=P+\frac{B^2}{2},$$

$$E=\rho e+\frac{\rho v^2}{2}+\frac{B^2}{2},{\mathbf{g}}_s={\Omega }_0^2(2qx\hat{\mathbf{x}}-z\hat{\mathbf{z}}),$$

$${\Pi }_{ij}=\nu \left[\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}-\frac{2}{3}\nabla \cdot \mathbf{v}{\delta }_{ij}\right].$$

$$\Lambda =\frac{\sigma T_e^4}{H},$$

$$T_e^4=\frac{4T^4}{3\Sigma \kappa }.$$

$$\frac{\partial \ln \alpha }{\partial \ln \Sigma }+\frac{1}{4}\frac{\partial \ln \tau }{\partial \ln \Sigma }+1<0,$$

$$\alpha P_{\mathrm{tot}}=⟨\rho (\delta v_r\delta v_{\varphi }-v_{Ar}{v}_{A\varphi })⟩,$$

$$\eta =\frac{5.6\times 10^{11}\ln {\Lambda }_{eH}}{T^{3/2}}{\mathrm{cm}}^2{\mathrm{s}}^{-1},{\nu }_{\mathrm{Rad}}=\frac{6.7\times 10^{-26}{T}^4}{\kappa {\rho }^2}{\mathrm{cm}}^2{\mathrm{s}}^{-1},$$

$$\mathrm{Pm}\simeq \frac{{\nu }_{\mathrm{Rad}}}{\eta }=1.9\times 10^{-38}\frac{T^{11/2}}{\kappa {\rho }^2},$$

$$\mathrm{Rm}=\frac{c_{s}H}{\eta },\mathrm{Re}=\frac{c_{s}H}{\nu },\mathrm{Pm}=\frac{\mathrm{Rm}}{\mathrm{Re}}.$$

$$\eta =\frac{c_{s0}{H}_0}{25600},\nu =\frac{c_{s0}{H}_0}{25600}{\left(\frac{T}{T_0}\right)}^{11/2}.{\left(\frac{\rho }{{\rho }_0}\right)}^{-2},$$

$$\frac{\sigma T_e^4}{H}=\frac{3}{2}\Omega \alpha (\mathrm{Pm})P_{\mathrm{tot}},$$

Full text

††pagerange: Demonstration of a magnetic Prandtl number disc instability from first principles–Demonstration of a magnetic Prandtl number disc instability from first principles††pubyear: 2017

1 Introduction

Accretion discs surrounding compact objects exhibit complex and dramatic cyclic changes in X-ray luminosity on a variety of timescales (e.g. Fender et al. 2004, Remillard & McClintock 2006 and Done et al. 2007). Black hole binaries (hereafter “BHB”) are observed to undergo flaring and quiescent cycles on month-long timescales in which the X-ray luminosity of the disc changes by orders of magnitude. This long-timescale cyclic flaring is thought to be caused by a thermal-viscous instability due to the sharp increase in disc opacity when hydrogen becomes ionised. This is known as the disc instability mechanism (DIM) (Faulkner et al. 1983, Hameury et al. 1998 and Lasota 2001). In the last decade, cyclic changes in the disc luminosity and X-ray spectrum have been observed in which the disc cycles between a radiatively efficient high soft state, resembling a classic composite blackbody and a radiatively inefficient low hard state, with a power law non-thermal hard X-ray spectrum. These changes in disc state are of particular interest because they are intimately linked to the production of a relativistic jet in the system; in the hard state a radio jet is observed, whilst in the soft state no jet is observed. However, the physical mechanism responsible for these state changes is not yet understood. In this paper we present simulations of a promising new thermal-viscous instability which could explain the hard/soft changes in disc state.

In a previous paper (Potter & Balbus, 2014) (hereafter PB), we put forth arguments for the existence of a new type of thermal-viscous disc instability triggered when the disc stress (usually parameterised by $\alpha$, Shakura & Sunyaev 1973) depends sensitively on the physical properties of the disc plasma, especially temperature (see also the instability calculations in Takahashi & Masada 2011). Simulations both of accretion disc turbulence (Fromang et al. 2007, Lesur & Longaretti 2007, Simon & Hawley 2009, Meheut et al. 2015), and of driven turbulent dynamos (Schekochihin et al. 2004) had shown that when it is near unity, the magnetic Prandtl number (the ratio of the microscopic viscosity to resistivity, ${\rm Pm}=\nu/\eta$), can have a significant effect on the strength and the maintenance of turbulent fluctuations. A possible explanation for the dependence of the disc $\alpha$ parameter on ${\rm Pm}$ is that the rate of small scale magnetic reconnection in turbulent flow is sensitive to the strength of the viscous stress in a reconnecting layer when Pm is near unity (Balbus & Henri, 2008).

Magnetic and kinetic energy are extracted from the differential rotation of the disc by the magnetorotational instability (MRI), predominantly on spatial scales set by the disc scale height Balbus & Hawley (1998). Once the MRI is fully non-linear, the gas hosts a turbulent cascade, transferring energy downwards to ever smaller length scales, until finally it is dissipated as heat. In the ${\rm Pm}>1$ regime, the viscosity is larger than the resistivity on a given scale, and so the viscous dissipation is correspondingly larger than resistive dissipation; the opposite is true when ${\rm Pm}<1$. Both the growth rate of the MRI and the magnitude of the turbulent disc stress (most of which is magnetic) are related to the RMS magnetic field strength of the disc. The saturated magnetic field strength in turn depends upon the rate of magnetic reconnection in the plasma. When ${\rm Pm}>1$, the viscous dissipation scale exceeds the resistive scale and so the viscosity damps out velocity fluctuations on the resistive scale. This leads to a lower rate of reconnection via dissipation, a build-up of the magnetic energy, and an increased value of $\alpha$. This is because the field dissipation process will generally require large velocity gradients over the resistive lengthscale in order to bring misaligned magnetic field lines together. If these large gradients produce correspondingly large viscous stresses, the dissipation will be inhibited. If on the other hand ${\rm Pm}<1$, the resistive length scale is larger than the viscous length scale, and velocity gradients on this resistive scale will not produce important dynamical stress. Magnetic dissipation unencumbered by viscosity will proceed apace, the cascade decreasing the saturated magnetic field strength, and with it the $\alpha$ stress.

The disc radius at which the transition between ${\rm Pm}>1$ and ${\rm Pm}<1$ occurs is $r\sim 100r_{\mathrm{s}}$ in a typical BHB accretion disc (PB). The dependence of $\alpha$ on ${\rm Pm}$ is typically a power law, $\alpha\propto\mathrm{Pm}^{n}$. PB showed that within the classic $\alpha$ formalism, when $n>0.5$, the disc is unstable. Many (though not all) MHD simulations of the dependence of $\alpha$ on ${\rm Pm}$ had found $n\sim 0.5-1$, suggesting that the instability might plausibly be present in astrophysical discs (Fromang et al. 2007 Lesur & Longaretti 2007 Simon & Hawley 2009 etc.). Using an idealised 1D dynamic thin disc approximation, PB demonstrated that the instability does indeed lead to the formation of a local unstable limit cycle. Dynamical heating and cooling fronts move throughout the disc, as in classical dwarf nova modelling (e.g. Hameury et al. 1998).

These preliminary results suggest two further avenues of exploration, which we follow here. The first is to investigate with some care the dependence of the disc stress upon the magnetic Prandtl number in a variety of initial conditions using 3D MHD isothermal shearing-box simulations. It is obviously critical to establish that the ${\rm Pm}$ dependence is physical, and not a numerical artefact. The second challenge is to conduct simulations of the disc thermal equilibrium curve using more realistic cooling, including a temperature and density dependent Pm. Does the suspected instability actually occur in a first-principle MHD simulation?

2 Numerical setup

We carry out this investigation with the PLUTO MHD code (Mignone et al., 2007), which is widely-used and publicly available. The dissipative MHD equations are solved in the local shearing-box approximation using a Godunov-type finite volume scheme with explicit dissipation terms (Bodo et al. 2008, Hawley et al. 1995, Mignone et al. 2012, Bodo et al. 2014). In conservation form, the equations are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho$ is the mass density, ${\bf v}$ the fluid velocity, ${\bf B}$ the magnetic field, $P_{\mathrm{tot}}$ the total pressure (thermal plus magnetic), $P$ the thermal pressure, ${\bf g}_{s}$ the effective gravity in the shearing-box approximation, $\Omega_{0}$ is the Keplerian angular velocity at the centre of the shearing-box, $\Pi$ the viscous stress tensor, $E$ the energy density, $e$ the internal energy per unit mass, $\eta$ the resistivity, $\nu$ the shear viscosity (we assume $\eta$ and $\nu$ are diagonal and isotropic tensors, i.e. scalars), ${\bf J}$ the electric current density and $\Lambda$ the energy loss per unit volume due to radiative cooling.

2.1 Cooling

We implement an effective cooling function in the simulation assuming that the opacity is dominated by electron scattering. This is likely to be the case in the inner regions of the disc, where we expect the ${\rm Pm}=1$ transition to occur (PB). Following Faulkner et al. 1983 and Latter & Papaloizou 2012, we use a bulk cooling function to stand in for the diffusive cooling of the form

[TABLE]

where $H$ is the disc scale height $H=c_{s}/\Omega$, $c_{s}$ is the thermal sound speed i.e. $P=\rho c_{s}^{2}$, $\Omega$ is the Keplerian angular velocity and $T_{e}$ is the effective surface temperature given by

[TABLE]

$T$ is the temperature of the simulated central disc plasma, $\kappa=0.4\mathrm{cm}^{2}\mathrm{g}^{-1}$ is the electron scattering opacity and $\Sigma=\rho H$ is the disc surface density. The cooling is then calculated using the temperature and density in each simulation cell.

2.2 Instability criterion and ${\rm Pm}$

In PB we derived a thermal-viscous instability criterion for a standard, radiatively efficient thin disc in which the $\alpha$ parameter is a variable depending on $\rho$ and $T$. This is

[TABLE]

where $\tau=\Sigma\kappa$ is the optical depth of the disc. Assuming a power-law dependence for $\alpha$ upon the magnetic Prandtl number, $\alpha\propto{\rm Pm}^{n}$, a parametisation suggested by simulations (Lesur & Longaretti 2007, Simon & Hawley 2009), the disc is unstable if $n>0.5$ (PB, equations 46 and 47). Here $\alpha$ is determined by the usual formula

[TABLE]

where $\delta v_{x}$ is residual velocity after subtraction of the local Keplerian circular velocity (i.e. $\delta v_{\phi}=v_{\phi}-\Omega R$, $\delta v_{r}=v_{r}$), $v_{A\,x}$ is the Alfvén velocity in the $x$-direction, $v^{2}_{A\,x}=B_{x}^{2}$, and the angle brackets denote a spatial average. For a typical BHB accretion disc environment, the radiative viscosity dominates the Coulomb viscosity, and one finds:

[TABLE]

[TABLE]

where $\kappa$ is the opacity of the plasma and $\ln\Lambda_{eH}$ is the electron-proton Coulomb collision factor, which we take to have a value of $\sqrt{40}\simeq 6.3$ (see PB section 2.1.3 for more detail).

3 How does $\alpha$ depend on dissipation and resolution?

In this section we wish to answer two questions, (i) to what extent do the viscous and resistive dissipation coefficients determine the strength of the saturated disc turbulence quantified by $\alpha$; and (ii) to what extent are these effects physical or numerical in origin? Recent simulations from several groups show that $\alpha$ depends on a whole array of physical and numerical parameters, such as the height of the simulating box, stratification, the presence of net magnetic fields, numerical resolution and convection, etc. (e.g., Simon et al. 2012, Hawley et al. 1995, Fromang et al. 2007, Hirose et al. 2014 and Ryan et al. 2017). Here we are interested in isolating the effect of dissipation coefficients. To avoid conflating this with other complications, we initially choose the simplest isothermal, unstratified local shearing-box. We wish to study the effect of dissipation coefficients, numerical resolution and the initial magnetic field configuration. In these isothermal simulations we maintain fixed dissipation coefficients throughout an individual run. (We shall later allow the viscosity and resistivity to become time-dependent functions of temperature and density.)

In Figures 1-3 and Table 1 we summarise the results of the isothermal simulations for a variety of values of resistivity, viscosity, net magnetic field and resolution. The resistivity and viscosity can be expressed in terms of the dimensionless Reynolds number, $\rm{Re}$, and magnetic Reynolds number Rm, given by.

[TABLE]

3.1 Initial conditions

Explicit parameters of our numerical simulations are provided in Table 1. All simulations start with a zero net field in the $z$-direction, $B_{z}=B_{0}\sin(2\pi x/H)$ and plasma beta, $\beta_{0}=100$, where $\beta_{0}=2P/B_{0}^{2}$ is the ratio of thermal to magnetic pressure. We shall subsequently quantify the strength of any net magnetic field as $\beta_{x}=2P/B^{2}_{x}$, the ratio of thermal pressure to the magnetic pressure of the net field in the $x$-direction. The dimensionless orbital frequency and thermal sound speed are chosen to be $\Omega=c_{s}=0.001$, so $H=1$, and box lengths are $L_{x}=H$, $L_{y}=4H$ and $L_{z}=H$. Simulations were initialised with random pressure perturbations of maximum amplitude $1\%$, as in Simon & Hawley (2009).

We first investigate disc turbulence for a variety of values of Rm and Re. Our findings are in basic agreement with previous investigations, which found that larger values of viscosity increase the value of the $\alpha$, whilst large resistivities decrease $\alpha$ (Fig. 1). There is clearly a complex interaction between the viscosity, resistivity, net B-field and resolution. Large values of the viscosity and resistivity, $\rm{Re}<800$, $\rm{Rm}<800$, result in sufficiently high dissipation on large lengthscales to prevent sustainable turbulence. These values are not representative of what would be expected in realistic astrophysical disc environments, and we avoid this regime. On the other hand, excessively small values of the dissipation coefficients (i.e. large Reynolds numbers) will produce no measurable effect in the simulations, since the relevant lengthscales on which the viscosity and resistivity become important will be below the grid resolution (and unavoidable numerical dissipation) of the simulations. Moreover, whilst an excessively large viscosity reduces the linear growth rate of the MRI, a large resistivity eliminates the linear MRI. We focus therefore on the intermediate asymptotic regime in which the resistivity is held constant at the lowest resolvable value, and allow only the viscosity to vary in order to study the effect of altering ${\rm Pm}$.

The results of varying ${\rm Pm}$ (with fixed Rm) are shown in Figures 2 and 3 and Table 1. The results, as expected, show an increase in the strength of the saturated MRI turbulence for values of ${\rm Pm}>1$. At small values of ${\rm Pm}$, $\alpha$ tends to a constant value determined by the initial net B-field in the simulation, $\rm{Rm}$, and the resolution of the simulation. The presence of an imposed initial net magnetic field is expected to increase the strength of saturated disc turbulence because it provides a sort of ‘backbone’ magnetic field which cannot be dissipated and is thus a permanent source of magnetic field. The initial net magnetic field is crucial in determining the minimum value of $\alpha$ at small ${\rm Pm}$. Net $B_{\phi}$ fields have a strong impact on the $\alpha$ value at ${\rm Pm}\lesssim 1$, with larger $B_{\phi}$ fields increasing the value of $\alpha$. Net $B_{z}$ fields increase $\alpha$ at all values of ${\rm Pm}$ (the strongly enhancing effect of a net $B_{z}$ field has been known for some time Hawley et al. 1995).

The size of the dissipation coefficients is of necessity unphysically large in simulations (due to a lack of dynamic range). To establish the dependence of $\alpha$ on ${\rm Pm}$ as a physical, rather than a numerical, effect, we need to show that variations in $\alpha$ exist which depend only upon the ratio of dissipation coefficients (i.e. ${\rm Pm}$) and not on their absolute values. This is done by comparing the results from simulations covering the same range of ${\rm Pm}$ with the same initial conditions, but as the numerical resolution of the simulation is increased, the absolute values of viscosity and resistivity are simultaneously decreased. Results are shown in Figure 3 where it can be seen that using the same initial setup, the simulations at higher numerical resolution (which at a given value of ${\rm Pm}$ have substantially decreased values of both viscous and resistive dissipation coefficients) obtain remarkably similar results to the lower resolution simulations. If the ${\rm Pm}$-$\alpha$ effect were due to artificially large dissipation lengthscales and insufficient scale separation, we would expect the effect to decrease with increased resolution and decreased values of the dissipation coefficients. Figure 3 shows this is clearly not the case. In fact, the ${\rm Pm}$-$\alpha$ effect is slightly more pronounced for ${\rm Pm}<1$ at higher resolution, even when the absolute dissipation coefficients are smaller. This is a clear indication that whatever residual numerical effects may be present, there appears to be a distinct and genuine physical ${\rm Pm}$ effect at work (Fig. 4).

The formal instability criterion is $\partial\ln\alpha/\partial\ln{\rm Pm}>0.5$ (PB). From Figure 3 this is satisfied in the ${\rm Pm}>1$ regime for the high resolution simulations with small mean magnetic field strengths.

4 Thermal equilibrium curve

Let us now investigate whether the dependence of $\alpha$ on ${\rm Pm}$ leads to the thermal-viscous instability outlined in PB. To test this, the local thermal-equilibrium curve of the disc is simulated. It is necessary to include both turbulent heating and radiative cooling (the gas is no longer isothermal), as well as a variable magnetic Prandtl number which has the correct dependence on temperature and density calculated in (13). To minimise the direct effect of a large resistivity on the linear MRI growth rate, we adopt a constant low value for the resistivity ($\rm{Rm}\simeq 25600$) and include the temperature and density dependence of ${\rm Pm}$ in the viscosity alone. The explicit resistivity and viscosity are given by

[TABLE]

where in the formula above $c_{s0}$ and $H_{0}$ are the initial values of the thermal sound speed and disc scale height. $T_{0}$ and $\rho_{0}$ are normalisation parameters which are chosen for convenience to be of order unity in Figure 5. The range of allowed ${\rm Pm}$ values is limited to $1/16>{\rm Pm}>16$ (or $409600>\rm{Re}>1600$). This range is chosen because when ${\rm Pm}<1/16$ the viscosity is completely sub-dominant to the numerical grid viscosity and when ${\rm Pm}>16$ the viscosity becomes unrealistically large, with the viscous lengthscale becoming comparable to the total size of the simulation. Instability, if present, will manifest itself as two sets of overlapping quasi-stable solutions (stable on the thermal timescale but evolving on the mass accretion timescale, hence the term thermal-viscous instability). To find such solutions, we carry out two simulations with nearly identical initial conditions, the sole difference being that one simulation starts with a fixed high ${\rm Pm}=16$ and one with a fixed low ${\rm Pm}=1/16$. The MRI is allowed to grow, become fully turbulent, and attain a dynamical-thermal equilibrium in the two cases. From Fig. 2 the high ${\rm Pm}$ simulation is expected to have stronger turbulent heating due to the larger $\alpha$ value, and so to reach a higher equilibrium temperature than the low ${\rm Pm}$ simulation. After equilibrium has been achieved, ${\rm Pm}$ is freed and its value is calculated via equation 15 for each grid cell. Thus ${\rm Pm}$ becomes a function of the local fluid temperature and density. The two simulations are then evolved until they relax to a new quasi-stable thermal equilibrium, if this in fact differs from the initial equilibrium. If the instability is present in the simulation, this should manifest itself as two stable thermal equlibria at the same surface density $\Sigma$, producing the characteristic unstable S-curve, as in the case of the hydrogen ionisation instability in dwarf novae (e.g. Lasota 2001, Frank et al. 2002 and PB).

Figure 5 shows the local thermal equilibrium curve of the instability calculated using a set of simulations with a weak net magnetic field in the azimuthal direction ($\beta_{\phi}=10^{4}$). It shows that two stable thermal equilibria exist at the same density in the unstable range of ${\rm Pm}$ (${\rm Pm}\approx 1-16$). This is the first self-consistent MHD simulation demonstrating the ${\rm Pm}$ instability from first principles. It shows that the dependence of $\alpha$ on ${\rm Pm}$ in a full MHD simulation is sufficiently sensitive to induce a thermal-viscous instability in the disc plasma. The black curve in the figure is the expected theoretical thermal equilibrium curve calculated by balancing turbulent heating and radiative cooling. This is in close agreement with the results of the simulation. The black curve is calculated from the thermal equilibrium equation

[TABLE]

where the LHS corresponds to radiative cooling and the RHS to turbulent heating (see e.g. Balbus & Hawley 1998). The value of alpha is determined as a function of ${\rm Pm}$ for the black curve by using the $\alpha$-${\rm Pm}$ dependence found in simulation IsoBphi6, which is the equivalent isothermal simulation. Due to the finite resolution of the simulation, turbulent fluctuations in the spatially averaged values of $\alpha$ and $T$ are much larger than would be expected in an actual disc (in part, because of Poisson noise from the finite number of grid points in the simulation) and so this precludes accurately simulating the solutions close to the edges or corners of the solid S-curve. Close to the edges of the upper and lower branches of solutions small temporal fluctuations in temperature are sufficient to move the average temperature of the plasma over the dashed line of solutions which are unstable to perturbations. The dashed line represents the watershed between states which will experience runaway heating (above the dashed line) or cooling (below the dashed line) until they reach an equilibrium solution on the solid curve. Examples of the evolution of thermally stable and unstable simulations are shown in Figure 6. The precise properties of the thermal-equilibrium curve are clearly dependent on the initial net magnetic field, which is likely to vary between different astrophysical systems.

5 Conclusion

In this paper we address two important questions - (i) to what extent does the turbulent disc stress depend on the magnetic Prandtl number, and is this effect numerical or physical? (ii) is the $\alpha$-${\rm Pm}$ dependence sufficiently sensitive to trigger a thermal-viscous instability in the disc?

In the first section of the paper we address (i) by conducting a suite of isothermal 3D MHD local shearing-box simulations to investigate the dependence of the turbulent disc stress $\alpha$ on the viscous and resistive dissipation coefficients. In agreement with previous studies we find that $\alpha$ depends on the magnetic Prandtl number (${\rm Pm}$), the ratio of viscosity to resistivity. In the regime $1<{\rm Pm}<16$, $\alpha$ increases with ${\rm Pm}$ and the dependence is sufficiently sensitive to induce an important new thermal-viscous instability, the magnetic Prandtl number disc instability (see PB). We investigate whether the $\alpha$-${\rm Pm}$ dependence is physical or numerical in origin by conducting a suite of simulations covering the same range in ${\rm Pm}$, but at different numerical resolutions and with different absolute values of viscosity and resistivity. If the $\alpha$-${\rm Pm}$ dependence were numerical and caused by artificially large dissipation coefficients or a lack of scale separation in the turbulent cascade, the dependence should decrease as the simulation resolution is increased and the dissipation coefficients are decreased. In fact, the $\alpha$-${\rm Pm}$ dependence persists as the simulation resolution is increased and the absolute values of the dissipation coefficients are decreased. This shows that the ${\rm Pm}$-$\alpha$ effect is not a numerical artefact caused by artificially large dissipation coefficients or a lack of scale separation. It is firm evidence that the $\alpha$-${\rm Pm}$ dependence is physical in origin.

To investigate whether the $\alpha$-${\rm Pm}$ dependence is sufficiently sensitive to induce the magnetic Prandtl number disc instability we conducted a further set of MHD simulations which include realistic turbulent heating and radiative cooling to map out the local thermal-equilibrium curve of the disc. In these simulations ${\rm Pm}$ is no longer held fixed; instead it is given a realistic physical dependence on the plasma density and temperature. Significantly, we find that the $\alpha$-${\rm Pm}$ dependence is sufficiently strong to trigger the ${\rm Pm}$ instability using physical conditions appropriate for an accretion disc. The thermal-equilibrium curve is found to have a characteristic S-shape forming an unstable limit cycle. In the unstable region corresponding to the regime $1<{\rm Pm}<16$, the instability manifests itself as two sets of overlapping stable solutions with different equilibrium temperatures but the same surface density. The two sets of stable solutions correspond to a hotter, high $\alpha$ branch of solutions with large ${\rm Pm}$ and a cooler, low $\alpha$ branch of solutions corresponding to ${\rm Pm}<1$ as predicted in PB. This is the first self-consistent MHD simulation demonstrating the existence of the ${\rm Pm}$ disc instability from first principles. It was shown in PB that the local ${\rm Pm}$ instability leads to the production of global heating and cooling fronts propagating through the disc, resulting in changes to the disc state on timescales substantially shorter than those of the DIM. This makes the magnetic Prandtl number disc instability a good potential candidate to explain the hard/soft changes in disc state observed on week long timescales. A detailed comparison of the expected timescale and spectral properties of the global instability will be the subject of future work.

6 Acknowledgements

WJP is supported by a Junior Research Fellowship from University College, University of Oxford. SAB acknowledges support from the Royal Society in the form of a Wolfson Research Merit Award. We acknowledge support from STFC in the form of a Consolidated Grant to Oxford Astrophysics. The simulations presented here were run on the DiRAC Complexity Cluster at Leicester and on Berg, the DiRAC facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Oxford. WJP is particularly grateful to Julien Faure for important numerical advice and to Sebastien Fromang for hosting a visit to the CEA Saclay. We would like to thank Sebastien Fromang, Julien Faure, Charles Gammie, Rob Fender, John Miller, Jean-Pierre Lasota and Alexander Schekochihin for helpful conversations and suggestions.