Particle-in-cell and weak turbulence simulations of plasma emission
Sang-Yun Lee, L. F. Ziebell, P. H. Yoon, R. Gaelzer, E. S. Lee
DDraft version November 7, 2018
Preprint typeset using L A TEX style AASTeX6 v. 1.0
PARTICLE-IN-CELL AND WEAK TURBULENCE SIMULATIONS OF PLASMA EMISSION
Sang-Yun Lee
School of Space Research, Kyung Hee University, Yongin, Gyeonggi 17104, Korea
L. F. Ziebell
Instituto de F´ısica, UFRGS, 91501-970 Porto Alegre, RS, Brazil
P. H. Yoon
School of Space Research, Kyung Hee University, Yongin, Gyeonggi 17104, KoreaInstitute for Physical Science and Technology, University of Maryland, College Park, MD 20742-2431, USAandKorea Astronomy and Space Science Institute, Daejeon 34055, Korea
R. Gaelzer
Instituto de F´ısica, UFRGS, 91501-970 Porto Alegre, RS, Brazil
E. S. Lee
School of Space Research, Kyung Hee University, Yongin, Gyeonggi 17104, Korea (Dated: November 7, 2018)
ABSTRACTThe plasma emission process, which is the mechanism for solar type II and type III radio burstsphenomena, is studied by means of particle-in-cell and weak turbulence simulation methods. Byplasma emission, it is meant as a loose description of a series of processes, starting from the solarflare associated electron beam exciting Langmuir and ion-acoustic turbulence, and subsequent partialconversion of beam energy into the radiation energy by nonlinear processes. Particle-in-cell (PIC)simulation is rigorous but the method is computationally intense, and it is difficult to diagnose theresults. Numerical solution of equations of weak turbulence (WT) theory, termed WT simulation,on the other hand, is efficient and naturally lends itself to diagnostics since various terms in theequation can be turned on or off. Nevertheless, WT theory is based upon a number of assumptions.It is, therefore, desirable to compare the two methods, which is carried out for the first time inthe present paper with numerical solutions of the complete set of equations of the WT theory andwith two-dimensional electromagnetic PIC simulation. Upon making quantitative comparisons it isfound that WT theory is largely valid, although some discrepancies are also found. The presentstudy also indicates that it requires large computational resources in order to accurately simulatethe radiation emission processes, especially for low electron beam speeds. Findings from the presentpaper thus imply that both methods may be useful for the study of solar radio emissions as they arecomplementary.
Keywords: methods: analytical – methods: numerical – plasmas – radiation processes: thermal –turbulence – waves [email protected] a r X i v : . [ a s t r o - ph . S R ] N ov INTRODUCTIONSolar radio bursts phenomena in the meter wavelengths were discovered and subsequently classified into five cate-gories, starting from the 1950s and are still being studied to this date (Wild & McCready 1950; Wild 1950a,b; Wildet al. 1954). The present work relates particularly to types II and III radio bursts. Both phenomena are related tohigh-energy beams of electrons produced during solar flares and Coronal Mass Ejections, which occur in the solarchromosphere and corona (Benz 2017), and are characterized by emissions in the VHF band, with frequencies rangingfrom 10 MHz to approximately 300 MHz (hence with wavelengths of the order of meters). Dynamic spectra of bothtypes II and III typically display negative frequency drift with time, but while the type III typically occurs in shortbursts lasting a few minutes, the type II is characterized by a much slower downward drift rate, which can last forseveral tens of minutes (Ergun et al. 1998; Bastian et al. 1998; Cane et al. 2002; Reiner et al. 2009; Hudson 2011; Reid& Ratcliffe 2014).The generally accepted mechanism by which types II and III solar radio bursts are generated was initially proposedby the pioneering work of Ginzburg & Zheleznyakov (1958) and has been greatly improved by several contributionssince then. According to the standard theory, disruptive events occurring in the solar chromosphere, such as solarflares, create bursts of energetic electrons that stream outwards from the Sun, along open magnetic fields lines. As thesebursts propagate through the corona, which already contains a denser population of particles in a quasi-thermalizedstate, a beam-plasma instability is triggered when the speed of the beam electrons exceeds the thermal velocity ofthe background electron population. This bump-on-tail instability then excites the exponential growth of dispersiveLangmuir waves ( L mode waves), starting from the initially low level of local thermal radiation field, which arepredominantly emitted in the same direction as the beam. As the Langmuir waves are convectively amplified alongtheir ray paths, the wave intensity becomes sufficiently high for nonlinear processes to start taking place. Accordingto the weak turbulence theory of plasmas, the most important nonlinear processes in this phase are the three-wavedecay/coalescence and nonlinear wave-particle scattering, which will first amplify the backscattered Langmuir modefrom thermal background, mostly due to the coupling of the forward-propagating L mode with the ion-acoustic( S ) mode. At later stages of the nonlinear turbulent processes, when the backward L mode reaches a sufficientlyhigh intensity level, subsequent wave-wave and wave-particle interactions promote its coupling with the forward L mode, thereby exciting the electromagnetic (transverse, T ) mode. Longer time-scale processes can finally lead to theexcitation of harmonics of the transverse mode. The entire series of events, which was broadly described above, isknown in the literature as the plasma emission process and is supported by several contributions published in theliterature (Tsytovich 1967; Kaplan & Tsytovich 1968; Melrose 1982; Goldman & DuBois 1982; Goldman 1983; Melrose1987; Cairns 1987a; Robinson & Cairns 1998a,b; Zlotnik et al. 1998; Kontar 2001a,b,c; Gosling et al. 2003; Li et al.2005a,b, 2006, 2008, 2011, 2012; Ratcliffe et al. 2012; Reid & Kontar 2017).The plasma emission mechanism has also been studied within the framework of the Generalized Weak Turbulence(GWT) theory (Yoon 2000, 2005, 2006; Yoon et al. 2012b). In this theory, the evolution of the turbulence levelin a thermal plasma is described by a set of kinetic equations for the velocity distribution functions of the plasmaspecies and for the spectral wave intensities of the various normal modes interacting with the particles. Examples ofapplications of the GWT theory to the study of the nonlinear evolution of the beam-plasma instability, taking intoaccount wave-wave and nonlinear wave-particle interactions, are given by Ziebell et al. (2001, 2008b); Gaelzer et al.(2008); Ziebell et al. (2008a, 2011a,b, 2012); Yoon et al. (2012a); Ziebell et al. (2014c,a,b, 2015, 2016). In particular,in Ziebell et al. (2014b) and Ziebell et al. (2015), all the steps involved in the plasma emission process are considered.In addition to the wave kinetic equations for the forward/backward L and S modes, containing all the usual linearand nonlinear interactions included in the traditional weak turbulence theory, the above-mentioned publications alsoincluded the kinetic equations for the forward/backward transverse modes. The full set of kinetic equations for particlesand waves is self-consistently solved during a sufficiently long evolution time that shows not only the excitation ofthe fundamental T mode, but the first two transverse harmonic modes as well. The results obtained by Ziebell et al.(2014b) and Ziebell et al. (2015) give additional support to the hypothesis that the underlying generating mechanismof the Type II/III solar radio bursts is the plasma emission process.In another theoretical front, some authors applied the technique of particle-in-cell (PIC) simulation to the questionof the plasma emission as the mechanism behind solar radio emissions. Some examples of contributions in this front are [email protected]@[email protected]@khu.ac.kr given by Kasaba et al. (2001), Karlick´y & Vandas (2007), Rhee et al. (2009a,b), Ganse et al. (2012a,b), and Thurgood& Tsiklauri (2015).Until now, the few attempts that have been made for a direct comparison, either qualitative or quantitative, betweenthe results obtained from the weak turbulence theory and from numerical simulations have all been of a limited scope.Rha et al. (2013) employed PIC simulation to test a theory for the generation of asymmetric superthermal tails inthe electronic distribution function. According to the theory (Yoon et al. 2012a), asymmetric tails can be created dueto the nonlinear interactions of electrons with L and S waves, when the ion-electron temperature ratio ( T i /T e ) varieswithin the range 0 . (cid:46) T i /T e (cid:46)
1. Ziebell et al. (2014c) employed simulation to test a theory for the generation ofa quasi-thermal electromagnetic radiation field, resulting as the time-asymptotic state of the nonlinear interaction ofthe transverse mode with the longitudinal modes and particles, in the absence of free energy sources ( e.g. , beams) andcollisions ( i.e. , no Brehmsstrahlung emission). On the other hand, Ratcliffe et al. (2014) performed a direct comparisonof the plasma emission mechanism with a PIC simulation. However, the weak turbulence formulation that was testedwas one-dimensional and the Langmuir wave kinetic equation contained only the three-wave decay term, besides theusual quasi-linear diffusion term. Moreover, the transverse mode kinetic equation was absent in their comparison.Here, we report for the first time a detailed comparison between the full two-dimensional (2D) plasma emissiontheory with a 2D particle-in-cell (PIC) simulation. The complete numerical solutions of 2D equations of generalizedweak turbulence (GWT) theory can be termed the weak turbulence (WT) simulation, which are directly comparedagainst numerical results obtained from the 2D PIC code simulation. Employing the same sets of physical parameters,the plasma emission process is simulated by both the 2D relativistic electromagnetic (EM) PIC simulation, employinga physical proton-to-electron mass ratio, and the GWT theory, employing the full set of 2D self-consistent kineticequations for the particles and for L , S and T waves. The results show a reasonably good quantitative agreementbetween the results from both approaches, with better agreement obtained for the low-beam speed regime.The organization of the paper is as follows: In section 2 a short overview of the generalized weak turbulence theoryis presented, with emphasis on the interpretation of the physical origin of different terms in wave-particle kineticequations. In section 3 a concise description of the PIC code is made. Then, in section 4 results from both approachesare presented and the comparison between them is carried out. Finally, in section 5 we summarize the major findingsand present our conclusions. WEAK TURBULENCE THEORYEssential theoretical developments of the generalized weak turbulence theory relevant to the study of plasma emissioncan be found in the papers by Yoon (2006); Yoon et al. (2012b); Ziebell et al. (2015). The theory describes weaklyturbulent nonlinear interactions of Langmuir ( L ), ion-sound ( S ), and transverse electromagnetic ( T ) waves, as well asthe electrons and protons. The longitudinal electric field wave energy density is given in spectral form by the sum ofLangmuir and ion-sound mode intensities, and is defined by (cid:68) δE (cid:107) (cid:69) k ,ω = (cid:88) σ = ± (cid:88) α = L,S I σα ( k ) δ ( ω − σω α k ) (1)where I σα ( k ) represents the individual mode intensity, σ = ± ω α k denotes the wave dispersionrelation. The electric and magnetic field wave energy densities for the transverse mode T is expressed in terms of the T mode wave intensity I σT ( k ) by (cid:10) δE ⊥ (cid:11) k ω = (cid:88) σ = ± I σT ( k ) δ ( ω − σω T k ) , (cid:10) δB (cid:11) k ω = (cid:88) σ = ± (cid:12)(cid:12)(cid:12)(cid:12) ckω T k (cid:12)(cid:12)(cid:12)(cid:12) I σT ( k ) δ ( ω − σω T k ) . (2)In the above σ = ± σ = +1 representing the forward propagation, while σ = − L , S , and T modes are given, respectively, by ω L k = ω pe (cid:18) k λ D (cid:19) , ω S k = kc S (1 + 3 T e /T i ) / (1 + k λ D ) / , ω T k = ( ω pe + c k ) / , (3)where ω pe = (cid:18) πn e e m e (cid:19) / , λ D = (cid:18) T e πn e e (cid:19) / = v th √ ω pe , v th = (cid:18) T e m e (cid:19) / c S = (cid:18) T e m i (cid:19) / . (4)In the above ω pe is the plasma frequency, n e , e , and m e being the electron number density, unit electric charge,and electron mass, respectively; λ D represents the Debye length, T e and T i being the electron and ion temperatures,respectively; and v th and c S stand for electron thermal speed and ion-sound speeds, respectively, m i being the protonmass.The equations of weak turbulence theory are made of kinetic equations that govern the dynamical evolution of theelectron velocity distribution function, F e ( v , t ), and spectral wave intensities, I ± L ( k , t ), I ± S ( k , t ), and I ± T ( k , t ), ∂F e ( v ) ∂t = πe m e (cid:88) σ = ± (cid:90) d k k k · ∂∂ v δ ( σω L k − k · v ) (cid:18) m e σω L k π k F e ( v ) + I σL ( k ) k k · ∂F e ( v ) ∂ v (cid:19) , (5) ∂I σL ( k ) ∂t = 4 πe m e k (cid:90) d v δ ( σω L k − k · v ) (cid:18) n e e F e ( v ) + πσω L k k · ∂F e ( v ) ∂ v I σL ( k ) (cid:19) + πe σω L k T e (cid:88) σ (cid:48) ,σ (cid:48)(cid:48) = ± (cid:90) d k (cid:48) µ k − k (cid:48) ( k · k (cid:48) ) k k (cid:48) | k − k (cid:48) | (cid:32) σω L k I σ (cid:48) L ( k (cid:48) ) I σ (cid:48)(cid:48) S ( k − k (cid:48) ) µ k − k (cid:48) − σ (cid:48) ω L k (cid:48) I σ (cid:48)(cid:48) S ( k − k (cid:48) ) µ k − k (cid:48) I σL ( k ) − σ (cid:48)(cid:48) ω L k − k (cid:48) I σ (cid:48) L ( k (cid:48) ) I σL ( k ) (cid:33) δ ( σω L k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω S k − k (cid:48) )+ σω L k e n e m e ω pe (cid:88) σ (cid:48) (cid:90) d k (cid:48) (cid:90) d v ( k · k (cid:48) ) k k (cid:48) (cid:26) n e e ω pe (cid:104) σω L k I σ (cid:48) L ( k (cid:48) ) − σ (cid:48) ω L k (cid:48) I σL ( k ) (cid:105) [ F e ( v ) + F i ( v )]+ πm e m i I σ (cid:48) L ( k (cid:48) ) I σL ( k )( k − k (cid:48) ) · ∂F i ( v ) ∂ v (cid:27) δ [ σω L k − σ (cid:48) ω L k (cid:48) − ( k − k (cid:48) ) · v ] , (6) ∂∂t I σS ( k ) µ k = 4 πµ k e m e k (cid:90) d v δ ( σω S k − k · v ) (cid:2) n e e [ F e ( v ) + F i ( v )]+ πσω L k (cid:18) k · ∂F e ( v ) ∂ v + m e m i k · ∂F i ( v ) ∂ v (cid:19) I σS ( k ) µ k (cid:21) + πe σω L k T e (cid:88) σ (cid:48) ,σ (cid:48)(cid:48) (cid:90) d k (cid:48) µ k [ k (cid:48) · ( k − k (cid:48) )] k k (cid:48) | k − k (cid:48) | (cid:16) σω L k I σ (cid:48) L ( k (cid:48) ) I σ (cid:48)(cid:48) L ( k − k (cid:48) ) − σ (cid:48) ω L k (cid:48) I σ (cid:48)(cid:48) L ( k − k (cid:48) ) I σS ( k ) µ k − σ (cid:48)(cid:48) ω L k − k (cid:48) I σ (cid:48) L ( k (cid:48) ) I σS ( k ) µ k (cid:19) δ ( σω S k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω L k − k (cid:48) ) , (7) ∂∂t I σT ( k )2 = πe σω T k m e ω pe (cid:88) σ (cid:48) ,σ (cid:48)(cid:48) (cid:90) d k (cid:48) ( k × k (cid:48) ) k k (cid:48) | k − k (cid:48) | (cid:32) k (cid:48) σ (cid:48) ω L k (cid:48) − | k − k (cid:48) | σ (cid:48)(cid:48) ω L k − k (cid:48) (cid:33) δ ( σω T k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω L k − k (cid:48) ) × (cid:18) σω T k I σ (cid:48) L ( k (cid:48) ) I σ (cid:48)(cid:48) LL ( k − k (cid:48) ) − σ (cid:48) ω L k (cid:48) I σ (cid:48)(cid:48) L ( k − k (cid:48) ) I σT ( k )2 − σ (cid:48)(cid:48) ω L k − k (cid:48) I σ (cid:48) L ( k (cid:48) ) I σT ( k )2 (cid:19) + πe σω T k T e (cid:88) σ (cid:48) ,σ (cid:48)(cid:48) (cid:90) d k (cid:48) µ k − k (cid:48) ( k × k (cid:48) ) k k (cid:48) | k − k (cid:48) | δ ( σω T k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω S k − k (cid:48) ) × (cid:32) σω T k I σ (cid:48) L ( k (cid:48) ) I σ (cid:48)(cid:48) S ( k − k (cid:48) ) µ k − k (cid:48) − σ (cid:48) ω L k (cid:48) I σ (cid:48)(cid:48) S ( k − k (cid:48) ) µ S k − k (cid:48) I σT ( k )2 − σ (cid:48)(cid:48) ω L k − k (cid:48) I σ (cid:48) L ( k (cid:48) ) I σT ( k )2 (cid:33) + πe σω T k m e (cid:88) σ (cid:48) ,σ (cid:48)(cid:48) (cid:90) d k (cid:48) | k − k (cid:48) | ( ω T k ) ( ω T k (cid:48) ) (cid:18) k · k (cid:48) ) k k (cid:48) (cid:19) δ ( σω T k − σ (cid:48) ω T k (cid:48) − σ (cid:48)(cid:48) ω L k − k (cid:48) ) × (cid:32) σω T k I σ (cid:48)(cid:48) L ( k − k (cid:48) ) I σ (cid:48) T ( k (cid:48) )2 − σ (cid:48) ω T k (cid:48) I σ (cid:48)(cid:48) L ( k − k (cid:48) ) I σT ( k )2 − σ (cid:48)(cid:48) ω L k − k (cid:48) I σ (cid:48) T ( k (cid:48) )2 I σT ( k )2 (cid:33) + σω T k e nm e ω pe (cid:88) σ (cid:48) (cid:90) d k (cid:48) (cid:90) d v ( k × k (cid:48) ) k k (cid:48) (cid:20) ˆ ne ω pe (cid:18) σω T k I σ (cid:48) L ( k (cid:48) ) − σ (cid:48) ω L k (cid:48) I σT ( k )2 (cid:19) [ F e ( v ) + F i ( v )]+ π m e m i I σ (cid:48) L ( k (cid:48) ) I σT ( k )2 ( k − k (cid:48) ) · ∂F i ( v ) ∂ v (cid:21) δ (cid:2) σω T k − σ (cid:48) ω L k (cid:48) − ( k − k (cid:48) ) · v (cid:3) , (8)where µ k = | k | λ De (cid:18) m e m i (cid:19) / (cid:18) T i T e (cid:19) / . (9)In (5)–(8), the dependence of each of the quantities, F e ( v , t ), I ± L ( k , t ), I ± S ( k , t ), and I ± T ( k , t ), on time t is implicit.Physical meanings of various terms in (5) – (8) have been expounded in the paper by Ziebell et al. (2015), which isbriefly recapitulated here. The electron particle kinetic equation (5) is given by the Fokker-Planck type of equationwhere the linear wave-particle resonance between the electrons and Langmuir turbulence is retained.The first term on the right-hand side of Langmuir wave kinetic equation (6), which is dictated by δ ( σω L k − k · v ),designates the spontaneous and induced emissions of L waves; the second term that contains the overall three-waveresonance condition δ ( σω L k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω S k − k (cid:48) ) describes the decay/coalescence processes involving two L modes andan S mode; the term dictated by nonlinear wave-particle resonance condition δ [ σω L k − σ (cid:48) ω L k (cid:48) − ( k − k (cid:48) ) · v ], depicts thespontaneous and induced scattering processes involving two Langmuir waves and the distribution of electrons as wellas stationary thermal protons.The first term on the right-hand side of S mode wave kinetic equation that contains the factor δ ( σω S k − k · v )corresponds to spontaneous and induced emissions of S waves; nonlinear terms with the factor δ ( σω S k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω S k − k (cid:48) )depicts the decay/coalescence involving an S mode and two L modes.For T mode we have ignored the linear wave-particle resonance term since it is impossible for the electrons tolinearly interact with the superluminal T mode. The k (cid:48) -integral term on the right-hand side of equation (8) thathas the three-wave resonance condition δ ( σω T k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω L k − k (cid:48) ) describes the coalescence of two L modes into a T mode, hence, this term is responsible for the harmonic emission at twice the plasma frequency, ω ∼ ω pe . The termassociated with the resonance factor δ ( σω T k − σ (cid:48) ω L k (cid:48) − σ (cid:48)(cid:48) ω S k − k (cid:48) ) represents the decay of L mode into an S mode anda T mode at the fundamental plasma frequency, ω ∼ ω pe , which is one of the processes that leads to the fundamentalemission. The k (cid:48) -integral with the factor δ ( σω T k − σ (cid:48) ω T k (cid:48) − σ (cid:48)(cid:48) ω L k − k (cid:48) ) describes the merging of a T mode and an L mode into the next higher harmonic T mode. This is known as the incoherent Raman scattering, and is responsiblefor higher-harmonic plasma emission, including the third harmonic emission, ω ∼ ω pe . The term with the condition δ [ σω T k − σ (cid:48) ω L k (cid:48) − ( k − k (cid:48) ) · v ] attached represents the spontaneous and induced scattering processes involving T and L modes as well as the charged particles. This process can be termed the transformation of L mode into the radiation,mediated by the particles, and is another mechanism responsible for the fundamental emission.Ziebell et al. (2015) numerically solve the set of equations (5) – (8), and by carefully turning various terms on oroff, they carried out the detailed diagnostics of each process. The purpose of the present analysis is not to repeat thedetailed tasks already carried out by Ziebell et al. (2015), but rather, the major aim of the present work is to validatethe efficacy of the WT simulation method by testing the numerical solutions against the PIC code simulation, whichis more rigorous. PARTICLE-IN-CELL SIMULATIONWe have carried out a series of two-dimensional relativistic and electromagnetic particle-in-cell (PIC) simulations.The simulation box size is L X = L Y = 64 . c/ω pe = 1024 v th /ω pe , where c is the speed of light, v th is the electronthermal speed and ω pe is the electron plasma frequency. The number of grids is N X = N Y = 1024 so that the gridsize is ∆ x = √ λ D . The simulation time step is ∆ t = ∆ x/ c , which satisfies the Courant-Friedrichs-Lewy (CFL)condition, where ∆ t is the time step.The number of particles is 200 per grid per species, and we used three kinds of species: protons, background electrons,and beam electrons. The proton-to-electron mass ratio is realistic, m p /m e = 1836, where m p is the proton mass and m e is the electron mass. The electron thermal speed is v th /c = 4 . × − . The electron temperature is 7 times higherthan proton temperature, T i /T e = 1 /
7, where T i is the ion (proton) temperature and T e is the electron temperature.The beam consists of 0 .
1% of the total electron content, namely, n b /n = 10 − , where n b is the beam numberdensity and n is the electron number density. The Maxwellian beam temperature is the same as the backgroundelectron temperature, T b = T e , where T b is the electron beam temperature. The plasma to cyclotron frequency ratiois ω pe / Ω ce = 100, where Ω ce is electron cyclotron frequency. The boundary conditions are periodic for both X - and Y -directions.The above described parameters are chosen as in the paper by Ziebell et al. (2015). We have carried out threesimulations where we have varied the average beam drift speed, V b /v th , from 6, to 8, to 10. Again, these choices arethe same as those considered by Ziebell et al. (2015). NUMERICAL ANALYSISWe solved the equations of GWT turbulence theory, that is, we have carried out the WT simulation, and also carriedout the PIC code simulation under the same set of initial conditions so that direct comparisons can be made, which isunprecedented. In the WT simulation we have taken the protons as a stationary background with the two-dimensionalvelocity distribution function given by F i ( v ) = m i πT i exp (cid:18) − m i v T i (cid:19) , (10)where T i is the proton thermal speed.In the PIC code, on the other hand, the protons are free to dynamically evolve, as there are no reasons to fix theprotons as a neutralizing background. The initial electron velocity distribution function is the same as that adopted inthe analysis by Ziebell et al. (2015), namely, the electrons are composed of a backward drifting background Maxwelliancomponent plus a tenuous forward drifting Maxwellian distribution of electron beam, which in two dimensions, is givenby F e ( v ,
0) = (cid:18) − n b n (cid:19) m i πT e exp (cid:18) − m e v ⊥ T e − m e ( v (cid:107) + V ) T e (cid:19) + n b n m e πT b exp (cid:18) − m e v ⊥ T b − m e ( v (cid:107) − V b ) T b (cid:19) . (11)We may define the thermal speeds associated with each component in their respective drifting frame, v th =(2 T e /m e ) / and v tb = (2 T b /m e ) / . Here, and V and V b are the drift velocities of the background and the for-ward beam, respectively, where the background drift velocity V is assumed in order to preserve the zero currentcondition in the proton frame, V = n b V b n − n b . (12)The PIC code simulation also initialized the electron configuration in the same manner as in the analytical initialvelocity distribution function. Initial spectral forms for L , S , and T mode intensities in the WT simulation are chosenas in Ziebell et al. (2015), I σL ( k ,
0) = T e π
11 + 3 k λ De ,I σS ( k ,
0) = T e π k λ De (cid:18) k λ De k λ De (cid:19) / (cid:82) d v δ ( σω S k − k · v )( F e + F i ) (cid:82) d v δ ( σω S k − k · v )[ F e + ( T e /T i ) F i ] ,I σT ( k ,
0) = T e π
11 + c k /ω pe . (13)Of course, in the PIC code simulation, the initial noise is present in the system so that one does not need to specifythe initial spectral form for the modes.The initial input parameters are taken to be the same as that considered by Ziebell et al. (2015). These are1 n e λ D = 5 × − , T e = T b , T i T e = 17 , v th c = 2 T e m e c = 4 . × − , n b n = 10 − . (14)Among the above input parameters, the plasma parameter, 1 / ( nλ D ) is strictly applicable only for the WT simulation.As we already explained in the previous section, the PIC code simulation was designed with the same specification as(14), except for the plasma parameter. In the PIC code, the number of particles per cell is loosely connected to theplasma parameter, as recently demonstrated by L´opez & Yoon (2018), but there is no precise way to strictly associatethe value of plasma parameter to the PIC code design parameters such as the dimensions, the number of particles percell, etc.We adopt the normalized wave vector, velocity, and time, k v th ω pe , v v th , ω pe t, (15)in plotting the numerical results. Of course, the proton-to-electron mass ratio is realistic, m i /m e = 1836. Ziebell et al.(2015) considered three cases of normalized forward beam speed, V b /v th = 6, 8 and 10. The PIC code simulations arecarried out with the same conditions so that direct comparisons can be made.In the following, we present the quantitative comparisons of the results obtained by solving the equations of WTtheory, that is, WT simulation, and PIC code simulation. In the WT analysis, where we revisited the earlier approach,we have considered a wider velocity space compared to that of Ziebell et al. (2015), with the objective of minimizing theeffect of boundary conditions. We present the results in the same format as in the paper by Ziebell et al. (2015), exceptthat we choose to present only the snapshots of electron velocity distribution function, Langmuir wave spectrum, andtransverse wave spectrum, or equivalently, the radiation pattern, as these quantities are dynamically important. Ionacoustic mode is excited in both PIC code simulation and WT theory, but as their wave intensities are generally low,we choose not to plot the results. In the PIC code simulation, in particular, fluctuations in the ion acoustic wavefrequency range is too weak to distinguish from numerical noise. It should also be noted that the theoretical ionacoustic mode turbulence intensity discussed in the paper by Ziebell et al. (2015) is also very low. We thus focus onthose quantities that may be meaningfully compared, namely, electron velocity distribution function, high-frequencylongitudinal electric field intensity in the Langmuir wave frequency range, and perturbed magnetic field spectrumconstructed from PIC simulation, and the corresponding theoretical electron velocity distribution, Langmuir waveintensity, and transverse radiation intensity computed on the basis of WT theory.4.1. Case 1: V b /v th = 6In Figure 1 we show the electron velocity distribution function (VDF) computed on the basis of weak turbulence(WT) simulation (top four panels) and the electron velocity phase space distribution constructed from the PIC codesimulation (bottom four panels), versus two dimensional normalized perpendicular and parallel speeds, v ⊥ /v th and v (cid:107) /v th , for the first case of V b /v th = 6 (case 1). Snapshots of the electron VDF at four different times correspondingto ω pe t = 500, 1000, 1500, and 2000, are shown.On the basis of the direct comparisons between the theory, or WT simulation, and PIC simulation, it can be seenthat there exists an overall quantitative, and even qualitative agreement between the two methods. Specifically, for ω pe t = 500, both methods produce similar results in that the beam has undergone a partial plateau formation, butsignificant positive gradient, ∂F e /∂v (cid:107) > ω pe t = 2000, the plateau formation is almost complete in the case of WT simulation. Bycontrast, for the PIC code case, a small but finite positive gradient still persists at ω pe t = 2000. Despite such a ratherinsignificant difference, the two methods agree quite well.The Langmuir turbulence spectrum is shown in Figure 2. In the top four panels, the Langmuir spectral intensity I L ( k ) is plotted versus k ⊥ v th /ω pe and k (cid:107) v th /ω pe . In computing for the forward and backward modes, we restrictedourselves to only positive k (cid:107) > σ = ± I +1 L ( k ) in k (cid:107) > I − L ( k ) in the negative k (cid:107) space. Thus, the portion of k space corresponding to k (cid:107) > I +1 L ( k ), while the other half space with k (cid:107) < L mode, I − L ( k ). For the bottom four panels, which correspond to the simulated electric field spectral intensity, itis not so easy to single out only the contribution from the eigenmode, which L mode is. Instead, we have filteredout the low- and high-frequency fluctuations, and focused only on the spectrum that roughly encompasses the plasmafrequency. Note, however, that the spectrum in the vicinity of plasma frequency not only contains the Langmuir modeintensity, which is the linear eigenmode of the plasma, but also nonlinear eigenmode – see, Rhee et al. (2009b). Thesimulated electric field spectrum does not distinguish the two. The total electric field may also contain the radiationmode, but in order to eliminate the transverse component as much as possible, we have selected only the electricfield component parallel to the beam propagation direction before implementing the fast Fourier transformation of thesimulated electric field. These subtle differences notwithstanding, the comparison between the theoretical intensityand simulated spectrum shows a rather striking resemblance to each other.Specifically, for early time, ω pe t = 500, the enhanced forward-propagating component (the primary L ) can be seento be excited in both WT and PIC simulations, which is the result of initial gentle bump-on-tail instability. Atthis relatively early time the weak backward propagating L mode is also evident. For ω pe t = 1000 and beyond, thebackscattering of primary L mode into the back-scattered L mode, via a combined three-wave decay process andnonlinear scattering off ions becomes increasingly more visible (Ziebell et al. 2001, 2008a, 2012, 2014b, 2015). In anoverall sense, the characteristic time scale of wave evolution and the spectral feature at each stage of time evolution,in particular, the formation of semi-ark shape spectra in both WT and PIC code simulations, are quite consistent.One minor difference is that, whereas in the WT calculation the growth of near k ∼ k (cid:39) Figure 1 . Case 1 ( V b /v th = 6): Time evolution of the electron velocity distribution function (VDF) F e ( v ) versus v ⊥ /v th and v (cid:107) /v th , for four different time steps corresponding to ω pe t = 500, 1000, 1500, and 2000. Top four panels correspond to WTsimulation, while the bottom four panels show results from PIC code simulation. Figure 2 . Case 1 ( V b /v th = 6): Time evolution of the Langmuir wave spectral intensity I L ( k ), in the case of WT simulation– top four panels, and total electric field intensity δE ( k ), in the case of PIC code simulation – bottom four panels, versus k x v th /ω pe and k z v th /ω pe , for ω pe t = 500, 1000, 1500, and 2000. Figure 3 . Case 1 ( V b /v th = 6): Time evolution of the radiation spectral intensity I T ( k ), in the case of WT simulation – top fourpanels, and total magnetic field intensity δB ( k ), in the case of PIC code simulation – bottom four panels, versus k x v th /ω pe and k z v th /ω pe , for ω pe t = 500, 1000, 1500, and 2000. T mode) spectrum computed on the basis of WT theoretical method (topfour panels), and the magnetic field spectrum constructed from PIC simulation (bottom four panels). According to WTcalculation, for early time ω pe t = 500, finite level of fundamental emission is first generated. Then at ω pe t = 1000 andbeyond, second harmonic mode gradually appears until at ω pe t = 2000, the fundamental/harmonic bi-model radiationpattern is established. In contrast, in the present case of weak beam speed, V b /v th = 6, the magnetic field spectrum isdominated by the noise almost throughout the entire simulated time domain. Only at the very end, ω pe t = 2000, doesa weak signature of fundamental/harmonic pair emission structure becomes barely discernible. This seems to indicatethe need for a higher number of particles per cell in order to reduce the numerical noise.4.2. Case 2: V b /v th = 8Ziebell et al. (2015) also considered the second case corresponding to V b /v th = 8 (case 2), which we consider next.Figure 4 displays the electron VDF in the same format as Figure 1. As with case 1, the time evolution of electronVDF is qualitatively similar for both WT and PIC simulations, especially for relatively early time, ω pe t = 500. Note,however, that as the system evolves in time, while the qualitative agreement is still maintained, upon close inspection,some quantitative discrepancies become noticeable. Most appreciable is the fact that while the two dimensionalstructure associated with the theoretically computed electron VDF is defined by elliptical outer contours for thebeam population, the PIC simulated electron VDF exhibits a distinct broadening of the beam along v ⊥ , while thehighest parallel velocity portion of the beam population does not suffer from such broadening. Another feature that isnoteworthy is that for ω pe t = 2000, the theoretical VDF still shows a mild positive parallel derivative associated withthe beam, whereas for the simulated VDF such a feature is almost completely gone. In spite of these, it is quite fairto say that the WT simulation features a good overall comparison against PIC simulation.The Langmuir turbulence spectrum as well as the simulated longitudinal electric field spectrum are plotted inFigure 5. Again, despite the subtle differences in the two quantities plotted, as already explained (namely, thetheoretical quantity is the Langmuir wave intensity while the simulated quantity is longitudinal electric field fluctuationcentered around the plasma oscillation frequency, which may contain both the Langmuir mode as well as the nonlineareigenmode), the overall agreement is rather remarkable. Upon direct comparison, it is seen that the contours for boththe theoretical spectrum and the simulated intensity evolve into more or less ring-like morphology in two dimensional k space. For the present case of higher beam speed, more free energy is available for the excitation of Langmuirinstability. Consequently, the entire instability and ensuing nonlinear processes develop much more rapidly. Thisresults in many back-and-forth decay and scattering processes, which leads to the multiply-peaked structures in thewave spectrum along k (cid:107) . This is particularly apparent in the WT simulation. In the PIC code simulation, on theother hand, owing to the inherent noise, the clear delineation of multiple-peak structure along k (cid:107) is not so evident.However, upon close examination, especially for ω pe t = 2000, the forward propagating longitudinal mode does indeedshow a faint evidence for the structure along parallel wave number. As with the first case of V b /v th = 6, the PIC codesimulation exhibits a rather robust Langmuir condensation phenomenon, while the WT simulation shows very weakor no evidence for Langmuir condensation.Moving on the transverse EM radiation ( T mode) for WT simulation, and the transverse magnetic field fluctuationspectrum for PIC simulation, Figure 6 plots these spectra. In the present case of V b /v th = 8, the fundamental plusweak harmonic emission already takes place at ω pe t = 500, in the case of WT calculation. The pair emission patterngradually and monotonically increases in intensity until the end of the computation, namely, ω pe t = 2000. In contrast,as with the previous case, PIC simulation shows no identifiable radiation emission pattern for relatively early times, ω pe t = 500 and 1000. Even at ω pe t = 1500, the pair emission pattern becomes barely visible with great difficulty.However, at the end of the simulation period, ω pe t = 2000, the radiation emission pattern now becomes quite discernibleover the background noise. Again, the present PIC code simulation study implies the difficulty in faithfully simulatingthe plasma emission radiation, and calls for a higher number of particles per cell and thus, quieter simulation, whichis, needless to say more computationally demanding, and is beyond the scope of the present work. In this regard, theWT simulation is advantageous, since such an approach is free from the noise issue.4.3. Case 3: V b /v th = 10The third case study is for V b /v th = 10 (case 3), which was considered by Ziebell et al. (2015). In their paper, theauthors speculated that the applicability of weak turbulence theory, and more specifically, the use of weak-growthrate formula inherent in the standard WT formalism might be suspect, but they did not have any standard to verifysuch a suspicion. For the present relatively high-beam speed, their numerical solution for the electron VDF featureda particularly undesirable aspect of the velocity plateau spreading widely until it reached the boundary, where the2 Figure 4 . Case 2 ( V b /v th = 8): Electron velocity distribution function (VDF) F e ( v ) versus v ⊥ /v th and v (cid:107) /v th , for four differenttime steps corresponding to ω pe t = 500, 1000, 1500, and 2000. Top four panels correspond to WT simulation, while the bottomfour panels show results from PIC code simulation. Figure 5 . Case 2 ( V b /v th = 8): Langmuir wave spectral intensity I L ( k ), in the case of WT simulation – top four panels, andtotal electric field intensity δE ( k ), in the case of PIC code simulation – bottom four panels, versus k x v th /ω pe and k z v th /ω pe ,for ω pe t = 500, 1000, 1500, and 2000. Figure 6 . Case 2 ( V b /v th = 8): Radiation spectral intensity I T ( k ), in the case of WT simulation – top four panels, and totalmagnetic field intensity δB ( k ), in the case of PIC code simulation – bottom four panels, versus k x v th /ω pe and k z v th /ω pe , for ω pe t = 500, 1000, 1500, and 2000. Figure 7 . Case 3 ( V b /v th = 10): Electron velocity distribution function (VDF) F e ( v ) versus v ⊥ /v th and v (cid:107) /v th , for four differenttime steps corresponding to ω pe t = 500, 1000, 1500, and 2000. Top four panels correspond to WT simulation, while the bottomfour panels show results from PIC code simulation. ω pe t = 2000, which appears to show that the beam has spread to the upper boundary of the figure, is actually nota problem since the actual velocity boundary is much wider than what is shown in the figure. In Figure 7 bottompanels, we display the results of PIC code simulation. As the comparison readily shows, the WT simulation enjoys atleast a quantitative agreement, but only in the relative early time periods corresponding to ω pe t = 500 and 1000. For ω pe t = 1500 and 2000, it is evident that the WT calculation exaggerates the velocity space diffusion of the beam. ForPIC code simulation, the beam is spread along v ⊥ as in case 2, but the peak velocity portion of the beam along v (cid:107) doesnot evolve much. It is interesting to note that for ω pe t = 2000, the PIC simulation indicates parallel acceleration ofelectrons in both positive and negative portions of v (cid:107) axis. This feature is absent in the WT calculation. In an overallsense, while there exist some discrepancies, the qualitative agreement is arguably present, especially for relatively earlytimes. This assessment notwithstanding, the present case 3 study implies that the beam speed of V b /v th = 10 may beat the limit of validity for WT theory, at least from the standpoint of electron VDF.However, as for the Langmuir turbulence, the agreement between WT method and PIC code simulation is not thatbad. Indeed, as Figure 8 demonstrates, the Langmuir turbulence spectrum and the simulated longitudinal electricfield spectrum are qualitatively similar. The overall morphologies of the two dimensional spectra for both methodssomehow produce rather consistent results, including the fact that even the WT simulation generates intense Langmuircondensation, which was missing in the previous two cases.Finally, moving on to the radiation emission pattern, jumping to the final state corresponding to ω pe t = 2000, shownin Figure 9 top and bottom panels, one may immediately appreciate the similarities in the WT versus PIC simulatedradiation pattern, which show fundamental/second-harmonic emissions, as well as weak third harmonic emission. Evenin the PIC code result, numerical noise notwithstanding, the third harmonic emission is easy to identify. Now, as forrelatively early time periods, especially for ω pe t = 500 and 1000, the PIC code simulation is still too noisy in order tovisually identify and discern clear radiation emission at the harmonics. This contrasts to the WT calculation, whichis free of noise problem. By the time the PIC code simulation is carried out to ω pe t = 1500, however, the radiationemission pattern begins to manifest itself, albeit, rather faintly.In the paper by Ziebell et al. (2015), the authors have analyzed the detailed physics of the plasma emission. Specif-ically, they discussed that the fundamental emission takes place as a result of combined processes of L mode decayinginto T and S modes, as well as the scattering involving the beating of L and T modes mediated by the particles. Theymention that both the decay and scattering mechanisms are governed by coupling coefficient of the form,( k × k (cid:48) ) k k (cid:48) ∝ sin ϑ, (16)where ϑ represents the angle between the two vectors k (cid:48) and k . They also argued that the fundamental emissionalong the direction specified by ϑ = 0 should be prohibited, thus resulting in the dipole radiation. In the numericalsimulation, whether it be based upon the WT theory, or it is by means of direct PIC code method, the dipolar patternassociated with the fundamental radiation is difficult to discern, since it involves a narrow region around k ∼ k × k (cid:48) ) k k (cid:48) (cid:16) k (cid:48) − | k − k (cid:48) | (cid:17) ∼ sin ϑ cos ϑ. (17)This implies a quadrupole pattern, but since the radiation emission generally involves multiple wave modes, the abovecoupling coefficient is to be integrated over k (cid:48) , or equivalently, ϑ , hence, the strict quadrupole emission is not evidentin reality.Ziebell et al. (2015) also confirmed the earlier theories of third- and higher harmonic emission (Zheleznyakov &Zlotnik 1974; Cairns 1987b; Kliem et al. 1992). The coupling coefficient for the higher-harmonic emission is given by (cid:18) k · k (cid:48) ) k k (cid:48) (cid:19) ∼ ϑ. (18)Ziebell et al. (2015) also reminded the readers that the third-harmonic plasma emission associated with the solar radiobursts is quite rare (Takakura & Yousef 1974; Zlotnik et al. 1998; Brazhenko et al. 2012).Finally, Ziebell et al. (2015) analyzed the details of the various emission mechanisms by artificially turning certainterms in the T wave kinetic equation (8) on or off in order to investigate the consequences thereof. By employing such7 Figure 8 . Case 3 ( V b /v th = ‘0): Langmuir wave spectral intensity I L ( k ), in the case of WT simulation – top four panels, andtotal electric field intensity δE ( k ), in the case of PIC code simulation – bottom four panels, versus k x v th /ω pe and k z v th /ω pe ,for ω pe t = 500, 1000, 1500, and 2000. Figure 9 . Case 3 ( V b /v th = 10): Radiation spectral intensity I T ( k ), in the case of WT simulation – top four panels, and totalmagnetic field intensity δB ( k ), in the case of PIC code simulation – bottom four panels, versus k x v th /ω pe and k z v th /ω pe , for ω pe t = 500, 1000, 1500, and 2000. SUMMARY AND CONCLUSIONSThe purpose of the present paper has been to investigate the plasma emission process, by making use of two differentapproaches, and to discuss the compatibility between the results obtained from these two approaches. The plasmaemission is generally acknowledged to be the fundamental radiation emission mechanism for solar type II and type IIIradio bursts phenomena.In one of the approaches, we have utilized a self-consistent system of coupled equations, obtained using the frameworkof weak turbulence (WT) theory, in order to study the time evolution of the velocity distribution function of theelectrons and of the spectra of electrostatic and transverse waves. The formulation incorporates the effects of differentphysical mechanisms. In such an approach the roles of different physical processes can be identified unambiguously byturning certain terms on and off, which has in fact been done by Ziebell et al. (2015). The WT formulation, thereforeis a convenient tool to test the validity of various theories proposed in the literature for the generation of plasmaemission. The complete set of equations can also be solved without a priori assumptions in order to quantitativelyanalyze the plasma emission process, as was done by Ziebell et al. (2015). On the other hand, the WT theory, as withany analytical theory, is based upon a series of assumptions, whose limits of applicability have not been clearly definedor tested.The second approach employed in the present investigation has been the direct numerical simulation, based uponthe particle-in-cell (PIC) paradigm. The PIC simulations rely on a smaller number of theoretical constraints than theapproach based on WT theory, as such an approach basically solves the Lorentz equation of motions for a collectionof charged particles, plus the Maxwell’s equation. Nonetheless, such a numerical approach is not without someshortcomings. For instance, the necessity of discretization and the finite grid size in velocity space and in wave numberspace, places some limitations on resolution of very small wavelengths. The simulation method also needs to take intoaccount a large number of particles inside a cell in order to reduce the numerical noise, with the consequent burdenon the requirements for computational resources. In addition, it is not so easy to make clear diagnostics about thephenomena which occur in the system, since all mechanisms act simultaneously, and cannot be arbitrarily turned onor off for verification of certain physical processes that operate in the system. In contrast, the WT approach naturallylends itself to such manipulations. The WT method requires far less computational resources.We have made use of these two approaches to study a plasma system containing one ion population and two electronpopulations, constituted by an initially Maxwellian background and a tenuous beam. We have considered the sameparameters, which have been utilized in a previous paper in which the WT equations have been employed in orderto discuss in detail the physical mechanisms involved in the plasma emission process (Ziebell et al. 2015). The WTanalysis of the present paper is essentially the same as that of Ziebell et al. (2015), except that we have considered awider velocity space in order to minimize the boundary effects. The the initial setup for the PIC code simulation isconsistent with that of WT analysis by Ziebell et al. (2015). Of course, the PIC code assumes additional parameters,such as the grid size and the number of particles per cell, etc., but the physical condition is consistent with that ofZiebell et al. (2015).The numerical results discussed in detail indicate that the results obtained with the WT approach and with the PICsimulations are largely compatible. Regarding the evolution of the electron distribution function, both approachesshow the formation of plateau in the beam region, within compatible time scales. The agreement between the WTand the PIC results is more noticeable for the case of lower beam velocity, which is not too surprising, since the WTtheory is based upon the assumption of weak wave growth and low wave energy density when compared to the particlethermal energy density. For the intermediate and high beam velocity cases considered in the present paper, we noticedthat after the formation of plateau , the high velocity part of the beam in the WT results is slightly wider along v ⊥ direction than in the PIC results. Moreover, the high velocity case (specifically for V b /v th = 10) WT calculationresulted in the formation of an extended tail along the forward direction, which is not seen in the PIC simulationresults. This appears to be an indication that the high beam velocity case of V b /v th = 10, when all other parametersare held constant, corresponds to the limit of applicability of the WT equation.0Regarding the spectra of waves, the WT approach singles out spectra for each eigenmode, which is built intothe theory. In contrast, for PIC simulation, the various eigenmodes, such as Langmuir or transverse waves mustbe carefully interpreted. For instance, in order to contruct the spectrum of Langmuir waves on the basis of PICsimulation, we ave considered only the electric field component along the direction parallel to the beam velocity. Thismay separate electrostatic waves from transverse waves, but may still retain some amount of intensity associatedto nonlinear harmonics of electrostatic waves, hence the interpretation of PIC simulation result, as far as Langmuirwaves are concerned, is not without some ambiguities. Nevertheless, the results obtained with the two approaches arelargely compatible. Both the WT results and the PIC results show the formation of the primary Langmuir wave peak,with comparable widths along the directions of parallel and perpendicular wave numbers, the growth of backwardpropagating waves, which are also characterized by consistent widths in wave number space, and also the spread ofthe primary peak towards the region of smaller wave numbers. A discrepancy remains, mostly in the cases of lowand intermediate beam velocities, namely, V b /v th = 6 and 8, in that whereas the PIC simulation shows the earlyappearance of waves for k (cid:39)
0, the so-called Langmuir condensation effect, in the WT results the region of very smallvalues of k is not quite attained, until the final computational time attained in our analysis. In the case of higherbeam velocity which has been considered, namely, V b /v th = 10, on the other hand, the WT results also show somegrowth of waves at k (cid:39)
0. The reason for this localized discrepancy between the two approaches is not yet completelyunderstood. In an overall sense, however, the agreement between WT and PIC methods are more consistent than theelectron velocity distribution, which is interesting.In order to obtain information about the spectra of electromagnetic waves in the case of PIC simulation results, wehave taken into account the total magnetic field intensity. The spectrum of the magnetic field fluctuations is thenused for comparison with the spectrum of transverse waves computed on the basis of WT method. The comparativeanalysis produced largely favorable results, although because of numerical noise associated with the PIC code, andbecause of the generally low level of radiation, the relatively early time results do not clearly show easily identifiableplasma emission. The compatibility between WT and PIC results was seen to improve with the increase of beamvelocity. Clear demonstration of the radiation emission was difficult to show with the limited numerical setup adoptedin the PIC code. This implies that the simulation of plasma emission requires large number of particles per cell inorder to reduce the numerical noise, which requires high computational resources.To conclude, the PIC code simulation is supposed to be more rigorous, but it necessitates computationally intenseefforts. In contrast, while the WT theory is a reduced approach, the comparative analysis presented herewith providedevidence that suggests that the use of WT theory can be reliable, if it is carefully applied to a parameter regime for whichthe theory is valid. The present investigation has focused on three examples, but more systematical statistical surveyof the parameter space could be carried out in order to further establish the region of validity of WT approach. Thepresent paper also indicates the possibility of improved agreement between the WT approach and the PIC simulationapproach if the numerical noise can be reduced in the PIC simulations. In short, we find that both the WT theoryand PIC code simulations are useful research tools in the fundamental study of solar radio bursts problem, as they aremutually complementary.The present research was supported in part by the BK21 plus grant from the National Research Foundation (NRF) ofthe Republic of Korea, to Kyung Hee University. L.F.Z. and R.G. acknowledge support provided by Conselho Nacionalde Desenvolvimento Cient´ıfico e Tecnol´ogico (CNPq), grants No. 304363/2014-6 and 307626/2015-6. This study wasalso financed in part by the Coordena¸c˜ao de Aperfei¸coamento de Pessoal de N´ıvel Superior - Brasil (CAPES) - FinanceCode 001. P.H.Y. acknowledges the grant from the GFT Charity, Inc., to the University of Maryland. This work waspartly carried out while P.H.Y. was visiting Korea Astronomy and Space Science Institute (KASI). E.L. acknowledgessupport by Space Core Technology Development Program through the NRF of Korea funded by Ministry of Scienceand ICT (NRF-2017M1A3A3A02016781). REFERENCES
Bastian, T. S., Benz, A. O., & Gary, D. E. 1998, ARA&A, 36,131–188, doi:
Benz, A. O. 2017, LRSP, 14, doi:
Brazhenko, A. I., Melnik, V. N., Konovalenko, A. A., et al. 2012,OAP, 25, 181–183 Cairns, I. H. 1987a, JPlPh, 38, 169,doi: —. 1987b, JPlPh, 38, 199, doi:
Cane, H. V., Erickson, W. C., & Prestage, N. P. 2002, JGR, 107,1315, doi: Ergun, R. E., Larson, D., Lin, R. P., et al. 1998, ApJ, 503,435–445, doi:
Gaelzer, R., Ziebell, L. F., Vi˜nas, A. F., Yoon, P. H., & Ryu,C. M. 2008, ApJ, 677, 676–682, doi:
Ganse, U., Kilian, P., Spanier, F., & Vainio, R. 2012a, ApJ, 751,145, doi:
Ganse, U., Kilian, P., Vainio, R., & Spanier, F. 2012b, SoPh,280, 551–560, doi:
Ginzburg, V. L., & Zheleznyakov, V. V. 1958, SvPhA, 2, 653Goldman, M. V. 1983, SoPh, 89, 403–442,doi:
Goldman, M. V., & DuBois, D. F. 1982, PhFl, 25, 1062,doi:
Gosling, J. T., Skoug, R. M., & McComas, D. J. 2003, GeoRL,30, 1697, doi:
Hudson, H. 2011, SSRv, 158, 5–41,doi:
Kaplan, S. A., & Tsytovich, V. N. 1968, SvPhA, 11, 956Karlick´y, M., & Vandas, M. 2007, P&SS, 55, 2336–2339,doi:
Kasaba, Y., Matsumoto, H., & Omura, Y. 2001, JGR, 106,18693, doi:
Kliem, B., Krueger, A., & Treumann, R. A. 1992, SoPh, 140,149–160, doi:
Kontar, E. P. 2001a, PPCF, 43, 589, doi: —. 2001b, SoPh, 202, 131–149, doi: —. 2001c, A&A, 375, 629–637,doi:
Li, B., Cairns, I. H., & Robinson, P. A. 2008, JGR, 113, A06104,doi: —. 2011, ApJ, 730, 20, doi: —. 2012, SoPh, 279, 173–196, doi:
Li, B., Robinson, P. A., & Cairns, I. H. 2006, PhPl, 13, 092902,doi:
Li, B., Willes, A. J., Robinson, P. A., & Cairns, I. H. 2005a,PhPl, 12, 012103, doi: —. 2005b, PhPl, 12, 052324, doi:
L´opez, R. A., & Yoon, P. H. 2018, JGR,doi:
Melrose, D. B. 1982, AuJPh, 35, 67–86, doi: —. 1987, SoPh, 111, 89–101, doi:
Ratcliffe, H., Bian, N. H., & Kontar, E. P. 2012, ApJ, 761, 176,doi:
Ratcliffe, H., Brady, C. S., Che Rozenan, M. B., & Nakariakov,V. M. 2014, PhPl, 21, 122104, doi:
Reid, H. A. S., & Kontar, E. P. 2017, A&A, 598, A44,doi:
Reid, H. A. S., & Ratcliffe, H. 2014, RAA, 14, 773–804,doi:
Reiner, M. J., Goetz, K., Fainberg, J., et al. 2009, SoPh, 259,255–276, doi:
Rha, K., Ryu, C.-M., & Yoon, P. H. 2013, ApJ, 775, L21,doi:
Rhee, T., Ryu, C.-M., Woo, M., et al. 2009a, ApJ, 694, 618–625,doi:
Rhee, T., Woo, M., & Ryu, C.-M. 2009b, JKPS, 54, 313–316,doi:
Robinson, P. A., & Cairns, I. H. 1998a, SoPh, 181, 363–394,doi: —. 1998b, SoPh, 181, 395–428, doi:
Takakura, T., & Yousef, S. 1974, SoPh, 36, 451–458,doi:
Thurgood, J. O., & Tsiklauri, D. 2015, A&A, 584, A83,doi:
Tsytovich, V. N. 1967, SvPhU, 9, 805,doi:
Wild, J. P. 1950a, AuSRA, 3, 399–408, doi: —. 1950b, AuSRA, 3, 541–557, doi:
Wild, J. P., & McCready, L. L. 1950, AuSRA, 3, 387–398,doi:
Wild, J. P., Murray, J. D., & Rowe, W. C. 1954, AuJPh, 7,439–459, doi:
Yoon, P. H. 2000, PhPl, 7, 4858, doi: —. 2005, PhPl, 12, 042306, doi: —. 2006, PhPl, 13, 022302, doi:
Yoon, P. H., Hong, J., Kim, S., et al. 2012a, ApJ, 755, 112,doi:
Yoon, P. H., Ziebell, L. F., Gaelzer, R., & Pavan, J. 2012b, PhPl,19, 102303, doi:
Zheleznyakov, V. V., & Zlotnik, E. Y. 1974, SoPh, 36, 443–449,doi:
Ziebell, L. F., Gaelzer, R., Pavan, J., & Yoon, P. H. 2008a,PPCF, 50, 085011 (15pp), doi:
Ziebell, L. F., Gaelzer, R., & Yoon, P. H. 2001, PhPl, 8, 3982,doi: —. 2008b, PhPl, 15, 032303, doi:
Ziebell, L. F., Petruzzellis, L. T., Yoon, P. H., Gaelzer, R., &Pavan, J. 2016, ApJ, 818, 61,doi:
Ziebell, L. F., Yoon, P. H., Gaelzer, R., & Pavan, J. 2012, PPCF,54, 055012, doi: —. 2014a, PhPl, 21, 012306, doi: —. 2014b, ApJL, 795, L32, doi:
Ziebell, L. F., Yoon, P. H., Pavan, J., & Gaelzer, R. 2011a, ApJ,727, 16, doi: —. 2011b, PPCF, 53, 085004, doi:
Ziebell, L. F., Yoon, P. H., Petruzzellis, L. T., Gaelzer, R., &Pavan, J. 2015, ApJ, 806, 237,doi:
Ziebell, L. F., Yoon, P. H., Sim˜oes, F. J. R., Gaelzer, R., &Pavan, J. 2014c, PhPl, 21, 010701, doi:10.1063/1.4861619