Cold-hot coupled waves in a flowing magnetized plasma

Nonlinear coupling of cold and hot waves in a flowing magnetized plasma is analyzed with the Vlasov equation. An analytical solution is obtained for cold waves of a small amplitude (weak flow) and a long wavelength. The distribution function is obtained by integrating the kinetic equation along a perturbed phase-space trajectory for a time-varying plasma flow. The kinetic description presents a generalized dispersion relation that involves resonances depending on cold and hot wave dispersions. Coherent fluid motion leads to radiation peaks in addition to the cyclotron harmonics, where the wavenumber of the cold wave determines the peak frequencies. The peaks appear narrow when the wave propagates perpendicular to the time-averaged flow while they become broad due to the Doppler effect when the wave propagates parallel to the flow. Fully kinetic particle-in-cell simulations corroborate the theoretical predictions. The dispersion relation and resulting wave spectra provide information about plasma parameters and flow properties.


Introduction
Rapid convective transport events are frequently observed in a wide variety of magnetically confined plasmas. In tokamaks, magnetohydrodynamic (MHD) instability modes often collapse with an abrupt transport of energy and particles. Notable examples are sawtooth crashes [1,2], the formation of tearing modes [3,4], and the burst of edge-localized modes (ELMs) [5][6][7]. On a larger scale, other examples include the ionospheric flows [8] and the sporadic explosion of the coronal loops on the surface of the Sun [9]. These events involve a time-dependent plasma flow of particles, which drives electrical currents at least in some localized regions Original Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. albeit not globally. The time-varying currents will induce electromagnetic (EM) waves via Maxwell's equations, which makes relevant the analysis of EM waves for the understanding of the explosive transport dynamics.
In this paper, we report the wave dynamics in magnetically confined plasmas with a macroscopic flow. As a basis for establishing the wave model, we adopt the cold wave analysis to describe the macroscopic flow and analyze the excitation of hot waves from the flow. When a plasma flow is induced in magnetically confined plasma, the macroscopic flow of charged particles can oscillate as a form of cold wave [10]. The frequency range of the cold wave is determined by the wavenumber of the stream. On the other hand, the hot waves are manifestations of the interactions among the individual particles and the EM fields [11], where the wave dispersion relation presents the cyclotron harmonic spectrum. Besides the cold and hot waves, nonlinearly coupled harmonic waves can also be generated via cold-hot wave interactions. Our analysis of cold and hot waves induced by macroscopic flow leads to a generalized dispersion relation of cold, hot, and coupled waves, bearing information about the fundamental plasma parameters such as the density, flow velocity, temperature, and electric and magnetic fields.
The EM waves associated with time-varying plasma flow may complement the long-standing interests in the wave instabilities driven by minority energetic particles as the primary energy source [12,13]. Modeling and experimental investigations on the magnetoacoustic cyclotron instability showed that fusion-born alpha-particles can excite superthermal ion cyclotron emissions [14,15]. As the intensity of the ion cyclotron harmonics depends on fusion reactivity and local plasma parameters [16], the radiation power measurement has been proposed for the diagnostics in future fusion devices such as ITER [17] and DEMO [18]. For the electron cyclotron frequency range, the anomalous Doppler instability driven by superthermal electrons was suggested to trigger the radiative relaxation during the ELM crash in tokamaks [19,20]. The electron cyclotron waves observed in planetary magnetospheres were interpreted as a result of the injection of energetic plasma particles [21][22][23]. These studies show that the kinetic waves are proxy to the dynamics of energetic particles. Similarly, the EM waves can carry essential information on the convective transport dynamics in magnetized plasmas.
This paper is structured as follows. In section 2 the dispersion relation is obtained by solving the Vlasov equation with the method of characteristics [24][25][26]. The waves propagating perpendicular to the magnetic field are investigated to understand the coupled wave harmonics. In section 3 the analytic results are corroborated by the particle-in-cell (PIC) simulation [27]. The numerical analysis not only demonstrates the analytical dispersion relation but also shows the behavior of wave intensity. The appendix provides the kinetic dispersion relation of the shifted bi-Maxwellian distribution for inhomogeneous plasma.

Overview
For the wave analysis, a distribution function f s (x, v, t) of a charged species s is decomposed into the zeroth order f s0 and the small deviation f s1 (f s0 ≫ f s1 ), Defining the lowest order flow velocity from f s0 Then the number density, flow velocity, and current density become to the first order where q s is the charge, The goal of this paper is to study the wave dynamics of J 1 and EM fields deriving the dispersion relation of a flowing plasma. Therefore f s0 in this work will be a general Maxwellian distribution depending on the zeroth order fields, while f s1 depends not only on f s0 but also on a non-Maxwellian distribution and the first order fields. Figure 1 illustrates an example of uniform time-varying flow in a Maxwellian plasma. A single electron velocity in figure 1(a) is given where v Te is the thermal speed and ω ce is the cyclotron frequency. Time evolution of the normal- is shown in figure 1(b).
To obtain J 1 in equation (9) we solve the Vlasov equation for f s1 with f s0 given using the method of characteristics [25]. The dispersion relation of a plasma wave is obtained by combining Faraday's and Ampère's laws in wavenumber space [28] nn · E − n 2 E + K · E = 0, where n = ck/ω is the refractive index, c is the speed of light, k is the wavenumber, and ω is the wave frequency. The effective dielectric tensor is defined as where σ is the conductivity tensor defined by J = σ · E and ε 0 is the vacuum permittivity.

Particle trajectory in magnetized flow
The current density J 1 is given by f s1 in equation (9) and f s1 can be obtained by solving the linearized Vlasov equation with where m s is the mass. Integrating the above equation along the characteristic curve gives , and v(t ′ ) are understood in the integrand. As seen in equation (13) the integration evaluates the change of f s1 along the zero-order trajectory x(t ′ ), v(t ′ ) determined by E 0 , B 0 [10,11]. The characteristic curve where the boundary conditions are x(t) = x and v(t) = v at t ′ = t. When solving this equation we decompose E 0 =Ē 0 + E 0 and B 0 =B 0 +B 0 . In our analysis, a bar denotes a constant field and a tilde denotes a fluctuating field of a cold wave.
When we perform the integration in equation (14) we . This can be justified by an approximation |k · ∆x| ≪ |ω∆t|, i.e. |k ·ũ s0 | ≪ |ω| in theū * s0 frame when the majority of particles are distributed near the mean velocity, so that most of the contribution to the integral comes from the vicinity ofũ s0 . Therefore the flow velocity of a cold wave will be approximate toũ s0 ∝ exp (−iωt ′ ) in our analysis. For multiple modes, we writeũ where the lth modeũ has the amplitudeũ ′ s0lr , the frequencyω l , and the phase shift ϕ slr in the direction ofr (r = x, y, z). Equation (22) describes not only the cold waves but also more general time-varying macroscopic motions satisfying |k ·ũ s0l | ≪ |ω l (k)| (e.g. externally driven waves, hot plasma waves, waves in inhomogeneous plasmas). Finally, the total particle velocity v(t ′ ) =ū * s0 +ũ is determined by equations (16)- (22) and can be used to evaluate f s1 in equation (14).
The classical analysis of the kinetic wave [10,29] used the unperturbed trajectory v =ū * s0 + v * withoutũ s0 . However, the kinetic analysis in this paper establishes the hybrid dynamical model that incorporates both the continuum and particle mechanics. Consideration of the flow-driven perturbation extends wave analysis to more general situations where time-varying plasma transport exists.

Dispersion relation
In this section we consider a uniform anisotropic plasma. A non-uniform plasma will be considered in the appendix. 2.3.1. General form. The distribution function of a uniform anisotropic plasma is described by a shifted bi-Maxwellian, where v Ts⊥ = (2T s⊥ /m s ) 1/2 and v Ts∥ = (2T s∥ /m s ) 1/2 (the subscript 0 will be omitted) are respectively the perpendicular and parallel thermal speeds with respect toB 0 . Note that the Maxwellian distribution function (equation (23)) and the cold wave motion (equations (19)- (22)) exactly satisfy the Vlasov equation. The integrals of the distribution function over the particle trajectory (equation (14)) and the velocity (equation (9)) are summarized in the appendix. As a consequence, the dielectric tensor in a uniform plasma can be expressed as where here Λ s = k 2 ⊥ r 2 Ls represents the finite Larmor radius effect, r Ls = v Ts⊥ / √ 2|ω cs | is the Larmor radius, J m lr = J m lr (k rũ ′ s0lr /ω l ) and I n = I n (Λ s ) are respectively the Bessel and the modified Bessel functions of the first kind, and I ′ n = ∂I n (Λ s )/∂Λ s . The notation In equation (25) the subscripts a and b can be either x, y, or z and Z = Z(ζ sm lr n ) with The plasma dispersion function Z(ζ) is defined as and it follows that Z ′ = dZ/dζ = −2(1 + ζZ). The variable u χs is given by which is the sum ofū s0 and the hybrid phase velocities ω l /k rr that depend on the cold wave frequencyω l and the hot wavenumber k. Note that u ds vanishes (ū s0 =ū * s0 ) for a uniform plasma in equation (29). The effects of nonvanishing u ds are considered for a non-uniform plasma in the appendix. The Z(ζ sm lr n ) terms in equation (25) show the resonances at where ω D = ω − k ·ū * s0 is the Doppler-shifted frequency.

Perpendicular propagation.
Since the hot wave harmonics are originated from the cyclotron motion, the essence of the cold-hot coupled harmonics can be understood by analyzing perpendicular propagation acrossB 0 (k ∥ → 0, |ζ sm lr n | ≫ 1). For convenience, let x be the direction of the perpendicular propagation, i.e. k ⊥ → k x . For ignorable parallel flows u s0∥ , K xz , K yz , K zx , and K zy can be ignored and then the dielectric tensor components in equations (24)-(25) reduce to where we have defined Φ = m l ϕ slx − (k xũ ′ s0lx /ω l ) sin ϕ slx and omitted the subscript x from m l because m ly and m lz terms vanish as k y , k ∥ → 0. The diamagnetic effect has been ignored, i.e. u s0 = u * s0 . The imaginary part of Z(ζ sm l n ), , has been neglected because of ζ sm lr n ≫ 1 and we examine only the stable hot waves using the real part of the dispersion. With D (ω, k) defined to be the determinant of the operator nn − n 2 I + K in equation (10), the dispersion relation for the perpendicular wave is Then the dispersion relation becomes, for the (E 1x , E 1y , 0) mode, and, for the (0, 0, E 1∥ ), Note that equations (33) and (34), respectively, describe extraordinary and ordinary hot waves similar to the cold waves of equation (20). While a constant-flow plasma has resonances at ω D = ±nω cs , the time-varying-flow plasma has resonances at ω D = ±nω cs + ∑ l m lωl due to the cold-hot coupled harmonics.
In the electrostatic limit (E 1x ≫ E 1y , E 1∥ ), the dispersion relation of equation (33) reduces to K xx = 0, that is, or, for smallũ s0 (t), This equation expresses the dispersion relation for the hot electrostatic wave in a flowing plasma. For a longitudinally static plasmaū s0x = 0, by replacing ω D → ω, the dispersion relation describes the Bernstein wave [28,29]. The harmonic resonance in a flowing plasma (equation (30)) is affected by the wavenumber because of the Doppler shift k ·ū * s0 . It suggests, by invoking the inverse Fourier transform (e. g. E(ω, k) → E(ω, x)), that the frequency spectrum of the EM field at a location x depends on k ·ū * s0 . Since the waves with different wavenumbers resonate at different shifted frequencies, the spectrum of ω(k) becomes broader when the magnitude of k ·ū * s0 increases. This analysis of the Doppler broadening should be applied to the spectral interpretation of plasma transport.

Example study of perpendicular propagation.
Now we exemplify the dispersion relation (33) of the extraordinary wave for the EM wave. Consider a deuterium plasma: m i /m e = 3670, n e0 = n i0 = 10 19 m −3 , T e = T i = 1keV (isotropic), and B 0 = 2ẑ T. The dispersion relation is obtained from equations (31) and (33) where k xūs0y is ignored for simplicity. For the infinite sums in equation (31), 1 ≤ n ≤ 100 and −100 ≤ m ≤ 100 are considered for convergent results. The cold wave dispersion relation is given by equation (33) with K replaced byK of equation (20) and the electrostatic wave dispersion is obtained from equation (36). Figure 2 presents both the cold and hot wave dispersion relations in the Dopplershifted frequency domain for the constant flow velocity. In figure 2, (a) shows both the cold and hot wave branches of the EM wave while (b) does not show the cold wave in the electrostatic limit, as expected. The top panel of figure  2(b) also shows the electron Bernstein wave whose frequency approaches the upper hybrid frequency ω uh = (ω 2 pe + ω 2 ce ) 1/2 as k x r Le → 0 [30]. Figure 3 shows the detail of figure 2 for cold and hot branches and more details in (a) and (b) at the corresponding points of the left panel.
Instead of directly solving D X (ω, k x ) = 0 with theũ s0lr terms, we obtain the solutions by adding the frequencies of 0.3mω ce and the roots of the time-independent relation (the top plot of figure 2(a)). Figure 5 shows that the approximate results agree well with the PIC simulation. Note that the negative roots ω(k x ) < 0 of the time-independent dispersion relation has been used since they become positive by the shift 0.3mω ce . For visibility, the analytical results are only depicted for |m| ≤ 3 in the left panel and but for |m| ≤ 9 in the zoomedin view of the right panel. The details of PIC simulation and the numerical verification of the cold-hot coupled harmonics are discussed in the next section.

PIC simulation
The analytic description of the cold-hot coupled wave can be verified by the 1D3V (one spatial dimension, three velocity dimensions) PIC simulation. We use the EPOCH code [27] to demonstrate the cold-hot coupled wave. The numerical convergence has been investigated and confirmed. The accuracy is ensured by setting the numerical setup based on the plasma parameters. The spatial grid is 0.05r Le , which is less than 4% of the Debye length. The total length is of the order of 10 3 r Le , such that several numbers of cold wavelength 2π/k x satisfy the periodic boundary condition in the x direction. The time step is restricted by the Courant-Friedrichs-Lewy condition. The number of superparticles per cell is 10 4 (corresponding to 6 × 10 5 per Debye length).
The particle distribution and the cold wave frequency have been carefully chosen to satisfy the zeroth order Vlasov  (23)) with the isotropic temperature. To generate a cold mode of wavenumberk x , we start the simulation with the same velocity profile u init = u ′ init cosk x x for electron and ion with E init = 0. Note that the same velocity profile has been used for a selfconsistent E × B drift.
To describe the simplest situation, we first consider a uniform flow, i.e.k x = 0. The initial flow is imposed to be u init = u ′ initŷ in figure 6(a) and u init = u ′ initx in figure 6(b), where u ′ init = 0.02v Te . For each case, the temporal Fourier transform of the cold wave simulation with E init = 0 (red solid line) is compared to that of the simulation with E init = B init × u init (black dotted line) which is chosen for constant E × B drift. In each figure 6(a) and (b), the top panel shows v x of a single particle with the initial velocity v init = u ′ init (x +ŷ). The middle panel plots the flow velocity u ex . In figure 6(b) the spectral peaks at ω = 0 indicate the existence of a constant E × B drift even for E init = 0 case. The initial flow can generate not only the perturbed fieldẼ 0 but also the constant fieldĒ 0 that satisfies u init ≃Ē 0 × B init /B 2 init [31,32]. The fieldĒ 0x generated by u init = u initŷ can be seen in the bottom panel of figure 6(a). The two peaks near the electron cyclotron frequency correspond to the cold wave frequencies 0.2ω ce and 1.2ω ce for k x r Le ≪ 1 in figure 3(a). The cold-hot coupled harmonics in the vicinity of ω = (n + 0.2m 1 + 1.2m 2 ) ω ce are observed in the u ex and E x plots in 6(a) but not in (b), because those peaks of different  wavenumbers overlap due to the Doppler effect in (b) wherē u e0 (=ū * e0 in a uniform plasma) is parallel to the wavenumber. Note that the cold-hot coupled wave still exists (see the bottom panel of figure 7(b)). Figure 7 shows the E x spectra in ω-k x space obtained from the simulations. The columns (a) and (b), respectively, correspond to the columns (a) and (b) of figure 6. The upper and lower plots present the hot wave dispersions without and with cold wave coupling respectively. The spectra clearly show the dispersion relations. The E y and B ∥ spectra show the same results. One can also see that the upper plots reproduce the same dispersions with correspondingū e0 in figure 2(a) top. In addition to the dispersions shown in the top panels, one can see the cold-hot coupled harmonics in the lower plots. Note that the Doppler shift is observed in figure 7(b) asū e0 is parallel to the wavenumber. This shift causes the overlap of spectral peaks and hence produces the broad frequency spectrum as in figure 6(b). So far we have analyzed a single-mode cold wave. Now we study a superposition of multiple waves. As an example study we consider four wave numbersk x r Le = 10 −3 , 0.01, 0.02, and 0.04 with E init = 0 and u ′ init = 0.02v Te . The corresponding cold wave frequencies in the electron cyclotron frequency range are, respectively,ω/ω ce = {0. 21 figure 2. The result for each wavenumber is presented in figure 8(a). Note that the cold-hot coupled harmonics forω/ω ce = 0.30 in figure 5 is clearly shown in the second panel. Figure 8(b) shows that the multiple cold waves excite the coupled harmonics at ω ≃ nω ce + ∑ l m lωl as shown by the kinetic analysis. For the multi-mode cold wave simulation the flow amplitude is reduced to u ′ init = 0.02v Te /4 to obtain the same order of amplitude for E x . The result also shows that the combination of four modes produces the broad spectrum and increase the EM field energy. In general, both the Doppler shift and the cold-hot coupled harmonics can broaden the frequency spectrum, and the  time-dependent flow-driven waves can increase the EM wave energy in the plasma. Figure 9 shows the oscillatory components of field and particle energies. The energy density is calculated as w Er = ε 0 E 2 r /2 for the electric field, w Br = B 2 r /2µ 0 for the magnetic field, and w Ks = ∑ m s v 2 /2 for the species s, and averaged over the spatial domain after removing the constant component. Figure 9(b)-(d) shows that the combination of the induced cold wave and the intrinsic hot wave (seen in figure  9(a)) drives the coupled waves shown in figures 5-8. For the performed setup the simulations do not exhibit a noticeable instability. A complete study of instability can be performed by thorough scanning of complex number space of frequency.

Summary and discussion
We have examined the generation of cold-hot coupled waves in a flowing plasma using the Vlasov formalism and the PIC simulation. To facilitate an analytic study on the relation between the flow and waves, we developed an astute kinetic model, which incorporates the time-varying macroscopic flow and microscopic perturbations. In particular, by representing the time-varying macroscopic flow as a cold wave, the kinetic model yielded the wave dispersion relations containing the cold waves, the hot waves, and their coupling.
We obtained an analytically manageable expression for the cold-hot wave dispersion relation for the particular case where the wave propagation is perpendicular to the background magnetic field, which is corroborated by PIC simulations. The PIC simulations confirm the analytic result that a narrow peak in the frequency spectrum appears at each coupled harmonic for the wave propagation perpendicular to the timeaveraged flow. In contrast, the Doppler effect results in a broad frequency spectrum for the wave propagation with a finite parallel component to the time-averaged flow. Furthermore, the simulations show that the time-dependent part of the plasma flow broadens the spectrum and increases the EM field intensity via the generation of the cold-hot coupled waves.
The analytic dispersion relation and the PIC simulation results based on our kinetic model illuminate the essential physics on the wave dynamics in a flowing magnetized plasma. This will improve our understanding of the transport processes in the laboratory and magnetospheric plasmas which are associated with the MHD instability [33,34], magnetic reconnection [35][36][37], turbulence [38], and the time-varying electric field [39]. The analysis of waves in a variety of time-varying flow may provide a theoretical background for plasma diagnostics, for example, the EM wave diagnostics [40] are useful tools for studying MHD-induced transport events and energetic particle physics.
Although we have investigated the specific cases where both the flow and wave propagation are perpendicular to the background magnetic fields, the method developed here can be extended to more general situations where the plasma  . Energy densities of electric field w E , magnetic field w B , and particles w K . The simulations correspond to (a) the reference (u init = 0, E init = 0) and the cold wave simulations (b)kx = 0 (figure 6(a)), (c)kx = 10 −2 /r Le (figure 8(a)), and (d) multi-mode ( figure  8(b)). Total energy conservation was confirmed.

Appendix A.
The distribution function with non-uniform density can be written as where u ds = −T s⊥ ∇n s0 ×B 0 /q s n s0B 2 0 and |k r n s0 /(∂n s0 /∂r)| ≫ 1 is assumed for the Fourier analysis. By using equations (16) and (22) in the particle trajectory v(t ′ ) = u * s0 (t ′ ) + v * (t ′ ) and x(t ′ ) =´v(t ′ )dt ′ with the boundary conditions v(t) = v and x(t) = x, the space-time dependence exp (ik · x(t ′ ) − iωt ′ ) of E 1 , B 1 in equation (14) becomes The exponential factor having the summations over l and r can be recast using the Bessel identity e iz sin θ = where J m lr = J m lr (k rũ ′ s0lr /ω l ) and the infinite series ∑ ∞ {m lr }=−∞ is define in equation (26). Inserting equations (A1)-(A4) into the linearized Vlasov equation (14) and then using Faraday's law B 1 = k × E 1 /ω gives where A x = a x cos ω cs τ − a y sin ω cs τ, A y = a y cos ω cs τ + a x sin ω cs τ, and, from now on, E 1 , B 1 ∝ exp (ik · x − iωt ′ ). In equation (A5), the approximation has been made that u ds ≪ v Ts⊥ so exp in ∂f s0 /∂v. The contribution of f s1 (x, v, t 0 ) when t 0 → −∞ is neglected [10]. We insert equation (A5) into equation (9), integrate over v = u * s0 + v * , and use equation (11) to obtain the general dielectric tensor where κ is given by equation (25).
where I n = I n (Λ s ), Z = Z(ζ sm lr n ), u * χs = u χs − u ds , ω ds = k · u ds is the drift frequency, and the subscripts a and b can be either x, y, or z. Equation (A7) together with the wave equation (10) provides a generalized dispersion relation of the collisionless Maxwellian plasma.
The form of equation (34) can include the dispersion relations for waves in an unmagnetized plasma by setting Λ s = 0. This is because the electric field is parallel to the background magnetic field. The dispersion relations for EM and electrostatic waves respectively are obtained by letting k ∥ → 0 and k ⊥ → 0.