Electron current drive by fusion-product-excited lower hybrid drift instability
We present first principles simulations of the direct collisionless coupling of the free energy of fusion-born ions into electron current in a magnetically confined fusion plasma. These simulations demonstrate, for the first time, a key building block of some “alpha channelling” scenarios for tokamak experiments. A fully self-consistent electromagnetic 1D3V particle-in-cell code is used to evolve a parallel drifting ring-beam distribution of 3MeV protons in a 10keV thermal deuterium-electron plasma with realistic mass ratio. Collective instability gives rise to electromagnetic field activity in the lower hybrid range of frequencies. These spontaneously excited obliquely propagating waves undergo Landau damping on resonant electrons, drawing out an asymmetric tail in the distribution of electron parallel velocities, which carries a current.
Optimal exploitation of the free energy of fusion products, for example the alpha-particles born at 3.5MeV in reactions between thermal deuterons and tritons at 10-20keV, is central to achieving fusion power through magnetic confinement of plasma. In the traditional framework, this energy is transferred through multiple collisions to the thermal electrons on a cumulative timescale 1s; electron heating by fusion alpha-particles has been observed in the TFTR1 and JET2 tokamaks. The electrons in turn sustain the temperature of the thermal ions to which they are collisionally coupled. It may be preferable, however, to use some of this fusion product free energy in “alpha channelling” scenarios. This term, coined by Fisch and Rax3 , refers to mechanisms by which rapid collisionless collective instabilities (natural or induced) of the fusion product population could directly benefit the plasma equilibrium, for example by helping sustain toroidal current4 ; 5 ; 6 . Here we report particle-in-cell (PIC) simulations of fusion-born protons in deuterium plasmas that demonstrate from first principles, for the first time, a key alpha channelling phenomenon for tokamak fusion plasmas. We focus on the collective instability of centrally born fusion products on banana orbits at the plasma edge, a population known to be responsible for observations of ion cyclotron emission in JET7 and TFTR8 . A fully self-consistent electromagnetic 1D3V PIC code evolves a parallel drifting ring-beam distribution of 3MeV protons in a 10keV thermal deuterium-electron plasma with realistic mass ratio. Collective instability gives rise to electromagnetic field activity in the lower hybrid range (LH) of frequencies. The spontaneously excited obliquely propagating waves undergo Landau damping on resonant electrons, drawing out an asymmetric tail in the distribution of electron parallel velocities, which carries a current. These simulations thus demonstrate a key building block of some alpha channelling scenarios: the direct collisionless coupling of fusion product free energy into a form which can help sustain the basic equilibrium of the tokamak plasma.
Spatially localized inversions of the velocity distribution of fusion-born ions can arise due to the particle energy and pitch angle dependence of particle orbits. Alpha-particles born with pitch angles just inside the trapped-passing boundary generate ion cyclotron emission7 ; 8 , for example. This motivates our model9 of the initial fusion product velocity distribution function as a ring travelling anti-parallel to the magnetic field with distribution function ; is the magnetic field aligned velocity and is the perpendicular velocity. The radial extent of the emitting region is of order a few energetic particle gyroradii in JET and TFTR7 ; 8 .
Our simulations use physical mass ratios to ensure that the physically relevant instability is excited. We select parameter values similar to those in the edge of a large tokamak, subject to computational resource constraints. The minority energetic fusion product protons at 3MeV have pitch angle and their number density is of that of the background (Maxwellian) 10keV deuterons . The electron number density of the initially quasineutral plasma is and the electron beta . The speed of the energetic protons is approximately half the Alfvén speed given the applied magnetic field . These parameter values imply a total energy of the energetic protons times that of the thermal deuterons and electrons combined. While this is times the value anticipated for next step fusion plasmas10 , it is necessary so as to drive the instability on an acceptable timescale computationally.
The collective instability of the energetic protons can be identified as a form of lower hybrid drift instability (LHDI), a topic of widespread relevance to space and astrophysical plasmas11 ; 12 , while LH waves are used for current drive13 ; 14 ; 15 in the fusion context and are associated with acceleration mechanisms in ionospheric plasmas16 . The treatment here extends LHDI theory in several respects: the energetic ion population does not contribute to the plasma equilibrium, unlike most space and astrophysical applications, where the ion beams are typically associated with currents and gradients central to the equilibrium; values of the key dimensionless parameters are guided by large tokamak edge plasma conditions; and the physically correct mass ratio of electrons to ions is used, which was not possible in some of the classic computational studies of LHDI.
Our fully self-consistent, relativistic kinetic simulations use a particle-in-cell code epoch1d based on the approach of Ref.17 . Computational macroparticles represent the particle distribution functions in full three dimensional velocity space. All three components of all vector quantities, that is, particle velocities and electromagnetic fields, are functions of a single spatial coordinate (, referred to as the direction of variation or the simulation domain) and time . Field and particle boundary conditions are periodic. Wavevectors are then parallel or antiparallel to the direction. To focus on obliquely propagating modes, we set the background field at an angle to the direction. The fields and bulk plasma parameters are resolved (in ) by computational cells of width , where is the electron Debye length. The simulation domain (in ) is optimized to capture the rather broad range of characteristic lengthscales of the energetic protons, background deuterons and electrons. Each of the cells corresponds to species gyroradii of , and . The simulation requires over 2 million computational particles to give good phase space resolution and is typically iterated over timesteps.
An overview of the energy flows in the simulation is shown in Figure 1. The total electric and magnetic energy in excited fields in the simulation are obtained by summing the (suitably normalized) energy densities over all grid cells giving electric and magnetic field energies where is the applied magnetic field. Both the electric field energy (light blue trace) and the magnetic field energy (magenta trace) rise with time. Four stages of the evolving simulation are indicated by vertical lines (i)-(iv) corresponding to times and lower hybrid oscillation periods respectively. Stage (i) corresponds to the onset of the linear phase of field growth, which is well developed by stages (ii) and (iii). By stage (iv) the wave amplitude is approaching its saturated level. The magnetic fluctuations contain nearly an order of magnitude less energy than the electric fluctuations, implying that the waves produced by the instability are largely electrostatic.
The total kinetic energy of each particle species is obtained by summing over all the computational particles in the system. From stage (iii) onwards, Figure 1 shows that the total kinetic energies of the protons and electrons approximately mirror each other, the protons losing energy as the electrons gain energy, whilst the deuterons show little change.
To obtain the relative gain, or loss, in the energy of each population we define a change in total kinetic energy of the ensemble of particles , and a total kinetic energy of fluctuations relative to the initial conditions, where indicates an average over all the particles of species . The time variation of the proton fluctuation energy is shown in Figure 1. It grows from the start of the simulation, ultimately increasing by three orders of magnitude, whereas the total proton kinetic energy declines by much less than one order of magnitude. This reflects the role of proton energy dispersion in the early phase of the instability. The electron kinetic energy rises with the electric field energy during the linear stage of the simulation, arising from electron participation in the principally electrostatic waves excited by the instability. The corresponding effect for the deuterons is also visible, but is much less due to their higher mass.
We now turn to the evolving distribution function of the electron parallel momentum shown in Figure 2 for times (i)-(iv). From stage (iii) onwards the electron distribution function develops an asymmetric tail in reflecting net directional electron acceleration. We infer Landau damping of the excited waves on resonant electrons which results in the flattening of the negative tail of the distribution function.
This is confirmed by the analysis shown in Figure 3 where we plot the spatio-temporal fast Fourier transform of the electric field component along the simulation direction , transformed from for the time interval ending at time (iv). The oblique cold-plasma normal mode of the background deuteron-electron plasma is marked by the black dashed trace. The presence of the additional energetic proton population modifies the nature and dispersion properties of the normal modes of the plasma and, through resonance, couple energy to these modified modes. The Figure shows peaks in intensity at which indicate resonance of the proton population with a modified deuteron-electron normal mode branch that lies on the surface between the lower extraordinary wave and the whistler wave in space. The locations in space where coupling takes place correspond to where the velocities of the peaks of the proton distribution function in approximately match the phase velocities of the normal mode. As the proton distribution function evolves, the positions in of the peaks in the distribution function move. These velocities at and at time (iv) are plotted on Figure 3 as black dash-dot traces. We note that there are harmonics of the dominant backward propagating (in negative direction) mode at .
On Figure 2 we indicate with vertical lines the two electron parallel momenta resonant with the dominant backward and forward electromagnetic structures shown in Figure 3. Taken together, Figures 2 and 3 indicate that the electrons are principally accelerated by the dominant wave at negative phase velocity, which is in turn excited by the backward-drifting fusion products.
To capture the coherent oscillatory features of the dominant excited fields, and hence their phase resonant characteristics, we decompose the electric field into component waves travelling with positive and negative phase velocities. These are and , where IFT denotes the inverse Fourier transform, which we plot in Figure 4. The sum of and recreates the original field. In Figure 4 , the amplitude of the backward travelling wave is approximately an order of magnitude greater than the forward travelling wave and is dominated by a single long wavelength component. The parallel phase velocity of this backward travelling wave is shown by the vertical line on the left of Figure 2.
The predominantly drift character of the instability is directly seen in the fully resolved velocity space of the PIC simulation. The velocity of each particle can be expressed in terms of a field aligned component , a component aligned with the simulation domain (and ) , and gyrophase . We test for velocity space patterns in these coordinates at snapshot (iii) when the wavefields are well established. We select protons from a narrow region in configuration space that is smaller in extent that the wavelength of the dominant wavemode (). The coordinates , , and the gyrophase for each particle are plotted in panel (a) of Figure 5, where gyrophase is along the abscissa, is along the ordinate and colour indicates . We see asymmetry in the oscillatory pattern in which is a function of : a characteristic of drift instability rather than gyroresonance. As such, these oscillatory perturbations in should vary as the wave amplitude at the point of resonance , and in panel (b) of Figure 5 we confirm that there is indeed resonance with the dominant LH wave identified above. Panel 5(b) plots on the ordinate the normalized wave amplitude , that is, the sine of the wave phase at the location and time of unperturbed test particles that satisfy the condition for resonance with the dominant wave . The point of resonance is obtained for test particles which initially are distributed uniformly in gyrophase and which follow unperturbed cycloidal orbits. To reconstruct a single snapshot in and (to compare with the simulation snapshot of panel (a)) these test particle orbits are then advanced in , and to a single point . We plot in panel (b), for a snapshot in space and time , the normalized wave amplitude experienced by the particle at the point of resonance against particle gyrophase with , represented in colour. We see that panel 5(b) closely tracks the large scale oscillations seen in the PIC code velocity space of panel 5(a). We refine this idea in panel 5(c) where we plot the effective wave amplitude arising from both the forward and backward waves identified in Figure 4; the weaker second wave can be seen to introduce fine scale oscillations visible in panel 5(a).
These first principles simulations demonstrate, for the first time, key physical elements of alpha channelling scenarios for future tokamak plasmas: the excitation of LH waves by a fusion product population, whose functional form and parameters are aligned with prior observations on JET and TFTR, combined with the subsequent Landau damping of the excited waves on resonant electrons, drawing out an asymmetric tail in their parallel velocity distribution, which carries a current. These simulations also contribute to the question of what instabilities may arise –- transiently, locally, or otherwise –- during the initiation and propagation of fusion burn in tokamak plasmas. Furthermore they deepen our understanding of LH drift instabilities, which may be widespread in space and astrophysical plasmas, by extending their study into parameter ranges that approximate (subject to computational resource constraints) edge plasma conditions in large tokamaks.
Acknowledgements.The authors acknowledge the EPSRC for support. This work was supported in part by EURATOM under a contract of association with the CCFE. The views expressed do not necessarily reflect those of the European Commission.
- (1) K. M. McGuire et al, Phys Plasmas 2, 2176 (1995); G. Taylor, J. D. Strachan, R. V. Budny and D. R. Ernst, Phys. Rev. Lett. 76, 2722 (1996)
- (2) P. R. Thomas et al, Phys. Rev. Lett. 80, 5548 (1998)
- (3) N. J. Fisch, and J-M. Rax, Phys. Rev. Lett. 69, 612 (1992); N. J. Fisch, Nucl Fusion 40, 1095 (2000)
- (4) M. C. Herrmann and N. J. Fisch, Phys. Rev. Lett. 79, 1495 (1997)
- (5) A. J. Fetterman and N.J Fisch, Phys. Rev. Lett. 101, 205003 (2008)
- (6) D. A. Gates, N. N. Gorelenkov and R. B. White, Phys. Rev. Lett. 87, 205003 (2001)
- (7) G. A. Cottrell and R. O. Dendy, Phys. Rev. Lett. 60, 33 (1988); G A Cottrell et al, Nucl Fusion 33, 1365 (1993); K. G. McClements, C. Hunt, R. O. Dendy and G. A. Cottrell, Phys. Rev. Lett. 82, 2099 (1999)
- (8) S. Cauffman et al, Nucl. Fusion 35, 1597 (1995)
- (9) R. O. Dendy, C. N. Lashmore-Davies, K. G. McClements and G. A. Cottrell, Phys. Plasmas 1, 1918 (1994)
- (10) S. V. Putvinski, Nucl. Fusion 38, 1275 (1998)
- (11) W. Daughton, G. Lapenta and P Ricci, Phys. Rev. Lett. 93, 105004 (2004)
- (12) T. A. Carter et al, Phys. Rev. Lett. 88, 015001 (2001)
- (13) S. Bernabei et al, Phys. Rev. Lett. 49, 1255 (1982)
- (14) N. C. Hawkes et al, Phys. Rev. Lett. 87, 115001 (2001)
- (15) G. T. Hoang et al, Phys. Rev. Lett. 90, 155002 (2003)
- (16) P. M. Kintner et al, Phys. Rev. Lett. 68, 2448 (1992)
- (17) M. Bonitz, G. Bertsch, V. S. Filinov and H. Ruhl, Introduction to Computational Methods in Many Body Physics, Cambridge University Press: Cambridge (2004)