Drift surface solver for runaway electron current dominant equilibria during the current quench

Lu Yuan(袁露) and Di Hu(胡地)

School of Physics,Beihang University,Beijing 100191,China

Keywords: tokamak,disruption,runaway electrons,equilibrium

1.Introduction

Plasma disruptions are the result of non-linearly unstable magneto-hydrodynamic(MHD)instabilities which destroy the global plasma confinement and result in rapid and localized energy deposition upon the plasma facing components(PFCs),potentially causing irreversible damage to future high performance devices.[1]Among the threats posed by disruptions,the runaway electrons and their localized deposition onto the first wall is one of the most feared consequence.[1–3]In the most disastrous scenario, significant portion of the original plasma current could be converted into runaway electron current,and all of their kinetic energy as well as a large portion of the magnetic energy it carries could be unleashed locally.[4]This should be avoided at all costs, and there have a lot of effort going on regarding the control and depletion of such runaway current.[5–9]

To address the confinement and loss of runaway electrons in the configuration space, however, one should have a firm grasp on the equilibrium solution,upon which instability studies could be carried out.The runaway electrons, on the other hand, will not exactly follow the magnetic field lines due to their high momentum.Once they become the main current carrier,the current density stops following the magnetic surface,and the current density(or it multiplied by the major radius in toroidal geometry) is no longer a flux function.Instead, the current will follow the drift surface of the runaway electrons,which is essentially the constant canonical angular momentum surface in the axisymmetric geometry.[10,11]Thus the conventional Grad–Shafranov equation[12]for solving plasma equilibrium needs to be revised to describe this new class of equilibria, and a new multi-fluid-like equilibrium equation[13–15]needs to be derived which solves the runaway electron drift surfaces directly.

We are satisfied with the consideration of a simplified runaway current model in this study.We obtain the force balance equation for the strongly passing runaway electron guiding center from its Lagrangian.[10]Such force balance equation is then combined with the toroidal symmetry to yield a Grad–Shafranov-like second order derivative equation for the effective flux (the canonical angular momentum), which can then be solved with given boundary condition.For simplicity,we have assumed that all runaway electrons have the same parallel momentum,although they can have a spatial distribution in density.We then solve the equation numerically assuming static magnetic flux boundary condition for a range of geometric parameters, runaway electron currents, runaway electron momentum and poloidal magnetic fields.Apparent deviation could be found between the calculated drift surfaces and the flux surfaces,depending on the runaway electron momentum,device geometry and the poloidal field strength.It is found that the current center will exhibit a horizontal displacement as the runaway electron momentum changes,even with constant total current and poloidal field,in consistent with results from previous simpler model.[16]The flattened current profile is more susceptible to such displacement with the same momentum change.Furthermore, in the presence of up–down asymmetry in the poloidal field, the runaway current exhibits vertical displacement accompanied by said horizontal displacement.Such displacement is purely caused by their parallel momentum change,and may have important implications for the stability of the current quench plasmas where runaway current is dominant,as well as for the runaway electron scrapping-off as a result of said displacement.

The rest of this paper is arranged as follows.In Section 2,we derive the new Grad–Shafranov-like equation for the runaway electrons from their guiding center Lagrangian.In Section 3, we perform a numerical parameter scan of numerical solutions to investigate the equilibrium evolution under different geometric parameters, runaway electron current, parallel momentum and poloidal field.Finally,the discussion and conclusion are shown in Section 4.

2.The equilibrium equation for runaway electron current

The Grad–Shafranov-like equation could be obtained from two ingredients: the force balance equation and a symmetric direction.For the latter,we are considering a simple 2D equilibrium in this paper,and the symmetric direction is only the toroidal direction.The force balance equation for the cold current quench(CQ)plasma should be simply asJ×B=0 in the limit that the current carrier does not have strong parallel momentum.However, for runaway electrons, the centrifugal force as a result of their high momentum must be considered,and this is balanced by theJ×Bterm.

A simple way of obtaining such force balance equation is to look at the description of the guide center trajectory from its Lagrangian,[10]and realize that such trajectory must automatically satisfy the force balance for the individual particle.In principle,one should then construct a series of moment equations for the runaway electron Boltzmann equation to render the runaway electron fluid equations,[17]from which the force balance equation is obtained.However,we will only consider a simplified case here as a first step that all runaway electrons share the same parallel momentum(δfunction distribution in momentum space for a given configuration space coordinate).In such case, the runaway current density could be expressed as

HerenREis the runaway electron density,qis its charge anduis the same with the individual guiding center velocity.

The Hamiltonian equations of motion for each coordinate of the runaway electron guiding center are[10,11]

whereXis the configuration space coordinates,p‖is the parallel momentum,μis the magnetic moment andθis the gyroangle,bis the magnetic field direction,γis the relativistic factor,qandmare the electron charge and rest mass,respectively.We have also defined the following effective fields:

whereAis the vector potential andΦis the electric potential.It should be noted that,if we want to consider higher parallel momentum scenario wherep‖/(qBR)∝εwithε ≡a/Rbeing the inverse aspect ratio, as is assumed in Ref.[11], then additional curvature drift terms need to be added in the above equations.In this study,however,we considerp‖/(qBR)∝ε2,so that the higher order curvature terms in Ref.[11]are considered as next order effects in this paper and are neglected.This ordering, although less ambitious, should already cover the major parameter range of interest during the current quench for strong magnetic field device such as the International Thermonuclear Experimental Reactor (ITER), since the radiative and collisional drag,as well as kinetic instabilities would tend to limit the upper bound of the runaway electrons’ parallel momentum.For strongly passing runaway electrons in quasineutral plasmas, we approximate ?Φ??0.In steady states,all partial time derivative terms also vanish.So that we have

The parallel momentum is thus a constant along the trajectory line,thus a function of the drift surface,and the runaway electron trajectory is determined by the effective magnetic fieldB?.Substituting into Eq.(1),we obtain

Note that Eq.(12)automatically satisfies the force balance in the cases we consider here, and the force balance equation simply yields

This is analogous to a force free equilibrium in an ordinary plasma[18]where we have

Substituting Eqs.(6) and (7) into (13) and compare with Eq.(14), the only difference is the additional termJ×.We will see later on that this term will enter the right hand side of our new equilibrium equation.In the axisymmetric system we considered here, the magnetic fieldBand effective magnetic fieldB?could be expressed as the reduced form

HereΨis the poloidal magnetic flux, the effective fluxΨ?is related to the toroidal component of the effective vector potentialA?which in turn is related to the canonical angular momentum of the runaway electrons.Meanwhile ˉFand ˉF?are related to the toroidal field and the effective toroidal field,respectively.The constantΨ?surfaces then define the drift surfaces of the runaway electrons, which are the new characteristic surfaces analogous to the flux surfaces in an ordinary plasma equilibrium.Combining Eqs.(12)and(16),one could find

Hence the runaway electron trajectory and thus the current density flow along the drift surfaces.Note that this meansnREcannot be a flux function ofΨ?in toroidal geometry due to the demand of current conservation, althoughnRERis.This must be compensated by the cold bulk electrons which have much higher density,so that the quasi-neutrality could be maintained.

With Eqs.(17)and(18),it is evident that the runaway current dominant equilibrium can be defined by the drift surfaces,just like the flux surfaces define the ordinary plasma equilibrium.We then proceed to solveΨ?.In our toroidal symmetric system,realizing that

Hence we get a new(cold)Grad–Shafranov-like equation for the runway electron current

Note that here we have used the fact that the toroidal magnetic field ˉFis a function of the effective fluxΨ?.Unfortunately,this equation could not be directly solved since the RHS is not a function ofΨalone, but are related toΨ?.To solve it, one could try to convert this equation into a second order differential equation forΨ?instead,and make sure that the RHS only contains the function ofΨ?and some spatial coordinates such asR.

To do this,we realize

We note that Eqs.(21)and(22)should be equivalent so that

which is a function ofΨ?.However,F?is not necessarily a function ofΨ?,and it is desirable to rewriteF?to separate the part which is a function ofΨ?and the part which is not.

To do this,we rewrite Eq.(16)as

whereζandrepresent the additional contribution from the mechanical momentum.From Eqs.(15) and (24) we must have

Realizing that under our current assumptions we have Eq.(11),so thatp‖is a constant on a given drift surface.In this paper,as a first step, we are satisfied with consideration of the simplest case thatp‖is also a constant over the whole volume.We can then further simplify and obtain

We could obtain the expression of ?ζin similar manner from Eq.(24)

Thus we recover the relationship between the effective fluxΨ?and the magnetic fluxΨ

which is consistent with Eq.(6).Note here the contour ofΨis our flux surface solution while the contour ofΨ?is the drift surface solution, hence for finitep‖there will be a deviation between the flux surfaces and the drift surfaces.Substituting Eqs.(30)and(32)into Eq.(22),we immediately get the relationship between??Ψ?and??Ψaccording to Eq.(32)

This then yields

Using the following approximation:

we could rewrite

So that we finally have

Compared with ordinary Grad–Shafranov equation and taking into account of Eq.(33),it can be seen that both our drift surface solutionΨ?and flux surface solutionΨwould differ from the flux surface solution of the ordinary Grad–Shafranov equation for finitep‖,while both of the first two solutions converge to the last in the limit ofp‖ →0.With given current density profile,p‖and boundary condition, the RHS of Eq.(40) is known and we could proceed to solve the equilibrium.

3.The numerical solution of the runaway electron current dominated equilibrium

In this section,we present the numerical investigation of the runaway electron current under different control variables such as the parallel momentum,the current density profile and the externally applied vertical magnetic field.Specifically, in Subsection 3.1,we introduce our numerical setup for the runaway electron current equilibrium solver.In Subsection 3.2,the runaway current equilibrium with various parameters are calculated and plotted for both the ITER-like plasma and a smaller MAST-like plasma with complete up–down symmetric boundary condition, and found the change in the parallel momentum would result in horizontal current center displacement even with all the others parameters intact.In Subsection 3.3,we demonstrate that such horizontal displacement is accompanied by a vertical displacement if the boundary condition exhibits some up–down asymmetry.

3.1.Numerical setup

Our numerical solver bases on the open source Grad–Shafranov solver for ordinary plasma equilibria FreeGS.[19]The equilibrium equation is changed to that of our runaway current equilibrium in Section 2, and the boundary condition is modified so that we have an ideally conducting boundary condition for the magnetic flux,plus whatever static externally applied poloidal magnetic field.Explicitly,the boundary condition for the magnetic flux in the up–down symmetric case is

This corresponds to the toroidal component of the magnetic vector potentialAφ=Bz0R,which represents a constant vertical fieldBz0.The corresponding boundary condition forΨ?is then

This boundary condition represents the scenario that the current displacement within the vacuum vessel is faster than the resistive time of the wall,hence all magnetic flux changes will be screened, although the external field, which penetrates on a much longer timescale,is frozen in the wall.For simplicity,we only consider a rectangular boundary in this study.For the so-called MAST-like geometry, we haveRmin= 0.1 m,Rmax=2.0 m,Zmin=?2.0 m, andZmax=2.0 m, while for the ITER-like geometry, we haveRmin=4 m,Rmax=8.5 m,Zmin=?4 m, andZmax=5.0 m.More detailed studies with realistic 2D first wall shape are left for future works.

In this section, we mainly scan the parallel momentump‖and a constant,uniform vertical external magnetic fieldBz0under different current density profiles.Such uniform vertical field of course is a simplification meant to provide a first estimation for our equilibrium solution.The realistic poloidal field or vertical field generated by the external coils is much more complicated in reality,[20]and we aim to include such realistic external coils in our future works.We choose the monotonic current density profile to have the following shape:

Here ˉψ?is the normalized effective flux.Several different parametersmandnare chosen to represent three kinds of distributions from the peaked profile to the flattened profile.Additionally, we consider one non-monotonic current profile to represent the possible hollowed runaway current spatial distribution[21–23]

All four kinds of current density profile are shown in Fig.1 for comparison, with different colour representing different choice ofmandnas well as the monotonic and non-monotonic profiles.Once the current profile shape and the total runaway current are chosen,the RHS of Eq.(40)is known.

Fig.1.The shape of current density profiles according to Eqs.(44) and(45).The blue,orange and green curves represent three different monotonic current profifiles,the red one is the non-monotonic one.

Apart from the above constant,uniformBz0case,we also investigate an additional case with up–down asymmetry in the magnetic field boundary condition,to mimic the up–down asymmetry in the poloidal field coil in real tokamaks.In this simple study we only consider two poloidal field coils located atRu=5 m,Zu=7.5 m andRd=6 m,Zd=?7.5 m, respectively,in addition to the uniform vertical fieldBz0.In this up–down asymmetric case,the boundary condition becomes

Here

is the total toroidal magnetic vector potential contribution from the external PF coils, whereaiis the radius of coili,riis the distance from the boundary point to the coil centre.Here we assumeri>ai.Further,θiis the poloidal angle of the boundary point related to the coil centre,Iiis the coil current andis the first order associated Legendre polynomial.The corresponding boundary condition forΨ?is

With those boundary conditions, we can proceed to solve the runaway electron current dominant equilibrium.

3.2.Numerical results for up–down symmetric boundary condition

We demonstrate here the runaway current equilibria calculated using Eq.(40)and the up–down symmetric boundary conditions described in Subsection 3.1.

To demonstrate the difference between our new equilibrium solution and the ordinary Grad–Shafranov solution, as well as their convergence in the limit ofp‖ →0, we hereby show two benchmark cases with FreeGS results using the ideal boundary MAST-like geometry withBz0=0 andIp=200 kA.The flux surface and the drift surface of our equilibrium solution as well as the flux surface solution from FreeGS are compared in Fig.2 for both thep‖=0mecand thep‖=20meccases.The black chained lines represent the Grad–Shafranov solution, the dark blue solid lines represent the new drift surface solution and the red dashed lines represent the new flux surface solution.In Fig.2(a) where we havep‖=0, exact agreement is found between all three surfaces as the contours of the three solutions overlap with each other.This demonstrates we could recover the ordinary force free Grad–Shafranov solution in the low electron momentum limit.As the electron momentum increase,however,deviation between the solutions becomes gradually more pronounced, as can be seen in Fig.2(b)where we havep‖=20mec.Here the Grad–Shafranov solution is unchanged, as the parallel momentum of the current carrier does not enter the equilibrium equation explicitly.Meanwhile,both the flux surfaces and the drift surfaces of our new equilibrium solution exhibit an outward displacement as the electron parallel momentum increases.The deviation between the drift surface and the flux surface solution also begins to manifest itself,as this deviation is proportional top‖as we discussed in Section 2.

Fig.2.The comparison between the flux surface (red dashed lines) and the drift surface(blue solid lines)of our new equilibrium solution as well as that of an ordinary Grad–Shafranov equation(black chained lines)with MAST-like geometry,Bz0=0 and Ip=200 kA for(a)γ=0 and(b)γ=20.

Fig.3.The runaway current equilibrium with Bz0=?0.045 T,and(a)γ=56,(b)γ=116,(c)γ=176 and(d)γ=200,respectively.The current profile is monotonic with m=1 and n=6 and total runaway electron current IRE=200 kA.The solid contours represent the drift surfaces and the dashed contours represent the flux surfaces.

We then look at the peaked current profile case described by Eq.(44)withm=1 andn=6 and total runaway electron currentIRE=200 kA but with varyingp‖.The equilibrium solution for several chosenp‖andBz0=?0.045 T as is shown in Fig.3,where the solid contours represent the drift surfaces and the dashed contours represent the flux surfaces.This choice ofBz0value is on the same order of magnitude as one would expect to maintain the horizontal position of the current centre for an ordinary plasma,[24]although it is not exactly the same due to the contribution from the runaway electron parallel momentum and the back reaction from the ideal wall.The orange vertical dashed line marks the horizontal position of the current center.The color bar for the right figure in each subplot represents the value ofJ‖Rwhich is a function ofΨ?.Note that such a centrally peaked current profile is likely to trigger some kink instability for the runaway current, but this is not the concern of our 2D equilibrium study here and we leave detailed stability studies for future works.From Figs.3(a)to 3(d), it can be seen that even with exactly the same vertical field and the same runaway electron current, the current center shows apparent displacement as the parallel momentum changes,in agreement with previous simplified model.[16]Such parallel momentum change could be caused by many mechanism,such as the collision with background species,[25]radiative drag,[26]kinetic instabilities[27]or toroidal electric field.In the smallp‖limit,the deviation between the drift surface and the flux surface is small, while such deviation consistently grows asp‖increases.In the largep‖limit,although there is only very little closed flux surface remains,significant closed drift surfaces still exist, and the runaway electrons remain well confined even in the open field line region.This feature cannot be captured with the ordinary Grad–Shafranov equation since the current flows along the flux surfaces, and there can not be closed current surfaces in open field line region.This runaway electron drift is closely associated with the aforementioned current center displacement,as the drifted runaway electrons are also the main current carrier, thus any runaway electron drift would also move the current channel,which in turn moves the magnetic flux, which causes further runaway electron drift.Thus we see inward current displacement for decreasing runaway electron momentum, and outward displacement for growing momentum.

As a comparison, we show another example of runaway current equilibrium with more flattened current profilem=2 andn=2, as is plotted in Fig.4.Similar horizontal current displacement caused by the change of the parallel momentum alone can be seen.Consistent with them=1 andn=6 case,closed runaway electrons drift surfaces exist where field lines are open.Compared with Fig.3, however, Fig.4 shows significant shrinking of the current channel as the horizontal current displacement occurs.This is due to the shrinking of our last closed drift surface as the runaway electron current moves towards the wall, which mimics the scraping-off of the edge runaway electron current and it redistributes to the core by the inductive electric field in real current displacement event.Comparing the value ofp‖at similar current center position in Figs.3 and 4, the more flattened current profile requires less change inp‖for the same current displacement, suggesting that their current center position is more sensitive to the parallel momentum change.This may be due to the fact that the flattened current profile case experiences more significant current shrinking compared with the peaked current profile case,which results in stronger change in their internal inductance.

We now take a look at the non-monotonic current profile case withm=1 andn=9 described by Eq.(45),the equilibrium solutions for several choices ofp‖are shown in Fig.5.Comparing Figs.5 and 4, this hollowed current profile case shows similar displacement with them=2 andn=2 monotonic current profile shape for the samep‖change.A similar scraping-off of the edge current density could be seen comparing Figs.5(a)and 5(b)with Fig.5(c).For our case studied here, the non-monotonic feature does not seem to have significant contribution to the susceptibility of the current center displacement.

A more straightforward comparison of the current center displacement as a function ofp‖andBz0for the four current profiles shown in Fig.1 is plotted in Fig.6.All cases are with the same constant total currentIRE=200 kA.There the colour of the pseudo-colour map represents the major radius of the current center, with deep blue representing the cases where the equilibrium can not be found due to too significant current center displacement.All figures in Fig.6 show the same trend that increasing the value of the vertical field (which points downward alongzdirection)pushes the current center inward,which is also expected for ordinary plasma equilibrium.Contrary to the ordinary equilibrium,changing the momentum of the runaway electrons without changing the total current or the current profile shape pushes the current center outward.Comparing Figs.6(a)–6(c), it is found that the current center displacement occurs more easily for the more flattened current profile with the same amount ofp‖change andBz0change,consistent with our observation above,and the boundary within which the equilibrium could be found becomes narrower for increasingly flattened current profile.Such trend can be more directly seen in Fig.7, where the different colours represent different current profiles andBz0=?0.045 T.It is apparent that the peaked profile (blue) case takes morep‖change to achieve the same current center displacement compared with the flattened profiles (orange and green).Again,the non-monotonic current profile does not seem to have significant impact of the displacement susceptibility.

Fig.4.The runaway current equilibrium with Bz0=?0.045 T,and(a)γ=56,(b)γ=116,(c)γ=176,and(d)γ=200,respectively.The current profile is monotonic with m=2 and n=2 and total runaway electron current IRE=200 kA.The solid contours represent the drift surfaces and the dashed contours represent the flux surfaces.

Fig.5.The runaway current equilibrium with Bz0 =?0.045 T, and (a) γ =56, (b) γ =116, (c) γ =176, and (d) γ =200, respectively.The current profile is non-monotonic with m=1 and n=9 and total runaway electron current IRE =200 kA.The solid contours represent the drift surfaces and the dashed contours represent the flux surfaces.

Fig.6.The current center major radius position as a function of p‖ and Bz0 for(a)the monotonic m=1,n=6 case,(b)the monotonic m=2,n=2 case, (c) the monotonic m=5, n=1 case and (d) the non-monotonic m=1, n=9 case.The dark blue colour represents where the equilibrium cannot be found due to too much current displacement.

Fig.7.The current centre displacement as a function of p‖ for Bz0 =?0.045 T for different current profile shapes and the same total current IRE=200 kA.

Furthermore, from Eq.(40), it can be deduced that the deviation between the drift surfaces and the flux surfaces is determined by the relative strength of the second term of the RHS compared with that of the first term.One would thus naturally expect that with stronger runaway electron current,thus stronger poloidal magnetic field,the aforementioned deviation would decrease; while for smaller and more compact device,such deviation will be more pronounced.

To investigate this, we first compare two cases with the sameIRE=200 kA and monotonic current profile withm=2 andn=2, one with the same ITER-like geometry we have been using, the other with the MAST-like geomery as is described in Subsection 3.1.The current center displacements as a function of their respectivep‖andBz0are shown in Fig.8,noting the different range of the axes and the color bar.For the same amount ofp‖,the MAST-like case,with its smaller major radius, experiences more significant horizontal drift compared with the ITER-like case.Mathematically, this is due to the 1/Rcontribution in the second term on the RHS of Eq.(40).Physically, this is due to the stronger centrifugal effect in the more compact device.So that we would expect the runaway electron drift behavior we show here to be more important for future small,compact machines.

We then move on to investigate the effect of the runaway electron current with two ITER-like geometry cases,one withIRE=200 kA, the other withIRE=1 MA.Both cases share the same monotonic current profile shape withm= 2 andn=2.Examples of their respective equilibrium, as well as the pseudo-colour plot of the current center displacement are shown in Fig.9.Comparing Figs.9(a) and 9(b), both cases usingp‖=200mecandBz0=?0.045 T, it can be seen that the deviation between the drift surfaces and the flux surfaces is weaker for theIRE=1 MA case and with less outward current center displacement, confirming our previous argument regarding such deviation depending on the relative strength of thep‖/(qR) term.Nevertheless, the runaway electron drift is still on the order of tens of centimeters for this high parallel momentum case,causing the runaway current equilibrium to deviate significantly from that of an ordinary equilibrium.Comparing Figs.9(c) and 9(d), note the different color bar range, we could see that the current center displacement susceptibility is less for theIRE=1 MA case for given amount ofp‖, consistent with our previous argument.Hence we would expect the above discussed runaway electron drift and the associated current displacement to be a significant effect mostly near the end of the runaway current depletion process when the total runaway current is small.

Fig.8.The current center major radius position as a function of p‖ and Bz0 for(a)the monotonic m=2,n=2 case,(b)the monotonic m=2,n=2 case.The two panels(a)and(b)depict the range of o point values under different devices,the panel(a)is taken from Fig.6 above when the o point value range is m=2, n=2 under the ITER device, and the panel (b) depicts the MAST, that is, the o point value range under the small device.All cases are with the same constant total current IRE=200 kA.The dark blue colours represent where the equilibrium cannot be found due to too much current displacement.

Fig.9.(a)The equilibrium solution for IRE=200 kA and(b)the equilibrium solution for IRE=1 MA.Both with m=2,n=2 monotonic current profile,Bz0=?0.045 T and γ=100.The current center major radius positions as a function of p‖and Bz0 for both cases are shown in panels(c)and(d),respectively.Here panel(c)is taken from Fig.6 and corresponds to IRE =200 kA and panel(d)corresponds to IRE =1 MA.Note the difference in the color bar range between panels(c)and(d).

3.3.Displacement of the guide center under asymmetric boundary conditions

Last, let us take a look at the interesting case of an up–down asymmetric boundary condition,as is introduced in Subsection 3.1.We set the upper coil current to beIu=85.94 kA and the lower coil current to beId=77.35 kA.The poloidal magnetic field generated by both coils would be reflected in the boundary condition.The equilibrium solutions under such boundary condition for several selectedp‖andBz0=?0.018 T are shown in Fig.10.Again,the total runaway electron current isIRE=200 kA and the current profile is monotonic withm=2 andn=2.

From Figs.10(a)to 10(d),one could clearly see that asp‖decrease, there is an inward horizontal current displacement just as in the previous cases, but there is also an accompanied vertical displacement due to the asymmetry in the external poloidal field.Since it is rather common to have some up–down asymmetry in the poloidal field coil in real tokamaks,one would expect that the momentum change in real runaway current would also induce a corresponding vertical displacement during the current quench, which should be considered as an additional contribution to the vertical displacement in the ordinary equilibrium solution.Here we are satisfied with merely demonstrating the existence of such vertical displacement as a result of the runaway electron momentum change.The more detailed study with more realistic poloidal coil and 2D wall will be left for future works.

Fig.10.The runaway current equilibrium with Bz0=?0.018 T,and(a)γ =48,(b)γ =36,(c)γ =24 and(d)γ =12,respectively.The current profile is non-monotonic with m=2 and n=2 and total runaway electron current IRE =200 kA.The position of the coil we add in this case is xu =5 m, xd =6 m, yu =7.75 m, and yd =?7.75 m.The magnitude of the current of the coil we add in this case is Iu =85.94 kA and Id=77.35 kA.The solid contours represent the drift surfaces and the dashed contours represent the flux surfaces.The red dashed line represents the separatrix.

4.Discussion and conclusion

A Grad–Shafranov-like equation for the simple case that all current is carried by the runaway electrons and all runaway electrons share the same parallel momentum being derived,solving directly for the runaway electron drift surfaces instead of the magnetic surfaces.A new force balance equation is obtained for the runaway electrons from their guide-center equation of motion, and combined with the toroidal symmetry to yield the aforementioned runaway current equilibrium equation,where an additional contribution from the runaway electron parallel momentum appears on the RHS.This additional contribution would result in the relative drift between the runaway electron drift surface and the flux surface.Since the runaway electrons are the current carriers,such drift naturally results in the current channel displacement and thus the flux surfaces themselves are displaced,which in turn incurs further drift of the runaway electrons.With given runaway current profile shape and given boundary condition,we are able to numerically solve this new equilibrium equation.

Numerical investigations are carried out using simple rectangular geometry with ITER-like and MAST-like geometric parameters.A number of different current profile shapes are tested.In general,it is found that with given external vertical fieldBz0and given total runaway electron currentIRE,the current center would exhibit an inward major radial displacement as the runaway electron momentum decreases,and likewise an outward displacement as their momentum increases.The flattened current profile shows more susceptibility on this displacement compared with the peak profile,probably due to the more significant current scraping-off for the former.The hollowed current profile does not show significant impact on such current displacement susceptibility.As discussed above,this current displacement is caused by the relative drift between the runaway electrons and the flux surfaces,and its effect depends on both the relative strength of the parallel momentum to the poloidal field, and the major radius.Indeed,it is found that in smaller device the runaway electron drift is more apparent and the current center is more easily displaced by changes in the parallel momentum.Meanwhile,if the total current is large and the poloidal field is strong, such drift is milder,so is the current displacement.Thus we conclude that the above effect is more pronounced near the end of the runaway current depletion process when the total current is small,and so is the case in future compact devices.Furthermore, it is found that with up–down asymmetry in the boundary condition,such as those caused by asymmetric poloidal field coil in real tokamaks, the horizontal current displacement is accompanied by a vertical displacement.Hence the change in the runaway electron momentum could also additionally contribute to the vertical displacement event.

Despite the result obtained above, there is still a lot of room for future improvement of this new equilibrium theory.For instance,one could consider the situation where the parallel momentum is no longer a constant for all runaway electrons, but have a distribution in the configuration space instead.This would result in some additional term concerning ?p‖, and is a natural future development of our theory.Moreover,one could also consider the scenario where the runaway electron distribution is not aδfunction in the momentum space,but could have a distribution instead.To obtain the force balance equation in such a scenario,one have to construct the moment equations from the runaway electron kinetic equation to yield the equation of motion for the runaway electron fluid.Then one could combine such force balance equation and the symmetry of the system to obtain the Grad–Shafranov-like equation for the runway electron fluid just like we did in this paper.This is also one of the future works we will pursue.Another possible expansion of current work is to consider a more realistic 2D wall and the associated poloidal field reflecting that of a real tokamak instead of the uniform vertical field we used here.Yet another natural future expansion of the current simplified model is to extend to the higher parallel momentum regime ofp‖/(qBR)∝εinstead ofp‖/(qBR)∝ε2in this paper.The curvature terms in Ref.[11]need to be included for this.

Nevertheless,the runaway current equilibrium theory presented in this study already provides a qualitative estimation on some of the important runaway current dynamic such as their relative drift to the flux surface and the current displacement as a result of parallel momentum change.Based on the equilibrium theory presented here and future works along similar lines, further stability analysis could be carried out for runaway electron dominated current and their consequential runaway electron transport or termination scenarios could be studied.These would help to achieve better runaway current control and depletion in the current quench phase of disruption.

Acknowledgments

The authors thank L.E.Zakharov for fruitful discussion.The authors thank Ben Dudson for the original open source equilibrium solver, upon which our numerical solver for the runaway current equilibrium is developed.Project supported by the National MCF Energy Research and Development Program of China(Grant No.2019YFE03010001)and the National Natural Science Foundation of China(Grant No.11905004).