Quasi-linear analysis of the extraordinary electron wave destabilized by runaway electrons

Runaway electrons with strongly anisotropic distributions present in post-disruption tokamak plasmas can destabilize the extraordinary electron (EXEL) wave. The present work investigates the dynamics of the quasi-linear evolution of the EXEL instability for a range of different plasma parameters using a model runaway distribution function valid for highly relativistic runaway electron beams produced primarily by the avalanche process. Simulations show a rapid pitch-angle scattering of the runaway electrons in the high energy tail on the $100-1000\;\rm \mu s$ time scale. Due to the wave-particle interaction, a modification to the synchrotron radiation spectrum emitted by the runaway electron population is foreseen, exposing a possible experimental detection method for such an interaction.


I. INTRODUCTION
Disruptions in tokamaks can lead to the generation of a high-current beam of highly energetic runaway electrons [1], which poses great challenges for the disruption mitigation system of future tokamaks [2]. The runaway electron beam has a strongly anisotropic velocity distribution and may destabilize high-frequency electromagnetic and electrostatic waves through a resonant interaction. Several high-frequency instabilities driven by runaway electrons have been considered before, using various models for the initial runaway distribution function [3][4][5][6][7][8]. In particular, the linear stability and the quasi-linear analysis of the whistler wave instability (WWI) has been investigated, and it was concluded that whistler waves may be destabilized by an avalanching runaway electron population [7,8]. The main motivation for that work was to investigate the possible effect of these waves on the runaway electron beam formation. If such an instability would lead to scattering of the runaway electrons in pitch-angle, resulting in higher synchrotron radiation losses, a passive mitigation mechanism limiting the detrimental effects of the runaway electrons would be provided. However, it was concluded that for the low temperatures characteristic of post-disruption plasmas, the collisional damping is likely to suppress the WWI and the effect of the instability on runaway beam formation is therefore small. On the other hand, the WWI may provide a diagnostic opportunity due to its sensitive dependence on the fast electron distribution function and the plasma parameters.
Recently, it has been shown that runaways can also destabilize so-called extraordinaryelectron (EXEL) waves at oblique propagation angles [9]. Compared to the WWI, it was found that significantly fewer energetic electrons are needed to destabilize the EXEL wave, which is therefore likely to be the most unstable wave [9]. The aim of this work is to determine the characteristics of the quasi-linear evolution of the EXEL instability and quantify its effects on the runaway electron beam. We also investigate the possibility of detecting signatures of the wave-particle interaction in the experimental infrared synchrotron emission data.
In large tokamak disruptions, where the principal source of runaway electrons is the secondary avalanche process [10], an analytical distribution function for the runaway electrons (in the absence of wave-particle interaction) can be obtained [7]. This distribution function has been benchmarked to the results of numerical simulations [11] and has been used in Ref. [8] as an initial runaway distribution function for the quasi-linear evolution of the WWI. In the present work we adopt a similar approach, extending the treatment to the EXEL wave.
One possible method of inferring the characteristics of the runaway population is to study the synchrotron radiation emitted by the energetic electrons. By calculating the integrated emission from the complete electron population [12], we show that the pitch-angle scattering of highly energetic runaway electrons due to the interaction with the EXEL wave causes a characteristic change in the synchrotron spectrum that could be detected in experiments.
The structure of the paper is as follows. In Sec. II the dispersion relation and the characteristics of the EXEL wave are described. In Sec. III we investigate the quasi-linear evolution of the EXEL wave and its effect on the distribution of fast electrons. Section IV completes the analysis with a study of the parametric dependencies of the process. The calculations of the synchrotron spectrum of the affected distribution, presented in Sec. V, provide guidelines for possible experimental detection of the instability. Finally, the results are discussed and summarized in Sec. VI.

II. EXCITATION OF THE EXTRAORDINARY ELECTRON WAVE
The characteristics of the EXEL wave can be derived from the wave dispersion relation in a homogeneous, magnetized plasma approximation [13]: Note that, to describe the EXEL wave, the frequently used electromagnetic approximation ǫ 33 ≫ (kc/ω) 2 cos θ sin θ has to be relaxed. Here k is the wave number, k and k ⊥ denote its components parallel and perpendicular to the static magnetic field, respectively, and cos θ = k /k. ω is the wave frequency, c is the speed of light, and ǫ ij are the elements of the dielectric tensor, consisting of the susceptibilities of the different plasma species: Here, the indices i and e denote the ion and thermal electron populations, respectively. We neglect the contribution of the runaway electron population to the real part of the frequency. In order to make the calculation of the instability growth rate more convenient, we rewrite the dispersion relation by introducing the cold plasma formulas for the ion and thermal electron susceptibility tensor elements in the high-frequency case of ω ≫ m e /m i ω ce . Equation (1) then becomes where C 1 (k, θ) = − 2k 2 c 2 + ω 2 ce + 3ω 2 pe , C 2 (k, θ) = k 4 c 4 + 2k 2 c 2 (ω 2 ce + 2ω 2 pe ) + ω 2 pe (ω 2 ce + 3ω 2 pe ), C 3 (k, θ) = − k 4 c 4 (ω 2 ce + ω 2 pe ) + k 2 c 2 ω 2 pe (3/2ω 2 ce + 2ω 2 pe + 1/2ω 2 ce cos 2θ) + ω 6 pe , C 4 (k, θ) = 1/2k 4 c 4 ω 2 ce ω 2 pe (1 + cos 2θ), ω pe is the electron plasma frequency and ω ce is the electron cyclotron frequency. Equation (2) is a fourth order equation for ω 2 giving four different branches of electromagnetic waves, as described in Ref. [9]. It has been shown in Ref. [5] that the two highest frequency branches cannot be destabilized by the runaway population. The remaining two branches, namely the electron-whistler and the EXEL wave, can be destabilized but the EXEL wave was shown to have a growth rate an order of magnitude higher than the electron-whistler wave for a runaway distribution function relevant for near-critical electric field [9].  A. Linear growth rate By taking into account the contribution of runaway electrons to the imaginary part of the frequency in the dispersion (2), the linear growth rate γ l of the EXEL wave is given by Ref. [9] as: where is the derivative of Eq. (2) with respect to ω.
The most important resonant interaction between the runaways and the EXEL wave occurs when the wave frequency ω and wave-number k are such that where v, γ and ω ce are the velocity, relativistic factor and the cyclotron frequency of the electrons taking part in the interaction, respectively. This resonance is called the anomalous Doppler resonance.
In the case of the EXEL wave, the anomalous Doppler resonance occurs with ultrarelativistic runaway electrons (p ≫ 1, where p = γv/c is the normalized momentum). In this region of the momentum space the distribution function of the runaway electrons is highly anisotropic. Meanwhile, the Cherenkov resonance ω − k v = 0 occurs with slightly relativistic runaways having significantly lower normalized momentum (p ≈ 1), for the same wave frequency and wave number vector. For other resonances, such as the Doppler resonance, the resonant momentum would be in the negative region (p < 0). Thus, for a velocity distribution which is sufficiently isotropic for low momentum, so that the Cherenkov resonance can be neglected, and anisotropic for higher momentum, the anomalous Doppler resonance will be dominant.
In the present analysis the effect of the Cherenkov resonance was neglected and a model for the ultra-relativistic runaway tail was used as initial electron distribution for the quasilinear analysis. The distribution is given by where the first part is the analytic secondary generation distribution derived in Ref. [7], valid for E ≫ 1. In the above equation, E = e|E |τ c /m e0 c is the normalized parallel electric field (assumed to be constant in time), m e0 is the electron rest mass, τ c = 4πǫ 2 0 m 2 e0 c 3 /n e e 4 ln Λ is the collision time for relativistic electrons, n e is the background electron density, c Z = 3(Z + 5)/π ln Λ, Z is the effective ion charge, α = (E − 1)/(Z + 1) and n r0 is the seed produced by primary generation. In Eq. (5) this form is supplemented by a Fermi function imposing a gradual cut-off at high momentum around p max with a width of σ p . This latter factor is necessary to account for the maximum energy the electrons typically reach, which is determined by the finite time duration of the accelerating electric field [14] and the energy loss due to close collisions [15]. In the present paper, p max = 30 (corresponding to an energy of 15 MeV) was chosen. This is the order of magnitude of the maximum runaway electron energies typically observed in experiments, see e.g. Figure 13 of Ref. [16]. The width was chosen to be σ p = 1. The runaway electron distribution in Eq. (5) is only valid for highly relativistic runaways, and as such can only be used to calculate the resonant interaction through the anomalous Doppler resonance.

B. Most unstable wave and stability thresholds
The linear growth rate of the EXEL wave is calculated by substituting the EXEL dispersion given by the second lowest frequency solution of Eq. (2) and the runaway electron susceptibility [7,13] into Eq. (3). It is positive in the whole wave number space, but the growth rate is highest in the high wave number region, where k c > ω, see Since the growth rate is increasing as we approach the k c = ω line (as seen on Fig. 2) and the closer we are to k c = ω, the higher is the resonant momentum p res , it follows that the resonant momentum of the most unstable EXEL wave is close to the chosen maximum momentum p res ≃ p max . However, we note that the exact value of the chosen p max does not have any drastic effect on either the growth rate or the parameters of the resonant wave.
The reason for that is that the resonant wave parameters are not changed significantly as p max changes (as seen on Fig. 3, the wave number k is almost the same whether we have e.g. p res = 30 or p res = 20). Also no divergence in the growth rate is observed when approaching the k c = ω. Therefore, the order of magnitude of the growth rate for the most unstable wave is the same, irrespective of the choice of p max . However, due to the fact that the line corresponding to p res = p max nearly coincides with the contour lines of the growth rate (as shown in Fig. 2), it is not trivial to find the exact parameters of the most unstable wave.
Fortunately, for precisely the same reason, the value of the growth rate or the number of runaways needed for the interaction are not affected significantly by the exact value of these parameters. Comparing the linear growth rate of the most unstable wave (γ l ) to the collisional damping ln Λ is the electron-ion collision time), and the convective damping rate γ v ≡ (∂ω/∂k ⊥ )/(4L r ) (where L r is the radius of the runaway beam [8]), gives the linear stability threshold -the number of runaway electrons needed for the destabilization of the wave. In the high electric field case studied in the present paper, the momenta of the resonant runaways is expected to be higher than in the corresponding near-critical case studied in Ref. [9], and both the most unstable EXEL and whistler waves therefore have lower frequencies. For the whistler wave this means that instead of the high-frequency electron-whistler approximation, the magnetosonic-whistler wave [7] -which also includes the ion susceptibilities in the dispersion relation -should be used. The stability thresholds for the EXEL and the magnetosonic-whistler wave are shown in Fig. 4. Note that the stability threshold for the EXEL wave is significantly lower in this high electric field case compared to the near-critical case studied in Ref. [9]. Here, the electric field was chosen to be 40 V/m.
where C ω = 0.92, C k = 0.011 and C θ = 0.35 around the most unstable wave in the reference scenario. The values of C ω and C θ tend to be quite robust with respect to changes in plasma density and magnetic field strength; a variation of only about 5% is observed for a change in the plasma parameters of roughly 20%. C k increases very strongly with increasing magnetic field, but remains almost insensitive to changes in the background electron density. This parameter gives a relatively small contribution to ω fit , so the fit is considered to be quite good in the close vicinity of the most unstable wave in the reference scenario. The fit also reproduces some of the dominant changes in the dispersion due to changes in the plasma parameters. However, in the region of interest (which is quite large due to the large spectral range of the waves destabilized in the quasi-linear interaction), ω fit deviates significantly from the exact dispersion. In the remainder of this paper, the exact dispersion given by the solution of Eq. (2), will be used.

EXTRAORDINARY-ELECTRON WAVE INSTABILITY
In the framework of quasi-linear theory, the evolution of the distribution function of the electrons is given by a diffusion equation in phase space, and the rate of growth of wave-energy is equal to the difference between the linear growth rate and the damping rates, The analysis of the dynamics of the interaction of runaway electrons with the EXEL wave can be performed similarly to that of the magnetosonic-whistler wave in Ref. [8]. Only the dispersion relation and the polarization of the wave are different in this case, but as we will show, this proves to have a significant effect on the temporal evolution of the instability. The evolution of the runaway distribution in the general case is given by [13]: where Ω = eB/m e = ω ce /γ,Π E kx , E ky , E kz are the components of the spatial Fourier transform of the electric field and J n (z) is the Bessel function of the first kind and of order n, with the argument z = k ⊥ p ⊥ c/ω ce .
Using the polarization for the EXEL wave (e x , e y , e z ) = 1 , −i where ω = ω EXEL (k, θ), equation (9) gives The wave instability is driven by the anisotropy of the runaway distribution via the anomalous Doppler resonance n = −1. For z ≪ 1 the Bessel function can be expanded as The quasi-linear equation for the runaway distribution becomes and if we assume k v ⊥ ∂f /∂p ≪ Ω∂f /∂p ⊥ , Eq. (13) simplifies to a diffusion equation with where W k (t) = ǫ 0 2 |E k (t)| 2 is the spectral energy of the wave. The assumption k v ⊥ ∂f /∂p ≪ Ω∂f /∂p ⊥ is valid when zk /k ⊥ ≪ (∂f /∂p ⊥ )/(∂f /∂p ), and is satisfied due to the ordering The time variation of the spectral energy of the wave is determined by the differential equation [8]: with the initial condition W k0 = W k (t = 0) = eT e /2, which is the thermal fluctuation level.

A. Numerical solution
Assuming a beam-like velocity distribution γ ≃ p and introducing all terms containing p in (14) into the diffusion operator we obtain a diffusion equation for f in whichD(p , t) is independent of p ⊥ . Introducing a dimensionless time the diffusion equation (14) takes the form: and with the initial condition (5) the solution according to [8] is where φ(p , t) = 2ατ (p , t) + p .
This formula gives the evolution of the runaway distribution as a function of the dimensionless time, τ (p , t). This enables us to create a numerical code which only has to solve for τ in each time step for every p value in a certain region. In order to calculate τ , we need to evaluate the integral in Eq. (17). Due to the azimuthal symmetry of the system, Eq. (17) can be written on the form where G(k, θ) is a function incorporating all dependences on the wave number and the propagation angle (other than the delta function and the Jacobian). We can evaluate the integral in k and arrive at where k res is the solution of the resonance condition ω(k, θ) + ω ce /p − kc cos θ = 0, which can be calculated numerically. From now on we will refer to k res as the resonant curve.
The numerical solution of the quasi-linear equations proceeds as follows.  An example, where the magnetic field strength was varied, is shown in Fig. 6. The figures for different magnetic fields are qualitatively similar; there are only two pronounced differences. The first is that a significantly larger number of runaway electrons is necessary for the destabilization of the EXEL wave for high magnetic field strengths, as indicated by the vertical axes of the distribution function plots. This agrees with the trend shown in Fig. 4. The other significant difference is the smaller angle of wave propagation at weak magnetic field, although the difference of about 0.2 radians is not particularly large. Upon closer inspection, the extent of isotropization due to the EXEL wave seems slightly larger for stronger magnetic fields, but the difference is not significant.
A quantitative analysis has also been performed regarding the change in the quasi-linear evolution due to changes in the plasma parameters. The results are summarized in Table I, which shows the value of the following characterizing parameters: • p m , the momentum resonant with the most unstable wave • n r1 , the runaway density at momentum p m integrated over p ⊥ at the time of the first wave destabilization. n r1 is thus the 'threshold linear density' • W max , the maximum wave energy over the 85.9 µs duration of the simulation • τ , the parameter characterizing the extent of velocity space diffusion (calculated according to Eq. 18) at p m at the end of the simulation • k m , the wave-number of the most unstable wave • θ m , the propagation angle of the most unstable wave From Tab. I we can infer that the most unstable momentum of the runaways is not affected significantly by the change of the plasma parameters through the quasi-linear evolution -it is close to the p max cut-off value introduced in the initial distribution function (Eq. 5).
The changes due to variations in the magnetic field already described in the discussion of Fig. 6 are also visible in the table, with the additional observation that larger magnetic fields shift the resonant EXEL waves towards larger wave numbers. On the other hand the maximum wave energy is an order of magnitude higher for lower magnetic field, and is obtained by the most unstable wave during the first phase of isotropization (whereas Fig. 6 shows a later time instant).  µs) after the first wave destabilization. In each column, only the parameter indicated by the column heading was changed -the remaining parameters where those of the reference scenario.
The dominant effect of a change in the density is a modification to the strength of the collisional damping. Accordingly, a higher density means a higher critical runaway density (n r1 ). On the other hand, the quasi-linear diffusion is significantly faster for high densities and the wave energies are significantly higher.
Increasing the accelerating electric field results in a decrease of the critical runaway density needed for the destabilization, which is explained by the increasing anisotropy of the runaway beam. There is no substantial effect on the other parameters in Tab. I.
The background plasma temperature enters through the collisional damping, so n r1 increases with decreasing temperature. It is a general observation that a higher threshold runaway density is accompanied by a higher wave energy. The only exception is modifications to the magnetic field strength, where the trend is the opposite.
In summary, the EXEL wave is expected to be destabilized in plasmas where the density and temperature are not too low, and where the magnetic field is weak. These conditions wave is expected among the most energetic runaways, but these are also the most strongly emitting particles in terms of synchrotron radiation. Therefore, the wave-particle interaction can result in a substantial change in the synchrotron spectrum [12].
The average synchrotron power emitted per runaway particle at a specific wavelength λ can be calculated as a convolution of the distribution function with the synchrotron emission from a single particle: where f is the momentum-space distribution of electrons, P describes the synchrotron emission and S RE is the runaway region in momentum space [12]. The synchrotron power radiated by a highly relativistic particle in a toroidal plasma was derived in Ref. [19] and is given by where c P = ce 2 /(ε 0 λ 3 γ 2 ), a = η/(1 + η 2 ), g(y) = y −1 + 2y, h(y) = 3ξ(y + y 3 /3)/2, R is the tokamak major radius, γ is the relativistic mass factor, J ν (x) is the Bessel function and J ′ ν (x) its derivative. Eq. (24) takes the drifts stemming from the curvature and gradient of the magnetic field into account and is valid when p ≫ p ⊥ , c/λ ≫ ω ce , and when the aspect ratio is large. It depends on the particle energy and pitch through the parameters γ and η, respectively. Due to their structure, the integrands in Eq. (24) are highly oscillatory and evaluating the integrals can be numerically demanding. Here we use a Matlab routine called SYRUP (SYnchrotron emission from RUnaway Particles), also used to obtain the results in Ref. [12], to calculate the synchrotron spectrum from the normalized distribution before and after the onset of the resonant interaction between the runaway distribution and the EXEL wave. SYRUP implements Eqs. (23) and (24).
The result of the synchrotron spectrum calculation for the reference JET-like scenario is shown in Fig. 7. As an initial distribution, Eq. (5) was used with R = 3 m. The distribution affected by the interaction was taken at 604 µs after the first destabilization of the most unstable wave (this is the distribution function plotted in Fig. 5). The runaway region in momentum space (S RE ) was defined by p ∈ [12,31] and p ⊥ ∈ [0, 3] (cf. Fig. 5). Due to the strong energy dependence of the synchrotron emission, and the exponential fall-off of the distribution with increasing p ⊥ , contributions from particles with lower parallel momentum or larger perpendicular momentum were negligible. The maximal parallel momentum was determined by the cut-off of the distribution function (20) at p max + σ p . The synchrotron spectra in Fig. 7 show that in 604 µs, the particle-wave interaction has a significant effect on the synchrotron spectrum, with the peak emission increasing by roughly a factor 2.5. The wavelength of peak emission is also shifted slightly towards shorter wavelengths. We emphasize that Fig. 7 shows the average emission per runaway, meaning that if the number of runaways can be considered fixed, the effect of the interaction with the EXEL wave is a significant increase of the total synchrotron emission.
As discussed in Ref. [12], however, there are several other factors that have a similar effect on the synchrotron spectrum. The spectrum is highly dependent on the properties of the runaway distribution, and is thus sensitive to plasma parameters such as temperature, density, impurity content and electric and magnetic fields. The effect of the EXEL wave interaction could only be discriminated by the characteristic sudden increase of the emission on the 100 − 1000 µs time scale.
The EXEL wave interaction is likely to produce an even larger effect than that shown in Figures 5 and 7, however, since as the assumption of beam-like distribution used in our modeling breaks down, we can not simulate the later stages of the interaction. Pitch-angle scattering might also increase radial transport and eventually lead to the mitigation of the runaway beam, but at an earlier stage our model predicts a burst of synchrotron radiation unaccompanied by macroscopic MHD activity.
For the distribution used in Fig. 7, the wavelength region of strong emission is in the farinfrared and sub-millimeter regions of the spectrum, implying that detection of the effect of the EXEL wave instability by means of synchrotron radiation would require an infrared camera sensitive to this wavelength range. The reason for the long wavelength emission is the cut-off of the runaway distribution (Eq. 5) at a particle energy of roughly 15 MeV. A runaway electron distribution extending to higher maximum energy would allow detection by ordinary near-infrared (or even visible light) cameras, but realistic modeling of the evolution of the distribution function in a disruption is out of the scope of this paper.

VI. CONCLUSIONS
Runaway electrons pose a significant threat to tokamaks. This is especially true for ITER, where the runaway current might be as high as 12 MA in disruptions, with the electron energy spectrum extending up to several tens of MeV [2]. In this paper the quasi-linear resonant interaction of the runaway population and the high-frequency obliquely propagat-ing extraordinary-electron (EXEL) wave, which leads to rapid pitch-angle scattering of the resonant runaways, was studied. The scattering occurs when the runaway density reaches a certain critical density of about 10 14 − 10 17 m −3 , depending on the plasma parameters. As soon as the EXEL wave is destabilized, it leads to a pitch-angle scattering of resonant electrons through quasi-linear diffusion in the velocity space on the 100 − 1000 µs time scale. In our simulations, the spectral energy of the destabilized EXEL wave did not exceed 10 −11 J in any of the scenarios considered, implying that direct experimental detection of the wave is likely to be difficult.
As the resonant interaction with the EXEL wave mainly affects the high energy runaways, which are the electrons characterizing the synchrotron radiation emitted by the whole population, the interaction causes a significant change in the synchrotron spectrum.
The interaction with the EXEL wave was shown to produce a burst of synchrotron radiation accompanied by a simultaneous shift of the spectrum towards shorter wavelengths, which might offer a possibility to detect the impact of the quasi-linear interaction in experiments.
By looking at a wide range of plasma parameters, we concluded that the characterizing quantities of the interaction (resonant runaway momentum, wave energy, critical runaway density, etc.) have a weak dependency on plasma parameters. We can therefore extend our conclusions to an ITER-like scenario. The intensity of the interaction, and the resulting change in the synchrotron spectrum, are expected to be qualitatively similar and of the same order of magnitude as for the investigated JET-like reference scenario. The only major difference in the ITER case is a slightly higher stability threshold (which could still easily be reached) due to he stronger magnetic field.
Our analysis shows that the EXEL wave is destabilized considerably more easily, as compared to the previously studied whistler wave [7,8]. A possibility for experimental confirmation of the results presented in this paper is offered through the predicted bursts in the far-infrared synchrotron emission. Our results also provide a basis for further theoretical work making use of more realistic kinetic simulations with advanced Fokker-Planck solvers, such as the LUKE code [20].