Vlasov Solver Research Articles

Enabling technology for global 3D+3V hybrid-Vlasov simulations of near-Earth space

We present methods and algorithms that allow the Vlasiator code to run global, three-dimensional hybrid-Vlasov simulations of Earth's entire magnetosphere. The key ingredients that make Vlasov simulations at magnetospheric scales possible are the sparse velocity space implementation and spatial adaptive mesh refinement. We outline the algorithmic improvement of the semi-Lagrangian solver for six-dimensional phase space quantities, discuss the coupling of Vlasov and Maxwell equations' solvers in a refined mesh, and provide performance figures from simulation test runs that demonstrate the scalability of this simulation system to full magnetospheric runs.

Global linear stability analysis of kinetic trapped ion mode (TIM) in tokamak plasma using the spectral method

Trapped ion modes (TIMs) belong to the family of ion temperature gradient (ITG) modes, which are one of the important ingredients in heat turbulent transport at the ion scale in tokamak plasmas. A global linear analysis of a reduced gyro-bounce kinetic model for trapped particle modes is performed, and a spectral method is proposed to solve the dispersion relation. Importantly, the radial profile of the particle drift velocity is taken into account in the linear analysis by considering the magnetic flux ψ dependency of the equilibrium Hamiltonian in both the quasi-neutrality equation and equilibrium gyro-bounce averaged distribution function . Using this spectral method, linear growth rates of TIM instability in the presence of different temperature profiles and precession frequencies of trapped ions, with an approximated constant Hamiltonian and the exact ψ dependent equilibrium Hamiltonian, are investigated. The growth rate depends on the logarithmic gradient of temperature κ T , density κ n and equilibrium Hamiltonian . With the exact ψ dependent Hamiltonian, the growth rates and potential profiles are modified significantly, compared to the cases with an approximated constant Hamiltonian. All the results from the global linear analysis agree with a semi-Lagrangian based linear Vlasov solver with good accuracy. This spectral method is very fast and requires much less computation resources compared to a linear version of the Vlasov-solver based on a semi-Lagrangian scheme.

A parallel low-rank solver for the six-dimensional Vlasov–Maxwell equations

Continuum Vlasov simulations can be utilized for highly accurate modelling of fully kinetic plasmas. Great progress has been made recently regarding the applicability of the method in realistic plasma configurations. However, a reduction of the high computational cost that is inherent to fully kinetic simulations would be desirable, especially at high velocity space resolutions. For this purpose, low-rank approximations can be employed. The so far available low-rank solvers are restricted to either electrostatic systems or low dimensionality and can therefore not be applied to most space, astrophysical and fusion plasmas. In this paper we present a new parallel low-rank solver for the full six-dimensional electromagnetic Vlasov–Maxwell equations that can utilize distributed memory architectures. Special care is taken to ensure the conservation of mass and a good representation of Gauss's law. The low-rank Vlasov solver is applied to standard benchmark problems of plasma turbulence and magnetic reconnection and compared to the full grid method. It yields accurate results at significantly reduced computational cost.

An Energy Conserving Vlasov Solver That Tolerates Coarse Velocity Space Resolutions: Simulation of MMS Reconnection Events

AbstractVlasov solvers that operate on a phase‐space grid are highly accurate but also numerically demanding. Coarse velocity space resolutions, which are largely unproblematic in particle‐in‐cell (PIC) simulations, can lead to numerical heating or oscillations in continuum Vlasov methods. To address this issue, we present a new dual Vlasov solver which is based on an established positivity preserving advection scheme for the update of the distribution function and an energy conserving partial differential equation solver for the kinetic update of mean velocity and temperature. The solvers work together via moment fitting during which the maximum entropy part of the distribution function is replaced by the solution from the partial differential equation solver. This numerical scheme makes continuum Vlasov methods competitive with PIC methods concerning computational cost and enables us to model large scale reconnection in Earth's magnetosphere with a fully kinetic continuum method. The simulation results agree well with measurements by the MMS spacecraft.

Conservative semi-Lagrangian kinetic scheme coupled with implicit finite element field solver for multidimensional Vlasov Maxwell system

In this paper, a conservative semi-Lagrangian (CSL) kinetic scheme coupled with an implicit finite element method (IFEM) is developed for multidimensional electromagnetic plasma simulations. The present method (CSL-IFEM) enjoys the respective advantages of the CSL and IFEM, including mass conservation, high-order spatial accuracy, as well as being free from the Courant-Friedrichs-Lewy (CFL) condition. In present CSL-IFEM, the CSL scheme with a general high order positivity-preserving limiter is utilized for the spatial discretization of the Vlasov equation, which enables the proposed method to remove the CFL limitation in phase space, exactly conserve mass, and preserve the positivity of the distribution function. The IFEM solver is developed for the spatial discretization of the Maxwell equation, which makes the current method free from CFL limitation induced by the speed of light and easily handles general field boundary conditions. To make the proposed CSL-IFEM more efficient in multidimensional simulations, the dimensional splitting procedure and sparse matrix technology are implemented in CSL and IFEM, respectively. Then the Vlasov solver CSL and Maxwell solver IFEM are coupled via the current density calculated from the general Ohm law. Finally, several numerical experiments, including electromagnetic wave (2d0v), particle gyromotion (0d2v), streaming Weibel instability (1d2v) and magnetic reconnection (2d3v), are performed to demonstrate the capabilities of the proposed method.

High-fidelity kinetic modeling of instabilities and gyromotion physics in nonuniform low-beta plasmas

A fourth-order accurate continuum kinetic Vlasov solver and a systematic method for constructing customizable kinetic equilibria are demonstrated to be powerful tools for the study of nonuniform collisionless low-beta plasmas. The noise-free methodology is applied to investigate two gradient-driven instabilities in 4D (x,y,vx,vy) phase space: the Kelvin–Helmholtz instability and the lower hybrid drift instability. Nonuniform two-species configurations where ion gyroradii are comparable to gradient scale lengths are explored. The approach sheds light on the evolution of the pressure tensor in Kelvin–Helmholtz instabilities and demonstrates that the associated stress tensor deviates significantly from the gyroviscous stress tensor. Even at high magnetization, first-order approximations to finite-gyromotion physics are shown to be inadequate for the Kelvin–Helmholtz instability, as shear scales evolve to become on par with gyromotion scales. The methodology facilitates exploring transport and energy partitioning properties associated with lower hybrid drift instabilities in low-beta plasma configurations. Distribution function features are captured in detail, including the formation of local extrema in the vicinity of particle-wave resonances. The approach enables detailed targeted investigations and advances kinetic simulation capability for plasmas in which gyromotion plays an important role.

Phase-space structure of protohalos: Vlasov versus particle-mesh

The phase-space structure of primordial dark matter halos is revisited using cosmological simulations with three sine waves and cold dark matter (CDM) initial conditions. The simulations are performed with the tessellation based Vlasov solver ColDICE and a particle-mesh (PM) N-body code. The analyses include projected density, phase-space diagrams, radial density ρ(r), and pseudo-phase space density: Q(r) = ρ(r)/σv(r)3 with σv the local velocity dispersion. Particular attention is paid to force and mass resolution. Because the phase-space sheet complexity, estimated in terms of total volume and simplex (tetrahedron) count, increases very quickly, ColDICE can follow only the early violent relaxation phase of halo formation. During the violent relaxation phase, agreement between ColDICE and PM simulations having one particle per cell or more is excellent and halos have a power-law density profile, ρ(r) ∝ r−α, α ∈ [1.5, 1.8]. This slope, measured prior to any merger, is slightly larger than in the literature. The phase-space diagrams evidence complex but coherent patterns with clear signatures of self-similarity in the sine wave simulations, while the CDM halos are somewhat scribbly. After additional mass resolution tests, the PM simulations are used to follow the next stages of evolution. The power law progressively breaks down with a convergence of the density profile to the well-known Navarro–Frenk–White universal attractor, irrespective of initial conditions, that is even in the three-sine-wave simulations. This demonstrates again that mergers do not represent a necessary condition for convergence to the dynamical attractor. Not surprisingly, the measured pseudo phase-space density is a power law Q(r) ∝ r−αQ, with αQ close to the prediction of secondary spherical infall model, αQ ≃ 1.875. However this property is also verified during the early relaxation phase, which is non-trivial.

Study of inertial electrostatic confinement fusion using a finite-volume scheme for the one-dimensional Vlasov equation.

While the majority of fusion energy research is focused on magnetic confinement, there have been several alternative confinement methods aimed at the development of smaller and less expensive reactors. A number of these alternative reactors are based on a spherically convergent beam of recirculating ions and include designs such as inertial electrostatic confinement (IEC), multigrid IEC, and the periodically oscillating plasma sphere concept. Here, a fully time-dependent GPU-based Vlasov solver was developed in order to study these spherically convergent devices. This code solves the Vlasov equation for a spherically symmetric system using a finite-volume method with a modified flux to account for electrode transparency. The solver accounts for secondary electron emission, interactions between the charged particles, and collisional effects such as ionization and charge exchange. This code was used to investigate a system similar to the ion-injected device described by Hirsch (see [R. L. Hirsch, J. Appl. Phys. 38, 4522 (1967)10.1063/1.1709162]), who had reported a neutron production rate for deuterium-deuterium reactions in the range of 10^{6}to10^{7} neutrons per second, which was attributed to the formation of a virtual electrode structure near the center of the chamber. Attempts to reproduce this experiment [B. J. Egle, Ph.D. thesis, 2010] yielded similar fusion rates, though the majority of the reactions were found not to occur near the center of the chamber. The results of this Vlasov solver, considering only beam-beam and beam-background fusion reactions, show that beam-background reactions would be dominant in such an ion-injected device. This result is consistent with work by Baxter and Stuart, who proposed a simplified steady-state Boltzmann model. However, the result of both models are inconsistent with the experimental results, which indicate a higher neutron production rate, and an inverse pressure scaling trend. It is shown that the higher experimental rates may be explained by beam-target fusion between the ion beam and deuterium embedded on the inner surface of the cathode.

Intrabunch motion

Impedance-driven (but not only) coherent beam instabilities are usually studied analytically with the linearized Vlasov equation, ending up with an eigenvalue system to solve. The eigenvalues describe the beam oscillation mode-frequency shifts, leading in particular to intensity thresholds defined by the longitudinal mode coupling instability in the longitudinal plane and by the transverse mode coupling instability in the transverse plane in the absence of chromaticity. This can be directly compared to measurements in particular for the lowest modes and in the absence of tune spread. In the presence of nonlinearities or when higher-order modes are involved, this becomes quite difficult, if not impossible, and the coupling between the modes cannot be directly measured (or simulated) anymore. Another important observable is the intrabunch motion, which can be also accessed analytically thanks to the eigenvectors. To the author's knowledge, until now, the intrabunch signal has only been explained theoretically for independent longitudinal or transverse beam oscillation modes, i.e., when the bunch intensity is sufficiently low compared to the mode coupling threshold. It was never explained theoretically in detail when two (or more) modes are involved. For instance, no answers were already given to these questions: is (are) there some fixed point(s) when the transverse mode coupling instability starts? If yes, where is it (are they)? And what happens in the presence of mode decoupling? Any number of modes can be treated with the general approach discussed in this paper, which is based on the galactic Vlasov solver (which was previously successfully benchmarked against the pyheadtail macroparticle tracking code as concerns the beam oscillation mode-frequency shifts). However, to be able to clearly see what happens when the bunch intensity is increased, the simple case of two modes is discussed in detail. The purpose of this paper is to describe the different regimes, below, at, above the transverse mode coupling instability and also after the mode decoupling (as it happens sometimes), using a simple analytical model (where two modes are considered together), which helps to really understand what happens at each step. Better characterizing an instability is the first step before trying to find appropriate mitigation measures and push the performance of a particle accelerator. The evolution of the intrabunch motion with intensity is a fundamental observable with high-intensity high-brightness beams.

Longitudinal and transverse mode coupling instability: Vlasov solvers and tracking codes

The study of collective effects in circular accelerators can be tackled by solving numerically the Vlasov equation or by using tracking codes. The two methods are obtained with different approaches: Vlasov solvers consider a continuous distribution function and describe the beam with coherent oscillation modes in frequency domain (ending up usually with an eigenvalue system to solve), while tracking codes use macroparticles and wakefields in time domain. In this paper we present two Vlasov solvers for the study of collective effects (from impedances/wakefields only) which evaluate the frequency shift of coherent oscillation modes and possible mode coupling instability in the single-bunch case for both longitudinal and transverse planes. In the longitudinal plane the Vlasov solver also takes into account the potential well distortion due to the wakefields under some conditions. In parallel to this theoretical approach, tracking codes, which include collective effects, have been used as benchmark. In particular, starting from their results, we also propose a new method to study the frequency shift of coherent modes and compare it with the output of the Vlasov solvers.

An alternative view of coherent synchrotron radiation induced microbunching development in multibend recirculation arcs

Abstract Coherent synchrotron radiation (CSR) induced microbunching instability has been one of the most challenging issues in the design of modern accelerators. We recently developed a semi-analytical Vlasov solver and applied it to the investigation of the physical processes associated with microbunching amplification for several types of recirculation arc lattices. In this study, we build on the concept of stage gain, proposed by Huang and Kim (2002), and develop an alternative view to characterize CSR microbunching development in terms of stage orders. The stage orders enable a fair and quantitative comparison of the impact of lattice optics on the microbunching development for different lattices under similar initial beam parameters. The multibend arcs considered in this study have a distinguishing feature of multistage amplification. The CSR microbunching gain grows as a six-stage amplification for the presented recirculation arcs with a total of 24 dipoles, while two-stage amplification was previously concluded for a three-dipole bunch compressor chicane. We also attempt to relate the lattice optics pattern with the obtained stage gain functions through a physical interpretation. The investigation in this paper may serve as a complementary discussion to our previous work (Tsai et al., 2017).

Effect of crab cavity high order modes on the coupled-bunch stability of High-Luminosity Large Hadron Collider

High frequency High Order Modes can significantly affect the coupled-bunch stability in a circular accelerator. With a large enough shunt impedance they may drive a coupled-bunch transverse instability. We have developed an analytical model and implemented it in the NHT Vlasov solver. The results show that the instability is characterized by the excitation of many azimuthal intra-bunch modes, which have similar growth rates, that make the traditional remedies such as a flat resistive feedback and chromaticity inefficient in suppressing it. For the High Luminosity Large Hadron Collider this could result in a significant increase of the stabilizing Landau octupole current, up to $\sim 100$ A ($\sim 20\%$ of the maximum available current). In order to limit the increase below 10 A ($\sim 2\%$), the transverse shunt impedance has to be kept below 1 M$\Omega$/m.

Kinetic Solvers with Adaptive Mesh in Phase Space for Low-Temperature Plasmas

We describe the implementation of 1d1v and 1d2v Vlasov and Fokker-Planck kinetic solvers with adaptive mesh refinement in phase space (AMPS) and coupling these kinetic solvers to Poisson equation solver for electric fields. We demonstrate that coupling AMPS kinetic and electrostatic solvers can be done efficiently without splitting phase-space transport. We show that Eulerian fluid and kinetic solvers with dynamically adaptive Cartesian mesh can be used for simulations of collisionless plasma expansion into vacuum. The Vlasov-Fokker-Planck solver is demonstrated for the analysis of electron acceleration and scattering as well as the generation of runaway electrons in spatially inhomogeneous electric fields.

Gyrokinetic simulation of ITG turbulence with toroidal geometry including the magnetic axis by using field-aligned coordinates

Simulation domain in field-aligned coordinates of the electrostatic gyrokinetic nonlinear turbulence global code, NLT, is extended to include the magnetic axis. The artificial boundary near the magnetic axis is replaced by the natural boundary. The singularity at the magnetic axis in Vlasov solver is treated by considering the spatial relation of fixed grid points in field-aligned coordinates. A new Poisson’s equation solver is developed, the coefficient matrix of algebraic equations is derived by using Gauss’s theorem. Linear ITG modes test, Rosenbluth–Hinton test and nonlinear relaxation test of the ITG turbulence with adiabatic electrons are performed. The gyrocenter conservation is much improved by including the magnetic axis in the simulation domain.

A Quasi-Conservative Dynamical Low-Rank Algorithm for the Vlasov Equation

Numerical methods that approximate the solution of the Vlasov--Poisson equation by a low-rank representation have been considered recently. These methods can be extremely effective from a computational point of view, but contrary to most semi-Lagrangian or Eulerian Vlasov solvers, they do not conserve mass and momentum, neither globally nor in respecting the corresponding local conservation laws. This can be a significant limitation for intermediate and long time integration. In this paper we propose a numerical algorithm that overcomes some of these difficulties and demonstrate its utility by presenting numerical simulations.

Phase-space structure analysis of self-gravitating collisionless spherical systems

In the mean field limit, isolated gravitational systems often evolve towards a steady state through a violent relaxation phase. One question is to understand the nature of this relaxation phase, in particular the role of radial instabilities in the establishment/destruction of the steady profile. Here, through a detailed phase-space analysis based both on a spherical Vlasov solver, a shell code, and a N-body code, we revisit the evolution of collisionless self-gravitating spherical systems with initial power-law density profiles ρ(r) ∝ rn, 0 ≤ n ≤ −1.5, and Gaussian velocity dispersion. Two sub-classes of models are considered, with initial virial ratios η = 0.5 (“warm”) and η = 0.1 (“cool”). Thanks to the numerical techniques used and the high resolution of the simulations, our numerical analyses are able, for the first time, to show the clear separation between two or three well-known dynamical phases: (i) the establishment of a spherical quasi-steady state through a violent relaxation phase during which the phase-space density displays a smooth spiral structure presenting a morphology consistent with predictions from self-similar dynamics, (ii) a quasi-steady-state phase during which radial instabilities can take place at small scales and destroy the spiral structure but do not change quantitatively the properties of the phase-space distribution at the coarse grained level, and (iii) relaxation to a non-spherical state due to radial orbit instabilities for n ≤ −1 in the cool case.

Solving the Vlasov equation in two spatial dimensions with the Schrödinger method

We demonstrate that the Vlasov equation describing collisionless self-gravitating matter may be solved with the so-called Schr\"odinger method (ScM). With the ScM, one solves the Schr\"odinger-Poisson system of equations for a complex wave function in d dimensions, rather than the Vlasov equation for a 2d-dimensional phase space density. The ScM also allows calculating the d-dimensional cumulants directly through quasi-local manipulations of the wave function, avoiding the complexity of 2d-dimensional phase space. We perform for the first time a quantitive comparison of the ScM and a conventional Vlasov solver in d=2 dimensions. Our numerical tests were carried out using two types of cold cosmological initial conditions: the classic collapse of a sine wave and those of a gaussian random field as commonly used in cosmological cold dark matter N-body simulations. We compare the first three cumulants, that is, the density, velocity and velocity dispersion, to those obtained by solving the Vlasov equation using the publicly available code ColDICE. We find excellent qualitative and quantitative agreement between these codes, demonstrating the feasibility and advantages of the ScM as an alternative to N-body simulations. We discuss, the emergence of effective vorticity in the ScM through the winding number around the points where the wave function vanishes. As an application we evaluate the background pressure induced by the non-linearity of large scale structure formation, thereby estimating the magnitude of cosmological backreaction. We find that it is negligibly small and has time dependence and magnitude compatible with expectations from the effective field theory of large scale structure.

Application of High Dimensional B-Spline Interpolation in Solving the Gyro-Kinetic Vlasov Equation Based on Semi-Lagrangian Method

AbstractThe computation efficiency of high dimensional (3D and 4D) B-spline interpolation, constructed by classical tensor product method, is improved greatly by precomputing the B-spline function. This is due to the character of NLT code, i.e. only the linearised characteristics are needed so that the unperturbed orbit as well as values of the B-spline function at interpolation points can be precomputed at the beginning of the simulation. By integrating this fixed point interpolation algorithm into NLT code, the high dimensional gyro-kinetic Vlasov equation can be solved directly without operator splitting method which is applied in conventional semi-Lagrangian codes. In the Rosenbluth-Hinton test, NLT runs a few times faster for Vlasov solver part and converges at about one order larger time step than conventional splitting code.

Electromagnetic cyclotron instabilities in bi-Kappa distributed plasmas: A quasilinear approach

Anisotropic bi-Kappa distributed plasmas, as encountered in the solar wind and planetary magnetospheres, are susceptible to a variety of kinetic instabilities including the cyclotron instabilities driven by an excess of perpendicular temperature T⊥>T∥ (where ∥,⊥ denote directions relative to the mean magnetic field). These instabilities have been extensively investigated in the past, mainly limiting to a linear stability analysis. About their quasilinear (weakly nonlinear) development, some insights have been revealed by numerical simulations using PIC and Vlasov solvers. This paper presents a self-consistent analytical approach, which provides for both the electron and proton cyclotron instabilities an extended picture of the quasilinear time evolution of the anisotropic temperatures as well as the wave energy densities.

Conditions for coherent-synchrotron-radiation-induced microbunching suppression in multibend beam transport or recirculation arcs

The coherent synchrotron radiation (CSR) of a high-brightness electron beam traversing a series of dipoles, such as transport or recirculation arcs, may result in beam phase space degradation. On one hand, CSR can perturb electron transverse motion in dispersive regions along the beam line and possibly cause emittance growth. On the other hand, the CSR effect on the longitudinal beam dynamics could result in microbunching instability. For transport arcs, several schemes have been proposed to suppress the CSR-induced emittance growth. Correspondingly, a few scenarios have been introduced to suppress CSR-induced microbunching instability, which however mostly aim for linac-based machines. In this paper we provide sufficient conditions for suppression of CSR-induced microbunching instability along transport or recirculation arcs. Examples are presented with the relevant microbunching analyses carried out by our developed semianalytical Vlasov solver [C.-Y. Tsai, D. Douglas, R. Li, and C. Tennant, Linear microbunching analysis for recirculation machines, Phys. Rev. ST Accel. Beams 19, (2016)]. The example lattices include low-energy ($\ensuremath{\sim}100\text{ }\text{ }\mathrm{MeV}$) and high-energy ($\ensuremath{\sim}1\text{ }\text{ }\mathrm{GeV}$) recirculation arcs, and medium-energy compressor arcs. Our studies show that lattices satisfying the proposed conditions indeed have microbunching gain suppressed. Beam current dependences of maximal CSR microbunching gains are also demonstrated, which should help outline a beam line design for different scales of nominal currents. We expect this analysis can shed light on the lattice design approach that aims to control the CSR-induced microbunching.