On the Measurements of Numerical Viscosity and Resistivity in Eulerian MHD Codes

We measured the numerical shear and bulk viscosity of the Aenus code using sound waves. We set the background density and pressure to ${\rho }_{0}={p}_{0}=1$, and imposed a perturbation of the form

$$\begin{eqnarray}&&{v}_{1x}(x,t=0)=\epsilon \sin ({kx}),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\rho }_{1}(x,t=0)=\displaystyle \frac{{v}_{1x}(x,0)}{{c}_{{\rm{s}}}}{\rho }_{0},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{p}_{1}(x,t=0)=\displaystyle \frac{{v}_{1x}(x,0)}{{c}_{{\rm{s}}}}{\rm{\Gamma }}{p}_{0},\end{eqnarray} \tag{ 30 }$$

where ${c}_{{\rm{s}}}=\sqrt{{\rm{\Gamma }}{p}_{0}/{\rho }_{0}}$ is the sound speed. The amplitude of the velocity perturbation is set to $\epsilon ={10}^{-5}$, which is small enough to prevent wave steepening (see Shore 2007) within the time of our simulations.

In the presence of (numerical) viscosity, the wave is damped with time. For a plane wave ${v}_{x1}(x,t)={\hat{v}}_{x1}\exp [i({kx}-\omega t)]$, one finds from the dispersion relation

$$\begin{eqnarray}&&\omega =\displaystyle \frac{-i(4\nu /3+\xi ){k}^{2}}{2}\pm {{kc}}_{{\rm{s}}}\sqrt{1-\displaystyle \frac{{k}^{2}{\rho }_{0}{(4\nu /3+\xi )}^{2}}{4{\rm{\Gamma }}{p}_{0}}}.\end{eqnarray} \tag{ 31 }$$

In the weak damping approximation, i.e., if

$$\begin{eqnarray}&&\displaystyle \frac{{k}^{2}{\rho }_{0}\,{(4\nu /3+\xi )}^{2}}{4{\rm{\Gamma }}{p}_{0}}\ll 1,\end{eqnarray} \tag{ 32 }$$

the phase velocity remains constant and the solution can be written as

$$\begin{eqnarray}&&{v}_{x}(x,t)={\hat{v}}_{1x}{e}^{-{{\mathfrak{D}}}_{s}t}{e}^{{ik}(x\mp {c}_{{\rm{s}}}t)},\end{eqnarray} \tag{ 33 }$$

where the sound damping coefficient ${{\mathfrak{D}}}_{s}$ is defined as

$$\begin{eqnarray}&&{{\mathfrak{D}}}_{s}=\displaystyle \frac{{k}^{2}}{2}\left(\displaystyle \frac{4}{3}\nu +\xi \right).\end{eqnarray} \tag{ 34 }$$

The sound wave propagates with a constant speed and its amplitude decreases with time. Performing simulations with different values of the physical shear and bulk viscosity, we found an excellent agreement between the analytical (Equation (33)) and the numerical solution (see Rembiasz 2013, for details).

With simulation series #S1, #S3, #S5, and #S6 (Table 1; upper left panel of Figure 1), we investigated the influence of the reconstruction scheme on the numerical dissipation. To keep the contribution of the time integration errors as small as possible, we set ${C}_{\mathrm{CFL}}=0.01$. For every simulation, we measure the damping rate ${{\mathfrak{D}}}_{{\rm{s}}* }$ from the decay of the kinetic energy, from which we compute the numerical dissipation of the code according to Equation (34) as

$$\begin{eqnarray}&&\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }=\displaystyle \frac{2}{{k}^{2}}{{\mathfrak{D}}}_{{\rm{s}}* },\end{eqnarray} \tag{ 35 }$$

where the right-hand side is given by the simulation results. Thus, in the case of sound waves, one cannot determine ${\nu }_{* }$ and ${\xi }_{* }$ separately from the value of ${{\mathfrak{D}}}_{{\rm{s}}* }$, but only a linear combination of both quantities. For every simulation series (i.e., reconstruction scheme), we fit the function

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }\right)=r\mathrm{ln}({\rm{\Delta }}x)+d,\end{eqnarray} \tag{ 36 }$$

where r is the measured order of convergence of the scheme. From the fit parameter

$$\begin{eqnarray}&&d=\mathrm{ln}\left[\left(\displaystyle \frac{4}{3}{{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}\right)\displaystyle \frac{{ \mathcal V }}{{{ \mathcal L }}^{r-1}}\right],\end{eqnarray} \tag{ 37 }$$

we can compute $\tfrac{4}{3}{{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}$ if both the characteristic speed and the length of the system are known. As we will show below, for this test ${ \mathcal L }=\lambda =1$ ($\lambda =2\pi /k$ being the wavelength) and ${ \mathcal V }={c}_{{\rm{s}}}=1.291$. The results presented in Table 1 and the upper right panel of Figure 1 show that all schemes have an exponent r close to the theoretical order ${r}_{\mathrm{th}}$ of the method.

Zoom In Zoom Out Reset image size

The results of the simulation series #S2–#S4 (Table 1) show that the LF, HLL, and HLLD Riemann solvers damp sound waves very similarly.

With the simulation series #S7–#S9, we determine the contribution of the RK3 (#S7 and #S8) and RK4 (#S9) time integrators to the numerical dissipation. To keep the contribution of the spatial discretization errors as low as possible, we use the MP9 reconstruction. We vary the time step either by varying the grid resolution (keeping ${C}_{\mathrm{CFL}}=0.5;$ series #S7 and #S9) or by varying the CFL factor (series #S8). According to Equation (19), both approaches should be equivalent. In both cases, the RK3 scheme performs very close to third-order accuracy, whereas the order of the RK4 integrator is higher than expected. We attribute the overperformance of RK4 in this test to the fact that it is not a TVD scheme since we have not computed the time-reversed operator $\tilde{L}$ as suggested in the Shu & Osher (1988) and in Suresh & Huynh (1997) for time integration schemes with orders of larger than three. However, the overperformance of RK4 in this test is very likely a fortunate coincidence (see Section 4.1.3, where it is not the case). The estimators $\tfrac{4}{3}{{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}t}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}t}$ are obtained by employing a fitting procedure analogous to the one for $\tfrac{4}{3}{{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}$ (see Table 1 and the upper right panel of Figure 1).

The natural characteristic speed of this flow problem should be the sound speed (${ \mathcal V }={c}_{{\rm{s}}}$). To test this hypothesis, we performed simulations varying the background pressure (series #cS1) or the density (#cS2) within the range given in Table 2. The results were fitted with the function

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }\right)=\alpha \mathrm{ln}({c}_{{\rm{s}}})+d.\end{eqnarray} \tag{ 38 }$$

From the fit, we obtain the value of α, expecting $\alpha =1$, and from the offset d, we determine ${{\mathfrak{N}}}_{\mathrm{tot}}^{{\rm{\Delta }}x}={{\mathfrak{N}}}_{\mathrm{tot}}^{{\rm{\Delta }}x}=(4/3){{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}$ (Table 2). The table and the lower left panel of Figure 1 clearly show that the sound speed is indeed the characteristic speed of the system.

Table 2.  Wave Damping Simulations II

Series	Wave	p0	${\rho }_{0}$	b0	λ	${{\mathfrak{N}}}_{\ \mathrm{tot}}^{{\rm{\Delta }}x}$	α
#cS1	sound	$1...{10}^{4}$	1	0	1	46.2 ± 2.3	0.993 ± 0.003
#cS2	sound	1	${10}^{-4}...10$	0	1	46.2 ± 2.3	0.993 ± 0.003
#cS3	sound	1	1	0	$0.1...20$	46.2 ± 2.4	0.9899 ± 0.0024
#cA1	Alfvén	10−1	1	${10}^{-3}...20$	1	34.8 ± 2.4	0.945 ± 0.015
#cA2	Alfvén	$2\times {10}^{-3}...{10}^{7}$	1	1	1	34.8 ± 2.4	0.945 ± 0.015
#cA3	Alfvén	$2\times {10}^{-3}$	${10}^{-4}...{10}^{4}$	1	1	34.8 ± 2.4	0.945 ± 0.015
#cA4	Alfvén	$2\times {10}^{-3}$	1	1	$0.1...10$	44 ± 2	−1.0003 ± 0.0003
#cMS1	magnetosonic	1	1	${10}^{-4}...{10}^{3}$	1	40 ± 3	0.997 ± 0.006
#cMS2	magnetosonic	${10}^{-4}...{10}^{4}$	1	1	1	40 ± 3	0.997 ± 0.006
#cMS3	magnetosonic	1	${10}^{-3}...{10}^{4}$	1	1	40 ± 3	0.997 ± 0.006

Note. Simulations performed with the MP5 reconstruction scheme, the HLL Riemann solver, and the RK3 time integrator to identify the characteristic velocity and length of the system. The CFL factor was set to 0.01 to guarantee negligible time integration errors. The columns give (from left to right) the series identifier, the wave type, the initial pressure, density and magnetic field, and the wavelength λ. ${{\mathfrak{N}}}_{\ \mathrm{tot}}^{{\rm{\Delta }}x}$ is defined as in Table 1 for the corresponding wave type. The exponent α is obtained for each simulation series by means of the fitting functions given by Equations (38), (39), (46), (49), and (61), respectively.

Download table as:  ASCIITypeset image

To determine the characteristic length of the system (the natural candidate being the wavelength λ), we performed the simulation series #cS3 varying the wavelength λ (and the size of the simulation domain accordingly, i.e., $L=\lambda$). The results were fitted with the function

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }\right)=\alpha \mathrm{ln}(\lambda )+d,\end{eqnarray} \tag{ 39 }$$

expecting again $\alpha =1$. Table 2 and the lower right panel of Figure 1 confirm our hypothesis. The figure shows that the numerical viscosity term is proportional to ${k}^{-1}$, i.e., ${D}_{s* }/{k}^{2}\propto {k}^{-1}$ (see Equation (34)).

With the help of Alfvén wave simulations, we determine a linear combination of the numerical shear viscosity and resistivity of the code. We set the background magnetic field and density to ${b}_{0x}={\rho }_{0}=1$, the pressure to ${p}_{0}=2\times {10}^{-3}$, and the transversal velocity to ${v}_{0y}={v}_{0}=0$. We imposed a perturbation of the form

$$\begin{eqnarray}&&{b}_{1y}(x,0)=\epsilon \sin ({kx}),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{v}_{1y}(x,0)=-\displaystyle \frac{{b}_{y1}}{\sqrt{{\rho }_{0}}}.\end{eqnarray} \tag{ 41 }$$

In ideal MHD, an Alfvén wave propagates with a constant amplitude at the Alfvén speed ${c}_{{\rm{A}}}={b}_{0x}/\sqrt{{\rho }_{0}}$. In the presence of viscosity and resistivity, the wave amplitude decreases with time. In the weak damping approximation, i.e., for ${k}^{4}{(\nu +\eta )}^{2}/(4{c}_{{\rm{A}}}^{2})\ll 1$, the velocity evolution reads (for the derivation, see Campos 1999)

$$\begin{eqnarray}&&{v}_{y}(x,t)={v}_{0}{e}^{-{{\mathfrak{D}}}_{{\rm{A}}}t}{e}^{{ik}(x\mp {c}_{{\rm{A}}}t)},\end{eqnarray} \tag{ 42 }$$

where the Alfvén damping rate is defined as

$$\begin{eqnarray}&&{{\mathfrak{D}}}_{{\rm{A}}}=\displaystyle \frac{{k}^{2}}{2}(\eta +\nu ).\end{eqnarray} \tag{ 43 }$$

We verified Equation (43) with the help of numerical simulations, and also checked that the bulk viscosity does not influence the damping coefficient (see Rembiasz 2013, for details).

In the simulation series #A2, #A4, and #A5 (see Table 1), we compared the influence of the MP5, MP7, and MP9 reconstruction schemes on the numerical shear viscosity ${\nu }_{* }$ and resistivity ${\eta }_{* }$. For every simulation, we measured the decrease of the kinetic energy, from which we determined a linear combination of the numerical shear viscosity and resistivity

$$\begin{eqnarray}&&{\nu }_{* }+{\eta }_{* }=\displaystyle \frac{2}{{k}^{2}}{{\mathfrak{D}}}_{{\rm{A}}* }.\end{eqnarray} \tag{ 44 }$$

The simulation results are fitted with the function

$$\begin{eqnarray}&&\mathrm{ln}({\nu }_{* }+{\eta }_{* })=r\mathrm{ln}({\rm{\Delta }}x)+d,\end{eqnarray} \tag{ 45 }$$

where r is the numerically measured order of accuracy of the reconstruction scheme. From the fit parameter d and Equations (17), (21), and (44), we determined ${{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$.⁴ Table 1 and the upper left panel of Figure 2 show that all methods have an order of convergence close to the theoretical expectation.

Zoom In Zoom Out Reset image size

According to the results of the simulation series #A1, #A2, and #A3 (Table 1), the numerical dissipation of the LF, HLL, and HLLD Riemann solvers are also very similar for Alfvén waves.

With the simulation series #A6 and #A7 (upper right panel of Figure 2), we assessed the contribution to the numerical dissipation of the RK3 and RK4 time integrators, respectively. We set ${C}_{\mathrm{CFL}}=0.8$ and changed the time step by varying the grid resolution. The results presented in Table 1 and the upper right panel of Figure 2 show that the RK3 time integrator performs at its theoretical order, whereas the order of the RK4 integrator is once again (like in the sound wave tests) higher than expected.

The characteristic velocity ${ \mathcal V }$ for the Alfvén wave test problem can be inferred from simulation series #cA1, #cA2, and #cA3, in which we varied the magnetic field, pressure, and density, respectively. We find that the logarithm of the numerical dissipation determined from these simulation data can be fitted (see the lower left panel of Figure 2) by the function

$$\begin{eqnarray}&&\mathrm{ln}({\nu }_{* }+{\eta }_{* })=\alpha \mathrm{ln}({c}_{\mathrm{ms}})+d,\end{eqnarray} \tag{ 46 }$$

where α and d are fitting parameters, and ${c}_{\mathrm{ms}}$ is the fast magnetosonic speed, which is defined as

$$\begin{eqnarray}&&{c}_{\mathrm{ms}}=\sqrt{\displaystyle \frac{1}{2}({c}_{{\rm{A}}}^{2}+{c}_{{\rm{s}}}^{2}+\sqrt{{({c}_{{\rm{A}}}^{2}+{c}_{{\rm{s}}}^{2})}^{2}-4{c}_{{\rm{A}}}^{2}{c}_{{\rm{s}}}^{2}{\cos }^{2}\theta })},\end{eqnarray} \tag{ 47 }$$

where θ is the angle between the perturbation wavevector and the background magnetic field. For a wavevector parallel to the background field ($\theta =0$),

$$\begin{eqnarray}&&{c}_{\mathrm{ms}}=\max \{{c}_{{\rm{A}}},{c}_{{\rm{s}}}\}.\end{eqnarray} \tag{ 48 }$$

The values of the fitting parameter α, which are given in Table 2, confirm that the characteristic velocity ${ \mathcal V }$ is the fast magnetosonic speed (not the Alfvén speed, as one could have presumed), both in the flow regime, where ${ \mathcal V }$ is dominated by the Alfvén speed and the sound speed.

The simulation series #cA4 (Table 2) shows that the characteristic length of the Alfvén wave simulations is, as for the sound wave test, the wavelength, and that the numerical dissipation can be fitted by (see Table 2 and the lower right panel of Figure 2)

$$\begin{eqnarray}&&\mathrm{ln}({\nu }_{* }+{\eta }_{* })=\alpha \mathrm{ln}\lambda +d.\end{eqnarray} \tag{ 49 }$$

Finally, to investigate whether the errors resulting from spatial discretization and time integration are additive, we performed simulations #A8, in which both types of errors should contribute non-negligibly to the numerical dissipation. Figure 3 shows the simulation results together with the expected numerical dissipation of the RK3 integrator (green), the MP5 scheme (blue), and the sum of both contributions (red). As the figure shows, the errors add linearly.

Zoom In Zoom Out Reset image size

From the simulations of magnetosonic waves, we determine the numerical resistivity and viscosity of the Aenus code. If not otherwise stated, the background pressure, density, and magnetic field strength are set to ${p}_{0}={\rho }_{0}={b}_{0y}=1$ and ${{\boldsymbol{b}}}_{0}={b}_{0y}\hat{{\boldsymbol{y}}}$, respectively. We perturb the background by a magnetosonic wave of the form

$$\begin{eqnarray}&&{v}_{x1}(x,0)=\epsilon \sin ({k}_{x}x),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{v}_{y1}(x,0)={v}_{x1}\displaystyle \frac{{k}_{x}^{2}{b}_{0x}{b}_{0y}}{{b}_{0x}^{2}{k}_{x}^{2}-{\rho }_{0}{\omega }^{2}},\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&{b}_{y1}(x,0)={v}_{x1}\displaystyle \frac{{k}_{x}{b}_{0y}\omega {\rho }_{0}}{{\rho }_{0}{\omega }^{2}-{b}_{0x}^{2}{k}_{x}^{2}},\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&{\rho }_{1}(x,0)={v}_{x1}\displaystyle \frac{{k}_{x}{\rho }_{0}}{\omega },\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&{e}_{1}(x,0)={v}_{x1}\displaystyle \frac{{k}_{x}{p}_{0}{\rm{\Gamma }}}{\omega ({\rm{\Gamma }}-1)},\end{eqnarray} \tag{ 54 }$$

where e₁ is the total specific energy of the wave. The velocity amplitude $\epsilon ={10}^{-5}$, and the wave's angular frequency ω is given by

$$\begin{eqnarray}&&{\omega }^{2}={k}^{2}{c}_{\mathrm{ms}}^{2}.\end{eqnarray} \tag{ 55 }$$

For $\theta =\pi /2$ (a value chosen in almost all simulations), the magnetosonic speed reads (see Equation (47))

$$\begin{eqnarray}&&{c}_{\mathrm{ms}}=\displaystyle \frac{\omega }{k}=\sqrt{\displaystyle \frac{{{{\boldsymbol{b}}}_{0}}^{2}+{\rm{\Gamma }}{p}_{0}}{{\rho }_{0}}}.\end{eqnarray} \tag{ 56 }$$

In the presence of viscosity or resistivity, the wave will be damped with time, i.e., the x component of the wave velocity will decrease as

$$\begin{eqnarray}&&{v}_{x}(x,t)=\epsilon \,{e}^{-{{\mathfrak{D}}}_{\mathrm{ms}}t}{e}^{{ik}(x\mp {c}_{\mathrm{ms}}t)},\end{eqnarray} \tag{ 57 }$$

where the damping coefficient for a fast magnetosonic wave propagating in the direction perpendicular to the background magnetic field is (for the derivation, see Campos 1999)

$$\begin{eqnarray}&&{{\mathfrak{D}}}_{\mathrm{ms}}=\displaystyle \frac{{k}^{2}}{2}\left(\displaystyle \frac{4}{3}\nu +\xi +\displaystyle \frac{\eta }{1+{c}_{{\rm{s}}}^{2}/{c}_{{\rm{A}}}^{2}}\right).\end{eqnarray} \tag{ 58 }$$

We verified this equation numerically (see Rembiasz 2013 for details).

We also performed simulations #MS1, #MS2, and #MS3 (Table 1) to investigate the influence of the MP5, MP7, and MP9 reconstruction schemes. From the measured kinetic energy damping, we determined a linear combination of the numerical resistivity, shear viscosity, and bulk viscosity, i.e., 

$$\begin{eqnarray}&&\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }+\displaystyle \frac{3}{8}{\eta }_{* }=\displaystyle \frac{2}{{k}^{2}}{{\mathfrak{D}}}_{\mathrm{ms}* }.\end{eqnarray} \tag{ 59 }$$

We fitted the simulation results with the function

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }+\displaystyle \frac{3}{8}{\eta }_{* }\right)=r\mathrm{ln}({\rm{\Delta }}x)+d,\end{eqnarray} \tag{ 60 }$$

where the fit parameter r is the numerically measured order of the reconstruction scheme. From the fit parameter d, and with Equations (17), (20), (21), and (59), we determined the combination of coefficients $\tfrac{4}{3}{{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}+\tfrac{8}{3}{{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$ (see Table 1 and the left panel of Figure 4).

Zoom In Zoom Out Reset image size

Using simulation series #MS4 and #MS5, we studied the contribution of the RK3 and RK4 time integrators to the numerical dissipation (see Table 1 and middle panel of Figure 4). We find that the RK4 integrator performs at a higher order than theoretically expected. Again, we point out that probably due to the non-TVD preserving property of our implementation of RK4, it overperformes in this test (see Section 4.1.2).

To determine the characteristic speed, we performed simulation series #cMS1, #cMS2, and #cMS3, varying background magnetic field strength, pressure, and density, respectively (Table 2). It is not surprising that the characteristic speed is the fast magnetosonic speed (see the bottom panel of Figure 4), which is confirmed quantitatively with the help of the fit function

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{4}{3}{\nu }_{* }+{\xi }_{* }+\displaystyle \frac{3}{8}{\eta }_{* }\right)=\alpha \mathrm{ln}({c}_{\mathrm{ms}})+d.\end{eqnarray} \tag{ 61 }$$

As expected, the fit (see Table 2) is consistent with $\alpha =1$ within the measurement errors. The value of d can be used to estimate $\tfrac{4}{3}{{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}+{{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}+\tfrac{3}{8}{{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$. In the asymptotic regime ${b}_{0}\ll {p}_{0}$, the numerical damping is independent of the magnetic field strength, while it is proportional to the field strength for ${b}_{0}\gg {p}_{0}$.

So far, we measured the numerical damping for three wave types separately. For each type of a wave, the damping coefficient depends on a linear combination of the resistivity, shear viscosity, and bulk viscosity (see Equations (34), (43), and (58)). This gives a system of three linearly independent equations with three unknowns, which has a unique solution.

If we consider a series of simulations in which time discretization errors are negligible (such as those in the upper left panels of Figures 1, 2, and 4), at a fixed grid resolution $k{\rm{\Delta }}x=\mathrm{const}.$ and for a given numerical method, the numerical viscosity and resistivity should be the same according to our ansatz (Equations (17), (20), and (21)). Therefore, the damping rates of the three propagating waves can (in principle) be used to compute the numerical viscosity and resistivity.

In Figure 5, we present the resistivity, shear viscosity, and bulk viscosity for three different reconstruction schemes (MP5, MP7, and MP9) as a function of grid resolution. In all simulations, we used the HLL Riemann solver, an RK4 time integrator, and ${C}_{\mathrm{CFL}}=0.01$. Fitting a power law to the data allows us to compute the coefficients entering in the ansatz given by Equations (17), (20), and (21). From the fit parameters, which can be found in Table 3, we compute the exponent r appearing in the ansatz independently for resistivity and viscosity. For all three reconstruction schemes (1) the values of r are close to the expected order of convergence, ${r}_{\mathrm{th}}$, except for a small deviation in the case of MP9, and (2) the value of the numerical viscosity is significantly larger than the absolute value of the numerical resistivity.

Zoom In Zoom Out Reset image size

Table 3.  Values of the Parameters in Our Ansatz for the Spatial Dependence of the Numerical Viscosity (Equations (17) and (20)) and Numerical Resistivity (Equation (21)), Based on a Fit of the Results Shown in Figure 5

Reconstruction	${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$	rη	${{\mathfrak{N}}}_{\nu }^{{\rm{\Delta }}x}$	rν	${{\mathfrak{N}}}_{\xi }^{{\rm{\Delta }}x}$	rξ
MP5	−7.0 ± 0.5	4.94 ± 0.02	49 ± 3	4.95 ± 0.02	51 ± 4	4.95 ± 0.02
MP7	−41 ± 3	6.81 ± 0.03	270 ± 50	6.80 ± 0.06	354 ± 18	6.881 ± 0.017
MP9	−3 ± 6	6.8 ± 0.6	300 ± 200	7.9 ± 0.3	200 ± 200	7.7 ± 0.3

Note. We obtain three different estimates of the exponent r from fits of resistivity (r_η), shear viscosity (r_ν), and bulk viscosity (r_ξ).

Download table as:  ASCIITypeset image

Our most striking result is that the value of the numerical resistivity ${\eta }_{* }$ is negative, its absolute value being about one order of magnitude smaller (and even two orders for MP9) than the value of the numerical viscosity (see Table 3). The value of the resistivity coefficient ${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$ obtained with MP9 (providing the most accurate result) is compatible with a non-negative (very small) numerical resistivity, while the values of the numerical shear viscosity and bulk viscosity are positive and very similar. The resulting damping rate, which is a combination of resistivity and viscosity, prevented an amplification of the wave amplitude in all three systems studied. Taken together, these facts suggest that the numerical viscosity must be considerably larger than the numerical resistivity of the code. Hence, we conjecture that there are large systematic uncertainties that prevent us from properly measuring the numerical resistivity of the code in all three wave propagation tests.

Given that the dissipation is dominated by viscosity rather than resistivity, we had to turn to a completely different system in order to study whether our ansatz for numerical resistivity is a valid one, and if true, whether the results are consistent with a positive value of the resistivity (see Section 4.3).

So far, in the wave damping simulations, we have set the background velocity to zero. To test whether a background velocity affects the numerical damping (i.e., by modifying the characteristic speed of the flow), we repeated the damping test for the sound waves, Alfvén, waves, and magnetosonic waves with a non-zero background velocity. All simulations were performed with the MP5 reconstruction scheme and the RK3 time integrator (with ${C}_{\mathrm{CFL}}=0.1$). For all three Riemann solvers and all types of waves, we observed the same behavior, i.e., the component of the background flow velocity that is perpendicular to the propagation of the wave (${v}_{0y}$) does not affect the numerical dissipation, whereas the parallel component (${v}_{0x}$) does. Figure 6 shows the numerical dissipation for some exemplary simulations of the sound wave damping test done with the HLL Riemann solver. Based on these simulations, we conclude that the characteristic speed of the flow is given by the sum of (the parallel component of) the flow velocity and the fast magnetosonic speed or, in other words, by the fast magnetosonic eigenvalue of the ideal MHD equations.

Zoom In Zoom Out Reset image size

In this section, we sketch how to analytically obtain a growth rate and an instability criterion of the TM instability, leaving out all technical details that can be found in Appendix A.1. Many of the results presented here were already obtained by FKR63, or they are different limits of expressions found in that work.

Consider a two-dimensional flow in the x–y plane of constant background density ${\rho }_{0}=1$ threaded by a magnetic field

$$\begin{eqnarray}&&{b}_{0x}={b}_{0}\tanh (y/a),\end{eqnarray} \tag{ 76 }$$

where b₀ is the magnetic field strength and a defines the shear length (see Figure 8). This magnetic field configuration gives rise to a current sheet at $| y/a| \lesssim 1$. To balance the resulting magnetic pressure gradient, one can introduce either a gas pressure gradient, so that ${{\rm{\nabla }}}_{y}(p+{b}_{0x}^{2}/2)=0$ (pressure equilibrium configuration), or an additional magnetic field component, so that ${{\rm{\nabla }}}_{y}({b}_{0x}^{2}/2+{b}_{0z}^{2}/2)=0$ (force-free configuration). Both equilibrium configurations are stable in ideal MHD, but are TM unstable in resistive MHD.

Zoom In Zoom Out Reset image size

FKR63 derived the instability criterion and the growth rate using the linearized resistive-viscous MHD equations in the incompressible limit, which read

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{b}}={\rm{\nabla }}\times ({{\boldsymbol{v}}}_{1}\times {{\boldsymbol{b}}}_{0})+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{b}},\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}{\rho }_{0}{\partial }_{t}{{\boldsymbol{v}}}_{1} & = & -{\rm{\nabla }}p+({\rm{\nabla }}\times {{\boldsymbol{b}}}_{1})\times {{\boldsymbol{b}}}_{0}+({\rm{\nabla }}\times {{\boldsymbol{b}}}_{0})\times {{\boldsymbol{b}}}_{1}\\ & & +{\rho }_{0}\nu {{\rm{\nabla }}}^{2}{{\boldsymbol{v}}}_{1},\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{v}}=0,\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{b}}=0,\end{eqnarray} \tag{ 80 }$$

where ${\boldsymbol{v}}={{\boldsymbol{v}}}_{0}+{{\boldsymbol{v}}}_{1}$, ${\boldsymbol{b}}={{\boldsymbol{b}}}_{0}+{{\boldsymbol{b}}}_{1}$, and we denote background and perturbed quantities with subscripts "0" and "1" respectively. In the incompressible limit, $| {{\boldsymbol{v}}}_{1}| \ll {c}_{s}$ holds, which was used to obtain the linearized equations. To simplify the notation, we omit hereafter the subscript "1" for the velocity perturbations, because the background flow is assumed to be at rest.

FKR63 solved the above equations using a WKB ansatz, i.e., 

$$\begin{eqnarray}&&{v}_{y}(x,y,t)=v(y){e}^{{ikx}+\gamma t},\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}&&{b}_{1y}(x,y,t)={b}_{1}(y){e}^{{ikx}+\gamma t},\end{eqnarray} \tag{ 82 }$$

where k is the wavevector in the x direction, and γ is the growth rate of the TM instability. This ansatz is justified only if the time dependence of the background magnetic field can be neglected. This is the case when the diffusion timescale is much larger than the instability timescale, i.e., ${a}^{2}/\eta \gg {\gamma }^{-1}$. The Alfvén crossing time must be sufficiently short too, i.e., $a/{c}_{{\rm{A}}}\ll {\gamma }^{-1}$, which is equivalent to considering instantaneous propagation of Alfvén waves through the system. Combining both conditions, we have

$$\begin{eqnarray}&&\displaystyle \frac{{a}^{2}}{\eta }\gg {\gamma }^{-1}\gg \displaystyle \frac{a}{{c}_{{\rm{A}}}}.\end{eqnarray} \tag{ 83 }$$

Among other cases, FKR63 also considered perturbations whose wavelengths in the x direction are comparable to (but smaller than) the shear width, i.e., 

$$\begin{eqnarray}&&k\lesssim {a}^{-1}.\end{eqnarray} \tag{ 84 }$$

For such perturbations, the wavevectors may differ from that of the fastest growing mode appreciably, and it is possible to set up a numerical test in which, for a given grid resolution, both the magnetic shear layer and the TM are well resolved.

FKR63 solved the TM problem in the limit (84) with a so-called boundary layer analysis (BLA; for details, see Appendix A.1). They define an inner region (see Figure 8), where resistive effects are important. We call this region the resistive layer of width ${L}_{{\rm{R}}}$ or, if the flow is also viscous, the resistive-viscous layer of width ${L}_{\mathrm{RV}}$. Far away from this layer, $| y| \gg {L}_{{\rm{R}}}$ or $| y| \gg {L}_{\mathrm{RV}}$, there is an outer region (see Figure 8), where resistivity can be ignored and the ideal MHD equations are valid. The layer width ${L}_{{\rm{R}}}$ or ${L}_{\mathrm{RV}}$ can be expressed in terms of the physical parameters of the system (see Appendix A). The velocity perturbation ${v}_{y}(x,y,t)$ of the WKB ansatz (Equation (81)) exhibits two extrema at $y=\pm {L}_{{\rm{R}}}$ or $y=\pm {L}_{\mathrm{RV}}$ (see Figure 8). Because the location of these extrema can be determined from our simulation results, the layer width ${L}_{{\rm{R}}}$ or ${L}_{\mathrm{RV}}$ is an appropriate quantity to compare simulation and theory.

In the BLA, the inner solution of the resistive MHD equations in a linearized background, i.e., 

$$\begin{eqnarray}&&{b}_{0x}(y)\approx {b}_{0}y/a.\end{eqnarray} \tag{ 85 }$$

(which holds for $| y| \ll a;$ blue dashed line in Figure 8), is matched with the ideal MHD solution in the outer region at some matching point $| {y}_{{\rm{m}}}|$. The coordinate ${y}_{{\rm{m}}}$ has to fulfill the condition ${L}_{{\rm{R}}}\ll | {y}_{{\rm{m}}}| \ll a$ (or ${L}_{\mathrm{RV}}\ll | {y}_{{\rm{m}}}| \ll a$), which implies that these transition regions can exist only if

$$\begin{eqnarray}&&{L}_{{\rm{R}}}\ll a\qquad \mathrm{or}\qquad {L}_{\mathrm{RV}}\ll a.\end{eqnarray} \tag{ 86 }$$

In the inviscid limit, analytic TM solutions of the resistive MHD equations can be obtained. For the background magnetic field given by Equation (76), TM will grow if

$$\begin{eqnarray}&&a\,k\lt 1\ \ (\mathrm{instability}).\end{eqnarray} \tag{ 87 }$$

In this case, the growth rate of the TM is

$$\begin{eqnarray}&&\gamma =1.82{\eta }^{3/5}{\left(\displaystyle \frac{{b}_{0}k}{\sqrt{{\rho }_{0}}}\right)}^{2/5}{a}^{-6/5}{\left(\displaystyle \frac{1}{ak}-ak\right)}^{4/5},\end{eqnarray} \tag{ 88 }$$

and the width of the resistive layer is according to our definition

$$\begin{eqnarray}&&{L}_{{\rm{R}}}=1.40{\eta }^{2/5}{\left(\displaystyle \frac{\sqrt{{\rho }_{0}}}{{b}_{0}k}\right)}^{2/5}{a}^{1/5}{\left(\displaystyle \frac{1}{ak}-ak\right)}^{1/5}\end{eqnarray} \tag{ 89 }$$

(see Equations (147) and (148), respectively).

The growth rate given by Equation (88) is of little value for our purpose, since it is obtained in the inviscid limit. However, as we have argued in Section 2, both numerical viscosity and resistivity are an unavoidable result of the discretization of the equations. Hence, if we want to use TMs to measure numerical resistivity, we have to use an approach that takes into account numerical viscosity as well. FKR63 also considered the resistive-viscous case for Prandtl numbers

$$\begin{eqnarray}&&{P}_{{\rm{m}}}\equiv \displaystyle \frac{\nu }{\eta }\gtrsim 1,\end{eqnarray} \tag{ 90 }$$

in the limit $k\ll {a}^{-1}$. Based on their approach, we obtained the TM growth rate for wavenumbers $k\lesssim {a}^{-1}$ (see Equation (153))

$$\begin{eqnarray}&&\gamma \approx 0.84{\eta }^{5/6}{\nu }^{-1/6}{\left(\displaystyle \frac{{b}_{0}k}{\sqrt{{\rho }_{0}}}\right)}^{1/3}{a}^{-4/3}\left(\displaystyle \frac{1}{a\,k}-a\,k\right),\end{eqnarray} \tag{ 91 }$$

and the width of the resistive-viscous layer

$$\begin{eqnarray}&&{L}_{\mathrm{RV}}\sim {(\eta \nu )}^{1/6}{\left(a\displaystyle \frac{\sqrt{{\rho }_{0}}}{{b}_{0}k}\right)}^{1/3}.\end{eqnarray} \tag{ 92 }$$

These expressions should be useful to set up a test to measure numerical resistivity. However, as we will show in the next sections, it is difficult to find a region in the numerical parameter space, where Equation (91) holds, i.e., where Equations (83), (84), (86), and (90) are fulfilled.

To demonstrate the possibility of using a TM simulation to measure numerical resistivity, we first set up the test using physical resistivity. This allows us to estimate the reconnection rate as a function of physical resistivity and viscosity.

Our numerical experiment is based on Landi et al. (2008), who were mainly interested in the nonlinear phase of TM, i.e., the formation of magnetic islands and the onset of turbulence. Since we want to study the exponential growth of a single TM in detail, we modified their setup for our purposes. We used a computational box of size $[-{L}_{x},{L}_{x}]\times [-{L}_{y},{L}_{y}]$, where ${L}_{x}={L}_{y}=\pi /3$, with periodic and open boundary conditions in x and y directions, respectively. We set the density and pressure to ${\rho }_{0}={p}_{0}=1$, and used an ideal-gas EOS with ${\rm{\Gamma }}=5/3$. In the expression for the background field, Equation (76), we set ${b}_{0}=1$ and a = 0.1.

We tested both the pressure equilibrium and force-free configurations and found that only the latter is suitable for our numerical experiments (see Rembiasz 2013 for details). To obtain the force-free configuration, we set

$$\begin{eqnarray}&&{b}_{0z}=\displaystyle \frac{{b}_{0}}{\cosh (y/a)}.\end{eqnarray} \tag{ 93 }$$

We note that our initial perturbations differ from those of Landi et al. (2008). Because those authors only perturbed the velocity, the TM instability is triggered promptly for high resistivities only ($\eta \geqslant {10}^{-5}$). Instead, we perturb both the velocity and the magnetic field based on an analytic solution of the TM:

$$\begin{eqnarray}&&{v}_{y}(x,y,t=0)=v(y)\sin ({kx}),\end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray}&&{b}_{1y}(x,y,t=0)={b}_{1}(y)\cos ({kx}).\end{eqnarray} \tag{ 95 }$$

The function $v(y)$ is given by Equations (132) and (145) for $| y| \geqslant {y}_{{\rm{m}}}$ and $| y| \leqslant {y}_{{\rm{m}}}$, respectively, where ${y}_{{\rm{m}}}$ is the matching point (typically ${y}_{{\rm{m}}}=2{L}_{\mathrm{RV}}$). Landi et al. (2015) used similar perturbations, i.e., eigenfunctions of the TM, in their studies of what they called "an ideal TM." This ideal case is a solution of the TM problem in a regime first studied by FKR63, but different from the one we consider here.

The function ${b}_{1}(y)$ in Equation (95) is given by Equation (131) for $| y| \geqslant {y}_{{\rm{m}}}$, and it is constant for $| y| \leqslant {y}_{{\rm{m}}}$, i.e., ${b}_{1}(y\leqslant {y}_{{\rm{m}}})={b}_{1}({y}_{{\rm{m}}})$ according to the "constant ψ approximation" (see Appendix A). The remaining perturbations ${v}_{1x}(x,y,t=0)$ and ${b}_{1x}(x,y,t=0)$ are determined from the divergence-free conditions ${\rm{\nabla }}\cdot {\boldsymbol{b}}={\rm{\nabla }}\cdot {\boldsymbol{v}}=0$. To reduce the computational cost, we chose the value of k in such a way that exactly one TM fits into the box, i.e., k = 3.

To compare the results of our TM simulations with the analytical predictions of Equations (91) and (92) for the TM growth rate and the width of the resistive-viscous layer, respectively, we must ensure that we are in the regime of applicability of these analytical predictions, i.e., Equations (83),(84), (86), and (90) should be fulfilled. The first condition (Equation (84)) is ensured by our choice of k and a. The other conditions can be written as

$$\begin{eqnarray}&&{{ \mathcal C }}_{1}\equiv \gamma {a}^{2}{\eta }^{-1}\gg 1,\end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray}&&{{ \mathcal C }}_{2}\equiv {a}^{-1}{c}_{{\rm{A}}}{\gamma }^{-1}\gg 1,\end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray}&&{{ \mathcal C }}_{3}\equiv {{aL}}_{\mathrm{RV}}^{-1}\gg 1,\end{eqnarray} \tag{ 98 }$$

$$\begin{eqnarray}&&{P}_{{\rm{m}}}\equiv \nu {\eta }^{-1}\gtrsim 1.\end{eqnarray} \tag{ 99 }$$

We plot iso-contours of these four quantities in the η-ν plane (see Figure 9) to locate the region where the analytical expressions are valid.

Zoom In Zoom Out Reset image size

We first discuss the results of a simulation with $\eta ={10}^{-5}$ and $\nu ={10}^{-4}$, which we call the reference model (#Rf) and which is marked by an asterisk in Figure 9. The first condition (Equation (96)) is only marginally satisfied for the reference model (${{ \mathcal C }}_{1}\approx 25$). To improve the situation, one should decrease the resistivity and viscosity, i.e., one should increase the grid resolution. This would place the model toward the lower left corner of Figure 9, where all the conditions are better satisfied. Therefore, we are limited here by the numerical resolution that we can afford. In the following, we present simulations with numerical resistivities and viscosities as low as 10⁻⁷, corresponding to values of ${{ \mathcal C }}_{1}$ in the range of $1\lt {{ \mathcal C }}_{1}\lt 100$, which are marginally consistent with Equation (96). As a result, the diffusion timescale is only about 10 times larger than the TM e-folding time, i.e., we observe diffusion of the background solution within the duration of the simulation. We circumvent this problem by solving, instead of the proper induction Equation (4), a modified (physically incorrect) version for a constant resistivity, η:

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{b}}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{b}})+\eta {{\rm{\nabla }}}^{2}({\boldsymbol{b}}-{{\boldsymbol{b}}}_{0}).\end{eqnarray} \tag{ 100 }$$

Thereby, the background magnetic field, ${{\boldsymbol{b}}}_{0}$, does not suffer from diffusion by resistivity.

The second condition (Equation (97)) yields ${{ \mathcal C }}_{2}\approx 403\gg 1$ for the reference model, i.e., the Alfvén crossing time is sufficiently small compared to the growth timescale of the TM instability.

The third condition (Equation (98)) is ${{ \mathcal C }}_{3}\approx 10$ for the reference model, i.e., it is only roughly fulfilled. This condition is the most challenging one to be met in numerical simulations, because the size of the resistive-viscous layer ${L}_{\mathrm{RV}}$ has to be much smaller than the width a of the magnetic shear layer. This can be achieved again by decreasing viscosity and resistivity, but since ${L}_{\mathrm{RV}}\propto {(\eta \nu )}^{-1/6}$ (Equation (92)) it is necessary to decrease $\eta \nu$ by six orders of magnitude to decrease ${L}_{\mathrm{RV}}$ by a factor of 10. Thus, if we aim for ${{ \mathcal C }}_{3}\approx 100$, we need $\sim {10}^{4}$ grid points for each box dimension to resolve the resistive-viscous layer with ∼10 grid points.

The fourth condition (Equation (99)) is satisfactorily fulfilled, since ${P}_{{\rm{m}}}=10\gtrsim 1$.

The reference simulation #Rf was performed with the HLL Riemann solver, an MP9 reconstruction scheme, and with a grid of 2048 × 2048 zones (Figure 10). We find that the initially imposed magnetic field and velocity perturbations do not evolve much with time, except for a growth of their amplitudes. This indicates that our initial perturbations, which are based on the TM solution in resistive-non-viscous MHD (in the constant ψ approximation), are very similar to the eigenfunctions of resistive-viscous TM.

The upper right panel of Figure 10 shows profiles in the y direction of the initial (black) and the evolved (at t = 40; green) magnetic field perturbations at $x=-0.5$, the latter being normalized to the ratio $\max | {b}_{y}(t=0)| /\max | {b}_{y}(t=100)|$. The corresponding velocity perturbation at x = 0 (bottom panels) exhibit two pronounced extrema surrounding the magnetic shear layer (marked by the two vertical green lines in the lower right panel, which is a zoom of the lower left panel), which are characteristic of TMs.

To measure the TM growth rate, we compute the evolution with time of the quantity

$$\begin{eqnarray}&&{\bar{E}}_{\mathrm{mag},y}\equiv {\int }_{-{L}_{x}}^{{L}_{x}}{\int }_{-a}^{a}\displaystyle \frac{{b}_{y}^{2}}{2}\ {dxdy},\end{eqnarray} \tag{ 101 }$$

where the integration is performed only up to $| y| \leqslant a$ to reduce a potential influence of boundary conditions. After an initial transient lasting up to 20 time units during which the initial perturbation adjusts to the analytic solution, ${\bar{E}}_{\mathrm{mag},y}(t)$ grows exponentially at a constant rate. Since ${b}_{y}\propto \exp (\gamma t)$,

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\mathrm{ln}({\bar{E}}_{\mathrm{mag},y})=\gamma t+{\rm{const.}},\end{eqnarray} \tag{ 102 }$$

where the constant depends on the initial perturbation amplitude and the box size. Using the above equation, we compute the instability growth rate by means of a simple linear regression. The black line in the upper left panel of Figure 10 shows the time evolution of ${\bar{E}}_{\mathrm{mag},y}$, while the green dashed line is the linear fit according to Equation (102). Note that both lines are almost indistinguishable after the initial transient time.

Zoom In Zoom Out Reset image size

To obtain the width of the resistive-viscous layer, we plot ${v}_{y}(x=0,y)$ for every simulation at t = 30 and measure the locations ${L}_{\mathrm{RV}}^{+}$ and ${L}_{\mathrm{RV}}^{-}$ of the two velocity extrema (see vertical green lines in the bottom right panel of Figure 10). To attribute a measurement error, we note that the extremum can be located anywhere inside of the corresponding computational zone of vertical size ${\rm{\Delta }}y$. Thus, the actual location of the extremum is uncertain up to an error $\pm {\rm{\Delta }}y/2$, i.e., the layer width is

$$\begin{eqnarray}&&{L}_{\mathrm{RV}}=\displaystyle \frac{{L}_{\mathrm{RV}}^{+}-{L}_{\mathrm{RV}}^{-}}{2}\pm \displaystyle \frac{{\rm{\Delta }}y}{2}.\end{eqnarray} \tag{ 103 }$$

The methodology explained above to measure the growth rate γ and the width of the resistive-viscous layer ${L}_{\mathrm{RV}}$ was applied to all TM simulations discussed below. To understand the dependence of these quantities on the different relevant parameters and to compare with the analytic results, we performed several series of simulations exploring the parameter space in the neighborhood of the reference model, by varying η, ν, b₀, and k. Details of these simulations can be found in Appendix A.

The main result extracted from this set of (numerically converged) simulations is the disagreement between the numerically obtained growth rates and the analytic ones given by Equation (91). The most likely explanation for the discrepancy is that the parameters of our TM simulations are outside of the regime of validity of the analytic results, particularly because of the difficulty to guarantee ${L}_{\mathrm{RV}}\ll a$. Unfortunately, this means that the analytic expressions (91) and (92) cannot be used to measure numerical resistivity. Thus, we decided to use an empirical approach to the problem.

Using the insight gained from the theoretical work of FKR63, we postulate an ansatz for the dependence of both the TM growth rate and the width of the resistive-viscous layer on the physical parameters, which we then calibrate using the series of numerical simulations mentioned above. The whole procedure is described in detail in Appendix A. The empirical expressions resulting for k = 3 and a = 0.1 are

$$\begin{eqnarray}&&\gamma (k=3,a=0.1)=34.56{\eta }^{4/5}{\nu }^{-1/5}{\left(\displaystyle \frac{{b}_{0}}{\sqrt{{\rho }_{0}}}\right)}^{2/5},\end{eqnarray} \tag{ 104 }$$

$$\begin{eqnarray}&&{L}_{\mathrm{RV}}(k=3,a=0.1)=0.634{(\eta \nu )}^{1/6}{\left(\displaystyle \frac{{b}_{0}}{\sqrt{{\rho }_{0}}}\right)}^{-1/3}.\end{eqnarray} \tag{ 105 }$$

Figure 11 shows the TM growth rates (upper panel) and the width of the resistive-viscous layer (lower panel) measured from two series of simulations done with a viscosity $\nu ={10}^{-4}$ (blue diamonds) and $\nu ={10}^{-5}$ (red asterisks) for different values of the resistivity. Solid lines represent the empirical expressions given by Equations (104) and (105), while the analytic results given by Equation (91) are plotted with dashed lines in the upper panel of Figure 11. The discrepancy in the growth rate between analytic and numerical results is obvious, whereas our empirical expression (105) for the width of the resistive-viscous layer is compatible with the analytical one (Equation (92)).

Zoom In Zoom Out Reset image size

With the knowledge acquired from the resistive-viscous simulations of the previous section, we can tackle the problem of estimating the numerical resistivity of the code. If we perform a simulation with $\eta =0$, the development of TM signals the presence of a non-zero numerical resistivity ${\eta }_{* }$, because TM are not present in ideal MHD.

For the numerical setup presented in the previous section and for a viscosity $\nu ={10}^{-4}$, the TM growth rate should be well described by Equation (104), if ${\eta }_{* }\lesssim {10}^{-5}$. In this case, we can determine ${\eta }_{* }$ using the expression

$$\begin{eqnarray}&&{\eta }_{* }={\left(\displaystyle \frac{\gamma (k=3,a=0.1)}{34.56}\right)}^{5/4}{\nu }^{1/4}{\left(\displaystyle \frac{\sqrt{{\rho }_{0}}}{{b}_{0}}\right)}^{1/2},\end{eqnarray} \tag{ 106 }$$

where we need to measure only the growth rate of the instability, γ, for a simulation with k = 3 and a = 0.1.

Alternatively, one could measure the resistive-viscous layer width (Equation (105)) to obtain ${\eta }_{* }$ from the expression

$$\begin{eqnarray}&&{\eta }_{* }={\left(\displaystyle \frac{{L}_{\mathrm{RV}}(k=3,a=0.1)}{0.634}\right)}^{6}{\nu }^{-1}{\left(\displaystyle \frac{{b}_{0}}{\sqrt{{\rho }_{0}}}\right)}^{2}.\end{eqnarray} \tag{ 107 }$$

This method is much less accurate, however, because measuring ${L}_{\mathrm{RV}}$ from a simulation is prone to rather large relative errors (of the order of 0.1), and because ${\eta }_{* }\propto {L}_{\mathrm{RV}}^{6}$.

To compute the numerical resistivity from Equation (106), we need to know the value of the viscosity ν. However, for a coarse numerical resolution, the value of the numerical viscosity can be of the same order. Therefore, we should require $\nu \gg {\nu }_{* }$. Expecting that the numerical resistivity and viscosity are of the same order, we had to choose a value of ν that is larger than the typical values of both numerical resistivity and numerical viscosity. On the other hand, ν must not be too large because the growth rate of the instability decreases with increasing ν, i.e., more expensive simulations are required. The size of the resistive-viscous layer also grows with ν and may become comparable to a, thus violating the condition ${L}_{\mathrm{RV}}\ll a$, i.e., Equation (104) no longer holds. As a compromise, we chose a value of $\nu ={10}^{-4}$ and performed all simulations with sufficiently high resolution to ensure that ${P}_{{\rm{m}}}\gg 1$.

Equation (104) was obtained from numerical simulations in which we removed the background field from the resistive term of the MHD equations (see Equation (100)) to prevent diffusion of the background field. The simulations to be discussed in the remainder of this section did not require this measure, because they were performed without physical resistivity. In spite of this difference, we can still apply the calibration obtained in the former series of simulations, because we find that the results of both series of simulations are consistent (TM develop in both cases, but the background magnetic field does not diffuse). Hence, numerical resistivity seems to act independently on large scales (diffusion of the background field across the box) and small scales (development of TM). This finding confirms our ansatz, which postulates that the value of numerical resistivity differs for phenomena occurring at different length scales, because ${\eta }_{* }$ depends on the typical length (${ \mathcal L }$) and velocity (${ \mathcal V }$) of the flow.

To determine the dependence of the numerical resistivity on the three Riemann solvers (LF, HLL, HLLD), we performed simulations with the MP5 reconstruction scheme, the RK3 time integrator, a CFL factor of 0.7, and grids of 128 × 128 to 1024 × 1024 zones. The default physical parameters were a = 0.1, k = 3, ${b}_{0}={\rho }_{0}=1$, $\nu ={10}^{-4}$, and $\eta =0$.

We find that TM are instigated by numerical resistivity for the LF solver. In the simulations performed with the HLL and HLLD Riemann solvers, no TM are observed, i.e., the numerical resistivity resulting from these solvers, though undetermined, must be much less than that of the LF solver. For the latter solver, as expected, the higher the grid resolution the smaller the instability growth rate, i.e., the lower the numerical resistivity (see Figure 12), since $\gamma \propto {\eta }^{4/5}$ (Equation (104)). Coarsening the grid resolution, the numerical resistivity eventually becomes so high that the width of the resistive-viscous layer is so large that Equation (104) is invalid, and we can no longer precisely measure the magnitude of the numerical resistivity. The resolution limit depends on the order of the reconstruction scheme, being 320, 256, and 224 zones per dimension for the MP5, MP7, and MP9 scheme, respectively.

Zoom In Zoom Out Reset image size

The results of an exemplary simulation (#Ex) without resistivity obtained with MP5 reconstruction scheme on a grid of 384 × 384 zones are shown in Figure 10. Like in the reference model #Rf (with $\eta ={10}^{-5};$ black dashed–dotted line in upper left panel), a TM grows exponentially with time in model #Ex (red dashed–dotted line in the panel), this time being driven by numerical resistivity (in this simulation ${\eta }_{* }=6.5\times {10}^{-6}$ marked with the third rightmost asterisk in Figure 12). The TM induced growth of the magnetic field (upper right panel of Figure 10) and velocity (bottom left panel) perturbations in model #Ex (without resistivity; red lines) are similar to those in model #Rf (with resistivity; green lines). This comparison clearly demonstrates that the behavior of the numerical resistivity closely resembles that of (real) physical resistivity.

We anticipate that the main contribution to the numerical resistivity comes from the y-direction. All variables exhibit a much stronger variation in the y-direction than in the x-direction. Hence, the characteristic length scales in our multi-D ansatz for the numerical resistivity (analogous one to ansatz (24) for numerical shear viscosity), are much larger in the x-direction than in the y-direction, ${{ \mathcal L }}_{x}\gg {{ \mathcal L }}_{y}$. More specifically, we can preliminarily estimate that ${{ \mathcal L }}_{x}\propto {k}^{-1}\propto {L}_{x}$ and ${{ \mathcal L }}_{y}\lesssim a$. Consequently, the total numerical errors coming from the x-direction will be negligible compared to the ones due to the discretization in y, allowing us to use the simpler one-dimensional ansatz in the following. Therefore, in the further discussion of TMs, we will refer to the 1D ansatz for numerical resistivity (Equation (21)) for the sake of simplicity. Furthermore, we will use ${\rm{\Delta }}x$ instead of ${\rm{\Delta }}y$, since they are equal in our simulations, having in mind, however, that the main contribution comes form errors proportional to ${({\rm{\Delta }}y/{{ \mathcal L }}_{y})}^{r}$. A similar situation occurred in the 2D simulation of sound waves (series #LS6–#LS11 from Table 4) in which ${({\rm{\Delta }}x/{{ \mathcal L }}_{x})}^{r}\gg {({\rm{\Delta }}y/{{ \mathcal L }}_{y})}^{r}$ and therefore the contribution of the y sweep to the numerical dissipation was negligible and one could equally well use 1D ansatzes for numerical dissipation.

The dependence of the numerical resistivity (which is determined from the measured growth rate of the instability using Equation (106)) on the grid resolution is shown in Figure 12. The results are fitted with the functions

$$\begin{eqnarray}\mathrm{ln}(\gamma ) & = & {\alpha }_{5}\mathrm{ln}({\rm{\Delta }}x)+{c}_{5},\\ \mathrm{ln}(\gamma ) & = & {\alpha }_{7}\mathrm{ln}({\rm{\Delta }}x)+{c}_{7},\\ \mathrm{ln}(\gamma ) & = & {\alpha }_{9}\mathrm{ln}({\rm{\Delta }}x)+{c}_{9},\end{eqnarray} \tag{ 108 }$$

where ${\alpha }_{i}$ and c_i are the coefficients of the MP reconstruction scheme of the ith order. Their values are (for a CFL factor in the range of 0.1–0.7)

$$\begin{eqnarray}{\alpha }_{5} & = & 2.35\pm 0.02,\quad {c}_{5}=9.85\pm 0.13,\\ {\alpha }_{7} & = & 2.740\pm 0.015,\quad {c}_{7}=11.06\pm 0.09,\\ {\alpha }_{9} & = & 2.88\pm 0.05,\quad {c}_{9}=11.1\pm 0.3.\end{eqnarray} \tag{ 109 }$$

If the numerical resistivity scales as ${\eta }_{* }\propto { \mathcal V }{ \mathcal L }\,{({\rm{\Delta }}x/{ \mathcal L })}^{r}$ as in Equation (21), one would naively expect that the growth rate scales as $\gamma \propto {({\rm{\Delta }}x)}^{4r/5}$, i.e., the expected theoretical values should be ${\alpha }_{5}=4$, ${\alpha }_{7}=5.6$, and ${\alpha }_{9}=7.2$, which do not agree with our results. As we explain below, this argumentation is wrong, however, because it fails to account for an (implicit) dependence of the quantities ${ \mathcal V }$ and ${ \mathcal L }$ on ${\eta }_{* }$.

To explain the apparent considerable reduction of the convergence order r of the MP reconstruction schemes in (109), we need to have a careful look at the ansatz (21) for the numerical resistivity (neglecting the contribution of time integration errors)

$$\begin{eqnarray}&&{\eta }_{* }={{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}{ \mathcal V }{ \mathcal L }{\left(\displaystyle \frac{{\rm{\Delta }}x}{{ \mathcal L }}\right)}^{r},\end{eqnarray} \tag{ 110 }$$

where ${ \mathcal V }$ and ${ \mathcal L }$ are the system's characteristic speed and length, respectively.

If we were to assume ${ \mathcal L }\propto a$ (which is constant), we would obtain ${r}_{i}=(5/4){\alpha }_{i}$. The conceptual mistake we have made here is that a is the correct choice for the characteristic length of the background magnetic field diffusion problem, but not for a TM whose length scale is much smaller than the shear width. It turns out (as we demonstrate below) that the characteristic length of the system is proportional to the width of the resistive-viscous layer, i.e., ${ \mathcal L }\propto {L}_{\mathrm{RV}}$. This seems logical because the current sheet can be described very well neglecting Ohmic dissipation everywhere outside the narrow resistive-viscous layer whose width is ${L}_{\mathrm{RV}}$ rather than a.

The value of ${ \mathcal L }$ is somewhat arbitrary, because the boundary of the resistive-viscous layer is (physically) not sharp. We defined its width to be set by the characteristic velocity peaks (see Figure 8), which is a useful convention for our purpose. In fact, there exists a transition region (marked in shaded yellow in the figure), where the ideal MHD equations can still approximately be applied, though one is already in the non-ideal regime.

For our applications, we found a useful definition based on the fact that resistive and viscous effects are largest in the vicinity of steep gradients of the MHD variables. The (in relative terms) most important gradient is that of the y-component of the velocity, which is very large between the two extrema close to the current sheet (see bottom left panel of Figure 10). Taking into account that ${L}_{\mathrm{RV}}$ is the half distance between the two extrema, which approximately corresponds to one-fourth of a wavelength of a sine function, we propose to use as the proper length scale

$$\begin{eqnarray}&&{ \mathcal L }=4{L}_{\mathrm{RV}}.\end{eqnarray} \tag{ 111 }$$

This choice is consistent with identifying ${ \mathcal L }$ with the wavelength in the wave damping tests. It further suggests that a similar reasoning based on local extrema may lead to the appropriate value in other systems as well. Combining Equations (111) and (110), we obtain

$$\begin{eqnarray}&&{\eta }_{* }={{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}{ \mathcal V }4{L}_{\mathrm{RV}}{\left(\displaystyle \frac{{\rm{\Delta }}x}{4{L}_{\mathrm{RV}}}\right)}^{r}.\end{eqnarray} \tag{ 112 }$$

On the other hand, from Equation (105), we have

$$\begin{eqnarray}&&{L}_{\mathrm{RV}}(k=3,a=0.1)=0.634{\eta }_{* }^{1/6}{\nu }^{1/6}{\left(\displaystyle \frac{{b}_{0}}{\sqrt{{\rho }_{0}}}\right)}^{-1/3}.\end{eqnarray} \tag{ 113 }$$

Note the explicit dependence of Equation (112) on ${L}_{\mathrm{RV}}$ and of Equation (113) on ${\eta }_{* }$. This dependence can be easily removed obtaining the expressions

$$\begin{eqnarray}&&{\eta }_{* }={\left[{2.536}^{(1-r)}{{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}{ \mathcal V }{\nu }^{(1-r)/6}{\left(\displaystyle \frac{{b}_{0}}{\sqrt{{\rho }_{0}}}\right)}^{-(1-r)/3}{({\rm{\Delta }}x)}^{r}\right]}^{6/(5+r)}\end{eqnarray} \tag{ 114 }$$

$$\begin{eqnarray}&&{L}_{\mathrm{RV}}={\left[\displaystyle \frac{0.2596}{{4}^{r}}\nu {\left(\displaystyle \frac{{b}_{0}}{\sqrt{{\rho }_{0}}}\right)}^{-2}{{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}{ \mathcal V }{({\rm{\Delta }}x)}^{r}\right]}^{1/(5+r)},\end{eqnarray} \tag{ 115 }$$

which are valid only for a = 0.1 and k = 3, and give the true dependence of ${\eta }_{* }$ and ${L}_{\mathrm{RV}}$ on the grid resolution. Consequently, the TM growth rate is expected to depend on ${\rm{\Delta }}x$ with an exponent ${\alpha }_{i}=24r/[5(5+r)]$, which allows us to compute the order of convergence from the numerical values ${\alpha }_{i}$ as

$$\begin{eqnarray}&&r=\displaystyle \frac{25{\alpha }_{i}}{24-5{\alpha }_{i}}.\end{eqnarray} \tag{ 116 }$$

Similarly, Equation (114) can be used to compute ${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$, resorting to the coefficient c_i from the fit to the growth rate and identifying ${ \mathcal V }$ as the magnetosonic speed, i.e., ${ \mathcal V }={c}_{\mathrm{ms}}\approx 1.63$ in this case ($\theta =\pi /2$ in Equation (47)).

In Table 5, we list the values of ${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$ and r computed with the procedure outlined above. The MP5 and MP7 schemes are almost fifth and seventh order accurate, whereas the MP9 scheme performs below the theoretical expectation. In other words, the higher the order of the reconstruction scheme, the higher the reduction of the convergence order. A possible explanation of this fact is the following. The function ${v}_{y}$ is proportional to ${y}^{-1}$ for $| y| \ll a$, i.e., outside of the resistive-viscous layer (see the discussion in Appendix A, in particular, Equation (129) and Figure 17). For this reason, all the derivatives of ${v}_{y}$ in the y-direction diverge. Thus, the neglected higher-order terms of the Taylor expansion in the reconstruction of ${v}_{y}$ for $y\approx 0$ can actually be dominant. Taking this into consideration, it is rather more surprising that the MP5 and MP7 schemes almost achieve their theoretical order of accuracy than the fact that the MP9 scheme performed below the theoretical expectation.

Table 5.  Resistivity Coefficient ${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$ (for the Definition, See Equation (21)) and Reconstruction Scheme Order r (Equation (116)) Determined from TM Simulations Performed with the LF Riemann Solver (see also Figure 15)

Reconstruction	${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$	r
MP5	16 ± 5	4.81 ± 0.09
MP7	142 ± 33	6.65 ± 0.08
MP9	170 ± 220	7.6 ± 0.6

Download table as:  ASCIITypeset image

According to Figure 13, the values of ${L}_{\mathrm{RV}}$ measured directly from the numerical simulations (see the previous section) agree well with those computed with Equation (115), where the values of r and ${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$ needed in this equation are extracted from the growth rate using Equations (109) and (116). This result shows that our assumptions are correct, which is far from being obvious, because we assumed that (1) numerical errors can be called "numerical resistivity," (2) this numerical resistivity can be treated as normal physical resistivity, and (iii) the same equations can be used to determine its magnitude or predict its influence on the system. Moreover, we also had to make use of ansatz (21) for the numerical resistivity.

Zoom In Zoom Out Reset image size

To test whether the magnetosonic speed is indeed the characteristic velocity of our TM setup, as expected from the arguments given in the discussion of the wave damping simulations, we ran several simulations with the same setup, but varying the fast magnetosonic speed from ${c}_{\mathrm{ms}}\approx 1$ to $\approx 39$. This was achieved by changing the background pressure from 0.01 to 900, keeping ${b}_{0}={\rho }_{0}=1$. The upper panel of Figure 14 shows that the numerical resistivity increases with ${c}_{\mathrm{ms}}$.

Zoom In Zoom Out Reset image size

Different from the wave damping tests, it is not straightforward, however, to compute the fast magnetosonic speed, because in TM simulations the perturbed fluid makes a "U-turn" in the vicinity of the magnetic shear layer (i.e., for $| y/a| \ll 1$). Therefore, determining the "correct" values of θ (which changes from 0 to $\pi /2$) and the background magnetic field strength (which changes from 1 for $| y/a| \gg 1$ to 0 for $| y/a| \approx 0$) is very error-prone. That is why we introduced ansatzes (17), (20), and (21) to have a simple way of estimating the code's numerical dissipation. Consequently, we took the maximum possible magnetosonic speed (obtained for $\cos \theta =0;$ Equation (47)) in our previous analysis, i.e., 

$$\begin{eqnarray}&&{c}_{\mathrm{ms}}=\sqrt{{c}_{{\rm{A}}}^{2}+{c}_{{\rm{s}}}^{2}}.\end{eqnarray} \tag{ 117 }$$

and put ${c}_{{\rm{A}}}=1$, which is well motivated by practical purposes (i.e., to keep the ansatzes as simple as possible). However, in the current analysis (simulations presented in Figure 14), we obtained a better fit with respect to ${c}_{\mathrm{ms}}$ with ${c}_{{\rm{A}}}=0.76$, which corresponds to the Alfvén speed at a distance of y = a, i.e., 

$$\begin{eqnarray}&&{c}_{\mathrm{ms}}=\sqrt{5/3{p}_{0}+{0.76}^{2}}.\end{eqnarray} \tag{ 118 }$$

We note that the precise choice of ${c}_{{\rm{A}}}$ is irrelevant in the high plasma β regime, while it only slightly affects the quality of the fits in the low plasma β regime (where the parameter $\beta \equiv {b}_{0}^{2}/(2{p}_{0})$).

From the measured growth rates γ, we determined the numerical resistivity in each simulation (${\eta }_{* }\propto {\gamma }^{5/4}$), and fitted the results with

$$\begin{eqnarray}&&\mathrm{ln}({\eta }_{* })=s\mathrm{ln}({c}_{\mathrm{ms}})+d\end{eqnarray} \tag{ 119 }$$

obtaining

$$\begin{eqnarray}s & = & 0.524\pm 0.002,\\ d & = & -\,12.9\pm 0.5.\end{eqnarray} \tag{ 120 }$$

From Equation (114) and Equation (104), we find that

$$\begin{eqnarray}&&{\eta }_{* }\propto {{ \mathcal V }}^{6/(5+r)},\end{eqnarray} \tag{ 121 }$$

and putting ${ \mathcal V }={c}_{\mathrm{ms}}$ in Equation (121), we finally obtain

$$\begin{eqnarray}&&s=\displaystyle \frac{6}{5+r}.\end{eqnarray} \tag{ 122 }$$

For MP5 reconstruction (${r}_{\mathrm{th}}=5$), the expected value is s = 0.6, which is close to the measured one. Using s from Equation (120), we determine the reconstruction scheme order to be $r=6.45\pm 0.05$, which is neither equal to 5 within the errors nor consistent with the value of $r=4.81\pm 0.09$ from Table 5. This discrepancy should not concern us, however, because we included only statistical errors in the measurement errors from the linear fit neglecting other errors, e.g., those originating from estimating the fast magnetosonic speed (which changes from zone to zone in the simulation). This implies that this way of determining the order of the reconstruction scheme is much less reliable than from the resolution studies.

In the bottom panel of Figure 14, we show the measured width of the resistive-viscous layer and its value (red curve) predicted from the fit Equation (119). Using the values of s and d in Equation (120), we determined r and ${{\mathfrak{N}}}_{\eta }^{{\rm{\Delta }}x}$, which we then inserted into Equation (115). This methodology demonstrates that our model is self-consistent, since the width of the resistive-viscous layer does depend on ${c}_{\mathrm{ms}}$ as expected.

At the beginning of this section, we made the assumption that the numerical resistivity is not changing the background magnetic profile, but affects only the flow in the resistive-viscous layer, where the TM grows. All consistency checks performed in this section seem to indicate that our assumption is correct. As a confirmation, we checked that there has not been a significant modification of the background profile in any of the simulations. This finding differs from the one obtained in the simulations with physical resistivity but without removing the background field from the induction equation (100). In those simulations, the background field started to diffuse during the simulations (which was the reason why we modified the induction equation (Equation (100)) in the first place).

Finally, in Figure 15, we present a comparison of the expected numerical dissipation in simulations of TM and magnetosonic wave damping based on our ansatz (see Equations (17), (20), and (21)) and the estimators from Tables 1 and 5. For the simulations, we set the characteristic velocities and lengths equal to one, i.e., ${ \mathcal V }={ \mathcal L }=1$. The box length is set to 1 too, hence "resolution" in the abscissa of Figure 15 refers to the number of zones per characteristic length. As we can see, the expected numerical dissipation based on calibration with the help of both types of simulations (MS waves and TM) is similar. This is an indication that our approach is presumably universal, and it makes us confident that it can be used to estimate dissipation coefficients for other flows.

Zoom In Zoom Out Reset image size

We consider perturbations of the system described in Section 4.3.1 whose wavelength in x direction is comparable to the shear width, i.e., $k\lesssim {a}^{-1}$ (Equation (84)). Note that one only needs to solve Equations (77) and (78) for ${v}_{y}$ and b_1y, respectively, while the other perturbation components can be easily determined from conditions (79) and (80). For the WKB ansatz (Equations (81) and (82)) to be justified, condition (83) must be satisfied.

Inserting the WKB ansatz into Equations (77) and (78) and taking the curl of the latter equation to eliminate the pressure, the y component of the induction equation and the derivative ${\partial }_{x}$ of the z component of the equation of motion read

$$\begin{eqnarray}&&\gamma {b}_{1y}={[{\rm{\nabla }}\times ({\boldsymbol{v}}\times {{\boldsymbol{b}}}_{0})]}_{y}+\eta (-{k}^{2}+{\partial }_{y}^{2}){b}_{1y},\end{eqnarray} \tag{ 124 }$$

$$\begin{eqnarray}&&\gamma {\rho }_{0}(-{k}^{2}+{\partial }_{y}^{2}){v}_{y}={\partial }_{x}\{{\rm{\nabla }}\times [({\rm{\nabla }}\times {{\boldsymbol{b}}}_{1})\\ &&\quad \times {{{\boldsymbol{b}}}_{0}+({\rm{\nabla }}\times {{\boldsymbol{b}}}_{0})\times {{\boldsymbol{b}}}_{1}+{\rho }_{0}\nu {{\rm{\nabla }}}^{2}{\boldsymbol{v}}]\}}_{z}.\end{eqnarray} \tag{ 125 }$$

After some algebra, we arrive at

$$\begin{eqnarray}&&\gamma {b}_{1y}={{ikv}}_{y}{b}_{0x}+\eta (-{k}^{2}+{\partial }_{y}^{2}){b}_{1y},\end{eqnarray} \tag{ 126 }$$

$$\begin{eqnarray}&&\gamma {\rho }_{0}(-{k}^{2}+{\partial }_{y}^{2}){v}_{y}={\rho }_{0}\nu ({k}^{4}-2{k}^{2}{\partial }_{x}^{2}+{\partial }_{x}^{4}){v}_{y}\\ &&\quad +{ik}[-{b}_{1y}{\partial }_{y}^{2}{b}_{0x}+{b}_{0x}(-{k}^{2}+{\partial }_{y}^{2}){b}_{1y}].\end{eqnarray} \tag{ 127 }$$

Because this system cannot be integrated analytically, FKR63 used the BLA method (see Section 4.3.1). They divided the domain into two regions, an outer one ($-{L}_{y}\leqslant y\lt -{y}_{\epsilon }$ and ${y}_{\epsilon }\lt y\leqslant {L}_{y}$, where ${y}_{\epsilon }$ is a small positive constant such that ${y}_{\epsilon }\ll a$) in which dissipative effects can be neglected, and an inner layer ($-{y}_{\epsilon }\lt y\lt {y}_{\epsilon }$) in which resistivity (and viscosity) are important (see Figure 8). The complete solution is determined by an interplay between both regions: resistivity acts in the inner region, but the rate at which magnetic field lines are reconnected also depends on the rate at which the field can be advected into and out of the inner region. Hence, the growth rate can be quite different for field profiles that are the same close to y = 0, but differ elsewhere.

Outer Layer—Condition (83) implies that there are negligible field gradients in the outer regions ($| y| \gt {y}_{\epsilon }$), i.e., resistive processes are slow there compared to the growth of TMs. On the other hand, $\gamma \gg {a}^{-2}\eta \sim {k}^{2}\eta$ (the second relation follows from assumption 84). Hence, to solve Equations (126) and (127) in the outer region the resistivity term in the induction Equation (126) can be neglected, i.e., 

$$\begin{eqnarray}&&{{ikv}}_{y}{b}_{0x}=\gamma {b}_{1y}.\end{eqnarray} \tag{ 128 }$$

Furthermore, from the second part of condition (83), we have $\gamma \ll {c}_{{\rm{A}}}/{L}_{y}\sim {c}_{{\rm{A}}}k$. This inequality together with Equation (128) implies that terms proportional to velocity (gradients) in Equation (127) are negligible, i.e., $| \gamma {\rho }_{0}{k}^{2}{v}_{y}| \ll | {{ik}}^{3}{b}_{0x}{b}_{1y}|$. Hence, TMs evolve so slowly that the plasma inertia (terms containing ${\rho }_{0}{v}_{y}$ in Equation (127)) can be neglected on the ideal MHD timescale, and Equations (126) and (127) simplify to

$$\begin{eqnarray}&&{{ikv}}_{y}=\displaystyle \frac{\gamma {b}_{1y}}{{b}_{0x}},\end{eqnarray} \tag{ 129 }$$

$$\begin{eqnarray}&&{ik}[-{b}_{1y}{\partial }_{y}^{2}{b}_{0x}+{b}_{0x}(-{k}^{2}+{\partial }_{y}^{2}){b}_{1y}]=0.\end{eqnarray} \tag{ 130 }$$

in the outer layers. So far, we have not yet made any assumption concerning the background magnetic field. For ${b}_{0x}(y)={b}_{0}\tanh (y/a)$, the solution of Equation (130) reads

$$\begin{eqnarray}&&{b}_{1y}(y)={b}_{1}{\left(\displaystyle \frac{1+\tanh (y/a)}{1-\tanh (y/a)}\right)}^{ak/2}\displaystyle \frac{\tanh (y/a)-{ak}}{{\rm{\Gamma }}(2-{ak})},\end{eqnarray} \tag{ 131 }$$

where b₁ is a constant (initial perturbation amplitude) and Γ denotes the Euler gamma function. The velocity perturbations can be easily determined combining Equations (129) and (131):

$$\begin{eqnarray}&&{v}_{1y}(y)=\displaystyle \frac{\gamma {b}_{1}}{{{ikb}}_{0}}{\left(\displaystyle \frac{1+\tanh (y/a)}{1-\tanh (y/a)}\right)}^{ak/2}\displaystyle \frac{\tanh (y/a)-{ak}}{{\rm{\Gamma }}(2-{ak})\tanh (y/a)}.\end{eqnarray} \tag{ 132 }$$

Note that Equation (129) has a singularity for $| y/a| \to 0$, i.e., for $| {b}_{0x}| \to 0$, which is removed (i.e., smearded out) by resistivity in the inner region.

Inner Layer—Resistive (and viscous) terms can no longer be neglected in the inner layer ($| y| \lt {y}_{\epsilon }$), and we have to solve Equations (126) and (127) simultaneously. Because in the inner region $| y/a| \ll 1$, we can approximate the background magnetic field (Equation (76)) as ${b}_{0x}(y)\approx {b}_{0}y/a$ (Equation (85)). In the inner region, perturbations in both velocity and magnetic field vary more strongly in the y direction than in the x direction, i.e., $| {k}^{2}{v}_{y}| \ll | {\partial }_{y}^{2}{v}_{y}|$ and $| {k}^{2}{b}_{1y}| \ll | {\partial }_{y}^{2}{b}_{1y}|$.⁵ Therefore, we can neglect the terms proportional to k² in Equations (126) and (127), and we obtain

$$\begin{eqnarray}&&\gamma {b}_{1y}={{ikv}}_{y}{b}_{0}{a}^{-1}y+\eta {\partial }_{y}^{2}{b}_{1y},\end{eqnarray} \tag{ 133 }$$

$$\begin{eqnarray}&&\gamma {\rho }_{0}{\partial }_{y}^{2}{v}_{y}={\rho }_{0}\nu (-2{k}^{2}{\partial }_{y}^{2}+{\partial }_{y}^{4}){v}_{y}+{{ikb}}_{0}{a}^{-1}y{\partial }_{y}^{2}{b}_{1y},\end{eqnarray} \tag{ 134 }$$

where we used also Equation (85). Combining both equations into a single one by eliminating terms that contain b_1y, we arrive at a sixth-order ordinary differential equation (ODE) for ${v}_{y}$ that reads

$$\begin{eqnarray}&&{v}_{y}^{(6)}[\nu \eta {y}^{2}]+{v}_{y}^{(5)}[-2\nu \eta y]+{v}_{y}^{(4)}[-(\gamma (\eta +\nu )\\ &&\quad +2\nu \eta {k}^{2}){y}^{2}+2\nu \eta ]+{v}_{y}^{(3)}[2\eta (\gamma +2\nu {k}^{2})y]\\ &&\quad +{v}_{y}^{(2)}[{k}^{2}{a}^{-2}{c}_{{\rm{A}}}^{2}{y}^{4}+\gamma (\gamma +2\nu {k}^{2}){y}^{2}-2\eta (\gamma +2\nu {k}^{2})]\\ &&\quad +{v}_{y}^{(1)}[2{k}^{2}{a}^{-2}{c}_{{\rm{A}}}^{2}{y}^{3}]=0,\end{eqnarray} \tag{ 135 }$$

where ${c}_{{\rm{A}}}^{2}\equiv {b}_{0}^{2}/{\rho }_{0}$ and ${v}_{y}^{(n)}\equiv {\partial }_{y}^{n}{v}_{y}$. This equation cannot be integrated analytically, i.e., some further approximations are necessary.

In the inner layer, but sufficiently far away from y = 0, velocity perturbations will have a solution of type ${v}_{y}\propto {y}^{-1}$, i.e., the terms with the highest order derivatives of ${v}_{y}$ dominate the solution of (135). Thus, we proceed neglecting lower-order derivatives in (135). Because the two terms with the highest order derivatives (${\partial }_{y}^{6}{v}_{y}$ and ${\partial }_{y}^{5}{v}_{y}$) depend on viscosity, we consider Equation (135) for two limiting cases, namely in the viscous and the inviscid regime.

In the inviscid case ($\nu =0$), Equation (135) reduces to a fourth-order differential equation that is still too complicated to be solved analytically. Therefore, we follow FKR63 and use their "constant ψ approximation".⁶

FKR63 noted that the function ${v}_{y}$ becomes singular in ideal MHD because ${v}_{y}\propto {y}^{-1}$ for $y\to 0$ (see Equation (129)), i.e., ${v}_{y}$ varies strongly in the limit of small resistivity. Resistivity regularizes this ill-behaved solution. The function b_1y varies less for $| y/a| \approx 0$ and can be approximated by a constant ${b}_{1y}(y)\approx {b}_{1y}(0)$. With this approximation, Equations (133) and (134) reduce to

$$\begin{eqnarray}&&\gamma \eta {\rho }_{0}{\partial }_{y}^{2}{v}_{y}-{k}^{2}{a}^{-2}{b}_{0}^{2}{v}_{y}{y}^{2}={ik}\gamma {b}_{0}{a}^{-1}{{yb}}_{1y}(0),\end{eqnarray} \tag{ 136 }$$

$$\begin{eqnarray}&&{\partial }_{y}^{2}{b}_{1y}=\displaystyle \frac{\gamma {\rho }_{0}}{{{ikb}}_{0}{a}^{-1}y}{\partial }_{y}^{2}{v}_{y}.\end{eqnarray} \tag{ 137 }$$

We solve this system of equations by first integrating Equation (136) to obtain ${v}_{y}$, which we insert into Equation (137) to get b_1y.

To express Equations (136) and (137) in dimensionless form, we introduce the dimensionless variables

$$\begin{eqnarray}&&s=y{\left(\displaystyle \frac{{k}^{2}{b}_{0}^{2}}{{a}^{2}\gamma \eta {\rho }_{0}}\right)}^{\tfrac{1}{4}}\equiv \displaystyle \frac{y}{{\epsilon }_{{\rm{R}}}},\end{eqnarray} \tag{ 138 }$$

$$\begin{eqnarray}&&{\rm{\Phi }}={{iv}}_{y}{\left(\displaystyle \frac{{b}_{0}^{2}{\rho }_{0}\eta {k}^{2}}{{a}^{2}{b}_{1y}^{4}(0){\gamma }^{3}}\right)}^{\tfrac{1}{4}},\end{eqnarray} \tag{ 139 }$$

$$\begin{eqnarray}&&\psi =\displaystyle \frac{{b}_{1y}}{{b}_{1y}(0)},\end{eqnarray} \tag{ 140 }$$

$$\begin{eqnarray}&&\lambda =\gamma {\left(\displaystyle \frac{{\rho }_{0}{a}^{2}}{{k}^{2}\eta {b}_{0}^{2}}\right)}^{\tfrac{1}{3}}\end{eqnarray} \tag{ 141 }$$

with the length scale

$$\begin{eqnarray}&&{\epsilon }_{{\rm{R}}}\equiv {\left(\displaystyle \frac{\gamma \eta {\rho }_{0}{a}^{2}}{{k}^{2}{b}_{0}^{2}}\right)}^{\tfrac{1}{4}}.\end{eqnarray} \tag{ 142 }$$

The physical interpretation of ${\epsilon }_{{\rm{R}}}$ is that resistive effects are important in the center (y = 0) up to $| y| \sim {\epsilon }_{{\rm{R}}}$. In these new variables, Equations (136) and (137) read

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\rm{\Phi }}}{{{ds}}^{2}}-{s}^{2}{\rm{\Phi }}=-s,\end{eqnarray} \tag{ 143 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}\psi }{{{ds}}^{2}}=-{\lambda }^{3/2}\displaystyle \frac{1}{s}\displaystyle \frac{{d}^{2}{\rm{\Phi }}}{{{ds}}^{2}}.\end{eqnarray} \tag{ 144 }$$

The solution of Equation (143) can be written as an integral over an auxiliary variable u (Goedbloed et al. 2010):

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \frac{s}{2}{\int }_{0}^{1}{(1-{u}^{2})}^{-1/4}{e}^{-{s}^{2}u/2}{du}.\end{eqnarray} \tag{ 145 }$$

The function Φ (black line in Figure 17) is positive for $s\gt 0$, and has a global maximum at ${s}_{\max }\approx 1.48$. Moreover, ${\rm{\Phi }}(s)\approx 1.06/\sqrt{\pi }\,s$ for $s\ll 1$, and ${\rm{\Phi }}(s)\approx 1/s$ for $s\gg 1$.

Zoom In Zoom Out Reset image size

To obtain the final form of the velocity perturbations, ${v}_{y}$, in the inner layer (given by Equations (139) and 145), we need to determine the TM growth rate γ. It can be calculated by matching the inner and outer solutions of ${v}_{y}$ (the former given by Equation (132)) at a certain point ${y}_{{\rm{m}}}$ in the region (marked in yellow in Figure 8), where both solutions are valid and overlap, i.e., where they give the same predictions for both the velocity and the magnetic field perturbations. The value of y_m must be large enough, so that ${\rm{\Phi }}(s)$ can be approximated as ${\rm{\Phi }}(s)\approx {s}^{-1}$, i.e., $s\gg 1$, yet small enough that the outer ideal MHD solution behaves like ${v}_{y}\propto {y}^{-1}$ (for $y\to 0$). Moreover, the inner resistive solution was found for such small values of y that ${b}_{0x}(y)$ can be approximated by Equation (85). Hence, ${y}_{{\rm{m}}}\ll a$ must hold. Recalling that s = 1 for $y={\epsilon }_{{\rm{R}}}$, we can combine the above conditions into

$$\begin{eqnarray}&&{\epsilon }_{{\rm{R}}}\ll {y}_{{\rm{m}}}\ll a.\end{eqnarray} \tag{ 146 }$$

The remaining part of the matching procedure is conceptually rather straightforward. Comparing Equations (145) and (132) at ${y}_{{\rm{m}}}$, we can determine the TM growth rate γ. We omit the details of these calculations⁷ and only give here the final expression for the TM growth rate in resistive (-inviscid) MHD:

$$\begin{eqnarray}&&\gamma ={\left(\displaystyle \frac{2}{2.12}\right)}^{4/5}{\eta }^{3/5}{\left(\displaystyle \frac{{b}_{0}k}{\sqrt{{\rho }_{0}}}\right)}^{2/5}{a}^{-6/5}{\left(\displaystyle \frac{1}{ak}-ak\right)}^{4/5}.\end{eqnarray} \tag{ 147 }$$

For $a\,k\gt 1$, the growth rate is complex (because of the last factor), i.e., the system is TM unstable only for perturbations with wavevectors $a\,k\lt 1$. On the other hand, for $a\,k\to 0$, the growth rate given by Equation (147) diverges, but in this regime the "constant ψ approximation" becomes invalid and Equation (147) no longer holds.

We note that the width of the resistive layer is somewhat arbitrary. Some authors (see Schnack 2009; Goedbloed et al. 2010) define that the layer extends to $| y| ={\epsilon }_{{\rm{R}}}$, whereas we define its boundary to be located at

$$\begin{eqnarray}&&| y| ={L}_{{\rm{R}}}=1.48{\epsilon }_{{\rm{R}}},\end{eqnarray} \tag{ 148 }$$

which corresponds to the maximum of ${\rm{\Phi }}(s)$ (Figure 17). This is a convenient definition because ${\rm{\Phi }}({s}_{\max })$ corresponds to the peaks of ${v}_{y}$, which can be determined from simulations rather easily (compare Figure 17 with the bottom panels of Figure 10).

Equations (132) and (145) for the velocity perturbation ${v}_{y}$, and Equation (131) for the magnetic field perturbation b_1y (supplemented by Equation (147) and the transformations (138) and (139) constitute a complete solution of the TM problem in resistive MHD. In the inner layer, the magnetic field perturbation b_1y is approximately constant, i.e., ${b}_{1y}(y\leqslant {y}_{{\rm{m}}})\approx {b}_{1y}({y}_{{\rm{m}}})$, and the perturbations ${v}_{x}$ and b_1x can be determined from conditions (79) and (80), respectively.

As a next step, we will generalize the above results to resistive-viscous MHD, because Equation (147) does not take into account (numerical) viscosity, which might be comparable to (or even larger than) numerical resistivity.

The derivation of the TM solution in resistive-viscous MHD by FKR63 is very similar to that of the non-viscous case. Therefore, we will only sketch their procedure (for more details, see FKR63). One integrates Equations (126) and (127) again separately in the outer and inner region. In the outer region, one can neglect dissipative terms and plasma inertia as in the inviscid case, i.e., the magnetic field perturbation b_1y and the velocity perturbation ${v}_{y}$ are still given by Equations (131) and (132), respectively. In the inner region, FKR63 simplified the sixth-order ODE (135) for ${v}_{y}$ using again the "constant ψ approximation" to a fourth-order ODE, but now including viscous terms. They further simplified the ODE by neglecting terms with lower-order derivatives, which is an acceptable approximation for magnetic Prandtl numbers ${P}_{{\rm{m}}}\equiv \nu /\eta \gtrsim 1$. Next, FKR63 introduced the dimensionless variables

$$\begin{eqnarray}&&\tilde{s}=\displaystyle \frac{y}{{\epsilon }_{\mathrm{RV}}},\end{eqnarray} \tag{ 149 }$$

$$\begin{eqnarray}&&\tilde{{\rm{\Phi }}}\propto {v}_{1y},\end{eqnarray} \tag{ 150 }$$

(where we used the tilde symbol to explicitly stress that functions 145 and 150 differ) with the length scale

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{RV}}={(\eta \nu )}^{1/6}{\left(\displaystyle \frac{a\sqrt{{\rho }_{0}}}{{b}_{0}k}\right)}^{1/3}.\end{eqnarray} \tag{ 151 }$$

The latter expression shows that both resistivity and viscosity affect the size of the region, where dissipative effects are important. We call this region the resistive-viscous layer. The width of this layer ${L}_{\mathrm{RV}}$ is proportional to ${\epsilon }_{\mathrm{RV}}$, and hence differs from the width ${L}_{{\rm{R}}}$ of the resistive (inviscid) layer, which is proportional to ${\epsilon }_{{\rm{R}}}$ (Equation (142)).

In the resistive-viscous case, matching condition (146) has to be replaced by

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{RV}}\ll {y}_{{\rm{m}}}\ll a,\end{eqnarray} \tag{ 152 }$$

and from matching of $\tilde{{\rm{\Phi }}}(\tilde{s})$ with Equation (132), we obtain the TM growth rate in resistive-viscous MHD:

$$\begin{eqnarray}&&\gamma \approx \displaystyle \frac{{2}^{4/3}}{3}{\eta }^{5/6}{\nu }^{-1/6}{\left(\displaystyle \frac{{b}_{0}k}{\sqrt{{\rho }_{0}}}\right)}^{1/3}{a}^{-4/3}\left(\displaystyle \frac{1}{a\,k}-a\,k\right).\end{eqnarray} \tag{ 153 }$$

The above expression differs from the result of FKR63 (see their Equation (H.8)) because these authors derived their equation in the $a\,k\ll 1$ limit (in our units), whereas we did not neglect terms proportional to $a\,k$. For the background magnetic field ${b}_{0x}={b}_{0}\tanh (y/a)$, we calculated the growth rate more accurately. Note that Equation (153) is only an approximation, because FKR63 only approximately solved ODE (135). Therefore, the above equation could be off by a small constant numerical factor.

Our numerical results differ from the analytic growth rates derived above. Restrictions of grid resolution and computing time prevented us from reaching the parameter regime where the derivation of the analytic growth rates holds. Therefore, in the following subsection, we present an "empirical" ansatz for the TM growth rate, which gives much better predictions for our simulation results presented in Section 4.3.