Particle-in-cell simulations of parametric decay instabilities near the upper hybrid layer

. The particle-in-cell (PIC) code EPOCH is used to simulate parametric decay instabilities (PDIs) converting a 105 GHz microwave X-mode pump wave into electrostatic daughter waves at the upper hybrid (UH) layer of a fusion plasma in 1D. The resulting ﬁelds are analyzed, identifying modes in f - and k -space, and estimating their spectral power as a function of pump wave intensity. Both linearly and nonlinearly converted modes are identiﬁed and their characteristics agree with literature. A dipole approximation employed in literature appears to be unjust.


Introduction
Observations of strong scattering shifted in frequency in several recent collective Thomson scattering (CTS) experiments [1] has led to investigations to quantify the scattering and uncover the mechanisms that cause it.If not anticipated, the scattering can be strong enough to damage microwave equipment in current reactors such as ASDEX Upgrade.CTS diagnostics are expected to be installed in ITER and a better understanding of the anomalous scattering is therefore necessary.The observed scattering occurs when the CTS gyrotron beam reaches the UH layer and it has characteristics associated with PDIs.An investigation of the daughter waves of such a PDI process may lead to efficient heating schemes and there is experimental evidence to suggest that it is possible to suppress the PDI by angling the CTS microwave beam in certain ways.Better knowledge of this PDI is expected to also benefit investigations of other PDIs such as two plasmon decay during 2nd harmonic ECRH, which may have a significant impact on its heating efficiency.Using a numerical code to simulate a warm fusion plasma where the anomalous scattering is expected to occur is a convenient tool to characterize the scattered modes on a number of parameters, not only in frequency.The PDIs of interest are expected to also rely on both wave and particle dynamics so a PIC code should be able recreate the observed scattering and test theoretical models describing the scattering as the product of a PDI.

Parametric decay theory
The PDIs of interest are described theoretically in [2].They are three wave interactions where a large amplitude * e-mail: mgse@fysik.dtu.dkelectromagnetic X-mode pump wave, of f 0 and k 0 , decays into a low frequency warm lower hybrid (LH) daughter wave, f 1 and k 1 , and a slightly downshifted electron Bernstein wave (EBW), f 2 and k 2 .Subsequent recombination of the pump wave and the low frequency daughter produces an upshifted daughter wave, f 3 and k 3 .These processes must conserve energy and momentum which translates to conservation of frequency and wave vector, i.e.
the decay and combination process respectively.Both f 2 and f 3 are then shifted by f 1 relative to f 0 .Using a dipole approximation, the electromagnetic X-mode pump wave is shown to decay into two electrostatic daughter waves which for perpendicular propagation to the magnetic field in 1D are described by the dispersion relations Here ω ce = eB/m e , ω pe = n e e 2 /( 0 m e ), ω pi = n i Z 2 i e 2 /( 0 m i ), ω UH = ω 2 ce + ω 2 pe and ω LH = |ω ce ω pi /ω UH | are the electron cyclotron, electron plasma, ion plasma, upper hybrid (UH) and LH angular frequency respectively, T i , T e are ion and electron temperature in electron volts, Z i is the ion atomic number, r Le = v te /( √ 2|ω ce |) is the electron Larmor radius, and v te = √ 2T e /m e is the thermal velocity.These waves are referred to as warm LH and UH waves as those are their frequencies in the cold limit.Above a pump wave amplitude threshold, the processes become unstable and increase nonlinearly with pump wave amplitude.The wave enhancement of an X-mode wave incident on the UH layer in a magnetically confined fusion plasma means that the threshold may be exceeded at CTS diagnostic power levels.This conversion process of the pump wave, however, is not the only one that is expected to occur near the UH layer.As the X-mode pump wave moves closer to the UH layer it gets increasingly electrostatic.Eventually coinciding with an EBW, the pump wave is linearly converted from X-mode to EBW [4] and is reflected.

Numerical Setup
The simulations are run using the PIC code EPOCH [3] version 4.9.0 and plasma parameters are chosen to resemble the region around the UH layer of a warm deuterium plasma.The numerical domain is a 1D domain in the xdirection, split into three regions.The left and right regions are completely homogeneous for easier mode identification and the middle region bridges the parameters of the other two regions.The background magnetic field and temperature are homogeneous and given by B 0 = B 0y = 3.35 T and T e = T i = 0.3 keV.The plasma is a deuterium plasma with a density profile, n(x) = n e (x) = n i (x), given by where x l = 4.00 cm, x r = 4.35 cm.An f 0 = 105 GHz beam is propagated from the left boundary located at x = 0 using a pseudo-current at the boundary.Particles leaving at the boundaries are replaced by new thermally distributed ones.It is currently only possible to excite a beam with transverse polarization in EPOCH.An X-mode beam is desired so the beam is polarized in the z-direction which then gives rise to an x-component also.Note that the Xmode beam is injected from the high density side as it is unable to propagate in the region between the UH layer and the R-cutoff.This region is found on the right side of the UH layer with the used parameters and launching an X-mode beam from this side would therefore not work.It seems, however, that exciting the X-mode in the described manner in EPOCH also generates what appears to be an unwanted EBW at the same frequency in addition to the X-mode wave at the boundary.The EBW has a much greater wave number and propagates significantly slower.This wave is likely excited to compensate for the missing x-component when exciting the X-mode.In that case, the unwanted EBW can be decreased to some degree by starting at an even higher density at the expense of a larger computational domain.To avoid unwanted effects arising from significant discontinuities in the fields, the beam power is ramped up from zero using an arctanfunction with a characteristic time scale of τ ramp = 0.20 ns.The amplitude of the pump beam varies between simulations to investigate the nonlinearity of the PDI. Figure 1a shows a continuous wavelet transform (CWT) of E x into k-space, covering the last 1 cm of the first homogeneous region as well as the entire inhomogeneous region which includes the UH layer.As expected, the pump wave is seen to propagate to the UH layer where it is amplified and waves of different k are scattered back into the homogeneous region.Some waves also seem to be present beyond the UH layer.In the homogeneous region, k 0 is seen to be much smaller than that of the back scattered waves, however, near the UH layer, the wave number of the pump wave increases while the wave number of the back scattered waves decrease and become comparable.This suggests that the dipole approximation used in [2] might not be valid.With the dipole approximation, the wave number of the daughter waves are expected to differ only by sign, i.e. k 2 ≈ −k 1 .Without the dipole approximation, the selection rules change and the frequency shift might not agree with theory then.

Frequency space
An FFT into f-space near the pump wave is shown in figure 1b.A large peak at the pump frequency f 0 = 105 GHz is clearly visible and around this peak, f 1 = 0.78 ± 0.03 GHz from it, are two smaller peaks.Since they are found symmetrically around the pump frequency and shifted approximately by the lower hybrid frequency, f 1 ∼ f LH = 0.71 GHz, the smaller peaks are suspected to be the daughter waves at f 2 ≈ 104.22 ± 0.03 GHz and f 3 ≈ 105.78 ± 0.03 GHz.

2D FFT into frequency and k-space
Plots showing selected parts of a 2D FFT of E x into f -and k-space are shown in figure 2 (a)-(b).Figure 2 (a) shows the squared spectrum near the pump frequency.Symmetrically around k = 0, lines are observed resembling the dispersion relations of warm X-mode connecting to EBWs at the UH layer.A number of peaks can be spotted, one coinciding with the X-mode pump wave and two other at the same frequency but on the EBW branches.The left EBW wave has a positive group velocity and is the unintended one created at the boundary whereas the right wave has negative group velocity and is the one originating from around the UH layer.Around the peak and along the EBW line, a number of smaller peaks separated by a similar distance in f and k can be seen.The closest two are likely to  be the daughter waves also observed in figure 1b.Using the observed frequency shift to calculate k 2 from equation (2) places it on the EBW branch as expected.Figure 2 (b) shows the squared spectrum near the warm LH frequency.Using the same frequency shift to calculate k 1 from equation (1) places it close to albeit not quite at the observed peak.Judging from the peaks and the dispersion curves observed in both plots and comparing their wave numbers with figure 1a, the backscattered waves at much greater wave number than k 0 seen in figure 1a are com-  prised mainly of a linearly converted EBW and the warm UH daughter wave from the nonlinear PDI.

Group velocity
By considering v g = ∂ω/∂k along the dispersion curves in figure 2 (a), the linearly converted EBW and the warm UH daughter wave are expected to propagate at a similar speed and tracking the position of the backscattered wavefront in time is expected to give a reasonable estimate of the group velocity for the warm UH daughter wave.The wavefront is tracked in the first homogeneous region by fitting a shifted arctan-function to |E x | averaged over 60 gridpoints.Fitting a first order polynomial to the wavefront position in time gives the group velocity v g = −(7.51± 0.01) × 10 5 m/s.The expected group velocity from equation (36) in [2] is v g,theo = 10.66 × 10 5 m/s.A similar approach gives v g,pump = (5.962± 0.006) × 10 6 m/s for the pump wave.Note that the very small errors are calculated as standard deviations which do not take systematic errors into account.The group velocity of the back scattered waves is much closer to the theoretical group velocity of the UH daughter than to the incoming X-mode pump beam.

Pump intensity dependence
As a measure of the power contained in a particular mode, the | Ẽx | 2 spectra shown in figures 2 (a)-(b) are numerically integrated around selected modes in boxes of dimension 940 m −1 along the k-axis and 0.22 GHz along the f -axis.The intensity of the pump wave is calculated as , where ê is the unit vector pointing in the direction of the electric field and M is the Maxwell operator as defined in equations (5.138)-(5.139) in [6].Values are taken from the first homogeneous region of the domain and M is assumed to be cold as it is not right at a resonance.This is done for all simulations and the integrals are normalized to their value at the lowest pump intensity simulation to give a relative power.The found values are shown in figure 3. The normalized power of the EBW mode at the pump frequency appears to follow that of the pump wave and grows linearly with pump intensity, implying that it is indeed linearly converted at the UH layer.The suspected PDI products all grow faster than linearly and whilst it is difficult to determine the shape of the curves with only three points, they are perhaps getting close to a saturation level.

Discussion
The PIC code EPOCH successfully recreated scattering peaks shifted in frequency relative to a microwave pump, resembling peaks observed experimentally.It was found that the scattering increases with a nonlinear dependence on the pump wave intensity as a PDI would.Whilst theoretical predictions about the scattering characteristics seemed to agree to the order or better, implying that it is indeed caused by a PDI at the UH layer, the dipole approximation employed in [2] to derive equations ( 1)-( 2) does not seem to hold as seen in figure 1a.Tracking the backscattered wavefront to get the group velocity of the warm UH wave yielded a result on the same order as the theoretical prediction under the assumption that it propagates at a similar speed as the linearly converted EBW since the amplitude of the warm UH wave is likely insignificant in comparison to the linearly converted EBW.This might be better to do at greater intensities.The box integration method for obtaining a measure of power in the modes is not ideal as it cuts off the tails of the peaks and is affected in particular by neighboring modes, dispersion lines etc. Assuming, however, that the peak is the most significant contribution to the integral and that nothing else changes in its vicinity as the pump intensity changes, the method does the job and will provide a measure in which growth of the modes can be observed.More simulations are needed to establish a threshold pump intensity and to explore the possibility of a saturation level of the daughter waves.
The beam is expected to reach the UH layer in the middle region at x UH = 4.27 cm where n e (x UH ) = 2.87 × 10 19 m −3 .The entire domain is 5.35 cm, and contains n x = 6600 grid points and n part = 60000n x pseudoparticles, of which half are electrons and the other half are deuterons.The time step is 0.95 of what is allowed by the CFL criterion, the output time step is δt = (10.3f 0 ) −1 = 9.25 × 10 −13 s and the total simulation time is almost 46.5 ns of which the first 10 ns is disregarded in most of the analysis as this time is spent propagating the pump from the boundary to the UH layer.
(a) CWT of E x leading up to the UH resonance at t = 30 ns, using the Morlet wavelet as base function.The dotted white line indicates x = x l .The solid yellow line is the solution to warm X-mode propagation[5], and the dashed green line is the solution to the back scattered warm UH mode from equation(2).Squared FFT of a 31.75 ns long E x signal for an input intensity of I pump = 16.7 MW/m 2 .Only the spectrum near the pump wave frequency and the expected UH daughter wave is shown.The plot is centered around the pump at f 0 = 105.00GHz and the two vertical dotted lines placed symmetrically around the pump wave mark the two PDI products f 2 and f 3 .

Figure 1 :
Figure 1: Plots showing transforms of the longitudinal electric field component, E x , into k-space (a) and f -space (b).
(a) Spectrum near the pump frequency.(b) Spectrum near the lower hybrid frequency.

Figure 2 :
Figure 2: 2D fast Fourier transform of E x in t ∈ [10.0, 46.7] ns and x ∈ [2, 4] cm into f-and k-space.The white circles are the expected location for the pump beam (a), the warm UH daughter wave (a) and the warm LH daughter wave (b).