We consider axisymmetric equilibrium of a tokamak plasma that includes current carried by relativistic runaway electrons (REs). Using a guiding center approach, a qualitative picture of the equilibrium of a pure RE beam is elucidated. In a hot thermal plasma, none of the classical drifts of charged particles contribute to the net field-perpendicular current density, which is purely due to magnetization current. In the case of a runaway beam, however, the curvature drift of REs provides the Lorentz force needed to maintain the centripetal acceleration associated with the relativistic toroidal motion. Two different equilibrium formulations are derived for the general case consisting of a mix of thermal and RE current. At higher RE energies, the shift between flux-surfaces and surfaces of constant generalized angular momentum of REs in such equilibria can exceed the radial extent of a typical magnetohydrodynamic mode such that its stability properties could be altered. Simplified one-dimensional governing equations are derived for the absolute and relative orbit shifts in the case of a circular tokamak, enabling quick estimates of parameter dependencies.

## I. INTRODUCTION

During the operation of a tokamak, large scale magnetohydrodynamic (MHD) plasma instabilities can sometimes result in the reorganization of the global magnetic field topology, leading to a disruption. This manifests as a stochastization of the magnetic field (Gorbunov and Razumova, 1964; Waddell , 1979), causing fast loss of thermal energy on a millisecond timescale. Such a thermal quench significantly increases the electrical resistivity and leads to the resistive decay of the plasma current. The large toroidal electric field associated with the increased resistivity can cause acceleration of suprathermal electrons to relativistic energies. Large-angle Coulomb collisions with thermal electrons can cause exponential growth of the small seed of relativistic electrons, eventually leading to a state wherein all the plasma current is carried by a beam of high energy runaway electrons (REs) (Jayakumar , 1993). In fusion-grade tokamaks like ITER and beyond, the final RE current is estimated to be a significant fraction of the pre-disruption current (Hender , 2007; Martín-Solís , 2017; and Vallhagen , 2020). An uncontrolled loss of such a large RE beam can cause deep melting of plasma facing components (Nygren , 1997; Reux , 2015; and Matthews , 2016) and, subsequently, long (and expensive) machine outages. In spite of considerable efforts to establish disruption mitigation schemes, which prevent the formation of an RE beam, it is still unclear whether this goal can be achieved in all circumstances (Vallhagen , 2020). Several schemes have been proposed to mitigate the risk associated with post-disruption REs in ITER, most of which rely on a distributed loss of REs triggered by some form of MHD instability (Hender , 2007; Breizman , 2019; Paz-Soldan , 2019; 2021; and Reux , 2021). To evaluate the efficacy of these schemes, it is important to understand the equilibrium and stability properties of an RE beam, which is the focus of the present work.

In a hot plasma in axisymmetric equilibrium, wherein all the plasma current is carried by thermal electrons, it is well known that particle drifts do not contribute to the net perpendicular current density and, hence, the minor-radial equilibrium. This is due to the fact that the contributions to the perpendicular current density from $ \u2207 B$ drift and curvature drift get exactly canceled by a part of the magnetization current density (Freidberg, 2007). It is the remaining part of the magnetization current density [equal to $ B \u2212 1 ( \u2207 p th \xd7 b \u0302 )$], which provides the Lorentz force $ J \xd7 B$ that balances the local pressure gradient to maintain equilibrium. In the above, $ b \u0302$ is the unit vector along the magnetic field ** B**, and

**and $ p th$ are the current density and thermal pressure-gradient, respectively.**

*J*Immediately after a thermal quench, the axisymmetric equilibrium can still be interpreted in the same manner as described above. Of course, in this specific case, due to negligible $ \u2207 p th$ after the thermal quench, the remnant magnetization current density is negligible as well and, hence, ** J** becomes nearly parallel to

**, which corresponds to a force-free equilibrium, i.e., $ J \xd7 B = 0$. However, when a significant portion of the plasma current is carried by REs, the equilibrium of the cold plasma cannot be interpreted anymore as force-free. The properties of the equilibrium differ significantly, more so at higher RE energies. An important difference is the major-radial shift between the flux surfaces and the RE orbit surfaces, which is proportional to the RE energy (Yoshida, 1990; Abdullaev , 2016). This means that at higher RE energies, the shift can be much larger than the radial extent of a tearing/kink mode, which can potentially alter the underlying physics of the MHD instability. To our knowledge, the special properties of an RE beam equilibrium are presently not accounted for in equilibrium reconstruction codes such as EFIT (Fitzgerald , 2013; Lao , 2022). In the recent years, motional Stark effect (MSE) measurements of the polarization of synchrotron radiation from REs (without any injection of diagnostic neutral beams) have been utilized to obtain useful information on RE pitch angle in post-disruption plasmas (Tinguely , 2019; Popovic , 2021). Nevertheless, pitch angle constraints from MSE being limited, have and might not be able to provide sufficient information for equilibrium reconstruction. Furthermore, the ability to have a self-consistent equilibrium (via reconstruction from experimental data or otherwise) allows more accurate predictions from MHD codes. This is the motivation of the present work.**

*B*Plasma equilibrium has been studied previously in topologically equivalent configurations starting with the Vlasov–Maxwell equations and transforming the distribution function into RE energy and canonical angular momentum space. In particular, this was examined in the context of E-layer equilibria formed by circulating REs in the ASTRON linear mirror device (Marder and Weitzner, 1970) and also in the context of stormtime ring current in earth's dipole field (Lackner, 1970). For RE beam equilibrium in the tokamak configuration, Yoshida (1989) presented a formulation based on an RE fluid velocity stream function. In this paper, we examine the equilibrium problem in further perspectives and provide alternative formulations for computing the equilibrium. In addition, a one-dimensional formulation is derived to estimate the total and relative major-radial shift of the angular momentum surfaces in a circular tokamak in the case of uniform RE energy. In the present treatment, we do not consider RE current decay or loss of REs from the plasma. In order to focus on the standalone effect of REs on the equilibrium and to have a simpler treatment, the effect of background plasma rotation is also not considered in the present study.

This paper is organized as follows: In Sec. II, a qualitative description of a pure RE beam equilibrium is presented. This is followed in Sec. III by the derivation of an equilibrium model for the general case when the plasma current contains both RE and thermal current. In Sec. IV, we derive simple one-dimensional governing equations to evaluate the total and relative RE orbit shifts in a circular tokamak. This is followed by closing comments in Sec. V.

## II. QUALITATIVE PICTURE OF A PURE RE BEAM EQUILIBRIUM

### A. RE fluid momentum balance

*ds*along the field line, $ e n \u0302$ is the unit vector along the inward-normal direction to the curve, and

*R*is the local radius of curvature of the curved motion. Also, the curvature vector $ \kappa \u0302$ is anti-parallel to the radius of curvature vector $ R c$ and is given by

_{c}*n*to be the RE number density, $ v r$ the RE velocity vector, and

_{r}*γ*the relativistic factor, the momentum balance for the RE species can be written as follows:

**and**

*E***are the electric and magnetic fields, respectively. We do not consider radiation losses in this analysis. Using Eq. (2.3), we can simplify the above to**

*B*The steady-state limit of (2.7) determines the equilibrium of the RE beam. There are a couple of ways to interpret this ( Appendix A). A better understanding of the same can be obtained by resorting to the guiding center picture, which we now turn to.

### B. RE equilibrium from guiding center formalism

**and evaluate the net current density due to REs through the area. The RE particle velocity due to the guiding center drifts can be written as follows:**

*B**e*is the electron charge. The three terms in the above equation correspond to the $ E \xd7 B , \u2009 \u2207 B$ and curvature drifts, respectively. It can be easily seen that the polarization drift is relatively negligible and, hence, is not considered. In principle, to compute the current density from the guiding center motion, integration using a distribution function is necessary. However, for the present purpose, as a simplifying approximation, the current density due to guiding center drifts is computed as $ J rg = \u2212 e n r v rg$, which takes the form

*r*from the edge contribute to this magnetization current (Freidberg, 2007). The current due to a single gyrating RE particle is $ I e = \u2212 e \omega c 2 \pi $, and the volume element (that goes all around the edge of the considered area) containing such particles is $ d r = \pi r L 2 b \u0302 \xb7 dl$. Using these, the magnetization current normal to the area due to all such particles in the volume element can be evaluated as follows:

_{L}*ω*is the gyro frequency. Therefore, the magnetization current density is simply

_{c}**and rearranging gives**

*B***is the total current density and $ p th$ is the thermal pressure. We include $ p th$ to keep the equations more general, which can later also cater to the more common special case of $ p th \u2248 0$ in post-thermal quench plasmas. Equations (2.14) and (2.15) constitute the force balance in a plasma with RE current. Adding the above two equations gives**

*J**Ze*and number density

*n*, $ J \u22a5 \u2260 Zen v drifts$ because the magnetization current contribution must be taken into account as well. However, in the specific case of REs, due to the fact that $ \u2207 B$ drift gets canceled and $ p \u22a5 \u226a p | |$, one can safely write

For this reason, one can conclude that in the context of equilibrium, effectively, the RE fluid perpendicular velocity is well approximated by just the curvature drift velocity. This in turn implies that the $ \u2207 B$ drift of REs has an insignificant role in the fluid treatment of REs and can be safely neglected in MHD codes.

## III. GENERAL AXISYMMETRIC EQUILIBRIUM FORMULATION IN THE PRESENCE OF REs

**and**

*B***:**

*J**ψ*is the poloidal magnetic flux, $ e \u0302 \varphi $ is the unit vector in the toroidal direction, and $ \Delta * \psi = \u2202 R R \psi \u2212 1 R \u2202 R \psi + \u2202 Z Z \psi $. The curvature vector $\kappa $ can be suitably decomposed into three components, along the directions $ e \u0302 R , \u2009 \u2207 \psi $, and $ e \u0302 \varphi $, where $ e \u0302 R$ is the unit vector in the major radial direction. It can be easily shown that the components along the $ \u2207 \psi $ and the $ e \u0302 \varphi $ directions are at least an order of magnitude smaller as compared to the major radial component (proof not shown here for the sake of brevity). Therefore, to leading order, only the major radial component of the curvature vector is important, and so $ \kappa = \u2212 b \u0302 \xd7 ( \u2207 \xd7 b \u0302 ) \u2248 \u2212 1 R e \u0302 R$. Additionally, from Eq. (2.15), the well known property follows that $ p th$ is a flux function

There are different possible choices for the form of the RE fluid velocity. Depending on the choice, we can describe the equilibrium via different sets of equations. We hereby briefly describe two formulations.

### A. Formulation-A

*u*is the stream function of the RE current density. The toroidal component of the RE current density is approximated by $ J \varphi , r \u2248 \u2212 e n r v \u2225$, and the parallel velocity is related to the relativistic factor via $ \gamma m = m e 0 1 \u2212 v \u2225 2 / c 2$.

*ψ*or in other words,

*u*. Also the function $ F ( \psi )$ includes the toroidal field from external coils and that from the poloidal thermal current.

*u*, that is $ \gamma = \gamma ( u )$. By taking the momentum balance for the REs

*A*is the generalized angular momentum of an RE particle about the major axis. Clearly, the shift between

*u*and

*ψ*surfaces increases with RE energy.

*u*has the advantage that some of the properties of the equilibrium are very apparent, and the inputs are in a convenient form. Equation (3.13) represents the force balance in the direction perpendicular to the flux-surfaces applied only to the thermal part of the current density. Physically, this representation is justified due to the fact that only the thermal part of the current density lie on the flux-surfaces as shown earlier. The last term on the right hand side of Eq. (3.13) along with the term on the left hand side represents the thermal part of the toroidal current density (and, hence, the Lorentz force due to its crossing with the poloidal magnetic field). The additional term $ e \mu 0 u$ (that appears in the second term on the right hand side of the equation) represents the modification to the toroidal magnetic field due to the RE contribution to the poloidal current density. On the other hand, Eq. (3.15) represents the major-radial force balance purely on the runaway part of the current. In Eq. (3.15), the right hand side term represents the centripetal force necessary to maintain the toroidal motion of the energetic REs, which is provided by the Lorentz force ( $ J r \xd7 B$) arising from RE motion as represented by the terms on the left hand side. The first term on the left hand side of (3.15) represents the Lorentz force due to the crossing of the RE current density in the vertical direction ( $ J Z , r$) and the toroidal magnetic field ( $ B \varphi $), while the second term on the left hand side represents the Lorentz force due to the crossing of the RE current density in the toroidal direction ( $ J \varphi , r$) and the vertical component of the magnetic field (

*B*). Finally, Eq. (3.14) that represents the definition of the generalized angular momentum

_{Z}*A*(

*u*) shows how the RE energy mediates a shift between the flux-surfaces and RE drift orbit surfaces. At low RE energies, one can see that

*A*(

*u*) tends toward $ e \psi $ and the separation between the flux-surfaces and surfaces of constant

*u*becomes negligible. In the special case when all the current is carried only by REs, $ p th = 0$ and $ F ( \psi )$ becomes a pure constant

*F*

_{0}.

### B. Formulation-B

*γ*,

*ψ*, and $ B \varphi $, the RE fluid velocity vector is completely determined. Since

*γ*is supposed to be an input parameter, only

*ψ*and $ B \varphi $ remain as unknowns in the problem. Suitable equations to evaluate them can be obtained as follows. From the background thermal plasma force balance, we get

In the special case when all the current is carried only by REs, $ p th = 0 , \u2009 F ( \psi )$ becomes a pure constant *F*_{0}, and the above two governing equations become decoupled. So one can solve for *ψ* first using $ \Delta * \psi = \mu 0 R ( e n r v | | )$ and then, subsequently, solve for $ B \varphi $.

This formulation (B) has the disadvantage that it is often difficult/impractical to provide the inputs $ \gamma ( R , Z )$ and $ n r ( R , Z )$ that actually correspond to an equilibrium situation. However, in the event of a known RE density and energy distribution (either from experiment or otherwise), this is the most straightforward way to compute the corresponding equilibrium magnetic field. Furthermore, in the special case when REs have a spatially uniform energy, this formulation reduces to much simpler forms that are very useful in practice. We now describe the uniform RE energy cases of formulation-B in the following.

*γ*, and $ N ( A )$:

*γ*) and all the current is carried only by REs. In this case, the two governing equations become decoupled. So with inputs

*γ*and $ N ( A )$, one can first solve for

*A*(

*R*,

*Z*) via Eq. (3.27) and then subsequently compute $ B \varphi ( R , Z )$ via

*R*= 6.2 m,

*a*= 2 m, and $ I p = I RE = 10$ MA. The computation was done using the (

*r*, $\theta $) co-ordinates, where

*r*and

*θ*represent the minor-radial co-ordinate and poloidal angle, respectively. In these co-ordinates, Eq. (3.27) can be expressed as follows:

*R*

_{0}is the major radius at the center of the circular domain. Runaway electron energies of 40 and 80 MeV have been chosen. The RE current density profile peaking is characterized via a spatial Gaussian function

*w*is a measure of the peaking for which we choose

*w*= 0.05. For the angular momentum

*A*, a fixed boundary condition $ A = \gamma m v | | R$ has been used that corresponds to

*ψ*= 0. Second order finite-differences are used for spatial discretization of (3.27) on a uniform polar grid of size 64 × 64. The solution of the resulting linear system of equations was obtained using LU decomposition that gives $ A ( r , \theta )$ at all grid points in the domain. Once

*A*is known, $ \psi ( r , \theta )$ is computed at all grid points directly using $ e \psi = \gamma m v | | R \u2212 A$.

The solution is shown in Figs. 1(a) and 1(b) via surfaces of constant *ψ* and constant *A*. As expected, the shift between *ψ* and *A* surfaces is higher at larger RE energy. The shift at the boundary is approximately $ \u223c 8$ cm for the 40 MeV case and $ \u223c 15$ cm for the 80 MeV case.

## IV. QUANTITATIVE ESTIMATION OF FLUX SURFACE AND RE ORBIT SHIFTS

*r*(circular surfaces) and a first order component that is also dependant on

*θ*. That is,

*A*. We also define the relative shift $ \delta ( r )$ between

*A*surfaces and

*ψ*surfaces as follows:

*r*=

*a*simplifies to

*δ*.

*δ*as follows:

Numerical solutions for Eqs. (4.11)–(4.14) have been obtained for an ITER-like circular plasma with $ I p 0 = I RE 0 = 6$ MA, $ R 0 = 6.2$ m, and *a* = 2 m and a leading-order current density profile $ J 0 ( r ) = J \u0302 0 ( 1 \u2212 r 2 a 2 ) \nu $. Runaway electron energies in the range 2–100 MeV and $ 2 \u2264 \nu \u2264 8$ (which translates to a normalized internal inductance $ 1.23 \u2264 l i \u2264 2.17$) have been considered. Figure 2(a) shows the minor-radial distribution of the flux-surface and RE orbit shifts for the case with the RE energy of 50 MeV and *ν* = 2. The relative shift *δ* clearly increases as we move minor-radially outward. Figures 2(b) and 2(c) show the dependence of $ \Delta S = \Delta ( r = 0 )$ vs the normalized internal inductance *l _{i}* and RE energy. From Fig. 2(b), it is seen that $ \Delta S$ increases with both the increase in the profile peaking (or

*l*) and RE energy. At very low RE energies, the solution tends to the limiting case of the classical Shafranov shift in a force-free equilibrium $ \Delta S FF$ indicated by the black curve in Fig. 2(b). The dependence of the minor-radially averaged relative shift $ \delta avg = a \u2212 1 \u222b 0 a \delta d r$ on RE energy and

_{i}*l*is shown in Fig. 2(d). It is interesting to note that an increase in peaking of the current density profile leads to a decrease in $ \delta avg$, qualitatively in line with previous numerical predictions (Yoshida, 1990). Similar qualitative behavior was observed from the 2D solution of Eq. (3.27).

_{i}## V. SUMMARY AND DISCUSSION

The problem of axisymmetric tokamak plasma equilibrium in the presence of relativistic electron current has been analyzed in detail. Using the guiding center formalism, it is shown that for the RE beam equilibrium, the current density due to the RE curvature drift provides the necessary centripetal acceleration to maintain the relativistic toroidal motion. This changes the equilibrium properties significantly at high RE energies, as compared to the force-free equilibrium of a cold plasma. Two different formulations (simpler than proposed in an earlier study) are provided for the general case of equilibrium consisting of both thermal and RE current. Unlike the case of a thermal plasma, in the presence of REs, $ R B \varphi $ is not a pure function of *ψ* anymore, which implies that lines of the RE current density do not lie on magnetic flux surfaces. Both the formulations are applicable for the general cases in which both thermal and RE current co-exist and when the spatial energy distribution of REs is non-uniform. In such general cases, further simplification is not possible, and solutions can be obtained by solving the equations numerically. Nevertheless, in certain special cases such as uniform RE energy and/or all current carried by REs, decoupling of equations as well as significant simplification of the equations can be achieved. For the case of a circular tokamak with all the current carried by monoenergetic REs, simple one-dimensional formulations for the total and relative RE orbit shifts have been presented. This enables quick and approximate estimates of the radial distribution of shifts as a function of the aspect ratio, normalized internal inductance, and RE energy.

The work presented here is useful in improving RE fluid models used in fusion MHD codes that often assume a force-free equilibrium even in the presence of REs (Munaretto , 2020; Zhao , 2020; and Bandaru , 2021). The RE fluid model (Bandaru , 2019) in the non-linear MHD code JOREK (Huysmans and Czarny, 2007; Hoelzl , 2021) has been extended to incorporate the correct equilibrium as outlined in this paper (details not shown in this paper for the sake of brevity). We consider this work to be also helpful for tokamak equilibrium reconstruction codes that currently do not take the relevant effects into account. An important consequence of the modified equilibrium with REs is its potential effect on the stability of MHD modes at high RE energies due to the fact that the RE orbit-surface shift can then be comparable or even larger than the MHD mode width. For example, in an ITER post-disruption plasma at $ T e \u223c 20$ eV without impurities, the resistive layer width would be of the order of a few centimeters and would be comparable to the average relative shift in a 20 MeV RE energy situation as shown in Fig. 2(d). The effect of the shifts on tearing mode stability is currently being studied and will be communicated in a future publication.

## ACKNOWLEDGMENTS

We acknowledge fruitful discussion with Karl Lackner. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. Part of this work was supported by the EUROfusion-Theory and Advanced Simulation Coordination (E-TASC) via the Theory and Simulation Verification and Validation (TSVV) Project on MHD Transients (2021–2025). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Vinodh Bandaru:** Conceptualization (equal); Formal analysis (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (equal). **Matthias Hoelzl:** Conceptualization (equal); Formal analysis (supporting); Methodology (supporting); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

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

### APPENDIX A: PHYSICAL INTERPRETATION OF THE STATE-STATE PERPENDICULAR MOMENTUM EQUATION

### APPENDIX B: DERIVATION OF EQ. (38) FROM RE MOMENTUM BALANCE

*e*give

*u*, the gradient of $ \gamma m v \varphi / n r$ will be parallel to $ \u2207 u$. So, it can be taken into the gradient operator along with

*R*. This gives

## REFERENCES

*Plasma Physics and Fusion Energy*

^{1}