Three-dimensional, global, radiative GRMHD simulations of a thermally unstable disc

In order to make connections with previous shearing box (esp. Jiang, Stone & Davis, 2013) and global, pseudo-Newtonian thin-disc (Reynolds & Miller, 2009) simulations, instead of starting from the Shakura-Sunyaev solution, we initialized a slab of gas of uniform thickness, orbiting everywhere at the local Keplerian frequency. We chose a black hole mass of $M_{\mathrm{BH}}=6.62M_{\odot}$, for close comparison with Jiang, Stone & Davis (2013). The initial density profile of the azimuthally symmetric disc is

where $\rho_{0}$ is the mid-plane density, $H$ is the initial height of the disc, and $R_{i}$ is the inner radius of the initial disc. We fix the initial inner radius to that of the innermost stable circular orbit (ISCO) and use an exponential cutoff to smooth the transition there. The initial disc structure is thus entirely governed by the chosen mid-plane density and disc height. We consider two cases, one with mid-plane density $\rho_{0}=10^{-3}~\mathrm{g}\,\mathrm{cm}^{-3}$ and height $H_{0}=0.4\,r_{g}$, where $r_{g}=GM/c^{2}$ is the gravitational radius, and the other with $\rho_{0}=10^{-6}~\mathrm{g}\,\mathrm{cm}^{-3}$ and $H_{0}=0.3\,r_{g}$. The first case matches the RSVET model in the shearing-box simulations of Jiang, Stone & Davis (2013). The Jiang, Stone & Davis (2013) simulations were done for a section of disc centered at $R=30\,r_{g}$, whereas we are performing global simulations, so the correspondence is imperfect. In our case, the disc is initially radiation pressure dominated, with $10\lesssim P_{\mathrm{rad}}/P_{\mathrm{gas}}\lesssim 1000$. We refer to this setup as RADPHR or RADPLR depending on the resolution used (see Table 1).

For RADPHR and RADPLR, we take the gas and radiation to be initially in local thermodynamic equilibrium (LTE), which means the total pressure is distributed between the two (magnetic pressure is initially negligible), as

where the total pressure, $P_{\mathrm{tot}}$, is initially taken to be

as in Reynolds & Miller (2009). The disc is thus initially in approximate hydrostatic equilibrium: on the one hand, the pressure value given by equation (3) exceeds the value required for vertical hydrostatic equilibrium in the Schwarzschild metric by a factor of $\approx 1.3$; on the other, until MRI develops, there is no heating in the disc. Substituting the expressions for radiation pressure (in local thermodynamic equilibrium) and ideal gas pressure into equation (2), we have

where $a_{\mathrm{R}}=4\sigma/c$ is the radiation constant. This is a quartic equation for the gas temperature with four possible roots, though only one is positive and real.

For the chosen values of the initial density and disc height, these simulations do not begin in thermal equilibrium. Fig. 1 shows a thermal equilibrium curve at $R=10\,r_{g}$ for an $\alpha=0.02$ thin disc. Our simulations start between the unstable, radiation-pressure-dominated and stable, gas-pressure-dominated branches.

The second case we consider is somewhat different in terms of its thermodynamic properties. This lower density disc is initially gas pressure dominated, with $10^{-11}\lesssim P_{\mathrm{rad}}/P_{\mathrm{gas}}\lesssim 10^{-7}$ for entire disc. We refer to this setup as GASPLR or GASPHR, again depending on the resolution (Table 1). GASPHR and GASPLR have an initial disc thickness $H_{0}=0.3r_{g}$. In this case, we do not assume LTE, initially. Instead, we define a uniform, initial radiation temperature, $T_{\mathrm{rad}}=10^{4}$ K, considerably below the initial gas temperature. We do this to ensure that the simulation starts out completely gas-pressure-dominated. In effect, we are artificially suppressing the initial radiation field. However, the surface density in this case is so low ($\leq 0.735~\mathrm{g~cm^{-2}}$) that the disc is effectively optically thin to electron scattering, so the coupling between the gas and radiation is weak.

The $\mathbf{M}_{1}$ closure scheme in Cosmos++ (Fragile, Olejar & Anninos, 2014) evolves the radiation energy density in the radiation rest frame, $E_{R}$, and the spatial components of the radiation rest frame 4-velocity, $u^{i}_{R}$ (see Section 2.5). However, it is easier to initialize our simulations by defining the radiation fields in the fluid frame. We already described in the previous section how we determine the initial radiation temperature, $T_{\mathrm{rad}}$. This temperature can be used to define the radiation energy density in the fluid frame

The initial radiative flux, $F^{i}$, is then calculated from the gradient of $E_{\mathrm{rad}}$. With these quantities, we can construct the contravariant radiation stress-energy tensor

where $h^{\alpha\beta}=g^{\alpha\beta}+u^{\alpha}u^{\beta}$ is the projection tensor and $u^{\alpha}$ is the fluid 4-velocity.

We can equate the radiation stress-energy tensor from equation (6) to the form written in terms of the radiation rest frame variables:

Following Sa̧dowski et al. (2013) and Fragile, Olejar & Anninos (2014), we start with the following two relationships, both of which come from equation (7):

and

Using these, we solve for the radiation energy density in the radiation rest frame, $E_{R}$, and the time component of the radiation rest frame four-velocity, $u^{t}_{R}$. With these, the remaining spatial components of the radiation rest frame 4-velocity, $u^{i}_{R}$, can easily be obtained from equation (7).

To seed the magneto-rotational instability (MRI) inside the disc, we impose a weak magnetic field ($\beta=P_{\mathrm{tot}}/P_{\mathrm{mag}}\geq 10$) on top of our hydrodynamical setup. The MRI is necessary to drive the accretion of matter into the black hole and transport angular momentum outwards (Balbus & Hawley, 1991). The turbulence of the MRI will also heat the disc, which is important to our goal of studying thermal stability.

It is important for the initialized magnetic field to be divergence free. The easiest way to accomplish this is to initialize the magnetic field starting from the vector potential. For our thin discs, we set $A_{r}=A_{z}=0$ and

where

and $R_{\mathrm{ISCO}}=6\,r_{g}$ is the radius of the ISCO. The effect is to create small magnetic field loops of roughly the same size as the disc thickness and alternating polarity centered along the mid-plane of the disc. The generalized curl of this vector potential

then gives the appropriate magnetic field components. The strength of the magnetic field is scaled to match the chosen $\beta$. In order to keep the magnetic field divergence free during the evolution, we use the staggered, constrained transport scheme described in Fragile et al. (2012).

In thin-disc simulations, a big challenge is to have enough resolution to capture the MRI. To help with this, we adopt a cylindrical, Kerr-Schild coordinate system computed by transformation from the usual spherical Kerr-Schild coordinate system ($R=r\sin\theta$ and $z=r\cos\theta$). To further improve the resolution near the disc mid-plane, we space the $n_{z}=160$ zones using a logarithmic coordinate

where $L_{z}=20H$ is the total box height and the sign of the expression comes from the sign of $z$. We also employ a logarithmic grid in the radial direction with $n_{R}=192$ zones spaced as

where $r_{\mathrm{BH}}=2\,r_{g}$ is the black hole radius. The radial domain runs over $4r_{g}\leq R\leq 40r_{g}$. A linear spacing is used along the azimuthal direction, $\phi$, with $n_{\phi}=32$ (low resolution) or $n_{\phi}=64$ (high resolution). To reduce the computational expense we only simulate the $0\leq\phi\leq 0.5\,\pi$ wedge in the azimuthal direction.

We use outflow boundary conditions in the inner radial and top and bottom vertical boundaries, which means all the gas variables except the normal component of the flow velocity are copied from the last active zones to the ghost zones. If the normal velocity component points outward, then it, too, is copied to the ghost zone. Otherwise, it is adjusted such that the normal component of the velocity is zero at the boundary face. The outer radial boundary uses a constant boundary condition, which means that all variables will retain their initial values in the outer radial ghost zones. The azimuthal boundary conditions are periodic.

We solve the coupled radiative magnetohydrodynamics equations using the $\textrm{M}_{1}$ closure scheme to handle both the optically thick and thin limits, as described in Fragile, Olejar & Anninos (2014). Along with the radiation stress-energy tensor, defined in equation (7), we need the MHD stress-energy tensor

where $\rho$ is rest mass density, $\varepsilon$ is specific internal energy, $P_{\mathrm{gas}}$ is gas pressure, defined using the ideal gas equation of state, $P_{\mathrm{gas}}=(\gamma-1)\rho\varepsilon$, with $\gamma=5/3$, and $b^{\alpha}$ is the contravariant magnetic 4-field, measured by an observer co-moving with fluid.

We aim to solve the following set of conservation equations for the mass

fluid stress-energy

and radiation stress-energy

together with the induction equation for the magnetic fields. In equations (17) and (18), $G_{\alpha}$ is the radiation 4-force density, which couples the fluid and radiation fields (Mihalas & Mihalas, 1984). This term includes normal scattering, absorption, and emission, as well as thermal Comptonization of the radiation. We do not include the relativistic corrections to thermal Comptonization, as we make the simplifying assumption that Compton scattering is symmetric in the fluid frame (hence, there is no associated momentum exchange). The form of $G^{\alpha}$ is

where

$\kappa_{\mathrm{s}}=0.34\,\mathrm{cm^{2}\,g^{-1}}$ is the opacity due to electron scattering, $\kappa^{\mathrm{a}}_{\mathrm{R}}=1.6\times 10^{21}T^{-7/2}\rho\,\mathrm{cm^{2}\,g^{-1}}$ is the Rosseland mean of absorption, $\kappa_{\mathrm{J}}^{\mathrm{a}}$ is the J-mean of absorption, and $\kappa_{\mathrm{P}}^{\mathrm{a}}$ is the Planck mean of absorption. We assume a nearly Planck spectrum, such that $\kappa_{\mathrm{J}}^{\mathrm{a}}=\kappa_{\mathrm{P}}^{\mathrm{a}}=6.4\times 10^{22}T^{-7/2}\rho\,\mathrm{cm^{2}\,g^{-1}}$. Here, $T_{\mathrm{rad}}$ and $T_{\mathrm{gas}}$ are the radiation and gas temperatures, respectively. Note that we do not solve independently for the temperature of the electrons, but simply assume it is equal to the temperature of the plasma. This should be sufficient in strongly coupled systems like the ones we simulate here, although in lower luminosity systems this does not hold (Ressler et al., 2015).

The form of the conservation equations that we actually solve can be written as

and

where $D=W\rho$ is the generalized fluid density, $W=\sqrt{-g}u^{t}$ is the generalized boost, $V^{i}=u^{i}/u^{t}$ is the fluid transport velocity, $\mathcal{E}=-\sqrt{-g}T^{t}_{t}$ is the total energy density, $\mathcal{S}=\sqrt{-g}T^{t}_{j}$ is the covariant momentum density, $\mathcal{R}=\sqrt{-g}R^{t}_{t}$ and $\mathcal{R}_{j}=\sqrt{-g}R^{t}_{j}$ are the conserved radiation energy density and momentum, respectively, and $\mathcal{B}^{j}=\sqrt{-g}B^{j}$ is the boosted magnetic field three-vector. The magnetic field, $B^{i}={}^{*}F^{\alpha i}$, is related to the co-moving field by

where ${}^{*}F^{\alpha\beta}$ is the dual of the Faraday tensor. These equations are solved using the explicit-implicit scheme described in Fragile, Olejar & Anninos (2014).

Each simulation is post-processed in order to extract the thermodynamic and geometric properties of the disc. We mainly use density-weighted shell- and time-averaged quantities, as well as spacetime diagrams to present our results. A general expression for the density-weighted shell-average of a quantity is given by the following expression:

where $dA_{R}$ is an area element normal to the radial direction. A time average of this quantity is computed as

For the height of the disc we use a density-squared weighting:

as it agrees better with our target height than other possible expressions, where $z_{0}=0$ represents the disc mid-plane.

We also track the photosphere of each disc, which we define as the $\tau=1$ surface. This is obtained by integrating the quantity $\kappa\rho$ from the lowest value of the $z$ coordinate on the grid, $z_{\mathrm{min}}$, to the height where $\tau=1$ and similarly from the highest value, $z_{\mathrm{max}}$, to the height where $\tau=1$:

where $\kappa=\kappa_{\mathrm{s}}+\kappa^{\mathrm{a}}_{\mathrm{R}}$. We emphasize here that the photosphere is always well within our simulation domain.

The local heating rate per unit surface area owing to turbulence caused by the MRI is computed within the volume enclosed by the photosphere as

where the integration is carried out within the photosphere and the integrand is azimuthally averaged, with $V^{\phi}=\Omega$ the azimuthal component of the three velocity and $W_{\hat{r}\hat{\phi}}$ the contravariant $r$-$\phi$ component of the MHD stress tensor in the co-moving frame. The radiative cooling is computed by tracking the radiative flux through the photosphere for each cylindrical shell:

where $F^{z}_{\mathrm{photo}}(R)=-4/3E_{R}u_{R}^{z}(u_{R})_{t}$ is the flux escaping through the top or bottom photosphere. As advective cooling is not important in these simulations, we ignore its contribution to $Q^{-}$. Also neglected is the contribution of the radial component of the radiative flux to cooling, which is appreciable only close to the ISCO.

The goal of this setup (RADPHR/RADPLR) is to study the thermal stability of a radiation-pressure-dominated disc. Before addressing the thermal stability, though, we present a general overview of this simulation.

Fig. 2 shows $R$-$z$ slices of the initial (left) vs. final (right) distribution of mass density (upper panel), temperature (middle panel), and the ratio of radiation pressure to gas pressure (lowest panel) for the RADPHR simulation. We clearly see in the upper-right panel that the disc collapses to a very thin structure. The middle-right panel shows that the disc cools down by about half an order of magnitude with respect to its initial gas temperature. This panel, especially, shows small, wave-like structures outside the main body of the disc. Those features are symptoms of primitive solver failures in the low-density background gas. In these cases, the internal energy that is used to calculate the temperature and gas pressure can introduce such artificial features. These errors do not appear to strongly influence the evolution of the disc itself. The lower-right panel shows that during the entire simulation the radiation pressure remains dominant within the disc.

In Fig. 3 we show additional $R$-$z$ slices of this simulation, now time averaged over the last 10 ISCO orbital periods ($t_{\mathrm{ISCO}}=92.3\,GM/c^{3}$). The top panel shows mass density again with streamlines of the azimuthally averaged fluid velocity. The photosphere is also shown with green solid curves. Again we see that the disc has collapsed to a very thin height. In this case, we terminated the simulation where we did, because we can no longer resolve the disc. The middle panel shows the radiation energy with streamlines of the azimuthally averaged radiation velocity (corresponding to the direction of the radiative flux vectors). We see in this panel that the radiation escapes mostly vertically from the disc, which is expected for geometrically thin discs. The lowest panel shows the ratio of the magnetic pressure to the sum of the radiation and gas pressures, $P_{\mathrm{mag}}/(P_{\mathrm{rad}}+P_{\mathrm{gas}})$, with lines of the azimuthally averaged magnetic field. There is a definite radial structure to the magnetic field lines in the background gas, while the field is more turbulent near the disc mid-plane. Magnetic pressure never dominates in the body of the disc during the entire simulation; thus, we do not expect the disc to be stabilized by the magnetic fields.

Fig. 4 shows the vertical acceleration components

in the coordinate frame. We can see that at the beginning of the simulation, as well as during the subsequent collapse, the disc is approximately in hydrostatic equilibrium. This demonstrates that the collapse is not caused by a loss of hydrostatic equilibrium. It is only after the collapse that hydrostatic equilibrium is violated (for $|z|<1$). We also note in the bottom panel that magnetic pressure dominates and even exceeds gravity above and below the disc. As we will see, this additional pressure support allows the photosphere to be outside the main body of the disc.

Fig. 5 shows the viscosity parameter, defined here as the density-weighted, height-averaged ratio of the (covariant) $\hat{r}$-$\hat{\phi}$ component of the stress tensor to the total pressure $\alpha\equiv\langle W_{\hat{r}\hat{\phi}}/P_{\mathrm{tot}}\rangle_{\rho}$. We find a nearly constant and uniform value of $\alpha=0.02.$

In Fig. 6 we show a spacetime plot of the azimuthally averaged radial profile of disc height, computed using equation (31). The figure shows that the height initially increases after one local MRI growth time. Subsequently, the disc collapses on a timescale comparable to the local equilibrium cooling time (solid white curve), which is estimated from thin disc theory to be $t_{\mathrm{cool}}(R)=2\pi/(\alpha\Omega)$. The disc collapse negates the initial expansion such that the final height is at least a factor of five smaller than the initial one. We show a similar spacetime plot for the height of our low-resolution, RADPLR case in Fig. 7. Comparing Figs. 6 and 7, we see that it takes somewhat longer for the disc to collapse in our high-resolution simulation. It could be that at even higher resolution we might find that the disc will take even longer to collapse. This point will have to await future simulations.

The vertical collapse of the disc can also be seen in Fig. 8, which shows a spacetime plot of the radially and azimuthally averaged density as a function of height, $z$. The narrowing of the high density region of the plot with time reflects the collapse of the disc. As the disc is collapsing, the photosphere (black curve) remains at a relatively constant height. This is because there is a magnetically supported atmosphere at $|z|\geq 1GM/c^{2}$ (see Fig. 4). Recall, too, that the photosphere is calculated by integrating the opacity from the domain boundaries toward the mid-plane.

Another clue to the thermal state of the disc comes from looking at the time evolution of the heating and cooling rates. In Fig. 9, we show a spacetime plot of the ratio of the heating rate to the cooling rate, computed using equations (33) and (34), respectively. The important takeaway is that cooling dominates over heating everywhere in the disc until after the disc has collapsed. This is even true at the outset of the simulation. As we noted, this simulation does not begin in thermal equilibrium. Heating does appear to balance cooling inside the ISCO, but this has little impact on our results.

As mentioned at the beginning of this section, the RADPHR and RADLR simulations are initialized with the disc being radiation pressure dominated. In Fig. 10 we show that it remains so throughout the simulation. Obviously the region inside the ISCO is strongly radiation pressure dominated. This is not surprising as this region is largely devoid of gas, but filled with radiation (compare the top and middle panels of Fig. 3).

To demonstrate a thermal instability, we need to show that the heating and cooling rates responding differently to changes in the total mid-plane pressure, which we do in Figs. 11 and 12. In particular, Fig. 11, shows that the cooling rate, $Q^{-}$, is proportional to the mid-plane pressure, as expected from the radiative diffusion equation. Note that $Q^{-}$ has been scaled in the figure with one power of total pressure. At early times ($t<2000\,GM/c^{3}$), and especially at small radii ($R<12\,GM/c^{2}$), we see in Fig. 12 that the heating rate, $Q^{+}$, scales as the square of the mid-plane pressure. Note that in this figure, $Q^{+}$ has been scaled by the square of total pressure. Thus, the dependence on pressure is steeper in the case of heating ($Q^{+}\propto P^{2}_{z_{0}}$) than in the case of cooling ($Q^{-}\propto P_{z_{0}}$), which is the hallmark of the thermal instability. So, even if this simulation had started in thermal equilibrium, it would have still been prone to a thermal runaway.

At later times the discussion of heating and cooling, and relating the numerical results to classic theory of accretion discs, is complicated by the presence of two different scale-heights in the disc. The density scale height $H$ illustrated in Figs. 6 and 7 is very different from the height of the photosphere in Fig. 8. We remind the reader that the cooling rate $Q^{-}(R)$ is the flux leaving through the photosphere at radius $R$, and the heating rate $Q^{+}(R)$ is computed in a cylindrical shell of the same radius from the lower to the upper photosphere.

In Fig. 13, we show a spacetime plot of the surface density, $\Sigma$, of the disc. First, we note that $\Sigma$ remains nearly constant throughout the simulation. This is important as it demonstrates that the disc collapse is not caused by matter being drained into the black hole. However, a closer inspection of Fig. 13 shows that the disc seems to collect into distinct rings after a time of $\sim 1000\,GM/c^{3}$. For instance, a very dense ring forms at $R\approx 7\,GM/c^{2}$. As seen in Fig. 14, these really are rings and not spiral structures. This is very suggestive of the viscous instability (Lightman & Eardley, 1974). Additional rings are seen forming at $R\approx 11,13$, and $15\,GM/c^{2}$. Thus, the spacing between the rings in the inner disc is roughly the same as the height of the disc ($\Delta R\sim H$), as expected in the early stages of the viscous instability (Lightman & Eardley, 1974). The anticipated growth timescale for the instability is roughly the viscous timescale multiplied by the relative size squared of the structures being formed, i.e. $t_{\mathrm{LE}}\sim(R/H)^{2}/(\alpha\Omega)\times(\Delta R/R)^{2}$. For the innermost ring, this is approximately $900\,GM/c^{3}$, roughly consistent with the first appearance of the ring in Fig. 13. If this collection of mass into rings is indeed owing to the viscous instability, it should also be accompanied by a rise in the viscous stress in the regions between the rings since $W_{\hat{r}\hat{\phi}}\propto\Sigma^{-1}$. Comparing Figs. 13 and 15, we do see that the prominent low-$\Sigma$ gap between $R=8$ and $10\,GM/c^{2}$ in Fig. 13 corresponds to the highest stress region in Fig. 15, providing additional evidence that we may be seeing the viscous instability in a simulation for the first time.

In the previous section we described a radiation-pressure-dominated thin disc which is unstable to collapse. In this section, we contrast that result with a baseline, gas-pressure-dominated simulation. The primary goal in performing this simulation is to ensure that the result of the previous section is not a numerical artifact and that, indeed, our numerical setup is able to evolve stable thin disc configurations for sufficiently long duration. Here we focus on the low-resolution simulation, GASPLR, as it ran for significantly longer than its high-resolution counterpart. Even so, it did not run long enough to achieve inflow equilibrium through much of the radial domain. If the goal was really to understand the behavior of this disc, then this simulation should be run even longer. However, we were able to run it for many thermal times, long enough to provide a convincing counter-example to our collapsing disc.

We start with a spacetime plot of the disc height (Fig. 16), comparable to Figs. 6 and 7. In it we see that the region inside of $R\approx 12\,GM/c^{2}$ maintains a nearly constant height following the onset of the MRI (dashed line). However, the GASPLR simulation does exhibit some thinning in its outer regions, an effect that was also seen in the radiation-pressure-dominated simulations. We suspect this is owing to the fact that the MRI is not sufficiently well resolved in those regions (see Sec. 3.4), which allows cooling to dominate.

In Fig. 17 we show a space-time plot of the ratio of heating to cooling in the GASPLR case. We see that heating and cooling are nearly in balance over time over most of the disc. However, in the outer regions and at late times, there are parts of the disc where cooling dominates. As with the collapsing scale height in these regions, we attribute this behavior to the under resolved MRI.

In Fig. 18 we show the vertical profile of the radially and azimuthally averaged mass density for the GASPLR simulation. As with Fig. 16, we find that the gas-pressure-dominated disc maintains its height much better than the radiation-pressure-dominated one. The growth in the mid-plane density at late times is dominated by the thinning of the outer regions ($R>16\,GM/c^{2}$) of the disc (also seen in Fig. 16).

One major concern in performing numerical simulations of very thin accretion discs is resolving the MRI. Multiple resolution studies have shown that nearly all global simulations done to date have been under-resolved in this regard (e.g. Hawley, Guan & Krolik, 2011; Hawley et al., 2013). A particular worry in our case is that under-resolving the MRI might cause the disc to collapse due to insufficient heating. Admittedly, as we show, some of our simulations fall below the ideal resolution. This is an unfortunate limitation of trying to perform global simulations, and doing so for thin discs only exacerbates the problem. However, an attempt to perform such thin disc simulation is necessary to make progress in understanding the physics of thin discs.

The standard measure of MRI resolution is the so-called $Q$ parameter, defined as

where $\lambda_{\mathrm{MRI},i}=2\pi v_{A,i}/|V^{\phi}|$ is the wavelength of fastest growing MRI mode, $\Delta x_{i}$ is a typical zone length, and $v_{A,i}=\sqrt{b^{i}b_{i}/\rho}$ is Alfven speed, all in a given direction, $i$. We checked both the vertical (Fig. 19) and azimuthal (Fig. 20) MRI $Q$ parameters for most of the models in Table. 1. For our radiation-pressure-dominated simulations, we actually capture the vertical MRI well, with $Q_{2}\gtrsim 10$ throughout most of the disc, especially in simulation RADPHR. The gas-pressure-dominated simulation is not nearly as well resolved, with $Q_{2}\lesssim 3$ at nearly all radii. The azimuthal direction is more problematic, with $Q_{3}\lesssim 5$, even in the high-resolution RADPHR case. One way to possibly improve this in future simulations without requiring even more computational resources would be to reduce the azimuthal extent of the domain, while keeping the number of zones fixed. For now, we are left to point to the similarity of our results in both the low- and high-resolution simulations as evidence that the poor azimuthal MRI resolution does not negate our results of a thermal collapse in our radiation-pressure-dominated disc.