Self-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas
SSelf-consistent kinetic simulations of lower hybriddrift instability resulting in electron current drivenby fusion products in tokamak plasmas
J W S Cook , S C Chapman and R O Dendy , Centre for Fusion, Space and Astrophysics, Department of Physics, WarwickUniversity, Coventry CV4 7AL, U.K. Euratom/CCFE Fusion Association, Culham Science Centre, Abingdon, OxfordshireOX14 3DB, U.K.
Abstract.
We present particle-in-cell (PIC) simulations of minority energetic protonsin deuterium plasmas, which demonstrate a collective instability responsible foremission near the lower hybrid frequency and its harmonics. The simulations capturethe lower hybrid drift instability in a regime relevant to tokamak fusion plasmas,and show further that the excited electromagnetic fields collectively and collisionlesslycouple free energy from the protons to directed electron motion. This results in anasymmetric tail antiparallel to the magnetic field. We focus on obliquely propagatingmodes under conditions approximating the outer mid-plane edge in a large tokamak,through which there pass confined centrally born fusion products on banana orbits thathave large radial excursions. A fully self-consistent electromagnetic relativistic PICcode representing all vector field quantities and particle velocities in three dimensionsas functions of a single spatial dimension is used to model this situation, by evolving theinitial antiparallel travelling ring-beam distribution of 3MeV protons in a background10keV Maxwellian deuterium plasma with realistic ion-electron mass ratio. Thesimulations thus demonstrate a key building block of alpha channelling scenarios forburning fusion plasmas in tokamaks. a r X i v : . [ phy s i c s . p l a s m - ph ] O c t elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas
1. Introduction
There is sustained interest in the collective instabilities of energetic ion populations inplasmas with confining magnetic fields that involve emission and absorption of radiofrequency (RF) waves. For fusion plasmas, the primary aspects of practical interestinclude the suprathermal ion cyclotron emission arising from the products of deuteriumand tritium fusion reactions, observed in both JET and TFTR, typically in the > ci (and its harmonics), and the dynamics of wavesinvolving collective motion of the electrons, for which the lowest natural frequency isthe lower hybrid frequency [3], ω LH , defined by ω LH = (cid:18) Ω ce Ω ci ce Ω ci /ω pi (cid:19) / (1)Here Ω ce , Ω ci and ω pi represent the electron and ion cyclotron frequencies and the ionplasma frequency respectively. Fusion alpha-particles are born at energies of 3.5MeVin deuterium-tritium fusion plasmas whose thermal energies are of order 10keV. Thiscreates the possibility of population inversion in velocity space, at least transientlyduring burn initiation before the alpha-particles have slowed through collisions withelectrons. Population inversion opens the door to rapid energy exchange throughcollective instability. It is already known that, in tokamak plasmas, spatially localisedinversions of the energy distribution of fusion-born ions can be sustained for long periodsas a result of the interplay between particle energy and the pitch angle-dependentcharacter of particle drift orbits. Specifically, ion cyclotron emission was excited in theouter mid-plane edge region of JET and TFTR as a consequence of the distinctive radialexcursions [4,5] of the drift orbits of fusion products born with pitch angles just inside thetrapped-passing boundary. This was first seen for fusion protons born in pure deuteriumplasmas in JET [4, 6], and subsequently for fusion alpha-particles in deuterium-tritiumplasmas in JET [7,8] and TFTR [1, 5,9], as well as injected beam ions [10] in TFTR andheated minority ions in JET [8]. For reviews and further theory, we refer to [1, 11, 12].These observations motivate the choice of model distribution for the energetic ions inthe present study. Given the variety of possible alpha channelling physics scenarios, itis necessary to find a focus for the representation of population inversion in the present elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas − energetic minority ions, thermal majorityions, and thermal electrons represented with realistic mass ratios. The results includedin this paper expand and elaborate upon the mechanism identified in Ref. [19] for thegeneration of a suprathermal tail in the distribution of electron velocities parallel tothe magnetic field and consequently provides a building block for the construction oftokamak-relevant alpha channelling scenarios.The linear phase of the instability underlying the phenomena seen in these fullynonlinear numerical simulations is consistent with a lower hybrid drift instability(LHDI). To simplify and generalise somewhat, the wave excited is predominantlyelectrostatic and propagates nearly perpendicular to the magnetic field, with a parallelphase velocity which (in the cases we consider) is sufficiently close to the electron thermalvelocity to enable significant electron Landau damping. The wave frequency is high incomparison to the ion cyclotron frequency, so that the action of the energetic ions onthe timescale of relevance is predominantly that of a quasi-unmagnetised beam.The scenario we consider differs fundamentally from that considered in much of theLHDI literature in one important aspect. Most previous work on LHDI has considereda free energy source residing in bulk ion drifts that are consistent with, and contributeto, the overall equilibrium of the plasma; notably, diamagnetic drifts. For the tokamakapplication, however, it is necessary to consider the LHDI of diffuse minority energeticion populations that do not contribute to the plasma equilibrium. This is the caseboth for the edge localised population that appears to drive ion cyclotron emission,and for many alpha channelling scenarios. This degree of freedom arises from theinterplay between particle energy and particle orbits in tokamak magnetic geometry,which ultimately derives from the combination of the particle toroidal velocity withthe toroidal component of the magnetic vector potential in the canonical toroidal elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas et al [23] considerthe LHDI using imposed drifts from field inhomogeneities for the strictly perpendicularcase. Hsia et al [22] includes finite k (cid:107) in their analytical approach. The more recentstudy by Silveira et al [32] builds on these studies, and other important works toonumerous to mention by name here. They present a unified formulation for instabilitiesin the lower hybrid range of frequencies due to field inhomogeneities. In these casesthe LHDI nevertheless is excited by inhomogeneity induced drifts that are slower, indimensionless terms, than the fusion product drifts considered in this paper. Thepioneering computational study by Chen and Birdsall [34] has cold fluid electrons.Further numerical work by Yoon et al [21] investigates the LHDI with gyro-averagedelectrons and kinetic ions, while work by Daughton et al [31] uses physical mass ratiokinetic particle simulations for ∇ n and ∇ B inhomogeneities which generate low driftvelocities.Finally we note that PIC codes have been used to address astrophysical problemsthat involve the same combination of effects that is considered here: beam-type minorityion distributions, resulting in RF wave excitation that leads on, through damping, toelectron acceleration [35–38]. PIC codes are also widely used to study other fundamentalinstabilities. Population inversion [39], kink instability [40], Weibel instability [41],Kelvin-Helmholtz instability [42], shocks [43, 44] and wave-particle interactions [45–49]are among the many avenues investigated using PIC.
2. Simulations
The 1D3V PIC code epoch1d, which is based on PSC [50], represents the fulldistributions of all the species in the simulation, using computational macroparticlesin three dimensional velocity space and in one configuration space dimension. These areevolved in time via the relativistic Lorentz force law in the absence of collisions. From elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas x ) and time.The spatial domain is split into N G cells of width ∆ x , each of which has a grid pointat the centre. Cell size, which is uniform across the domain, is set at 1 /
10 the electronDebye length and resolves the gyroradii of all species. The phase space density function f of each species is represented by a large (682670) set of computational macroparticleswhich are initially distributed evenly in space amongst 2048 cells. Although allcomponents of velocity are evolved, only the velocity component in the direction of thesimulation domain causes the particles’ position to change. Each velocity componentof all particles contributes to the respective current component on the grid points innearby cells. The field and particle advance timestep, ∆ t , divided by the cell size, ∆ x ,resolves the speed of light, and thus resolves cyclotron frequencies. Periodic boundaryconditions are used for particles and fields and any wavemodes present propagate parallelor antiparallel to the direction of the simulation domain.The trajectories of background thermal electrons and deuterons, both of whichhave non-drifting Maxwellian distributions with temperature of 10keV, are evolved intime alongside the minority energetic protons (initially at 3MeV) which initially have aring-beam velocity distribution modelled by f p = 12 πv r δ ( v || − u ) δ ( v ⊥ − v r ) (2)see for example Dendy et al. [13] and McClements et al. [9]. Here u is the velocityalong the magnetic field and v r is the perpendicular velocity of the ring. As mentionedpreviously, this choice of model is motivated by multiple studies of ion cyclotron emissionin JET and TFTR. For this reason the initial pitch angle of the ring distribution (anglebetween u and v r ) is initialised at 135 ◦ with respect to the background magnetic fieldat the start of the simulation. The initial velocity distributions of all the species in thesimulation are uniform in space. The angle of inclination of the direction of spatialvariation x , of our simulation to the initial background magnetic field is θ = 84 ◦ .We simulate a quasi-neutral plasma with an electron number density n e = 10 m − ,and an applied magnetic field of B = 3 T . The electron beta is β e = 0 . v p ( t = 0) ≡ ( u + v r ) / , to the Alfv´en speed V A , is v p ( t = 0) /V A (cid:39) .
51. This set of model parameters is selected to enableus to address quasi-perpendicular propagating modes under conditions approximatingthe edge plasma of a large tokamak, subject to the necessary simplifications. Itoptimizes the grid size with respect to the gyroradii of the particles, so as to resolvethe essential physics with reasonable computational resources. The gyroradius of theenergetic protons, and the gyroradii at thermal speeds of electrons and deuterons, areresolved by 965, 1.8 and 111 grid cells respectively. The ratio of protons to deuterons n p /n d = 10 − . The 3MeV ring of protons is not replenished as the distribution function elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas ∼
10 the value anticipated for various energetic ion populations in next stepfusion plasmas, see for example Table 1 of S. Putvinski (Nucl. Fus. , 1275 (1998)).This value is necessary to drive the instability on an acceptable timescale given finitecomputational resources, and brings this study closer to fusion relevant regimes thanmuch of the LHDI literature.
3. Flow of energy between particles and fields
In order to examine the instability produced by the free energy in the energetic protonpopulation we first consider the time evolution of the electric ε E and magnetic ε B totalfield energies in the simulation: ε E = ε N G (cid:88) i =1 E i ∆ x ε B = 12 µ N G (cid:88) i =1 ( B i − B )∆ x (3)where B is the applied magnetic field. Figure 1(a) shows that the electric field energy, ε E , (light blue trace) starts rising exponentially after a period of ∼ τ LH . The magneticfield energy, ε B (magenta trace), rises at approximately the same rate as the electricfield, but contains nearly an order of magnitude less energy, implying that the wavesproduced by the instability are largely electrostatic. The total kinetic energy of species i is defined by ε i = w i N P,i (cid:88) j =1 ( (cid:113) p j c − m i c − m i c ) (4)where w i is the weighting factor to map the energy of each computational macroparticle, j , to the energy of the number of physical particles it represents. Particle weight isdefined as w i = n ,i L/N
P,i where n ,i is the species number density, L is the lengthof the spatial domain, and N P,i is the number of computational macroparticles thatrepresent the species in the simulation. The essential features of figure 1(a) indicate alinear phase of field growth during 10 (cid:46) t/τ LH (cid:46)
15, accompanied by a correspondingsmall decline in total proton energy. Field amplitudes remain approximately constantafter t (cid:38) τ LH , while proton energy continues to decline and is matched by an increasein electron energy.We now define a quantity which better shows the variations in kinetic energy of theproton population, which we denote as fluctuation energy ˜ ε summed over the j = 1 − N P,i computational macroparticles of each i th species:˜ ε i ( t ) = N P,i (cid:88) j =1 ( | ε j ( t ) − < ε i ( t = 0) > | ) (5) elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas ε is the sum of the modulus of the change in kinetic energy of each particlefrom the initial rest frame mean energy value of the species, < ε i ( t = 0) > . Thus ˜ ε i captures kinetic energy dispersion regardless of sign, whereas a time evolving change intotal kinetic energy, ∆ ε i = ε i ( t ) − ε i ( t = 0) where ε i ( t = 0) = (cid:80) N P,i j =1 ε j ( t = 0), onlyrepresents changes in the combined energy of the ensemble of particles. In particular weuse ˜ ε to quantify the initial energy dispersion of the protons during the early phases ofthe instability, see figure 1(b) and the discussion in section 5, because in some respectsit is less susceptible to noise.Figure 1(b) shows in greater detail the species energy dynamics. The values ofproton fluctuation energy ˜ ε p , and of the changes ∆ ε e and ∆ ε d in total electron anddeuterium kinetic energy are plotted alongside changes in ε E and ε B . The value of ˜ ε p grows, ultimately increasing by three orders of magnitude, whereas figure 1(a) showsthat the total proton kinetic energy declines by much less than one order of magnitude.This reflects the role of the proton population as the source of free energy. The changein electron kinetic energy ∆ ε e rises in close alignment with electric field energy ε E during the linear stage of the simulation, 10 < t/τ LH <
15. This is evidence of theelectron oscillation within the largely electrostatic waves excited by the instability. Thecorresponding effect for the deuterons is also visible in figure 1(b) but is much less dueto their higher mass. −2 Time [2 π / ω LH ] E ne r g y [ J ] ε e ε d ε p ε E ε B (i) (ii) (iv)(iii) (a) −6 −4 −2 Time [2 π / ω LH ] E ne r g y [ J ] ∆ε e ∆ε d ε p ε E ε B (iv)(iii)(ii)(i) ~ (b) Figure 1.
Panel (a): Time evolution of the kinetic energy, ε i defined by equation4, in each plasma species, and of the energy of electric field ε E , and magnetic field ε B defined by equation 3. The energy of the magnetic field excludes the appliedcomponent. Panel (b): Time evolution of the proton fluctuation energy, ˜ ε i defined byequation 5, change in kinetic energy of electrons and deuterons, ∆ ε i = ε i ( t ) − ε i ( t = 0),and of the electric and magnetic field energies ε E and ε B . Both panels show verticallines labelled ( i ) − ( iv ) which denote the four snapshots in time at which particle datais shown in other figures. Time is measured in units of the lower hybrid period. elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas −1 −0.5 0 0.5 110 Momentum [p p (t=0)] P r obab ili t y D en s i t y p x p y p z, ⊥ p || (p ⊥ ,12 +p ⊥ ,22 ) (a) −6 −4 −2 0 2 410 −5 −4 −3 −2 −1 Parallel Momentum [p e (t=0)] P r obab ili t y den s i t y (iv)(iii)(ii)(i) (i)(iv)(i) (iv) (b) τ LH ] S k e w ne ss , γ ∆ ε [ J ] ∆ε e γ e (i) (iv)(ii) (iii) (c) τ LH ] J || [ A r b . U n i t s .] J e (v || > −2.5v e (t=0))J e (v || < −2.5v e (t=0)) (iv)(iii)(ii)(i) (d) Figure 2.
Panel ( a ): Initial proton momentum PDF: along the simulation domain, p x ;in the plane created by the magnetic field and the simulation domain, p y ; perpendicularto both magnetic field and simulation domain p z ≡ p ⊥ ; parallel to the magnetic field p (cid:107) ; and total perpendicular p ⊥ ≡ ( p ⊥ , + p ⊥ , ) / . Panel (b): Electron PDFs inrelativistic parallel momentum p (cid:107) , for the four snapshots in time ( i ) − ( iv ) defined inFig. 1, also marked by vertical lines in panels ( c ) and ( d ). Vertical lines indicatethe phase velocities of the forward and backward travelling dominant waves. Panel(c): Skewness γ e of the electron v (cid:107) PDF (black circles) and change in electron kineticenergy ∆ ε e = ε e ( t ) − ε e ( t = 0), (red crosses). Panel (d): Time evolving parallelcurrent in the electron bulk v x > − . v e ( t = 0) (red triangles), and in the electron tail v x < − . v e ( t = 0) (black stars). Momentum axes are in units of the initial relativisticmomentum p i ( t = 0) = (cid:104) (1 /N P,i )Σ N P,i j =1 p j . p j ( t = 0) (cid:105) / of the species i , and time axesare in units of the lower hybrid oscillation period. Figure 2(a) shows projections in momentum space of the PDF of the protonpopulation at t = 0 defined by (2). The ring-beam distribution, as aligned with respectto the magnetic field, is projected into the cartesian co-ordinates of the simulation where x is along the simulation domain, z is perpendicular to both the simulation domain andthe magnetic field and y completes the orthogonal set. The initial PDF is a deltafunction in both p (cid:107) and p ⊥ ≡ ( p ⊥ , + p ⊥ , ) / . The three momentum components areprojected into the simulation domain, in which the ring distribution becomes a more elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas p x (red circled trace) are not equal and opposite.This has consequences for the wave-particle interaction which are dealt with in the nextsection.Figure 2(b) shows electron parallel momentum at four different snapshots in time:( i ) ≡ . τ LH , ( ii ) ≡ . τ LH , ( iii ) ≡ . τ LH , ( iv ) ≡ . τ LH . Henceforth, thesesnapshots in time will be referred to by their Roman numeral.Figures 2(b) and 2(c) show that for t > τ LH the electron distribution developsan asymmetric tail in v (cid:107) with corresponding skewness, γ ( t ) = 1 /N P,i Σ N P,i j =1 (cid:32) v (cid:107) ,j ( t ) − < v (cid:107) ,i ( t ) >σ v (cid:107) ,i ( t ) (cid:33) (6)where < v (cid:107) ,i ( t ) > and σ v (cid:107) ,i ( t ) are respectively the mean and standard deviation of electron v (cid:107) at time t .Figure 2(d) shows the time evolving parallel electron current summed acrossconfiguration space. We plot the current for the bulk ( v (cid:107) > − . v e ( t = 0)) and tail( v (cid:107) < − . v e ( t = 0)) of the electrons separately. We see that approximately equaland opposite currents are created by the bulk and tail electron populations in thiscollisionless simulation. This is a consequence of the periodic boundary conditions usedin the code. The rate of change of current across the simulation domain under periodicboundary conditions follows from Amp`ere’s and Faraday’s laws: µ ∂∂t (cid:90) L J dx = − (cid:90) L ∇×∇× E dx − µ (cid:15) ∂ ∂t (cid:90) L E dx (7)where in one dimension ∇ = ( ∂/∂x, ,
0) and ∇ × ∇× ≡ (0 , − ∂ ∂x , − ∂ ∂x ). Periodicboundary conditions enforce F (0 , t ) = F ( L, t ) and ∇ F (0 , t ) = ∇ F ( L, t ) where F is anycomponent of the electromagnetic field. Equation (7) is a wave equation in the plasmamedium. For current density integrated over the box (cid:18) c ∇×∇× − ∂ ∂t (cid:19) (cid:90) L J dx = 1 /(cid:15) ∂∂t σ (cid:90) L J dx (8)where J = σ E and σ is the plasma conductivity tensor. Consequently, any currentdensity integrated across the domain may oscillate about zero but cannot grow in timefor these simulations. We also note that collisions would differentially affect the bulkand tail electron populations, and hence their relative drift; however our PIC code doesnot incorporate collisions.The asymmetry in the electron parallel PDF, reflecting net directional electronacceleration, continues to grow in the period beyond t = 20 τ LH , during which waveenergy quantified by ε E and ε B is approximately constant, see figure 1(a). We inferLandau damping of the excited waves on resonant electrons, which results in theflattening of the electron parallel momentum PDF in figure 2(b). The vertical lines inthis panel indicate the electron nonrelativistic momentum that corresponds to the phasevelocity of the dominant relevant features in the ω, k transform of the electromagneticfield, to which we now turn. elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas
4. The nature of the electromagnetic fields
Figure 3 displays the spatio-temporal fast Fourier transform (FFT) of the electricfield component along the simulation domain E x , for three epochs of the simulation,showing temporal evolution of the ω, k structure. Referring to figure 1, the threeepochs correspond to the following phases of electromagnetic field evolution: initial,5 ≤ t/τ LH ≤
10; early linear growth, 10 ≤ t/τ LH ≤
15; and late linear growth,15 ≤ t/τ LH ≤
20. The peaks in intensity at k (cid:39) ± ω pe /c and ω (cid:39) ω LH in allthree panels of figure 3 reflect the resonant interaction of the energetic initial protonpopulation with the normal mode of the background plasma to which they most easilycouple, namely the electron-deuteron lower extraordinary-Whistler wave mode whoseanalytical characteristics are outlined below. elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas Wavenumber [ ω pe /c] F r equen cy [ ω L H ] −10 −5 0 5 100102030 −50510 (a) Wavenumber [ ω pe /c] F r equen cy [ ω L H ] −10 −5 0 5 100102030 −50510 (b) Wavenumber [ ω pe /c] F r equen cy [ ω L H ] −10 −5 0 5 100102030 −50510 (c) Figure 3.
Spatio-temporal FFTs of the electric field component in the direction ofthe simulation domain E x , integrated over all x and the following time periods. Panel(a): 5 ≤ t/τ LH ≤
10. Panel (b): 10 ≤ t/τ LH ≤
15. Panel (c): 15 ≤ t/τ LH ≤
20. Overplotted on all panels is the cold plasma dispersion relation (dash trace).Overplotted on panel (a) and panel (b) are velocities corresponding to the peaks in theinitial proton v x PDF (dot-dash trace). Overplotted on panel (c) are the velocities ofthe peaks and extremes in the proton v x PDF at snapshot ( iv ) (dot-dash trace). Thefrequency axes are in units of lower hybrid frequency and wavenumber axes in inverseelectron skin depths ω pe /c . Figure 3(a) shows the waves present in the simulation during its initial phasebefore the onset of linear growth of electric and magnetic fields, in the time interval5 < t/τ LH <
10. Overplotted is the cold plasma electron-deuteron dispersion relation,implicitly defined by (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε − N − iε iε ε − N cos ( θ ) N sin ( θ ) cos ( θ )0 N sin ( θ ) cos ( θ ) ε − N sin ( θ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 (9)This is the determinant of the cold plasma dielectric tensor for arbitrary propagationdirection θ with respect to the magnetic field, where the refractive index N = ck/ω .The components of the dielectric tensor are elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas ε = 1 + ω pe Ω ce − ω + ω pd Ω cd − ω ε = Ω ce ω ω pe Ω ce − ω + Ω cd ω ω pd Ω cd − ω ε = 1 − ω pe ω (10)where ω pi and Ω ci are the plasma and cyclotron frequencies of species i .The cold plasma dispersion relation shows predominantly perpendicular propagat-ing normal wavemodes; ordinary and extraordinary. However some properties of par-allel propagating waves (Whistler) are additionally present, because the angle of wavepropagation is oblique. The result is an extraordinary-Whistler branch with a cut-offfrequency > ω LH . The instability is still essentially a lower hybrid drift-type instability,although it is active at frequencies above ω LH .Figure 3(b) is obtained using data from the interval 10 < t/τ LH <
15. Figure1(a) shows that this corresponds to the early linear stage of electric field wave growth,during which the changes in energy of electrons and deuterons also rise in tandem,reflecting their coherent participation in wave motion. The trace corresponding tosnapshot ( ii ) in figure 2(b) shows that a monotonically decreasing electron tail is drawnout during this epoch, with skewness becoming apparent (figure 2(c)) at its end. Figure3(b) enables us to identify the ω, k character of the waves excited in the linear growthphase of the simulation. Both forward and backward propagating waves are present,at 6 ω LH (cid:46) ω (cid:46) ω LH and 2 ω pe /c (cid:46) | k | (cid:46) ω pe /c . Figure 3(b) also exhibits twopeaks in electric field intensity which appear to be backward propagating harmonicsof the fundamental activity at frequency ω (cid:39) ω LH . During the late linear phase15 ≤ t/τ LH ≤
20, these multiple higher harmonics become established.Figure 3(c) shows data from 15 < t/τ LH <
20 which encompasses the end of thelinear growth phase during which electron acceleration arises (figures 2(b) and 2(c)). Thedominant Fourier structures of the electric field are similar to the linear phase but withmore pronounced backward propagating harmonics. Traces corresponding to snapshots( iii ) and ( iv ) in figure 2(b) show that the electron p (cid:107) PDF becomes increasingly flattenedaround the resonant momentum during this stage.From figure 3(b) we see that the peaks in energy in ω, k space correspond to theregion around ω = 6 ω LH and k = ± ω pe /c . The precise values of ω, k of the dominantwaves correspond approximately to where the phase velocity of the mode matches thevelocities of the peaks in the proton v x PDF. Thus we can infer the parallel phasevelocity of forward v + , and backward v − , travelling waves by reading off coordinates ofthe peaks in power from the ω, k plots. Thus the electrons that are in resonance withthese waves have forward and backward parallel momenta of p (cid:107) , + = − v + m e / cos θ (cid:39) . p e ( t = 0) elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas p (cid:107) , − = − v − m e / cos θ (cid:39) − . p e ( t = 0) (11)in the nonrelativistic limit respectively. These momenta are indicated in figure 2(b) byvertical lines. These show that the instability driven by energetic protons, for thesesimulation parameters, generates waves with phase velocity resonant with the tail of theinitial Maxwellian distribution of electron velocities, where both the gradient and thesmall number of particles are conducive to efficient Landau damping. This gives rise toelectron acceleration and a strongly asymmetric velocity tail: hence current drive.
5. Proton velocity space dynamics
The particle-in-cell simulation presented here provides a full kinetic description of thespecies evolution in configuration space and velocity space. This differs from, andcomplements, analytical approaches. In the linear phase, analytical treatments aretypically based on a small amplitude single harmonic perturbation. The linear phasehere involves, as we have seen, fields mediating subtle interactions between the energeticprotons, background deuterons, and electrons, in ways that would be difficult to capturein detail analytically. In the present section we therefore focus on the kinetics of theenergetic proton population that drives the system, thereby obtaining a fresh perspectiveon the character of the instability during its lifetime.Figure 4 provides snapshots of the proton distributions in velocity space at twoinstants: ( ii ), during the linear phase (upper panels); and ( iv ), at the end of the linearphase (lower panels). The two left panels in figure 4 focus on velocity space. They showthe proton velocity distribution plotted with respect to three axes defined by v (cid:107) and twoorthogonal velocity components v ⊥ , and v ⊥ , in the directions perpendicular to B . Dataplotted are for protons from a distinct region in configuration space as described below.The two right hand panels focus on configuration space behaviour. Each panel comparestwo subsets of protons. Each subset is localised in a distinct region of configuration spaceof extent δx , where δx (cid:28) λ and λ = 2 π/ (2 ω pe /c ) is the wavelength of the dominantelectromagnetic wave. The two distinct regions are separated in configuration space by ∼ λ/
2. To achieve this, protons are binned into subsets in configuration space such thatthe ratio of the size of each bin to the wavelength of the dominant electric field wave is δx/λ (cid:39) . δx = 0 . λ .The upper panels of figure 4 show early linear phase behaviour of the proton ringbeam. There is a distinct wave-like structure oscillating in v (cid:107) around the gyro-orbit,shown in panel (b) for two subsets of protons. The two datasets show velocity spacewaves that are in anti-phase with each other, consistent with resonance with a wave atlocations λ/ elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas −0.71 −0.705−0.500.5−0.500.5 v || [v p (t=0)]v ⊥ ,1 [v p (t=0)] v ⊥ , [ v p ( t = ) ] −2 0 2−0.71−0.704 Gyrophase α [rad] v || [ v p ( t = ) ] −0.72 −0.7−0.500.5−0.500.5 v || [v p (t=0)]v ⊥ ,1 [v p (t=0)] v ⊥ , [ v p ( t = ) ] α [rad] v || [ v p ( t = ) ] ∆ x < 16160 Left panels: Velocities in magnetic field-aligned coordinates of protonsfound in the configuration space region 0 < x/ ∆ x < 16, where ∆ x is the spatial extentof one grid cell. Colour indicates speed in units of v p ( t = 0). Right panels: Proton v (cid:107) asa function of gyrophase α = arctan v ⊥ , /v ⊥ , for protons from two configuration spaceregions: 0 < x/ ∆ x < 16 (blue dots) and 160 < x/ ∆ x < 176 (red pluses). Top panelsshow data from snapshot ( ii ), bottom panels from snapshot ( iv ). Velocity axes are inunits of the proton characteristic velocity, v p ( t = 0) = (cid:104) (1 /N P,p )Σ N P,p j =1 v j . v j ( t = 0) (cid:105) / ,where N P,p is the number of computational macroparticles representing the protons. The two bottom panels of figure 4 show data from the end of the linear phase atsnapshot ( iv ). The proton v (cid:107) structure as a function of gyrophase from the same twoconfiguration space bins as in panel (b) is shown for contrast in the panel (d). Thewaves in v (cid:107) are beginning to break nonlinearly.The predominantly drift, as distinct from gyroresonant, character of the instabilitycan be seen in the fully resolved velocity space of the protons. However it is challengingto extract the nature of the interaction between protons and the excited waves, becausesnapshots in time of proton phase space contain the history of all previous interactions.We now present a method for modelling the effect of previous interactions using onlyunperturbed proton trajectories which, of course, are easily calculable.In this linear phase the interaction between the protons and the waves is weak:the protons are not trapped but follow trajectories through phase space close to thosefollowed by gyro-orbits in the unperturbed fields. We can therefore directly verify thatthe large scale phase space structure in figure 4(b) arises from interaction with the elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas ω (cid:39) ω LH , k = − ω pe /c . To do this, we will reconstruct a phase spacesnapshot to compare with figure 4, assuming that: protons move on unperturbed gyro-orbits; and the perturbations in v (cid:107) are proportional to the wave amplitude experiencedby the proton when last at resonance. Thus we can see the history of previous wave-particle interactions in the phase space of protons not via perturbations to the trajectory,but in the history of the amplitude of the wavefield E ( x, t ) at the last point of resonance,which is calculated as follows.The components of parallel and perpendicular velocities of a single proton gyro-orbit projected onto v x can be written, referring back to (2), in the form v x ( t ) = u cos θ − v r ( t ) sin θ (12)Here θ is the angle between the background magnetic field and the x -direction, and v r ( t ) = v r sin( α ( t )), where the gyrophase is α ( t ) = Ω p t + α , and α is an initialcondition. We can calculate the time of resonance t R for this proton by noting that atresonance ω/k = v x , so that v x ( t R ) = ω/k = u cos θ − v r sin( α ( t R )) sin θ (13)The position x ( t ) of the proton is given by x ( t ) = tu cos θ + v r cos( α ( t )) sin θ/ Ω p + x (14)where x is the initial position of the proton. We eliminate the dependence on the initialcondition x by calculating the distance between the position of resonance x R and theposition at time t , x ( t ) − x R = ( t − t R ) u cos θ + v r sin θ (cos( α ( t )) − cos( α ( t R ))) / Ω p (15)The normalized amplitude of the wave at the point of resonance is E ( x, t ) = sin( kx R ( α , x, t ) − ωt R ( α )) (16)We now use (15) and (16) to generate a phase space diagram for protons distributedacross a continuum of initial gyrophase. This is shown in Figure 5, which plots thenormalized wave field E ( x R , t R ) as a function of proton gyrophase α ( t ) at points ( t , x )and ( t , x + λ/ ω (cid:39) ω LH , k = − ω pe /c . Thismethod recreates the large scale structure in Figure 4(b). The fine scale structure canbe attributed in part to the interaction of the protons with the relatively weak forwardtravelling wave, as shown in Figure 5(b) which combines the effects of both waves. Here E ( x, t ) = sin ( k b x Rb − ω b t Rb ) + 0 . sin ( k f x Rf − ω f t Rf ) is plotted at points ( t , x ) and( t , x + λ/ 2) against the gyrophase of the protons, where ( x b , t b ) and ( x f , t f ) are thepoints of resonance for the backward propagating wave ( ω b , k b ) and the weaker forwardpropagating wave ( ω f , k f ) respectively. In this case we take ω f /k f to equal the maximumin proton v x . elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas −2 0 2−101 Gyrophase [rad] F i e l d a m p li t ude [ a r b . u .] −101 x x + λ /2 (a)(b) Figure 5. Panel (a): Normalized wave amplitude seen by protons at resonance withthe dominant backward wave, plotted as a function of gyrophase for two regions inconfiguration space separated by half the wavelength of the dominant wave (see text).Panel (b): As panel (a), for the sum of the backward propagating wavemode and therelatively weak forward wavemode. 6. Conclusions We have studied a collective instability arising from spatially localised energetic ionpopulations that can arise in tokamaks, and the phenomenology of the resultant energytransfer to the thermal background plasma, notably the electrons. We employ a self-consistent relativistic electromagnetic 1D3V particle-in-cell code to evolve minorityenergetic protons and background thermal electrons and deuterons. The 1D3V PICcode provides 3D velocity and electromagnetic field vectors along one spatial directionand time. The single spatial direction restricts the propagation direction of waves toan oblique angle to the magnetic field. Plasma parameters are chosen to approximatewhere possible those at the outer mid-plane edge inside the last closed flux surface in alarge tokamak plasma. The velocity distribution of confined fusion products is modelledby an antiparallel travelling ring beam, motivated by observations and interpretation ofion cyclotron emission in JET and TFTR.In using a PIC code, the initial perturbation is seeded by thermal noise thathas broad spectral content, whereas analytical studies rely on eikonal approaches.Computation of the proton fluctuation energy is valuable in demonstrating the existenceof structured early time perturbations to the ensemble of particles, which are otherwisemasked by noise. The initial perturbation to the proton distribution function generatesa collective instability, resulting in growth of waves in which the electric field issubstantially larger than the magnetic: predominantly electrostatic waves in the lowerhybrid range of frequencies. Snapshots of the proton v (cid:107) distribution function enable usto link the multiharmonic mode structure of the electric field to oscillations in proton v (cid:107) as a function of gyrophase.Fast Fourier transforms of the electric field along the simulation domain before thelinear phase, and during its early and late stages, capture proton-excited ω, k structure. elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas 7. Acknowledgements This work was partly funded by the UK Engineering and Physics Sciences ResearchCouncil under grant EP/G003955 and by The European Communities under thecontract of association between Euratom and CCFE. The views and opinions expressedherein do not necessarily represent those of the European Communities. The authorswould like to thanks Christopher Brady and the epoch development team for their helpwith the code. 8. References [1] R. O. Dendy, K. G. McClements, C. N. Lashmore-Davies, G. A. Cottrell, R. Majeski, andS. Cauffman. Ion cyclotron emission due to collective instability of fusion products and beamions in TFTR and JET. Nuclear Fusion , 35:1733–1742, Dec 1995.[2] Nathaniel J. Fisch and Jean-Marcel Rax. Interaction of energetic alpha particles with intenselower hybrid waves. Phys. Rev. Lett. , 69(4):612–615, Jul 1992.[3] T. H. Stix. Waves in Plasmas , page 29. Springer-Verlag New York, Inc., New York, 1992.[4] G. A. Cottrell and R. O. Dendy. Superthermal radiation from fusion products in JET. Phys. Rev.Lett. , 60(1):33–36, Jan 1988.[5] S. Cauffman, R. Majeski, K.G. McClements, and R.O. Dendy. Alfv´enic behaviour of alpha particledriven ion cyclotron emission in TFTR. Nuclear Fusion , 35(12):1597–602, 1995.[6] P.Schild, G.A.Cottrell, and R.O.Dendy. Sawtooth oscillations in ion-cyclotron emission from JET. Nuclear Fusion , 29(5):834–839, 1989.[7] G.A. Cottrell, V.P. Bhatnagar, O. Da Costa, R.O. Dendy, J. Jacquinot, K.G. McClements, D.C.McCune, M.F.F. Nave, P. Smeulders, and D.F.H. Start. Ion cyclotron emission measurementsduring JET deuterium-tritium experiments. Nuclear Fusion , 33(9):1365–1387, 1993.[8] K. G. McClements, C. Hunt, R. O. Dendy, and G. A. Cottrell. Ion cyclotron emission from JETD-T plasmas. Phys. Rev. Lett. , 82(10):2099–2102, Mar 1999.[9] K. G. McClements, R. O. Dendy, C. N. Lashmore-Davies, G. A. Cottrell, S. Cauffman, and elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas R. Majeski. Interpretation of ion cyclotron emission from sub-Alfv´enic fusion products in theTokamak Fusion Test Reactor. Physics of Plasmas , 3(2):543–553, 1996.[10] R. O. Dendy, K. G. McClements, C. N. Lashmore-Davies, R. Majeski, and S. Cauffman. Amechanism for beam-driven excitation of ion cyclotron harmonic waves in the Tokamak FusionTest Reactor. Physics of Plasmas , 1(10):3407–3413, 1994.[11] R O Dendy. Interpretation of ion cyclotron emission from fusion and space plasmas. PlasmaPhysics and Controlled Fusion , 36(12B):B163–B172, 1994.[12] R.O. Dendy, C.N. Lashmore-Davies, and K.F. Kam. A possible excitation mechanism for observedsuperthermal ion cyclotron emission from tokamak plasmas. Physics of Fluids B: PlasmaPhysics , 4(12):3996–4006, 1992.[13] R. O. Dendy, C. N. Lashmore-Davies, K. G. McClements, and G. A. Cottrell. The excitation ofobliquely propagating fast Alfv´en waves at fusion ion cyclotron harmonics. Physics of Plasmas ,1(6):1918–1928, 1994.[14] K. G. McClements and R. O. Dendy. Ion-cyclotron harmonic wave generation by ring protons inspace plasmas. Journal of Geophysical Research-Space Physics , 98(A7):11689–11700, July 1993.[15] K. G. McClements, R. O. Dendy, and C.N. Lashmore-Davies. A model for the generation ofobliquely propagating ULF waves near the magnetic equator. Journal of Geophysical Research-Space Physics , 99(A12):23685–23693, Dec 1994.[16] K G McClements, R O Dendy, L O’C Drury, and P Duffy. Excitation of ion cyclotron harmonicwaves in cosmic ray shock precursors. Monthly Notices of the Royal Astronomical Society ,280(219-226), 1996.[17] K. G. McClements, R. O. Dendy, R. Bingham, J. G. Kirk, and L. O. Drury. Acceleration ofcosmic ray electrons by ion-excited waves at quasi-perpendicular shocks. Monthly Notices of theRoyal Astronomical Society , 291:241–249, October 1997.[18] J G Kirk and R O Dendy. Shock acceleration of cosmic rays - a critical review. Journal of PhysicsG: Nuclear and Particle Physics , 27(7):1589, 2001.[19] J W S Cook, S C Chapman, and R O Dendy. Electron current drive by fusion-product-excitedlower hybrid drift instability. arXiv:1009.1041v1 [physics.plasm-ph].[20] Peter H. Yoon and Anthony T. Y. Lui. Effects of magnetized ions on the lower-hybrid-driftinstability. Physics of Plasmas , 10(11):4260–4264, 2003.[21] P. H. Yoon, Y. Lin, X. Y. Wang, and A. T. Y. Lui. Theory and simulation of lower-hybrid driftinstability for current sheet with guide field. Physics of Plasmas , 15(11):112103, 2008.[22] J. B. Hsia, S. M. Chiu, M. F. Hsia, R. L. Chou, and C. S. Wu. Generalized lower-hybrid-driftinstability. Physics of Fluids , 22(9):1737–1746, 1979.[23] R. C. Davidson, N. T. Gladd, C. S. Wu, and J. D. Huba. Effects of finite plasma beta on thelower-hybrid-drift instability. Physics of Fluids , 20(2):301–310, 1977.[24] Zoltan Dobe, Kevin B. Quest, Vitali D. Shapiro, Karoly Szego, and Joseph D. Huba. Interactionof the solar wind with unmagnetized planets. Phys. Rev. Lett. , 83(2):260–263, Jul 1999.[25] Peter H. Yoon, Anthony T. Y. Lui, and Chia-Lie Chang. Lower-hybrid-drift instability operativein the geomagnetic tail. Physics of Plasmas , 1(9):3033–3043, 1994.[26] J. F. Drake, P. N. Guzdar, A. B. Hassam, and J. D. Huba. Nonlinear mode coupling theory ofthe lower-hybrid-drift instability. Physics of Fluids , 27(5):1148–1159, May 1983.[27] N T Gladd. The lower hybrid drift instability and the modified two stream instability in highdensity theta pinch environments. Plasma Physics , 18(1):27–40, 1976.[28] T. A. Carter, H. Ji, F. Trintchouk, M. Yamada, and R. M. Kulsrud. Measurement of lower-hybriddrift turbulence in a reconnecting current sheet. Phys. Rev. Lett. , 88(1):015001, Dec 2001.[29] S. D. Bale, F. S. Mozer, and T. Phan. Observation of lower hybrid drift instability in the diffusionregion at a reconnecting magnetopause. Geophysical Research Letters , 29(24):33–36, 2002.[30] J. D. Huba and K. Papadopoulos. Nonlinear stabilization of the lower-hybrid-drift instability byelectron resonance broadening. Physics of Fluids , 21(1):121–123, 1978.[31] William Daughton, Giovanni Lapenta, and Paolo Ricci. Nonlinear evolution of the lower-hybrid elf-consistent kinetic simulations of lower hybrid drift instability resulting in electron current driven by fusion products in tokamak plasmas drift instability in a current sheet. Phys. Rev. Lett. , 93(10):105004, Sep 2004.[32] O. J. G. Silveira, L. F. Ziebell, R. Gaelzer, and Peter H. Yoon. Unified formulation forinhomogeneity-driven instabilities in the lower-hybrid range. Phys. Rev. E , 65(3):036407, Feb2002.[33] A K Gwal, M S Tiwari, and K D Misra. Lower hybrid drift instability in the ionosphere. PhysicaScripta , 19(5-6):533–537, 1979.[34] Yu-Jiuan Chen and C. K. Birdsall. Lower-hybrid drift instability saturation mechanisms in one-dimensional simulations. Physics of Fluids , 26(1):180–189, 1983.[35] M E Dieckmann, K G McClements, S C Chapman, R O Dendy, and L O’C Drury. Electronacceleration due to high frequency instabilities at supernova remnant shocks. Astronomy andAstrophysics , 356(377), 2000.[36] K. G. McClements, M. E. Dieckmann, A. Ynnerman, S. C. Chapman, and R. O. Dendy.Surfatron and stochastic acceleration of electrons at supernova remnant shocks. Phys. Rev.Lett. , 87(25):255002, Nov 2001.[37] H.Schmitz, S.C.Chapman, and R.O. Dendy. The influence of electron temperature and magneticfield strength on cosmic-ray injection in high Mach number shocks. Astrophysical Journal ,570:637–646, May 2002.[38] H. Schmitz, S. C. Chapman, and R. O. Dendy. Electron preacceleration mechanisms in the footregion of high Alfv´enic mach number shocks. The Astrophysical Journal , 579(1):327–336, 2002.[39] Jae Koo Lee and C. K. Birdsall. Velocity space ring-plasma instability, magnetized, Part II:Simulation. Physics of Fluids , 22(7):1315–1322, 1979.[40] William Daughton. Nonlinear dynamics of thin current sheets. Physics of Plasmas , 9(9):3668–3678, 2002.[41] Tsunehiko N. Kato. Saturation mechanism of the Weibel instability in weakly magnetized plasmas. Physics of Plasmas , 12(8):080705, 2005.[42] K. Theilhaber and C. K. Birdsall. Kelvin-Helmholtz vortex formation and particle transport in across-field plasma sheath. Phys. Rev. Lett. , 62(7):772–775, Feb 1989.[43] R. E. Lee, S. C. Chapman, and R. O. Dendy. Numerical Simulations of Local Shock Reformationand Ion Acceleration in Supernova Remnants. The Astrophysical Journal , 604:187–195, March2004.[44] R. E. Lee, S. C. Chapman, and R. O. Dendy. Ion acceleration processes at reforming collisionlessshocks. Physics of Plasmas , 12(1):012901, 2005.[45] P. E. Devine, S. C. Chapman, and J. W. Eastwood. One- and two-dimensional simulations ofwhistler mode waves in an anisotropic plasma. Journal of Geophysical Research-Space Physics ,100:17189–17204, September 1995.[46] M. Oppenheim, D. L. Newman, and M. V. Goldman. Evolution of electron phase-space holes ina 2D magnetized plasma. Phys. Rev. Lett. , 83(12):2344–2347, Sep 1999.[47] M. E. Dieckmann, S. C. Chapman, A. Ynnerman, and G. Rowlands. The energy injection intowaves with a zero group velocity. Physics of Plasmas , 6:2681–2692, July 1999.[48] P. C. Birch and S. C. Chapman. Detailed structure and dynamics in particle-in-cell simulationsof the lunar wake. Physics of Plasmas , 8(10):4551–4559, 2001.[49] P. C. Birch and S. C. Chapman. Particle-in-cell simulations of the lunar wake with high phasespace resolution. Geophys. Res. Lett. , 28(2):219, 2001.[50] M. Bonitz, G. Bertsch, V. S. Filinov, and H.Ruhl. Introduction to Computational Methods inMany Body Physics . Cambridge University Press, Cambridge, 2004.[51] R O Dendy, S C Chapman, and M Paczuski. Fusion, space and solar plasmas as complex systems. Plasma Physics and Controlled Fusion , 49(5A):A95, 2007.[52] S.V. Putvinski. Physics of energetic particles in ITER.