Verification of particle simulation of radio frequency waves in fusion plasmas

Verification of particle simulation of radio frequency waves in fusion plasmas


Radio frequency (RF) waves can provide heating, current and flow drive, as well as instability control for steady state operations of fusion experiments. A particle simulation model has been developed in this work to provide a first-principles tool for studying the RF nonlinear interactions with plasmas. In this model, ions are considered as fully kinetic particles using the Vlasov equation and electrons are treated as guiding centers using the drift kinetic equation. This model has been implemented in a global gyrokinetic toroidal code (GTC) using real electron-to-ion mass ratio. To verify the model, linear simulations of ion plasma oscillation, ion Bernstein wave, and lower hybrid wave are carried out in cylindrical geometry and found to agree well with analytic predictions.


I Introduction

The importance of radio frequency (RF) waves as a source for heating and current drive has been recognized from the early days of magnetically confined plasma researchFisch (1987); Gormezano et al. (2007). The RF waves provide one of the very few options for steady state operation of the burning plasma experiment ITER, the crucial next step in the quest for the fusion energy. First, the RF waves in ITER will be used to deliver sufficient central heating power to access the H-mode confinement regime and to control the plasma temperature. Secondly they can provide a non-inductive central current drive and an off-axis current drive capability for the current profile control. Thirdly they will be used for the control of magnetohydrodynamic (MHD) instabilities in ITER. It has also been proposedRostoker (2011) that the RF waves can be used for driving plasma flows and current in the field reversed configurationBinderbauer et al. (2010). To effectively utilize the RF power we need a better understanding of the key physics of RF waves in plasmas, e.g., wave-particle interactionKominis (2008); Lee et al. (2012); Hellsten et al. (2012), mode conversionLin et al. (2011); Tsujii et al. (2012) and nonlinear effectsTripathi et al. (1979); Liu et al. (1984); Kuley et al. (2010a, b); Cesario et al. (2011); Goniche et al. (2010); Gao et al. (2011)

Two computational methods have been widely used to study wave-particle interactions in fusion plasmas. The first solves the wave equation derived from the linearized Vlasov-Maxwell system (the full wave model). This approach has been used in the eigenvalue solvers like TORICBrambilla (1999) and AORSAJaeger et al. (2001) to study high frequency waves such as the lower hybrid wave and the ion Bernstein wave. However, this method does not capture the crucial nonlinear physics. The second method is the initial value simulation in which a kinetic equation is integrated in time, retaining all nonlinearities. Such an approach has been taken by gyrokinetic (GK) simulation codes, which have revolutionized studies of turbulent transport driven by low frequency drift wavesLee (1983); Lin et al. (1998). Nonlinear phenomena of the RF waves have been studied in the slab geometries with particle codes such as GeFiQi et al. (2013), VorpalJenkins et al. (2013) and G-gaugeYu and Qin (2009).

For waves in the intermediate frequency range, between the ion and electron cyclotron frequencies (e.g., lower hybrid wave, ion Bernstein wave, etc.), the GK model is not valid, but a fully kinetic model for both ions and electrons is inefficient due to the small electron-to-ion mass ratio. These waves often play important roles in the kinetic processes of magnetized plasmas, e.g., particle acceleration, current drive, plasma heating and spectral cascade of turbulence from long to short wavelength. In this work, we develop a simulation model for these waves, which uses fully kinetic (FK) ions but treats electrons in the drift kinetic approximation (DK). We will study only waves with wavelength longer than the electron gyroradius, so that the electron GK equation reduces to the DK equation. The current FK/DK hybrid simulation modelChen and Parker (2009) can be regarded as a reduced version of the FK/GK modelLin et al. (2005), which overcomes the difficulty associated with the small electron mass by analytically removing the high frequency modes (electron cyclotron frequency and electron plasma frequency). Our goal is to develop a new nonlinear toroidal particle simulation model, which is the most effective approach to study the nonlinear physics in the RF heating and current drive.

Realistic RF simulations for fusion plasmas also require the global toroidal geometry and massively parallel computing due to multiple temporal and spatial scales. The current work utilizes the gyrokinetic toroidal code (GTC)Lin et al. (1998) to take advantage of its existing physics capability, toroidal geometry and computational power. GTC has been extensively applied to study turbulent transport in fusion plasmas including ion and electron temperature gradient turbulence,Lin et al. (1999, 2002, 2007) collisionless trapped electron mode turbulence,Xiao and Lin (2009) energetic particle turbulence and transportZhang et al. (2008, 2011, 2012); Deng et al. (2012) and neoclassical transportLin et al. (1997). As a first step in developing this nonlinear toroidal particle simulation model, the verification of the linear physics of lower hybrid wave (LHW) and ion Bernstein wave (IBW) in cylindrical geometry are presented in this paper.

The paper is organized as follows: the fully kinetic ion and drift kinetic electron simulation model is described in Sec. II, Sec III gives the verification of the GTC simulation of the electrostatic normal modes in uniform plasmas. Sec. IV summarizes this work.

Ii Formulation of fully Kinetic ion and Drift kinetic electron simulation model

ii.1 Formulation of FK ion and DK electron model

The FK ion and DK electron simulation model treats the ion with the fully kinetic (FK) model and the electron with the drift kinetic (DK) approximation. For the FK ion, the dynamics is described by the six dimensional Vlasov equation


where is the ion distribution function, is the ion charge and is the ion mass. is the equilibrium magnetic field. In the current simulation we use the cylindrical coordinates , where is the radial position, is the poloidal angle and is the length of the cylinder with circular cross section. The evolution of the ion distribution function can be described by the Newtonian equation of motion in the presence of self-consistent electromagnetic field as follows


In the fully kinetic version of the GTC code we use for the velocity space, where and are the parallel and perpendicular velocity, respectively, and is the gyro phase angle. This model retains full finite Larmor radius effects and wave frequencies larger than , where is the ion gyrofrequency.

Electron dynamics is described by the drift kinetic equation using guiding center position , perpendicular and parallel velocity as a set of independent variables


where is the guiding center distribution function. The evolution of the electron distribution function can be described by the following equations of guiding center motionBrizard and Hahm (2007)


where (by definition), . The above Eq.(4) is valid only for uniform magnetic field. This electron model is suitable for the dynamics with the wave frequency and , where is perpendicular to the magnetic field, is the electron cyclotron frequency and is the electron gyroradius.

The electrostatic potential can be found from the Poisson’s equation


assuming to suppress the undesirable high frequency electron plasma oscillation along the magnetic field line. Second term on the left hand side corresponds to the electron density due to its perpendicular polarization drift of the electrostatic field. The number densities are defined as the fluid moments of the corresponding distribution function,


Eqs. (2)-(6) are implemented using both non-perturbative (full- and perturbative methods in GTC. We use the simulation for the fully kinetic ion to reduce the particle noise in this work. In the current linear simulation, we assume that the background plasma is uniform in density and temperature. We decompose the ion distribution function into its equilibrium and perturbed part , where . By defining the particle weight for the linear simulation, we can rewrite the Vlasov equation for ion as follows

Figure 1: Coordinate system on the poloidal cross section of a cylinder.

where the second and third terms on the right hand side arise due to the change in the perpendicular energy and the correction of the gyro frequency, respectively. By considering the background plasma as a Maxwellian with the temperature , one can further simplify the weight equation as follows


Similarly the weight equation for the electron in a uniform Maxwellian background with the temperature can be written asHolod et al. (2009)


where for the linear simulation. and are the equilibrium and perturbed distribution function, respectively. Eqs. (8) and (9) are valid only for uniform density and temperature. The parallel component of the electric field can be written as


With a uniform magnetic field one can write down the change in the perpendicular energy as follows


where the particle equations of motion in cylindrical coordinates are


In the fully kinetic version of the GTC code the perpendicular component of the velocity and the gyro phase angle can be calculated from Eq. (2) using the Boris push method Boris (1970); Birdsall and Langdon (2005) . In the following section we will discuss the implementation of the Boris push technique in GTC. However, for the calculation of we use conventional Runge-Kutta method.

Figure 2: Schematic diagram for Boris push method. The first step indicates the addition of the first half of the electric field impulse to the velocity. The red color defines the rotation of the velocity vector in the second step. In the third step, we add the second half of the electric field impulse to the rotated velocity component.

ii.2 Boris push implementation in GTC

The particle push is an important part of the simulation process. Eq. (2) is basically Newton’s second law with the force being the Lorentz force. It is numerically challenging to integrate the particle velocity in the presence of the magnetic field. This problem can be overcome by defining the velocity as suggested by Boris Boris (1970); Birdsall and Langdon (2005). This explicit algorithm is simple to implement, with second order accuracy. It is symmetric to the time reversal, i.e., it preserves the canonical invariantsPenn et al. (2003). The Boris push process can be summarized in the following three steps as described in Fig. 2.

In the cylindrical geometry with magnetic field in the z direction, we decompose the velocity components in the direction perpendicular and parallel to the magnetic field. In the first step we add the first half of the electric field impulse to the velocity vector to obtain a new as


We use (x, y) coordinates to represent (see Fig.1). From Eq. (2) we get




In the second step we consider the rotation of the velocity vector . The vector form of this rotation is given by




where , is also a form of rotation vector T scaled to satisfy that the magnitude of the velocity should remain unchanged during the rotation. Eqs. (16) and (17) together give the rotation of the velocity vector as shown by the red color in the Fig. 2.

In the third step we add the remaining half of the electric field impulse to the rotated vector to obtain


Now we can write down the new and gyro phase angle from and


where is chosen to vary in the range of .

Figure 3: (a) Ion plasma oscillation frequency as a function of normalized wavelength , and its verification with the analytical theory (cf. Eq. (22)), (b) comparison of the electrostatic potential of the ion plasma wave as a function of the normalized radius between analytical theory and GTC simulation.

Iii Verification of Normal modes

In this section we will discuss the electrostatic normal modes with in uniform plasmas and uniform magnetic field. The corresponding dispersion relation can be written as


By considering the uniform Maxwellian background plasma using Eqs. (1) and (2), one can write down the susceptibility asLiu and Tripathi (1986)


where , , , are the electron, ion cyclotron frequencies\textcolorred, respectively\textcolorred. and are the electron and ion Larmor radius\textcolorred, respectively. There are only three electrostatic normal modes in the uniform plasma for , e.g., ion plasma oscillation, lower hybrid wave and ion Bernstein wave.

iii.1 Ion plasma oscillation

Unmagnetized ions and magnetized electrons support the normal mode called ion plasma oscillation when . In the massless electron limit, the ion and electron contributions to the susceptibility can be written as


To verify the fully kinetic ion model, we carried out simulation\textcolorreds for different equilibrium plasma density (i.e., varying the ion Debye length ). Fig. 3(a) demonstrate that for small value of , we can recover , the ion plasma oscillation. In the presence of the finite ion temperature the ion plasma wave will be damped after a few oscillations because of ion Landau damping. During this process the electric field can penetrate up to the ion Debye length. GTC simulation of the ion Debye shielding effect agrees well with the analytic theory [Fig. 3(b)]. In the simulations the boundary conditions for the electrostatic potential are at the inner boundary and constant at the outer boundary. These one-dimensional simulations are carried out using the full- method. The system length is about 10 ion Debye lengths. The number of grid points in radial, poloidal, and parallel direction is Nx=100, Ny=100, and Nz=32, respectively. A total of 4000 particles per cell are used. Initially the particles are loaded uniformly with a Maxwellian velocity distribution. The initial fluctuations are due to the random noise.

iii.2 Lower Hybrid waves

Figure 4: (a) Time history of m=4 lower hybrid wave amplitude excited by the antenna (b) Poloidal mode structure of electrostatic potential and (c) Radial profile of .
Figure 5: Comparison of ion Bernstein wave dispersion relation between the analytical solution of Eq.(27) and the GTC simulations with for the first and second harmonics.

Lower hybrid waves are space-charge waves in the frequency range . In this limit ion motion can be taken to be unmagnetized and the ion susceptibility become\textcolorredsLiu and Tripathi (1986)


Now we consider the finite mass of the electron. For such normal modes in the magnetized plasma with , is dominated by term as


which arises due to the guiding center polarization drift. We implement the electron polarization term in GTC similar to the ion polarization term calculated in the gyrokinetic simulation.

By using Eq. (20) in the limit of , the frequency of the lower hybrid wave is


We use an artificial antenna to excite these modes and to verify the mode structure and frequency in our simulation. The antenna is implemented for the electrostatic potential as follows Zhang et al. (2010)


To find the eigenmode frequency of the system, we carry out the scan with different antenna frequencies and find out the frequency in which the mode has the maximum growth of the amplitude. That frequency is then identified as the eigenmode frequency of the system.

In our simulation the background plasma density is uniform with a uniform temperature. The simulations are all linear and electrostatic. We apply a poloidal mode filter to select only the mode. In this simulation , and . Fig. 4(a) is the time evolution of the (m=4) LHW excited with an antenna frequency , which gives the maximal growth of the wave amplitude. Fig. 4(b) is the poloidal mode structure of the electrostatic potential. The simulation result of the LHW frequency agrees well with the analytical result of (cf. Eq. (25)).

iii.3 Ion Bernstein waves

An important kinetic feature for the normal modes of magnetized ion plasma is the finite Larmor radius effect, which modifies the cold plasma mode with frequency close to the harmonics of ion cyclotron frequency, known as ion Bernstein waves (IBW). Using Eq. (20) the dispersion relation of the IBW becomes


Fig. (5) shows the dispersion relation of the ion Bernstein wave obtained by solving the Eq. (27) analytically with and for the first and second harmonics. To compare our GTC simulation with analytical results we carried out our simulations in different wavelengths and for different harmonics and . Fig. (5) demonstrates a good agreement between the analytical and GTC simulation results of the IBW frequency. These simulations are carried out using the method. The number of grid points in radial, poloidal, and parallel direction is Nx=100, Ny=200, and Nz=32, respectively. A total of 90 particles per cell are used. We have carried out the convergence study of the real frequency as a function of the number of particles per cell (cf. Fig. 6). The simulation results do not depend sensitively on the number of particles, as the grid numbers per wavelength is sufficiently large (100 in this case). In our simulation we have , where is the frequency of the normal mode and is the time step. So, we have more than 600 time steps per wave period. To measure the wave frequency we count the number of time steps in several wave periods from the time history of the wave amplitude. The uncertainty in measuring the frequency is defined as the inverse of number of time steps. This provides a better accuracy than the FFT in measuring the frequency.

Iv Discussions

In summary, with the implementation of the fully kinetic ion and drift kinetic electron model, GTC is particularly applicable to problems in which the electrostatic normal mode frequency ranges from ion Bernstein wave to lower hybrid waves. This new simulation model should have wide applications in the areas of radio frequency heating and current drive, control of MHD instability, and other nonlinear phenomenon. The model is more efficient for the physical process with , , and can handle the realistic electron-to-ion mass ratio, by removing the fast electron gyro motion from the wave dynamics. The LHW and IBW excitation by artificial antenna provides the verification of the mode structure, and the frequency using the predicted by linear theory. Our initial verification suggests that the present simulation model is promising and can incorporate a broad range of realistic issues (toroidal geometry, electromagnetic effects, nonlinear kinetic effects, nonlinear ion Landau damping, parametric instabilities, and ponderomotive effects).

Figure 6: Convergence study of real frequency as a function of number of particles per cell. Red line represents the fit to the data points.
This work is supported by the Trialpha Energy Inc., U. S. Department of Energy (DOE), and China National Magnetic Confinement Fusion Science Program, Grant No. 2013GB111000. Simulations were performed using supercomputers at ORNL, NERSC and NSCC-TJ.


  1. preprint: AIP/123-QED


  1. N. J. Fisch, Rev. Mod. Phys. 59, 175 (1987).
  2. C. Gormezano, A. C. C. Sips, T. C. Luce, S. Ide, A. Becoulet, A. Litaudon, A. Isayama, J. Hobrik, M. R. Wade, T. Oikawa, R. Prater, A. Zvonkov, B. Llyod, T. Suzuki, E. Barbato, P. Bonoli, C. K. Phillips, V. Vdovin, E. Joffrin, T. Casper, J. Ferron, D. Mazon, D. Moreau, R. Bundy, C. Kessel, A. Fukuyama, N. Hayashi, F. Imbeaux, M. Murakami, A. R. Polevoi,  and H. E. StJohn, Nucl. Fusion 47, S285 (2007).
  3. N. Rostoker, Private Comm. (2011).
  4. M. W. Binderbauer, H. Y. Guo, M. Tuszewski, S. Putvinski, L. Sevier, D. Barnes, N. Rostoker, M. G. Anderson, R. Andow, L. Bonelli, F. Brandi, R. Brown, D. Q. Bui, V. Bystritskii, F. Ceccherini, R. Clary, A. H. Cheung, K. D. Conroy, B. H. Deng, S. A. Dettrick, J. D. Douglass, P. Feng, L. Galeotti, E. Garate, F. Giammanco, F. J. Glass, O. Gornostaeva, H. Gota, D. Gupta, S. Gupta, J. S. Kinley, K. Knapp, S. Korepanov, M. Hollins, I. Isakov, V. A. Jose, X. L. Li, Y. Luo, P. Marsili, R. Mendoza, M. Meekins, Y. Mok, A. Necas, E. Paganini, F. Pegoraro, R. Pousa-Hijos, S. Primavera, E. Ruskov, A. Qerushi, L. Schmitz, J. H. Schroeder, A. Sibley, A. Smirnov, Y. Song, X. Sun, M. C. Thompson, A. D. Van Drie, J. K. Walters,  and M. D. Wyman (the TAE Team), Phys. Rev. Lett. 105, 045003 (2010).
  5. Y. Kominis, Phys. Rev. E 77, 016404 (2008).
  6. J. Lee, I. F. Parra, R. R. Parker,  and P. T. Bonoli, Plasma Phys. Control. Fusion 54, 125005 (2012).
  7. T. Hellsten, T. J. Johnson, D. V. Eester, E. Lerche, Y. Lin, M. Mayoral, J. Ongena, G. Calabro, K. Crombé, D. Frigione, C. Giroud, M. Lennholm, P. Mantica, M. F. F. Nave, V. Naulin, C. Sozzi, W. Studholme, T. Tala, T. Versloot,  and J.-E. Contributors, Plasma Phys. Control. Fusion 54, 074007 (2012).
  8. Y. Lin, J. E. Rice, S. J. Wukitch, M. L. Reinke, M. Greenwald, A. Hubbard, Y. Marmar, E.S. Podpaly, M. Porkolab, T. N.,  and the Alcator C-Mod team, Nucl. Fusion 51, 063002 (2011).
  9. N. Tsujii, M. Porkolab, P. T. Bonoli, Y. Lin, J. C. Wright, S. J. Wukitch, E. F. Jaeger, D. L. Green,  and R. W. Harvey, Phys. Plasmas 19, 082508 (2012).
  10. V. K. Tripathi, C. S. Liu,  and C. Grebogi, Phys. Fluids 22, 301 (1979).
  11. C. Liu, V. Tripathi, V. Chan,  and V. Stefan, The Physics of fluids 27, 1709 (1984).
  12. A. Kuley, C. Liu,  and V. Tripathi, Physics of Plasmas 17, 072506 (2010a).
  13. A. Kuley, C. S. Liu,  and V. K. Tripathi, Phys. Plasmas 17, 072506 (2010b).
  14. R. Cesario, L. Amicucci, C. Castaldo, M. Kempenaars, S. Jachmich, J. Mailloux, O. Tudisco, A. Galli, A. Krivska,  and J. contributors, Plasma Phys. Control. Fusion 53, 085011 (2011).
  15. M. Goniche, L. Amicucci, Y. Baranov, V. Basiuk, G. Calabro, A. Cardinali, C. Castaldo, R. Cesario, J. Decker, D. Dodt, A. Ekedahl, L. Figini, J. Garcia, G. Giruzzi, J. Hillairet, G. T. Hoang, A. Hubbard, E. Joffrin, K. Kirov, X. Litaudon, J. Mailloux, R. Oosako, T. Parker, V. Pericoli Ridolfini, Y. Peysson, P. Platania, F. Rimini, P. K. Sharma, C. Sozzi,  and G. Wallace, Plasma Phys. Control. Fusion 52, 124031 (2010).
  16. Z. Gao, N. Fisch,  and H. Qin, Physics of Plasmas 18, 082507 (2011).
  17. M. Brambilla, Plasma Phys. Control. Fusion 41, 1 (1999).
  18. E. F. Jaeger, L. A. Berry, E. D’Azevedo, D. B. Batchelor,  and M. D. Carter, Phys. Plasmas 8, 1573 (2001).
  19. W. W. Lee, Physics of Fluids 26, 556 (1983).
  20. Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang,  and R. B. White, Science 281, 1835 (1998).
  21. L. Qi, X. Wang,  and Y. Lin, Physics of Plasmas 20, 062107 (2013).
  22. T. G. Jenkins, T. M. Austin, D. N. Smithe, J. Loverich,  and A. H. Hakim, Phys. Plasmas 20, 012116 (2013).
  23. Z. Yu and H. Qin, Phys. Plasmas 16, 032507 (2009).
  24. Y. Chen and S. E. Parker, Phys. Plasmas 16, 052305 (2009).
  25. Y. Lin, X. Wang, Z. Lin,  and L. Chen, Plasma Phys. Contr. Fusion 47, 657 (2005).
  26. Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang,  and P. H. Diamond, Phys. Rev. Lett. 83, 3645 (1999).
  27. Z. Lin, S. Ethier, T. S. Hahm,  and W. M. Tang, Phys. Rev. Lett. 88, 195004 (2002).
  28. Z. Lin, I. Holod, L. Chen, P. H. Diamond, T. S. Hahm,  and S. Ethier, Phys. Rev. Lett. 99, 265003 (2007).
  29. Y. Xiao and Z. Lin, Phys. Rev. Lett. 103, 085004 (2009).
  30. W. Zhang, Z. Lin,  and L. Chen, Phys. Rev. Lett. 101, 095001 (2008).
  31. W. Zhang, Z. Lin,  and L. Chen, Phys. Rev. Lett. 107, 239501 (2011).
  32. H. S. Zhang, Z. Lin,  and I. Holod, Phys. Rev. Lett. 109, 025001 (2012).
  33. W. Deng, Z. Lin, I. Holod, Z. Wang, Y. Xiao,  and H. Zhang, Nucl. Fusion 52, 043006 (2012).
  34. Z. Lin, W. M. Tang,  and W. W. Lee, Phys. Rev. Lett. 78, 456 (1997).
  35. A. J. Brizard and T. S. Hahm, Rev. Mod. Phys. 79, 421 (2007).
  36. I. Holod, W. L. Zhang, Y. Xiao,  and Z. Lin, Phys. Plasmas 16, 122307 (2009).
  37. J. Boris, in Proceedings of the fourth international conference on numerical simulation of plasmas (NRL, 1970) pp. 3–67.
  38. C. K. Birdsall and A. B. Langdon, Plasma physics via computer simulation (Institute of Physics, New York, 2005).
  39. G. Penn, P. H. Stoltz, J. R. Cary,  and J. Wurtele, J. Phys. G: Nucl. Part. Phys. 29, 1719 (2003).
  40. C. S. Liu and V. K. Tripathi, Physics Reports 130, 143 (1986).
  41. H. S. Zhang, Z. Lin, I. Holod, X. Wang, Y. Xiao,  and W. L. Zhang, Phys. Plasmas 17, 112505 (2010).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description