Evolution of Nonthermal Particle Distributions in Radio Frequency Heating of Fusion Plasmas

1. Evolution of Nonthermal Particle Distributions in Radio Frequency Heating of Fusion Plasmas Paul Bonoli on behalf of the SciDAC Center for Simulation of Wave-Plasma Interactions SciDAC 2007 Conference Boston, Massachusetts June 25-28, 2007

2. Participants in the Center for Simulation of Wave – Plasma Interactions L.A. Berry, D.B. Batchelor, E.F. Jaeger, E. D`Azevedo, M. Carter D. Smithe C.K. Phillips, E. Valeo N. Gorelenkov, H. Qin P.T. Bonoli, J.C. Wright H. Kohno R.W. Harvey, A.P. Smirnov N.M. Ershov M. Brambilla R. Bilato Politecnico di Torino R. Maggiora V. Lancellotti M. Choi D. D’Ippolito, J. Myra - Lodestar Research

3. Background and Motivation for Work

4. The role of RF power in fusion plasmas has evolved over the years • Early applications involved bulk plasma heating and non-inductive maintenance of the entire plasma current: • Ion cyclotron resonance heating (ICRH), lower hybrid heating (LHH), and electron cyclotron resonance heating (ECRH). • Lower hybrid current drive (LHCD) was especially successful and efficient. • More recent applications have utilized RF waves for localized control of the plasma current and pressure profiles: • LH and EC current drive for sawtooth control, NTM control. • Mode converted ICRF waves for shear flow generation and current generation. • Current profile control for access to high confinement (‘advanced tokamak”) regimes with high bootstrap current fraction (> 70%). • Most of the RF applications listed above are planned for ITER: • A predictive capability is needed to insure success in ITER applications.

5. So what are some of the challenges ? • Complete description of the wave-particle interaction involves integrating several “separate” calculations: • Antenna coupling (linear and nonlinear response) • Wave propagation • Wave absorption • Multiple spatial scales can occur: • ICRF Mode conversion and current drive • RF sheath formation in the edge • RF waves can interact in a nonlinear fashion with thermal and nonthermal particles in the plasma (long timescale): • Nonthermal electron and ion distributions produced by the RF wave itself. • Nonthermal ions already present in the plasma due to neutral beam injection (NBI) and fusion reactions (alpha particles).

6. ICRF Heating Involves Two Important Processes:Ion Cyclotron Absorption with Nonthermal Ion Tail Production Conversion of Fast Wave to Short Wavelength Modes

7. Physics of Nonthermal Ion Tail Evolution • If a “minority” ion species (5% hydrogen) is present in a “majority” ion species plasma (95% deuterium) then an RF wave at the cyclotron frequency of the minority ion will have an electric field component with the same polarization as the minority ion: • Secular interaction - wave damps its power via cyclotron absorption on the minority ion. • Nonthermal, anisotropic minority ion tail is generated that slows down and heats background electrons (drag) and majority ions (collisions). • But this process is greatly complicated by a phenomenon known as “mode conversion”.

8. Dissection of the “minority” heating scheme in a controlled fusion device the tokamak:Simulation and detection of mode converted ion cyclotron RF waves

9. Mode Conversion Changes Polarization of Incoming Fast Wave, Thus Modifying Resonant Ion Cyclotron Absorption

10. Two full-wave solvers have been advanced within the RF SciDac Center Blowup region Slow ion cyclotron wave • All Orders Spectral Algorithm (AORSA) – 1D, 2D & 3D (Jaeger) • Spectral in all 3 dimensions • Cylindrical coordinates (x, y, ) • Includes all cyclotron harmonics • No approximation of small particle gyro radius r compared to wavelength l • Produces huge, dense, non-symmetric, indefinite, complex matrices • TORIC – 2D (Brambilla/Bonoli/Wright) • Mixed spectral (toroidal, poloidal), finite element (radial) • Flux coordinates • Up 2nd cyclotron harmonic • Expanded to 2nd order in r/l • Sparse banded matrices, with dense blocks Electrostatic ion Bernstein wave

11. Understanding the mode conversion aspect of ICRF heating required theory, experiment, and computation • Initial observations of mode converted ICRF waves in Alcator C-Mod presented a scientific conundrum: • Waves were detected on the tokamak LFS and at kR≈7 cm-1 • This was the “wrong” location and wavenumber to be the anticipated ion Bernstein wave (IBW) • But our full-wave simulations (TORIC and AORSA) also revealed the presence of these waves at the “wrong” location and wavenumber.

12. Experimental Observation of a “new” type of mode conversion- the ion cyclotron wave (ICW) A. Mazurenko, PhD thesis, MIT (2001). Nelson-Melby et al, PRL 90 (2003) 155004

13. Both the TORIC and AORSA Solvers also predicted the ICW wave field feature ICW TORIC at 240Nr x 255 Nm AORSA at 230Nx x 230 Ny • In fact the mode converted ICW had been predicted to exist years ago [F. W. Perkins, Nuclear Fusion (1977)] but had been forgotten. IBW FW

14. Predicted RF Electric Field from the TORIC field solver has been used in a synthetic diagnostic code for the PCI TORIC E field is used in synthetic PCI on C-Mod that measures the perturbed density due to mode converted ICRF waves. Y. Lin, A. Parisot, J. Wright; PoP, 2005, PPCF, 2005 • Shear flow generation for pressure profile control is possible in fast wave to ion cyclotron wave mode conversion (Myra, Jaeger et al Phys Rev Lett, 2003 and Melby et al Phys Rev Lett. 2003)

15. Calculations on the Cray XT3 Jaguar have allowed the first simulations of mode conversion in ITER ITER with D:T:HE3 = 20:20:30 with NR = NZ = 350, f = 53 MHz, n =2.5x1019 m-3 (4096 processors for 1.5 hours on the Cray XT-3) E_perp Blowup (E_parallel) Mode converted Ion Cyclotron Wave (ICW) Calculation of flow or current drive for full antenna spectrum in ITER requires petaflop capability

16. Scaling of Full-wave ICRF solvers to > 20,000 processors demonstrated for ICW Mode Conversion in ITER in preparation for ITER mode conversion studies ITER with D:T:HE3 = 20:20:30 with NR = NZ = 500, f = 53 MHz, n =2.5x1019 m-3

17. Dissection of the “minority” heating scheme in a controlled fusion device the tokamak:Simulation and measurements of velocity space structure of nonthermal ion distributions.

18. Wave propagation and the plasma response are governed by the Maxwell-Boltzmann system of equations For time harmonic (rapidly oscillating) wave fields E with frequency ω, Maxwell’s equations reduce to the Helmholtz wave equation: Wave Solvers (AORSA) (TORIC) The plasma current (Jp) is a non-local, integral operator (and non-linear) on the rf electric field and conductivity kernel: SIGMAD Module gives  (f0,s) The long time scale response of the plasma distribution function is obtained from the bounce averaged Fokker-Planck equation: Plasma Response (CQL3D) 0 where Need to solve this nonlinear, integral set of equations for wave fields and velocity distribution function self-consistently. This requires an iterative process to attain self-consistency.

19. Calculation for C-Mod minority H, NR = 128, NZ = 128,[256 processors for 3 hrs on Cray XT3 – ORNL]

20. 2D Field and Dissipation Contours Show that Heating is Concentrated at Ion Turning Points on the Minority Resonance Chord at r/a  0.45 Wave fields Heating (H)

21. H+ ~0.1% B Experimental measurements of the energetic ion tail on C-Mod have been made using a compact neutral particle analyzer Good agreement between simulated and measured tail temperature Tion 70 keV Courtesy of V. Tang, PhD Thesis, MIT (2006);

22. 3D (r, V , V//) distribution function from CQL3D – AORSA reproduces CNPA measurements using a synthetic code diagnostic Building this synthetic diagnostic required a close collaboration between theory and experiment (V. Tang and R. Harvey) Courtesy of V. Tang, PhD Thesis, MIT (2006); also PPCF, 49, 873 (2007).

23. Outstanding challenges for the near future:Finite ion drift orbit effects

24. Simulations of high harmonic fast wave (HHFW) – fast ion beam interaction in DIII-D are still unresolved • DIII–D high density L-mode Stronger Beam Interactions at 4D (60 MHz) Than at 8D (116 MHz) Observed in DIII-D CQL3D-AORSA predicts increased absorption as frequency was raised – in disagreement with expt. Monte Carlo ORBIT code (ORBIT-RF) combined with an RF operator (using fields from TORIC solver) does reproduce the experimental trend. 1014 (/s) Neutron reaction rate Power Sn: neutron enhancement factor 217-05/MC/jy

25. Orbit topology modifies wave-particle resonance • Shown at right are trajectories for 12 particles in the C-Mod case: • 4 equi-spaced || velocities • 3 equi-spaced  velocities • 409,600 complete poloidal orbits • Particle cyclotron resonances and strong quasilinear diffusion occur in roughly vertical planes in zero-orbit width description. • But orbit topology can move particles away from (or towards) resonances that would be sampled (not sampled) in full-wave solver.

26. We are investigating finite ion drift orbit effects using two approaches: • The diffusion coefficient (D) has been evaluated by a direct orbit integration using electric fields from AORSA: • The “DC” code computes averages of the changes in velocity, pitch angle, and radial position over a complete bounce orbit, to obtain a set of RF induced diffusion coefficients. • Diffusion Coefficient calculations done on CRAY XT3 (ORNL) using 256 processors @ 10 min. • The Monte Carlo code ORBIT RF has been combined with the TORIC ICRF solver: • Self-consistent iteration not yet carried out. • ORBIT RF code ported to JAGUAR where good scaling to > 1000 processor cores has been demonstrated. • We are now examining best way to pass statistical distribution from ORBIT RF to TORIC & AORSA to do self-consistent iteration.

27. Monte Carlo ORBIT Code has been coupled to the TORIC full-wave solver through an RF Operator Finite orbit effects can be studied quantitatively using this approach: QL Diffusion Operator Formulated in terms of Multi-Fourier Poloidal Modes from the TORIC ICRF Solver and used to compute increment in magnetic moment due to the ICRF interaction: Numerical distribution function from ORBIT-RF not yet coupled back to TORIC or AORSA – this is an important and difficult next step !

28. Beam pressure computed with f(E) from ORBIT-RF agrees qualitatively with experiment – but iteration is still needed between the ORBIT code and full-wave solver Particle distribution : f(E) ORBIT-RF 104 4th harmonic NB only 8th harmonic 116MHz (1.7MW) Beam injection energy (80keV) Sn(ORBIT-RF)=1.2 60 MHz (0.8MW) Sn(ORBIT-RF)=1.9 Beam Ion Pressure (N/m 2 ) r/a Energy (keV) 217-05/MC/jy

29. Outstanding challenges for the near future:Nonlinear effects at the RF antenna-tenuous edge plasma

30. ICRF Launchers in Contact with Plasma are Subject to Nonlinear Effects, Leading to Parasitic Power Losses. RF sheaths can form due to mismatch between equilibrium B and the antenna structure (E//), resulting in power dissipation: Alcator C-Mod Dipole Antenna Electrons are preferentially accelerated out of the sheath region. A DC voltage Vrfis set up to maintain ambipolarity.

31. Capability to efficiently compute 3D wave fields will be important for assessing antenna – edge interaction, especially in weak single pass damping regime NSTX ITER Scenario 2 NSTX simulation summed over 81 toroidal modes.ITER simulation summed over 169 toroidal modes)[AORSA run on JAGUAR using 2048 processors for 8 hrs]

32. Two Approaches are Being Pursued to Study the Nonlinear ICRF antenna – edge Interaction • Implementation of RF sheath boundary conditions in full-wave solver (spectral solution): • Start with linear field response from a coupled full-wave field solver (TORIC) and 3D electromagnetic antenna code (TOPICA). • Modify metal wall BC in field solver to include sheath dissipation and then iterate with antenna code. • Approach will quantify how much ICRF power is coupled to the plasma. • Time domain simulations using 3D EM field solver - VORPAL • Fully implicit time domain dielectric response module has been implemented for electrons and ions. • Use PIC-treatment in future for ion response (fully nonlinear).

33. VORPAL Time Domain Simulation of Antenna LoopD. Smithe, Poster at this Conference (Tuesday – Evening)

34. VORPAL Time Domain Simulation of Antenna Loop Surface E-field on Loop antenna RF B-field Surface E-field Wavefronts Radiation Pattern Radiation from Behind

35. Summary • Combined full-wave and Fokker Planck solvers (CQL3D-SIGMAD-AORSA) have been used to simulate minority ICRF heating in present day tokamaks and in ITER: • Full coupling uses self-consistent nonthermal ion distributions and quasilinear diffusion coefficient in differential form. • Comparison with synthetic diagnostic (CNPA) and experiments validates the use of this simulation to predict these experiments in burning plasmas such as ITER. • Both AORSA and TORIC can simulate ICW/IBW mode conversion in present day tokamaks and in ITER: • Comparison of synthetic diagnostic (PCI) with experiment validates predictive capability of simulation: • Can now assess use of mode converted waves for pressure and current profile control in present day tokamaks and for ITER.

36. Summary • Fully self-consistent coupling of our full-wave solvers to a Monte Carlo orbit code (ORBIT RF) is underway: • Full coupling will use statistical nonthermal ion distributions and quasilinear diffusion coefficient in differential form • A self-consistent treatment of the RF antenna – edge plasma is underway: • Linear antenna coupling problem is substantially completed using the TORIC – TOPICA suite. • Boundary conditions for sheaths are being implemented in our full-wave solver (TORIC) • Proof of principle time domain simulations of sheath formation have been done using VORPAL and will now be extended to the nonlinear regime using PIC treatment for ions.