The first fully kinetic particle-in-cell (PIC) simulations of sheared flow stabilized Z-pinch plasmas show the suppression of the sausage instability by shear, ∂rvz ≠ 0, with flow Mach numbers , consistent with experimental observations. Experimental investigations of sheared-flow stabilized Z-pinches demonstrated stability for 10 s of microseconds, over 1000 Alfvén radial transit times, in quasi steady-state plasmas that are an intermediate between conventional inertial and magnetic confinement systems. The observed stability coincides with the presence of radial shear in axial flow profiles with peak speeds less than Mach 1, and experiments are underway to validate scaling this design to fusion conditions. The experimentally observed stability agrees with models of m = 1 kink mode suppression by sheared flows, but existing models of the m = 0 sausage mode underestimate the efficacy of sheared flow stabilization. These models rely on fluid approximations and find that stabilization requires flows ranging from Mach 1.7 to 4.3, and in some cases, stabilization is not reproduced in the models. This is faster than the measured flows in long-lived plasmas and would necessitate substantial energy convection out of the Z-pinch and the need to drive and sustain supersonic flows in future devices. The MHD models typically used in the literature are invalid in the high-temperature, high-current environments desirable for many Z-pinch applications, and they ignore large Larmor radius effects and viscous dissipation which are known to impact Z-pinch stability. PIC simulations can capture all these effects as well as kinetic instabilities that could influence the performance of high-temperature sheared flow stabilized Z-pinch plasmas. The PIC simulations presented here show the suppression and damping of m = 0 modes by sheared flows ∂rvz = 0.75vA/r0 with flow Mach numbers . Equivalent stability occurs under plasma conditions ranging from the limits of present-day experimental capabilities to the projected conditions of a sheared flow stabilized Z-pinch reactor.
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.1 ICF uses a convergent implosion to compress plasmas which heat to multi-keV temperatures and densities1 up to 1025 cm−3. The plasma confinement terminates on hydrodynamic time scales of several nanoseconds, which requires high input power,1 up to 1014 W, and thus large-scale drivers. MCF plasmas are maintained at low density, 1014 cm−3, for seconds or more by a magnetic field that enhances alpha energy deposition and inhibits thermal conduction.1 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,5 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,6 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-pinches8,9 and dense plasma foci (DPF).10 However, these instabilities prevent sustained thermonuclear fusion and the possibility of energy gain due to their rapid growth on time scales of r0/vA over small volumes with axial extent λ ∼ r0. Here, r0 is the characteristic pinch radius, λ is the axial wavelength of the instabilities, and is a representative Alfvén velocity taking the peak magnetic field B0 and density n0. A promising approach to stabilize these MHD instabilities in the Z-pinch is the use of sheared axial flows, ∂rvz ≠ 0.11–20 Experimental studies of the SFS Z-pinch11–14 report the suppression of m = 1, 2, 3 fluctuations to a low amplitude and stable Z-pinch operation for over 1000r0/vA. 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 current13,14 generating peak plasma temperatures up to 1 keV at densities up to 1017 cm−3. 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 1020 cm−3 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 flows15 and stabilization of m = 0 modes through the plasma pressure profile as described by the Kadomtsev criterion.22 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 ∂rvz ≥ 0.1kvA, where k is the wavevector of the mode.15 This stabilization model has been validated in experiments measuring m = 1 activity and shear11,12 for a reasonable choice of the dominant mode, k = π/r0, at multiple plasma currents. Efforts to produce a similar sheared-flow stabilization model for the m = 0 mode16–20 find that stabilization requires shear ranging from 0.2 to 0.56kvA, which necessitates supersonic flows ranging from Mach 1.7 to Mach 4.3. This is stronger than measured flows11–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.22 Simulations with a kinetic ion treatment have shown23,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-pinches9 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.30 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 conditions5,31 of a SFS Z-pinch reactor with peak flow speeds Mach 1, consistent with the experiment.
II. PLASMA EQUILIBRIUM AND SIMULATION SETUP
The Z-pinches are described by the Bennett32 profile
Here, B0 is the maximum magnetic field, n0 is the axial density of ions, T and J/n are constant, and r0 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, . The plasma extends to a conducting boundary at 4r0 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 10r0 and 4r0 has good agreement, consistent with a similar test by Arber,23 so a boundary at 4r0 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 solver36 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, , and up to 3 × 104 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 10r0/vA simulation time scale. The total energy and current conservation are 97% and 99%, respectively. Chicago offers two collision models: a cumulative binary collision37 model and a reduced scattering model that assumes drifting-Maxwellian particle distributions.38 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−2 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 is the linear density averaged over the axial length L. The sheared flows, vz(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 δnk(r, t) before integrating radially.
III. SIMULATION BENCHMARK AND KINETIC EFFECTS IN SHEAR-FREE Z-PINCHES
The PIC simulations are benchmarked against shear-free m = 0 growth rate calculations from a Vlasov-fluid simulation study24 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 results24 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−2, 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 × 1017 cm−1.
|Device .||I (kA) .||n0 (cm−3) .||T (keV) .||r0 (mm) .||νiir0/vA .|
|FuZE||300||4.25 × 1018||1.27||0.91||0.087|
|Reactor||1500||1.41 × 1020||31.7||0.158||0.0008|
|Device .||I (kA) .||n0 (cm−3) .||T (keV) .||r0 (mm) .||νiir0/vA .|
|FuZE||300||4.25 × 1018||1.27||0.91||0.087|
|Reactor||1500||1.41 × 1020||31.7||0.158||0.0008|
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 kr0 = 5.
A rough estimate39 of the shear required for stability is ∂rvz ≥ γ, 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 = kr0/5.8 ≃ 1, where ρLi0 is the ion Larmor radius in the maximum B field. This is consistent with observations41 of macroscopic sausage modes with 2π/k ∼ r0, leading to a Z-pinch current disruption rather than large-k modes with 2π/k ≪ r0. 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.
IV. SHEARED-FLOW STABILIZATION
A linear flow profile, vz(r) = v0r/r0, is given to the ions and electrons in simulations seeded with the strongest mode in the shear-free equilibrium, kr0 = 5. The plasma develops a radial electric field when the sheared axial flows are initialized in the plasma, . 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 v0 = 0.75vA in the bottom row. Here, vA is defined based on the peak equilibrium magnetic field and density, . The m = 0 mode amplitudes in these simulations, given by Eq. (3), are shown in Fig. 5 along with a stronger shear profile, v0 = vA.
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 v0 = vA 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 v0 = 0.75vA 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.
V. STABILITY IN LOCALIZED FLOW PROFILES
The shear in the stabilizing flow profile, vz(r) = 0.75vAr/r0, extends from the axis to the maximum simulation radius, shown in Fig. 7(bottom, purple dashed) normalized to the sound speed . The initial flow profile is well preserved over the simulation time scale aside from a flattening near the axis due to the higher viscosity42 in the unmagnetized region. The flattening occurs over a much faster time scale than the m = 0, and the flow profile at t = 5r0/vA 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 r0 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 = 5r0/vA 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 experiment12 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 v0 = 0.75vA flow profile show faster and stronger overall damping in short wavelength modes, kr0 > 5. The long-wavelength instabilities damped to around the noise level in the v0 = vA profile, but the limiting stabilizing behavior of the v0 = 0.75vA could not be reproduced over the extended time and spatial scales of the long-wavelength modes. If the sheared flows become degraded by Kelvin Helmholtz43 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 ∼ mec/eB0 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 simulations11,21,44 also report improved stability in nonuniform shear profiles, e.g., , 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, r0, 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 4r0. 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 validation11–13 of the sheared-flow stabilization of m = 1 modes when15 ∂rvz ≥ 0.1kvA. The suppression of the strongest m = 0 modes in the present study requires more shear, ∂rvz = 0.15kvA, 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 ∂rvz = 0.75vA/r0 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.