A fusion neutron source (FNS) based on the gas-dynamic trap (GDT, Budker Institute, Novosibirsk) is considered for confinement of two-species plasma heated by neutral beam injection in a regime where the fast ion distribution function is far from Maxwellian. Kinetic instabilities are expected to develop in this regime, and in this paper we investigate the ion-cyclotron instability evolving in moderate densities of pure hydrogen and mixed deuterium–hydrogen target plasmas. The properties of the studied unstable mode, such as its azimuthal wavenumbers, propagation direction and its being affected by changes in the bulk plasma density and composition, allow us to identify it as the drift cyclotron loss cone (DCLC) instability. This mode scatters fast ions and thereby leads to drops in diamagnetic flux signals and increases longitudinal energy and particle losses, with the average energy of the lost ions estimated to be far above the temperature of warm Maxwellian ions. Our interpretation is that the unstable wave grows due to interaction with the fast ions located near the loss cone in the velocity space and scatters them. Applying the method of suppressing the DCLC instability by filling the loss cone with warm plasma, we have determined the values of plasma density and deuterium percentage that allow us to suppress the DCLC instability in the GDT. These findings justify using mixed bulk plasmas in fusion neutron source operation.

The propagation of electromagnetic waves in a linearly varying index of refraction is a fundamental problem in wave physics, being relevant in fusion science for describing certain wave-based heating and diagnostic schemes. Here, an exact solution is obtained for a given incoming wavefield specified on the boundary transverse to the direction of inhomogeneity by performing a spectral, rather than asymptotic, matching. Two case studies are then presented: a Gaussian beam at oblique incidence and a speckled wavefield at normal incidence. For the Gaussian beam, it is shown that when the waist $W$ is sufficiently large, oblique incidence manifests simply as rigid translation and focal shift of the corresponding diffraction pattern at normal incidence. The destruction of the hyperbolic umbilic caustic (corresponding to a critically focused beam) as $W$ is reduced is then demonstrated. The caustic disappears once $W \lesssim \delta _a \sqrt {L}$ ($L$ being the medium length scale normalized by the Airy skin depth $\delta _a$), at which point the wave behaviour is increasingly described by Airy functions, but experiences less focusing as a result. To maximize the intensity of a launched Gaussian beam at a turning point, one should therefore minimize the imaginary part of the launched complex beam parameter while having the real part satisfy critical focusing. For the speckled wavefield, it is shown that the transverse speckle pattern only couples to the Airy longitudinal pattern when the coupling parameter $\eta = \sqrt {L}/f_{\#}$ is large, with $f_{\#}$ being the f-number of the launching aperture. When $\eta \ll 1$, a reduced description of the total wavefield can be obtained by simply multiplying the incoming speckle pattern with the Airy swelling.

In the radiation hydrodynamic simulations used to design inertial confinement fusion (ICF) and pulsed power experiments, nonlinear radiation diffusion tends to dominate CPU time. This raises the interesting question of whether a quantum algorithm can be found for nonlinear radiation diffusion which provides a quantum speedup. Recently, such a quantum algorithm was introduced based on a quantum algorithm for solving systems of nonlinear partial differential equations (PDEs) which provides a quadratic quantum speedup. Here, we apply this quantum PDE (QPDE) algorithm to the problem of a non-equilibrium Marshak wave propagating through a cold, semi-infinite, optically thick target, where the radiation and matter fields are not assumed to be in local thermodynamic equilibrium. The dynamics is governed by a coupled pair of nonlinear PDEs which are solved using the QPDE algorithm, as well as two standard PDE solvers: (i) Python's py-pde solver; and (ii) the KULL ICF simulation code developed at Lawrence-Livermore National Laboratory. We compare the simulation results obtained using the QPDE algorithm and the standard PDE solvers and find excellent agreement.

In magnetic confinement fusion devices close to axisymmetry, such as tokamaks, a key element is the winding profile of the magnetic field lines, or its inverse, the safety profile $q=q_{\boldsymbol {B}}$. A corresponding profile, $q_{\boldsymbol {J}}$, can be defined for the current density field lines. Ampère's law relates any mode of current perturbation $\delta \boldsymbol {J}_{m,n}$ with a mode of magnetic perturbation $\delta \boldsymbol {B}_{m,n}$. It is shown that the knowledge of the pair $(q_{\boldsymbol {B}},q_{\boldsymbol {J}})$ allows us then to characterize the resonant, or non-resonant, nature of the modes for both the magnetic and current density field lines. The expression of $q_{\boldsymbol {J}}$ in the flux coordinate is derived. Including this calculation in real-time Grad–Shafranov equilibrium reconstruction codes would yield a comprehensive view of the magnetics. The monitoring of the pair $(q_{\boldsymbol {B}},q_{\boldsymbol {J}})$ would then allow us to investigate the role played by the resonant modes for the current density, that are current filamentary modes, in the plasma small-scale turbulence. By driving the magnetic and current density profiles apart so that the images of $q_{\boldsymbol {B}}$ and $q_{\boldsymbol {J}}$ are disjoint, these filamentary modes would not impact the magnetic field topology, being not associated with magnetic islands but with non-resonant magnetic modes. It remains to be explored to what extent such a configuration, where the spectrum of tiny current density filaments produces a spectrum of magnetic modes that has practically no effect on heat transport, is beneficial.

Free normal-flow (NF) conditions at the plasma boundary are shown to be essential for three-dimensional magnetohydrodynamic (MHD) simulations to agree with linear stability theory. A comparative verification study is presented between two different formulations of the boundary conditions (BCs) for velocity perturbations: (i) fully consistent free NF implementation and (ii) fixed NF formulation, neglecting flow perturbations at the numerical boundary. Numerical results are compared with consolidated figures of merit from the linear theory of external kink modes. We consider two classes of initial equilibria presenting different numerical challenges: a uniform current channel surrounded by pure vacuum and a shaped Wesson-like equilibrium, with smooth (polynomial) radial dependency. Only the fully consistent free NF formulation is invariably accurate in modelling the plasma interface at the numerical boundary, without the need of enforcing a pseudovacuum region at the edge of the simulation domain, as in most analogous past studies. This study employs the cylindrical code SPECYL (Cappello & Biskamp, Nucl. Fusion, vol. 36, no. 5, 1996, p. 571) that solves a full-MHD model without pressure gradients, whose fully consistent resistive wall module with free NF BCs was recently successfully verified against the independent code PIXIE3D (Spinicci et al., AIP Adv., vol. 13, no. 9, 2023, p. 095111).

We derive an advanced surrogate model for predicting turbulent transport at the edge of tokamaks driven by electron temperature gradient (ETG) modes. Our derivation is based on a recently developed sensitivity-driven sparse grid interpolation approach for uncertainty quantification and sensitivity analysis at scale, which informs the set of parameters that define the surrogate model as a scaling law. Our model reveals that ETG-driven electron heat flux is influenced by the safety factor $q$, electron beta $\beta _e$ and normalized electron Debye length $\lambda _D$, in addition to well-established parameters such as the electron temperature and density gradients. To assess the trustworthiness of our model's predictions beyond training, we compute prediction intervals using bootstrapping. The surrogate model's predictive power is tested across a wide range of parameter values, including within-distribution testing parameters (to verify our model) as well as out-of-bounds and out-of-distribution testing (to validate the proposed model). Overall, validation efforts show that our model competes well with, or can even outperform, existing scaling laws in predicting ETG-driven transport.

Quasilinear treatments are widely used for tokamaks to evaluate radio frequency (rf) heating and current drive. Even though the core of a tokamak plasma is weakly collisional, the solution of the linearized kinetic equation is evaluated using unperturbed collisionless trajectories while often treating successive poloidal circuits of the passing (and trapped) particles as uncorrelated or nearly so. In addition, the most important effect of tokamak geometry, the mirror force, is usually mistreated or ignored when obtaining the solution. These concerning aspects of rf treatments are clarified by considering lower hybrid heating and current drive to illustrate that the electrons in resonance with the applied rf are enclosed by narrow collisional boundary layers, and that tokamak geometry makes it necessary to retain poloidal variation when solving a weakly collisional linearized kinetic equation. Other aspects such as collisional boundary layers at the trapped–passing boundary, cyclotron resonances, and the limitations of quasilinear theory are also considered. The new insights lead to a fundamentally different formulation and interpretation of the solution of the linearized Fokker–Planck equation used for rf quasilinear theory in a tokamak, while retaining many of the features that have contributed to its successful application to rf heating and current drive.

Hot electrons produced by parametric instabilities can negatively affect the thermonuclear gain of a direct-drive inertial confinement fusion capsule. A Monte Carlo transport is coupled to the hydrodynamics MULTI-IFE code in order to study the hot electron transport in a high-gain target implosion. The thermonuclear energy produced by the directly driven implosion of a reference target drops from 27 MJ without hot electrons to 13.5 kJ when the hot electrons are taken into account. A clear relationship is established between the level of hot electrons produced and the degradation of the thermonuclear yield produced by the target. It is shown that adding a relatively short laser spike (shock ignition) to the main pulse restores the thermonuclear energy released by the implosion. Different shock-ignition windows are presented depending on the duration of the spike. It appears that the longer the duration of the spike, the lower the power required for shock ignition. Finally, a 300 ps-duration spike requires only 50 TW to restore a thermonuclear energy output of 20 MJ.

The paper considers the electron transport in toroidal systems taking into account relativistic effects for electrons. The treatment is based on the relativistic drift-kinetic equation with the thermodynamic equilibrium given by the relativistic Maxwell–Jüttner distribution function. The definition of relativistic fluxes is given in a classic-like form using the same set of thermodynamic forces as in the classical (non-relativistic) approach. Such a formulation allows us to apply the currently used non-relativistic solvers for calculation of relativistic mono-energetic transport coefficients. As an example, the procedure for calculating electron fluxes is proposed, in which relativistic effects are taken into account using the DKES code. The model can be easily implemented in various transport codes, developed for the non-relativistic limit, making them accurate also for hot plasmas with non-negligible relativistic effects.

In this work we consider the problem of optimizing a stellarator subject to hard constraints on the design variables and physics properties of the equilibrium. We survey current numerical methods for handling these constraints, and summarize a number of methods from the wider optimization community that have not been used extensively for stellarator optimization thus far. We demonstrate the utility of new methods of constrained optimization by optimizing a quasi-axisymmetric stellarator for favourable physics properties while preventing strong shaping of the plasma boundary, which can be difficult to create with external current sources.

We have derived a necessary condition for the achievement of breakeven in axisymmetric plasmas with zero toroidal field and confined by a dipole magnetic field (B. Lehnert, Nature, vol. 181, 1958, p. 4605; A. Hasegawa, Comments Plasma Phys. Control. Fusion, vol. 11, no. 3, 1987). Excellent MHD stability, high values of $\beta$ (up to $26$ %) and good confinement properties awaken the interest of private investors after years of neglect due to lack of public funding and competing alternative lines of research like the tokamak. Starting from a requirement of self-consistency between the balances of momentum and energy in a dipole-confined, two-species plasma and assuming a Maxwellian distribution function for ions and electrons, we derive a necessary condition for breakeven. This condition is more stringent than the Lawson criterion because of the lack of a stabilizing toroidal field. For a given current flowing across the toroidal coil internal to the plasma, the crucial factor at stake is the ratio between the radius of the main toroidal coil and the radius of the vacuum chamber.

Electrostatic gyrokinetic instabilities and turbulence in the Wendelstein 7-X stellarator are studied. Particular attention is paid to the ion-temperature-gradient (ITG) instability and its character close to marginal stability [Floquet-type turbulence (Zocco et al., Phys. Rev. E, vol. 106, 2022, p. L013202) with no electron temperature gradient]. The flux tube version of the $\delta f$ code stella (Barnes et al., J. Comput. Phys., vol. 391, 2019, pp. 365–380) is used to run linear and nonlinear gyrokinetic simulations with kinetic electrons. The nature of the dominant instability depends on the wavelength perpendicular to the magnetic field, and the results are conveniently displayed in stability diagrams that take this dependence into account. This approach highlights the presence of universal instabilities, which are less unstable but have longer wavelengths than other modes. A quasi-linear estimate of the heat flux suggests they are relevant for transport. Close to the stability threshold, the linear eigenmodes and turbulence form highly extended structures along the computational domain if the magnetic shear is small. Numerical experiments and diagnostics are undertaken to assess the resulting radial localisation of the turbulence, which affects the interaction of the latter with zonal flows. Increasing the amplitude of the magnetic shear (e.g.through current drive) has a stabilising effect on the turbulence and, thus, reduces the nonlinear energy transport.

A surrogate model of the runaway electron avalanche growth rate in a magnetic fusion plasma is developed. This is accomplished by employing a physics-informed neural network (PINN) to learn the parametric solution of the adjoint to the relativistic Fokker–Planck equation. The resulting PINN is able to evaluate the runaway probability function across a broad range of parameters in the absence of any synthetic or experimental data. This surrogate of the adjoint relativistic Fokker–Planck equation is then used to infer the avalanche growth rate as a function of the electric field, synchrotron radiation and effective charge. Predictions of the avalanche PINN are compared against first principle calculations of the avalanche growth rate with excellent agreement observed across a broad range of parameters.

Microinstabilities drive turbulent fluctuations in inhomogeneous, magnetised plasmas. In the context of magnetic confinement fusion devices, this leads to an enhanced transport of particles, momentum and energy, thereby degrading confinement. In this work, we describe an application of the adjoint method to efficiently determine variations of gyrokinetic linear growth rates on a general set of external parameters in the local $\delta f$-gyrokinetic model. We then offer numerical verification of this approach. When coupled with gradient-based techniques, this methodology can facilitate the optimisation process for the microstability of the confined plasmas across a high-dimensional parameter space. We present a numerical demonstration wherein the ion-temperature-gradient instability growth rate in a tokamak plasma is minimised with respect to flux surface shaping parameters. The adjoint method approach demonstrates a significant computational speed-up compared with a finite-difference gradient calculation.

GX is a code designed to solve the nonlinear gyrokinetic system for low-frequency turbulence in magnetized plasmas, particularly tokamaks and stellarators. In GX, our primary motivation and target is a fast gyrokinetic solver that can be used for fusion reactor design and optimization along with wide-ranging physics exploration. This has led to several code and algorithm design decisions, specifically chosen to prioritize time to solution. First, we have used a discretization algorithm that is pseudospectral in the entire phase space, including a Laguerre–Hermite pseudospectral formulation of velocity space, which allows for smooth interpolation between coarse gyrofluid-like resolutions and finer conventional gyrokinetic resolutions and efficient evaluation of a model collision operator. Additionally, we have built GX to natively target graphics processors (GPUs), which are among the fastest computational platforms available today. Finally, we have taken advantage of the reactor-relevant limit of small $\rho _*$ by using the radially local flux-tube approach. In this paper we present details about the gyrokinetic system and the numerical algorithms used in GX to solve the system. We then present several numerical benchmarks against established gyrokinetic codes in both tokamak and stellarator magnetic geometries to verify that GX correctly simulates gyrokinetic turbulence in the small $\rho _*$ limit. Moreover, we show that the convergence properties of the Laguerre–Hermite spectral velocity formulation are quite favourable for nonlinear problems of interest. Coupled with GPU acceleration, which we also investigate with scaling studies, this enables GX to be able to produce useful turbulence simulations in minutes on one (or a few) GPUs and higher fidelity results in a few hours using several GPUs. GX is open-source software that is ready for fusion reactor design studies.

A systematic theory of the asymptotic expansion of the magnetohydrostatics (MHS) equilibrium in the distance from the magnetic axis is developed to include arbitrary smooth currents near the magnetic axis. Compared with the vacuum and the force-free system, an additional magnetic differential equation must be solved to obtain the pressure-driven currents. It is shown that there exist variables in which the rest of the MHS system closely mimics the vacuum system. Thus, a unified treatment of MHS fields is possible. The mathematical structure of the near-axis expansions to arbitrary order is examined carefully to show that the double-periodicity of physical quantities in a toroidal domain can be satisfied order by order. The essential role played by the leading-order Birkhoff–Gustavson normal form in solving the magnetic differential equations is highlighted. Several explicit examples of vacuum, force-free and MHS equilibrium in different geometries are presented.

Filament-based coil optimizations are performed for several quasi-helical stellarator configurations, beginning with the one from Landreman & Paul (Phys. Rev. Lett., vol. 128, 2022, 035001), demonstrating that precise quasi-helical symmetry can be achieved with realistic coils. Several constraints are placed on the shape and spacing of the coils, such as low curvature and sufficient plasma–coil distance for neutron shielding. The coils resulting from this optimization have a maximum curvature 0.8 times that of the coils of the Helically Symmetric eXperiment (HSX) and a mean squared curvature 0.4 times that of the HSX coils when scaled to the same plasma minor radius. When scaled up to reactor size and magnetic field strength, no fast particle losses were found in the free-boundary configuration when simulating 5000 alpha particles launched at $3.5\,\mathrm {MeV}$ on the flux surface with a normalized toroidal flux of $s=0.5$. An analysis of the tolerance of the coils to manufacturing errors is performed using a Gaussian process model, and the coils are found to maintain low particle losses for smooth, large-scale errors up to amplitudes of approximately $0.15\,\mathrm {m}$. Another coil optimization is performed for the Landreman–Paul configuration with the additional constraint that the coils are purely planar. Visual inspection of the Poincaré plot of the resulting magnetic field-lines reveal that the planar modular coils alone do a poor job of reproducing the target equilibrium. Additional non-planar coil optimizations are performed for the quasi-helical configuration with $5\,\%$ volume-averaged plasma beta from Landreman et al. (Phys. Plasma, vol. 29, issue 8, 2022, 082501), and a similar configuration also optimized to satisfy the Mercier criterion. The finite beta configurations had larger fast-particle losses, with the free-boundary Mercier-optimized configuration performing the worst, losing approximately $5.5\,\%$ of alpha particles launched at $s=0.5$.

When simulating runaway electron dynamics in tokamak disruptions, fluid models with lower numerical cost are often preferred to more accurate kinetic models. The aim of this work is to compare fluid and kinetic simulations of a large variety of different disruption scenarios in ITER. We consider both non-activated and activated scenarios; for the latter, we derive and implement kinetic sources for the Compton scattering and tritium beta decay runaway electron generation mechanisms in our simulation tool Dream (Hoppe et al., Comput. Phys. Commun., vol. 268, 2021, 108098). To achieve a diverse set of disruption scenarios, Bayesian optimization is used to explore a range of massive material injection densities for deuterium and neon. The cost function is designed to distinguish between successful and unsuccessful disruption mitigation based on the runaway current, current quench time and transported fraction of the heat loss. In the non-activated scenarios, we find that fluid and kinetic disruption simulations can have significantly different runaway electron dynamics, due to an overestimation of the runaway seed by the fluid model. The primary cause of this is that the fluid hot-tail generation model neglects superthermal electron transport losses during the thermal quench. In the activated scenarios, the fluid and kinetic models give similar predictions, which can be explained by the significant influence of the activated sources on the runaway dynamics and the seed.

The problem of the assimilation of a cryogenic fuel pellet injected into a hot plasma is considered. Due to the transparency to ambient particles of the plasmoid, the localised region of high-density plasma created by ionisation of the ablated pellet material, electrons reach a ‘quasiequilibrium’ (QE) state which is characterised by a steady-state on the fastest collisional time scale. The simplified electron kinetic equation of the QE state is solved. Taking a velocity moment of the higher-order electron kinetic equation, which is valid on the expansion time scale, permits a fluid closure, yielding an evolution equation for the macroscopic parameters describing the QE distribution function. In contrast to the Braginskii equations, the closure does not require that electrons have a short mean free path compared with the size of density perturbations, and permits an anisotropic and highly non-Maxwellian distribution function. As the QE distribution function accounts for both trapped and passing electrons, the self-consistent electric potential that causes the expansion can be properly described, in contrast to earlier models of pellet plasmoid expansion with an unbounded potential. The plasmoid expansion is simulated using both a Vlasov model and a cold-fluid model for the ions. During the expansion plasmoid ions and electrons obtain nearly equal amounts of energy; as hot ambient electrons provide this energy in the form of collisional heating of plasmoid electrons, the expansion of a pellet plasmoid is expected to be a potent mechanism for the transfer of energy from electrons to ions on a time scale shorter than that of ion–electron thermalisation.

We present a stellarator configuration optimized for a large threshold (‘critical gradient’) for the onset of the ion temperature gradient (ITG) driven mode, which achieves the largest critical gradient we have seen in any stellarator. Above this threshold, gyrokinetic simulations show that the configuration has low turbulence levels over an experimentally relevant range of the drive strength. The applied optimization seeks to maximize the drift curvature, leading to enhanced local-shear stabilization of toroidal ITG modes, and the associated turbulence. These benefits are combined with excellent quasi-symmetry, yielding low neoclassical transport and vanishingly small alpha particle losses. Analysis of the resulting configuration suggests a trade-off between magnetohydrodynamic (MHD) and ITG stability, as the new configuration possesses a vacuum magnetic hill.