Beams of energetic runaway electrons are generated during disruptions in tokamaks, and fluid models are used to study their effects on macroscale dynamics. Linear computations of a massless, runaway electron beam coupled to MHD plasma show that resistive hose instabilities grow faster than tearing modes at large resistivity. Eigenvalue results with reduced models of the resistive hose instability are compared with results from the full MHD and beam system, showing that the resistive hose decouples from any plasma response. An estimate of plasma temperature at which growth of the resistive hose dominates tearing for post-disruption DIII-D plasma parameters is in a physically relevant regime.

The tokamak is the leading candidate for a magnetic confinement fusion reactor. It confines plasma in a toroidal configuration with an externally generated toroidal magnetic field and a poloidal magnetic field supported by toroidal electrical current that is driven through the plasma. Tokamaks are prone to disruptions—sudden dynamic events in which confinement of the plasma is lost and the plasma current is terminated. A typical disruption begins with a rapid loss of plasma thermal energy, the thermal quench (TQ), followed by a slower current quench (CQ). The decrease in temperature during the TQ increases the electrical resistivity of the plasma, increasing the electric field along magnetic field-lines and accelerating a small population of electrons to high energies. Since the effective collision frequency decreases at higher energies, the resulting beam-like population of supra-thermal electrons is largely unimpeded by collisional drag.^{1} These runaway electrons (REs) can carry a substantial fraction of the initial plasma current and have caused considerable damage to plasma-facing components. Understanding the dynamics of REs during disruptions is critical for the development of tokamak fusion reactors. For the ITER experiment, which will carry a plasma current of 15 MA, mitigating RE-induced damage will be essential.^{2} To that end, the development and characterization of computational tools to simulate RE dynamics in tokamaks is required.

*γ*is the relativistic factor associated with the runaway parallel velocity: $ \gamma r = 1 / 1 \u2212 ( c r / c ) 2$. Since the RE current is essentially resistance free,

_{r}^{1}the typical resistive MHD Ohm's law is modified,

Because electrical resistance is modified, a natural question is the effect of the inclusion of REs on resistive MHD instabilities. Prior work with this model has explored the effect of a fluid runaway current on the linear behavior of the tearing mode in a sheet pinch,^{3} and tearing and resistive kink modes in a large aspect ratio cylinder.^{4–6} Helander *et al.* additionally consider nonlinear saturation of the tearing mode.^{7} Li *et al.* perform linear tearing mode calculations in toroidal geometry and include effects related to non-zero plasma pressure.^{8}

The stability of high-energy, relativistic particle beams propagating through resistive, neutralizing background plasma has also been studied. The earliest analytical description of an instability in such a system is due to Rosenbluth.^{9} Rosenbluth identifies a kink-type instability of a self-pinched beam of relativistic particles that couples the motion of the beam particles with the resistive response of magnetic field. Weinberg generalizes the work of Rosenbluth first to a spatially modulated beam,^{10} where the instability is dubbed a “hose” mode, and later, to a more general dispersion relation for the self-pinched equilibrium.^{11} In each case, the dynamics of the beam particles are treated via a relativistic kinetic description, and the background plasma is assumed to interact only by providing an Ohmic current density, $ \sigma E$. Rosenbluth uses a “rigid-beam” approximation for low-frequency *m* = 1 modes, where the entire beam density at any *z*-value is displaced from the magnetic axis without distorting. Assuming also that the resistive skin-depth is large compared to the beam radius, $ \mu 0 \sigma \omega r 2 \u226a 1$, allows an analytic derivation of a dispersion relation for arbitrary initial beam profiles. Weinberg's treatment extends to the more general case, which recovers the rigid-beam limit for low-frequency perturbations. Closed-form dispersion relations are given for the case of a flat beam current profile. The hose mode has long been recognized in the accelerator community as an impediment to the propagation of beams, and controlling its growth with properly tailored beam profiles has been a subject of study.^{12}

In the resistive hose case, it is assumed that the plasma column surrounding the beam is generated by ionization of a dense gas by the beam particles. The beam is self-pinched, meaning that the self-induced magnetic field is azimuthal, and the equilibrium force balance condition requires that this field is sufficiently strong to provide the centripetal acceleration of the beam particle orbits around the axis. Some relatively recent work includes the effect of an axial magnetic field, which modifies the equilibrium orbits of the beam particles.^{13} In contrast, in the case of the RE beam generated in a tokamak, the background plasma carries current density before the RE population is generated. There is an existing MHD equilibrium with both axial (toroidal) and azimuthal (poloidal) magnetic field, and after its generation, the RE beam carries some substantial fraction of the equilibrium current density.

Stability of runaway electron beams in tokamaks has been extensively analyzed with models of the background plasma that neglect resistivity.^{14–16} Although the model fidelity varies, a consistent result is that the beam modifies the marginal stability of the system to kink modes, slightly, and introduces a real frequency when viewed in the lab frame. The neglect of a mechanism of current decay in the plasma means that these analyses do not describe resistive hose modes.

In this Letter, we draw the connection between the earlier results on resistive hose instabilities of particle beams and resistive instabilities in the context of RE modeling in tokamaks. It has previously been noted that resistive hose instabilities are possible in a tokamak, but they have been regarded as unimportant since the resistivity of tokamak plasmas is usually very low.^{17} However, in the context of a post-TQ runaway beam, the high resistivity of the cold background plasma suggests it may be worthwhile to revisit the relevance of hose modes.

The dynamics of the RE beam in this work are treated under the assumption that the electrons have negligible mass and are guided by the tokamak magnetic fields in contrast to the particle-beam analysis, where inertia is important for the particle orbits. No analysis of the resistive hose in the literature has, heretofore, explored the massless limit. To check its relevance, we estimate that the resistive hose mode that appears in this work would grow faster than the tearing mode at post-TQ plasma temperatures of $ T e \u2248 1.5$ eV or less in DIII-D.^{18} We do not suspect that the hose instability alone is responsible for the final loss of RE current, but it may couple with other modes, nonlinearly.

The model consists of a background plasma whose dynamics are governed by single-fluid, resistive MHD. In the present discussion, the effects of pressure are presumed to be negligible, as appropriate for post-TQ conditions. The MHD equations are augmented with a second fluid species that represents the collective low-frequency behavior of REs. The model for REs is an electron beam with a large constant speed parallel to the magnetic field. In the present work, the $ E \xd7 B$ and curvature drifts are neglected. The RE beam dynamics are coupled to the resistive MHD equations via the modification to the Ohm's law noted earlier. Quasi-neutrality between the species is assumed so that there are no space-charge effects from the beam.

**b**, the bulk plasma velocity,

**v**, and the runaway density,

*n*. The variables with capital letters and the mass density

_{r}*ρ*represent the equilibrium fields, and the lowercase variables are the perturbations.

The model in Eqs. (1)–(9) is intended to capture the effect of runaway electron current on resistive MHD instabilities. An early analytic application of this model to tearing modes in slab geometry found that in the limit of zero inertia, stability of the tearing mode is unaffected.^{3} More recent results using reduced MHD in cylindrical geometry also suggest that the runaway current does not impact the stability of the tearing mode, but nonlinear saturation is affected.^{7} The analytic and numerical work in Ref. 4 finds an additional correction to the linear dispersion relation that results in overstability of the tearing mode.

We have implemented Eqs. (1)–(9) in the NIMROD extended MHD code.^{19} The linear equations are Fourier analyzed in the axial direction, and a high-order spectral element representation is used in the poloidal $ ( r \u2212 \theta )$ plane. For a single axial Fourier harmonic, *e ^{ikz}*, the initial value problem is solved by numerically integrating the linear equations in time from an arbitrary initial condition. As time $ t \u2192 \u221e$, the solution approaches the most unstable mode from which the growth rate and frequency are determined.

For cylindrical problems, the equations can be reduced to a set of coupled ordinary differential equations in the radial coordinate by also assuming a Fourier representation in the azimuthal angle, $ e i m \theta $, and a complex exponential dependence in time $ e \gamma t e \u2212 i \omega t$. To verify NIMROD results, a 1D eigenvalue code is used to numerically solve the resulting ordinary differential equations in *r* to determine *ω* and *γ* given *k* and *m*. The eigenvalue solver also uses spectral elements, albeit for the single spatial dimension *r*.^{20}

^{21}The safety factor profile is $ q ( r ) = q 0 ( 1 + q 2 ( r / a ) 2 )$. In the cases presented, $ q 0 = 1.15$ and $ q 2 = 1.54$, so that $ 1 < q ( r ) < 3$ and

*q*increases monotonically in

*r*. Thus, the only

*n*= 1 MHD instabilities expected in the absence of runaway effects are tearing modes with $ m \u2265 2$. To lowest order in

*a*/

*R*, the shape of the parallel current profile is given by

The cylinder is considered periodic in the axial coordinate *z*. It has a unit minor radius, *a* = 1; there is a perfectly conducting wall placed at *r* = *a*; and the RE beam radius is equal to the plasma radius. The axial field strength at *r* = 0 is $ B 0 = 1$. The density *ρ* is chosen so that the Alfvén velocity defined by the axial field at *r* = 0, $ V A \u2261 B 0 / \mu 0 \rho $ is unity when $ B 0 = 1$. A particular value of $ c r = 20 V A$ is chosen as a representative example for $ c r \u226b V A$. Lengths are normalized by *a*, resistivity is normalized to $ \mu 0 a V A$, and the frequencies and growth rates are normalized to $ V A / a$. In the results that follow, the normalized axial wavenumber is $ k = \u2212 0.1$.

Figure 1 plots the growth rate and frequency of the fastest growing mode from NIMROD calculations of Eqs. (1)–(9) as a function of varying resistivity values. Also plotted are the growth rates of the fastest growing mode in the absence of the RE beam, growth rates and frequencies computed with the M3D-C1 code,^{4} and growth rates and frequencies computed with our 1D eigenvalue code with *m* = 2. For values of $ \eta \u2272 10 \u2212 4$, the fastest growing mode is the *m* = 2 tearing mode, which displays the expected localization around the *q* = 2 surface. For $ \eta > 10 \u2212 4$, the fastest growing mode computed in NIMROD is a hose-like mode with *m* = 1. The appearance of this hose mode accounts for the distinct growth rate and frequency at $ \eta = 10 \u2212 3$ reported by NIMROD. Figure 2 plots contours of the radial and azimuthal components of the solution for $ \eta = 10 \u2212 2$.

The transition from the tearing mode to the hose mode is explicated by calculating each *m* separately with the eigenvalue code. Figure 3 plots the frequencies and growth rates for *m* = 1 and *m* = 2 as a function of resistivity. The shaded region denotes the range of *η* values where the hose mode grows faster than the tearing mode. It is clear from the eigenvalue code results that the hose mode scales linearly with resistivity for small *η*, whereas the tearing mode scales as $ \eta 3 / 5$. The difference in scaling suggests that for equilibria that are unstable to both modes, there is a value of *η* at which the fastest growing mode transitions from tearing to the hose. The exact details of the transition depend on the equilibrium current profile, and quantifying this will require an analytic expression for the hose dispersion relation in this model. Deriving such an expression remains a task for future work.

**v**. This is a good approximation when the background plasma responds on a much longer timescale than the hose mode dynamics. In making this simplification, we also redefine the RE density variables as

*V*is not a relevant parameter.

_{A}*ψ*and

*λ*,

*ψ*. As shown in Fig. 4, numerical solutions of these equations confirm that they are sufficient to reproduce the features of the hose instability observed in the solutions of Eqs. (1)–(9).

The agreement between these three calculations confirms that the flow of the background plasma is unimportant and that the polarization of the magnetic field is predominantly transverse. Moreover, the form of Eq. (16) is equivalent to the field equation derived by Rosenbluth and Weinberg to describe the resistive hose mode.^{9,11}

The hose mode in these calculations does not extend beyond a given radius. Figure 5 shows radial profiles of the hose mode. It is clear that there is some radius $ r \u0303$ outside which the profiles of the solution for *λ* and *ψ* drop to zero. The value of $ r \u0303$ appears to be determined by the solution to $ \omega + c r F ( r \u0303 ) = 0$, which is plotted in Fig. 6. Although there is no *m* = 1 rational surface in the tearing sense of *F*(*r*) = 0, the equation for *λ* becomes singular at the radial location where $ \omega + c r F ( r ) = 0$. Unlike tearing where *ω* = 0, the location of the singularity is not known *a priori*.

The shape of the equilibrium beam current profile also affects the behavior of this mode. Varying the *q*_{2} parameter of Eq. (10) in the eigenvalue calculations shows that the growth rate and frequency of the mode increases as the current profile becomes more peaked despite a decrease in the total beam current. Changes in the shape of the Λ profile also affects the *c _{rF}* term, and the eigenfunction in higher

*q*

_{2}case is more localized near

*r*= 0. Further work to elucidate the precise influence of the equilibrium profile on the instability characteristics is ongoing.

We have shown that resistive hose-type modes associated with a RE beam interacting with a thermal background plasma grow faster than the typical resistive MHD tearing mode at large resistivity. The scaling of the hose mode is linear with the resistivity in the low resistivity regime, and the frequency of the mode is much larger than the growth rate. Its existence in the absence of plasma flow distinguishes it from any MHD instability. Since the hose mode is dominant when the background plasma resistivity is high, it may appear in post-TQ tokamak plasmas.

For the RE beam profile considered here, the hose mode grows faster than the tearing mode at a normalized resistivity of $ \eta \u2248 2 \xd7 10 \u2212 4$. Ignoring for now whether or not this current density profile is representative of the experiment, we may use the plasma parameters from a DIII-D discharge reported in Ref. 18 to estimate the temperature where this transition occurs. Using $ a = 0.5 , V A \u2248 6.6 \xd7 10 6$, and *Z _{eff}* = 3, the Spitzer formula produces a temperature estimate of $ T e \u2248 1.5$ eV. This is in the regime of temperatures reported for the post-TQ plasmas observed in DIII-D.

^{18}At these temperatures, recombination of the cold background plasma would lead to a further increase in resistivity from bulk electron–neutral collisions. Experiments in DIII-D using neutral deuterium injection observe recombination of the background plasma that precedes an instability of the RE beam.

^{22}Even in cases where the hose mode is subdominant to resistive MHD instabilities, it can introduce an

*m*= 1 component into the fluctuation spectrum that would otherwise be unexpected from resistive MHD alone in cases where the minimum

*q*-value is larger than unity.

Future work will aim to find an analytic dispersion relation to gain insight on stability thresholds and the dependence on the equilibrium profiles. We also intend to address the nonlinear dynamics of the hose mode in tokamaks to determine whether it may enhance transport of REs out of the plasma, possibly in combination with MHD instabilities.

The authors thank Dr. Chang Liu for providing data on the tearing mode results with runaways. This work was supported by the U.S. DOE through Grant No. DE-SC00180001.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Alexandre Paul Sainterme:** Investigation (equal); Software (equal); Writing – original draft (lead). **C. R. Sovinec:** Project administration (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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