Experimental Runaway Electron Current Estimation in COMPASS Tokamak

Runaway electrons present a potential threat to the safe operation of future nuclear fusion large facilities based on the tokamak principle (e.g., ITER). The article presents an implementation of runaway electron current estimations at COMPASS tokamak. The method uses a theoretical method developed by Fujita et al., with the difference in using experimental measurements from EFIT and Thomson scattering. The procedure was explained on the COMPASS discharge number 7298, which has a significant runaway electron population. Here, it was found that at least 4 kA of the plasma current is driven by the runaway electrons. Next, the method aws used on the set of plasma discharges with the variable electron plasma density. The difference in the plasma current was explained by runaway electrons, and their current was estimated using the aforementioned method. The experimental results are compared with the theory and simulation. The comparison presented some disagreements, showing the possible direction for the code development. Additional application on runaway electron energy limit is also addressed.


Introduction
In the last decade, nuclear fusion started to shift from science to industry, where the tokamak-based ITER is playing the leading role.Runaway electrons (REs) present a potential threat to the safe operation of the future nuclear fusion power plants based on the tokamak principle.Namely, in the latest issue of ITER Physics Basis [1], REs are considered as the second-highest priority for the ITER disruption mitigation.
COMPASS tokamak [2] is a suitable and low-risk test bed for studying RE disruption physics and RE model benchmarking.However, one of the common missing RE parameters for comparison of experiment and simulation is RE current or RE density.This article addresses this issue for certain cases.
The basic idea behind the used model was reported on ORMAK tokamak by Fujita [3].In contrast to Fujita's full-theoretical approach, we implemented experimental data from EFIT [4] and Thomson Scattering [5] to the calculation.The obtained results are in the acceptable range and their application on density scans and RE localization is presented.
The article starts with a description of an experimental observation of the particular discharge.Theoretical models and detailed description of the settings used in the models for comparison with the experiments are then reported.An estimation of the RE current I RE is the first experimental analysis.Subsequently, the influence of the REs on the current ramp-up phase of the discharge is evaluated.Following that, a principle on how to localize RE beam is done in the next section, using the knowledge of the RE beam current and its maximum energy.Finally, discussions and conclusions are given in the last section with an outlook towards expected future works.

Experimental Setup: COMPASS Tokamak
The COMPASS tokamak is a compact experimental fusion facility with a major radius R 0 = 0.56 m and a minor radius a p = 0.23 m.The toroidal magnetic field B tor is in the 0.9-1.6T range (typically set to 1.15 T), coming from 16 toroidal field coils.The plasma current I p can reach up to 400 kA using an air-cored transformer.The range for the electron densities is flexible and is typically in the 10 19 -10 20 m −3 domain.Plasma shaping varies from circular and elliptical to single-null D-shaped ITER-like plasmas.When circular, the plasma is limited by a carbon HFS wall.The regular pulse length is ∼ 0.4 s, although low current circular discharges with RE can last almost 1 s [6,7].
The Thomson Scattering diagnostic system in COMPASS uses two Nd:YAG lasers of wavelength 1064 nm with an energy of 1.5 J for each laser.The whole system is oriented vertically and the scattered light is recorded in the radial direction.Each laser has 30 Hz repetition rate, which offers a ∼ 16.7 ms time resolution if the two lasers are operated equidistantly.
The EFIT reconstruction provides relevant quantities on a 2D mesh cross-section of the poloidal plasma plane solving the Grad-Shafranov equation for plasma MHD equilibrium, constrained by the magnetic diagnostics at the COMPASS tokamak (current measured in plasma and field coils, pick-up coils).Additionally, in the analysis addressed here, a toroidal loop voltage in plasma V plasma for computing a maximum RE energy W max is calculated using METIS [8], a Multi-Element Tokamak-oriented Integrated Simulator.

Experiment: High RE Current Discharge
The very first magnetic observations of high RE current in the COMPASS tokamak are reported in Ref. [9] for discharge #7298.The discharge had circular-shaped limited plasma, with toroidal magnetic field B tor = 1.15 T and plasma current I p = 120 kA.The plasma density n e was lower than 1.5 × 10 19 m −3 in the plasma core as measured by the Thomson Scattering system.
The COMPASS tokamak is well-supplied with the general plasma measurements, which provided us with the following list of observations in the discharge:

•
Loop voltage measured plasma voltage V plasma bellow 1 V, while the typical range is 1.5-2 V; • Electron temperature from Thomson Scattering was 1 order of magnitude colder than ordinary tokamak plasma; • Hard X-ray scintillators observed a significant amount of RE losses; • Poloidal internal pick-up coils detected inner structure that is moving outwards; • Vertical field for feedback of horizontal plasma position was increasing during the discharge, even though plasma current was constant; • EFIT reconstruction of poloidal cross-section depicted plasma shrinkage towards inner side, due to the increasing magnetic pressure coming from the above-mentioned vertical field; • Plasma pressure in terms of its ratio to the poloidal magnetic field pressure (i.e., poloidal beta β pol ), estimated from EFIT showed unrealistic large values.
All aforementioned information provides us with enough evidence to conclude that RE beam co-existed with the bulk plasma.For more detailed discussion on the measurements, the interested reader is referred to Section 5.1 in Ref. [10].

Method: RE Current Calculation
The overestimate of the β pol is exploited here for the estimate of the RE current I RE .The analysis is based on work done by Fujita et al. [3].Even though the basic principle used here is the same as Fujita's, the final approach to the calculation is different.While the original paper estimates all necessary parameters theoretically, we use direct measurements from the Thomson scattering and calculations from EFIT reconstruction.In this section, we only address our changes to the method.The main idea behind the Fujita's model is an estimate of RE pressure P RE from the difference between total plasma pressure p tot (coming from EFIT β pol ) and pure bulk plasma pressure p pl (coming from the Thomson Scattering measured radial T e profile) where r is the radial position along the minor radius a p .Examples of p tot V and p pl V are presented in Figure 1.The P RE = p tot V − p pl V is calculated as an average from obtained P RE values corresponding to results when following conditions p tot V > 1.4 p pl V and p tot V > 2 p pl V are satisfied simultaneously, to reduce the propagation of error given by this approach.The first condition specifies the accuracy of Thomson Scattering measurements and thereby defines the minimal threshold necessary to distinguish between p tot V and p pl V and eliminates the error of Thomson Scattering measurements.The second condition is always valid in the case of non-negligible RE pressure and it is artificially taken as a conservative limit.Furthermore, P RE is connected to the RE density where represents average in both spatial and momentum coordinates, γ is relativistic Lorentz factor, and v and v ⊥ denote the RE velocity in the parallel and perpendicular direction to the magnetic field B, respectively.By estimating the RE energy, one can estimate the RE density n RE , from which RE current I RE can be approximately derived using the following expression: where A RE is the area of the RE beam cross-section and β is normalized relativistic velocity.

RE Energy Calculation
The model for maximum obtainable RE energy W max is based on Ref. [11], where radiation losses are subtracted from the acceleration.It is a 0D model, as there is no dimensional dependency, but only time evolution.Particularly, the lost power by synchrotron radiation P synch and bremsstrahlung radiation P brems are subtracted from the power gained by the RE thanks to the electric field P E dW max dt = P E − P synch (5) to obtain the time dependence of W max .The power gain by the toroidal electric field E tor is simply calculated as with R p standing for the plasma major radius.The synchrotron power loss where r e and m e are the classical electron radius and mass, respectively.The curvature radius R c is calculated using Equation ( 9) from Ref. [12]: where 2/π at the end of the denominator comes from the average taken over the gyro-motion angle, and where η = v ⊥ /v dr and v dr is the drift velocity with Ω being the fundamental gyro-frequency.Finally, the angle brackets are making the whole calculation dependent on the RE distribution function (REDF).

RE Distribution Functions
With mentioning REDF, we need to address its influence on the I RE calculation.The REDFs used by Fujita et al. and utilized here are monoenergetic f mono , uniform f uni , and linear f lin .The definitions of the first three distributions as functions of RE energy w are as follows: with δ() being the Dirac delta function.The coefficients come from setting the REDF maximum value to 0 at W max and then normalizing the REDF to 1. Furthermore, the exponential REDF is introduced to present a distribution with a steep decrease that could be important if the avalanche mechanism were to dominate.Additionally, two more REDFs are considered here: skewed Gaussian f sG [13] with s = (w − l)/ζ and Maxwell-Jüttner f MJ [14] f MJ (w; distributions.The (negatively) skewed Gaussian represents a more realistic case of f mono , as it is not a delta-function.In other words, the negatively skewed Gaussian represents REDF where almost all RE tend towards W max but they never reach it.The Maxwell-Jüttner distribution function f MJ is selected as one of the most commonly relativistic distribution functions used in the literature.As all three latter functions have asymptotic behavior towards zero and/or W max , normalization has to be different than for the former three.Their parameters are set in such a way that asymptotic edge is smaller than 10 −5 , while f sG maximum is located just over w = 0.95W max .

Pitch Angle
Equation (3) indicates that the pressure separation through the velocities v and v ⊥ implies knowing the pitch angle θ.Fujita et al. omitted the γv 2  ⊥ term and approximated the other one γv 2 by c 2 γβ 2 .However, we consider the θ influence to be consistent with the estimation of W max reported in Section 4.1, where θ plays a role through the parameter η ∝ β ⊥ / γ β 2 .
Table 1 shows the influence of different REDFs on essential parameters for I RE calculation.Firstly, it is important to state that averages presented in Table 1 are done for energy span from 0 to W max using the mean of the variable, i.e., neglecting the RE limit at lower energies.The error is negligible due to the high W max reached in the COMPASS discharges.The average energy w is calculated for illustration, while β can be seen in Equation (4).From Equation (4), one would expect that the exponential REDF has the smallest I RE for the same P RE , as it has the far smallest β .However, due to the term β 2 γ + 0.5 β 2 ⊥ γ in Equation ( 3) used for the n RE estimation (by knowing P RE in Equation ( 3)), the result is opposite.Namely, difference in 1 of an order of magnitude creates much larger difference in estimated n RE between REDFs than β does in Equation (4).Therefore, even though the exponential has by far the smallest β , it has the largest I RE .
Furthermore, even though two extreme values are taken for θ to investigate maximum possible effect from the corresponding variations, the difference due to β 2 γ is smaller than a few percent for every REDF.Furthermore, the aforementioned η is calculated for θ = {0.1,0.3} rad, as it is zero at θ = 0 rad.Here, a significant influence of the pitch angle on the η is observed.Henceforth, θ is implemented in η that is relevant for W max calculation.Another expected dependence that can be noted from Table 1 is that the most similar REDF to the monoenergetic one is the skewed Gaussian f sG REDF, as their parameters are close to each other.

Result: COMPASS Discharge 7298
The pitch angle θ is taken to be 0.2 rad as a mean value found in Ref. [15].The RE minor radius a RE for the RE area A RE calculation is assumed to be constant and equal to 5 cm.The METIS (with the effective ion charge Z e f f = 2.5) data for the electric field and plasma pressure p pl V from Thomson scattering are interpolated and extrapolated to the EFIT time scale.As EFIT measures before the first and after the last Thomson scattering measurement, the first and the last point of the extrapolated p pl V are approximated to be 95% (arbitrarily taken) of the first and the last Thomson scattering measurement, respectively.The main result is presented in Figure 2. The estimate of the RE current I RE for different RE distribution functions are presented in Figure 2c and their maximum values are listed in Table 2.The plasma current I p is plotted in Figure 2b for comparison purposes.In the same figure, Shielded HXR and Standard HXR are also plotted.The comparison between the W max estimate coming from different RE distribution functions described in Section 4.2 is shown in Figure 2a and their maximum values are tabulated in Table 2.
In Figure 2a, one can see that only W max from the exponential f exp and the Maxwell-Jüttner f MJ distribution function are significantly different, as could be expected from the η values in Table 1.
Even though β pol from EFIT increases significantly during the ramp-down phase of the discharge, a relatively strong decrease in RE current I RE can be seen in Figure 2c.This is a consequence of the current drop in the total plasma pressure p tot V calculation from EFIT in Equation (1).It is interesting that the I RE drop coincides with the RE losses seen from Shielded HXR signal, even though RE losses are not implemented in the calculation.Henceforth, the I RE decrease can be seen as numerical artifact.
However, physical interpretation could be an additional RE energy limit due to the RE beam drift, which is reported in Section 7. Next, notice that I RE values span over two orders of magnitude-from a few kA to tens of kA.
Even though f mono is the most unrealistic of all the used RE distribution functions, it gives minimum I RE as β, γ and η have extreme values.Therefore, we conclude that at least a few kA of current should definitely be driven by the REs in the COMPASS discharge #7298.However, it is not possible to determine which REDF gives the most accurate RE current I RE estimate from Figure 2.
The above experimental results are compared here with the Kruskal-Bernstein theory [16] and NORSE simulation [17] in Figure 3. Kruskal-Bernstein theory is the analytical solution for growth rate of the Dreicer mechanism of RE generation, whose last shape was derived by Connor and Hastie [16].The Kruskal-Bernstein equation is a strong function of the E tor .NORSE calculateS only Dreicer mechanism of RE generation and n e , T e and E tor are taken as time-varying parameters, while Z e f f and B tor are constant.Finally, the n RE is computed by multiplying the total electron density from measurements n e,c with the RE fraction calculated by NORSE.
When applied to the COMPASS discharge #7298, the Kruskal-Bernstein theory predicts RE density n RE larger by 1-3 orders of magnitude than the total n e .This is a non-physical result of course and we consider this theory not relevant for this particular case.Furthermore, NORSE reaches slide-away regime (i.e., when all electrons runaway) at the very beginning of the discharge.Such early slide-away regime is not supported by the experimental results.The possible reasons for the overestimating results from Kruskal-Bernstein theory and NORSE simulation could result from overestimating the electric field E tor by METIS or from too high sensitivity of the theories on the E tor parameter.Moreover, RE energy calculation used here is also simplistic and EFIT results could be misleading in the presence of the RE.

Application: RE Influence on Ramp-Up
REs are frequently generated during the current ramp-up phase in the COMPASS tokamak.Here, the influence of the RE generation on the plasma current I p at the end of the ramp-up phase is investigated using the reported model for the RE current I RE .
The density scan in the range 1-5 × 10 19 m −3 was done during the second RE campaign.The scan consisted of 10 COMPASS discharges from #8552 to #8561.Plasma current I p was feedback controlled to 130 kA during the flat-top phase and the toroidal magnetic field B tor was 1.15 T for all discharges.
Figure 4a presents an estimate of n RE , which is more suitable for comparison of the method reported here with Kruskal-Bernstein and NORSE theories.The three discharges were chosen to cover the standard tokamak discharge (#8553), the slide-away regime (#8559) and the transient between the previous two (#8555).To avoid redundant lines, n RE was calculated only for the two most extreme REDF-the monoenergetic f mono and the exponential f exp ones.
The discharge #8553 has n RE lower by an order of magnitude during the ramp-up phase than the other two discharges (see Figure 4).Surprisingly, discharges #8555 and #8559 have similar n RE at the beginning of the discharge.However, n RE drops for the former discharge as expected from the missing β p rise from EFIT and density drop, as observed by the Thomson Scattering system.Finally, the corresponding RE current I RE at the end of the ramp-up phase for discharges #8553, #8555 and #8559 are 0.3-2.0kA, 2.3-17.8kA and 2.4-19.6 kA, respectively.Therefore, the expected trend is observed.Note that the non-continuous line for discharge #8553 comes from the defined thresholds of the method.To complete the analyses, Kruskal-Bernstein theory and NORSE code were used to theoretically estimate n RE .The results are presented in Figure 4b.
First to notice from Figure 4 is that both theories have the expected trend in density of RE rising from #8553 to #8559.Differently from discharge #7298 in Figure 3, the Kruskal-Bernstein theory does not give n RE over n e , while n RE do rise towards n e at the end of discharge for #8555 and #8559.For these two discharges, the Kruskal-Bernstein theory does not show a significant difference in n RE .On the other hand, NORSE estimates that those two discharges are quite different-according to this code #8559 reaches slide-away regime already around 1060 ms, which is probably too early.
Interestingly, in the case of standard discharge #8553, both codes seem to approximately agree in the order of magnitude with the experimental n RE (between 10 14 and 10 15 m −3 ).The same is valid when the NORSE simulation is compared with experimental n RE using f exp for discharge #8555.However, the NORSE code is not made to predict the drop of n RE , clearly observed in Figure 4a and from β pol from the EFIT reconstruction.

Application: RE Localisation
In this section, we investigate the RE energy limit from the RE outward drift through the analytically derived calculation reported by Zehrfeld [18].
Zehrfeld's analysis provides the W dri f t as a function of the normalized minor radius ρ = r/a p .One of the main messages from this analysis is that W dri f t profile can have (depending on machine and plasma parameters) maximum value inside plasma (i.e., ρ max < 1-dashed coinciding colored lines in Figure 5), which indicates that REs of certain energy can loose confinement before they reach the plasma edge (i.e., the last closed flux surface) that is represented by R out in Figure 5. Namely, peaking parameter of the plasma current profile m1 proved to be the most significant parameter for RE confinement in respect to the drift.Therefore, knowledge of the W dri f t (ρ) profile can be used for localization of the RE beam.Knowing W max at a given time of the discharge, one can find the minimum minor radius ρ min (solid colored lines in Figure 5) where W dri f t (ρ min ) = W max , corresponding to the minor radius below which REs with energy W max cannot be confined.Generally, the major radius of plasma geometrical center R geom and bulk plasma barycenter R sh do not coincide in the tokamak plasma due to the Shafranov shift [19].Additionally, there is RE beam present in the observed discharges shifting the plasma current barycenter R bc more outwards.All these radii R geom , R sh and R bc are taken from EFIT and shown in Figure 5 for orientation purposes.
EFIT predicts m = 1.5 while Zehrfeld's method assumes m to be an integer.Accordingly, we used both rounding integers values: m = 1 and m = 2.
One can notice from Figure 5 that RE orbits are predicted in minor radius of 10-15 cm, which is not very limiting knowing that the COMPASS minor plasma radius is around 20 cm.RE orbits are localized more inside the plasma for higher m, as could be expected due to the higher peaking of W dri f t profile.For the same reason, later significant losses of high energy RE are expected for higher m, which can be seen by slightly delayed equalization of the RE radius with R out (for 5-15 ms) for factor of current profile m = 2 than in the m = 1 case.However, both timings for significant high energy RE losses are a few tens of milliseconds after the Shielded HXR observes the losses, showing that these losses are probably due to the RE outward drift as I p is decreasing and connection of strong Shielded HXR signal and I p ramp-down phase is indeed observed in Figure 2b.On the other hand, this does not explain the interaction on the high-field side seen by the visible cameras.

Conclusions
Understanding RE physics gives us better prediction towards understanding RE generation and mitigation in ITER.One basic RE parameter to compare model and experiment is the RE current.Therefore, we have presented here one method for estimating the RE current in the COMPASS tokamak.
The pitch angle showed to be relevant for the RE energy calculation, which plays an internal role in the RE current estimate method.Implementation of the method on the COMPASS discharge #7298 shows that at least a few kA of the total plasma current is driven by the RE.The method is then used for n e -scan and showed the different amount of RE current present in the plasma at the expected trend.Namely, discharges with lower density have higher RE current.Commonly used Kruskal-Bernstein theory does not appear suitable for high RE current discharges.NORSE, a code built for such discharges, showed to be more sensitive to plasma parameters, but further development is necessary for better prediction of the RE dynamics.Applying the knowledge of RE current to the RE drift energy limit, a crucial role of the limit at the final stage of the COMPASS plasma discharge, is reported.
At the first instance, an implementation of the RE current into the equilibrium simulation is left for future work on the COMPASS tokamak.Furthermore, theories used here are also the most commonly used model basis in the RE research community.However, they obviously overestimate the number of REs.Their further development is far beyond the scope of this article.The results obtained from theories are here to show the issue is oversimplified and demonstrate the need for a better modeling of the RE generation.

Figure 2 .
Figure 2. (a) The maximum kinetic energy W max of discharge #7298 obtained for the different RE distribution functions: mono-energetic f mono (black solid), uniform f uni (black dotted), linear f lin (black dashed), exponential f exp (red), skewed Gaussian f sG (green) and Maxwell-Jüttner f MJ (blue).(b) Time traces of the plasma current I p (blue), Shielded HXR (black) and Standard HXR (red).(c) Estimated RE current I RE corresponding to each RE distribution function with the same labeling as on (a).

Figure 3 .
Figure 3.Time traces of total electron density (black) measured with interferometer, RE density from Kruskal-Bernstein theory (green), RE density from NORSE simulation (cyan) and RE densities corresponding to RE currents from Figure 2c (dashed lines).

Figure 4 .
Figure 4. Time traces of the estimated n RE for the three typical discharges: #8553 (blue, standard discharge), #8555 (black, intermediate case) and #8559 (red, slide-away regime).(a) Estimates from measurements where for each discharge the two most extreme REDF: the monoenergetic f mono (full line) and the exponential f exp (dashed line).(b) Estimates from Kruskal-Bernstein theory (solid line) and NORSE simulation (dashed line).

Figure 5 .
Figure 5.Time traces for lower ρ min (solid colored lines) and upper ρ max (dashed coinciding colored lines) limits of RE major radii for all six REDFs are presented for discharge #7298.For orientation, major radius of plasma geometrical center R geom (black dash-dotted line), outer plasma major radius R out (black solid line), plasma current barycenter R bc (black dashed line) and theoretical bulk plasma barycenter from Shafranov shift R sh (black dotted line) are added.Beside all major radii, the Shielded HXR signal (green points) is presented.Figure is plotted for two of current profile factor m values: (a) m = 1; and (b) m = 2.

Table 2 .
The maxima of W max and I RE over the time domain of the whole discharge, plotted in Figure2.