Kinetic simulations of sheared flow stabilization in high-temperature Z-pinch plasmas

Conventional strategies to reach fusion relevant conditions encompass a broad range of plasma densities and confinement times from high-density short-lived plasmas in inertial confinement fusion (ICF) schemes to low-density magnetic confinement fusion (MCF) plasmas with long confinement times.¹ ICF uses a convergent implosion to compress plasmas which heat to multi-keV temperatures and densities¹ up to 10²⁵ cm⁻³. The plasma confinement terminates on hydrodynamic time scales of several nanoseconds, which requires high input power,¹ up to 10¹⁴ W, and thus large-scale drivers. MCF plasmas are maintained at low density, 10¹⁴ cm⁻³, for seconds or more by a magnetic field that enhances alpha energy deposition and inhibits thermal conduction.¹ The need to mitigate thermal losses over these long periods has motivated large-scale device designs. Magneto-inertial confinement fusion (MICF) utilizes compressional heating and thermal transport inhibition by magnetic fields in plasmas with densities and confinement times that fall between the extremes of ICF and MCF. There are several promising MICF approaches,^1–4 and this paper investigates an MICF approach with a simple linear design and favorable scaling,⁵ the sheared-flow stabilized Z-pinch (SFS Z-pinch). The SFS Z-pinch could reduce the reactor size and cost of an MICF design relative to conventional designs by avoiding the prolonged thermal insulation requirements of MCF and the power input requirements of ICF. This design is scaling to fusion conditions and this paper describes the first fully kinetic simulations of the SFS Z-pinch which reproduce the observed stability at present-day experimental conditions and predict equivalent stability at reactor scales.

Z-pinches were among the first fusion reactor candidates,⁶ but they consistently form MHD instabilities which disrupt the current flow and terminate the plasma compression and confinement, e.g., m = 0 sausage and m = 1 kink modes where m is the azimuthal mode number.^5,7 These instabilities can facilitate beam-target fusion by accelerating ion beams which create short, bright neutron sources in gas-puff Z-pinches^8,9 and dense plasma foci (DPF).¹⁰ However, these instabilities prevent sustained thermonuclear fusion and the possibility of energy gain due to their rapid growth on time scales of r₀/v_A over small volumes with axial extent λ ∼ r₀. Here, r₀ is the characteristic pinch radius, λ is the axial wavelength of the instabilities, and $vA=B02/4πmin0$ is a representative Alfvén velocity taking the peak magnetic field B₀ and density n₀. A promising approach to stabilize these MHD instabilities in the Z-pinch is the use of sheared axial flows, ∂_rv_z ≠ 0.^11–20 Experimental studies of the SFS Z-pinch^11–14 report the suppression of m = 1, 2, 3 fluctuations to a low amplitude and stable Z-pinch operation for over 1000r₀/v_A. The stable period coincides with the presence of shear in axial flow measurements,^12,21 and this design has been scaled to 200 kA of pinch current^13,14 generating peak plasma temperatures up to 1 keV at densities up to 10¹⁷ cm⁻³. The Fusion Z-pinch Experiment (FuZE) at the University of Washington is currently scaling up to 300 kA of plasma current, and if stability persists up to 1–1.5 MA of pinch current, the resulting multi-keV plasmas at 10²⁰ cm⁻³ could reach positive fusion gain with pinch lengths of 0.1–1 m. Up to this point, SFS Z-pinch models have used fluid approximations and explain the stability as an m = 1 stabilization by sheared flows¹⁵ and stabilization of m = 0 modes through the plasma pressure profile as described by the Kadomtsev criterion.²² These fluid models are invalid in the high-temperature collisionless fusion plasmas of interest, and controlling the plasma pressure profile is a restrictive design constraint for future devices. This paper describes the first fully kinetic simulations of the SFS Z-pinch which show that the m = 0 modes are suppressed by sheared flows consistent with previous experiments.

An MHD model of SFS Z-pinch stability found that m = 1 modes are stabilized in sheared-flow profiles satisfying ∂_rv_z ≥ 0.1kv_A, where k is the wavevector of the mode.¹⁵ This stabilization model has been validated in experiments measuring m = 1 activity and shear^11,12 for a reasonable choice of the dominant mode, k = π/r₀, at multiple plasma currents. Efforts to produce a similar sheared-flow stabilization model for the m = 0 mode^16–20 find that stabilization requires shear ranging from 0.2 to 0.56kv_A, which necessitates supersonic flows ranging from Mach 1.7 to Mach 4.3. This is stronger than measured flows^11–13 in experimental plasmas with dramatically extended lifetimes, suggesting that the m = 0 suppression is due to the plasma pressure profile as described by the Kadomtsev criteria.²² Simulations with a kinetic ion treatment have shown^23,24 substantial changes in Z-pinch stability relative to the predictions of fluid theory and have revealed species separation effects and neutron production by nonmaxwellian distributions in ICF plasmas.^25–28 In gas puff Z-pinches⁹ and DPFs,^8,10,29 fully kinetic simulations have been used to predict the development of instabilities and ion beams for neutron source applications in optimized full-device simulations using a binary neutron algorithm.³⁰ Rather than full-device modeling, as is done in DPFs or gas puff Z-pinches, the present study uses small-scale, high-resolution, low-noise simulations of assembled SFS Z-pinch plasmas to produce a stability model for the SFS Z-pinch in kinetic regimes. In the multi-keV conditions of a fusion reactor, electrons are also collisionless and require a fully kinetic treatment. The simulations described here are appropriate for these conditions and show the efficacy of shear-flow m = 0 suppression under conditions ranging from the limits of present-day experimental capabilities to the projected conditions^5,31 of a SFS Z-pinch reactor with peak flow speeds $≲$ Mach 1, consistent with the experiment.

The Z-pinches are described by the Bennett³² profile

Here, B₀ is the maximum magnetic field, n₀ is the axial density of ions, T and J/n are constant, and r₀ is the characteristic radius of the pinch. This effective pinch radius encloses half of the total current and mass of the plasma and corresponds to the peak density scale length, $|n/(dn/dr)|$⁠. The plasma extends to a conducting boundary at 4r₀ with a periodic boundary condition along z. The radial boundary condition can substantially impact the plasma stability and m = 0 profile,^19,33 particularly when comparing free-boundary and internal m = 0 modes, i.e., a vacuum vs conducting boundary at the maximum radius. In SFS Z-pinch experiments, the plasma is surrounded by a cylindrical outer electrode with a radius over 10 times larger than the plasma pinch radius, as determined by digital holographic interferometry (DHI) diagnostics.^12,13 Simulation with a conducting boundary at 10r₀ and 4r₀ has good agreement, consistent with a similar test by Arber,²³ so a boundary at 4r₀ is used in the present study to reduce cost.

The simulations are carried out using the particle-in-cell (PIC) code Chicago.^34,35 The direct-implicit solver³⁶ option in Chicago is used to allow a much larger time step and cell size than a typical explicit scheme, while preserving low frequency oscillations like m = 0. The simulations contain 100 cells along r, dr ∼ ρ_Li0/4, axial resolution based on the wavelength of interest, 0.06ρ_Li0 < dz < 0.12ρ_Li0, time step dictated by the electron cyclotron frequency, $dt<0.7ωce0−1$⁠, and up to 3 × 10⁴ virtual particles per cell. Here, ρ_Li0 is the smallest ion Larmor radius in the plasma and ω_ce0 is the largest electron cyclotron frequency in the plasma which changes as the magnetic field perturbations grow. The large particle number reduces PIC density noise and lowers the numerical electron kinetic energy losses to 20% over the 10r₀/v_A simulation time scale. The total energy and current conservation are 97% and 99%, respectively. Chicago offers two collision models: a cumulative binary collision³⁷ model and a reduced scattering model that assumes drifting-Maxwellian particle distributions.³⁸ In the high temperature environments of this study, there is a negligible difference in the m = 0 behavior with the different collision models, and so the latter is adopted to reduce cost.

The Z-pinch is initialized with an m = 0 perturbation seed with a single wavelength along z, L = 2π/k,

where ϵ is the initial perturbation amplitude, which is 10⁻² in the present study. The growth of the m = 0 modes is determined by Fourier analyzing the density

and finding the exponential growth of

where $N=∫n 2πr drdz/L$ is the linear density averaged over the axial length L. The sheared flows, v_z(r), convect growing perturbations in the z direction, which can transform the bottlenecking density profile seen in static Z-pinches into a collection of low-density rings that are shifted axially based on the flow speed at the given radius. These axial shifts appear through the phase of the Fourier amplitudes

which are removed when assessing the overall mode growth by taking the complex magnitude of δn_k(r, t) before integrating radially.

The PIC simulations are benchmarked against shear-free m = 0 growth rate calculations from a Vlasov-fluid simulation study²⁴ of a shear-free Z-pinch using a continuum kinetic ion treatment with warm fluid electrons, shown in Fig. 1. The PIC simulations have good agreement with the kinetic results²⁴ and capture the suppression of large-k modes which does not occur in MHD models.^6,40 The mode amplitude from the PIC simulations is shown in Fig. 2 along with the unseeded second harmonic with wavelength λ = L/2. The perturbations are seeded with an initial amplitude of ϵ = 10⁻², which is near the low amplitude limit set by the unseeded harmonics of the simulation, shown in Fig. 2. These harmonics consist of PIC density noise that results from asymmetries in the random velocity sampling at the beginning of the simulation. As seen in Fig. 2, the noise and thus the harmonics decrease in amplitude as the number of particles increases, but the increased computational expense does not justify the noise improvement around 16 000 to 32 000 particles per cell. Although a short linear growth period reduces the confidence in growth rate calculations, the kinetic treatment is preferable for the collisionless or weakly collisional plasmas of interest and an investigation of the stabilization of m = 0 modes in sheared axial flows does not necessitate a long linear growth time. After verifying the m = 0 growth rates over the amplitude range accessible in the PIC simulations, the plasma equilibrium conditions relevant to the FuZE are examined.

The plasma equilibrium conditions from the FuZE and the projected conditions of a sheared-flow stabilized reactor are given in Table I. The projected reactor equilibrium is found by scaling the FuZE equilibrium from 300 kA to 1500 kA, with a fixed linear density, N = 1.1 × 10¹⁷ cm⁻¹.

The m = 0 growth rates from PIC simulations of the projected FuZE equilibrium are shown in Fig. 3 along with the MHD prediction. The PIC results and the MHD results in Fig. 3 use the Bennett plasma profile given by Eq. (1). These simulations have no sheared flows, and each point represents an independent simulation of one wavelength of the desired mode. The FuZE and reactor equilibria have comparable normalized growth rates with a peak growth rate at kr₀ = 5.

A rough estimate³⁹ of the shear required for stability is ∂_rv_z ≥ γ, which is maximum in the limit k → ∞ in MHD.^6,40 The finite Larmor radius stabilization captured through the kinetic treatment of the ions shows the suppression of large-k modes relative to the MHD prediction and reveals the strongest m = 0 modes with kρ_Li0 = kr₀/5.8 ≃ 1, where ρ_Li0 is the ion Larmor radius in the maximum B field. This is consistent with observations⁴¹ of macroscopic sausage modes with 2π/k ∼ r₀, leading to a Z-pinch current disruption rather than large-k modes with 2π/k ≪ r₀. With this result, sheared flow stability criteria for the strongest m = 0 instabilities can be determined by examining the range of wavelengths which mark the large-k bounds of MHD validity.

A linear flow profile, v_z(r) = v₀r/r₀, is given to the ions and electrons in simulations seeded with the strongest mode in the shear-free equilibrium, kr₀ = 5. The plasma develops a radial electric field when the sheared axial flows are initialized in the plasma, $Er(r)=−v0rr0Bθ(r)$⁠. For the conditions of these simulations, this electric field causes negligible charge separation, and the plasma pressure and current profiles retain the Bennett profile shape given by Eq. (1). Ion density contours are shown in Fig. 4, where the different rows show the increasing suppression of m = 0 growth with increasing shear, up to v₀ = 0.75v_A in the bottom row. Here, v_A is defined based on the peak equilibrium magnetic field and density, $vA=B0/μ0min0$⁠. The m = 0 mode amplitudes in these simulations, given by Eq. (3), are shown in Fig. 5 along with a stronger shear profile, v₀ = v_A.

The apparent saturation of the shear-free m = 0 amplitude, given by the black curve in Fig. 5, is a consequence of the normalization. As seen in the top row of Fig. 4, the radial confinement in the shear-free FuZE equilibrium continues to degrade through this saturation period. In the v₀ = v_A profile in Fig. 5, the m = 0 damps to the noise level of the simulation but cannot damp further. The noise level is given by the maximum amplitude among the unseeded axial modes, 2π/k= L/2, L/3, …. The amplitude of these harmonics increases over the simulation time scale, but no coherent mode appears at these wavelengths. The harmonics have small radial structures 2 to 4 cells wide along r, 2dr–4dr, and remain under-resolved in simulations with a higher radial resolution. When the harmonic wavelengths are seeded above the noise level of the simulation, they exhibit the suppression and damping behavior shown in Fig. 5. The minimum shear required for stability is difficult to determine over the limited time scales available in PIC simulations, but the mode suppression in the v₀ = 0.75v_A flow profile, which damps to roughly the initial seed amplitude, is most similar to a limiting shear in this survey of shear strengths. The shear stabilization scan was repeated in the projected conditions of a fusion reactor given by Table I and is plotted in Fig. 6. The similarity between the mode damping in FuZE and reactor conditions suggests that the sheared flow stabilization in the environment of a reactor has equivalent efficacy to the stabilization in the ongoing FuZE.

The shear in the stabilizing flow profile, v_z(r) = 0.75v_Ar/r₀, extends from the axis to the maximum simulation radius, shown in Fig. 7(bottom, purple dashed) normalized to the sound speed $cs2=(Ti+5Te/3)/mi$⁠. The initial flow profile is well preserved over the simulation time scale aside from a flattening near the axis due to the higher viscosity⁴² in the unmagnetized region. The flattening occurs over a much faster time scale than the m = 0, and the flow profile at t = 5r₀/v_A is shown in Fig. 7(bottom, solid purple).

In the Bennett profile, the plasma extends to the maximum simulation radius, and the m = 0 fluctuation amplitude, δn(r) peaks at the pinch radius r₀ with a relatively narrow radial localization as seen in Fig. 5(top, gray shaded). Mode stabilization does not require shear at large radii beyond the localization of the m = 0 fluctuations, and equivalent m = 0 suppression occurred in the flow profile shown in black in Fig. 7(bottom, black). As with the linear flow profiles, the dashed curve shows the initial profile, and the solid curve shows the flow profile at t = 5r₀/v_A which has flattened on axis. This localized flow profile stabilizes the plasma with a significantly lower peak flow speed than the linear flow profile and is more consistent with measured flow profiles in previous experiments. A flow profile measured with ion spectroscopy in the ZaP experiment¹² is plotted in Fig. 7(bottom, blue), normalized to the plasma radius and sound speed. The measurement is taken along chords with positive and negative impact parameters, relative to the peak plasma density, and the simulated profiles are mirrored across the axis for comparison. Although a lower limit on the sheared flow necessary for stability was not extensively searched, the simulated profile (bottom, black) reproduces the stabilization observed in earlier experiments in flows comparable to the measured flows with peak flow speeds $≲$ Mach 1.

Simulations of other wavelengths in the v₀ = 0.75v_A flow profile show faster and stronger overall damping in short wavelength modes, kr₀ > 5. The long-wavelength instabilities damped to around the noise level in the v₀ = v_A profile, but the limiting stabilizing behavior of the v₀ = 0.75v_A could not be reproduced over the extended time and spatial scales of the long-wavelength modes. If the sheared flows become degraded by Kelvin Helmholtz⁴³ instabilities (KHI) or microturbulence, the Z-pinch stability would be compromised. No KHI activity is seen in these PIC simulations, although simulations of longer wavelength modes may reveal KHI activity. The simulations have axial resolution up to dz = ρ_Li0/16.4, and temporal resolution dt ∼ m_ec/eB₀ capable of resolving drift instabilities and micro-instabilities. No microturbulence or drift-wave activity is seen, and the stabilizing flow profile in Fig. 7(black) has diminishing shear in the low β regions where these modes are operative.

MHD simulations^11,21,44 also report improved stability in nonuniform shear profiles, e.g., $vz=v0r2/r02$⁠, when the peak shear is localized in the radial vicinity of the peak amplitude of m = 0 fluctuations, δn(r). These studies use slab and parabolic profiles, in which the plasma pressure gradient or plasma pressure goes to zero beyond the pinch radius, r₀, and the peak m = 0 fluctuation amplitude is located near the pinch radius. For the diffuse Bennett profile used in the present study, the plasma pressure and m = 0 fluctuation amplitude extend smoothly beyond the pinch radius, where the mode peaks, to the simulations boundary at 4r₀. This increases the radial width of the modes and the necessary radial extent of the shear. Despite the larger mode width, the localized shear profile given by the black curve in Fig. 7 suppresses the m = 0 with peak flow speeds $≲$ Mach 1.

There is substantial experimental validation^11–13 of the sheared-flow stabilization of m = 1 modes when¹⁵ ∂_rv_z ≥ 0.1kv_A. The suppression of the strongest m = 0 modes in the present study requires more shear, ∂_rv_z = 0.15kv_A, but as seen in the measured flow profile, such flows are within the device capabilities. This is in contrast to previous fluid models, which predict stability in sheared flows with supersonic speeds up to Mach 4.3 or in plasma pressure profiles satisfying the Kadomtsev criteria. In summary, PIC simulations reveal the wavelengths of the dominant m = 0 modes, and the stabilization of these modes in sheared flows with ∂_rv_z = 0.75v_A/r₀ and peak flows $≲$ Mach 1. Stability occurs in flow profiles consistent with experimental measurements with similar efficacy at scales ranging from present experimental capabilities up to a reactor.

The information, data, or work presented herein was funded in part by the Advanced Research Projects Agency - Energy (ARPA-E), U.S. Department of Energy, under Award No. DE-AR-0000571. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Computing support was provided by the LLNL Computing Grand Challenge Program. U. Shumlak gratefully acknowledges support of the Erna and Jakob Michael Visiting Professorship at the Weizmann Institute of Science and as a Faculty Scholar at the Lawrence Livermore National Laboratory. Special thanks to the experimental team for the valuable discussions and input, Elliot L. Claveau, Ellie G. Forbes, James M. Mitrani, Anton D. Stepanov, Tobin R. Weber, and Yue Zhang.