The effects of kinetic instabilities on the electron cyclotron emission from runaway electrons
Chang Liu Lei Shi Eero Hirvijoki Dylan P. Brennan Amitava Bhattacharjee Carlos Paz-Soldan Max E. Austin





Princeton Plasma Physics Laboratory, Princeton, New Jersey 08540, USA
cliu@pppl.gov

University of California, Irvine, California 92697, USA	

Princeton Plasma Physics Laboratory, Princeton, New Jersey 08540, USA

Princeton University, Princeton, New Jersey 08544, USA

Princeton Plasma Physics Laboratory, Princeton, New Jersey 08540, USA
Princeton University, Princeton, New Jersey 08544, USA

General Atomics, San Diego, California 92186, USA

University of Texas, Austin, Texas 78712, USA


In this paper we show that the kinetic instabilities associated with runaway electron beams play an essential role for the production of high-level non-thermal electron-cyclotron-emission (ECE) radiation. Most of the non-thermal ECE comes from runaway electrons in the low-energy regime with large pitch angle, which are strongly scattered by the excited whistler waves. The power of ECE from runaway electrons is obtained using a synthetic diagnostic model based on the reciprocity method. The electron distribution function is calculated using a kinetic simulation model including the whistler wave instabilities and the quasilinear diffusion effects. Simulations based on DIII-D low-density discharge reproduces the rapid growth of the ECE signals observed in DIII-D experiments. Unlike the thermal ECE where radiation for a certain frequency is strongly localized inside the resonance region, the non-thermal ECE radiation from runaway electrons is nonlocal, and the emission-absorption ratio is higher than that of thermal electrons. The runaway electron tail is more significant for ECE with higher frequencies, and the ECE spectrum becomes flatter as RE population grows. The nonlinear behavior of the kinetic instabilities is illustrated in the osculations of the ECE waves. The good agreement with the DIII-D experimental observations after including the kinetic instabilities clearly illustrate the significance of the scattering effects from wave-particle interactions, which can also be important for runaway electrons produced in disruptions.  



The effects of kinetic instabilities on the electron cyclotron emission from runaway electrons
    Max E. Austin
    December 30, 2023
==============================================================================================





§ INTRODUCTION


In recent tokamak studies, runaway electrons (RE) have drawn a lot of attention due to their importance in disruptions<cit.>. The classic theory of runaway electrons is based on the dominance of a strong electric field force over the collisional force in the high-energy regime, where an electron an electron can be accelerated to extremely high energy if a threshold in energy is exceeded<cit.>. Based on this theory, it is predicted that in a typical disruption event in ITER, the strong inductive electric field can cause an avalanche growth<cit.> of a high-energy RE population up to tens of MeV, which can bring severe damage to the plasma facing material<cit.>. Given the importance of RE physics, many experiments have been conducted in current tokamaks<cit.>, including flattop scenarios<cit.> and disruption cases<cit.> , to validate the theoretical model and provide predictions for ITER. In these experiments, the dynamics of runaway electrons are diagnosed mainly through the radiation emitted, including the hard X-ray (HXR) and gamma rays from bremsstrahlung<cit.>, visible and infrared lights from synchrotron emission (SE)<cit.>, and radio-frequency waves from  electron cyclotron emission (ECE). Analysis of the diagnostic results show that, in addition to the electric field force and the collision force, the radiation reaction force<cit.> and the kinetic instabilities like whistler waves<cit.> are also very important in the RE dynamics, making the runaway electron distribution very different from what classic theory predicts.

The ECE has been widely used in tokamak experiments to diagnose the electron temperature (T_e)<cit.>. The diagnostic is based on the radiation of electron cyclotron waves (ECWs), as a consequence of the gyro-motion of electrons around magnetic field lines. ECWs can be emitted and absorbed by the electrons when the resonance condition is satisfied. Depending on the absorption coefficients, ECWs with different frequencies and polarizations have different optical thicknesses with respect to the plasma. According to the Kirchhoff's law of thermal radiation, if the plasma is locally in thermal equilibrium and optically thick so that the wave absorption in the resonance region is strong, then the radiated power from the resonance region is proportional to the local electron temperature. Based on this mechanism, the ECE-imaging (ECEI) approach, which scans the ECE signals with various frequencies and propagation angles, has been used to obtain the spatial profile and fluctuations of electron temperature.

The high-energy runaway electrons can also produce ECE during gyro-motions. In many RE experiments in tokamaks including both the flattop<cit.> and disruptions scenarios<cit.>, ECE signals have been observed that are much more intense than the ECE from thermal electrons.  Compared to HXR or gamma rays which mainly depend on the electron energy, the ECE power depends on both the energy and the pitch angle of the runaway electrons, similar to SE. On the other hand, the ECE has a wide-angle distribution which is rather different from the strong forward-beaming of SE. 

Given that the RE tail is not in thermal equilibrium and the absorption of ECWs is much smaller than that of thermal electrons, the ECE radiation power is not primarily determined by the electron temperature in the manner specified by Kirchoff's law. The ECE from suprathermal electrons has been studied in lower hybrid current drive experiments<cit.> using the Doppler broadening of the ECE spectrum, but for RE tail with a wide energy distribution, a synthetic diagnostic model is required to interpret the more complicated non-thermal ECE signals that elucidate the electron distribution function.

In this paper, we present an ECE synthetic diagnostic model for arbitrary electron distributions. The model is applied to study relativistic RE tail including the kinetic instabilities. Unlike the forward method used in previous codes like HORACE<cit.>, in our model we use a backward method to calculate the ECE power based on the reciprocal theorem<cit.>, an innovative approach which separates the wave emission and absorption by introducing an artificial wave propagating against the radiation direction along the same ray path. The reflection loss and the polarization scrambling<cit.> of the waves are included in the reciprocal calculations. The RE distribution function is calculated using a newly-developed kinetic simulation model, coupled to a kinetic instability calculation framework using the quasilinear diffusion model<cit.>. The results show that, when kinetic instabilities are taken into account, the synthetic diagnostic shows good agreement with experiments for both X-wave and O-wave, including the growth rate and peak amplitude of radiation power<cit.>, the evolution of the ECE spectrum shape, and the oscillatory behavior of ECE signals observed in RE decaying stage<cit.>. The radiation power is further expressed using the ECE weight function, which shows that most of the non-thermal ECE radiation comes from runaway electrons with low-energy and large pitch angle, and the excitation of whistler waves by REs plays an essential role in producing this radiation.

This paper is organized as follows. In Sec. <ref> we introduce the ECE synthetic diagnostic model, which is based on the reciprocity theorem in electrodynamics. In Sec. <ref> the kinetic simulation which provides the RE distribution for the synthetic diagnostic is introduced. Then in Sec. <ref> the synthetic ECE signals obtained using these two models based on a DIII-D low-density flattop RE experiment is shown, which reproduces features observed in experiments. In Sec. <ref>, we explain the physics mechanism causing the strong non-thermal ECE from RE in our simulation results, and calculate the ECE weight function in this case. In Sec. <ref> we summarize. 



§ ABSORPTION AND EMISSION OF ELECTRON CYCLOTRON WAVE



To calculate the wave power radiated by the plasma and received by a specified antenna, one can either use the forward or the backward method. In the forward method, the radiation intensity at every position in plasma is solved through the radiation transport equation<cit.>, which is derived from the wave kinetic equation following the Wentzel-Krammers-Brillouin (WKB) method. In solving the transport equation, one needs to calculate an integral along the ray path, and both the emission and the absorption of the wave needs to be taken into account. On the other hand, in the backward method the absorption and emission of the wave are treated separately, which can facilitate the computation. In this method, one needs to calculate an artificial wave which propagates from the antenna back into the plasma, to get the dependence function of the radiation power on the emission of every point, based on the so-called reciprocal theorem<cit.>. It can be proved mathematically that the two methods are equivalent, and for a 1D radiation problem of the kind studied in the present paper, the calculations are not much different. However, for more complicated radiation problems such as ECEI, the backward method is more computationally advantageous since it does not need to include calculations of the radiation not entering the antenna.

We now focus on the application of the backward method to the ECE. According to the reciprocal theorem, to get the dependence function we need to calculate the electric fields of an artificial wave. This artificial wave has the same frequency and polarization but  propagates in a direction opposite to that of the ECWs. The fields of the artificial wave in plasma can be solved according to the Maxwell equation,

    ∇×∇×𝐄+ω^2/c^2ϵ·𝐄=4π/c𝐣,

where 𝐄 is the wave electric fields, ω is the wave frequency, and c is the speed of light. ϵ is the plasma dielectric tensor. ϵ can be separated into a Hermitian part ϵ^H and an anti-Hermitian part ϵ^A. The Hermitian part determines the dispersion relation and the polarization of the wave, and the anti-Hermitian part determines the damping. In this work we use the cold electron dielectric tensor for the Hermitian part of ϵ<cit.> ignoring the ion contribution and the thermal effect, assuming their correction to the dispersion relation is small. The anti-Hermitian part includes both the absorption of waves due to particle resonance (ϵ^Ar) and damping from collisions(ϵ^Ac). The resonance absorption depends sensitively on the distribution function near the resonance point. The collision damping can be calculated from the electron collision frequency and the wave dispersion relation<cit.>, which is independent of the distribution function. The current 𝐣 on the right-hand-side of Eq. (<ref>) represents the emission of the wave from the fluctuation current, which can be ignored in the calculation of the artificial wave and will be addressed later.

For simplicity, we consider the artificial wave that is propagating along the major radius direction only, propagating from the low-field-side to the high-field. Set x as the direction of the wave propagation, and z as the magnetic field direction, the polarization of the X-wave is

    E_x=iω_ceω_pe^2/ω(ω^2-ω_ce^2-ω_pe^2)E_y-ω^2-ω_ce^2/ω^2-ω_ce^2-ω_pe^2(ϵ^A·𝐄)_x,   E_z=0,

where ω_ce and ω_pe are the electron cyclotron frequency (we choose ω_ce<0) and plasma frequency obtained from the local density and magnetic field. Eq. (<ref>) can then be simplified as

    ∂^2 E_y/∂ x^2+(ω^2-ω_pe^2)-ω^2ω_ce^2/(ω^2-ω_pe^2-ω_ce^2)c^2E_y+ω^2/c^2iω_ceω_pe^2/ω(ω^2-ω_ce^2-ω_pe^2)(ϵ^A·𝐄)_x+ω^2/c^2(ϵ^A·𝐄)_y=0.

Note in Eqs. (<ref>) and (<ref>), the ϵ^A  terms are much smaller than the rest of the terms. Taking the first order approximation, Eq. (<ref>) can be solved using the WKB approximation, and the result is

    E_y=C/√(k(x))exp[∫ dx [i k(x)+k_i(x)]],

where k(x)=√([(ω^2-ω_pe^2)^2-ω^2ω_ce^2]/[(ω^2-ω_pe^2-ω_ce^2)c^2]) is the real component of the wave vector. k_i is the imaginary component,

    k_i(x)=ω^2/c^2(𝐄^†·ϵ^A·𝐄)/2 k(x) |E_y|^2,

which describes the absorption of the X-wave due to particle resonances.

For O-wave, the polarization is

    E_x=0,   E_y=0,

and E_z satisfies the differential equation,

    ∂^2 E_z/∂ x^2+ω^2-ω_pe^2/c^2E_z+ω^2/c^2(ϵ^A·𝐄)_z=0.

We can similarly derive the WKB equation for E_z,

    E_z=C/√(k(x))exp[∫ dx [i k(x)+k_i(x)]],

where

    k(x)=√(ω^2-ω_pe^2/c^2),


    k_i(x)=ω^2/c^2(𝐄^†·ϵ^A·𝐄)/2 k(x) |E_z|^2.

For the artificial wave, constant C is arbitrary and here we choose C=1. The value of 𝐄 along the ray can then be calculated by integrating Eq. (<ref>).

When reaching the wall of tokamak, the artificial wave gets reflected into the opposite direction of the incoming ray. The intensity of the wave (I_X=E_y^2k for X-wave and I_O=E_z^2k for O-wave) will be multiplied by the reflection factor α_r<1 due to the reflection power loss. In addition, the intensity of the two polarizations will transfer between each other due to the polarization scrambling, which can be characterized by a polarization transfer fraction α_p. The intensity after a reflection I' will thus follow the equations<cit.>,

    I'_X=α_r [(1-α_p)I_X+α_p I_O]


    I'_O=α_r [(1-α_p)I_O+α_p I_X]


After obtaining the artificial wave electric fields 𝐄, the total radiation power of the ECE with the same frequency can be calculated according to the reciprocity theorem, and sum over X and O waves,

    P(ω)=∑_X,O∫ dx 𝐄^†·K·𝐄,

where K is the correlation tensor for the fluctuating current density, which represents the emission of ECWs by the resonant particles (the term on the right-hand-side of Eq. (<ref>)). The reciprocal theorem proves that the power obtained from Eq. (<ref>) is the same as the power calculated using the forward method, but in the backward method the absorption and emission are separately calculated and only coupled through Eq. (<ref>).

The values of ϵ^Ar and K can be calculated from the electron distribution function<cit.>,

    ϵ^Ar=ω_pe^2/ω∫ d^3𝐩∑_n=-∞^∞δ(ω-k_∥v_∥-nω_cb)[1/v∂ f/∂ p+ωcosθ-k_∥v/ω p v_⊥] T_n,


    K=ω_pe^2/π∫ d^3𝐩∑_n=-∞^∞δ(ω-k_∥v_∥-nω_cb) f T_n,

where


    =10pt
    T_n=([   n^2ω_cb^2/k_⊥^2J_n^2  i nω_cbv_⊥/k_⊥J_nJ'_n      nω_cbv_∥/k_⊥J_n^2; -i nω_cbv_⊥/k_⊥J_nJ'_n            v_⊥^2J'^2_n      -i v_⊥ v_∥J_nJ'_n;      nω_cbv_∥/k_⊥J_n^2       i v_⊥ v_∥J_nJ'_n            v_∥^2 J_n^2 ]),
k_∥ is the wave vector parallel to magnetic field (in our ECE calculations k_∥=0), and k_⊥ is the perpendicular component.  p is electron momentum, v is the velocity, and θ is the pitch angle. v_∥=vcosθ, v_⊥=vsinθ. ω_cb=ω_ce/γ is the cyclotron frequency for a relativistic electron (γ is the relativistic factor). f is the electron distribution function normalized to ∫ d^3𝐩 f=1, and J_n is the nth order Bessel function, with argument k_⊥ v_⊥/ω_cb.


According to Eqs. (<ref>) and (<ref>), both the wave emission and absorption happens with the electrons satisfying the resonance condition with the wave,

    ω-k_∥v_∥-nω_cb=0.

For unrelativistic electrons k_∥=0, the resonance condition can be simplified as ω=nω_ce with n<0. The wave-particle
resonance only happens in the regions where the value of B satisfies the resonance condition, which is called the resonance region. Regarding the absorption efficiency in the resonance region, the plasma can be categorized as optically thick or optically thin for different waves. In the optically thick case, the wave will be strongly damped in the resonance region (the characteristic length of damping L=1/k_i is much smaller than the scale length of the resonance region), which can be regarded as a black body with radiation power only depending on the temperature. In the optically thin case, the wave will be weakly absorbed in the resonance region, and the radiation power depends on density and length of the resonance region as well<cit.>. Typically for X-wave with n≥ -2, the plasma can be treated as optically thick, and for O-wave or smaller n the plasma is optically thin.


For suprathermal electrons like REs, the wave-particle resonance can happen in regions
outside the resonance regions due to the relativistic factor. These electrons
can provide emission and absorption in addition to the thermal radiation. However, the emission of ECWs also depends on the pitch angle of the resonant electrons, through the Bessel function J_n in Eq. (<ref>). The pitch angle requirement is particularly important for high-energy electrons and high frequency ECWs, where |n| is large, and the Bessel function J_n at higher |n| requires larger argument (comparable to |n|) to get a significant output value. This condition can limit the ECE from the runaway electron tail, since the pitch angle distribution of RE is usually strongly collimated due to the weak collisional scattering in the high-energy regime.

On the other hand, if the suprathermal electrons can have have large pitch angle, then the ECE from them will differentiate from thermal electrons not only in frequencies and resonances, but also in polarizations. According to the polarization of X-wave and O-wave, we know that the radiation of X-wave depends on the xx, yy and xy component of K, whereas the radiation of O-wave depends on the zz component. Given Eq. (<ref>) and taking the leading order term of the Bessel functions, we find that the X-wave emission is proportional to (k_⊥ v_⊥/ω_cb)^2n-2, while the O-wave emission is proportional to (k_⊥ v_⊥/ω_cb)^2n. For thermal electrons, it is known that k_⊥ v_⊥/ω_cb≪ 1, thus the X-wave emission dominates the O-wave. However, for superthermal electrons with large pitch angle, v_⊥∼ c, we find that the O-wave emission can be comparable to X-wave, which leads to an ECE with a combined polarization.

According to Eqs. (<ref>) and (<ref>), for a Maxwellian electron distribution with temperature T_e, ϵ^Ar and K satisfy

    K_ij/ϵ^Ar_ij=ω T/π

for every term in the tensor. This result is an outcome of the fluctuation-dissipation theorem, which states that for a system in thermal equilibrium, the fluctuation amplitude level is always proportional to the efficiency of the corresponding dissipative process, with the ratio proportional to the system temperature.

We can also define the normalized emission term,

    For X-wave,     U=πω/c^2𝐄^†·K·𝐄/2 k |E_y|^2


    For O-wave,     U=πω/c^2𝐄^†·K·𝐄/2 k |E_z|^2

If we only consider k_i from the resonance absorption ϵ^Ar,  we find U/k_i=T_e according to Eq. (<ref>)(<ref>). The radiation power is then

    P=c^2/πω∫ dx 2U(x)  exp(-2∫ dx k_i(x)).


Assuming all the power of the artificial wave is absorbed by the plasma with uniform temperature T_e, the ECE radiation power can be calculated as

    P=c^2T_e/πω∫ dx 2k_iexp(-2∫ dx k_i)=c^2T_e/πω

thus the radiation power is proportional to the plasma temperature. This result is consistent with Kirchhoff's radiation law. In the later text, we will use the effective radiation temperature T_𝐞𝐟𝐟=(πω/c^2)P to represent the ECE radiation power, which characterizes the plasma temperature in the resonance region for the optically thick case. 

In addition to resonance absorption, the collision effect (ϵ^Ac) can also contribute to the wave damping and increase k_i. This will results in U/k_i<T_e and the final radiation power will be smaller than that from Kirchhoff's radiation law. In the normal tokamak experimental parameters like what will be discussed in Sec. <ref>, where T_e∼keV, the collisional damping is much smaller than the resonance absorption and can be ignored. However, for post-disruption scenarios where T_e∼ a few eV, the resonance absorption from thermal electrons is very weak and the collisional damping is significantly enhanced. In this case, the contribution from collisions dominates the damping of ECWs.

In the numerical model, the integrals in Eq. (<ref>) are calculated using the  finite difference method. For each reciprocal calculation, 10 reflections are accounted in total, with α_r=0.76 and α_p=0.20. The profile of the magnetic field B is set as B(x)=B_0R/x where R is the major radius and B_0 is the magnetic field on the axis. In calculating k_i and U at every positions, the resonance harmonics n are calculated from -2 to -20n_s, where n_s=ω/ω_ce0 and ω_ce0=eB_0/mc. The electron distribution function f is obtained from a kinetic simulation, which is discussed in Sec. <ref>.



§ SIMULATION OF RE MOMENTUM SPACE DISTRIBUTION INCLUDING KINETIC INSTABILITIES



In order to model the non-thermal ECE radiation, the electron distribution function in momentum space needs to be modeled correctly. Especially for runaway electrons, the distribution function has significant deviations from a Maxwellian and cannot be simply modeled as a bi-Maxwellian. Here the distribution function is calculated by solving the spatially homogeneous kinetic equation. The coordinates for momentum space are (p,ξ), where p is the electron momentum normalized to mc (m is the electron mass and c is the speed of light), and ξ=p_∥/p is the cosine of the pitch angle. For runaway electrons, the important forces affecting momentum space dynamics are the electric force, the collisions, and the radiation reaction forces. Knock-on collisions which are ignored in the Fokker-Planck form of the collision operator provide an additional source for RE generation. The kinetic equation we solve can be written as

    ∂ f/∂ t+eE_∥/mc(ξ∂ f/∂ p+1-ξ^2/p∂ f/∂ξ)+C[f]+∂/∂𝐩·(𝐅_radf)+D[f]=S_A[f],

where E_∥ is the parallel electric field and e is the electron charge, C[…] represents the test particle collision operator<cit.>, 𝐅_rad is the synchrotron radiation reaction force<cit.>, D[…] is the quasilinear diffusion operator from the excited waves, and S_A[…] is the source term for the knock-on collisions.

In this model we take into account the kinetic instabilities associated with the RE beam, and its back-reaction on the RE distribution function using the quasilinear diffusion model in magnetized plasmas<cit.>. The RE distribution determined by all the kinetic forces has a wide distribution in energy, but is strongly forward beamed and thus has a very narrow distribution in pitch angle (see Fig. <ref> (a)). This anisotropic distribution can be susceptible to various kinds of kinetic instabilities, most notably one driven by its interaction with the whistler waves<cit.>. In addition, a bump-on-tail distribution formed by the energy dissipation from radiation reaction<cit.> can also driven plasma waves. The excited plasma waves can scatter the resonant electrons strongly in pitch angle when satisfying the resonance condition (Eq. (<ref>)), and isotropize the distribution, called "fan instabilities"<cit.>. Given that the ECE depends sensitively on the pitch angle of electrons, these kinetic instabilities can have a remarkable impact on the ECE radiation.





The impact of kinetic instabilities on the RE distribution is illustrated in Fig. <ref> and Fig. <ref>, where we show a comparison of the RE distribution functions obtained from the kinetic simulation, with quasilinear diffusion operator turned on and off. The plasma parameters are the same as the simulation in Sec. <ref> from a DIII-D flattop experiment. As shown in Fig. <ref>, the excited modes lead to an increase of RE pitch angle in both low-energy and high-energy regimes. This scattered electron distribution affected by kinetic instabilities (Fig. <ref> (b)) is quantitatively different from the unscattered distribution (Fig. <ref> (a)).  Especially in the low-energy regime, where the resonant electrons can be scattered to p_∥≤0, which can significantly enhance the emission of ECWs. In addition to scattering, interactions with the whistler waves also cause a higher avalanche growth rate and a larger RE population, as shown in Fig. <ref>. Affected by the waves, the electrons tend to accumulate in low-energy regime and stop going to very high energy. The value of the critical electric field is also altered by the kinetic instabilities<cit.>.




Knock-on collisions between highly energetic REs and low-energy thermal electrons cause an exponential growth of the RE population, called "runaway electron avalanche". This effect can be accounted for using a source term in the kinetic equation derived from the Boltzmann collision operator using a Møller cross section<cit.>. In previous studies of the source term<cit.>, the pitch angle distribution of the seed REs is ignored based on the assumption that the RE distribution is strongly anisotropic. Given the new distribution from the simulation including kinetic instabilities, this assumption is not valid. In terms of this, we apply a new knock-on collision source term<cit.> in Eq. (<ref>) which takes into account both the energy and the pitch angle distribution of the seed REs. This source term is calculated by transforming  the distribution function to a spectrum representation using Legendre polynomials<cit.>, in which the knock-on collision operator can be simplified.

Note that in the current model the distribution function is only solved in momentum space. In order to do a synthetic diagnostic of ECE, we need the electron distribution functions along the ray path. Based on the knowledge of RE generation and current diagnostic results<cit.>, it is known that most of the generated runaway electrons live in the region close to the core. Thus for simplicity, we choose a tokamak cross-section with r<0.5a as the RE region, where a is the minor radius. Inside this region, f is homogeneous in space, initialized as a Maxwellian with T_e=T_e0, and evolves with the kinetic simulation. Outside this region, f is assumed to be a Maxwellian with temperature profile T_e=(T_core-T_edge)(1-(r/a)^2)^2+T_edge, and unchanged with time. Plasma density is assumed to have a similar profile n_e=n_e0(1-(r/a)^2)^2.



§ ECE SIGNALS FROM RUNAWAY ELECTRON SIMULATIONS



We now use the ECE synthetic diagnostic model to calculate the ECE signals from the RE simulations. The simulation is set based on the low-density flattop RE experiments on DIII-D<cit.>. The experimental discharge contains two stages. In the first stage, the plasma density is low and the runaway electron tail is driven by a parallel electric field from the Ohmic coils. In the second stage, the plasma density increases due to gas puffing and the value of E/E_CH decreases (E_CH is the Connor-Hastie critical electric field<cit.>). An experimental example of the two stages divided by gas puffing can be found in <cit.>. The parameters we used are close to the numbers from the tokamak core diagnostic. For stage 1, density n_e0 is 0.6× 10^19m^-3, temperature T_e0 is 1.3keV, T_core is 2.0keV and T_edge is 0.2keV. The electron distribution f is initialized as a Maxwellian. Electric field E_∥=0.055V/m. For stage 2, n_e0 is 0.8× 10^19m^-3. The temperature profile is the same as stage 1 and the electron distribution is initialized from the final step of stage 1. The electric field E_∥ is 0.045V/m. For the dimension of the device, the major radius R is 2.0m and the minor radius a is 0.5m. The magnetic field on axis in two stages is B_0=1.45T.

The results are summarized as follows. In Fig. <ref> we show the effective radiation temperature T_eff of ECE signals calculated from the synthetic diagnostic, for two different polarizations and at two different frequencies in stage 1. The frequencies of ECE are chosen to be 2ω_ce0 and 3ω_ce0. The solid lines are calculated from the scattered RE distribution including the wave-particle interactions (Fig. <ref> (b)). The dashed lines are the results from the unscattered RE distribution function (Fig. <ref> (a)).



At the beginning of stage 1, the value of the T_eff for the X-mode with ω=2ω_ce0 is very close to T_e0, which is consistent with Kirchhoff's law. The value of the ECE signal for ω=3ω_ce0 is much smaller, since the plasma is optically thin for this high frequency wave. The O-mode signals mostly come from the polarization scrambling, given its own radiation is small compared to X-mode. During the RE population growth in stage 1, the ECE signals barely grow until 3.0s. After 3.0s, the ECE signals for both polarizations from the scattered RE distribution grow abruptly, whereas the ECE signals from the unscattered distribution remain dormant. An interesting feature of ECE signals in this abrupt growing phase is that for the X-mode, the ω=3ω_ce0 signal surpasses the ω=2ω_ce0 signal. Comparing with the experiments (Fig. 5 (d) in <cit.>), we find that the X-mode ECE results with the scattered distribution has much better agreement with the experiments, including the growth rate, peak amplitudes and the overpassing behavior. This is strong evidence that the kinetic instabilities are present and play an important role in the DIII-D QRE experiments<cit.>.

The different behaviors of the ECE signals at two frequencies inspire us to study the change of the whole ECE spectrum during RE growth. In Fig. <ref> we show the spectrum of the X-mode ECE signals at different times in stage 1, for the unscattered and the scattered distributions. At the beginning when f is close to a Maxwellian, the spectrum is like a step function, which reflects the uniform radial profile of the RE distribution function. At a later time, the spectrum for the unscattered distribution changes little, whereas for the scattered distribution the ECE signals grow at all frequencies and the spectrum becomes flatter, meaning that REs give larger contributions to ECE at higher frequency ECE than the lower frequencies. This flattening behavior is consistent with the ECE spectrum observed in experiments<cit.> (Fig. <ref>).





Fig. <ref> shows the X-mode ECE signals at two frequencies in stage 2, with and without the effects of kinetic instabilities. The initial states of f are both from the scattered distribution functions in stage 1. When the kinetic instabilities are turned off, the ECE signals drop to its base values very quickly (in 0.2s), and then stay at these low levels, despite the fact that there remains a large population of high-energy REs. One the other hand, with the wave diffusion turned on, the ECE signals stay at a high level and drop gradually, which is consistent with the evolution of RE density and agrees with experiments. This further confirms that the high level of non-thermal ECE signals observed is a reflection of not only the large population of runaway electrons generated, but also the intensity of the excited kinetic instabilities.



In addition to the trends of ECE signals within the discharge, we also study the behavior of the ECE signals on short timescales. Examination of the simulation results shows that the X-mode ECE signals in stage 2 has oscillatory behaviors, as shown in Fig. <ref>, with periods of about 0.002s. In each period, the ECE signals experiences a fast growth phase, followed by a slower decaying phase. This behavior is consistent with the "inverse-sawtooth" behavior of ECE signals observed in QRE experiments<cit.> and in post-disruption scenarios<cit.>.





§ EXPLANATIONS OF SIMULATION RESULTS



To have a better understanding of the results obtained in Sec. <ref>, here we show the details of the calculation using the reciprocity method.

We first look at the results of the emission and absorption profiles. Fig. <ref> shows the profiles of k_i (Eq. (<ref>)) and U/T_e0 (Eq. (<ref>)) in the plasma (from R=1.5m to 2.5m) in stage 1, for X-mode ECWs with ω=2ω_ce0. At the beginning of stage 1 when the electron distribution is close to a Maxwellian (Fig. <ref> (a)), both the emission and the absorption are strongly localized near the magnetic axis, where the resonance region for this wave is located. The value of U/T_e0 is almost equal to k_i in the resonance region, consistent with Eq. (<ref>).



Fig. <ref> (b) shows the emission and absorption profiles after a significant population of RE generated, including effects of wave-particle interactions (Fig. <ref> (b)). We can see that the RE tail increases the amplitudes of emission and absorption outside the resonance region, as discussed in Sec. <ref>. As shown by the orange and blue lines outside the resonance region, the ratio between the emission and absorption from RE tail is not equal to T_e0 but much larger. This means that REs contribute to strong emission of ECWs but weaker absorption, compared to thermal electrons. This difference originates from the shapes of distribution function f. Note that in Eqs. (<ref>) and (<ref>), the emission coefficient depends on the value of f satisfying the resonance condition, whereas the absorption coefficient depends on the gradients of f. As shown Fig. <ref> and Fig. <ref>, the RE tail distribution has a much smoother shape compared to the Maxwellian part, resulting in a weaker absorption for the ECWs. The wave-particle interaction can further broaden the difference by making the RE tail distribution smoother.
The result of artificial wave calculation, which only depends on the absorption, is almost identical in the two cases, as shown by the green line in Fig. <ref>.


Fig. <ref> shows the emission and absorption profiles for X-mode ECW with ω=3ω_ce0, and the O-mode ECWs with ω=2ω_ce0, both calculated from the RE tail distribution. We can see that for ECWs with higher frequency or different polarization, both the absorption and the emission in the resonance region is much weaker compared to Fig. <ref>. This leads to a difference in the optical thickness. As shown by the green lines in Fig. <ref> and Fig. <ref>, the damping of the artificial waves are very different in the resonance regions.
Unlike the thermal electrons, the level of ECE from runaway electrons outside resonance region (the orange lines) at different  frequencies and polarization are close, meaning that the RE tail takes a more important role in the ECE radiation power. This explains that the ECE spectrum will become flatter when the contributions from RE dominate. It also illustrate that the growth of O-wave ECE found in Sec. <ref> is not only due to the polarization scrambling at walls, but also comes from the enhanced emission of O-wave from the runaway electrons.
The weak damping of ω=3ω_ce0 wave and the O-wave means that the plasma cannot be considered as a black body and the received ECE power will thus be smaller than T_e0 even for a Maxwellian distribution. In addition, it also means that the emission of ECWs from REs in the whole domain can propagate to the receiver, whereas the for the X-wave with ω=2ω_ce0, the emission from the low-field-side will be largely absorbed in the resonance region. Thus the ECE power for ω=3ω_ce0 can be larger than 2ω_ce if the emission from RE dominates.



To illustrate the dependence of the ECE power on runaway electrons at different momentum, we introduce the weight function for ECE radiation. Note that in Eq. (<ref>), the absorption of the artificial wave k_i mainly depends on the thermal electrons and the wall reflection, since the absorptions by REs are weak. The emission term U depends on the both the Maxwellian part and the RE tail through Eqs. (<ref>) and (<ref>). In terms of that, we can substitute a delta function δ(𝐩-𝐩_0) for f into Eq. (<ref>) to obtain U, and  calculate the ECE weight function W(𝐩_0) with k_i from a Maxwellian distribution. The radiation effective temperature can then be expressed as

    T_eff=∫ d^3𝐩_0 W(𝐩_0) f(𝐩_0).


Fig. <ref> shows the values of W(𝐩_0) and W_𝐩_0f(𝐩_0) in RE momentum space, for X-mode ECWs with ω=3ω_ce0 and a scattered runaway electron distribution. We can see that only runaway electrons with large pitch angles can give significant contributions to ECE power. As shown in Fig. <ref> (b), for the scattered RE distribution the ECE power comes from both thermal electrons (the small circular region near p_∥=p_⊥=0) and the runaway electrons with p<5 and pitch angle larger than π/6. The low-energy runaway electrons with large pitch angles are a result of the wave-particle interaction from the high frequency whistler waves<cit.>. This explains why the unscattered RE distribution gives little ECE radiation power, since they have little population in this region (Fig. <ref> (a)).



According to the ECE weight function, the ECE signals mainly characterize the RE density in the low-energy regime (<2.5MeV), which is different from the SE and gamma ray radiation which mainly comes from very high-energy runaway electrons. The ECE power is also very sensitive to the pitch angle distribution. Given the shape of the weight function, a high level of ECE from RE can be regarded as a signal of anomalous pitch angle scattering of REs in the low-energy regime, and an indication of RE kinetic instabilities.

The oscillatory behavior of the ECE signals shown in Fig. <ref> is related to the nonlinear dynamics of the wave-particle interactions. This behavior has been observed in DIII-D QRE experiments, and is found to be correlated with the excitation and damping of whistler waves with frequency 100-200MHz<cit.>. As shown by the weight function in Fig. <ref>, the high level ECE signals are sensitive to the population of REs with large pitch angles. Given that the excitation of whistler waves requires a strong gradient of f, and the quasilinear diffusion tends to  smoothen the gradient, the oscillatory behavior of ECE represents the scatter-replenish cycle of the resonant electrons the kinetic forces and wave-particle interactions, and the unstable-stable cycle of the whistler waves. This nonlinear behavior of the whistler waves and RE distribution will be studied in future.



§ SUMMARY



ECE is a powerful tool to diagnose the runaway electron dynamics in tokamaks, but the interpretation of the signal is much more complicated than the black body case, and a synthetic diagnostic model is required. In this paper an ECE synthetic diagnostic model for arbitrary electron distributions is presented, and the conditions of runaway electrons producing significant non-thermal ECE radiation are discussed. According to the shape of the RE tail distribution, it was found that runaway electrons can give remarkable emissions of ECWs, but little contribution to the absorption. In addition, the emission power depends on both the energy and the pitch angle distribution of RE. Using the experimental parameters from DIII-D low-density flattop scenarios, we find that the RE tail can give strong non-thermal ECE radiation in both X-mode and O-mode, and make the ECE spectrum flatter. These results by the first time give solid explanations of these observations in RE experiments, and have good agreements with DIII-D diagnostic results. Further analysis of the weight function of ECE suggests that most of the high frequency ECE comes from REs in the low-energy regime with large pitch angles, a unique feature from other diagnostics. In addition to runaway electron studies, the synthetic diagnostic model can also be used to study ECE signals from other scenarios such as cases with strong current drive through wave particle interaction.


The kinetic instabilities associated with the non-monotonic and anisotropic RE beam play an essential role in producing the non-thermal ECE radiation.
In our simulations of the QRE experiments, we find that the ECE power is weak unless kinetic instabilities are excited and scatter the low-energy runaway electrons. In addition, the ECE from the unscatterd RE tail, with scattering from kinetic instabilities turned off, is not much different from a Maxwellian distribution, and only the scattered RE tails show the interesting properties of non-thermal ECE described above. Therefore, the significant non-thermal ECE radiation, which is ubiquitous in RE experiments, can be used as an indication of RE kinetic instabilities, which can be applied to other tokamak devices and ITER. The scattering effects from kinetic instabilities can also be important in the study of runaway electrons generated in disruptions, especially in ITER.

The model can be further improved in several ways. In the calculation of ECW emissions and absorptions, we use a uniform profile for the RE distribution near the core. This part can be improved by extending the kinetic simulation model from 2D to 3D, and applying a bounce-average kinetic simulation code like CQL3D or LUKE. Furthermore, in the calculation of ECE power using the reciprocal method, we only trace a single ray along the major radius direction, and ignore  the diffraction effects. This shortcoming can be improved by using an ECEI model incorporated in a recently-developed Synthetic Diagnostic Platform (SDP)<cit.>, which calculates ray tracing within a 2D wavefront plane including the refraction and diffraction of the wave.



Chang Liu wants to thank Robert Harvey, Ernest Valeo, John Krommes, Donald Spong, William Heidbrink and Yong Liu for fruitful discussions. This work has received funding from
the Department of Energy 
under Grant No. DE-SC0016268 and DE-AC02-09CH11466.

iopart-num

		
