We simulate the dynamics, including laser cooling, of three-dimensional (3-D) ion crystals confined in a Penning trap using a newly developed molecular dynamics-like code. The numerical integration of the ions’ equations of motion is accelerated using the fast multipole method to calculate the Coulomb interaction between ions, which allows us to efficiently study large ion crystals with thousands of ions. In particular, we show that the simulation time scales linearly with ion number, rather than with the square of the ion number. By treating the ions’ absorption of photons as a Poisson process, we simulate individual photon scattering events to study laser cooling of 3-D ellipsoidal ion crystals. Initial simulations suggest that these crystals can be efficiently cooled to ultracold temperatures, aided by the mixing of the easily cooled axial motional modes with the low frequency planar modes. In our simulations of a spherical crystal of 1000 ions, the planar kinetic energy is cooled to several millikelvin in a few milliseconds while the axial kinetic energy and total potential energy are cooled even further. This suggests that 3-D ion crystals could be well suited as platforms for future quantum science experiments.

We present a quantum algorithm based on repeated measurement to solve initial-value problems for nonlinear ordinary differential equations (ODEs), which may be generated from partial differential equations in plasma physics. We map a dynamical system to a Hamiltonian form, where the Hamiltonian matrix is a function of dynamical variables. To advance in time, we measure expectation values from the previous time step and evaluate the Hamiltonian function classically, which introduces stochasticity into the dynamics. We then perform standard quantum Hamiltonian simulation over a short time, using the evaluated constant Hamiltonian matrix. This approach requires evolving an ensemble of quantum states, which are consumed each step to measure the required observables. We apply this approach to the classic logistic and Lorenz systems, in both integrable and chaotic regimes. Our analysis shows that the solutions’ accuracy is influenced by both the stochastic sampling rate and the nature of the dynamical system.

This paper presents a novel approach for simulating plasma instabilities in tokamak plasmas using the piecewise field-aligned finite element method in combination with the particle-in-cell method. Our method traditionally aligns the computational grid, but defines the basis functions in piecewise field-aligned coordinates to avoid grid deformation while naturally representing the field-aligned mode structures. This scheme is formulated and implemented numerically. It also applies to the unstructured triangular meshes in principle. We have conducted linear benchmark tests, which agree well with previous results and traditional schemes. Furthermore, multiple-$n$ simulations are also carried out as a proof of principle, demonstrating the efficiency of this scheme in nonlinear turbulence simulations within the framework of the finite element method.

We analyse distributions of the spatial scales of coherent intermittent structures – current sheets – obtained from fully kinetic, two-dimensional simulations of relativistic turbulence in a collisionless pair plasma using unsupervised machine-learning data dissection. We find that the distribution functions of sheet length $\ell$ (longest scale of the analysed structure in the direction perpendicular to the dominant guide field) and curvature $r_c$ (radius of a circle fitted to the structures) can be well-approximated by power-law distributions, indicating self-similarity of the structures. The distribution for the sheet width $w$ (shortest scale of the structure) peaks at the kinetic scales and decays exponentially at larger values. The data shows little or no correlation between $w$ and $\ell$, as expected from theoretical considerations. The typical $r_c$ depends linearly on $\ell$, which indicates that the sheets all have a similar curvature relative to their sizes. We find a weak correlation between $r_c$ and $w$. Our results can be used to inform realistic magnetohydrodynamic subgrid models for plasma turbulence in high-energy astrophysics.

Two-dimensional particle-in-cell (PIC) simulations explore the collisionless tearing instability developing in a Harris equilibrium configuration in a pair (electron–positron) plasma, with no guide field, for a range of parameters from non-relativistic to relativistic temperatures and drift velocities. Growth rates match the predictions of Zelenyi & Krasnosel'skikh (Astron. Zh., vol. 56, 1979, pp. 819–832) modified for relativistic drifts by Hoshino (Astrophys. J., vol. 900, issue 1, 2020, p. 66) as long as the assumption holds that the thickness $a$ of the current sheet is larger than the Larmor radius $\rho _L$, with the fastest growing mode at $ka \approx 1/\sqrt {3}$. Aside from confirming these predictions, we explore the transitions from thick to thin current sheets and from classical to relativistic temperatures. We show that for thinner current sheets ($a< \rho _L$), the growth rate matches the prediction for the case $a=\rho _L$. We also explore the nonlinear evolution of the modes. While the wavenumber with the fastest growth rate initially matches the prediction of Zelenyi & Krasnosel'skikh (1979), these modes saturate moving the dominant mode to lower wavenumbers (especially for thick current sheets with low growth rates). Furthermore, at a late, nonlinear stage, the growth rate (initially following the growth rate prediction proportional to $(\rho _L/a)^{3/2} < 1$) increases faster than exponentially, reaching a maximum growth rate equivalent to the linear growth rate prediction at $\rho _L/a = 1$, before eventually saturating.

Quantum algorithms have been proposed to accelerate the simulation of the chaotic dynamical systems that are ubiquitous in the physics of plasmas. Quantum computers without error correction might even use noise to their advantage to calculate the Lyapunov exponent by measuring the Loschmidt echo fidelity decay rate. For the first time, digital Hamiltonian simulations of the quantum sawtooth map, performed on the IBM-Q quantum hardware platform, show that the fidelity decay rate of a digital quantum simulation increases during the transition from dynamical localization to chaotic diffusion in the map. The observed error per CNOT gate increases by $1.5{\times }$ as the dynamics varies from localized to diffusive, while only changing the phases of virtual RZ gates and keeping the overall gate count constant. A gate-based Lindblad noise model that captures the effective change in relaxation and dephasing errors during gate operation qualitatively explains the effect of dynamics on fidelity as being due to the localization and entanglement of the states created. Specifically, highly delocalized states that are entangled with random phases show an increased sensitivity to dephasing and, on average, a similar sensitivity to relaxation as localized states. In contrast, delocalized unentangled states show an increased sensitivity to dephasing but a lower sensitivity to relaxation. This gate-based Lindblad model is shown to be a useful benchmarking tool by estimating the effective Lindblad coherence times during CNOT gates and finding a consistent $2\unicode{x2013}3{\times }$ shorter $T_2$ time than reported for idle qubits. Thus, the interplay of the dynamics of a simulation with the noise processes that are active can strongly influence the overall fidelity decay rate.

The heat conductivity of a plasma is usually much higher along the magnetic field than across it, and, as a result, the presence of a magnetic island can significantly affect the temperature profile in its vicinity. Radiation energy losses, which depend sensitively on temperature, are thus strongly affected by magnetic islands. This phenomenon is explored in a simple mathematical setting, and it is shown that the presence of a magnetic island greatly enhances a plasma's capacity to radiate energy. In the limit of highly anisotropic heat conductivity, the steady-state heat conduction equation can be reduced to an ordinary differential equation. Although this equation operates in one dimension, the topology is not that of the real line, but corresponds to a rod with a cooling fin. As parameters such as the incoming heat flux or the radiation amplitude are varied, the radiation has a tendency to linger around the island, in particular in the region of the separatrix, and the total radiated energy is then significantly increased. The island acts as a ‘cooling fin’ to the plasma. Furthermore, the solutions exhibit bifurcations, where the location of the radiation zone suddenly changes.

Electron cyclotron resonance ion thrusters (ECRITs) have the potential to be used for space gravitational wave detection due to their wide thrust range. However, an unclear understanding of dynamic processes of ECRITs with strongly coupled multi-operating parameters limits further improvements on thrust noise and response velocity by feedback control systems. An integrative mathematical model considering the non-Maxwell electron energy distribution function for ECRITs is validated by experiments and used to study the influence of operating parameters on the dynamic processes of thrusters, which provides a new simplified grid model. Simulation results show the response processes with microwave (MW) power can be divided into two stages. The characteristic times of the first and second stages are respectively several microseconds and 10 ms, which are respectively dominated by plasma motion and the volume effect. The overshoot of screen grid (SG) current decreases, while its response time remains unchanged when the response time of MW power is prolonged. The response time of SG current with a step increase of flow rate is approximately 10 ms, consistent with the volume effect. The SG current decreases with rise of flow rate for high flow rate operations due to the small increment of ion density limited by low electron temperature, the decrease of ion Bohm velocity and reduction of sheath extraction area. The influence of grid voltage on the dynamic process of the SG current depends on variation ranges of extraction capabilities. When variations of sheath extraction area are limited, the response time is 5 μs, consistent with plasma response time. It is prolonged to 0.5 ms if sheath extraction area variations are large because they cause obvious variations of plasma parameters in the discharge chamber. These dynamic results can not only facilitate designing feedback controllers of micro-propulsion systems for high-precision space missions, but also provide guidance for ion sources to generate highly stable or rapid-response ion beam.

Confinement quality in fusion plasma is influenced significantly by the presence of heavy impurities, which can lead to radiative heat loss and reduced confinement. This study explores the clustering of heavy impurity, i.e. tungsten in edge plasma, using high-resolution direct numerical simulations of the Hasegawa–Wakatani equations. We use the Stokes number to quantify the inertia of impurity particles. It is found that particle inertia will cause spatial intermittency in particle distribution and the formation of large-scale structures, i.e. the clustering of particles. The degrees of clustering are influenced by the Stokes number. To quantify these observations, we apply a modified Voronoi tessellation, which assigns specific volumes to impurity particles. By determining time changes of these volumes, we can calculate the impurity velocity divergence, which allows the clustering dynamics to be assessed. To quantify the clustering statistically, several approaches are applied, such as probability density function (PDF) of impurity velocity divergence and joint PDF of volume and divergence.

MagNetUS is a network of scientists and research groups that coordinates and advocates for fundamental magnetized plasma research in the USA. Its primary goal is to bring together a broad community of researchers and the experimental and numerical tools they use in order to facilitate the sharing of ideas, resources and common tasks. Discussed here are the motivation and goals for this network and details of its formation, history and structure. An overview of associated experimental facilities and numerical projects is provided, along with examples of scientific topics investigated therein. Finally, a vision for the future of the organization is given.

Turbulence self-organization is studied in the flux-driven regime by means of the reduced model Tokam1D. Derived in the electrostatic and isothermal limit but keeping finite electron and ion temperatures, it features two instabilities that are suspected to dominate turbulent transport at the edge of L-mode tokamak plasmas: interchange (a reduced version of the resistive ballooning modes) and collisional drift waves, governed respectively by an effective gravity parameter $g$ and the adiabaticity parameter $C$. The usual properties of these two instabilities are recovered in the linear regime. The nonlinear study focuses on the self-organization of collisional drift-wave turbulence at $g=0$. It is found that the energy stored in zonal flows (ZFs) decreases smoothly at small $C$ due to the reduction of both electric and diamagnetic stresses. Conversely to gradient-driven simulations, no sharp collapse is observed due to the self-consistent evolution of the equilibrium density profile. The ZFs are found to structure into staircases at small and large $C$. These structures exhibit a rich variety of dynamics but are found to be robust to large perturbations. Their nucleation is found to be critically governed by the phase dynamics. Finally, the staircase structures are lost in the gradient-driven regime, when the system is prevented from storing turbulent energy into the equilibrium density (pressure) profile.

Isolated, undamped geodesic-acoustic-mode (GAM) packets have been demonstrated to obey a (focusing) nonlinear Schrödinger equation (NLSE) (E. Poli, Phys. Plasmas, 2021). This equation predicts susceptibility of GAM packets to the modulational instability (MI). The necessary conditions for this instability are analysed analytically and numerically using the NLSE model. The predictions of the NLSE are compared with gyrokinetic simulations performed with the global particle-in-cell code ORB5, where GAM packets are created from initial perturbations of the axisymmetric radial electric field $E_r$. An instability of the GAM packets with respect to modulations is observed both in cases in which an initial perturbation is imposed and when the instability develops spontaneously. However, significant differences in the dynamics of the small scales are discerned between the NLSE and gyrokinetic simulations. These discrepancies are mainly due to the radial dependence of the strength of the nonlinear term, which we do not retain in the solution of the NLSE, and to the damping of higher radial spectral components $k_r$. The damping of the high-$k_r$ components, which develop as a consequence of the nonlinearity, can be understood in terms of Landau damping. The influence of the ion Larmor radius $\rho _i$ as well as the perturbation wavevector $k_\text {pert}$ on this effect is studied. For the parameters considered here the aforementioned damping mechanism hinders the MI process significantly from developing to its full extent and is strong enough to stabilize some of the (according to the undamped NLSE model) unstable wavevectors.

A quasi-linear reduced transport model is developed from a database of high-$\beta$ electromagnetic nonlinear gyrokinetic simulations performed with spherical tokamak for energy production (STEP) relevant parameters. The quasi-linear model is fully electromagnetic and accounts for the effect of equilibrium flow shear using a novel approach. Its flux predictions are shown to agree quantitatively with predictions from local nonlinear gyrokinetic simulations across a broad range of STEP-relevant local equilibria. This reduced transport model is implemented in the T3D transport solver that is used to perform the first flux-driven simulations for STEP to account for transport from hybrid kinetic ballooning mode turbulence, which dominates over a wide region of the core plasma. Nonlinear gyrokinetic simulations of the final transport steady state from T3D return turbulent fluxes that are consistent with the reduced model, indicating that the quasi-linear model may also be appropriate for describing the transport steady state. Within the assumption considered here, our simulations support the existence of a transport steady state in STEP with a fusion power comparable to that in the burning flat top of the conceptual design, but do not demonstrate how this state can be accessed.

Using two counter-propagating ultra-intense laser interactions with a solid target, we conducted a study on the generation of electron-positron pairs via the multi-photon Breit–Wheeler (BW) process and trident process. These processes were simulated using the particle-in-cell (PIC) code EPOCH. Our proposed scheme involves irradiating two targets with two counter-propagating lasers. High-energy photons are produced when hot electrons collide with the reflected laser pulse at the target's front, leading to electron and positron pair production. In the single-target scenario, electron bunches are extracted from the target by the p-polarized laser electromagnetic field and accelerated by the laser ponderomotive force before colliding with the counter-propagating laser. However, using two targets enhances pair creation compared with the single-target set-up. We observed that in two-target configurations, the increased number of high-energy gamma-rays contributes to higher-energy electron–positron generation. Additionally, the generation of hot electrons is also more pronounced in this scheme. Consequently, the laser demonstrates higher efficiency in generating gamma photons and positrons in the dual-target set-up, which is beneficial for investigating high-energy pair production and gamma-ray emission. The generated positrons exhibit a density of the order of $10^{27}\,\text {m}^{-3}$ and can be accelerated to energies of 1.5 GeV. The involvement of hot electrons in the target is crucial for generating high-energy photons and positrons. The maximum pair yield reaches $8 \times 10^9$ for the BW process and $10^8$ for the trident process. Notably, the total laser energy conversion efficiencies to electrons, $\gamma$-rays and positrons show improvement in the dual-target configuration. Specifically, the laser energy absorbed by positrons increases from 11.62 % in Case A to 13.12 % in Case B. These enhancements in conversion efficiency and electron/positron density have significant practical implications in experimental set-ups. In both the BW and trident processes, the two-target set-up dominates, highlighting its effectiveness. We also compared the strengths of both approaches, suggesting that these simple models of implementing two targets can be used in experiments as well.

Magnetic reconnection leads to the formation of island-shaped magnetic structure(s). Due to disagreement between theoretical evaluations of the characteristic reconnection time and observations, it is commonly accepted that the collisionality (or resistivity) is too low to explain magnetic reconnection phenomena in fusion plasmas. Thus, magnetic reconnection still raises many open questions. The work presented here aims to improve the fundamental knowledge about ‘the life of a magnetic island’. Here, in the light of the many works of the last 70 years, a new paradigm for understanding magnetic reconnection in fusion plasmas is proposed. The life of a magnetic island (whatever its scale) follows three phases: the origin, the growth and the saturation. The possible physical mechanisms at play in these three phases will be investigated. First, for the island origin, typical time scales in link with magnetic reconnection will be evaluated for three tokamaks of different sizes (TCV, WEST and JET) to verify if magnetic reconnection is such an unexplained phenomenon in fusion plasmas. Second, for the island drive, the richness of possible mechanisms leading to ‘rapid’ magnetic island growth in fusion devices will be presented for small and large scales. Third comes the island saturation step. Results on the prediction of a large island width at saturation are presented and discussed.

The regulation of electron heat transport in high-$\beta$, weakly collisional, magnetized plasma is investigated. A temperature gradient oriented along a mean magnetic field can induce a kinetic heat-flux-driven whistler instability (HWI), which back-reacts on the transport by scattering electrons and impeding their flow. Previous analytical and numerical studies have shown that the heat flux for the saturated HWI scales as $\beta _e^{-1}$. These numerical studies, however, had limited scale separation and consequently large fluctuation amplitudes, which calls into question their relevance at astrophysical scales. To this end, we perform a series of particle-in-cell simulations of the HWI across a range of $\beta _e$ and temperature-gradient length scales under two different physical set-ups. The saturated heat flux in all of our simulations follows the expected $\beta _e^{-1}$ scaling, supporting the robustness of the result. We also use our simulation results to develop and implement several methods to construct an effective collision operator for whistler turbulence. The results point to an issue with the standard quasi-linear explanation of HWI saturation, which is analogous to the well-known $90^{\circ }$ scattering problem in the cosmic-ray community. Despite this limitation, the methods developed here can serve as a blueprint for future work seeking to characterize the effective collisionality caused by kinetic instabilities.

Ion cyclotron resonance heating is a versatile heating method that has been demonstrated to be able to efficiently couple power directly to the ions via the fast magnetosonic wave. However, at temperatures relevant for reactor grade devices such as DEMO, electron damping becomes increasingly important. To reduce electron damping, it is possible to use an antenna with a power spectrum dominated by low parallel wavenumbers. Moreover, using an antenna with a unidirectional spectrum, such as a travelling wave array antenna, the parallel wavenumber can be downshifted by mounting the antenna in an elevated position relative to the equatorial plane. This downshift can potentially enhance ion heating as well as fast wave current drive efficiency. Thus, such a system could benefit ion heating during the ramp-up phase and be used for current drive during flat-top operation. To test this principle, both ion heating and current drive have been simulated in a DEMO-like plasma for a few different mounting positions of the antenna using the FEMIC code. We find that moving the antenna off the equatorial plane makes ion heating more efficient for all considered plasma temperatures at the expense of on-axis heating. Moreover, although current drive efficiency is enhanced, electron damping is reduced for lower mode numbers, thus reducing the driven current in this part of the spectrum.

Fully relativistic particle-in-cell (PIC) simulations are crucial for advancing our knowledge of plasma physics. Modern supercomputers based on graphics processing units (GPUs) offer the potential to perform PIC simulations of unprecedented scale, but require robust and feature-rich codes that can fully leverage their computational resources. In this work, this demand is addressed by adding GPU acceleration to the PIC code Osiris. An overview of the algorithm, which features a CUDA extension to the underlying Fortran architecture, is given. Detailed performance benchmarks for thermal plasmas are presented, which demonstrate excellent weak scaling on NERSC's Perlmutter supercomputer and high levels of absolute performance. The robustness of the code to model a variety of physical systems is demonstrated via simulations of Weibel filamentation and laser-wakefield acceleration run with dynamic load balancing. Finally, measurements and analysis of energy consumption are provided that indicate that the GPU algorithm is up to ${\sim }$14 times faster and $\sim$7 times more energy efficient than the optimized CPU algorithm on a node-to-node basis. The described development addresses the PIC simulation community's computational demands both by contributing a robust and performant GPU-accelerated PIC code and by providing insight into efficient use of GPU hardware.

We construct a mean-field model that describes the nonlinear dynamics of a spin-polarized electron gas interacting with fixed, positively charged ions possessing a magnetic moment that evolves in time. The mobile electrons are modelled by a four-component distribution function in the two-dimensional phase space $(x,v)$, obeying a Vlasov–Poisson set of equations. The ions are modelled by a Landau–Lifshitz equation for their spin density, which contains ion–ion and electron–ion magnetic exchange terms. We perform a linear response study of the coupled Vlasov–Poisson–Landau–Lifshitz (VPLL) equations for the case of a Maxwell–Boltzmann equilibrium, focussing in particular on the spin dispersion relation. Conditions of stability or instability for the spin modes are identified, which depend essentially on the electron spin polarization rate $\eta$ and the electron–ion magnetic coupling constant $K$. We also develop an Eulerian grid-based computational code for the fully nonlinear VPLL equations, based on the geometric Hamiltonian method first developed by Crouseilles et al. (J. Plasma Phys., vol. 89, no. 2, 2023, p. 905890215). This technique allows us to achieve great accuracy for the conserved quantities, such as the modulus of the ion spin vector and the total energy. Numerical tests in the linear regime are in accordance with the estimations of the linear response theory. For two-stream equilibria, we study the interplay of instabilities occurring in both the charge and the spin sectors. The set of parameters used in the simulations, with densities close to those of solids (${\approx }10^{29}\ {\rm m}^{-3}$) and temperatures of the order of 10 eV, may be relevant to the warm dense matter regime appearing in some inertial fusion experiments.

This paper presents a theoretical model of plasma equilibrium in the diamagnetic confinement mode in an axisymmetric mirror device with neutral beam injection. The hot ionic component is described within the framework of the kinetic theory, since the Larmor radius of the injected ions appears to be comparable to or even larger than the characteristic scale of the magnetic field inhomogeneity. The electron drag of the hot ions is taken into account, while the angular scattering of the hot ions due to Coulomb collisions is neglected. The background warm plasma, on the contrary, is considered to be in local thermal equilibrium, i.e. has a Maxwellian distribution function and is described in terms of magnetohydrodynamics. The density of the hot ions is assumed to be negligible compared with that of the warm plasma. Both the conventional gas-dynamic loss and the non-adiabatic loss specific to the diamagnetic confinement mode are taken into account. In this work, we do not consider the effects of the warm plasma rotation as well as the inhomogeneity of the electrostatic potential. A self-consistent theoretical model of the plasma equilibrium is constructed. In the case of the cylindrical bubble, this model is reduced to a simpler one. The numerical solutions in the limit of a thin transition layer of the diamagnetic bubble are found. Examples of the equilibria corresponding to the gas-dynamic multiple-mirror trap device are considered.