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.

Topics
Magnetohydrodynamics, Plasma dynamics, Runaway electrons, Tokamaks

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 ∇ 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 − 1 ( ∇ p th × b ̂ )], which provides the Lorentz force J × B that balances the local pressure gradient to maintain equilibrium. In the above, b ̂ is the unit vector along the magnetic field B, and J and p th are the current density and thermal pressure-gradient, respectively.

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 ∇ p th after the thermal quench, the remnant magnetization current density is negligible as well and, hence, J becomes nearly parallel to B, which corresponds to a force-free equilibrium, i.e., J × 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.

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.

The acceleration corresponding to the parallel motion of an RE particle along a field line can be expressed as follows:
d d t ( v ∥ b ̂ ) = d v ∥ d t b ̂ + v ∥ d b ̂ d θ d θ d s d s d t = d v ∥ d t b ̂ + v ∥ e n ̂ 1 R c v ∥ .
(2.1)
Here, v | | is the component of the RE velocity along the magnetic field vector, d θ is the angle covered by the particle as it traverses a distance ds along the field line, e n ̂ is the unit vector along the inward-normal direction to the curve, and Rc is the local radius of curvature of the curved motion. Also, the curvature vector κ ̂ is anti-parallel to the radius of curvature vector R c and is given by
κ ̂ = R c R c 2 = R c ̂ R c ,
(2.2)
and e n ̂ = − R c ̂⁠, which leads us to e n ̂ / R c = − R c ̂ / R c = κ ̂⁠. Using this result, Eq. (2.1) can be simplified as follows:
d d t ( v ∥ b ̂ ) = b ̂ d v ∥ d t + v ∥ 2 κ ̂ .
(2.3)
If we choose nr to be the RE number density, v r the RE velocity vector, and γ the relativistic factor, the momentum balance for the RE species can be written as follows:
d d t [ γ m n r v r ] = − e n r [ E + v r × B ] ,
(2.4)
where the symbol γ m = γ m e , 0 ,   m e 0 denotes the electron rest mass and E and B 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
γ m n r ( b ̂ d v ∥ d t + v ∥ 2 κ ̂ + d v ⊥ d t ) + v r d d t ( γ m n r ) = − e n r [ E + v r × B ] .
(2.5)
The momentum equation (2.5) can be decomposed into a parallel component and a perpendicular component as follows:
d d t ( γ m n r v ∥ ) = − e n r E ∥ ,
(2.6)
d d t ( γ m n r v ⊥ ) + γ m n r v ∥ 2 κ ̂ = − e n r E ⊥ + j r × B .
(2.7)

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.

We consider a small area for which the normal vector is perpendicular to B 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:
v rg = E × B B 2 − γ m v ⊥ 2 2 e B × ∇ B B 3 + γ m v ∥ 2 e κ × B B 2 ,
(2.8)
where v ⊥ and v ∥ are the perpendicular and parallel RE particle velocity, respectively, and e is the electron charge. The three terms in the above equation correspond to the E × B ,   ∇ 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 = − e n r v rg⁠, which takes the form
J rg = − e n r E × B B 2 + p ⊥ B × ∇ B B 3 − p ∥ κ × B B 2 ,
(2.9)
where the variables p ⊥ = n r γ m v ⊥ 2 2 ,   p ∥ = n r γ m v ∥ 2⁠, denote the dynamic pressure in the directions perpendicular and parallel to the magnetic field.
We now turn to the current density due to gyro motion. It is well known that gyro motion can contribute to a net perpendicular current density only close to the edge of the small area considered. Specifically, particles that have their guiding centers within a distance of the gyro radius rL from the edge contribute to this magnetization current (Freidberg, 2007). The current due to a single gyrating RE particle is I e = − e ω c 2 π⁠, and the volume element (that goes all around the edge of the considered area) containing such particles is d r = π r L 2 b ̂ · dl⁠. Using these, the magnetization current normal to the area due to all such particles in the volume element can be evaluated as follows:
I r m = ∮ I e n r d r = ∮ − e ω c 2 π × n r × π γ m 2 v ⊥ 2 e 2 B 2 b ̂ · dl = − ∮ n r γ m v ⊥ 2 2 B b ̂ · dl = − ∮ p ⊥ B b ̂ · dl ,
where ωc is the gyro frequency. Therefore, the magnetization current density is simply
J rm = − ∇ × ( p ⊥ B b ̂ ) ,
which can be rewritten as follows:
J rm = − p ⊥ B × ∇ B B 3 + p ⊥ κ × B B 2 − 1 B ( ∇ p ⊥ × b ̂ ) .
(2.10)
The net current density would then be
J r = J rg + J rm = − e n r E × B B 2 + ( p ∥ − p ⊥ ) B × κ B 2 + 1 B ( b ̂ × ∇ p ⊥ ) .
(2.11)
Note that the current density due to the ∇ B drift canceled exactly, but the current density due to curvature drift does not cancel due to the strong anisotropy in the RE velocity. This is the fundamental difference of the RE beam current as compared to the classical thermal current. Doing a cross-product with B and rearranging gives
( p ∥ − p ⊥ ) κ = − e n r E ⊥ + J r × B − ∇ ⊥ p ⊥ .
(2.12)
Experiments and kinetic studies indicate that the RE pitch angle (⁠ v ⊥ / v⁠) is typically very small (Breizman , 2019) (less than 10° in most case, which is p ⊥ / p ∥   ≲   3 %⁠). There are exceptions, however, to this. RE-wave interactions can lead to higher pitch-angles. In this treatment, however, we discount such high pitch-angle scenarios by assuming p ∥ ≫ p ⊥⁠, which simplifies the above equation to
p ∥ κ = − e n r E ⊥ + J r × B .
(2.13)
This is exactly the same as the steady-state limit of Eq. (2.7). Furthermore, the term involving E ⊥ is of little relevance to equilibrium as all particles have the same E × B drift and does not contribute to the net current density. So the RE equilibrium is mainly determined by
p ∥ κ = J r × B .
(2.14)
The main point here is that the net perpendicular RE current density J r comes from the curvature drift motion of REs. This provides the necessary J r × B force to maintain the centripetal acceleration along the major-radial direction needed for the relativistic toroidal motion of REs.
For the background plasma, the equilibrium is determined by
0 = ( J − J r ) × B − ∇ ⊥ p th ,
(2.15)
where J 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 ≈ 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
p ∥ κ = J × B − ∇ ⊥ p th ,
(2.16)
for the whole plasma (all species included).
Another implication of the above analysis is as follows: In general, for any plasma species with charge Ze and number density n, J ⊥ ≠ 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 ∇ B drift gets canceled and p ⊥ ≪ p | |⁠, one can safely write
J r , ⊥ = − e n r v curv . drift .
(2.17)

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 ∇ B drift of REs has an insignificant role in the fluid treatment of REs and can be safely neglected in MHD codes.

We now aim to determine the magnetic field that corresponds to the axisymmetric equilibrium in the presence of REs. We include the possibility of a significant thermal pressure-gradient in the plasma, although in most cases of post-thermal-quench plasmas, the thermal pressure gradient can be neglected. For the magnetic field, axisymmetry plus solenoidal constraint easily leads to the following well known forms for B and J:
B = 1 R ∇ ψ × e ̂ ϕ + B ϕ e ̂ ϕ ,
(3.1)
J ϕ = − 1 μ 0 R Δ * ψ ,
(3.2)
J pol = 1 μ 0 R ∇ ( R B ϕ ) × e ̂ ϕ ,
(3.3)
where ψ is the poloidal magnetic flux, e ̂ ϕ is the unit vector in the toroidal direction, and Δ * ψ = ∂ R R ψ − 1 R ∂ R ψ + ∂ Z Z ψ⁠. The curvature vector κ can be suitably decomposed into three components, along the directions e ̂ R ,   ∇ ψ⁠, and e ̂ ϕ⁠, where e ̂ R is the unit vector in the major radial direction. It can be easily shown that the components along the ∇ ψ and the e ̂ ϕ 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 κ = − b ̂ × ( ∇ × b ̂ ) ≈ − 1 R e ̂ R⁠. Additionally, from Eq. (2.15), the well known property follows that p th is a flux function
p th = p th ( ψ ) .
(3.4)
In the classical Grad–Shafranov equilibrium, the toroidal field can be represented as R B ϕ = F ( ψ ) and J · ∇ ψ = 0⁠. That is, current density lines lie on flux surfaces. In the presence of REs, this is not the case anymore due to the dependence of the RE current density on the curvature drift.

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.

From number conservation ∇ · ( v r n r ) = 0 and the fact that ∂ ( … ) / ( ∂ ϕ ) = 0⁠, we can write the RE velocity vector as follows:
v r = 1 n r [ 1 R ∇ u × e ϕ − J ϕ , r e e ϕ ] ,
(3.5)
where u is the stream function of the RE current density. The toroidal component of the RE current density is approximated by J ϕ , r ≈ − e n r v ∥⁠, and the parallel velocity is related to the relativistic factor via γ m = m e 0 1 − v ∥ 2 / c 2⁠.
Taking the toroidal component of Eq. (2.15) gives
[ ( J − J r ) × B ] · e ϕ = 0 , [ { ( 1 μ 0 R ∇ ( R B ϕ ) + e R ∇ u ) × e ϕ } × 1 R ( ∇ ψ × e ϕ ) ] · e ϕ = 0 , − [ ( 1 μ 0 R ∇ ( R B ϕ ) + e R ∇ u ) · ( e ϕ × 1 R ∇ ψ ) ] e ϕ · e ϕ = 0 , e ϕ R 2 · [ ∇ ( R B ϕ μ 0 + e u ) × ∇ ψ ] = 0 .
This means that R B ϕ / μ 0 + e u is a pure function of ψ or in other words,
R B ϕ = − e μ 0 u + F ( ψ ) .
(3.6)
This implies that the thermal part of the current density would still lie on the flux-surfaces, while the RE current density lines lie on surfaces of constant u. Also the function F ( ψ ) includes the toroidal field from external coils and that from the poloidal thermal current.
We now assume that the RE energy is a pure function of u, that is γ = γ ( u )⁠. By taking the momentum balance for the REs
v | | 2 κ = − e γ m v r × B ,
(3.7)
approximating v | | ≈ − J ϕ e n r⁠, doing a cross-product with ∇ u and some algebra ( Appendix B for details) leads to
γ m R ( − J ϕ , r e n r ) − e ψ = A ( u ) ,
(3.8)
where 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.
Likewise, taking the major-radial projection of Eq. (3.7) gives easily
( e μ 0 u − F ) R ∂ u ∂ R + n r v ∥ ∂ ψ ∂ R = γ m n r v ∥ 2 e .
(3.9)
We now take the background plasma force-balance equation and do the operation [ ( J − J r ) × B ] · ∇ ψ = ∇ p th · ∇ ψ⁠, which simplifies as follows:
∇ ψ · [ ( 1 μ 0 R ∇ ( R B ϕ ) + e R ∇ u ) × e ϕ × ( B ϕ e ϕ ) + ( − Δ * ψ μ 0 R − J ϕ , r ) e ϕ × 1 R ( ∇ ψ × e ϕ ) ] = p ′ th ( ∇ ψ ) 2 , ∇ ψ · [ − B ϕ ( 1 μ 0 R ∇ ( R B ϕ ) + e R ∇ u ) + ∇ ψ R ( − Δ * ψ μ 0 R − J ϕ , r ) ] = p ′ th ( ∇ ψ ) 2 .
(3.10)
Since R B ϕ = − e μ 0 u + F ( ψ )⁠, we have
∇ ( R B ϕ ) μ 0 = − e ∇ u + F ′ μ 0 ∇ ψ .
(3.11)
Using this in Eq. (3.10), we get
− B ϕ R [ − e ∇ u + F ′ μ 0 ∇ ψ + e ∇ u ] · ∇ ψ + ( ∇ ψ ) 2 R ( − Δ * ψ μ 0 R − J ϕ , r ) = p ′ th ( ∇ ψ ) 2 .
Multiplying throughout by − R 2 μ 0 ( ∇ ψ ) 2 and rearranging, we get
Δ * ψ = − μ 0 R 2 p ′ th − μ 0 R J ϕ , r − R B ϕ F ′ .
In the above, using R B ϕ = − e μ 0 u + F and J ϕ , r = − e n r v ∥ gives
Δ * ψ = − μ 0 R 2 p ′ th − ( F − e μ 0 u ) F ′ + μ 0 R ( e n r v | | ) .
(3.12)
To summarize, given the inputs p ′ th ( ψ ) ,   F ( ψ ) ,   γ m ( u )⁠, and A ( u )⁠, one needs to solve for ψ ( R , Z ) ,   u ( R , Z )⁠, and n r ( R , Z )⁠, the equations
Δ * ψ = − μ 0 R 2 p ′ th − ( F − e μ 0 u ) F ′ + μ 0 R ( e n r v | | ) ,
(3.13)
e ψ = γ m v | | R − A ,
(3.14)
( e μ 0 u − F ) R ∂ u ∂ R + n r v ∥ ∂ ψ ∂ R = γ m n r v ∥ 2 e .
(3.15)
After this, B ϕ is easily evaluated from the expression R B ϕ = − e μ 0 u + F⁠. Clearly, the above set of equations have to be solved together iteratively. This formulation (A) with the use of the RE current density streamfunction 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 μ 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 × 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 ϕ⁠), 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 ϕ , r⁠) and the vertical component of the magnetic field (BZ). Finally, Eq. (3.14) that represents the definition of the generalized angular momentum 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 ψ 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 ( ψ ) becomes a pure constant F0.
Using the guiding center formalism, we have previously shown that the RE fluid perpendicular velocity is effectively the curvature drift in the context of axisymmetric equilibrium. In this respect, we could take the runaway fluid velocity to be
v r = v | | R B ϕ [ ( ∇ ψ − γ m v | | e ∇ R ) × e ϕ ] + v | | e ϕ .
(3.16)
One can easily show that the above form satisfies the density conservation ∇ · ( n r v r )⁠. This is apparent that the above form of the RE fluid velocity can also be obtained by recasting the streamfunction form [i.e., Eq. (3.5)] by simplifying the RE force balance as follows:
v | | 2 κ = − e γ m v r × B , v | | 2 R e R = e γ m R n r [ − B ϕ ∇ u − J ϕ e ∇ ψ ] .
Using e R = ∇ R ,   J ϕ = − e n r v | | and multiplying by γ m / e give
γ m v | | 2 e R ∇ R = − B ϕ R 1 n r ∇ u + v | | R ∇ ψ ,
which on rearranging and doing a cross-product with e ϕ gives
1 n r ( 1 R ∇ u × e ϕ ) = v | | R B ϕ [ ( ∇ ψ − γ m v | | e ∇ R ) × e ϕ ] .
(3.17)
So both the forms of the RE fluid velocity (formulation-A and B) are fully consistent with each other.
It can be seen that knowing γ, ψ, and B ϕ⁠, the RE fluid velocity vector is completely determined. Since γ is supposed to be an input parameter, only ψ and B ϕ 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
e ϕ · [ ( J − J r ) × B ] = ∇ p th · e ϕ , e ϕ · [ [ ( 1 μ 0 R ∇ ( R B ϕ ) + e n r v | | R B ϕ ( ∇ ψ − γ m v | | e ∇ R ) ) × e ϕ ] × ( 1 R ∇ ψ × e ϕ ) ] = 0 , e ϕ R 2 · [ ( ∇ ( R B ϕ ) + μ 0 e n r v | | B ϕ ( ∇ ψ − γ m v | | e ∇ R ) ) × ∇ ψ ] = 0 .
Therefore, we can write
∇ ( R B ϕ ) + μ 0 e n r v | | B ϕ ( ∇ ψ − γ m v | | e ∇ R ) = F ′ ∇ ψ ,
(3.18)
where F = F ( ψ )⁠. Multiplying throughout by R B ϕ⁠, we get
R B ϕ ∇ ( R B ϕ ) + μ 0 R e n r v | | ( ∇ ψ − γ m v | | e ∇ R ) = R B ϕ F ′ ∇ ψ .
(3.19)
Furthermore, we do the operation [ ( J − J r ) × B ] · ∇ ψ = ∇ p th · ∇ ψ⁠, which simplifies as follows:
∇ ψ · [ { 1 μ 0 R ∇ ( R B ϕ ) + e n r v | | R B ϕ ( ∇ ψ − γ m v | | e ∇ R ) } × e ϕ × ( B ϕ e ϕ ) ] + ∇ ψ · [ ( − Δ * ψ μ 0 R − J ϕ , r ) e ϕ × 1 R ( ∇ ψ × e ϕ ) ] = p ′ th ( ∇ ψ ) 2 , ∇ ψ · [ − B ϕ { 1 μ 0 R ∇ ( R B ϕ ) + e n r v | | R B ϕ ( ∇ ψ − γ m v | | e ∇ R ) } + ∇ ψ R ( − Δ * ψ μ 0 R − J ϕ , r ) ] = p ′ th ( ∇ ψ ) 2 , Δ * ψ = − μ 0 R 2 p ′ th − μ 0 R J ϕ , r − ∇ ψ | ∇ ψ | 2 [ R B ϕ ∇ ( R B ϕ ) + μ 0 R e n r v | | ( ∇ ψ − γ m v | | e ∇ R ) ] .
Now, replacing the term in the square brackets with R B ϕ F ′ ∇ ψ [Eq. (3.19)] gives
Δ * ψ = − μ 0 R 2 p ′ th − R B ϕ F ′ + μ 0 R ( e n r v | | ) .
(3.20)
In summary, with the given inputs p ′ th ( ψ ) ,   F ( ψ ) ,   γ ( R , Z )⁠, and n r ( R , Z )⁠, one needs to solve the following two coupled equations for the unknowns ψ ( R , Z ) and B ϕ ( R , Z )⁠:
∇ ( R B ϕ ) + μ 0 e n r v | | B ϕ ( ∇ ψ − γ m v | | e ∇ R ) = F ′ ∇ ψ ,
(3.21)
Δ * ψ = − μ 0 R 2 p ′ th − R B ϕ F ′ + μ 0 R ( e n r v | | ) .
(3.22)

In the special case when all the current is carried only by REs, p th = 0 ,   F ( ψ ) becomes a pure constant F0, and the above two governing equations become decoupled. So one can solve for ψ first using Δ * ψ = μ 0 R ( e n r v | | ) and then, subsequently, solve for B ϕ⁠.

This formulation (B) has the disadvantage that it is often difficult/impractical to provide the inputs γ ( 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.

When REs have a spatially uniform energy, that is when γ m and v ∥ are pure constants, formulation-B can be recast into a practically useful form (in terms of input needed) due to the fact that the RE fluid velocity given in Eq. (3.16) can now be rewritten as follows:
v r = − v | | e R B ϕ ( ∇ A × e ϕ ) + v | | e ϕ ,
(3.23)
with A = γ m v | | R − e ψ⁠. To the leading order (under the assumption that R B ϕ is approximately constant), one can see that ∇ · ( v r / R ) = 0⁠. This along with the fact that ∇ · ( v r n r ) = 0 leads easily to ∇ · [ R n r ( v r / R ) ] = 0⁠, and hence, R n r = N ( A )⁠. The solution can then be obtained by solving the following equations for ψ ( R , Z ) and B ϕ ( R , Z ) with inputs p ′ th ( ψ ) ,   F ( ψ )⁠, γ, and N ( A )⁠:
∇ ( R B ϕ ) − μ 0 v | | N R B ϕ ( γ m v | | ∇ R − e ∇ ψ ) = F ′ ∇ ψ ,
(3.24)
Δ * ψ = − μ 0 R 2 p ′ th − R B ϕ F ′ + μ 0 e v | | N .
(3.25)
Further simplification occurs, when REs have a spatially uniform energy (constant γ) 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 ϕ ( R , Z ) via
∇ ( R B ϕ ) − μ 0 v | | N R B ϕ ∇ A = 0 ,
(3.26)
Δ * A = − μ 0 e 2 v | | N − γ m v | | R .
(3.27)
After that, it is trivial to compute ψ ( R , Z ) from A = γ m v | | R − e ψ⁠.
Numerical solution for this specific case (all current carried by REs with uniform energy) has been obtained for an ITER-like circular plasma with R = 6.2 m, a = 2 m, and I p = I RE = 10 MA. The computation was done using the (r, θ⁠) 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:
( 1 r ∂ ∂ r r ∂ ∂ r + 1 r 2 ∂ 2 ∂ θ 2 ) A − 1 R 0 + r   cos   θ ( cos   θ ∂ ∂ r −   sin   θ r ∂ ∂ θ ) A = − μ 0 e 2 v | | N ( A ) − γ m v | | R 0 + r   cos   θ ,
(3.28)
where R0 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
N ( A n ) = R n r 0 e − A n 2 w ,
(3.29)
where A n = A / e is the normalized angular momentum and w is a measure of the peaking for which we choose w = 0.05. For the angular momentum A, a fixed boundary condition A = γ 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 , θ ) at all grid points in the domain. Once A is known, ψ ( r , θ ) is computed at all grid points directly using e ψ = γ m v | | R − 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 ∼ 8 cm for the 40 MeV case and ∼ 15 cm for the 80 MeV case.

FIG. 1.
Contours of the surfaces of constant ψ (dotted lines) and constant A (solid lines) for the equilibrium solution in the case of the pure RE current and a spatially constant γ in an ITER-sized circular plasma. The cases correspond to (a) 40 MeV, w = 0.05; (b) 80 MeV, w = 0.05, where w is a measure of the profile broadness. In both the cases, the total RE current (which is also the total plasma current) has been kept constant at Ip = 10 MA.
View largeDownload slide

Contours of the surfaces of constant ψ (dotted lines) and constant A (solid lines) for the equilibrium solution in the case of the pure RE current and a spatially constant γ in an ITER-sized circular plasma. The cases correspond to (a) 40 MeV, w = 0.05; (b) 80 MeV, w = 0.05, where w is a measure of the profile broadness. In both the cases, the total RE current (which is also the total plasma current) has been kept constant at Ip = 10 MA.

FIG. 1.
Contours of the surfaces of constant ψ (dotted lines) and constant A (solid lines) for the equilibrium solution in the case of the pure RE current and a spatially constant γ in an ITER-sized circular plasma. The cases correspond to (a) 40 MeV, w = 0.05; (b) 80 MeV, w = 0.05, where w is a measure of the profile broadness. In both the cases, the total RE current (which is also the total plasma current) has been kept constant at Ip = 10 MA.
View largeDownload slide

Contours of the surfaces of constant ψ (dotted lines) and constant A (solid lines) for the equilibrium solution in the case of the pure RE current and a spatially constant γ in an ITER-sized circular plasma. The cases correspond to (a) 40 MeV, w = 0.05; (b) 80 MeV, w = 0.05, where w is a measure of the profile broadness. In both the cases, the total RE current (which is also the total plasma current) has been kept constant at Ip = 10 MA.

Close modal
In this section, we derive differential equations that describe the major-radial shift in the magnetic flux surfaces and the RE drift orbit surfaces. Note that even the shift in the flux surfaces can be very different from what might be estimated from a classical thermal equilibrium. Here, we consider the case of a circular tokamak cold plasma, wherein all the current is carried by REs of uniform energy. As discussed in Sec. III B, the following equations describe the state of axisymmetric equilibrium in such a scenario:
Δ * ψ = μ 0 e v | | N ( A ) ,
(4.1)
Δ * A = − μ 0 e 2 v | | N ( A ) − γ m v | | R .
(4.2)
In the ( r , θ⁠) cylindrical co-ordinates, the above equations easily transform to
( 1 r ∂ ∂ r r ∂ ∂ r + 1 r 2 ∂ 2 ∂ θ 2 ) ψ − 1 R 0 + r   cos   θ ( cos   θ ∂ ∂ r −   sin   θ r ∂ ∂ θ ) ψ = μ 0 e v | | N ( A ) ,
(4.3)
( 1 r ∂ ∂ r r ∂ ∂ r + 1 r 2 ∂ 2 ∂ θ 2 ) A − 1 R 0 + r   cos   θ ( cos   θ ∂ ∂ r −   sin   θ r ∂ ∂ θ ) A = − μ 0 e 2 v | | N ( A ) − γ m v | | R 0 + r   cos   θ ,
(4.4)
where Eq. (4.4) is identical to Eq. (3.28) in Sec. III B. Similar to the treatment for the classical Shafranov shift (Wesson, 2004), we consider the solution for ψ ( r , θ ) and A ( r , θ ) to consist of a leading order component, which is a pure function of r (circular surfaces) and a first order component that is also dependant on θ. That is,
ψ ( r , θ ) = ψ 0 ( r ) + ψ 1 ( r , θ ) = ψ 0 − Δ ( r ) d ψ 0 d R = ψ 0 − Δ ( r )   cos   θ d ψ 0 d r ,
(4.5)
A ( r , θ ) = A 0 ( r ) + A 1 ( r , θ ) = A 0 − Δ RE ( r ) d A 0 d R = A 0 − Δ RE ( r )   cos   θ d A 0 d r .
(4.6)
Here, Δ ( r ) represents the major-radial shift in the flux-surfaces, and Δ RE ( r ) is the corresponding shift in the surfaces of constant A. We also define the relative shift δ ( r ) between A surfaces and ψ surfaces as follows:
δ ( r ) = Δ RE ( r ) − Δ ( r ) .
(4.7)
Using such a decomposition, the leading order and first order balance of Eq. (4.4) can be respectively written as follows:
1 r d d r ( r d A 0 d r ) = − μ 0 e 2 v | | N ( A 0 ) − γ m v | | R 0 ,
(4.8)
( 1 r ∂ ∂ r r ∂ ∂ r + 1 r 2 ∂ 2 ∂ θ 2 ) A 1 −   cos   θ R 0 d A 0 d r = − d d r ( μ 0 e 2 v | | N ( A 0 ) ) d r d A 0 A 1 + γ m v | | R 0 ( r R 0 cos   θ ) .
(4.9)
Using A 1 = − Δ RE ( r )   cos   θ d A 0 d r⁠, Eqs. (4.8) and (4.9) can be simplified to obtain
d d r [ r ( d A 0 d r ) 2 d Δ RE d r ] = r R 0 [ − γ m v ∥ R 0 r d A 0 d r − ( d A 0 d r ) 2 ] .
Now using d A 0 d r = − e d ψ 0 d r and B θ 0 ≈ − R 0 − 1 d ψ 0 d r⁠, one can rewrite the above equation as follows:
d d r [ r B θ 0 2 d Δ RE d r ] = r R 0 [ − γ m v ∥ e r R 0 B θ 0 − B θ 0 2 ] ,
(4.10)
which upon integration gives us the governing differential equation for Δ RE in the form
d Δ RE d r = 1 R 0 B θ 0 2 [ − γ m v ∥ e R 0 1 r ∫ 0 r r 2 B θ 0 d r − 1 r ∫ 0 r r B θ 0 2 d r ] .
(4.11)
Note that unlike Δ, here Δ RE can have a finite value on the boundary, i.e., Δ RE ( a ) ≠ 0⁠. The necessary boundary condition for Δ RE can be obtained from the angular momentum definition as follows:
A = γ m v ∥ R − e ψ , A 0 − Δ RE ( r )   cos   θ d A 0 d r = γ m v ∥ ( R 0 + r   cos   θ ) − e [ ψ 0 − Δ ( r )   cos   θ d ψ 0 d r ] ,
in which the first order balance becomes
− Δ RE ( r )   cos   θ d A 0 d r = γ m v ∥ r   cos   θ + e Δ ( r )   cos   θ d ψ 0 d r .
Using the relations Δ ( a ) = 0 ,   d A 0 d r = − e d ψ 0 d r⁠, and B θ 0 ≈ − R 0 − 1 d ψ 0 d r⁠, the above equation at r = a simplifies to
Δ RE ( a ) = − ( γ m v ∥ e ) a R 0 1 B θ 0 , a .
(4.12)
Equations (4.11) and (4.12) can be used to fully determine the radial distribution of Δ RE⁠. We now proceed to obtain a governing equation for δ.
The first order balance of Eq. (4.3) reads
( 1 r ∂ ∂ r r ∂ ∂ r + 1 r 2 ∂ 2 ∂ θ 2 ) ψ 1 −   cos   θ R 0 d ψ 0 d r = d d r ( μ 0 e v | | N ( A 0 ) ) d r d A 0 A 1 .
Using the leading order balance of Eq. (4.3), i.e., 1 r d d r ( r d ψ 0 d r ) = μ 0 e v | | N ( A 0 )⁠, the above equation simplifies to
d d r ( r B θ 0 2 d δ d r ) + r B θ 0 d d r ( 1 r d d r ( r B θ 0 ) ) δ = d d r ( r B θ 0 2 d Δ RE d r ) + r R 0 B θ 0 2 = − ( r R 0 ) 2 γ m v ∥ e B θ 0 .
The last simplification (on the right hand side above) comes from the use of Eq. (4.10). We, therefore, obtain a second-order ordinary differential equation for δ as follows:
r B θ 0 2 ∂ 2 δ ∂ r 2 + d d r ( r B θ 0 2 ) d δ d r + r B θ 0 d d r ( 1 r d d r ( r B θ 0 ) ) δ = − ( r R 0 ) 2 γ m v ∥ e B θ 0 ,
(4.13)
with the boundary conditions
( d δ d r ) r = 0 = 0 ,   δ ( a ) = − ( γ m v ∥ e ) a R 0 1 B θ 0 , a .
(4.14)

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 ̂ 0 ( 1 − r 2 a 2 ) ν⁠. Runaway electron energies in the range 2–100 MeV and 2 ≤ ν ≤ 8 (which translates to a normalized internal inductance 1.23 ≤ l i ≤ 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 Δ S = Δ ( r = 0 ) vs the normalized internal inductance li and RE energy. From Fig. 2(b), it is seen that Δ S increases with both the increase in the profile peaking (or li) 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 Δ S FF indicated by the black curve in Fig. 2(b). The dependence of the minor-radially averaged relative shift δ avg = a − 1 ∫ 0 a δ d r on RE energy and li 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 δ avg⁠, qualitatively in line with previous numerical predictions (Yoshida, 1990). Similar qualitative behavior was observed from the 2D solution of Eq. (3.27).

FIG. 2.
Estimates for a circular tokamak plasma with I p = I RE = 6 MA, R 0 = 6.2 m, and a = 2 m. (a) Minor-radial variation of the shift Δ in ψ, shift Δ RE in A and the relative shift δ = Δ RE − Δ for the case when the RE energy is 50 MeV and ν = 2 (or l i = 1.23); (b) major-radial shift of the magnetic axis ( Δ S) vs the normalized internal inductance (li) at different values of RE energy. The corresponding Shafranov shift for a force-free equilibrium is Δ FF. (c) Shift of magnetic axis ( Δ S) vs li and RE energy; (d) average relative shift δ avg vs li and RE energy.
View largeDownload slide

Estimates for a circular tokamak plasma with I p = I RE = 6 MA, R 0 = 6.2 m, and a = 2 m. (a) Minor-radial variation of the shift Δ in ψ, shift Δ RE in A and the relative shift δ = Δ RE − Δ for the case when the RE energy is 50 MeV and ν = 2 (or l i = 1.23⁠); (b) major-radial shift of the magnetic axis (⁠ Δ S⁠) vs the normalized internal inductance (li) at different values of RE energy. The corresponding Shafranov shift for a force-free equilibrium is Δ FF⁠. (c) Shift of magnetic axis (⁠ Δ S⁠) vs li and RE energy; (d) average relative shift δ avg vs li and RE energy.

FIG. 2.
Estimates for a circular tokamak plasma with I p = I RE = 6 MA, R 0 = 6.2 m, and a = 2 m. (a) Minor-radial variation of the shift Δ in ψ, shift Δ RE in A and the relative shift δ = Δ RE − Δ for the case when the RE energy is 50 MeV and ν = 2 (or l i = 1.23); (b) major-radial shift of the magnetic axis ( Δ S) vs the normalized internal inductance (li) at different values of RE energy. The corresponding Shafranov shift for a force-free equilibrium is Δ FF. (c) Shift of magnetic axis ( Δ S) vs li and RE energy; (d) average relative shift δ avg vs li and RE energy.
View largeDownload slide

Estimates for a circular tokamak plasma with I p = I RE = 6 MA, R 0 = 6.2 m, and a = 2 m. (a) Minor-radial variation of the shift Δ in ψ, shift Δ RE in A and the relative shift δ = Δ RE − Δ for the case when the RE energy is 50 MeV and ν = 2 (or l i = 1.23⁠); (b) major-radial shift of the magnetic axis (⁠ Δ S⁠) vs the normalized internal inductance (li) at different values of RE energy. The corresponding Shafranov shift for a force-free equilibrium is Δ FF⁠. (c) Shift of magnetic axis (⁠ Δ S⁠) vs li and RE energy; (d) average relative shift δ avg vs li and RE energy.

Close modal

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 ϕ 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 ∼ 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.

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.

The authors have no conflicts to disclose.

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).

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

The steady state limit of Eq. (2.7), which is
γ m n r v ∥ 2 κ ̂ = − e n r E ⊥ + j r × B ,
(A1)
can also be interpreted in the following manner. For this, it is useful to make the following ansatz for the RE fluid velocity:
v r = v ∥ b ̂ − γ m v ∥ 2 e B ( b ̂ × κ ̂ ) + E × B B 2 .
(A2)
If we now note that ( b ̂ × κ ̂ ) × B = ( b ̂ · B ) κ ̂ − ( κ ̂ · B ) b ̂ = B κ ̂⁠, along with the fact that ( E × B ) × B = ( E · B ) B − ( B · B ) E = E ∥ B B − B 2 E = − B 2 ( E − E ∥ b ̂ ) = − B 2 E ⊥⁠, it can be easily seen that
j r × B = − e n r [ v ∥ b ̂ − γ m v ∥ 2 e B ( b ̂ × κ ̂ ) + E × B B 2 ] × B , = γ m n r v ∥ 2 κ ̂ + e n r E ⊥ .
The above relation in combination with Eq. (A1) shows that the J r × B force arising from the curvature drift provides the centripetal acceleration necessary for the fast toroidal motion. This is in line with the conclusion from the guiding center formalism.
Taking the reduced momentum balance for the REs
v | | 2 κ = − e γ m v r × B ,
(B1)
assuming that κ ≈ − 1 R e ̂ R and just retaining the non-toroidal components from the RHS gives
v | | 2 R e ̂ R = e γ m R n r [ − B ϕ ∇ u − J ϕ , r e ∇ ψ ] .
(B2)
Now again approximating v | | ≈ − J ϕ , r e n r gives us
γ m J ϕ , r 2 e 2 n r e ̂ R = e [ − B ϕ ∇ u − J ϕ , r e ∇ ψ ] .
Doing a cross product of the above equation with ∇ u and multiplying by e give
( γ m J ϕ , r e n r e ̂ R + e ∇ ψ ) × ∇ u = 0 , ( γ m J ϕ , r e n r ∇ R + ∇ ( e ψ ) ) × ∇ u = 0.
Now since γ m and v | | [and, hence, J ϕ , r / ( e n r )] are pure functions of u, the gradient of γ m v ϕ / n r will be parallel to ∇ u⁠. So, it can be taken into the gradient operator along with R. This gives
[ ∇ ( γ m R J ϕ , r e n r ) + e ψ ] × ∇ u = 0 ,
And, hence,
γ m R ( − J ϕ , r e n r ) − e ψ = A ( u ) .
(B3)
The assumptions κ ≈ − 1 R e ̂ R and v | | ≈ − J ϕ , r e n r were essential in proving the above.
1.
Abdullaev
,
S. S.
,
Finken
,
K. H.
,
Wongrach
,
K.
, and
Willi
,
O.
, “
Theoretical and experimental studies of runaway electrons in the TEXTOR tokamak
,” Technical Report No. FZJ-2016-03329 (
Forschungszentrum Jülich
,
2016
).
Google Scholar
 
2.
Bandaru
,
V.
,
Hoelzl
,
M.
,
Artola
,
F. J.
,
Papp
,
G.
, and
Huijsmans
,
G. T. A.
, “
Simulating the nonlinear interaction of relativistic electrons and tokamak plasma instabilities: Implementation and validation of a fluid model
,”
Phys. Rev. E
99
,
063317
(
2019
).
https://doi.org/10.1103/PhysRevE.99.063317
Google Scholar
Crossref
Search ADS
PubMed
 
3.
Bandaru
,
V.
,
Hoelzl
,
M.
,
Reux
,
C.
,
Ficker
,
O.
,
Silburn
,
S.
,
Lehnen
,
M.
,
Eidietis
,
N.
,
JOREK Team, and JET Contributors
. “
Magnetohydrodynamic simulations of runaway electron beam termination in JET
,”
Plasma Phys. Controlled Fusion
63
,
035024
(
2021
).
https://doi.org/10.1088/1361-6587/abdbcf
Google Scholar
Crossref
Search ADS
 
4.
Breizman
,
B. N.
,
Aleynikov
,
P.
,
Hollmann
,
E. M.
, and
Lehnen
,
M.
, “
Physics of runaway electrons in tokamaks
,”
Nucl. Fusion
59
(
8
),
083001
(
2019
).
https://doi.org/10.1088/1741-4326/ab1822
Google Scholar
Crossref
Search ADS
 
5.
Fitzgerald
,
M.
,
Appel
,
L. C.
, and
Hole
,
M. J.
, “
EFIT tokamak equilibria with toroidal flow and anisotropic pressure using the two-temperature guiding-centre plasma
,”
Nucl. Fusion
53
(
11
),
113040
(
2013
).
https://doi.org/10.1088/0029-5515/53/11/113040
Google Scholar
Crossref
Search ADS
 
6.
Freidberg
,
J.
,
Plasma Physics and Fusion Energy
(
Cambridge University Press
,
2007
).
Google Scholar
Crossref
Search ADS
 
7.
Gorbunov
,
E. P.
 and
Razumova
,
K. A.
, “
Effect of a strong magnetic field on the magnetohydrodynamic stability of a plasma and the confinement of charged particles in the tokamak machine
,”
J. Nucl. Energy, Part C
6
(
5
),
515
–
525
(
1964
).
https://doi.org/10.1088/0368-3281/6/5/310
Google Scholar
Crossref
Search ADS
 
8.
Hender
,
T. C.
,
Wesley
,
J. C.
,
Bialek
,
J.
,
Bondeson
,
A.
,
Boozer
,
A. H.
,
Buttery
,
R. J.
,
Garofalo
,
A.
,
Goodman
,
T. P.
,
Granetz
,
R. S.
,
Gribov
,
Y.
,
Gruber
,
O.
,
Gryaznevich
,
M.
,
Giruzzi
,
G.
,
Günter
,
S.
,
Hayashi
,
N.
,
Helander
,
P.
,
Hegna
,
C. C.
,
Howell
,
D. F.
,
Humphreys
,
D. A.
,
Huysmans
,
G. T. A.
,
Hyatt
,
A. W.
,
Isayama
,
A.
,
Jardin
,
S. C.
,
Kawano
,
Y.
,
Kellman
,
A.
,
Kessel
,
C.
,
Koslowski
,
H. R.
,
Haye
,
R. L.
,
Lazzaro
,
E.
,
Liu
,
Y. Q.
,
Lukash
,
V.
,
Manickam
,
J.
,
Medvedev
,
S.
,
Mertens
,
V.
,
Mirnov
,
S. V.
,
Nakamura
,
Y.
,
Navratil
,
G.
,
Okabayashi
,
M.
,
Ozeki
,
T.
,
Paccagnella
,
R.
,
Pautasso
,
G.
,
Porcelli
,
F.
,
Pustovitov
,
V. D.
,
Riccardo
,
V.
,
Sato
,
M.
,
Sauter
,
O.
,
Schaffer
,
M. J.
,
Shimada
,
M.
,
Sonato
,
P.
,
Strait
,
E. J.
,
Sugihara
,
M.
,
Takechi
,
M.
,
Turnbull
,
A. D.
,
Westerhof
,
E.
,
Whyte
,
D. G.
,
Yoshino
,
R.
,
Zohm
,
H
, and
the ITPA MHD, Disruption and Magnetic Control Topical Group,
“
Chapter 3: MHD stability, operational limits and disruptions
,”
Nucl. Fusion
47
(
6
),
S128
(
2007
).
https://doi.org/10.1088/0029-5515/47/6/S03
Google Scholar
Crossref
Search ADS
 
9.
Hoelzl
,
M.
,
Huijsmans
,
G. T. A.
,
Pamela
,
S. J. P.
,
Becoulet
,
M.
,
Nardon
,
E.
,
Artola
,
F. J.
,
Nkonga
,
B.
, et al “
The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas
,”
Nucl. Fusion
61
(
6
),
065001
(
2021
).
https://doi.org/10.1088/1741-4326/abf99f
Google Scholar
Crossref
Search ADS
 
10.
Huysmans
,
G. T. A.
 and
Czarny
,
O.
, “
MHD stability in X-point geometry: Simulation of ELMs
,”
Nucl. Fusion
47
(
7
),
659
(
2007
).
https://doi.org/10.1088/0029-5515/47/7/016
Google Scholar
Crossref
Search ADS
 
11.
Jayakumar
,
R.
,
Fleischmann
,
H. H.
, and
Zweben
,
S. J.
, “
Collisional avalanche exponentiation of runaway electrons in electrified plasmas
,”
Phys. Lett. A
172
(
6
),
447
–
451
(
1993
).
https://doi.org/10.1016/0375-9601(93)90237-T
Google Scholar
Crossref
Search ADS
 
12.
Lackner
,
K.
, “
Deformation of a magnetic dipole field by trapped particles
,”
J. Geophys. Res.
75
(
16
),
3180
–
3192
, https://doi.org/10.1029/JA075i016p03180 (
1970
).
Google Scholar
Crossref
Search ADS
 
13.
Lao
,
L. L.
,
Kruger
,
S.
,
Akcay
,
C.
,
Balaprakash
,
P.
, et al “
Application of machine learning and artificial intelligence to extend EFIT equilibrium reconstruction
,”
Plasma Phys. Controlled Fusion
64
(
7
),
074001
(
2022
).
https://doi.org/10.1088/1361-6587/ac6fff
Google Scholar
Crossref
Search ADS
 
14.
Marder
,
B.
 and
Weitzner
,
H.
, “
A bifurcation problem in E-layer equilibria
,”
Plasma Phys.
12
,
435
–
445
(
1970
).
https://doi.org/10.1088/0032-1028/12/6/003
Google Scholar
Crossref
Search ADS
 
15.
Martín-Solís
,
J. R.
,
Loarte
,
A.
, and
Lehnen
,
M.
, “
Formation and termination of runaway beams in ITER disruptions
,”
Nucl. Fusion
57
(
6
),
066025
(
2017
).
https://doi.org/10.1088/1741-4326/aa6939
Google Scholar
Crossref
Search ADS
 
16.
Matthews
,
G. F.
,
Bazylev
,
B.
,
Baron-Wiechec
,
A.
,
Coenen
,
J.
,
Heinola
,
K.
,
Kiptily
,
V.
,
Maier
,
H.
,
Reux
,
C.
,
Riccardo
,
V.
,
Rimini
,
F.
,
Sergienko
,
G.
,
Thompson
,
V.
, and
Widdowson
,
A.
, “
Melt damage to the JET ITER-like wall and divertor
,”
Phys. Scr.
T167
,
014070
(
2016
).
https://doi.org/10.1088/0031-8949/T167/1/014070
Google Scholar
Crossref
Search ADS
 
17.
Munaretto
,
S.
,
Chapman
,
B. E.
,
Cornille
,
B. S.
,
DuBois
,
A. M.
,
McCollam
,
K. J.
,
Sovinec
,
C. R.
,
Almagri
,
A. F.
, and
Goetz
,
J. A.
, “
Generation and suppression of runaway electrons in MST tokamak plasmas
,”
Nucl. Fusion
60
,
046024
(
2020
).
https://doi.org/10.1088/1741-4326/ab73c0
Google Scholar
Crossref
Search ADS
 
18.
Nygren
,
R.
,
Lutz
,
T.
,
Walsh
,
D.
,
Martin
,
G.
,
Chatelier
,
M.
,
Loarer
,
T.
, and
Guilhem
,
D.
, “
Runaway electron damage to the Tore Supra Phase III outboard pump limiter
,”
J. Nucl. Mater.
241–243
,
522
–
527
(
1997
).
https://doi.org/10.1016/S0022-3115(96)00557-0
Google Scholar
Crossref
Search ADS
 
19.
Paz-Soldan
,
C.
,
Eidietis
,
N. W.
,
Liu
,
Y. Q.
,
Shiraki
,
D.
,
Boozer
,
A. H.
,
Hollmann
,
E. M.
,
Kim
,
C. C.
, and
Lvovskiy
,
A.
, “
Kink instabilities of the post-disruption runaway electron beam at low safety factor
,”
Plasma Phys. Controlled Fusion
61
(
5
),
054001
(
2019
).
https://doi.org/10.1088/1361-6587/aafd15
Google Scholar
Crossref
Search ADS
 
20.
Paz-Soldan
,
C.
,
Liu
,
Y.
,
Reux
,
C.
,
Del-Castillo-Negrete
,
D.
, et al “
A novel path to runaway electron mitigation via deuterium injection and current-driven kink instability
,”
Nucl. Fusion
61
(
11
),
116058
(
2021
).
https://doi.org/10.1088/1741-4326/ac2a69
Google Scholar
Crossref
Search ADS
 
21.
Popovic
,
Z.
,
Hollman
,
E. M.
,
del Castillo-Negrete
,
D.
, et al “
Polarized imaging of visible synchrotron emission from runaway electron plateaus in DIII-D
,”
Phys. Plasmas
28
,
082510
(
2021
).
https://doi.org/10.1063/5.0058927
Google Scholar
Crossref
Search ADS
 
22.
Reux
,
C.
,
Paz-Soldan
,
C.
,
Aleynikov
,
P.
,
Bandaru
,
V.
,
Ficker
,
O.
,
Silburn
,
S.
,
Hoelzl
,
M.
,
Jachmich
,
S.
,
Eidietis
,
N.
,
Lehnen
,
M.
,
Sridhar
,
S.
, and
JET Contributors
. “
Runaway electron beam suppression using impurity flushing and large magnetohydrodynamic instabilities
,”
Phys. Rev. Lett.
126
,
175001
(
2021
).
https://doi.org/10.1103/PhysRevLett.126.175001
Google Scholar
Crossref
Search ADS
PubMed
 
23.
Reux
,
C.
,
Plyusnin
,
V.
,
Alper
,
B.
,
Alves
,
D.
,
Bazylev
,
B.
,
Belonohy
,
E.
,
Boboc
,
A.
,
Brezinsek
,
S.
,
Coffey
,
I.
,
Decker
,
J.
,
Drewelow
,
P.
,
Devaux
,
S.
,
de Vries
,
P. C.
,
Fil
,
A.
,
Gerasimov
,
S.
,
Giacomelli
,
L.
,
Jachmich
,
S.
,
Khilkevitch
,
E. M.
,
Kiptily
,
V.
,
Koslowski
,
R.
,
Kruezi
,
U.
,
Lehnen
,
M.
,
Lupelli
,
I.
,
Lomas
,
P. J.
,
Manzanares
,
A.
,
Martin De Aguilera
,
A.
,
Matthews
,
G. F.
,
Mlynář
,
J.
,
Nardon
,
E.
,
Nilsson
,
E.
,
Perez von Thun
,
C.
,
Riccardo
,
V.
,
Saint-Laurent
,
F.
,
Shevelev
,
A. E.
,
Sips
,
G.
, and
Sozzi
,
C.
, “
Runaway electron beam generation and mitigation during disruptions at JET-ILW
,”
Nucl. Fusion
55
(
9
),
093013
(
2015
).
https://doi.org/10.1088/0029-5515/55/9/093013
Google Scholar
Crossref
Search ADS
 
24.
Tinguely
,
R. A.
,
Hoppe
,
M.
,
Granetz
,
R. S.
,
Mumgaard
,
R. T.
, and
Scott
,
S.
, “
Experimental and synthetic measurements of polarized synchrotron emission from runaway electrons in alcator c-mod
,”
Nucl. Fusion
59
(
9
),
096029
(
2019
).
https://doi.org/10.1088/1741-4326/ab2d1d
Google Scholar
Crossref
Search ADS
 
25.
Vallhagen
,
O.
,
Embreus
,
O.
,
Pusztai
,
I.
,
Hesslow
,
L.
, and
Fülöp
,
T.
, “
Runaway dynamics in the DT phase of ITER operations in the presence of massive material injection
,”
J. Plasma Phys.
86
(
4
),
475860401
(
2020
).
https://doi.org/10.1017/S0022377820000859
Google Scholar
Crossref
Search ADS
 
26.
Waddell
,
B. V.
,
Carreras
,
B.
,
Hicks
,
H. R.
, and
Holmes
,
J. A.
, “
Nonlinear interaction of tearing modes in highly resistive tokamaks
,”
Phys. Fluids
22
(
5
),
896
–
910
(
1979
).
https://doi.org/10.1063/1.862685
Google Scholar
Crossref
Search ADS
 
27.
Wesson
,
J.
,
Tokamaks
,
3rd ed.
(
Clarendon Press
,
Oxford
,
2004
).
Google Scholar
 
28.
Yoshida
,
Z.
, “
A self-consistent equilibrium model of plasma-beam systems
,”
Phys. Fluids B
1
,
1702
(
1989
).
https://doi.org/10.1063/1.858949
Google Scholar
Crossref
Search ADS
 
29.
Yoshida
,
Z.
, “
Numerical analysis of runaway tokamak equilibrium
,”
Nucl. Fusion
30
,
317
(
1990
).
https://doi.org/10.1088/0029-5515/30/2/010
Google Scholar
Crossref
Search ADS
 
30.
Zhao
,
C.
,
Liu
,
C.
,
Jardin
,
S. C.
, and
Ferraro
,
N. M.
, “
Simulation of MHD instabilities with fluid runaway electron model in M3D-C1
,”
Nucl. Fusion
60
,
126017
(
2020
).
https://doi.org/10.1088/1741-4326/ab96f4
Google Scholar
Crossref
Search ADS
 
© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).