We investigate the effects of runaway electron current on the dispersion relation of resistive magnetohydrodynamic modes in tokamaks. We present a new theoretical model to derive the dispersion relation, which is based on the asymptotic analysis of the resistive layer structure of the modes. It is found that in addition to the conventional resistive layer, a new runaway current layer can emerge whose properties depend on the ratio of the Alfvén velocity to the runaway electron convection speed. Due to the contribution from this layer, both the tearing mode and kink mode will have a real frequency in addition to a growth rate. The derived dispersion relation has been compared with numerical results using both a simplified eigenvalue calculation and a M3D-C^{1} linear simulation, and good agreement is found in both cases.

## I. INTRODUCTION

In tokamak discharges, suprathermal electron populations can be generated through the so-called runaway mechanism during both quiescent operation^{1} and disruption events.^{2} These high-energy runaway electrons pose a serious threat to the successful operation of large tokamaks^{3} and remain one of the outstanding challenges for ITER.^{4–6} During a disruption, runaway electrons can be formed through the Dreicer process,^{7–9} the hot-tail mechanism,^{10} or the avalanche mechanism.^{11,12} The avalanche growth rate of runaway electrons will be much larger in ITER disruptions compared to current tokamak devices, given the larger predisruption current and the induced toroidal electric field, as well as longer wall time for electric field diffusion. It is predicted that in a typical disruption event in ITER, a significant fraction of predisruption current will be converted to be carried by runaway electrons with energies of MeVs.^{13,14} In addition, the impurities injected to mitigate other disruption consequences such as wall forces during vertical displacement events (VDEs) can even worsen the runaway electron avalanche issue through partial screening effects.^{6,15–17} Uncontrolled loss of REs can lead to melting of plasma facing components;^{18} thus, it is important to study the confinement of RE beams during disruption events.

The generation of runaway electrons is very sensitive to the local plasma temperature and the parallel electric field.^{19} In addition, the generated runaway electrons can only survive long enough in regions with closed field lines to give rise to avalanche generation. These lead to more peaked current profiles than predisruption,^{20} which alters the magnetohydrodynamic (MHD) stability of the plasma. In recent disruption experiments at DIII-D, significant MHD instabilities have been observed coincident with a large generated runaway current.^{21} The MHD instabilities can lead to deconfinement of both runaway electrons and bulk plasma, with significant oscillations observed in ECE and interferometer signals.^{22} To understand the effect of runaway electrons on MHD instabilities, effort has been taken on both theoretical analysis and numerical simulations.^{23–26} In these works, it is found that the linear growth rate of MHD instabilities is approximately the same as in a plasma without runaway electrons, but the nonlinear saturation level is different.^{23} However, the perturbed current profile is much more peaked near the rational surface.^{25} These studies are conducted using a fluid description of runaway electrons and a current coupling scheme into the MHD equations, which can save computation time compared to electromagnetic kinetic simulation.^{27} Recently, this method has been implemented in the M3D-C^{1} code, which is a fully 3D nonlinear MHD code based on the finite element method.^{28,29} Both *m* > 1 (*m* is the poloidal mode number) tearing modes and resistive kink modes in the presence of runaway electrons have been studied using M3D-C^{1}, and several new interesting phenomena have been found in simulations.^{30}

One of the newly found phenomena observed in our M3D-C^{1} numerical simulation of MHD tearing instabilities in the presence of runaway electrons is that the mode has a real frequency (overstability) driven by the runaways and is propagating in the plasma frame. When real frequencies are driven in resistive MHD instabilities, it is critically important to understand the physics driving them and their effect on the mode, especially when studying experimental phenomena. Overstability can be highly impactful to both the stability of resistive MHD modes and the physics of their interaction with the plasma and resistive wall. The real frequency of the mode has been studied as an effect of plasma pressure gradient,^{31–34} two-fluid effect,^{35–38} polarization drifts,^{39} energetic ion interaction,^{40,41} effects of favorable average curvature,^{42} and coupling with other modes.^{43} It has been used to explain the Maxwell torque and plasma rotation as mode locking happens.^{44,45} With multiple effects simultaneously present, a single frequency is generally found for the mode depending on the combined influences. Under different conditions, each individual driving mechanism can be particularly important to the stability of the mode.

In this paper, we present a new mechanism causing the mode rotation that can explain the real frequency observed, which is attributable to the convection of runaway electrons. This work is based on the theoretical derivations in Ref. 23, where asymptotic analysis has been utilized to study reduced MHD equations^{46} plus a runaway electron convection equation in a Cylindrical plasma. In addition to the classical resistive layer theory, we find that runaway electrons can cause a secondary layer structure to form inside the resistive layer, whose width depends on the value of $vA/c$ (*v _{A}* is the Alfvén velocity and

*c*is the speed of light) and can be much narrower than the resistive layer. In addition to the

*m*> 1 tearing mode, we also use the model to study the $m=1,n=1$ resistive kink mode with runaway electrons based on Ref. 47, whose eigenfunction of perturbed flux behaves very differently than the

*m*> 1 tearing mode.

This paper is organized as follows. In Sec. II, the resistive tearing mode in the presence of runaway electrons with *m* > 1 is studied, and a new dispersion relation is derived, which gives both the growth rate and the real frequency. In Sec. III, a similar method is applied to study the *m* = 1 resistive kink mode. The results obtained in both sections are compared with numerical calculations using both a 1D eigenvalue code and the M3D-C^{1} MHD simulation code, which is presented in Sec. IV. We summarize in Sec. V.

## II. RESISTIVE TEARING MODE WITH *m* **>** 1

We first study a tearing mode with *m* > 1 in a periodic straight cylinder with its axis in the *z* direction. In the framework of reduced MHD, the magnetic field and the plasma flow can be conveniently represented as

where *ψ* is the normalized magnetic flux ($\psi =Az/(BTa)$, with *a* being the minor radius of the cylinder) and $\varphi $ is the normalized electric potential ($\varphi =\Phi /(BTavA)$). The linearized equations of motion of the magnetic fields, the perpendicular velocity fields (represented by the electric field since $E\u22a5=v\u22a5\xd7Bz$), and the convection equation of runaway electrons can be written as^{23}

Here, we ignore the effects of plasma pressure, given that the plasma *β* is close to zero in post-disruptions. *γ* is the linearized growth rate normalized to Alfvén time $\tau A=a/vA$, *j _{RE}* is the perturbed runaway electron current density,

*J*

_{0}is the equilibrium current density, $JRE0$ is the equilibrium runaway electron current density,

*η*is the normalized resistivity ($\eta =\eta \u2225/(\mu 0vAa)$), and

*m*and

*n*are the poloidal and toroidal mode numbers. In addition,

where *R* is the major radius and $q=rBT/(RBP)$ is the safety factor, with $BP=|\u2207\psi |$ being the poloidal field. Note that in Eq. (5), the runaway electrons are assumed to have a convection velocity equal to *c* on the opposite direction of the magnetic field line (so the runaway electron current is positive), plus the $E\xd7B$ drift, given that most runaway electrons are relativistic with $v\u2225\u226bv\u22a5$.

Note that in Eqs. (3)–(5), the gradient and curvature drift of runaway electrons are ignored. For high energy runaway electrons, the curvature drift can be important and comparable to the $E\xd7B$ drift. However, the magnitude of these drift terms depends on runaway electron energy, which is missing in the fluid description of runaway electrons. In addition, the effect of curvature drift can be canceled by vertical fields from external coils to maintain the force balance.^{48} The pressure of runaway electrons is also absent, since we are using a current coupling scheme rather than pressure coupling for runaway electrons. The effects of these terms on MHD stability will be studied in the future.

Given that $\gamma \u226a1,\eta \u226a1$, for the solution far away from the rational surface (outer region), terms proportional to *η* and *γ* can be ignored. The solution $\psi out$ can be used to determine the boundary condition of the inner region. For the solution close to the rational surface, $k\u2225\u21920$ and the resistive term becomes important. We introduce the coordinate variable $x=r\u2212rs$, where *r _{s}* is the minor radius of the rational surface. The equations in the inner layer can be rewritten as

where $x\u226ars$ is assumed, $\xi =ikc\varphi /\gamma $ is the plasma displacement, $kc=nq\u2032/(Rq)$, and^{49}

Note that the resistive term is only important in the layer $|x|<|\delta 1|$, where

The last two terms on the right-side of Eq. (9) come from the right-side of Eq. (4), which are ignored in Ref. 46 in the inner layer calculation. By including the first order contribution of these two terms, we can obtain a corrected value of the resistive tearing mode growth rate, as discussed in Ref. 49.

The perturbed runaway electron current in Eq. (5) can be rewritten as follows in the inner layer:

Assuming $JRE0=J0$ and ignoring the terms proportional to $vA/c$ in Eq. (13), we find that the runaway electron current term can cancel the last two terms in Eq. (9), leaving the equation exactly the same as the approximate equation in Ref. 46. Thus, the growth rate will follow the 3/5 power law with respect to *η* in Ref. 46, as discussed in Ref. 23. For $JRE0\u2260J0$, the last two terms in Eq. (9) will then come only from the thermal electron current $J0\u2212JRE0$, and the mode growth rate will be reduced by a certain factor.

We now consider the effect of finite value of $vA/c$. Given $vA\u226ac$, the term proportional to $vA/c$ in the numerator in Eq. (13) can be ignored. However, the second term in the denominator leads to an imaginary correction to *j* for $x\u21920$. Defining

we find that $|\delta 2|\u226a|\delta 1|$; thus, this correction is only significant in a sublayer inside the inner layer. To satisfy Eq. (8), we can, therefore, assume $\psi \u2033=\u2212jRE$ inside this sublayer. Following the constant-*ψ* approximation, it is found that the two terms on the right hand side of Eq. (8) are subdominant compared to $\psi \u2033$, given $|\delta 2|\u226a\gamma /\eta $. Therefore,

Here, we apply the Sokhotski–Plemelj theorem to calculate the complex integral and ignore the principle value, which is similar to the calculation of the imaginary term in a Landau contour integral. The sign of the integration result will be the same as the sign of the term $\gamma vA/(kcc)$. Note that the real part of *ψ* will also be affected by *j _{RE}* and behave like $\psi =\u2212\psi 0[1+mJ\u2032RE0/(kcrs)x\u2009ln\u2009|x|]$ near

*x*= 0, as discussed in Ref. 23. This part will also contribute to the matching condition through the last two terms on the right hand side of Eq. (9).

To match to the ideal MHD solution in the outer region, the solution needs to satisfy

where $\Delta \u2032=(\psi \u2032(x=0+)\u2212\psi \u2032(x=0\u2212))/\psi 0$, where $\psi \u2032(x=0+)$ and $\psi \u2032(x=0\u2212)$ are taken from the solution in the outer region. Here, we take the “constant-*ψ*” approximation in the layer by assuming $\psi =\psi 0$.

For the region $|\delta 2|\u226a|x|\u226a|\delta 1|$, the correction caused by finite value of $vA/c$ is very small. In addition, as discussed above, the runaway electron current term can cancel the last two terms in Eq. (9), leaving the equation exactly the same as the approximate equation in Ref. 46. We can, thus, assume that the functional form of the solution of *ξ* is little affected by the existence of the sublayer. Note that the amplitude of the solution, including both real and complex parts, will be affected by the sublayer as shown below.

Introduce $z=x/\delta 1,\u2009\chi =\delta 1\xi /\psi 0$, then in the region $|z|\u226b|\delta 2|/|\delta 1|$, and the equation can be simplified as

The solution can be obtained by doing a Fourier-Laplace transformation^{50}

and the integral required for matching condition is

Putting this result into Eqs. (15) and (16), we obtain the result of the integral in the region inside the resistive layer but outside the sublayer,

where the higher order terms are ignored. The last term is calculated by substituting the real part of *ψ* near *x* = 0 into Eq. (9). Note that the last term is small compared to $\Delta \u2032$ and can be ignored here. The matching condition that determines the growth rate is

where we assume $Re(\gamma )>0$. Comparing this result with Ref. 46, we find that the effect of runaway electron current on the resistive tearing mode can be treated by adding an imaginary term to $\Delta \u2032$. This correction term is of the same order of the original integral of $\psi \u2033$ in the resistive layer. This will make *γ* a complex number. The real part is the mode growth rate and the imaginary part means that the tearing mode will have a real frequency and propagate in the plasma frame. The amplitude of this imaginary term is proportional to $JRE0\u2032$ at the rational surface; thus, for the case with $JRE0\u2032=0$ (like the tearing mode in slab geometry with a symmetric current profile or a flat runaway electron current profile across the rational surface), the correction is zero. This imaginary part does not have a strong dependence on the value of $vA/c$, as long as $|\delta 2|\u226a|\delta 1|$. The sign of this term is determined by the sign of $JRE0\u2032$ and the streaming direction of the runaway electrons.

## III. RESISTIVE KINK MODE WITH *m* **=** 1

For a resistive mode with *m* = 1, the equation in the inner region is the same as Eqs. (8) and (9). However, the solution will match that of the internal kink mode in the outer region, which is $\xi =\xi 0$ for $r<rs$ and *ξ* = 0 for $r>rs$ for plasma *β* close to zero. This is a different boundary condition from the *m* > 1 tearing mode. In addition, the perturbed flux *ψ* will also be zero in the outer region, which means that the constant *ψ* approximation cannot be used to calculate the inner region solution.

In previous studies of the *m* = 1 resistive kink mode without runaway electron current, by ignoring the last two terms in Eq. (9), an exact solution is found,^{47,50}

where

The outer region gives the matching condition as

Thus, the growth rate satisfies

By including the runaway electron current term, as in Sec. II, the runaway electron current will cancel the last two terms in Eq. (9), making the growth rate close to the value given in Eq. (26) in both the small and large *η* cases. Like in Sec. II, we find that *j _{RE}* term can cancel the last two terms in Eq. (9), leaving the combined equation of

*ξ*exactly the same as the approximate equation in Ref. 47. The solution of

*ξ*will also be the same as in Eq. (22). The form of

*ψ*, however, will be affected by the presence of the last two terms in Eq. (9). Here, we use the form in Eq. (23) as the leading order solution of

*ψ*together with Eq. (9) to calculate the next order solution of $\psi \u2033$.

Like the derivation in Sec. II, the introduction of runaway electron current can give rise to a sublayer of *ψ* near *x* = 0 with the width $|\delta 2|$ in Eq. (14). This leads to an additional imaginary term in the integral of $\psi \u2033$ across the sublayer. Using the result of $\psi (x=0)$ calculated from Eq. (23), we get

In addition, the second term in Eq. (23) can also contribute to the integral of $\psi \u2033$ in the region inside the resistive layer but outside the sublayer, through the last two terms on the right side of Eq. (9),

In addition, there is also a contribution from the real part of *ψ* near *x* = 0 like the last term in Eq. (20). However, the calculation shows that this term is proportional to $\delta \u030212\u2009ln\u2009\delta \u03021$, which is a higher order effect compared to Eq. (28).

Combining these two additional terms with the integral of $\xi \u2033$, we obtain the new matching condition

Note that the last term in Eq. (29) does not depend on *J _{RE}*, which is simply a correction term due to the last two term in Eq. (9), and will lead to a change in the growth rate rather than real frequency. Given $\eta 1/2/\gamma 1/2\u223c\eta 1/3\u226a1$, the correction term coming from runaway electron current is of higher order compared to the original solution, which is different from the results in Eq. (21). However, for plasma with large resistivity like in a post-disruption scenario discussed in Sec. IV, the correction can still be important. In addition, the runaway electron current will also contribute to the real part of the dispersion relation, which will change the value of the growth rate and make it deviate from the 1/3 power law for large

*η*cases.

## IV. COMPARE WITH NUMERICAL RESULT AND M3D-C^{1} SIMULATION RESULT

In this section, we compare the MHD instability growth rates and real frequencies obtained in Secs. II and III with both an eigenvalue code and M3D-C^{1} simulation result. The eigenvalue code solves Eqs. (3)–(5) using the finite difference method. For the M3D-C^{1} code, we use the two-field linear version, which only solves for $\varphi $ for the velocity field and *ψ* for the magnetic field, in order to benchmark with reduced MHD results. The MHD equations and runaway electron equations are solved on a 2D circular plane, and the toroidal derivative terms are represented with a spectral method. To resolve the runaway current sublayer structure, fine mesh packing is applied near the rational surface where the mode is excited. The results of linear simulations with full MHD equations are discussed in Ref. 30.

Before proceeding with numerical results, it is useful to make some simple estimates on the value of normalized resistivity *η* on which the growth rates and the real frequencies will depend. For a runaway electron experiment in quiescent operation in DIII-D, where ions are all deuterium and *B* = 2T, electron density $ne=1021$ m^{−3}, electron temperature *T _{e}* = 1 keV, and

*a*= 0.67 m, the normalized resistivity

*η*is only about $2\xd710\u22129$ and $c/vA=30$. However, for a post-disruption scenario in ITER, with

*B*= 5T, $ne=1021$ m

^{−3},

*T*= 5 eV,

_{e}*a*= 2 m, $\eta \u224810\u22125$, and $c/vA\u2248123$. In addition, if the plasma is composed of argon instead of deuterium due to impurity injection, the values of

*η*and $c/vA$ can increase to $4\xd710\u22125$ and 550.

We first compare the results of *m* > 1 resistive tearing modes from these different methods. The MHD equilibrium in a cylinder is set up with $q=1.15(1+r2/0.6561)$. For the case with runaway electron current, all equilibrium current is assumed to be carried by runaway electrons ($J0=JRE0$). In Fig. 1, we show the values of *γ* for an $m=2,\u2009\u2009n=1$ mode, for cases both with runaway electron current and without runaway electron current. We can see that in the case without runaway electron current, the growth rate from both the eigenvalue code and the M3D-C^{1} simulation deviate from the 3/5 power law^{46} for large *η* but agree with the result in Ref. 49. After including the runaway electron term in the equation, the growth rates from both numerical codes agree well with Eq. (21) except for very large *η*. This is consistent with the numerical calculation in Ref. 23. The real frequencies of the modes from both the numerical codes also agree well with Eq. (21) and are of the same order as the growth rate for all values of *η*. Since in all cases $|\delta 2|\u226a|\delta 1|$, it is found that both the growth rates and real frequencies are insensitive to the value of $c/vA$, consistent with the results reported in Ref. 30. The deviation at very large *η* is due to the fact that the large value of the tearing layer width *δ* can become comparable to *a* in these cases, which is the nonasymptotic regime. Thus, the assumptions we took in the derivations in Sec. II, such as constant-*ψ* approximation and $x/rs\u226a1$, are no longer valid.

Figure 2 shows the real and imaginary parts of *ψ* and $\psi \u2032$ for the (2,1) modes near the rational surface, from the eigenvalue code, for $\eta =10\u22125$ and $c/vA=40$. We can see clearly the existence of a sublayer structure. Due to the existence of large *j _{RE}* near

*x*= 0, the imaginary part of $\psi \u2032$ will behave like a step function. The change of $\psi \u2032$ inside the sublayer will then be compensated at $|\delta 2|<|x|<|\delta 1|$ to satisfy the boundary condition at the outer bounds. The real part of $\psi \u2032$ is also affected by

*j*and behaves like $\u223c\u2009ln\u2009|x|$ near

_{RE}*x*= 0, as discussed in Sec. II.

We next compare *γ* of (1,1) resistive kink modes. The safety factor profile is now set to be $q=0.9(1+x2/2)$ to incorporate a *q* = 1 surface. The results for the growth rates and real frequencies are shown in Fig. 3. The analytical results are calculated from Eq. (29) using an iterative method. Similar to those for (2,1) modes, the growth rates deviate from the 1/3 power law for large *η* without runaway electron current, but with runaway electron current, the results follow the 1/3 power law except for very large *η*. This is due to the correction to the real part in Eq. (29). These results are consistent with those in Ref. 25. The real frequencies are much smaller than the growth rate for small values of *η*, but for large *η*, the real frequency can be even larger than the growth rate. The results for very large *η* deviate from Eq. (26) because the very large layer width breaks down the assumption used in the calculations in Sec. III.

## V. CONCLUSION

In this paper, we establish a new theoretical model to understand the linear properties of both resistive tearing and resistive kink modes in the presence of runaway electron current. The model is based on the asymptotic analysis method for the reduced MHD equations in cylindrical geometry as presented in previous studies of runaway electron current's effect.^{23} Extending from the analysis in Ref. 23, we studied the sublayer structure of the mode formed by the convection equation of runaway electron current and its effect on the mode dispersion relation. It is found that the sublayer can have a strongly peaked perturbed current, whose width can be much smaller than the resistive layer and depends on the value of $vA/c$. On the other hand, this runaway electron sublayer can give rise to an imaginary term in the dispersion relation whose amplitude does not depend on the value of $vA/c$, which is similar to the pole contribution in a Landau contour. The growth rates and real frequencies of modes calculated using this dispersion relation are compared with the results obtained using an eigenvalue calculation and M3D-C^{1} simulations.

The real frequency of the mode found in the paper is proportional to the gradient of runaway electron current $JRE\u2032$ at the rational surface. This current gradient will break the symmetry of the profile across the rational surface, making an odd part in the perturbed flux function. As shown in Ref. 49 and Eq. (29), this gradient can also change the mode growth rate. In addition, as discussed in Ref. 51, the break of symmetry can affect the saturation level of the tearing mode. It is believed that such an effect can be important for disruptive plasma since runaway electron current can have a more peaked profile than thermal plasma.^{20}

With the simplified model presented here, including a mode resonant at a single surface and a perfectly conducting wall, the effect of the runaways on the layer solution, frequency, and growth rate are clearly important. Given multiple rational surfaces and/or a surrounding resistive wall and bulk plasma flow, the contributions from the layer response at each surface, and the currents driven in the wall, will again combine to select a single frequency of the mode. The mode frequency determines the rate of propagation in the plasma frame and, via a Doppler shift, relative to the resistive wall. The stability and structure of the mode will then depend on the conditions at all the surfaces, the plasma flow, and the properties of the resistive wall. Given a significant runaway current and gradient across surfaces, ultimately the stability of the runaway electron confining plasma depends on this real frequency, and thus so does our ability to keep the runaways confined. Clearly, including this physics into MHD stability calculations of experiments with runaway electron generation is critically important.

The result of this work can help future numerical and theoretical simulation research on runaway electrons. The very thin layer of perturbed current inside the resistive layer poses a significant challenge to simulation of MHD modes with runaway electrons, especially for *m* > 1 tearing modes. To resolve the mode structure and get the correct value of mode real frequency, either a fine mesh packing is required near the rational surface or one has to use a reduced unphysical value of $c/vA$ like in Ref. 24. In addition, the width of this sublayer can be comparable to the electron skin depth ($de=c/\omega pe$, where *ω _{pe}* is the plasma oscillation frequency). This means that inside the sublayer, the electron inertial term ($de2dj/dt)$ can become important and should be included in Ohm's law. By including this term, we find that for the ITER disruption case ($\eta \u223c10\u22125$), the inertial term is still a subdominant term compared to the resistive term (

*ηj*). However, for the case with smaller

*η*like the quiescent runaway electron experiments, the inertial term can be important. The inertial term has been shown to play a role in the fast Hamiltonian reconnection

^{52–54}and gives rise to a singular current spike. The nonlinear evolution of the modes can also be strongly affected by the existence of this sublayer and current spike, the study of which we leave for future work.

## ACKNOWLEDGMENTS

Chang Liu would like to thank Per Helander, Guo-Yong Fu, Hank R. Strauss, and Akinobu Matsuyama for fruitful discussion. This work was supported by the Simulation Center of electrons (SCREAM) SciDAC center by Office of Fusion Energy Science and Office of Advanced Scientific Computing of U.S. Department of Energy, under Contract Nos. DE-SC0016268 and DE-AC02–09CH11466. This research used the high-performance computing cluster at the Princeton Plasma Physics Laboratory and the Eddy cluster at Princeton University.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

*C*

^{1},”