Energy spectrum of tearing mode turbulence in sheared background field
Abstract
The energy spectrum of tearing mode turbulence in a sheared background magnetic field is studied in this work. We consider the scenario where the nonlinear interaction of overlapping largescale modes stirs up a broad spectrum of smallscale modes, generating tearing mode turbulence. The spectrum of such turbulence is of interest since it is relevant to the smallscale backreaction to the largescale field. The turbulence we discuss here differs from traditional MHD turbulence mainly in two aspects, one being the existence of many linearly stable smallscale modes which cause an effective damping during energy cascade; the other being the scale independent anisotropy induced by the large scale modes tilting the sheared background field, as opposed to the scale dependent anisotropy frequently encountered in traditional weak turbulence or critically balanced turbulence theory. These two differences result in the deviation of energy spectrum from a simple power law behavior, taking the form of a power law multiplied by an exponential falloff. Numerical simulation is carried out using viscous resistive MHD equations to verify our theoretical prediction, and reasonable agreement is found between numerical result and our model.
While visiting at ]PPPL, Princeton, New Jersey
I Introduction
The generation of a spectrum of smallscale tearing modes or resistive kink modes by their largescale counterparts is a very relevant issue both in magnetically confined devices such as reverse field pinch (RFP) or tokamak, and in astrophysical plasma. For RFP, the constant interaction of tearing modes and resistive interchange modes keeps the plasma in a perpetual turbulent state RFP . For tokamak, the nonlinear excitation and overlapping of a spectrum of tearing modes or resistive kink modes may break the nested flux surfaces and lead to disruption Diamond1984 ; Strauss1986 ; Craddock1991 . In astrophysical plasma, the secondary plasmoid turbulence is found to play crucial role during magnetic reconnection both in kinetic Daughton2011 and in resistive MHDHuang16APJ investigations.
A very interesting aspect of these problems is the backreaction of smallscale field in a tearing turbulence to their largescale counterparts. A well known example of such backreaction is the hyperresistivity applied to the mean field, which has been the subject of intensive studies in the past decades Strauss1986 ; Craddock1991 ; Boozer1986 ; AB1986 ; Hameiri1987 . To understand this problem, however, knowledge regarding the structure of tearing turbulence spectrum is necessary Diamond1984 ; Strauss1986 ; Craddock1991 ; AB1986 ; Hameiri1987 . Hence, in this paper, we will try to construct a model describing the structure of an instability driven tearing turbulence spectrum in sheared strong magnetic field. This sheared and strongly magnetized case is most relevant to the laboratory plasmas, although it will also provide some insight to a more general consideration regarding any instability driven turbulence.
Two arguments are commonly invoked when studying the spectrum of MHD turbulence. One is the inertial range argument, which states that there exists a selfsimilar region in space between the stirring scale and dissipation scale where energy is conservatively transferred from one scale to another, directly result in a power law FrischBook ; DiamondBook . The other is the scale dependent anisotropy which indicates the ratio between the parallel and perpendicular length scale depends on . For weak turbulence in homogeneous magnetic field, threewave interaction results in no cascade along the parallel direction Ng96ApJ ; Ng96POP ; Galtier2000 ; Lithwick2003 . Hence is independent of , yielding the energy spectrum . For strong turbulence, assuming no scale dependent alignment, the frequently invoked critical balance assumption states that the nonlinear term and linear term are of the same order, where is the Alfvén speed of background field, is the velocity at a given perpendicular scale GS1995 ; GS1997 . This critical balance combined with the inertial range argument will yield the scale dependent anisotropy , corresponding to the energy spectrum DiamondBook . With scale dependent alignment, the balance between linear and nonlinear terms becomes , leading to the anisotropy , and energy spectrum Boldyrev2006 .
However, recent development in kinetic turbulence theory has pointed out the possibility that stable eigenmodes which are nonlinearly stirred up by unstable modes can act as an effective damping mechanism HatchPRL2011 ; HatchPOP2011 . This is equally true for tearing turbulence we are concerned with here. Unlike the externally driven turbulence in a homogeneous system commonly discussed, instability driven turbulence usually has a bunch of stable modes along with a few unstable modes which provide energy for the rest of the spectrum. The effective damping caused by those stable modes will break the previous conservative transfer of energy between scales and thus alter the structure of spectrum. As mentioned above, it can be expected that the resulting spectrum will deviate from the traditional power law structure , and take the form of a power law multiplied by an exponential fall Terry2009 ; Terry2012 . Here, , and are constant coefficients. Furthermore, recent resistive MHD simulation concerning plasmoid turbulence in sheared magnetic field has found discrepancy from the scale dependent anisotropy picture, and produced a somewhat scale independent anisotropy for strong turbulence cases where perturbation strength is comparable with the background field Huang16APJ , raising doubt on the scale dependence of anisotropy for tearing mode turbulence in a magnetically sheared system.
In this paper, we will study the impact to tearing turbulence spectrum caused by the violation of the scale dependent anisotropy assumption as well as the inertia range argument. On one hand, the presence of large scale perturbations in a sheared guide field is found to introduce a scale independent anisotropy in the small scale eddies. On the other hand, the effective damping of the turbulence will be calculated from linear stability of high modes, wherein the effective damping scales like , with and in the inviscid and viscous regime respectively. This power dependence is considerably lower than that of the classical dissipation, which generally scales like . An analytical model will be carried out for turbulence under such scale independent anisotropy and effective damping. Based on this model, the modified spectrum will be obtained by considering the local energy budget in space. This analytical spectrum will then be compared with resistive MHD simulation. Reasonable agreement is found between analytical predictions and numerical results.
The rest of the paper is arranged as follows. In Section II, the system of interest will be described and the basic resistive MHD equations will be introduced. In Section III, the theoretical model regarding the damped tearing turbulence will be given, as well as the modified turbulence spectrum. This new spectrum will be checked with simulation results in Section IV, and spectral properties and structure functions of the turbulence will be discussed. The issue of scale independence of the turbulence will be studied analytically as well as numerically. Further, this scale independent anisotropy will be checked for strong turbulence cases. Discussion about the implication of this new form of spectrum to future studies and a conclusion will be presented in Section V.
Ii System of interest
We will consider the standard compressible MHD equation with viscosity and resistivity included.
(1) 
(2) 
(3) 
(4) 
Here, Eq. (1) is the continuity equation, Eq. (2) is the equation of motion, Eq. (3) represents the equation of state and Eq. (4) is the Ohm’s law. is the plasma density, is the velocity, is the total magnetic field, is the current density and is the pressure. The vacuum permeability has been absorbed into and . Further, here is the adiabatic index, and should not be confused with the growth rate of the modes. The constant dissipation coefficients and stand for classic viscosity and resistivity respectively. An additional artificial fourthorder dissipation is also implemented to damp smallscale fluctuations at grid size. This should not be confused with the real hyperdissipation selfconsistently generated by the nonlinear terms Strauss1986 ; Craddock1991 .
In this study, we will consider a simple slab system with coordinates , and the boundary conditions are assumed to be periodic at all sides. The size of the system is , and , with the geometric center of the system set as . The three components of the magnetic field are as follows:
(5) 
Here, and are constants to be specified later. The system is initially in force free, with the pressure assumed to be constant and set equal to unity. The corresponding initial current profile is then:
(6) 
An artificial constant electric field along direction is implemented to sustain the initial current profile against resistive diffusion. Although the assumed global geometry is simple, it is sufficient to capture the fundamental physical process of the dynamics of smallscale tearing fluctuations, and the qualitative features of the theory are not expected to change in more realistic global geometry.
We can define a “safety factor” analogous to that of tokamak:
(7) 
and the “rotational transform”:
(8) 
The corresponding profile is then a function of . In the region , the profile is shown in Fig. 1, with , , , , and corresponding minimum safety factor . Numerical observation indicates that several largescale modes, such as , and modes, are unstable for this magnetic profile. The nonlinear growth and interaction of these modes will generate a spectrum of smallscale modes.

As the turbulence grows in strength, it will have a backreaction on the mean background field, leading to selfconsistent evolution of the background fields. The mean current profile will tend to relax under the turbulence spreading Diamond1984 ; Strauss1986 ; Craddock1991 , and it is observed that substantial profile flattening would occur over time after the turbulence has been fully established. Ultimately, the relaxation would reach a point that there is no free energy anymore, and the tearing turbulence would then gradually decay away. However, it would be shown in Section IV.1 that the characteristic time scale of such relaxation is much longer than the slowest nonlinear turnover time of eddies, thus the turbulence can be viewed as attaining a quasisteadystate within a certain time period.
Iii Analytical model for tearing turbulence
The structure of tearing turbulence spectrum will be discussed in this section. In our model, we consider a sheared background field, and the tilting of field lines caused by the presence of large scale perturbations. It will be seen that the small scale turbulence exhibits scale independent anisotropy under such perturbed sheared magnetic field. We further assume that the energy transfer across different scales is local in space as it is in traditional turbulence theories FrischBook ; DiamondBook , thus the energy flow from largescale unstable modes to medium and smallscale stable modes can be seen as a forward energy cascade process. The local energy budget in space under the effective damping of stable modes then yields the spectrum of such a damped turbulence.
iii.1 Effective damping caused by linearly stable modes
We consider the effective damping under the weak nonlinearity assumption, that is, the nonlinear interaction does not change the outer region solution much. Hence we can still use the linear theory to consider the mode structure, and the effective damping rate can be estimated by the negative linear growth rate.
The justification of using the linear growth rate to estimate effective damping is as follows. The effective island width for a given Fourier component of the magnetic perturbation has the following dependence on mode numbers and the magnetic perturbation strength WhiteRMP ; Di2015 :
(9) 
where is the perturbed oblique flux
(10) 
with being the oblique direction defined by given and . Also, is the component of corresponding magnetic perturbation and is the second order derivative of background oblique flux taken at the resonant surface. Further, is the magnetic shear length and is the magnetic shear. Meanwhile, in the inviscid limit, the tearing layer width scales as Rutherford1973 ; Coppi1966 ; Glasser1975 :
(11) 
And in the viscous regime we have Finn2005 :
(12) 
with the magnetic Prandtl number . The following two factors justify the use of linear stability analysis. First, the perturbation amplitudes of high modes are orders of magnitude smaller than that of low modes, thus the effective width of a high island will also be much smaller than that of a low island. Second, the effective island width will shrink faster than the tearing layer width for increasing , as the power dependence on for the former is greater than that of the latter. Simple estimation using the turbulence spectrum obtained later in Section IV indicates that, in our case of weak turbulence, the island width will be smaller than the tearing layer width when for our weak turbulence case. Further, the contribution from the hyperresistivity is also ignored since it is proportional to the driven mode width to the fourth, making their contribution less important for very small scale modes Craddock1991 .
The ideal linear eigenequation for slab geometry can be simply rewritten as Furth1963 ; BiskampBook ; Baalrud2012 :
(13) 
Here, , and is the wave number perpendicular to the oblique direction . It should be noted that we have due to as a result of the localized small scale mode structure. For straight tearing modes with , remains finite even at the resonant surface where . If , then the eigenstructure has the following form near resonant surface :
(14) 
Hence for high modes which are linearly stable, we have:
(15) 
The minus and plus signs here denote the left and the right side of the resonant surface. For oblique modes, there is a logarithmic singularity in the derivative of the ideal solution due to being singular near the resonant surface Furth1963 ; Di2015 . However, the contribution of this logarithmic singularity to is even in parity near resonance, and thus does not contribute to . Hence the of high oblique modes should have the same form as that of straight modes as shown in Eq. (15). Numerical solution of Eq. (13) confirms this statement Baalrud2012 .
The linear growth rate for oblique tearing modes in the inviscid limit is given by BiskampBook ; Baalrud2012 :
(16) 
while in the viscous regime we have Finn2005
(17) 
Here, is the gradient taken at resonance . These dispersion relations have the mode structure and parity corresponding to the driving modes.
As has been mentioned in Section II, the background magnetic field is constantly evolving throughout the course of turbulence, hence we need to track the evolution of numerically as the turbulence evolves. We define the following characteristic length scale of variation
(18) 
so that, in the inviscid limit, we find the following dependence for :
(19) 
while in the viscous regime, we have
(20) 
Thus the effective damping in space can be written as:
(21) 
with in the inviscid limit and in the viscous limit. The Lundquist number is defined as , with and , while is the Alfvén speed considering the guide field. Further, is the effective damping coefficient with dimension of . Combining Eq. (16) or Eq. (17) with Eq. (21), we then have
(22) 
in the inviscid limit, and
(23) 
in the viscous regime. This damping rate has a lower power dependence on than the classical dissipation, making the distinction between the inertial range and dissipation range hard to define. Thus, the present physical situation, in which damping appears to be important at all scales, does not permit a strict delineation of an inertial range in tearing turbulence.
iii.2 Scale independent anisotropy in sheared background field
The scale dependence of the ratio between the parallel and the perpendicular length scales of eddies are of great interest since they will directly affect the nonlinear turnover rate, thus further influencing the forward energy cascade rate of turbulence. The nonlinear turnover rate for MHD turbulence can be modeled as DiamondBook :
(24) 
Here, represent kinetic perturbation at scale.
For weak turbulence generated by oppositely propagating Alfvén waves with straight background field lines, the threewave interaction will preserve the space structure of the beating waves, thus preventing any energy cascade along the direction parallel to the background magnetic field Ng96ApJ ; Ng96POP . A simple way to see this is by considering the resonant condition of wave number and frequency for threewave interaction Sridhar2010 . We have:
(25) 
Here, and represent forward propagating and backward propagating Alfvén waves respectively. The oppositely propagating waves indicate that either or must be zero for both the wave number resonance and frequency resonance to be satisfied. Hence there is no transfer of energy along . And the nonlinear turnover rate scales as as a result.
For a spectrum of modes in a sheared guide field, the parallel length scale where is the dispersion in parallel mode number, and perpendicular length scale where is the perpendicular wave number. With truly weak turbulence in such a guide field where the average field can be seen as unperturbed, we will find a similar independence between and in the turbulence spectrum, that is, we have , although, the physical mechanism are somewhat different as described above. However, it can be seen that the inclusion of finite large scale perturbation will introduce an additional relationship between and in the spectrum, so long as we have where is the random large scale perturbation, is the shear length of background plasma and is the guide field along the ignorable direction. It is important to note that here is the perpendicular wave number of the small scale modes rather than the large scale perturbation, thus the aforementioned criterion should not be confused with the Kubo number of the large scale perturbation, which is defined as the ratio between the nonlinear and linear terms .
We assume the perturbation has the following form: , where and . For the unperturbed background field, we have:
(26) 
Again, , , and . It is worthwhile to emphasize that, and here do not contain the contribution of large scale perturbation. The length represents the distance to the resonant surface for a given . Meanwhile, the dispersion in parallel wave number is
(27) 
Here, represents averaging quantity over for a given , as well as averaging it across the plane for a given . This average over plane is needed because the smallscale mode structures are very localized, and we are looking at the spectrum for a specific x.
It will be shown later on in Section IV.1 that the characteristic length scale of turbulence strength envelope is much larger than in cases we are interested in, thus can be assumed to be independent of for a given . Then Eq. (14) suggests that:
(28) 
Note that, we are integrating over rather than because so long as . For the denominator, we have:
(29) 
Thus we obtain:
(30) 
Realising that the localized mode structure also implies all the small scale modes which can be “seen” from have similar , we can approximately write:
(31)  
(32) 
Therefore we can write down
(33) 
Substituting the above relationship into Eq. (30), the parallel length scale for small scale perturbations is found to be independent of
(34) 
This is similar to the weak turbulence limit discussed previously by Ref. [Ng96ApJ, ] and Ref. [Ng96POP, ], although the underlying physics is quite different.
Now, let’s consider the effect of a large scale perturbation on the small scale anisotropy. We consider the summation of several large scale modes as a random magnetic perturbation strong enough to twist the field “seen” by the small scale modes. Let be the perturbation component in plane. Thus the parallel wave number for each mode is now:
(35) 
Recalling Eq. (32), for a given , we have:
(36) 
For simplicity, let’s define the following parameters:
(37) 
An important feature of these two parameters is that the contribution from the large scale perturbation vanishes upon taking the plane average since vanishes under such spatial average, although does not.
Carrying out the same method as shown in Eq. 28, we again have:
(38) 
as well as
(39)  
Thus, so long as , we have
(40) 
resulting in a scale independent anisotropy . Here, we would like to emphasize that can still be small for the condition to be valid due to the largeness of , with being the wave number of small scale modes.
iii.3 Local energy budget in space
The impact of nonnegligible dissipation on the structure of the spectrum has been studied by considering the local energy budget in space Terry2009 ; Terry2012 . We will follow this methodology here in our analysis, albeit in the context of turbulence with scale independent anisotropy instead of the critically balanced one, as discussed in Section III.2.
Under the local interaction assumption, the local energy budget in space naturally arises from considerations of the effective damping and classical resistive diffusion Terry2009 ; Terry2012 :
(41) 
where , and is the energy forward transfer rate at scale and is the “energy density” in space. Here, we consider a priori the equipartition of magnetic and kinetic energy for medium to high ; we will check the validity of this assumption a posteriori.
The energy budget Eq. (41) can be solved to yield the energy spectrum if can be represented as a function of . The traditional scaling for MHD turbulence without scale dependent alignment indicates that DiamondBook :
(42) 
Here, represents the kinetic perturbation at scale, and again is the Alfvén speed considering the guide field. Due to the equipartition of kinetic and magnetic energy, also represents the forward cascade of magnetic energy as with being the magnetic perturbation at scale.
Using the scale independent anisotropy discussed before, the forward transfer rate is now:
(43) 
Here, is a constant characterizing the scale independent anisotropy, the value of which will be extracted from numerical observation of established tearing turbulence. We follow the closure technique used in Ref. [Terry2009, ]  Ref. [TLBook, ], and write the forward energy transfer rate
(44)  
Here, we have used the closure . This closure technique effectively build the inertial power law behavior into the energy spectrum as an asymptote in low damping limit, therefore, as can be seen later in this section, the spectrum approaches a simple power law when the effective damping vanishes.
Substituting Eq. (44) into Eq. (41), we obtain a linear first order ODE for the energy spectrum:
(45) 
The tearing turbulence spectrum where linear stabilities act as effective damping is then:
(46)  
with and in the inviscid and viscous regime respectively, while effective damping coefficient is given by Eq. (22) and Eq. (23). From Eq. (46), it can be seen that the primary impact of effective damping is an exponential multiplier on the original power law. In the limit of small effective damping, the spectrum recovers the simple power law behavior predicted by the assumption of an inertial range. Also, the power law behavior in the no damping limit tends to be due to the scale independent anisotropy .
Iv Simulation result
In this section, the analytical result from Section III will be tested against resistive MHD simulation using the same set of equations (1)  (4) described in Section II. Especially, we are concerned with the structure of the energy spectrum and the dependence of on . We will look into both the strong guide field case where and , and the comparable guide field case where , and the largescale perturbed field is only one order of magnitude smaller than the guide field . Especially in the latter case, the Kubo number for the large scale perturbation is , corresponding to the regime where turbulent shearing is comparable with parallel propagation. Here, the Kubo number is equivalent to the used by Goldreich and Sridhar in Ref. [GS1997, ]. Numerical observation of the magnetic shear length indicates we have when for both of the above cases. The numerical algorithm will follow that of Ref. [Huang16APJ, ] and Ref. [Guzdar1993, ]. Fivepoint finite difference scheme is used to calculate derivatives, and trapezoidal leapfrog is used for time stepping scheme. The resistivity is set to be , and the magnetic Prandtl number is , and thus we are in the viscous regime discussed above.
iv.1 Strong guide field case


(a)  (b) 


(c)  (d) 
In this subsection, we will compare the tearing turbulence in a strong guide field with our previous theoretical model. Let . , , and .
At the beginning of the simulation, a small initial perturbation with , and is seeded. The resonant surfaces corresponding to those modes lie in the central region of the system , as can be seen from Fig. 1. This is also the region where the turbulent perturbation is strongest later when the turbulence is established. The quadratic form of magnetic and kinetic perturbation is averaged across the  plane, providing us the summation of perturbation energy over the whole spectrum and :
(47) 
(48) 
These perturbation energies as a function of are plotted in Fig. 2 for different times. Initially, the dynamic is dominated by the few largescale unstable modes, and the envelope of their mode structure determines the perturbation energy profile, as can be seen from Fig. 2 (a) and Fig. 2 (b). Later, those initial islands overlap with each other and generate a whole bunch of smallscale modes, and the perturbation energy profile becomes smooth in the core region as can be seen from Fig. 2 (c) and Fig. 2 (d). By the time tearing turbulence has been established, both the magnetic and kinetic energy perturbation are confined within the region , and their profile are almost flattened within the central region. The kinetic perturbation seems to be much smaller than the magnetic one, raising doubt regarding our energy equipartition assumption. However, as will be seen later in Section IV.1, this is because equipartition is established not at the scale of the largescale instabilities driving the turbulence but at the small scales. Meanwhile, the spectrum for high modes actually agrees rather well with the equipartition assumption.
The alignment of the turbulent eddies to the local mean field is also of interest. That is, we wish to know whether or not the eddies have elongated structure along the mean field direction, as would be expected from highly magnetized MHD turbulence. Here, we look at the local property of magnetic perturbation for a given position, and perform Fourier decomposition along and direction for all components of magnetic field. We take the zeroth order harmonic as the local mean field for the given position, while all the other harmonics corresponds to modes with various scales. The alignment of those modes to the direction of local mean field line can be represented by looking at . We once again write down:
(49) 
Such alignment of smallscale tearing modes can then be checked by looking at the distribution of the 2D perturbed energy spectrum and in space. The result for is shown in Fig. 3, where the logarithm of perturbed energy is plotted as a function of mode number and . The black dashed line represents the contour of , with the one originated from point corresponds to . It can be seen that the tearing spectrum strongly aligns with the local mean field, indicating a strongly anisotropic structure. It is noteworthy that, for a magnetically sheared system, this localization in space directly corresponds to the localization of mode structures near their respective resonant surfaces in the configuration space. Due to this localized mode structure, the amplitude of the modes drop away rapidly as we move away from its resonant surface. Thus only modes which are near resonance (corresponds to low ) can be seen from the spectrum shown in Fig. 3, resulting in observed localization in space. This is especially true for high modes. The red dashed lines represent the contours of . The strong alignment behavior of smallscale perturbation in the presence of strong guide field indicates that we have , and it follows that . This confirms our previous assumptions.



(a)  (b) 
With the alignment of eddies known, we now look at how this highly anisotropic turbulence establishes itself. From Fig. 2 (c) and (d), it can be seen that the turbulence strength is rather flat in the central region, this implies that we can use the local spectrum for a given to represent the evolution of tearing turbulence as a whole. Here, we choose to look at the spectrum evolution at . The energy density in space can be obtained by integrating over the red dashed lines in Fig. 3. The spectrum of for several different times is presented in Fig. 4. The logarithm of magnetic energy perturbation is shown as a function of the logarithm of the perpendicular scale . It can be seen that at there are only several largescale unstable modes. Then the interaction of those largescale modes gradually stir up smallscale modes. At a later time, the tearing turbulence reaches a quasisteady state as can be seen in Fig. 4. The structure of spectrum changes very little from to while the longest nonlinear turnover time of the eddies is on the order of .

The comparison between magnetic and kinetic energy spectrum after the turbulence reached the quasisteady state is very important, as we have assumed the equipartition of energy in Section III.3. An example kinetic and magnetic spectrum for quasisteady state tearing turbulence is shown in Fig. 5 for . It can bee seen that for high modes the kinetic and magnetic energy are approximately the same, while at the largest scale there is a departure from this equipartition. This departure does not have a large impact on our theoretical analysis in Section III.3, since we are primarily concerned with smallscale modes which are linearly stable rather than the unstable largescale modes. As a side note, the fact that the magnetic perturbation is one order of magnitude larger than the kinetic perturbation for largest scale modes is also consistent with the observation in Fig. 2, as the total magnetic perturbation energy is also one order of magnitude larger than the total kinetic perturbation energy.

The next important property is the structure function of turbulence, which provides us information regarding the scale dependence of said anisotropy and thus has significant impact on the energy transfer rate and consequently the turbulence spectrum. We follow the procedure detailed in Ref. [Huang16APJ, ] and Ref. [Cho2000, ], and define the following twopoint structure functions:
(50) 
(51) 
Here, is the position of a random point in the configuration space, and is another random vector. Thus and define a random pair of point in configuration space. The bracket here indicates ensemble average over a number of such random pairs. Due to the strong localization of mode structure demonstrated in Fig. 3, we look at a 2D version of the structure function in our study. That is, we take the random pairs within the  plane for a given instead of considering the full 3D space. The parallel and perpendicular component of is defined by the local mean field direction, which is calculated by averaging the magnetic field at two points. We randomly average over pairs of points, and obtain the structure function for both kinetic and magnetic perturbation as functions of and . The contours of this structure function in then reflect the anisotropy of the eddy at different scales. The anisotropy thus obtained is plotted in Fig. 6, with two scalings and plotted as black dashed lines. It can be seen that the anisotropic behavior of simulation result largely agrees with our analytical model, and is mostly scale independent. There is some discrepancy between the length scale of kinetic and magnetic perturbations, which might be the consequence of their different distribution width in space as can be seen in Fig. 3. The ratio between parallel and perpendicular length scale ultimately deviates from the scale independent scaling at the very small scale where classical dissipation kicks in. Lastly, it is observed that is two orders of magnitude larger than , thus we hereby take as a reasonable estimation. This estimation also agrees with our prediction by Eq. (40) since we also have .

With the characteristic structure known, we can finally check our analytical model given by Eq. (46) against the simulation result. The magnetic perturbation spectrum for tearing turbulence is shown in Fig. 7 for and . The simulation result is compared with three analytical models: our damped turbulence model as shown in Eq. (46), a simple power law as a result of traditional inertial range argument, and the spectrum produced by Eq. (46) if only the resistivity is included as damping. From numerical observation, we evaluate the characteristic length scale to be near the central flattened region where the resonant surfaces of the concerned modes lie. Here can be larger than since it only serves as an indication of the local magnetic field gradient. The only free parameter in Eq. (46) is then the energy injection rate , which will be used to fit the simulation result. On the other hand, the power in the simple power law will also be used as a free parameter to fit the numerical result. Those fittings yield and , with fixed parameters and . The fitted curves are shown in Fig. 7, too. It can be seen that our analytical model gives a better agreement with the simulation result than both the simple power law and the resistivity only spectrum. It is noteworthy that although the final decay of the turbulence spectrum is due to the influence of resistivity, as indicated by a similar decay seen in the curve corresponding to resistivity only case, the actual curve deviates from the inertial range behavior long before that as a result of effective damping. This gradual deviation creates the impression that there exists a “inertia range” with a steeper slope as represented by the blue dashed line, while in fact this is merely a misimpression caused by a slow exponential decay.

iv.2 Weaker guide field case
The strong guide field case has been investigated in the previous subsection. Reasonable agreement has been found between the simulation result and our theoretical prediction provided in Section III. The perturbation is found to be two order of magnitude smaller than the guide field, corresponding to a typical weak turbulence scenario. However, we are also interested in the cases where the guide field is weaker, and the Kubo number . Again, here is equivalent to the used in Ref. [GS1997, ]. Note that, in this stronger turbulence case the perturbed field is still smaller than the guide field, although the Kubo number may exceed unity due to the anisotropy.
The initial magnetic fields are now and . To maintain a similar initial safety factor profile with the one shown in Fig. 1, the system size is now , and . We are mainly concerned with the anisotropic behavior and the energy spectrum of the turbulence, and we wish to see whether or not the distinctive features exhibited in our weak turbulence simulation persist for this stronger turbulence case.
We first examine the anisotropy. Again, we look at the contours of structure function for both the kinetic and magnetic perturbation as described in Section IV.1. The position is chosen at , and time , when the turbulence has already reached the quasisteady state. The anisotropy is shown in Fig.8, with the two scaling and plotted as black dashed lines. It can be seen that this stronger turbulence case still follows the scale independent anisotropy behavior described in Section III.2, only deviate from it at very small scale. This stronger turbulence anisotropy exhibts even better independence of the perpendicular scale comparing to that is shown in Fig. 7, possibly due to a stronger large scale perturbation.

We then look at the structure of the turbulence spectrum. The magnetic perturbation spectrum for and is shown in Fig.9. Once again, the simulation result is compared with a simple power law and the spectral form predicted by our model as described by Eq. (46), with fixed parameters and . The fitting result returns . Reasonable agreement is again found between the numerical result and our prediction, with the gradual departure from the original scaling well before entering the resistive dissipation scale.

A noteworthy feature of this stronger turbulence case is that the Kubo numbers for the largest perturbations exceed unity. From Fig. 8, it can be seen that the scale independent anisotropy is approximately . At the same time, numerical observation from Fig. 9 indicates the largest scale perturbation has . Hence, for the large scale perturbations, we have , while for smaller scale perturbation the Kubo number steadily decreases as the perturbation strength decreases. This is different from the critical balance scenario where the Kubo number remains on the order of unity across all scales. This deviation from critical balance is very similar to that discussed by Huang et al. in Ref. [Huang16APJ, ]. Thus we conclude that our analysis can also be applied to the case where the nonlinear mixing is stronger than the linear parallel propagation, such as those reported in plasmoid turbulence, where the critical balance condition is frequently assumed to be true.
V Discussion and conclusion
Instability driven tearing turbulence in sheared magnetic field is studied in this work. The turbulence consists of several largescale unstable modes and a whole spectrum of smallscale linearly stable modes which are stirred up by their largescale counterparts. It is found that the linearly stable modes will act as an effective damping mechanism which has a low power dependence on the perpendicular wave number , which goes as and for inviscid and viscous regime respectively. This low power dependence indicates that this damping mechanism will manifest itself long before we reach the resistive dissipation scale, thus a clean inertial range can not be found, and damping has to be considered in all scales instead. Further, we argue that the tilting of sheared background field by large scale perturbations will impose a scale independent anisotropy for small scale modes. This anisotropic behavior then determines the scale dependence of forward energy cascade rate.
With the knowledge of effective damping rate and energy cascade rate at hand, the structure of this damped turbulence can be obtained by considering local energy budget in space. The key idea is that the difference of energy forward transfer rate between the two ends of any interval in space corresponds to the damping within that interval. The resulting spectrum features a power law multiplied by exponential falloff, as opposed to the pure power law spectrum obtained by using inertial range argument.
The above analytical result is checked against resistive viscous MHD simulation. The turbulence is found to be highly anisotropic and tends to align with the strong local mean field direction. The two point structure function is calculated to demonstrate anisotropic property for different scales, and a scale independent anisotropy is found, confirming our argument. Further, the equipartition between kinetic and magnetic energy is found to be valid for the turbulence in question. The magnetic perturbation spectrum is then fitted with our analytical result, a simple power law and the resistivity only version of our model. The numerical result agrees well with our analytical model. On the other hand, fitting with the simple power law does not produce a physically meaningful result, and the resistivity only result fails in catching the additional departure from the original power law before resistive scale. Thus we conclude that the spectrum of such instability driven turbulence can be well described by our effective damping model.
The behavior of a stronger turbulence, where the Kubo number exceeds unity for certain scales, is also investigated. We find the scale independent anisotropy and the energy spectrum continues to hold for the stronger turbulence case, indicating that our analysis is still applicable even for the scenario where the perpendicular turbulence shearing is stronger than the parallel propagation.
With this knowledge regarding spectrum structure, the next step would be considering the smallscale backreaction to the large scale. This involves a sum of quadratic form of perturbed quantities over the whole spectrum, which requires knowledge regarding the form of spectrum given by our study here. An example is the smallscale spreading of mean field described by hyperresistivity, as is studied in Ref. [Strauss1986, ] and Ref. [Craddock1991, ]. Hence our analysis here provides a solid base for future study along those lines, and the smallscale backreaction to not only the mean field, but also unstable largescale modes can be investigated. This is left to future work.
Acknowledgments
The authors thank P. H. Diamond, X.G. Wang, H.S. Xie and L. Shi for fruitful discussion. This work is partially supported by the National Natural Science Foundation of China under Grant No. 11261140326 and the China Scholarship Council. A. Bhattacharjee and Y.M. Huang acknowledge support from NSF Grants AGS1338944 and AGS1460169. Simulations were performed with supercomputers at the National Energy Research Scientific Computing Center.
References
 (1) H. A. B. Bodin, “The reversed field pinch”. Nucl. Fusion 30 1717 (1990);
 (2) P. H. Diamond, R. D. Hazeltine, Z. G. An, B. A. Carreras and H. R. Hicks, “Theory of anomalous tearing mode growth and the major tokamak disruption”. Phys. Fluids 27 1449 (1984);
 (3) H. R. Strauss, “Hyperresistivity produced by tearing mode turbulence”. Phys. Fluids 29 3668 (1986);
 (4) G. G. Craddock, “Hyperresistivity due to densely packed tearing mode turbulence”. Phys. Fluids B 3 316 (1991);
 (5) W. Daughton, V. Roytershteyn, H. Karimabadi, L. Yin, B. J. Albright, B. Bergen and K. J. Bowers, “Role of electron physics in the development of turbulent magnetic reconnection in collisionless plasmas”. Nature Physics 7 539542 (2011);
 (6) Y.M. Huang, A. Bhattacharjee, “Turbulent magnetohydrodynamic reconnection mediated by the plasmoid instability”. Ap. J. 818 20 (2016);
 (7) A. H. Boozer, “Ohm’s law for mean magnetic field”. J. Plasma Phys. 35 133139 (1986);
 (8) A. Bhattacharjee, E. Hameiri, “Selfconsistent dynamolike activity in turbulent plasma”. Phys. Rev. Lett. 57 206 (1986);
 (9) E. Hameiri, A. Bhattacharjee, “Turbulent magnetic diffusion and magnetic field reversal”. Phys. Fluids 30 1743 (1987);
 (10) U. Frisch, “Turbulence  The legacy of A. N. Kolmogorov”. (Cambridge University Press, Cambridge, 1995), p. 72;
 (11) P. H. Diamond, S.I. Itoh and K. Itoh, “Physical kinetics of turbulent plasmas”. Modern Plasma Physics, Vol. 1 (Cambridge University Press, Cambridge, 2010), p. 52 & p. 350;
 (12) C. S. Ng and A. Bhattacharjee, “Interaction of shearAlfvén wave packet: implication for weak magnetohydrodynamic turbulence in astrophysical plasmas”. Ap. J. 465 845 (1996);
 (13) C. S. Ng and A. Bhattacharjee, “Scaling of anisotropic spectra due to the weak interaction of shearAlfvén wave packets”. Phys. Plasmas 4 605 (1996);
 (14) S. Galtier, S. V. Nazarenko, A. C. Newell and A. Pouquet, “A weak turbulence theory for incompressible magnetohydrodynamics”. J. Plasma Phys. 63 447488 (2000);
 (15) Y. Lithwick and P. Goldreich, “Imbalanced weak magnetohydrodynamic turbulence”. Ap. J. 582 12201240 (2003);
 (16) P. Goldreich and S. Sridhar, “Toward a theory of interstellar turbulence. 2: Strong Alfvénic turbulence”. Ap. J. 458 763 (1995);
 (17) P. Goldreich and S. Sridhar, “Magnetohydrodynamic turbulence revisited”. Ap. J. 485 680 (1997);
 (18) S. Boldyrev, “Spectrum of magnetohyrodynamic turbulence”. Phys. Rev. Lett. 96 115002 (2006);
 (19) D. R. Hatch, P. W. Terry, F. Jenko, F. Merz and W. M. Nevins, “Saturation of gyrokinetic turbulence through damped eigenmodes”. Phys. Rev. Lett. 106 115003 (2011);
 (20) D. R. Hatch, P. W. Terry, F. Jenko, F. Merz, M. J. Pueschel, W. M. Nevins and E. Wang, “Role of subdominant stable modes in plasma microturbulence”. Phys. Plasmas 18 055706 (2011);
 (21) P. W. Terry and V. Tangri, “Magnetohydrodynamic dissipation range spectra for isotropic viscosity and resistivity”. Phys. Plasmas 16 082305 (2009);
 (22) P. W. Terry, A. F. Almagri, G. Fiksel, C. B. Forest, D. R. Hatch, F. Jenko, M. D. Nornberg, S. C. Prager, K. Rahbarnia, Y. Ren and J. S. Sarff, “Dissipation range turbulent cascade in plasmas”. Phys. Plasmas 19 055906 (2012);
 (23) H. Tennekes and J. L. Lumley, “A first course in turbulence”. (MIT press, Cambridge, 1972), p. 268;
 (24) R. B. White, “Resistive reconnection”, Rev. Mod. Phys. 58 209 (1986);
 (25) D. Hu and L. E. Zakharov, “Quasilinear perturbed equilibria of resistive unstable current carrying plasma”, J. Plasma Phys. 81 515810602 (2015);
 (26) P. H. Rutherford, “Nonlinear growth of the tearing mode”, Phys. Fluids 16 1903 (1973);
 (27) B. Coppi, J. M. Greene, J. L. Johnson, “Resistive Instabilities in a diffuse linear pinch”. Nucl. Fusion 6 101 (1966);
 (28) A. H. Glasser, J. M. Greene, J. L. Johnson, “Resistive instabilities in general toroidal plasma configuration”. Phys. Fluids 18 875 (1975);
 (29) J. M. Finn, “Hyperresistivity due to viscous tearing mode turbulence”. Phys. Plasmas 12 092313 (2005);
 (30) H. P. Furth, J. Killeen and M. N. Rosenbluth, “Finiteresistivity instability of a sheet pinch”. Phys. Fluids 6 459 (1963);
 (31) D. Biskamp, “Magnetic reconnection in plasmas”. (Cambridge University Press, Cambridge, 2000), p. 81;
 (32) S. D. Baalrud, A. Bhattacharjee and Y.M. Huang, “Reduced magnetohydrodynamic theory of oblique plasmoid instabilities”. Phys. Plasmas 19 022101 (2012);
 (33) S. Sridhar, “Magnetohydrodynamic turbulence in a strongly magnetized plasma”. Astron. Nachr., 331 93 (2010);
 (34) P. N. Guzdar, J. F. Drake, D. McCarthy, A. B. Hassam and C. S. Liu, “Threedimensional fluid simulations of the nonlinear driftresistive ballooning modes in tokamak edge plasmas”. Phys. Fluids B 5 3712 (1993);
 (35) J. Cho and E. T. Vishniac, “The anisotropy of magnetohydrodynamic Alfvén turbulence”. Ap. J. 539 273 (2000);