Guiding Center Simulations Research Articles

Kinetic particle simulations in a global toroidal geometry

The gyrokinetic toroidal code has been upgraded for global simulations by coupling the core and scrape-off layer regions across the separatrix with field-aligned particle-grid interpolations. A fully kinetic particle pusher for high frequency waves (ion cyclotron frequency and beyond) and a guiding center pusher for low frequency waves have been implemented using cylindrical coordinates in a global toroidal geometry. The two integrators correctly capture the particle orbits and agree well with each other, conserving energy and canonical angular momentum. As a verification and application of this new capability, ion guiding center simulations have been carried out to study ion orbit losses at the edge of the DIII-D tokamak for single null magnetic separatrix discharges. The ion loss conditions are examined as a function of the pitch angle for cases without and with a radial electric field. The simulations show good agreement with past theoretical results and with the experimentally observed feature in which high energy ions flow out along the ion drift orbits and then hit the divertor plates. A measure of the ion direct orbit loss fraction shows that the loss fraction increases with the ion energy for DIII-D in the initial velocity space. Finally, as a further verification of the capability of the new code, self-consistent simulations of zonal flows in the core region of the DIII-D tokamak were carried out. All DIII-D simulations were performed in the absence of turbulence.

Collisional enhancement of energetic particle Alfvénic resonance width in tokamaks

The phase-space structure of resonances between fast ions and an Alfvénic mode and the associated modification of density profiles in tokamaks are studied as a function of particle collisions. Guiding-center simulations in a realistic tokamak equilibrium are employed to address the resonance broadening parametric dependencies with respect to changes in the pitch-angle scattering rate. The rate of collisional replenishment, along with resonance strength, given by the combination of eigenmode and resonance structures and equilibrium parameters, determines saturation amplitudes for a given damping rate. As seen from the distribution function flattening, collisions have an effect of broadening the resonances, while the absolute value of δf decreases with increasing collisionality. It is observed that the collisional broadening can be comparable to the collisionless resonance width due to the mode amplitude alone. The resonance broadening coefficients are compared with the existing theory based on analytically expected saturation levels, showing fair agreement. The results can be useful in assisting reduced kinetic models, such as quasilinear models, when prescribing the effective resonance phase-space width, i.e., the mode-particle interaction platform, due to collisional or turbulent processes.

Guiding center simulations on curvilinear grids

Semi-Lagrangian guiding center simulations are performed on sinusoidal perturbations of cartesian grids, thanks to the use of a B-spline finite element solver for the Poisson equation and the classical backward semi-Lagrangian method (BSL) for the advection. We are able to reproduce the standard Kelvin-Helmholtz instability test on such grids. When the perturbation leads to a strong distorted mesh, we observe that the solution differs if one takes standard numerical parameters that are used in the cartesian reference case. We can recover good results together with correct mass conservation, by diminishing the time step.

Fast-ion losses induced by ELMs and externally applied magnetic perturbations in the ASDEX Upgrade tokamak

Phase-space time-resolved measurements of fast-ion losses induced by edge localized modes (ELMs) and ELM mitigation coils have been obtained in the ASDEX Upgrade tokamak by means of multiple fast-ion loss detectors (FILDs). Filament-like bursts of fast-ion losses are measured during ELMs by several FILDs at different toroidal and poloidal positions. Externally applied magnetic perturbations (MPs) have little effect on plasma profiles, including fast-ions, in high collisionality plasmas with mitigated ELMs. A strong impact on plasma density, rotation and fast-ions is observed, however, in low density/collisionality and q95 plasmas with externally applied MPs. During the mitigation/suppression of type-I ELMs by externally applied MPs, the large fast-ion bursts observed during ELMs are replaced by a steady loss of fast-ions with a broad-band frequency and an amplitude of up to an order of magnitude higher than the neutral beam injection (NBI) prompt loss signal without MPs. Multiple FILD measurements at different positions, indicate that the fast-ion losses due to static 3D fields are localized on certain parts of the first wall rather than being toroidally/poloidally homogeneously distributed. Measured fast-ion losses show a broad energy and pitch-angle range and are typically on banana orbits that explore the entire pedestal/scrape-off-layer (SOL). Infra-red measurements are used to estimate the heat load associated with the MP-induced fast-ion losses. The heat load on the FILD detector head and surrounding wall can be up to six times higher with MPs than without 3D fields. When 3D fields are applied and density pump-out is observed, an enhancement of the fast-ion content in the plasma is typically measured by fast-ion D-alpha (FIDA) spectroscopy. The lower density during the MP phase also leads to a deeper beam deposition with an inward radial displacement of ≈2 cm in the maximum of the beam emission. Orbit simulations are used to test different models for 3D field equilibrium reconstruction including vacuum representation, the free boundary NEMEC code and the two-fluid M3D-C1 code which account for the plasma response. Guiding center simulations predict the maximum level of losses, ≈2.6%, with NEMEC 3D equilibrium. Full orbit simulations overestimate the level of losses in 3D vacuum fields with ≈15% of lost NBI ions.

Guiding center equations of high accuracy

Guiding center simulations are an important means of predicting the effect of resistive and ideal magnetohydrodynamic instabilities on particle distributions in toroidal magnetically confined thermonuclear fusion research devices. Because saturated instabilities typically have amplitudes of δB/B of a few times 10−4 numerical accuracy is of concern in discovering the effect of mode particle resonances. We develop a means of following guiding center orbits which is greatly superior to the methods currently in use. For full implementation, the method requires some breaking of axisymmetry, either through toroidal field ripple or magnetohydrodynamic instabilities. In the presence of ripple or time dependent magnetic perturbations both energy and canonical momentum are conserved in a time step to better than one part in 1014, an improvement of nine orders of magnitude over standard Runge–Kutta integration, and the relation between changes in canonical momentum and energy is also conserved to very high order.

Guiding center equations for ideal magnetohydrodynamic modes

Guiding center simulations are routinely used for the discovery of mode-particle resonances in tokamaks, for both resistive and ideal instabilities and to find modifications of particle distributions caused by a given spectrum of modes, including large scale avalanches during events with a number of large amplitude modes. One of the most fundamental properties of ideal magnetohydrodynamics is the condition that plasma motion cannot change magnetic topology. The conventional representation of ideal magnetohydrodynamic modes by perturbing a toroidal equilibrium field through δB→=∇×(ξ→×B→), however, perturbs the magnetic topology, introducing extraneous magnetic islands in the field. A proper treatment of an ideal perturbation involves a full Lagrangian displacement of the field due to the perturbation and conserves magnetic topology as it should. In order to examine the effect of ideal magnetohydrodynamic modes on particle trajectories, the guiding center equations should include a correct Lagrangian treatment. Guiding center equations for an ideal displacement ξ→ are derived which preserve the magnetic topology and are used to examine mode particle resonances in toroidal confinement devices. These simulations are compared to others which are identical in all respects except that they use the linear representation for the field. Unlike the case for the magnetic field, the use of the linear field perturbation in the guiding center equations does not result in extraneous mode particle resonances.

Guiding-center simulations on curvilinear meshes

The purpose of this work is to design simulation tools for magnetised plasmas in the ITER project framework. The specific issue we consider is the simulation of turbulent transport in the core of a Tokamak plasma, for which a 5D gyrokinetic model is generally used, where the fast gyromotion of the particles in the strong magnetic field is averaged in order to remove the associated fast time-scale and to reduce the dimension of 6D phase space involved in the full Vlasov model. Very accurate schemes and efficient parallel algorithms are required to cope with these still very costly simulations. The presence of a strong magnetic field constrains the time scales of the particle motion along and accross the magnetic field line, the latter being at least an order of magnitude slower. This also has an impact on the spatial variations of the observables. Therefore, the efficiency of the algorithm can be improved considerably by aligning the mesh with the magnetic field lines. For this reason, we study the behavior of semi-Lagrangian solvers in curvilinear coordinates. Before tackling the full gyrokinetic model in a future work, we consider here the reduced 2D Guiding-Center model. We introduce our numerical algorithm and provide some numerical results showing its good properties.

Modification of particle distributions by MHD instabilities I

The modification of particle distributions by low amplitude magnetohydrodynamic modes is an important topic for magnetically confined plasmas. Low amplitude modes are known to be capable of producing significant modification of injected neutral beam profiles, and the same can be expected in burning plasmas for the alpha particle distributions. Flattening of a distribution due to phase mixing in an island or due to portions of phase space becoming stochastic is a process extremely rapid on the time scale of an experiment but still very long compared to the time scale of guiding center simulations. Thus it is very valuable to be able to locate significant resonances and to predict the final particle distribution produced by a given spectrum of magnetohydrodynamic modes. In this paper we introduce a new method of determining domains of phase space in which good surfaces do not exist and use this method for quickly determining the final state of the particle distribution without carrying out the full time evolution leading to it.

Review of modeling of losses and sources of relativistic electrons in the outer radiation belt I: Radial transport

In this paper, we focus on the modeling of radial transport in the Earth's outer radiation belt. A historical overview of the first observations of the radiation belts is presented, followed by a brief description of radial diffusion. We describe how resonant interactions with poloidal and toroidal components of the ULF waves can change the electron's energy and provide radial displacements. We also present radial diffusion and guiding center simulations that show the importance of radial transport in redistributing relativistic electron fluxes and also in accelerating and decelerating radiation belt electrons. We conclude by presenting guiding center simulations of the coupled particle tracing and magnetohydrodynamic (MHD) codes and by discussing the origin of relativistic electrons at geosynchronous orbit. Local acceleration and losses and 3D simulations of the dynamics of the radiation belt fluxes are discussed in the companion paper [Shprits, Y.Y., Subbotin, D.A., Meredith, N.P., Elkington, S.R., 2008. Review of modeling of losses and sources of relativistic electrons in the outer radiation belt II: Local acceleration and loss. Journal of Atmospheric and Solar-Terrestrial Physics, this issue. doi:10.1016/j.jastp.2008.06.014].

Guiding-center simulations of nonlocal and negative inertia effects on rotation in a tokamak

The magnitude of the radial electric field (Er), resulting from nonambipolar fluxes of neoclassical origin, is evaluated using Monte Carlo guiding-center simulations for a low-current plasma corresponding to the FT-2 tokamak {Fisichiskii Tokamak-2, Ioffe Institute, St. Petersburg [S. I. Lashkul, V. N. Budnikov, E. O. Vekshina et al., Plasma Phys. Rep. 27, 1001 (2001)]}. The Er-values are found to significantly exceed those given by the standard neoclassical theory, based on thin-orbit assumption, when the plasma current is sufficiently low and the pressure gradient is sufficiently high. Strong Er-structures are found to form in the same low plasma-current range where enhanced confinement is reported in the FT-2 tokamak. In the simulations, the physics behind the strong increase in the field values is intimately related to the poloidal Mach-number, together with the wide ion orbits.

Asymmetric self-phase modulation and compression of short laser pulses in plasma channels.

A relativistically intense femtosecond laser pulse propagating in a plasma channel undergoes dramatic photon deceleration while propagating a distance on the order of a dephasing length. The deceleration of photons is localized to the back of the pulse and is accompanied by compression and explosive growth of the ponderomotive potential. Fully explicit particle-in-cell simulations are applied to the problem and are compared with ponderomotive guiding center simulations. A numerical Wigner transform is used to examine local frequency shifts within the pulse and to suggest an experimental diagnostic of plasma waves inside a capillary.

Bifurcation of the radial electric field in the presence of edge polarization in tokamaks

In guiding-center simulations of a tokamak plasma edge, bifurcations of both the radial electric field and the electrode current are observed for a threshold value of the electrode bias voltage, in qualitative agreement with experiments. The simulations are neoclassical and Monte Carlo based, also containing the perpendicular viscosity effect. Above the voltage threshold, an inward moving, solitary radial electric field structure is generated with a concomitant drop in the electrode current. The position, width, and height of the soliton depend on the magnitude and radial variation of the conductivity (and, thus, on local plasma parameters), finite orbit effects, and also on the perpendicular viscosity and neutral density. The obtained Er structures are narrower and higher than the ones measured in the edge polarization experiments in the Jülich tokamak TEXTOR-94 [S. Jachmich et al., Plasma Phys. Controlled Fusion 40, 1105 (1998)].

Simulation of the outer radiation belt electron flux decrease during the March 26, 1995, magnetic storm

In this paper we study the variation of the relativistic electron fluxes in the Earth's outer radiation belt during the March 26, 1995, magnetic storm. Using observations by the radiation environment monitor (REM) on board the space technology research vehicle (STRV‐Ib), we discuss the flux decrease and possible loss of relativistic electrons during the storm main phase. In order to explain the observations we have performed fully adiabatic and guiding center simulations for relativistic equatorial electrons in the nonstationary Tsygarienko96 magnetospheric magnetic field model. In our simulations the drift of electrons through the magnetopause was considered as a loss process. We present our model results and discuss their dependence on the magnetospheric magnetic and electric field model, as well as on the prestorm fluxes used in the simulations.

Energetic particle orbits in the National Spherical Tokamak Experiment

The trajectories of neutral beam injected energetic ions in spherical tokamaks are examined. The large poloidal magnetic field in the outboard region of spherical tokamaks causes neutral beam injected ions to be born on trapped orbits even with cotangential injection. Numerical solutions to the equations for particle motion and for guiding center drifts are compared in several magnetic equilibria for a range of particle initial conditions. Even when rL/a∼1/4 the guiding center orbits closely resemble the path of the instantaneous center of gyration of the particle motion; exceptions occur primarily for orbits near the trapped/passing boundary. Finite Larmor radius effects are included in guiding center simulations of prompt orbit loss in the National Spherical Tokamak Experiment (NSTX) [J. Spitzer et al., Fusion Technol. 30, 1337 (1996)]. Orbit loss in the NSTX is caused primarily by collisions with the close fitting conducting shell, and severe losses would be expected with counter directed injection. While most orbits are similar to those found in conventional tokamaks, additional orbit types are possible in spherical tokamaks.

Observations of iron, silicon, and other heavy ions in the geostationary altitude region during late March 1991

Following the great sudden storm commencement of March 24, 1991, observations were made aboard CRRES of Fe, Mg, Si, and other heavy ions in the energy range of several tens to a few hundred keV/n. The charge states of these heavy ions were observed to change abruptly. For example, the dominant Fe charge state changed at one point from q = 9+ to q = 16+. We suggest that ions from the solar corona, retaining the charge state frozen in near the Sim, were convected inward to the location of CRRES from the plasma sheet within minutes. The ions were energized in the convection process, and the ionic charge state was unaltered by charge exchange. The very large magnetospheric electric fields associated with the highly disturbed geomagnetic conditions were required to enable the ions to reach CRRES. Guiding center simulations were carried out which support this hypothesis. At least two ion populations were observed corresponding to temperatures at two different coronal locations, and the observed changes in heavy‐ion charge state were caused by a sequence of solar wind parcels sweeping over the Earth's magnetosphere. The signature of a change from low to high temperature seems to fit well the suggestions of Tsurutani and Gonzales [1994], of a fast solar wind stream tamped by a slower one, as a generator of great storms.

Bounce‐averaged Hamiltonian for charged particles in an axisymmetric but nondipolar model magnetosphere

In order to facilitate bounce‐averaged guiding center simulations of geomagnetically trapped particles, we express the kinetic energy of a particle with magnetic coordinates (L,ϕ) as an analytic function of the first two adiabatic invariants (M, J) and the L value of the field line. The magnetic field model is axisymmetric, consisting of a dipolar B field plus a uniform southward magnetic field parallel to the dipole moment μE. This model magnetosphere is surrounded by a circular equatorial neutral line whose radius b is an adjustable parameter. The L value of a field line is (by definition) inversely proportional to the flux enclosed by the corresponding magnetic shell of equatorial radius r0, and the L value at the neutral line (r0 = b) is denoted L*. The azimuthal coordinate ϕ measures magnetic local time. The best functional representation found for the normalized difference (L³a³/μE)(Bm ‐ B0) between mirror‐point field Bm and equatorial field B0 along any field line is a 5‐term expansion in powers (2/3 through 6/3) of the quantity X ≡ (La/μE)1/2K, where K ≡ (J²/8m0M)1/2 is an adiabatically conserved quantity independent of particle energy, m0 is the rest mass of the particle, and a is the radius of the Earth. This functional form is motivated by results for limiting cases in which particles mirror very near and very far from the magnetic equator. Expansion coefficients corresponding to various powers of X are obtained from least squares fits to numerically computed results for X as a function of L and Bm. These are accurately expressible as fourth‐order polynomials in (r0/b)³, hence indirectly as functions of L/L* = 3La/2b. This representation, which leads (except for a manageably small region of parameter space) to better than 1% accuracy in the specification of Bm as a function of K and L, allows bounce‐averaged guiding center simulations to be performed without actually tracing the bounce motions of individual particles. Bounce‐averaged drifts L (meridional) and ϕ (azimuthal) are proportional to derivatives of the Hamiltonian H (sum of kinetic and potential energies) with respect to ϕ and L, respectively. Our formulation thus provides a computationally efficient method for tracing the bounce‐averaged adiabatic motion (conserving all three invariants) and nonadiabatic transport (violating the third invariant while conserving the first two invariants) of geomagnetically trapped particles in the model magnetosphere.

Simulations of phase space distributions of storm time proton ring current

We use results of guiding‐center simulations of ion transport to map phase space densities of the stormtime proton ring current. We model a storm as a sequence of substorm‐associated enhancements in the convection electric field. Our pre‐storm phase space distribution is an analytical solution to a steady‐state transport model in which quiet‐time radial diffusion balances charge exchange. This pre‐storm phase space spectra at L∼2‐4 reproduce many of the features found in observed quiet‐time spectra. Using results from simulations of ion transport during model storms having main phases of 3, 6, and 12 hr, we map phase space distributions from the pre‐storm distribution in accordance with Liouville's theorem. We find stormtime enhancements in the phase space densities at energies E∼30‐160 keV for L∼2.5‐4. These enhancements agree well with the observed stormtime ring current. For storms with shorter main phases (∼3 hr), the enhancements are caused mainly by the trapping of ions injected from open night side trajectories, and diffusive transport of higher‐energy (≳ 160 keV) ions contributes little to the stormtime ring current. However, the stormtime ring current is augmented also by the diffusive transport of higher‐energy ions (E≳ 160 keV) during storms having longer main phases (≳ 6 hr). In order to account for the increase in Dst associated with the formation of the stormtime ring current, we estimate the enhancement in particle‐energy content that results from stormtime ion transport in the equatorial magnetosphere. We find that transport alone cannot account for the entire increase in |Dst| typical of a major storm. However, we can account for the entire increase in |Dst| by realistically increasing the stormtime outer boundary value of the phase space density relative to the quiet‐time value. We compute the magnetic field produced by the ring current itself and find that radial profiles of the magnetic field depression resemble those obtained from observational data.