TEMPERATURE FLUCTUATIONS DRIVEN BY MAGNETOROTATIONAL INSTABILITY IN PROTOPLANETARY DISKS

We performed simulations of unstratified, magnetized, shearing boxes with net vertical field including Ohmic resistivity and radiative thermal diffusion using the Pencil Code,⁴ a sixth-order, central difference in space method, (Brandenburg & Dobler 2002), as well as confirming our major results with Athena⁵ (Stone et al. 2008; Stone & Gardiner 2010), a constrained-transport, Godunov method. Here we use Athena configured with orbital advection, linear reconstruction, and the HLLD Riemann solver.

The field variables evolved in the Pencil Code were density, velocity, magnetic vector potential, and thermal energy density. The Pencil Code requires diffusion operators for stability, so we employed a grid-scaled hyperdiffusion in velocity and density (which have no physically resolved diffusion effect in the model), along with shock diffusion on all fields so that all shocks are captured at the grid scale. Because the Pencil Code is sixth order accurate in space and third order accurate in time, it resolves small scale subsonic flow structures almost as well as a spectral code. However, it does not conserve energy exactly. To confirm our results, we replicated the configuration in Athena, a locally energy-conservative, finite volume code. In such a method, no explicit stabilizing diffusion on the momentum and density fields is required. Furthermore, the conversion of energy from kinetic and magnetic energy in shocks and current sheets to thermal energy is fully conservative.

We included two explicit physical diffusive effects in our model, resistivity η, and radiative thermal diffusion. We allowed neither η nor the opacity κ to vary with temperature, unlike in an actual disk. The form of our Ohmic resistivity is

$$\begin{eqnarray} &&\eta \equiv \frac{c^2}{4 \pi \sigma } \end{eqnarray} \tag{ 2 }$$

in Gaussian cgs units, where σ is the conductivity. The evolution of the magnetic field ${\boldsymbol {B}} {}$ is then described by the induction equation

$$\begin{equation} \frac{\partial {\boldsymbol {B}} {}}{\partial t} = \nabla \times \left(\left[{\boldsymbol {v}} {}- \frac{3}{2} \Omega _0 x \hat{\boldsymbol {y}} {}\right] \times {\boldsymbol {B}} {}\right) - \frac{4\pi }{c} \nabla \times [\eta {\boldsymbol {J}} {}], \end{equation} \tag{ 3 }$$

where the current

$$\begin{equation} {\boldsymbol {J}} {}\equiv \frac{c}{4 \pi } \nabla \times {\boldsymbol {B}} {}, \end{equation} \tag{ 4 }$$

and the shearing box fluctuation velocity is ${\boldsymbol {v}} {}$. The $(3/2) \Omega _0 x \hat{\boldsymbol {y}} {}$ term is the orbital shear.

So long as the resistivity remains constant in space,

$$\begin{equation} \frac{\partial {\boldsymbol {B}} {}}{\partial t} = \nabla \times \left(\left[{\boldsymbol {v}} {}- \frac{3}{2} \Omega _0 x \hat{\boldsymbol {y}} {}\right] \times {\boldsymbol {B}} {}\right) -\eta \nabla \times [\nabla \times {\boldsymbol {B}} {}]. \end{equation} \tag{ 5 }$$

The resistive term was integrated with the same scheme as the other operators. The Pencil Code evolves an analogous equation for the magnetic vector potential ${\boldsymbol B} = \nabla \times {\boldsymbol A}$.

The thermal energy evolution includes a diffusive term based on radiative transfer in the optically thick limit with Rosseland mean opacities, as well as a thermal relaxation term to model the large scale energy transfer away from the midplane of the disk. The full thermal energy equation is thus

$$\begin{eqnarray} \frac{\partial e}{\partial t} =& - \nabla \left(\left[{\boldsymbol {v}} {}- \frac{3}{2} \Omega _0 x \hat{\boldsymbol {y}} {}\right] e\right)- p_\mathrm{th}\nabla \cdot {\boldsymbol v} \nonumber \\ &+ \frac{4 \pi \eta }{c^2} J^2 -\nabla \cdot {\boldsymbol F} - C_v \frac{T-T_0}{\tau _0}, \end{eqnarray} \tag{ 6 }$$

where $p_\mathrm{th}=(\gamma -1)e=\rho k_B T/\bar{m}$ is the thermal pressure, $C_v=(\gamma -1)^{-1} k_B/\bar{m}$ is the heat capacity at constant volume, T is the gas temperature, T₀ the reference temperature, τ₀ the thermal relaxation time, k_B the Boltzmann constant, $\bar{m}$ the mean mass per particle, and with the Rosseland mean opacity thermal diffusion flux ${\boldsymbol F}$ given by

$$\begin{equation} {\boldsymbol F} = -\frac{16 \sigma _B T^3}{3 \rho \kappa } {\boldsymbol {\nabla }} {}T, \end{equation} \tag{ 7 }$$

where σ_B is the Stefan–Boltzmann constant, ρ the density and κ the Rosseland mean opacity.

The thermal diffusion treatment is motivated by the high optical depth expected over length scales of the MRI unstable wavelength expected for the outer edge of the MRI active region in the dusty inner disk. In this region λ_MRI ≃ H, so the optical depth of an MRI wavelength is comparable to that of the disk itself. For our parameters, λ_MRI has an optical depth of τ > 15,000. We choose the thermal relaxation timescale to be an orbital period τ₀ = 2π/Ω₀ (Table 1). We make this choice to probe the coolest reasonable disk state, which corresponds to one that can cool so quickly that it approaches marginal gravitational stability. The thermal relaxation rate is chosen to be fast, so that it strongly limits the heating resulting from the energy input by the background shear. Additionally, we can place an upper limit on the physically plausible cooling rate by recognizing that the model is aimed to produce the environment at the inner edge of a dead zone. If the cooling rate of the neighboring dead zone was much faster than $\tau _0=\Omega _0^{-1}$, the region would become gravitationally unstable (Rice et al. 2014), and would thus be prone to generating strong spiral waves, strongly altering the global disk structure. As such, we can argue that this region of the disk must, to be consistent with an MRI-unstable quasi-steady state at the inner edge of the dead zone, have a cooling timescale longer than $\tau _0=\Omega _0^{-1}$. Therefore, our choice of $\tau _0=\Omega _0^{-1}$ should be a conservative one, and should result in a model that underestimates the temperature fluctuations.

We chose the parameters listed in Table 1 to approximate a disk that is marginally MRI active at 1 AU. Conceptually, these parameters can be thought of as corresponding to a region at the inner edge of a dead zone. Unusually for astrophysical magnetohydrodynamics, effects of the microphysical resistivity are resolvable at the inner edge of the dead zone (unlike the unresolvable microphysical viscosity) because there the MRI is marginally super-critical and numerical simulations can clearly show the growth of magnetic fields in an MRI dynamo.

The marginal instability criterion for the MRI is that the most unstable wavelength of the MRI be comparable to the disk scale height. In our case, the MRI wavelength λ_MRI and the effective scale height of the disk $H \equiv c_s/\sqrt{\gamma } \Omega$ are within 10% of each other. This means that even though we have performed unstratified numerical simulations, we can link our length scales to those of stratified simulations near MRI criticality. With these parameters, the initial net vertical magnetic field corresponds to an initial plasma β₀ = 750. This is low compared to many existing models with zero resistivity. However, in this very resistive regime the wavelength of the fastest growing MRI mode increases with β₀. Accordingly, a relatively strong initial field is needed to make the most unstable length scale match the nominal disk scale height. The wavenumbers corresponding to the vertical heights of the shearing box domains used in this study are noted in Figure 1 along with the dispersion relation curve for the MRI.

While temperature is not strictly diffused, for the parameters in Table 1, extrapolating from Equation (7) and using $T=(\gamma -1)\bar{m} e/\rho k_B$, an effective thermal temperature diffusion coefficient can be estimated as

$$\begin{equation} \mu \equiv \frac{16 (\gamma -1) \bar{m} \sigma _B T_0^3}{3 \rho _0^2 \kappa k_B} \simeq 1.8 \times 10^{14} \ \frac{{\rm cm}^2}{{\rm s}}, \end{equation} \tag{ 8 }$$

where σ_B is the Stefan–Boltzmann constant, and the other parameters are as defined in Table 1. When measured with respect to the MRI wavelength and orbital timescale of our setup, we have

$$\begin{equation} \mu \sim \frac{\lambda _{{\rm MRI}}^2 \Omega _0}{1000}. \end{equation} \tag{ 9 }$$

Accordingly, the radiative cooling time for a structure with the size of the most unstable MRI wavelength is approximately 130 yr, far longer than the imposed cooling time of 1 yr. This means that large scale thermal structures are dissipated by the thermal relaxation cooling term, but small scale structures (less than 0.1λ_MRI) cool radiatively.

We used the Pencil Code to perform two sequences of runs. The first sequence uses a cubic volume with resolutions from 64³ to 512³, with 5λ_MRI per box length. Our highest two resolutions are equivalent to 50 and 100 zones per scale height, respectively. The higher resolution runs were initialized from the next highest resolution run. The 64³ simulation was run for t = 105 orbits. At t = 45 orbits, it was remeshed to a resolution of 128³ and this was also run to t = 105 orbits. At t = 60 orbits, the 128³ run was remeshed to 256³ and this was run to t = 105 orbits. Finally, at t = 75 orbits, the 256³ run was remeshed to 512³ and this was run to t = 105 orbits. The second sequence of runs repeats the same parameters with a slab volume one quarter the height, that is with aspect ratio L_x:L_y: L_z of 1: 1: (1/4). Thus, these runs use a shearing box with a height of approximately one scale height. These runs preserved the same cell-size as the first sequence of runs, so the highest resolution was 512² × 128. Otherwise the procedures and analysis of this sequence of runs follows the cubical runs.

The comparison runs done with Athena used the cubical grid and did not use remeshing. Both started from t = 0 and ran for 50 orbits. Resolutions of 128³ and 256³ are reported here. The setup is otherwise identical to that used with the Pencil Code.

We first consider the behavior of the global average quantities in our models. In the top panels of Figures 2 and 3, we show time series of the volume average plasma beta

$$\begin{equation} \overline{\beta }=\frac{8 \pi \langle p_\mathrm{th} \rangle _V}{\langle B^2 \rangle _V}, \end{equation} \tag{ 10 }$$

where 〈⋅⋅⋅〉_V denotes volume averaging, for our different resolution runs. The volume average of the magnetic field was taken separately to avoid having this average skewed by very high values of β within demagnetized reconnecting regions. Our cubical models settle at values of $\overline{\beta } \sim 10$, while the slab models only reach $\overline{\beta } \sim 50$. Similarly, the bottom panel of both figures shows the volume averaged Shakura & Sunyaev (1973) α turbulent viscosity parameter:

$$\begin{equation} \alpha =\frac{1}{p_\mathrm{th}} \left(\rho v_x v_y - \frac{B_x B_y}{4\pi } \right). \end{equation} \tag{ 11 }$$

The volume averaged value of this parameter settles at a value of about 〈α〉_V = 4 × 10⁻² in the cubical case and 〈α〉_V = 3 × 10⁻² in the slab case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 4 and 5, we show temperature and heating data from the Pencil Code runs. Two facts jump out. First, there is an order unity difference between the tenth and ninetieth percentile temperatures at all times. Second, the dominant energy dissipation mechanism is resistive rather than compressive or shock dissipation, because when averaged over the volume the adiabatic heating and cooling approximately cancel. This occurs even though the compressive heating has much stronger peak values than the resistive heating, as we show in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To more closely investigate the source of the temperature variations, we turn to a study of the morphology of current sheets. We find that the heating is usually dominated by one or two major sheets at any given time. In Figure 6, we show a set of slices at constant azimuth through our highest resolution cubical Pencil Code run. In the slice chosen, there is a current sheet lying in the radial–azimuthal plane, perpendicular to the slice plane, and visible in all variables except the total pressure, although only barely visible in the adiabatic heating and cooling panel. This figure makes it clear that, even though the volume average $\overline{\beta } \simeq 8$ at the orbit in question, the local minimum ratio associated with the peak magnetic field is β_p ≃ 1.

Zoom In Zoom Out Reset image size

Furthermore, the largest temperature variation generally traces the highest resistive dissipation in the current sheet structure, although there is enough difference between the two to make clear that it is time-averaged, rather than instantaneous, heating that determines the temperature perturbation. Although the compressive heating is high along the weak shocks filling the domain, these are accompanied by large regions of expansion cooling. Hence, as was seen in the volume average shown in the upper panel of Figure 4, the heating is dominated by the resistive dissipation.

Similarly, Figure 7 shows a set of slices at constant azimuth through our highest resolution slab run. The qualitative pattern of dominant azimuthal field bundles is similar, and a strong, hot current sheet structure can be seen. Consistent with the results for the globally averaged quantities, the values are less extreme in this vertically restricted domain than in the cubic one: the volume averaged $\overline{\beta } \simeq 20$ and the local minimum β associated with the peak magnetic field is β_p ≃ 1.

Zoom In Zoom Out Reset image size

This difference between the local β_p and volume averaged $\overline{\beta }$ is straightforward. The strongest magnetic field structures of the MRI are generated by orbital shear, which stretches radial magnetic field lines into nearly azimuthally constant bundles of strong, azimuthally directed magnetic flux. If the magnetic field is dominantly contained in azimuthally constant structures varying along a single large scale wave-vector $\boldsymbol {k}$, we would expect $\beta _p \simeq \overline{\beta }/4$ because we need to average over sin (k_xx)²sin (k_zz)². If β_p is of order unity, its value is expected to further drop because the magnetic field bundles will have low thermal pressure compared to the volume averaged thermal pressure (because the magnetic pressure pushes fluid out of the high magnetic field regions). In Figure 8, we show a vertical cut through the current sheet identified in Figure 6 that demonstrates that this inverse magnetic and thermal pressure correlation occurs in large current sheets found in the highest resolution aspect ratio cubical model.

Zoom In Zoom Out Reset image size

These figures also show that even though the simulation has β_p ∼ 1 and hence drives mildly supersonic flows, the large temperature variations are due to resistive heating. Regions heated by the hydrodynamical shocks quickly reexpand and adiabatically cool, leaving little imprint on the temperature structure of the gas.

Figures 2 and 3 compare the time variation of $\overline{\beta }$, 〈α〉_V, in runs of varying resolution for the two aspect ratios. The average value for the highest resolution runs are shown by the dashed lines to guide the eye. Similarly Figures 4 and 5 compare the temperature and heating rates. Episodes of very low field strength are common for the lower resolution runs, which are less well-resolved, and hence subject to stronger numerical resistivity. Our simulations appear to have converged to at least the intrinsic time variation in these averaged quantities for resolutions of 128³ or larger. However, the temperature and the resistive heating rate converge only for the two highest resolutions (256³ and 512³, or 50 and 100 zones per λ_MRI or equivalently per scale height).

The motivation for the thinner domain is to avoid the cycles of strong channel flow and intermittency that have been observed to be exacerbated by the use of a cubic domain (Bodo et al. 2008). Indeed, the variation on an ∼10 orbit timescale appears less in these runs than in the cubical ones. However, channel mode-like pairs of azimuthal field bundles still dominate the magnetic field, as can be seen in the top left panel of Figure 7. The heating is less extreme for the thinner box, but the largest MRI active scales are also truncated by the vertical height (see Figure 1). Accordingly, the volume average $\bar{\beta }$ is on the order of five times higher in these runs.

We explain the increase in volume averaged heating rate with resolution by reference to the theory described in Section 5.2. Magnetic fields are current sources: ${\boldsymbol {J}} {}\propto {\boldsymbol {\nabla }} {}\times {\boldsymbol {B}} {}$. The MRI often produces oppositely directed azimuthal magnetic field bundles. Given two such bundles with strength B₀/2 separated by a distance L, the current flowing between them scales as J ∼ B₀/L. However, resistive dissipation scales with the square of the current, $\eta |{\boldsymbol {J}} {}|^2$, so the net resistive heating scales as $\eta B_0^2/L^2$. If the current sheet is not resolved, the heating rate will be understated because the actual value of L, set by the resolution, is larger than the physical one set by η.

The high resolution required to resolve the current sheets, even at the high resistivity edge of the dead zone, can be understood from dimensional considerations. Numerical resistivity η_n∝u × Δx, where u is the turbulent velocity. If the imposed resistivity

$$\begin{equation} \eta \gg \eta _n, \end{equation} \tag{ 12 }$$

then the physical resistivity η dominates over the numerical resistivity η_n in Equation (3) and resistive effects are resolved.

Assuming equipartition between the turbulent kinetic and turbulent magnetic energies, $u^2 \sim v_A^2$, so η_n ∼ v_A × Δx, where v_A is the Alfvén speed of the saturated magnetic field. The ratio of η_n for the saturated state to η_n0 of the initial state is η_n/η_n0 ≃ v_A/v_A0 = (β_i/β_s)^1/2, where β_i is the plasma beta of the initial field, and β_s is that of the saturated field. This ratio is often an order of magnitude or larger.

This can cause issues even near the relatively high resistivity edge of the dead zone where Equation (12) is best satisfied. The MRI criticality condition is that the initial Elsasser number $\Lambda _0 =v_{A0}^2/\eta \Omega$ be of order unity, where v_A0 is the Alfven speed of the seed field. The magnetic Reynolds number of the saturated state is

$$\begin{equation} {\rm Re}_M\equiv u^2 t_t/\eta, \end{equation} \tag{ 13 }$$

where t_t is the turbulent turnover time, t_t ≃ Ω⁻¹, and taking u to be the maximum shearing-box fluctuation velocity. In our case, the initial Elsasser number is Λ₀ = 0.5 and at late times the turbulent magnetic Reynolds number is  Re_M ∼ 100, which is a non-trivial value of  Re_M to resolve numerically. Note that  Re_M/Λ₀ = (η_n/η_n0)², so the grid resolution needed to resolve the turbulent flow in the saturated state is a factor of over 14 times more stringent than at the initial time. All our runs, including the lowest resolution, 64³ ones, resolved the initial state according to Equation (12). During the highest energy episodes, however, even the highest resolution 512³ simulations did not quite satisfy Equation (12), though only by a factor of about two. To understand the actual convergence behavior, we need to move beyond that criterion's simple dimensional analysis.

To better probe the resolution required to resolve resistive dissipation, we examine the structure of the current density. Our current sheets show multiple lengthscales: the full box length in azimuth, λ_MRI in width and a far shorter length perpendicular to the sheet. This rules out a simple Fourier analysis of the current sheet structures. Instead, to examine the convergence of the complex current sheet structures, we decomposed $|{\boldsymbol {J}} {}|^2$ with a multiresolution decomposition (Harten 1994).

We successively coarse-grained the data, and subtracted the coarse-grained data from the original. We call the results of the subtraction residuals, although they are described as error coefficients in other contexts. They measure how much structure there is on a given scale and hence at what scale the dominant features of $|{\boldsymbol {J}} {}|^2$ exist. The decomposition is conservative, in that the sum of the residuals on every level over the entire volume is zero. To detect how much structure in the field resides on each level, we examine the average absolute value of the residuals on each grid level.

The transform used to calculate the residuals is formally

$$\begin{equation} v^L \longleftrightarrow \left\lbrace e^L,\ldots,e^1, v^0 \right\rbrace. \end{equation} \tag{ 14 }$$

The function v is on the finest mesh v^L, the residuals on each level N are e^N, and v⁰ is the representation of the function v on the coarsest mesh: a single point. The transform uses a decimation operator to move the representation of the function v from mesh level N to a mesh with half as many points in each direction on level N − 1 by

$$\begin{equation} v^{N-1}_{i,j,k} = \frac{1}{8} \sum _{l,m,n=0}^1 v^N_{i+l,j+m,k+n} \end{equation} \tag{ 15 }$$

which is a simple average over a cube of eight neighboring points. In this way, the volume integral of v is conserved in the successively coarser representations. Thus the most coarse representation v⁰ is a scalar value that is equal to the volume average of v^L. The residuals on each level are defined by the difference between the representation of the function v on level N and level N − 1 as

$$\begin{eqnarray} e^N_{2i+l,2j+m,2k+n} &=&v^N_{2i+l,2j+m,2k+n} - v^{N-1}_{i,j,k} \nonumber \\ && l,n,m \in \lbrace 0,1\rbrace . \end{eqnarray} \tag{ 16 }$$

The mean absolute value of the residual e^N on level N, denoted as 〈|e^N|〉 is then a measure of the total volume-integrated changes between v^{N − 1} and v^N.

If we examine the set of values {〈|e¹|〉, ..., 〈|e^L|〉}, the location of the largest residual corresponds to the grid scale where the original function v^L changes most. The single |v⁰| value measures the volume averaged $|{\boldsymbol {J}} {}|^2$. To account for the temporal variations of the system, we compute this decomposition once per orbit, and report averages of the absolute values of the coefficients, normalized to the time average of |v⁰| from the highest resolution simulation.

The results of this analysis are shown in Figures 9 and 10. The time averages were performed over the interval t = 80 to t = 105. The |v⁰| values show that the volume integrated $|{\boldsymbol {J}} {}|^2$ only approaches convergence at resolution of 256³ to 512³. Both the cubical and slab simulations show the same convergence behavior. The averaged residuals e^N consistently peak at the scale corresponding to level 5; as resolution increases, the finest scales show decreasing e^N coefficients. The strength of the residuals converges well at the highest resolutions. Note that e⁵ corresponds to the structure of  J² on a mesh with resolution 32³, i.e., to a scale of approximately λ_MRI/6. Hence, even though $|{\boldsymbol {J}} {}|^2$ shows structure on all scales, our highest resolution simulations are well able to resolve the dominant scales of current sheet dissipation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify that the results from the Pencil Code resolution study are not dominated by any possible local, nonconservative, energy dissipation, we have performed a smaller set of runs with Athena, which conserves energy to numerical accuracy. These runs all start from t = 0 and are run to t = 50 orbits. As expected, they display a different convergence behavior and should be compared to the highest resolution Pencil Code result. The volume average $\overline{\beta }$ at 128³ and 256³ shown in Figure 11 gives reasonable agreement with the Pencil Code result, generally varying in the range $\overline{\beta } =$ 5–10. The 〈α〉_V viscosity also agrees reasonably, varying in the range 0.05–0.1, as shown in Figure 12. Finally, the variation in temperatures seen in Athena in Figure 13 also agrees reasonably, with the mean temperature typically ∼1200 K and the spread between the 10th and 90th percentile values being typically ∼200 K. As the Godunov method employed in Athena does not capture shocks through an explicit shock viscosity, it is difficult to decompose the heating rates as was done for the Pencil Code runs. Altogether, the Athena results are entirely consistent with the major result that current sheets can cause substantial local heating under the conditions we consider.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The comparatively short perpendicular length scale of these current sheets implies that as long as β ⩾ 1, structures that are smaller than a scale height H and survive for at least an orbital timescale Ω⁻¹ should be in near pressure equilibrium with the bounding magnetic field bundles. Accordingly, we can estimate the temperature variation by assuming that adiabatic compression maintains pressure equilibrium between the current sheet and the surrounding magnetic field bundles. Outside of the current sheet, the pressure is altered by both the growth of the magnetic field during the formation of the bundles and adiabatic expansion into the current sheet, while inside the current sheet the pressure is altered solely by adiabatic compression.

If we assume a sinusoidally varying magnetic field, and the maximum value of the magnetic pressure is small enough that we can expand linearly, then the adiabatic compression in the current sheets is equal to the adiabatic expansion outside the current sheets. Thus, if ρ₀, T₀ and p₀ are the density, temperature, and pressure in the disk before the generation of the field bundles, the pressure p₁ of the system after magnetic field growth and pressure equilibration occurs will be

$$\begin{equation} p_1=p_0+\delta p_{{\rm th}}=p_0-\delta p_{{\rm th}} + \delta p_m, \end{equation} \tag{ 17 }$$

where δp_m ≃ p₀/β_p is the maximum magnetic pressure increase, β_p is the minimum plasma β associated with the peak magnetic field, and δp_th is the amplitude of the thermal pressure variation. This implies that

$$\begin{equation} \delta p_{{\rm th}}= \frac{p_0}{2 \beta _p}. \end{equation} \tag{ 18 }$$

The temperature variation between the magnetized regions and the current sheets is

$$\begin{equation} \frac{\delta T}{T_0} = \left(\frac{p_0+\delta p_{{\rm th}}}{p_0-\delta p_{{\rm th}}}\right)^{(\gamma -1)/\gamma } -1 \simeq \frac{\gamma -1}{\gamma } \frac{1}{\beta _p}. \end{equation} \tag{ 19 }$$

For our value of γ = 1.5, this reduces to δT/T₀ ≃ 1/(3β_p).

Equation (19) implies that strong magnetic fields will always generate significant point-to-point temperature fluctuations through adiabatic compression alone, but this effect alone clearly is insufficient to explain the simulated behavior if $\beta _p \sim \overline{\beta } \sim 10$. However, there are two further effects we must consider. Thermal diffusion will reduce the temperature variation, but resistive heating in the current sheet will increase it.

In this section, we analyze the time-dependent, resistive heating of current sheets for magnetic fields whose peak magnetic field is weak, with a minimum value of the magnetic to thermal pressure ratio $\beta _p^{-1} \ll 1$. We expect that the azimuthal magnetic field bundles that bound the sheets form from the amplification of radial field variations by orbital shear. Such a mechanism should result in anticorrelation of radial and azimuthal magnetic fields. The top panel of Figure 8 demonstrates that this is indeed the case in our simulation.

We can therefore approximate the system as an initial radial field that varies only as a function of z and decays resistively. We do not include any variation in the x direction as the current sheets are thin, extended structures. As the azimuthal field bundles that appear in the simulations span the full box in azimuth, we also do not include any variation in the y direction. However, such y-direction variations could be an interesting study, as they would sharpen with the action of the background shear. That radial field is sheared by a flow of u_y = Sx where S = −3Ω₀/2 is the Keplerian shear. The shear then generates an azimuthal field that also decays resistively. Under these conditions, the induction Equation (5) becomes (Brandenburg & Subramanian 2005):

$$\begin{eqnarray} &\frac{\partial B_x}{\partial t} = \eta \frac{\partial ^2 B_x}{\partial z^2} \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} &\frac{\partial B_y}{\partial t} = -\frac{3\Omega _0}{2} B_x +\eta \frac{\partial ^2 B_y}{\partial z^2}. \end{eqnarray} \tag{ 21 }$$

Assuming an initial radial field of B_x = B₀sin (kz), and no initial azimuthal field, the time evolution of the magnetic field is

$$\begin{eqnarray} &B_x (t) = B_0 \exp ({-t/\tau }) \sin (kz) \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &B_y (t) = -\frac{3\Omega _0 t}{2} B_x = -B_0 \left(\frac{3\Omega _0 t}{2}\right) \exp ({-t/\tau }) \sin (kz), \nonumber\\ \end{eqnarray} \tag{ 23 }$$

where τ ≡ 1/(ηk²) is the resistive decay time associated with the wavenumber k.

The peak azimuthal field occurs at time t = τ. As long as the ratio of the resistive timescale τ to the shear timescale S⁻¹, namely 3Ω₀τ/2 ≫ 1, B_y dominates the magnetic field, with a peak strength of

$$\begin{equation} B_p = \frac{3\Omega _0 \tau }{2 \exp (1)} B_0. \end{equation} \tag{ 24 }$$

In what follows, we consider only the heating due to this dominant azimuthal field, neglecting B_x for all purposes other than feeding B_y.

We normalize the peak azimuthal field strength to the initial thermal pressure by defining the plasma-β value for this peak azimuthal field with respect to the initial thermal pressure:

$$\begin{equation} \beta _p \equiv \frac{8 \pi p_0}{B_p^2}. \end{equation} \tag{ 25 }$$

We assume that β_p ≫ 1, so that we can linearize in $\beta _p^{-1}$. In this limit, defining β_p with respect to p₀ rather than the thermal pressure at the peaks of the perturbed magnetic field, where sin (kz) = ±1 is equivalent to expanding first order in $\beta _p^{-1}$ because the thermal pressure perturbation does contribute at lowest order. We can use Equations (24) and (25) to rewrite the initial seed field strength in terms of β_p:

$$\begin{equation} B_0^2 = \frac{32 \pi \exp (2) p_0}{9 \Omega _0^2 \tau ^2 \beta _p}. \end{equation} \tag{ 26 }$$

In Gaussian cgs units, the Ohmic resistive energy dissipation is

$$\begin{equation} \frac{4 \pi \eta }{c^2} J^2 = \frac{\eta }{4\pi } \left(\nabla \times {\boldsymbol {B}} {}\right)^2 \,. \end{equation} \tag{ 27 }$$

Following the approximation that the magnetic field is dominated by the B_y component yields

$$\begin{equation} \frac{4 \pi \eta }{c^2} J^2 = \frac{2 \exp (2) p_0}{\beta _p \tau } \left(\frac{t}{\tau }\right)^2 \exp ({-2t/\tau }) \cos (kz)^2. \end{equation} \tag{ 28 }$$

Note that in the case of Ohmic resistivity, the magnetic energy dissipation is exactly 90° out of phase spatially with the magnetic energy. This is neither the case with ambipolar diffusion (which scales with ${\boldsymbol {J}} {}\times {\boldsymbol {B}} {}$), nor with numerical dissipation.

Between time t = 0 and t, Ohmic resistivity deposits an energy density

$$\begin{equation} E = \int _0^t dt^{\prime } \frac{\exp (2) p_0}{\beta _p \tau ^3} t^{\prime 2} \exp ({-2t^{\prime }/\tau }) \left(1+\cos (2kz)\right) \end{equation} \tag{ 29 }$$

into the gas. Using the relation for the thermal energy density e = p/(γ − 1), we can relate the total energy density E deposited resistively to the resulting temperature fluctuation to first order in $\beta _p^{-1}$: δT = [(γ − 1)E/p₀]T₀.

5.2.1. Adiabatic Limit

In the adiabatic limit, where the thermal diffusion coefficient μ ≪ η, and to first order in $\beta _p^{-1}$, there is neither advective nor diffusive energy transport. Therefore, the highest point-to-point temperature variation occurs after all the magnetic energy has resistively dissipated. Then

$$\begin{eqnarray} &&\delta T = 2 \int _0^{\infty } dt^{\prime } \frac{(\gamma -1) T_0 \exp (2)}{\beta _p \tau ^3} t^{\prime 2} \exp ({-2t^{\prime }/\tau }) \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &&= \frac{(\gamma -1) \exp (2) }{2\beta _p} T_0, \quad \frac{\mu }{\eta }\ll 1\;\;\,\quad\qquad \end{eqnarray} \tag{ 31 }$$

where the factor of sin (π/2) − sin (3π/2) = 2 comes from taking the difference between the hottest and coldest points.

5.2.2. Intermediate τ_E = τ/2 Case

An analytical solution is possible in the particular case that the thermal diffusion has a timescale of τ_E = 1/(4μk²) = τ/2, then Equation (29) gives

$$\begin{equation} \frac{\partial \delta T}{\partial t} =\frac{2 (\gamma -1) T_0 \exp (2)}{\beta _p \tau ^3} t^2 \exp ({-2t/\tau }) -\frac{\delta T}{\tau _E}. \end{equation} \tag{ 32 }$$

In computing this solution note that the temperature fluctuation has wavenumber 2k and Equation (32) can be solved analytically, becoming

$$\begin{equation} \delta T = \frac{ (\gamma -1) \exp (2) t^3 \exp ({-t/\tau _E})}{12 \tau _E^3 \beta _p} T_0, \end{equation} \tag{ 33 }$$

which has a maximum for t = 3τ_E = 1.5τ of

$$\begin{equation} \delta T = \frac{9(\gamma -1)}{4 \exp (1) \beta _p} T_0. \end{equation} \tag{ 34 }$$

This underestimates the final result because we have ignored the adiabatic expansion of the magnetic field bundle due to the magnetic pressure, which also generates temperature variations that are first order in $\beta _p^{-1}$ and are spatially in phase with the resistive heating. However, that signal would have modestly diffused away by t = 3τ_E.

5.2.3. Fast Cooling Limit

Finally, if the temperature diffusion is very fast, the largest temperature variation will occur at the time of fastest heating, or t = τ, when B = B_p, and the thermal energy perturbation will be just (4π/c²)ηJ²τ_E. This means

$$\begin{eqnarray} \delta T &=& \frac{2 (\gamma -1) }{\beta _p }\frac{\tau _E}{\tau } T_0 \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} \;\;\quad\qquad\qquad &=& \frac{(\gamma -1)}{2\beta _p} \frac{\eta }{\mu } T_0, \qquad \frac{\mu }{\eta }\gg 1. \end{eqnarray} \tag{ 36 }$$

Here, the heating simply scales by the ratio of magnetic diffusion to thermal diffusion η/μ.

In our simulation, we use temperature independent values of the resistivity η ∼ 9 × 10¹⁴ cm² s⁻¹ (Table 1), and the thermal diffusivity μ ∼ 1.8 × 10¹⁴ cm² s⁻¹ (Equation (8)), so

$$\begin{eqnarray} & \tau _E k^2 = \frac{1}{4 \mu } = 1.4\times 10^{-15} \mathrm{\ s\ cm^{-2}}, \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} &\tau k^2 = \frac{1}{\eta } = 1.1 \times 10^{-15} \mathrm{\ s\ cm^{-2}}. \end{eqnarray} \tag{ 38 }$$

This means that, for any k, the timescale τ_E ∼ τ, so we have slightly slower cooling than the regime of Equation (34), even though we are not in the $\beta _p^{-1} \ll 1$ regime. However, if we do apply Equation (34) to our simulations where β_p ≈ 1 in both geometries, then the predicted ΔT/T₀ = 0.42. This agrees reasonably well with the results obtained in Figures 4 and 5, given the differing β_p regimes, especially considering that, as discussed in Section 5.2.2, Equation (34) neglects the heating in the current sheet due to the adiabatic expansion of the magnetic field bundles. This latter will be particularly important when the system enters the β_p ≈ 1 regime during particularly energetic phases. Due to this, one would expect that Equation (34) underestimates the heating that occurs in our simulations, as is the case.

The generalizability of this result beyond the regime of our unstratified, net vertical field simulation does remain to be confirmed. In particular, including stratification could allow for the export of magnetic energy vertically through buoyancy; while net-vertical field MRI is relatively violent, so relaxing these constraints might lead to weaker current sheets, providing less heating.

We have further demonstrated that resolving the current sheet structure of the MRI is a major undertaking. Even though we plausibly resolved the time–volume averaged magnetic field strength and stresses with only 128³ simulations (about 26 zones/scale height), the heating terms were only beginning to be resolved at 256³ (or about 50 zones/scale height). Furthermore, this resolution requirement was for nearly the highest resistivity that would still allow the MRI with the fastest growing wavelength contained in the thin aspect ratio box height. Lower levels of resistivity will require even higher resolutions.

We employed a unstratified shearing box model, which limits the manner in which the large scale cooling of the disk can be taken into account. Within this framework, we chose the fastest consistent value for the cooling timescale. A vertically global model with full radiative transfer may have effectively slower thermal relaxation at the midplane, and therefore the temperature fluctuations produced in our model may be underestimates.

A significant difference between the case treated in our simulations and an actual protoplanetary disk is that we allow neither η nor κ to vary with temperature. Accordingly, because our simulations already show strong temperature fluctuations, we can state that any MRI simulations with temperature-independent resistivity and opacity in the temperature range we consider are not self-consistent. Our one-dimensional results on the short-circuit instability (McNally et al. 2013) suggest that substantially stronger fluctuations would occur if a temperature-dependent resistivity and opacity are included.

For our simulations to progress, we must choose initial conditions that are MRI unstable. Awkwardly, as the MRI develops it heats the box. Were we to combine physically motivated, temperature-dependent resistivity with an initial condition that was marginally unstable to the MRI, the rising temperatures and falling resistivities would take the MRI away from marginal instability. Indeed, initially MRI-stable regions have been seen to become unstable as the turbulent disk heats (Latter & Balbus 2012; Faure et al. 2014). Such radial variations cannot be captured by a local shearing box model, so the method used in this paper has a limited ability to treat the dynamically evolving region at the potentially migrating, certainly self-determining inner edge of the dead zone.

Our result further supports the hypothesis that magnetic dissipation in general, and the short-circuit instability in particular may be responsible for the deduced thermal histories of calcium–aluminum-rich inclusions (CAIs) and chondrules, as suggested by Hubbard et al. (2012) and McNally et al. (2013).

Several classes of igneous CAIs, particularly compact Type A, Type B, and Type C CAIs, appear to have been briefly heated back over their melting temperatures at some time after their initial condensation (Stolper & Paque 1986; Scott & Krot 2005). Being formed of the most refractory minerals, CAIs may trace the inner, warm, and hence MRI-active, regions of the protoplanetary disk where temperature fluctuations generated by current sheets are unavoidable. Chondrules are made of less refractory material than CAIs, and so their thermal histories are often less extreme than those of remelted CAIs (e.g., Ebel 2006). However, they still require heating to temperatures where thermal ionization is easily capable of allowing the MRI to act.

Even in regions where the MRI is active but where the ionization does not sharply increase with temperature as required for the short-circuit instability, the temperature fluctuations we have found could provide a route to annealing silicates at more modest temperatures than are required to melt chondrules or CAIs (Sargent et al. 2009).

The scale of the temperature variations seen suggests that significant hysteresis may be possible in MRI activity. Once MRI begins, it may be able to provide the ionization to sustain itself even if the equilibrium disk becomes too neutral for it to continue (Latter & Balbus 2012; Faure et al. 2014). The strongly time varying accretion rates of protoplanetary disks, for example as seen in FU Orionis type events (Hartmann & Kenyon 1996), predict, and can be explained by, the inner active zone extending radially outward (Armitage et al. 2001). If MRI activity is adequate to generate local hot regions with low enough resistivity to support the MRI, the MRI might be active even in regions where the median resistivity is too high to allow MRI.

This is a similar mechanism to that suggested by Inutsuka & Sano (2005), although we favor current-driven heating resulting in thermal ionization rather than relying on direct non-thermal ionization by electrons accelerated in the current sheets. However, the end result, that current sheets may have a significantly lower resistivity, was explored by Muranushi et al. (2012) (see also erratum Muranushi et al. 2013) who found that some regions of what would otherwise be the dead zone can maintain an MRI active state by such a mechanism. Unlike the electron acceleration mechanism of Inutsuka & Sano (2005), resistively heated current sheets can be expected to be most effective at maintaining MRI when the disk is massive, as this regime has the lowest thermal diffusion.

Another consequence of temperature variations is the broadening of ice lines. If the background temperature varies with radius as R^−1/2, then 20% temperature variations will broaden the ice-line to ∼40% of its "mean" position. A significant amount of material would then be repeatedly evaporated and recondensed (Ros & Johansen 2013). This effect is significant even with quite small temperature variations: 4% temperature variations will still leave 8% broadening, producing an annulus with a width of a few local scale heights. The accretion crossing timescale for such an annulus exceeds α⁻¹, long enough for material within it to be processed and reprocessed. Multiple rounds of evaporation and condensation are unlikely to result in the same grain porosity, and hence collision resilience, as direct collisional growth, so larger grain sizes could be reached than those expected from collisional growth, as computed, for example, by Windmark et al. (2012).