Guiding-center Model Research Articles

Testing the conservative character of particle simulations: II. Spurious heating of guiding centers and full orbits subject to fluctuations expressed in terms of E and B

For an axisymmetric tokamak plasma, Hamiltonian theory predicts that the orbits of charged particles must stay on invariant tori of conserved energy in the moving frame of reference of a wave that propagates along the torus with fixed angular phase velocity, amplitude, and shape. The mode structure in the poloidal plane is arbitrary if the fluctuations are expressed in terms of potentials Φ and A, which satisfy Faraday's law and the solenoidal condition by definition. Consequently, smoothing operations (such as gyroaveraging and noise suppression) do not violate the conservative laws. However, this is not guaranteed for models expressed in terms of the physical fields E and B. Here, we demonstrate that manipulations of E and B in the poloidal (R, z) plane can cause spurious heating that is independent of time steps or numerical methods, but can be sensitive to geometry. In particular, we show that secular acceleration is enhanced when one imposes nonnormal modes that possess strong up–down asymmetry instead of the usual in–out asymmetry of normal toroidal (eigen)modes. We compare full gyro-orbit and guiding center models and find similar behavior. We also examine the effect of ad hoc N-point gyroaveraging in a guiding center model, as is done in some simulation codes. If one uses Faraday's law to (re)compute B(t) after gyroaveraging E, the guiding center motion remains conservative. Otherwise, spurious heating should be expected and monitored, but it may be tolerable when normal modes dominate.

Testing the conservative character of particle simulations: I. Canonical and noncanonical guiding center model in Boozer coordinates

The guiding center (GC) Lagrangian in Boozer coordinates for toroidally confined plasmas can be cast into canonical form by eliminating terms containing the covariant component BΨP of the magnetic field vector with respect to the poloidal flux function ΨP. In an unperturbed plasma, BΨP can be eliminated via exact coordinate transformations, but, in general, one relies on approximations, assuming that the effect of BΨP is small. Here, we are interested in the question whether Hamiltonian conservation laws are still satisfied when BΨP is retained in the presence of fluctuations. Considering fast ions in the presence of a shear Alfvén wave field with fixed amplitude, fixed frequency, and a single toroidal mode number n, we show that simulations using the code ORBIT with and without BΨP yield practically the same resonant and nonresonant GC orbits. The numerical results are consistent with theoretical analyses (presented in the appendix), which show that the unabridged GC Lagrangian with BΨP retained yields equations of motion that possess two key properties of Hamiltonian flows: (i) phase space conservation and (ii) energy conservation. As counter-examples, we also show cases where energy conservation (ii) or both conservation laws (i) and (ii) are broken by omitting certain small terms. When testing the conservative character of the simulation code, it is found to be beneficial to apply perturbations that do not resemble normal (eigen) modes of the plasma. The deviations are enhanced and, thus, more easily spotted when one inspects wave-particle interactions using nonnormal modes.

Calculation of collisionless pitch-angle scattering of runaway electrons with synchrotron radiation via high-order guiding-centre equation

Recently, the collisionless pitch-angle scattering for relativistic runaway electrons (REs) in toroidal geometries such as tokamaks was discovered through a full orbit simulation approach (Liu et al., Nucl. Fusion, vol. 56, 2016, p. 064002), and it was then theoretically investigated that a new expression for the magnetic moment, including the second-order corrections, could essentially reproduce the so-called collisionless pitch-angle scattering process (Liu et al., Nucl. Fusion, vol. 58, 2018, p. 106018). In this paper, with synchrotron radiation, extensive numerical verification of the validity of the high-order guiding-centre theory is given for simulations involving REs by incorporating such an expression for the magnetic moment into our particle tracing code. A high-order guiding-centre simulation approach with synchrotron radiation (HGSA) is applied. Synchrotron radiation plays an essential role in the life cycle of REs. The energy of REs first increases and then becomes saturated until the electric field acceleration is balanced by the radiation dissipation. Unfortunately, the process cannot be simulated accurately with the standard guiding-centre model, i.e. the first-order guiding-centre model. Remarkably, it is found that the HGSA can effectively produce the fundamental process of REs. Since the time scale of the energy saturation of REs is close to seconds, the computational cost becomes significant. In order to save costs, it is necessary to estimate the time of energy saturation. An analytical estimate is derived for the time it takes for synchrotron drag to balance an accelerating electric field and the provided formula has been numerically verified. Test calculations reveal that HGSA is favourable for exploiting the dynamics of REs in tokamak plasmas.

Accelerating the estimation of collisionless energetic particle confinement statistics in stellarators using multifidelity Monte Carlo

In the design of stellarators, energetic particle confinement is a critical point of concern which remains challenging to study from a numerical point of view. Standard Monte Carlo (MC) analyses are highly expensive because a large number of particle trajectories need to be integrated over long time scales, and small time steps must be taken to accurately capture the features of the wide variety of trajectories. Even when they are based on guiding center trajectories, as opposed to full-orbit trajectories, these standard MC studies are too expensive to be included in most stellarator optimization codes. We present the first multifidelity Monte Carlo (MFMC) scheme for accelerating the estimation of energetic particle confinement in stellarators. Our approach relies on a two-level hierarchy, in which a guiding center model serves as the high-fidelity model, and a data-driven linear interpolant is leveraged as the low-fidelity surrogate model. We apply MFMC to the study of energetic particle confinement in a four-period quasi-helically symmetric stellarator, assessing various metrics of confinement. Stemming from the very high computational efficiency of our surrogate model as well as its sufficient correlation to the high-fidelity model, we obtain speedups of up to 10 with MFMC compared to standard MC.

A High Order Bound Preserving Finite Difference Linear Scheme for Incompressible Flows

We propose a high order finite difference linear scheme combined with a high order bound preserving maximum-principle-preserving (MPP) flux limiter to solve the incompressible flow system. For such problem with highly oscillatory structure but not strong shocks, our approach seems to be less dissipative and much less costly than a WENO type scheme, and has high resolution due to a Hermite reconstruction. Spurious numerical oscillations can be controlled by the MPP flux limiter. Numerical tests are performed for the Vlasov-Poisson system, the 2D guiding-center model and the incompressible Euler system. The comparison between the linear and WENO type schemes will demonstrate the good performance of our proposed approach.

Analysis of the isotropic and anisotropic Grad–Shafranov equation

Pressure anisotropy is a commonly observed phenomenon in tokamak plasmas, due to external heating methods such as neutral beam injection and ion-cyclotron resonance heating. Equilibrium models for tokamaks are constructed by solving the Grad–Shafranov equation; such models, however, do not account for pressure anisotropy since ideal magnetohydrodynamics assumes a scalar pressure. A modified Grad–Shafranov equation can be derived to include anisotropic pressure and toroidal flow by including drift-kinetic effects from the guiding-centre model of particle motion. In this work, we have studied the mathematical well-posedness of these two problems by showing the existence and uniqueness of solutions to the Grad–Shafranov equation both in the standard isotropic case and when including pressure anisotropy and toroidal flow. A new fixed-point approach is used to show the existence of solutions in the Sobolev space $H_0^1$ to the Grad–Shafranov equation, and sufficient criteria for their uniqueness are derived. The conditions required for the existence of solutions to the modified Grad–Shafranov equation are also constructed.

Inner radiation belt simulations of the proton flux response in the South Atlantic Anomaly during the Geomagnetic Storm of 15 May 2005

A test particle simulation code was developed to simulate the inner proton belt response during the intense geomagnetic storm of 15 May 2005. The guiding center model was implemented to compute the proton trajectories with an energy range of 70–180 MeV. The time-varying magnetic field model implemented in the simulations was computed by the Tsyganenko model TS05 with the associated inductive electric field. One of the most important features of the low-earth orbit (LEO) environment is the South Atlantic Anomaly, which imposes a dangerous radiation load on most LEO missions. This research aims to investigate the proton flux variations in the anomaly region with respect to space weather conditions. The results showed that during the main phase of the geomagnetic storm, the proton flux in the SAA decreased, whereas, throughout the initial and recovery phases, the proton flux was increased at most of the altitudes. Satellite measurements confirmed numerical results.

An efficient trajectory tracking algorithm for the backward semi-Lagrangian method of solving the guiding center problems

In this paper, we develop an effective numerical algorithm that tracks the trajectory needed to solve guiding center models using the backward semi-Lagrangian method. In terms of numerical calculations, two appreciably fast algorithms for the departure points are designed. One is a completely explicit formula for numerical solutions of the discrete system for each Cauchy problem. This formula is characterized by numerical factors that are less than half the multiplication number used in the usual Gaussian elimination. The other finds the required departure points with an interpolation method that is at least 30% less expensive for heavy Cauchy problems. Last, we propose a method to modify the solution to improve estimation of physical quantities, such as conservation of mass, which can be lost in interpolation solutions calculated at the departure points. It turns out that the proposed method not only saves a great deal of computation time, but also preserves physical quantities such as mass and total kinetic energy much better than conventional methods. To demonstrate the numerical evidence, we use the proposed method to simulate several problems such as the incompressible Euler equation, Kelvin-Helmholtz instability, Diocotron instability and a three-dimensional guiding center model.

Solving hyperbolic-elliptic problems on singular mapped disk-like domains with the method of characteristics and spline finite elements

A common strategy in the numerical solution of partial differential equations is to define a uniform discretization of a tensor-product multi-dimensional logical domain, which is mapped to a physical domain through a given coordinate transformation. By extending this concept to a multi-patch setting, simple and efficient numerical algorithms can be employed on relatively complex geometries. The main drawback of such an approach is the inherent difficulty in dealing with singularities of the coordinate transformation.This work suggests a comprehensive numerical strategy for the common situation of disk-like domains with a singularity at a unique pole, where one edge of the rectangular logical domain collapses to one point of the physical domain (for example, a circle). We present robust numerical methods for the solution of Vlasov-like hyperbolic equations coupled to Poisson-like elliptic equations in such geometries. We describe a semi-Lagrangian advection solver that employs a novel set of coordinates, named pseudo-Cartesian coordinates, to integrate the characteristic equations in the whole domain, including the pole, and a finite element elliptic solver based on globally C1 smooth splines (Toshniwal et al., 2017 [27]). The two solvers are tested both independently and on a coupled model, namely the 2D guiding-center model for magnetized plasmas, equivalent to a vorticity model for incompressible inviscid Euler fluids. The numerical methods presented show high-order convergence in the space discretization parameters, uniformly across the computational domain, without effects of order reduction due to the singularity. Dedicated tests show that the numerical techniques described can be applied straightforwardly also in the presence of point charges (equivalently, point-like vortices), within the context of particle-in-cell methods.

A semi-Lagrangian finite difference WENO scheme for scalar nonlinear conservation laws

For a nonlinear scalar conservation law in one-space dimension, we develop a locally conservative semi-Lagrangian finite difference scheme based on weighted essentially non-oscillatory reconstructions (SL-WENO). This scheme has the advantages of both WENO and semi-Lagrangian schemes. It is a locally mass conservative finite difference scheme, it is formally high-order accurate in space, it has small time truncation error, and it is essentially non-oscillatory. The scheme is nearly free of a CFL time step stability restriction for linear problems, and it has a relaxed CFL condition for nonlinear problems. The scheme can be considered as an extension of the SL-WENO scheme of Qiu and Shu (2011) [2] developed for linear problems. The new scheme is based on a standard sliding average formulation with the flux function defined using WENO reconstructions of (semi-Lagrangian) characteristic tracings of grid points. To handle nonlinear problems, we use an approximate, locally frozen trace velocity and a flux correction step. A special two-stage WENO reconstruction procedure is developed that is biased to the upstream direction. A Strang splitting algorithm is used for higher-dimensional problems. Numerical results are provided to illustrate the performance of the scheme and verify its formal accuracy. Included are applications to the Vlasov–Poisson and guiding-center models of plasma flow.

Optimization of particle-in-cell simulations for Vlasov-Poisson system with strong magnetic field

We study the dynamics of charged particles under the influence of a strong magnetic field by numerically solving the Vlasov-Poisson and guiding center models. By using appropriate data structures, we implement an efficient (from the memory access point of view) particle-in-cell method which enables simulations with a large number of particles. We present numerical results for classical one-dimensional Landau damping and two-dimensional Kelvin-Helmholtz test cases. The implementation also relies on a standard hybrid MPI/OpenMP parallelization. Code performance is assessed by the observed speedup and attained memory bandwidth.

Solving the guiding-center model on a regular hexagonal mesh

This paper introduces a Semi-Lagrangian solver for the Vlasov-Poisson equations on a uniform hexagonal mesh. The latter is composed of equilateral triangles, thus it doesn't contain any singularities, unlike polar meshes. We focus on the guiding-center model, for which we need to develop a Poisson solver for the hexagonal mesh in addition to the Vlasov solver. For the interpolation step of the Semi-Lagrangian scheme, a comparison is made between the use of Box-splines and of Hermite finite elements. The code will be adapted to more complex models and geometries in the future.

KINETIC EVOLUTION OF CORONAL HOLE PROTONS BY IMBALANCED ION-CYCLOTRON WAVES: IMPLICATIONS FOR MEASUREMENTS BY SOLAR PROBE PLUS

We extend the kinetic guiding-center model of collisionless coronal hole protons presented in Isenberg and Vasquez to consider driving by imbalanced spectra of obliquely propagating ion-cyclotron waves. These waves are assumed to be a small by-product of the imbalanced turbulent cascade to high perpendicular wavenumber, and their total intensity is taken to be 1% of the total fluctuation energy. We also extend the kinetic solutions for the proton distribution function in the resulting fast solar wind to heliocentric distances of 20 solar radii, which will be attainable by the Solar Probe Plus spacecraft. We consider three ratios of outward-propagating to inward-propagating resonant intensities: 1, 4, and 9. The self-consistent bulk flow speed reaches fast solar wind values in all cases, and these speeds are basically independent of the intensity ratio. The steady-state proton distribution is highly organized into nested constant-density shells by the resonant wave-particle interaction. The radial evolution of this kinetic distribution as the coronal hole plasma flows outward is understood as a competition between the inward- and outward-directed large-scale forces, causing an effective circulation of particles through the (v{sub ∥}, v{sub ⊥}) phase space and a characteristic asymmetric shape to the distribution. These asymmetries are substantial and persistmore » to the outer limit of the model computation, where they should be observable by the Solar Probe Plus instruments.« less

An Iteration Free Backward Semi-Lagrangian Scheme for Guiding Center Problems

In this paper, we develop an iteration free backward semi-Lagrangian method for nonlinear guiding center models. We apply the fourth-order central difference scheme for the Poisson equation and employ the local cubic interpolation for the spatial discretization. A key problem in the time discretization is to find the characteristic curve arriving at each grid point which is the solution of a system of highly nonlinear ODEs with a self-consistency imposed by the Poisson equation. The proposed method is based on the error correction method recently developed by the authors. For the error correction method, we introduce a modified Euler's polygon and solve the induced asymptotically linear differential equation with the midpoint quadrature rule to get the error correction term. We prove that the proposed iteration free method has convergence order at least $3$ in space and $2$ in time in the sense of the $L_{2}$-norm. In particular, it is shown that the proposed method has a good performance in computational...

A high order time splitting method based on integral deferred correction for semi-Lagrangian Vlasov simulations

Semi-Lagrangian schemes with various splitting methods, and with different reconstruction/interpolation strategies have been applied to kinetic simulations. For example, the order of spatial accuracy of the algorithms proposed in Qiu and Christlieb (2010) [29] is very high (as high as ninth order). However, the temporal error is dominated by the operator splitting error, which is second order for Strang splitting. It is therefore important to overcome such low order splitting error, in order to have numerical algorithms that achieve higher orders of accuracy in both space and time. In this paper, we propose to use the integral deferred correction (IDC) method to reduce the splitting error. Specifically, the temporal order accuracy is increased by r with each correction loop in the IDC framework, where r=1,2 for coupling the first order splitting and the Strang splitting, respectively. The proposed algorithm is applied to the Vlasov–Poisson system, the guiding center model, and two dimensional incompressible flow simulations in the vorticity stream-function formulation. We show numerically that the IDC procedure can automatically increase the order of accuracy in time. We also investigate numerical stability of the proposed algorithm via performing Fourier analysis to a linear model problem.

Gyroaverage effects on chaotic transport by drift waves in zonal flows

Finite Larmor radius (FLR) effects on E × B test particle chaotic transport in the presence of zonal flows is studied. The FLR effects are introduced by the gyro-average of a simplified E × B guiding center model consisting of the linear superposition of a non-monotonic zonal flow and drift waves. Non-monotonic zonal flows play a critical role on transport because they exhibit robust barriers to chaotic transport in the region(s) where the shear vanishes. In addition, the non-monotonicity gives rise to nontrivial changes in the topology of the orbits of the E × B Hamiltonian due to separatrix reconnection. The present study focuses on the role of FLR effects on these two signatures of non-monotonic zonal flows: shearless transport barriers and separatrix reconnection. It is shown that, as the Larmor radius increases, the effective zonal flow profile bifurcates and multiple shearless regions are created. As a result, the topology of the gyro-averaged Hamiltonian exhibits very complex separatrix reconnection bifurcations. It is also shown that FLR effects tend to reduce chaotic transport. In particular, the restoration of destroyed transport barriers is observed as the Larmor radius increases. A detailed numerical study is presented on the onset of global chaotic transport as function of the amplitude of the drift waves and the Larmor radius. For a given amplitude, the threshold for the destruction of the shearless transport barrier, as function of the Larmor radius, exhibits a fractal-like structure. The FLR effects on a thermal distribution of test particles are also studied. In particular, the fraction of confined particles with a Maxwellian distribution of gyroradii is computed, and an effective transport suppression is found for high enough temperatures.

A semi-Lagrangian AMR scheme for 2D transport problems in conservation form

In this paper, we construct a semi-Lagrangian (SL) Adaptive-Mesh-Refinement (AMR) solver for 1D and 2D transport problems in conservation form. First, we describe the à-la-Harten AMR framework: the adaptation process selects a hierarchical set of grids with different resolutions depending on the features of the integrand function, using as criteria the point value prediction via interpolation from coarser meshes, and the appearance of large gradients. We integrate in time by reconstructing at the feet of the characteristics through the Point-Value Weighted Essentially Non-Oscillatory (PV-WENO) interpolator. We propose, then, an extension to the 2D setting by making the time integration dimension-by-dimension thanks to a Strang splitting. We discuss the quality of the results and the speedup with respect to a Fixed Mesh (FM) strategy through the following benchmark tests: in 1D, constant and variable-coefficient advections; in 2D, the 1D Vlasov–Poisson system (2D in the phase space) for the case of constant-coefficient advections, and, for the case of variable-coefficient advections, the deformation flow and the guiding-center model.

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.

Fluid Moments in the Reduced Model for Plasmas with Large Flow Velocity

Representation of particle fluid moments in terms of fluid moments in the modified guiding-centre model for flowing plasmas with large E × B velocity [N. Miyato et al., J. Phys. Soc. Jpn. 78, 104501 (2009)] is derived from the formal exact representation by a perturbative expansion in the subsonic flow case. It is similar to that in the standard gyrokinetic model in the long wavelength limit, except it has an additional flow term. The flow term has no effect on the representation for particle density, leading to the same representation as the standard one formally. In the conventional guiding-centre models for flowing plasmas, on the other hand, the representation for particle density is different from the standard one. This is due to the difference in the transformation for the guiding-centre position. Although the exact representation usually used in the standard gyrokinetic model has a different form from that in the modified guiding-centre case, correspondence between the two models is shown by considering the alternative form of exact representation in the standard gyrokinetic case. The representation for particle density is also obtained from the single particle Lagrangian by a variational method which is used to derive the representation in the transonic case.

Numerical Solution of the Gyroaverage Operator for the Finite Gyroradius Guiding-Center Model

In this work, we are concerned with numerical approximation of the gyroav-erage operators arising in plasma physics to take into account the effects of the finite Larmor radius corrections. Several methods are proposed in the space configuration and compared to the reference spectral method. We then investigate the influence of the different approximations considering the coupling with some guiding-center models available in the literature.