The confinement of energetic particles in nuclear fusion devices is studied in the presence of an oscillating radial electric field and an axisymmetric magnetic equilibrium. It is shown that, despite the poloidal and toroidal symmetries, initially integrable orbits turn into chaotic regions that can potentially intercept the wall of the tokamak, leading to particle losses. It is observed that the losses exhibit algebraic time decay different from the expected exponential decay characteristic of radial diffusive transport. A dynamical explanation of this behaviour is presented, within the continuous time random walk theory. The central point of the analysis is based on the fact that, contrary to the radial displacement, the poloidal angle is not bounded and a proper statistical analysis can therefore be made, showing for the first time that energetic particle transport can be super-diffusive in the poloidal direction and characterised by asymmetric poloidal displacement. The connection between poloidal and radial positions ensured by the conservation of the toroidal canonical momentum, implies that energetic particles spend statistically more time in the inner region of the tokamak than in the outer one, which explains the observed algebraic decay. This indicates that energetic particles might be efficiently slowed down by the thermal population before leaving the system. Also, the asymmetric transport reveals a new possible mechanism of self-generation of momentum.

Understanding particle drifts in a non-symmetric magnetic field is of primary interest in designing optimized stellarators in order to minimize the neoclassical radial loss of particles. Quasisymmetry and omnigeneity, two distinct properties proposed to ensure radial localization of collisionless trapped particles in stellarators, have been explored almost exclusively for magnetic fields with nested flux surfaces. In this work, we examine radial particle confinement when all field lines are closed. We then study charged particle dynamics in the special case of a non-symmetric vacuum magnetic field with closed field lines obtained recently by Weitzner & Sengupta (Phys. Plasmas, vol. 27, 2020, 022509). These magnetic fields can be used to construct magnetohydrodynamic equilibria for low pressure. Expanding in the amplitude of the non-symmetric fields, we explicitly evaluate the omnigeneity and quasisymmetry constraints. We show that the magnetic field is omnigeneous in the sense that the drift surfaces coincide with the pressure surfaces. However, it is not quasisymmetric according to the standard definitions.

Following up on earlier work which demonstrated an improved numerical stellarator coil design optimization performance by the use of stochastic optimization (Lobsien et al., Nucl. Fusion, vol. 58 (10), 2018, 106013), it is demonstrated here that significant further improvements can be made – lower field errors and improved robustness – for a Wendelstein 7-X test case. This is done by increasing the sample size and applying fully three-dimensional perturbations, but most importantly, by changing the design sequence in which the optimization targets are applied: optimization for field error is conducted first, with coil shape penalties only added to the objective function at a later step in the design process. A robust, feasible coil configuration with a local maximum field error of 3.66 % and an average field error of 0.95 % is achieved here, as compared to a maximum local field error of 6.08 % and average field error of 1.56 % found in our earlier work. These new results are compared to those found without stochastic optimization using the FOCUS and ONSET suites. The relationship between local minima in the optimization space and coil shape penalties is also discussed.

Accelerated statistical computation of collisionless fusion alpha particle losses in stellarator configurations is presented based on direct guiding-centre orbit tracing. The approach relies on the combination of recently developed symplectic integrators in canonicalized magnetic flux coordinates and early classification into regular and chaotic orbit types. Only chaotic orbits have to be traced up to the end, as their behaviour is unpredictable. An implementation of this technique is provided in the code SIMPLE (symplectic integration methods for particle loss estimation, Albert et al., 2020b, doi:10.5281/zenodo.3666820). Reliable results were obtained for an ensemble of 1000 orbits in a quasi-isodynamic, a quasi-helical and a quasi-axisymmetric configuration. Overall, a computational speed up of approximately one order of magnitude is achieved compared to direct integration via adaptive Runge–Kutta methods. This reduces run times to the range of typical magnetic equilibrium computations and makes direct alpha particle loss computation adequate for use within a stellarator optimization loop.

In this work, we highlight an issue that may reduce the accuracy of many local nonlinear gyrokinetic simulations – turbulent self-interaction through the parallel boundary condition. Given a sufficiently long parallel correlation length, individual turbulent eddies can span the full domain and ‘bite their own tails’, thereby altering their statistical properties. Such self-interaction is only modelled accurately when the simulation domain corresponds to a full flux surface, otherwise it is artificially strong. For Cyclone Base Case parameters and typical domain sizes, we find that this mechanism modifies the heat flux by approximately 40 % and it can be even more important. The effect is largest when using kinetic electrons, low magnetic shear and strong turbulence drive (i.e. steep background gradients). It is found that parallel self-interaction can be eliminated by increasing the parallel length and/or the binormal width of the simulation domain until convergence is achieved.

A gyrokinetic model is presented that can properly describe large and small amplitude electromagnetic fluctuations occurring on scale lengths ranging from the electron Larmor radius to the equilibrium perpendicular pressure gradient scale length, and the arbitrarily large deviations from thermal equilibrium that are present in the plasma periphery of tokamak devices. The formulation of the gyrokinetic model is based on a second-order accurate description of the single charged particle dynamics, derived from Lie perturbation theory, where the fast particle gyromotion is decoupled from the slow drifts assuming that the ratio of the ion sound Larmor radius to the perpendicular equilibrium pressure scale length is small. The collective behaviour of the plasma is obtained by a gyrokinetic Boltzmann equation that describes the evolution of the gyroaveraged distribution function. The collisional effects are included by a nonlinear gyrokinetic Dougherty collision operator. The gyrokinetic model is then developed into a set of coupled fluid equations referred to as the gyrokinetic moment hierarchy. To obtain this hierarchy, the gyroaveraged distribution function is expanded onto a Hermite–Laguerre velocity-space polynomial basis. Then, the gyrokinetic equation is projected onto the same basis yielding the spatial and temporal evolution of the Hermite–Laguerre expansion coefficients. A closed set of fluid equations for the lowest-order coefficients is presented. The Hermite–Laguerre projection is performed accurately at arbitrary perpendicular wavenumber values. Finally, the self-consistent evolution of the electromagnetic fields is described by a set of gyrokinetic Maxwell equations derived from a variational principle where the velocity integrals are explicitly evaluated.

The dependence of confinement on input power for a tokamak plasma with regions having a stiff temperature profile is explored. The resilience of the confinement of the core energy to increasing power loss by core radiation from impurities in such situations, as it is anticipated will be required in a demonstration fusion reactor (DEMO) design, is examined.

The concept of the available energy of a collisionless plasma is discussed in the context of magnetic confinement. The available energy quantifies how much of the plasma energy can be converted into fluctuations (including nonlinear ones) and is thus a measure of plasma stability, which can be used to derive linear and nonlinear stability criteria without solving an eigenvalue problem. In a magnetically confined plasma, the available energy is determined by the density and temperature profiles as well as the magnetic geometry. It also depends on what constraints limit the possible forms of plasma motion, such as the conservation of adiabatic invariants and the requirement that the transport be ambipolar. A general method based on Lagrange multipliers is devised to incorporate such constraints in the calculation of the available energy, and several particular cases are discussed for which it can be calculated explicitly. In particular, it is shown that it is impossible to confine a plasma in a Maxwellian ground state relative to perturbations with frequencies exceeding the ion bounce frequency.

We present an energy-conserving discontinuous Galerkin scheme for the full-$f$ electromagnetic gyrokinetic system in the long-wavelength limit. We use the symplectic formulation and solve directly for $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$, the inductive component of the parallel electric field, using a generalized Ohm’s law derived directly from the gyrokinetic equation. Linear benchmarks are performed to verify the implementation and show that the scheme avoids the Ampère cancellation problem. We perform a nonlinear electromagnetic simulation in a helical open-field-line system as a rough model of the tokamak scrape-off layer using parameters from the National Spherical Torus Experiment (NSTX). This is the first published nonlinear electromagnetic gyrokinetic simulation on open field lines. Comparisons are made to a corresponding electrostatic simulation.

Mechanisms setting the density decay in the scrape-off layer (SOL) at the outer midplane of a tokamak plasma are disentangled using two-fluid numerical simulations in a double-null magnetic configuration and analytical estimates. Typical experimental observations are retrieved, in particular increasing intermittency of the turbulence going from the near to the far SOL, which is reflected in two different density decay lengths. The decay length of the near SOL is well described as the result of transport driven by a nonlinearly saturated ballooning instability, while in the far SOL, the density decay length is described using a model of intermittent transport mediated by blobs. The analytical estimates of the decay lengths agree well with the simulation results and typical experimental values and can therefore be used to guide tokamak design and operation.

A direct construction of equilibrium magnetic fields with toroidal topology at arbitrary order in the distance from the magnetic axis is carried out, yielding an analytical framework able to explore the landscape of possible magnetic flux surfaces in the vicinity of the axis. This framework can provide meaningful analytical insight into the character of high-aspect-ratio stellarator shapes, such as the dependence of the rotational transform and the plasma beta limit on geometrical properties of the resulting flux surfaces. The approach developed here is based on an asymptotic expansion on the inverse aspect ratio of the ideal magnetohydrodynamics equation. The analysis is simplified by using an orthogonal coordinate system relative to the Frenet–Serret frame at the magnetic axis. The magnetic field vector, the toroidal magnetic flux, the current density, the field line label and the rotational transform are derived at arbitrary order in the expansion parameter. Moreover, a comparison with a near-axis expansion formalism employing an inverse coordinate method based on Boozer coordinates (the so-called Garren–Boozer construction) is made, where both methods are shown to agree at lowest order. Finally, as a practical example, a numerical solution using a W7-X equilibrium is presented, and a comparison between the lowest-order solution and the W7-X magnetic field is performed.

Plasma terminating disruptions in tokamaks may result in relativistic runaway electron beams with potentially serious consequences for future devices with large plasma currents. In this paper, we investigate the effect of plasma elongation on the coupled dynamics of runaway generation and resistive diffusion of the electric field. We find that elongated plasmas are less likely to produce large runaway currents, partly due to the lower induced electric fields associated with larger plasmas, and partly due to direct shaping effects, which mainly lead to a reduction in the runaway avalanche gain.

The shape gradient is a local sensitivity function defined on the surface of an object which provides the change in a characteristic quantity, or figure of merit, associated with a perturbation to the shape of the object. The shape gradient can be used for gradient-based optimization, sensitivity analysis and tolerance calculations. However, it is generally expensive to compute from finite-difference derivatives for shapes that are described by many parameters, as is the case for typical stellarator geometry. In an accompanying work (Antonsen, Paul & Landreman J. Plasma Phys., vol. 85 (2), 2019), generalized self-adjointness relations are obtained for magnetohydrodynamic (MHD) equilibria. These describe the relation between perturbed equilibria due to changes in the rotational transform or toroidal current profiles, displacements of the plasma boundary, modifications of currents in the vacuum region or the addition of bulk forces. These are applied to efficiently compute the shape gradient of functions of MHD equilibria with an adjoint approach. In this way, the shape derivative with respect to any perturbation applied to the plasma boundary or coil shapes can be computed with only one additional MHD equilibrium solution. We demonstrate that this approach is applicable for several figures of merit of interest for stellarator configuration optimization: the magnetic well, the magnetic ripple on axis, the departure from quasisymmetry, the effective ripple in the low-collisionality $1/\unicode[STIX]{x1D708}$ regime $(\unicode[STIX]{x1D716}_{\text{eff}}^{3/2})$ (Nemov et al. Phys. Plasmas, vol. 6 (12), 1999, pp. 4622–4632) and several finite-collisionality neoclassical quantities. Numerical verification of this method is demonstrated for the magnetic well figure of merit with the VMEC code (Hirshman & Whitson Phys. Fluids, vol. 26 (12), 1983, p. 3553) and for the magnetic ripple with modification of the ANIMEC code (Cooper et al. Comput. Phys. Commun., vol. 72 (1), 1992, pp. 1–13). Comparisons with the direct approach demonstrate that, in order to obtain agreement within several per cent, the adjoint approach provides a factor of $O(10^{3})$ in computational savings.

A novel derivation of the parallel ion velocity, and the bootstrap and Pfirsch–Schlüter currents in an imperfectly optimized (that is, almost omnigenous) stellarator magnetic field, $\boldsymbol{B}$, is presented that somewhat more generally recovers expressions completely consistent with previous analytic results. However, it is also shown that, when the conventional radially local form of the drift kinetic equation is employed, the flow velocity and the bootstrap current acquire a spurious contribution proportional to $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D708}$, where $\unicode[STIX]{x1D714}$ denotes the $\boldsymbol{E}\times \boldsymbol{B}$ rotation frequency (due to the radial electric field $\boldsymbol{E}$) and $\unicode[STIX]{x1D708}$ the collision frequency. This contribution is particularly large in the $\sqrt{\unicode[STIX]{x1D708}}$ regime and at smaller collisionalities, where $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D708}\gtrsim 1$, and is presumably present in most numerical calculations, but it disappears if a more accurate drift kinetic equation is used.

A non-parametric Gaussian process regression model is developed in the three-dimensional equilibrium reconstruction code V3FIT. A Gaussian process is a normal distribution of functions that is uniquely defined by specifying a mean function and covariance kernel function. Gaussian process regression assumes that an unknown profile belongs to a particular Gaussian process and uses Bayesian analysis to select the function the give the best fit to measured data. The implementation in V3FIT uses a hybrid representation where Gaussian processes are used to infer some of the equilibrium profiles and standard parametric techniques are used to infer the remaining profiles. The implementation of the Gaussian process is tested using both synthetic data and experimental data from multiple machines.

Integrated modelling of electron runaway requires computationally expensive kinetic models that are self-consistently coupled to the evolution of the background plasma parameters. The computational expense can be reduced by using parameterized runaway generation rates rather than solving the full kinetic problem. However, currently available generation rates neglect several important effects; in particular, they are not valid in the presence of partially ionized impurities. In this work, we construct a multilayer neural network for the Dreicer runaway generation rate which is trained on data obtained from kinetic simulations performed for a wide range of plasma parameters and impurities. The neural network accurately reproduces the Dreicer runaway generation rate obtained by the kinetic solver. By implementing it in a fluid runaway-electron modelling tool, we show that the improved generation rates lead to significant differences in the self-consistent runaway dynamics as compared to the results using the previously available formulas for the runaway generation rate.

A method is given to rapidly compute quasisymmetric stellarator magnetic fields for plasma confinement, without the need to call a three-dimensional magnetohydrodynamic equilibrium code inside an optimization iteration. The method is based on direct solution of the equations of magnetohydrodynamic equilibrium and quasisymmetry using Garren & Boozer’s expansion about the magnetic axis (Phys Fluids B, vol. 3, 1991, p. 2805), and it is several orders of magnitude faster than the conventional optimization approach. The work here extends the method of Landreman et al. (J. Plasma Phys., vol. 85, 2019, 905850103), which was limited to flux surfaces with elliptical cross-section, to higher order in the aspect-ratio expansion. As a result, configurations can be generated with strong shaping that achieve quasisymmetry to high accuracy. Using this construction, we give the first numerical demonstrations of Garren and Boozer’s ideal scaling of quasisymmetry breaking with the cube of the inverse aspect ratio. We also demonstrate a strongly non-axisymmetric configuration (vacuum rotational transform $\unicode[STIX]{x1D704}>0.4$) in which symmetry-breaking mode amplitudes throughout a finite volume are ${<}2\times 10^{-7}$, the smallest ever reported. To generate boundary shapes of finite-minor-radius configurations, a careful analysis is given of the effect of substituting a finite minor radius into the near-axis expansion. The approach here can provide analytic insight into the space of possible quasisymmetric stellarator configurations, and it can be used to generate good initial conditions for conventional stellarator optimization.

We show analytically that for $\unicode[STIX]{x1D704}$-profiles similar to the one of the Wendelstein 7-X stellarator, where $\unicode[STIX]{x1D704}$ is the rotational transform of the equilibrium magnetic field, a highly conducting toroidal plasma is unstable to kinetically mediated pressure-driven long-wavelength reconnecting modes, of the infernal type. The modes are destabilized either by the electron temperature gradient or by a small amount of current, depending on how far from unity the average value of $\unicode[STIX]{x1D704}$ is, which is assumed to be slowly varying. We argue that, for W7-X, a broad mode with toroidal and poloidal mode numbers $(n,m)=(1,1)$ can be destabilized due to the strong geometric side-band coupling of the resonant kinetic electron response at locations where $\unicode[STIX]{x1D704}$ is rational for harmonics that belong to the mode family of the $(n,m)=(1,1)$ mode itself. In many regimes, the growth rate is insensitive to the plasma density, thus it is likely to persist in high performance W7-X discharges. For a peaked electron temperature, with a maximum of $T_{e}=5~\text{keV}$, larger than the ion temperature, $T_{i}=2.5~\text{keV}$, and a density $n_{0}=10^{19}~\text{m}^{-3}$, instability is found in regimes which show plasma sawtooth activity, with growth rates of the order of tens of kiloHertz. Frequencies are either electron diamagnetic or of the ideal magnetohydrodynamic type, but sub-Alfvénic. The kinetic infernal mode is thus a good candidate for the explanation of sawtooth oscillations in present-day stellarators and poses a new challenge to the problem of stellarator reactor optimization.

The magnetic presheath is a boundary layer occurring when magnetized plasma is in contact with a wall and the angle $\unicode[STIX]{x1D6FC}$ between the wall and the magnetic field $\boldsymbol{B}$ is oblique. Here, we consider the fusion-relevant case of a shallow-angle, $\unicode[STIX]{x1D6FC}\ll 1$, electron-repelling sheath, with the electron density given by a Boltzmann distribution, valid for $\unicode[STIX]{x1D6FC}/\sqrt{\unicode[STIX]{x1D70F}+1}\gg \sqrt{m_{\text{e}}/m_{\text{i}}}$, where $m_{\text{e}}$ is the electron mass, $m_{\text{i}}$ is the ion mass, $\unicode[STIX]{x1D70F}=T_{\text{i}}/ZT_{\text{e}}$, $T_{\text{e}}$ is the electron temperature, $T_{\text{i}}$ is the ion temperature and $Z$ is the ionic charge state. The thickness of the magnetic presheath is of the order of a few ion sound Larmor radii $\unicode[STIX]{x1D70C}_{\text{s}}=\sqrt{m_{\text{i}}(ZT_{\text{e}}+T_{\text{i}})}/ZeB$, where e is the proton charge and $B=|\boldsymbol{B}|$ is the magnitude of the magnetic field. We study the dependence on $\unicode[STIX]{x1D70F}$ of the electrostatic potential and ion distribution function in the magnetic presheath by using a set of prescribed ion distribution functions at the magnetic presheath entrance, parameterized by $\unicode[STIX]{x1D70F}$. The kinetic model is shown to be asymptotically equivalent to Chodura’s fluid model at small ion temperature, $\unicode[STIX]{x1D70F}\ll 1$, for $|\text{ln}\,\unicode[STIX]{x1D6FC}|>3|\text{ln}\,\unicode[STIX]{x1D70F}|\gg 1$. In this limit, despite the fact that fluid equations give a reasonable approximation to the potential, ion gyro-orbits acquire a spatial extent that occupies a large portion of the magnetic presheath. At large ion temperature, $\unicode[STIX]{x1D70F}\gg 1$, relevant because $T_{\text{i}}$ is measured to be a few times larger than $T_{\text{e}}$ near divertor targets of fusion devices, ions reach the Debye sheath entrance (and subsequently the wall) at a shallow angle whose size is given by $\sqrt{\unicode[STIX]{x1D6FC}}$ or $1/\sqrt{\unicode[STIX]{x1D70F}}$, depending on which is largest.

A gyrokinetic Coulomb collision operator is derived, which is particularly useful to describe the plasma dynamics at the periphery region of magnetic confinement fusion devices. The derived operator is able to describe collisions occurring in distribution functions arbitrarily far from equilibrium with variations on spatial scales at and below the particle Larmor radius. A multipole expansion of the Rosenbluth potentials is used in order to derive the dependence of the full Coulomb collision operator on the particle gyroangle. The full Coulomb collision operator is then expressed in gyrocentre phase-space coordinates, and a closed formula for its gyroaverage in terms of the moments of the gyrocentre distribution function in a form ready to be numerically implemented is provided. Furthermore, the collision operator is projected onto a Hermite–Laguerre velocity space polynomial basis and expansions in the small electron-to-ion mass ratio are provided.