The generation and radial structure of zonal flows are studied in competing collisional drift waves and interchange turbulence using the reduced flux-driven nonlinear model Tokam1D. Zonal flows are generated in both the interchange dominated and adiabatic regimes with the former favoring radially structured flows and avalanche transport. The distance to the instability threshold proves to be key, with a more stable radial flow structure emerging near the threshold and increased energy stored in the flows for interchange turbulence. The avalanches are shown to perturb zonal flow structures in drift-wave turbulence and to reactivate them in the interchange regime. Finally, the ExB staircases with radially structured, stable in time zonal flows are proved beneficial for the overall confinement.

We present two-dimensional global simulations of mitigated and vertically unstable disruptions in ITER in the presence of runaway electrons (REs). An elongated plasma in free-boundary equilibrium is subjected to an artificial thermal quench (TQ) and current profile flattening, followed by one or more massive material injections and a RE avalanche. Scenarios of major disruptions as well as upward and downward vertical displacement events are considered. Results provide important insights into the effects of RE formation, post TQ current profile, injection quantities and timings, and impurity flushout on the overall evolution of disruption and the plasma vertical motion thereof. Interplay between the various effects offers scope for potentially beneficial RE mitigation scenarios.

The linear collisionless plasma response to a zonal-density perturbation in quasisymmetric stellarators is studied, including the geodesic-acoustic-mode oscillations and the Rosenbluth–Hinton residual flow. While the geodesic-acoustic-mode oscillations in quasiaxisymmetric configurations are similar to tokamaks, they become non-existent in quasi-helically symmetric configurations when the effective safety factor in helical-angle coordinates is small. Compared with concentric-circular tokamaks, the Rosenbluth–Hinton residual is also found to be multiplied by a geometric factor $\mathcal {C}$ that arises from the flux-surface-averaged classical polarization. Using the near-axis-expansion framework, we derive an analytic expression for $\mathcal {C}$, which varies significantly among different configurations. These analytic results are compared with numerical simulation results from the global gyrokinetic particle-in-cell code GTC, and good agreement with the theoretical Rosenbluth–Hinton residual level is achieved when the quasisymmetry error is small enough.

Mercier’s criterion is typically enforced as a hard operational limit in stellarator design. At the same time, past experimental and numerical studies have shown that this limit may often be surpassed, though the exact mechanism behind this nonlinear stability is not well understood. This work aims to contribute to our current understanding by comparing the nonlinear evolution of Mercier unstable Wendelstein stellarators with that of nonlinearly stable quasi-interchange modes in tokamaks. A high mirror, very low $\iota$, W7-X-like configuration is first simulated. Broad flow structures are observed, which produce a similar magnetohydrodynamic (MHD) dynamo term to that in hybrid tokamak discharges, leading to flux pumping. Unlike in tokamaks, there is no net toroidal current to counterbalance this dynamo, and it is unclear if it can be sustained to obtain a similar quasistationary nonlinear state. In the simulation, partial reconnection induced by the overlap of multiple interchange instabilities leads to a core temperature crash. A second case is then considered using experimental reconstructions of intermediate $\beta$ W7-AS discharges, where saturated low-n modes were observed experimentally, with sustained MHD signatures over tens of milliseconds. It is shown that these modes do not saturate in a benign quasistationary way in current simulations even in the presence of background equilibrium $\boldsymbol{E} \times \boldsymbol{B}$ flow shear. This leads to a burst of MHD behaviour, inconsistent with the sustained MHD signatures in the experiment. Nevertheless, the (1, 2) mode is observed at the experimental Spitzer resistivity, and its induced anomalous transport can be overcome using an experimentally relevant heat source, reproducing these aspects of the dynamics. The possible reasons for the discrepancies between experiment and simulation, and the observation of partial reconnection in contrast to flux pumping are discussed, in view of reproducing and designing for operation of stellarators beyond the Mercier stability limit.

This work provides an optimisation mechanism to ensure the compatibility of non-planar stellarator coils with ReBCO (rare-earth barium copper oxide) high-temperature superconducting (HTS) tape. ReBCO coils enable higher field strengths and/or operating temperatures for the magnet systems of future fusion reactors, but they are sensitive to mechanical strain due to their brittle, ceramic functional layer. To ensure that non-planar coils can be wound without damage, we have introduced into the stellarator optimisation framework SIMSOPT a penalty on the binormal curvature and torsion of the tape. This metric can be used to optimise the tape winding orientation along a given coil filament or the coil filament itself can also be free to vary as part of the strain optimisation. We demonstrate the strain optimisation in three examples. For the EPOS (electrons and positrons in an optimised stellarator) design, we combined the strain penalty with an objective for quasisymmetry into a single-stage optimisation; this enables us to find a configuration with excellent quasisymmetry at the smallest possible size compatible with the use of ReBCO tape. For CSX (Columbia stellarator experiment), in addition to HTS strain, we added a penalty to prevent full turn tape rotation, so as to ease the coil winding process. For a coil at reactor scale, we found a considerable variation of the binormal and torsional strain over the cross-section of the large winding pack (54 cm x 54 cm); by exploiting the overall orientation of the winding pack as a degree of freedom, we were able to reduce strains below limits for all of the ReBCO stacks in the pack.

In magnetized, stratified environments such as the Sun's corona and solar wind, Alfvénic fluctuations ‘reflect’ from background gradients, enabling nonlinear interactions that allow their energy to dissipate into heat. This process, termed ‘reflection-driven turbulence’, likely plays a key role in coronal heating and solar-wind acceleration, explaining a range of detailed observational correlations and constraints. Building on previous works focused on the inner heliosphere, here we study the basic physics of reflection-driven turbulence using reduced magnetohydrodynamics in an expanding box – the simplest model that can capture local turbulent plasma dynamics in the super-Alfvénic solar wind. Although idealized, our high-resolution simulations and simple theory reveal a rich phenomenology that is consistent with a diverse range of observations. Outwards-propagating fluctuations, which initially have high imbalance (high cross-helicity), decay nonlinearly to heat the plasma, becoming more balanced and magnetically dominated. Despite the high imbalance, the turbulence is strong because Elsässer collisions are suppressed by reflection, leading to ‘anomalous coherence’ between the two Elsässer fields. This coherence, together with linear effects, causes the growth of ‘anastrophy’ (squared magnetic potential) as the turbulence decays, forcing the energy to rush to larger scales and forming a ‘$1/f$-range’ energy spectrum in the process. Eventually, expansion overcomes the nonlinear and Alfvénic physics, forming isolated, magnetically dominated ‘Alfvén vortices’ with minimal nonlinear dissipation. These results can plausibly explain the observed radial and wind-speed dependence of turbulence imbalance (cross-helicity), residual energy, fluctuation amplitudes, plasma heating and fluctuation spectra, as well as making a variety of testable predictions for future observations.

Magnetic geometry has a significant effect on the level of turbulent transport in fusion plasmas. Here, we model and analyse this dependence using multiple machine learning methods and a dataset of ${\gt}200\,000$ nonlinear gyrokinetic simulations of ion-temperature-gradient turbulence in diverse non-axisymmetric geometries. The dataset is generated using a large collection of both optimised and randomly generated stellarator equilibria. At fixed gradients and other input parameters, the turbulent heat flux varies between geometries by several orders of magnitude. Trends are apparent among the configurations with particularly high or particularly low heat flux. Regression and classification techniques from machine learning are then applied to extract patterns in the dataset. Due to a symmetry of the gyrokinetic equation, the heat flux and regressions thereof should be invariant to translations of the raw features in the parallel coordinate, similar to translation invariance in computer vision applications. Multiple regression models including convolutional neural networks (CNNs) and decision trees can achieve reasonable predictive power for the heat flux in held-out test configurations, with highest accuracy for the CNNs. Using Spearman correlation, sequential feature selection and Shapley values to measure feature importance, it is consistently found that the most important geometric lever on the heat flux is the flux surface compression in regions of bad curvature. The second most important geometric feature relates to the magnitude of geodesic curvature. These two features align remarkably with surrogates that have been proposed based on theory, while the methods here allow a natural extension to more features for increased accuracy. The dataset, released with this publication, may also be used to test other proposed surrogates, and we find that many previously published proxies do correlate well with both the heat flux and stability boundary.

From the near-Earth solar wind to the intracluster medium of galaxy clusters, collisionless, high-beta, magnetized plasmas pervade our universe. Energy and momentum transport from large-scale fields and flows to small-scale motions of plasma particles is ubiquitous in these systems, but a full picture of the underlying physical mechanisms remains elusive. The transfer is often mediated by a turbulent cascade of Alfvénic fluctuations as well as a variety of kinetic instabilities; these processes tend to be multi-scale and/or multi-dimensional, which makes them difficult to study using spacecraft missions and numerical simulations alone. Meanwhile, existing laboratory devices struggle to produce the collisionless, high ion beta ($\beta _i \gtrsim 1$), magnetized plasmas across the range of scales necessary to address these problems. As envisioned in recent community planning documents, it is therefore important to build a next generation laboratory facility to create a $\beta _i \gtrsim 1$, collisionless, magnetized plasma in the laboratory for the first time. A working group has been formed and is actively defining the necessary technical requirements to move the facility towards a construction-ready state. Recent progress includes the development of target parameters and diagnostic requirements as well as the identification of a need for source-target device geometry. As the working group is already leading to new synergies across the community, we anticipate a broad community of users funded by a variety of federal agencies (including National Aeronautics and Space Administration, Department of Energy and National Science Foundation) to make copious use of the future facility.

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.

The generation of an autoresonantly phase-locked high-amplitude plasma waves to the chirped beat frequency of two driving lasers is studied in two dimensions using particle-in-cell simulations. The two-dimensional plasma and laser parameters correspond to those that optimized the plasma wave amplitude in one-dimensional simulations. Near the start of autoresonant locking, the two-dimensional simulations appear similar to one-dimensional particle-in-cell results (Luo et al., Phys. Rev. Res., vol. 6, 2024, p. 013338) with plasma wave amplitudes above the Rosenbluth–Liu limit. Later, just below wave breaking, the two-dimensional simulation exhibits a Weibel-like instability and eventually laser beam filamentation. These limit the coherence of the plasma oscillation after the peak plasma wave field is obtained. In spite of the reduction of spatial coherence of the accelerating density structure, the acceleration of self-injected electrons in the case studied remains at $70\,\%$ to $80\,\%$ of that observed in one dimension. Other effects such as plasma wave bowing are discussed.

Near-future experiments with Petawatt class lasers are expected to produce a high flux of gamma-ray photons and electron–positron pairs through strong field quantum electrodynamical processes. Simulations of the expected regime of laser–matter interaction are computationally intensive due to the disparity of the spatial and temporal scales, and because quantum and classical descriptions need to be accounted for simultaneously (classical for collective effects and quantum for nearly instantaneous events of hard photon emission and pair creation). We study the stochastic cooling of an electron beam in a strong, constant, uniform magnetic field, both its particle distribution functions and their energy momenta. We start by obtaining approximate closed-form analytical solutions to the relevant observables. Then, we apply the quantum-hybrid variational quantum imaginary time evolution to the Fokker–Planck equation describing this process and compare it against theory and results from particle-in-cell simulations and classical partial differential equation solvers, showing good agreement. This work will be useful as a first step towards quantum simulation of plasma physics scenarios where diffusion processes are important, particularly in strong electromagnetic fields.

Bounce-averaged theories provide a framework for simulating relatively slow processes, such as collisional transport and quasilinear diffusion, by averaging these processes over the fast periodic motions of a particle on a closed orbit. This procedure dramatically increases the characteristic time scale and reduces the dimensionality of the modelled system. The natural coordinates for such calculations are the constants of motion (COM) of the fast particle motion, which by definition do not change during an orbit. However, for sufficiently complicated fields – particularly in the presence of local maxima of the electric potential and magnetic field – the COM are not sufficient to specify the particle trajectory. In such cases, multiple domains in COM space must be used to solve the problem, with boundary conditions enforced between the domains to ensure continuity and particle conservation. Previously, these domains have been imposed by hand, or by recognising local maxima in the fields, limiting the flexibility of bounce-averaged simulations. Here, we present a general set of conditions for identifying consistent domains and the boundary condition connections between the domains, allowing the application of bounce-averaged theories in arbitrarily complicated and dynamically evolving electromagnetic field geometries. We also show how the connections between the domains can be represented by a directed graph, which can help to succinctly represent the trajectory bifurcation structure.

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 driven by gyrokinetic instabilities is largely responsible for transport in magnetic fusion devices. To estimate this turbulent transport, integrated modelling codes often use mixing length estimates in conjunction with reduced models of the linearized gyrokinetic equation. One common method of formulating and solving the linearized gyrokinetic eigenvalue problem equation uses a Ritz variational principle, particularly in the local collisionless limit. However, the variational principle as typically stated in the literature is mathematically incorrect. In this work, we derive a mathematically correct form of the variational principle that applies to local linear collisionless gyrokinetics in general geometry with electromagnetic effects. We also explicitly derive a weak form of the gyrokinetic field equations suitable for numerical applications.

We solve ‘half’ the problem of finding three-dimensional quasisymmetric magnetic fields that do not necessarily satisfy magnetohydrostatic force balance. This involves determining which hidden symmetries are admissible as quasisymmetries, and then showing explicitly how to construct quasisymmetric magnetic fields given an admissible symmetry. The admissibility conditions take the form of a system of overdetermined nonlinear partial differential equations involving second derivatives of the symmetry's infinitesimal generator.

In one of the leading theories for the origin of the solar wind, photospheric motions launch Alfvén waves (AWs) that propagate along open magnetic-field lines through the solar atmosphere and into the solar wind. The radial variation in the Alfvén speed causes some of the AWs to reflect, and counter-propagating AWs subsequently interact to produce Alfveńic turbulence, in which AW energy cascades from long wavelengths to short wavelengths and dissipates, heating the plasma. In this paper we develop a one-dimensional two-fluid solar-wind model that includes Alfvénic turbulence, proton temperature anisotropy and a novel method for apportioning the turbulent heating rate between parallel proton heating, perpendicular proton heating and electron heating. We employ a turbulence model that accounts for recent observations from NASA’s Parker Solar Probe, which find that AW fluctuations in the near-Sun solar wind are intermittent and less anisotropic than in previous models of anisotropic magnetohydrodynamic turbulence. Our solar-wind model reproduces a wide range of remote observations of the corona and in-situ measurements of the solar wind, and our turbulent heating model consists of analytic equations that could be usefully incorporated into other solar-wind models and numerical models of more distant astrophysical plasmas.

The Floquet exponents of periodic field lines are studied through the variations of the magnetic action on the magnetic axis, which is assumed to be elliptical. The near-axis formalism developed by Mercier, Solov'ev and Shafranov is combined with a Lagrangian approach. The on-axis Floquet exponent is shown to coincide with the on-axis rotational transform. A discrete solution suitable for numerical implementation is introduced, which gives the Floquet exponents as solutions to an eigenvalue problem. This discrete formalism expresses the exponents as the eigenvalues of a $6\times 6$ matrix.

Motivated by explosive releases of energy in fusion, space and astrophysical plasmas, we consider the nonlinear stability of stratified magnetohydrodynamic equilibria against two-dimensional interchanges of straight magnetic-flux tubes. We demonstrate that, even within this restricted class of dynamics, the linear stability of an equilibrium does not guarantee its nonlinear stability: equilibria can be metastable. We show that the minimum-energy state accessible to a metastable equilibrium under non-diffusive two-dimensional dynamics can be found by solving a combinatorial optimisation problem. These minimum-energy states are, to good approximation, the final states reached by our simulations of destabilised metastable equilibria for which turbulent mixing is suppressed by viscosity. To predict the result of fully turbulent relaxation, we construct a statistical mechanical theory based on the maximisation of Boltzmann's mixing entropy. This theory is analogous to the Lynden-Bell statistical mechanics of collisionless stellar systems and plasma, and to the Robert–Sommeria–Miller theory of two-dimensional vortex turbulence. Our theory reproduces well the results of our numerical simulations for sufficiently large perturbations to the metastable equilibrium.

We investigate self-consistent, steady-state axisymmetric solutions of an incompressible tokamak plasma using a visco-resistive magnetohydrodynamic model. A key contribution of this work is the formulation of Poisson’s equation that governs the pressure profile. Our analysis reveals that the current modelling fails to produce realistic pressure levels. To overcome this limitation, we introduce additional non-inductive current drives, akin to those generated by neutral beam injection or radio frequency heating, modelled as modifications to the toroidal current. Numerical simulations validate our enhanced model, showing significant improvements in pressure profile characteristics. In the cases examined, the effect of these current drives on the velocity profiles is moderate, except when the non-inductive current drives induce reversals in the total toroidal current density, leading to non-nested flux surfaces with internal separatrices.

The use of gyrokinetics, wherein phase-space coordinate transformations result in a phase-space dimensionality reduction as well as the removal of fast time scales, has enabled the simulation of microturbulence in fusion devices. The state-of-the-art gyrokinetic models used in practice are parallel-only models wherein the perpendicular part of the vector potential is neglected. Such models are inherently not gauge-invariant. We generalise the work of Burby & Brizard (2019 Phys. Lett. A vol. 383, no. 18, pp. 2172–2175) by deriving a sufficient condition on the gyrocentre coordinate transformation that ensures gauge invariance. This leads to a parametrised family of gyrokinetic models for which we motivate a specific choice of parameters that results in the smallest gyrocentre coordinate transformation for which the resulting gyrokinetic model is consistent, gyro-phase independent, gauge-invariant and has an invariant magnetic moment. Due to gauge invariance, this model can be expressed directly in terms of the electromagnetic fields rather than the potentials, and the gyrokinetic model thereby results in the macroscopic Maxwell’s equations. For the linearised model, it is demonstrated that the shear and compressional Alfvén waves are present with the correct frequencies. The fast compressional Alfvén wave can be removed by making use of a Darwin approximation. This approximation retains the gauge invariance of the proposed model.