Z Pinch Kinetics II - A Continuum Perspective: Betatron Heating and Self-Generation of Sheared Flows

Doubly adiabatic CGL theory describes collisionless compression of magnetized plasma, and predicts different degrees of parallel and perpendicular heating due to magnetic flux compression Chew, Goldberger, and Low (1956). The classic CGL model consists of two invariants, one for each pressure component

	$\displaystyle\frac{d}{dt}\Big{(}\frac{p_{\perp}}{\rho B}\Big{)}$	$\displaystyle=$	$\displaystyle 0,$		(4a)
	$\displaystyle\frac{d}{dt}\Big{(}\frac{p_{\parallel}B^{2}}{\rho^{3}}\Big{)}$	$\displaystyle=$	$\displaystyle 0,$		(4b)

in contrast to the singly adiabatic (conserved specific entropy) invariant $\frac{d}{dt}\big{(}\frac{p}{\rho^{\gamma}}\big{)}=0$.

Equations 4 hold in the cylindrical Z pinch for fluid elements which are made up primarily of cyclotron trajectories. This is demonstrated simply as follows. Suppose that a plasma fluid element is characterized by three adiabatic invariants, namely specific magnetic moment $\mu$, specific angular momentum $L_{\theta}$, and specific magnetic flux $s_{m}$ (from Alfvén’s law)

	$\displaystyle\frac{d\mu}{dt}=0$	$\displaystyle\quad\implies\quad$	$\displaystyle\frac{d}{dt}\Big{(}\frac{v_{\mathrm{th},\perp}^{2}}{B_{\theta}}\Big{)}=0,$		(5a)
	$\displaystyle\frac{dL_{\theta}}{dt}=0$	$\displaystyle\quad\implies\quad$	$\displaystyle\frac{d}{dt}(rv_{\mathrm{th},\parallel})=0,$		(5b)
	$\displaystyle\frac{ds_{m}}{dt}=0$	$\displaystyle\quad\implies\quad$	$\displaystyle\frac{d}{dt}\Big{(}\frac{B_{\theta}}{\rho r}\Big{)}=0$		(5c)

where $v_{\mathrm{th},\parallel}=\sqrt{k_{B}T_{\parallel}/m}$ is the parallel thermal velocity and $v_{\mathrm{th},\perp}=\sqrt{k_{B}T_{\perp}/m}$ the perpendicular one. Then substituting the parallel and perpendicular pressures ($p_{\parallel}=nkT_{\parallel}$ and $p_{\perp}=nkT_{\perp}$) for the thermal velocities and combining Eqs. 5 yields exactly Eqs. 4, the CGL equations.

Equations 5 clarify the double-adiabatic CGL model in cylindrical geometry as consisting of the magnetic moment plus the angular momentum invariant, versus the Cartesian picture in which the second invariant is the mirror bounce action $J_{\parallel}\sim v_{\parallel}\ell$ along a field-aligned path of length $\ell$. The singly adiabatic model (invariant specific entropy) arises in the collisional limit $\omega_{ci}\tau\ll 1$. Observe that the CGL model is truly triply adiabatic considering the necessity of Alfvén’s law to obtain the CGL equations and the thermodynamic character of the specific magnetic flux Crews et al. (2024).

Compression in the Cartesian slab model and in the cylindrical pinch are distinguished by the Jacobian appearing in Alfvén’s law Turchi and Shumlak (2023). Namely, the Cartesian slab has $\frac{d}{dt}\big{(}\frac{B}{\rho}\big{)}=0$ while the cylinder has $\frac{d}{dt}\big{(}\frac{B_{\theta}}{\rho r}\big{)}=0$, a subtle difference which can lead to opposite expectations. That is, combining Eqs. 4 with these frozen flux invariants produces the following equations,

	Slab,	$\displaystyle\frac{d}{dt}\Big{(}\frac{p_{\parallel}}{p_{\perp}}B\Big{)}=0,$			(6a)
	Cylinder,	$\displaystyle\frac{d}{dt}\Big{(}\frac{p_{\parallel}}{p_{\perp}}B_{\theta}r^{2}\Big{)}=0.$			(6b)

Cartesian slab anisotropy is simply inversely proportional to the flux compression ratio, $\frac{p_{\parallel,2}}{p_{\perp,2}}=C_{m}^{-1}$ (with $C_{m}\equiv\frac{B_{2}}{B_{1}}$) between an isotropic initial state $(\cdot)_{1}$ and an anisotropic compressed state $(\cdot)_{2}$. In the same vein, anisotropy in the cylinder follows from Eq. 6b as

$$\frac{p_{\parallel,2}}{p_{\perp,2}}=\Big{(}\frac{I_{1}}{I_{2}}\Big{)}\Big{(}\frac{r_{1}}{r_{2}}\Big{)}$$

(7)

where $\mu_{0}I\equiv 2\pi rB_{\theta}$ is the axial current enclosed up to radius $r$.

Thus, the anisotropy produced by flux compression of cyclotron trajectories in a Z pinch is the product of two factors: the current amplification ratio and the ratio of radii (alternatively, the product of the flux and area compression ratios). Flux compression in the Cartesian slab always leads to the anisotropy orientation $p_{\parallel}<p_{\perp}$, but the Z pinch will easily display the orientation $p_{\parallel}>p_{\perp}$ if current does not increase by a greater fraction than the radius. This azimuthal over-heating can happen because angular momentum conservation competes with magnetic moment conservation, and has important consequences for the resulting radial force balance.

To summarize the preceding, fluid elements composed of cyclotron trajectories display two adiabatic invariants (Eqs. 5), magnetic moment and angular momentum, and an additional frozen-in flux invariant. In the cylinder, magnetic moment is the invariant associated with radial and axial bounce of the cyclotron orbit, and the angular momentum with periodic azimuthal orbits.

Idealized betatron trajectories display invariants closely related to those of the ideal cyclotrons with a crucial difference: the lines of electric current are the principal symmetry axis rather than the magnetic field lines. Ergo, fluid elements composed of betatron trajectories should follow the same CGL equations as the cyclotrons, but the notions of parallel and perpendicular refer to the electrical current axis ($\hat{z}$) instead of the magnetic field lines, which we demonstrate below.

The adiabatic invariants of an ideal betatron orbit consist of its radial bounce action $J_{r}$, its axial bounce action $J_{z}$, and its azimuthal angular momentum $L_{\theta}$ (or action of the periodic azimuthal motion). These invariants are inferred from the ideal betatron orbit (cf. the prequel Crews, Meier, and Shumlak (2024)), which we present here for self-completeness.

The radial potential energy of a charged particle in the vicinity of the Z-pinch axis is

$$U=\frac{(P_{z}-q_{s}A_{z})^{2}}{2m_{s}}\approx\frac{P_{z}^{2}}{2m_{s}}-P_{z}\frac{q_{s}A_{z}}{m_{s}}$$

(8)

with $P_{z}$ the z-component of canonical momentum, provided that $P_{z}\gg-\frac{1}{2}q_{s}A_{z}$. Expanding the vector potential to $\mathcal{O}(r^{2})$ around the Z-pinch axis yields a harmonic potential $U\approx\frac{1}{2}kr^{2}$ describing a radial bounce with spring constant $k$. The axial bounce motion follows from conservation of canonical momentum. These radial and axial bounce motions are expressed as

	$\displaystyle r(t)$	$\displaystyle=$	$\displaystyle r_{0}\sin(\omega_{\beta}t),$		(9a)
	$\displaystyle mv_{z}(t)$	$\displaystyle=$	$\displaystyle P_{z}-\kappa\sin^{2}(\omega_{\beta}t).$		(9b)

The axial momentum fluctuation amplitude $\kappa\equiv q_{s}A_{z}(r_{0})$ is the potential momentum at the radial turning point $r_{0}$ and $\omega_{\beta}$ is the betatron frequency with $\omega_{\beta}^{2}\equiv\frac{q_{s}B_{\theta}}{m_{s}}\frac{v_{z}}{r}$. The radial and axial bounce motions contribute substantially to temperature in the species drift frame.

The radial bounce invariant follows as the area of the $(r,v_{r})$ phase-plane orbit,

$$J_{r}=\pi m_{s}r_{0}^{2}\omega_{\beta}\sim v_{r0}r_{0},$$

(10)

or the product of the radial bounce amplitude $r_{0}$ and the radial velocity amplitude $v_{r0}$. The radial bounce action is similar in form to the angular momentum or mirror bounce invariants, i.e. velocity times length, as a simple harmonic motion.

The axial bounce invariant is found by transforming to the bounce-averaged ($\langle\cdot\rangle_{T}$) axial drift frame, i.e. to the oscillation-center coordinates $(\delta z,\delta v_{z})$. The oscillation-center trajectory is

	$\displaystyle\langle m_{s}v_{z}\rangle_{T}$	$\displaystyle=$	$\displaystyle P_{z}-\frac{\kappa}{2},$		(11a)
	$\displaystyle\delta(m_{s}v_{z})\equiv m_{s}v_{z}-\langle m_{s}v_{z}\rangle_{T}$	$\displaystyle=$	$\displaystyle-\kappa\cos(2\omega_{\beta}t),$		(11b)

and the axial bounce invariant as the area of the oscillation-center phase-plane ellipse,

$$J_{z}=\pi\frac{\kappa^{2}/2m_{s}}{\omega_{\beta}}.$$

(12)

Equations 10 and 12 hold approximately for non-ideal betatron orbits, in a similar way that magnetic moment invariance is taken to hold for cyclotron orbits, although it is truly conserved only up to guiding-center corrections. We also consider axis-encircling betatron orbits to be approximately described by Eqs. 10 and 12 with the addition of an angular momentum invariant, $J_{\theta}=m_{s}r_{0}v_{\theta}$, because betatron orbits with azimuthal momentum may be approximated as an in-plane bounce orbit that rotates around the $\hat{z}$-axis Crews, Meier, and Shumlak (2024).

We now specialize to ion trajectories because their orbits are much larger than the electrons. The characteristic betatron frequency $\omega_{\beta}$ follows from $\omega_{\beta}^{2}\equiv\omega_{ci}\frac{v_{z}}{r}$ by setting $v_{z}$ equal to the ion axial drift $v_{i}$ (in the perfectly neutral Lorentz frame). Three useful forms of this frequency are

$$\omega_{\beta}=\frac{v_{\mathrm{th},i}}{r_{p}}=\frac{v_{i}}{\lambda_{i}}=\omega_{ci}\mathfrak{B}_{i}^{-1/2}$$

(13)

where $v_{\mathrm{th},i}$ is ion thermal velocity, $r_{p}$ is pinch radius, $\lambda_{i}=c/\omega_{pi}$ is ion inertial length, $\omega_{ci}$ is characteristic ion cyclotron frequency, and $\mathfrak{B}_{i}$ the ion Budker parameter Crews, Meier, and Shumlak (2024).

Applying Eq. 13 to Eq. 12 and supposing $\mathfrak{B}_{i}=\text{const.}$, we infer three Lagrangian invariants associated to an ion fluid element threaded by betatron orbits, in analogy to Eqs. 5, to be

	$\displaystyle\frac{d}{dt}\Big{(}\frac{v_{\mathrm{th},z}^{2}}{B_{\theta,0}}\Big{)}=0,$		(14a)
	$\displaystyle\frac{d}{dt}(rv_{\mathrm{th},r\theta})=0,$		(14b)
	$\displaystyle\frac{d}{dt}\Big{(}\frac{B_{\theta,0}}{\rho r}\Big{)}=0$		(14c)

having assumed $(r,\theta)$-isotropy, corresponding to the thermal axial invariant (Eq. 14a), radial invariant (Eq. 14b) and azimuthal angular momentum invariant (Eq. 14b). In Eqs. 14, $B_{\theta,0}$ is the characteristic flux density in the vicinity of the magnetic O-point about which the betatron orbits are bouncing. Combination of Eqs. 14 obtains the CGL model for an ion betatron fluid element,

	$\displaystyle\frac{d}{dt}\Big{(}\frac{p_{r\theta}}{\rho B_{\theta,0}}\Big{)}$	$\displaystyle=$	$\displaystyle 0,$		(15a)
	$\displaystyle\frac{d}{dt}\Big{(}\frac{p_{z}B_{\theta,0}^{2}}{\rho^{3}}\Big{)}$	$\displaystyle=$	$\displaystyle 0.$		(15b)

The electrical current axis assumes the principal anisotropy axis of the betatron fluid, and is otherwise modeled by the same CGL equations as a cyclotron plasma fluid. The radial and azimuthal directions display harmonic motions essentialized by a characteristic distance and a characteristic velocity (similar to the longitudinal mirror invariant), and the harmonic motion of the axial direction is essentialized by characteristic kinetic energy and the characteristic cyclotron frequency (like the magnetic moment of the cyclotrons). Anisotropic heating of the radial and axial directions is an intuitive consequence of the inductive axial electric field associated with flux compression splitting its energy between: 1) acceleration of the betatron axial drift, and 2) particle heating through increasing the energy of the betatron’s radial, axial, and azimuthal periodic motions.

Equations 14 come with important caveats. These invariants have been inferred based on the properties of “linear” betatron orbits, for which $P_{z}\gg-q_{s}A_{z}/2$, which translates to a very small trapping parameter in the prequel’s terminology Crews, Meier, and Shumlak (2024). This analysis based on the linear orbit’s invariants parallels invariance of the magnetic moment for the cyclotron trajectories only to lowest order in the finite Larmor radius (FLR) parameter, viz. trapping parameter. The action integrals of a particular orbit are truly nonlinear functions of its trapping parameter.

Indeed, for the Bennett distribution, nonlinearity is non-negligible for many near-axis orbits. Sufficiently close to the magnetic O-point, these adiabatic invariants are expressed in elliptic integrals. Such nonlinear expressions for the radial betatron action were given explicitly by Sonnerup Sonnerup (1971) in the context of a reconnecting current sheet. Sonnerup’s analysis identified a significant fact: adiabatic compression of the center of the current sheet causes all orbits (both cyclotrons and betatrons) to approach the same elliptic modulus (trapping parameter) $\alpha=2^{-1/4}\approx 0.84$ (at constant canonical momentum $P_{z}$). This trapping parameter $\alpha<1$ represents a nonlinear betatron oscillation Crews, Meier, and Shumlak (2024). Thus, a corollary of Sonnerup’s result, originally identified for adiabatic compression of a reconnecting current sheet, is that sufficient flux compression of the Z pinch drives all near-axis particles into a nonlinear betatron orbit.

Sonnerup’s result assumes that the non-adiabatic cyclotron-betatron separatrix crossing is adiabatic-in-mean, as argued by Janicke Janicke (1975). In reality, the separatrix crossing is chaotic, as explored by Cary et al. Cary, Escande, and Tennyson (1986) and Büchner Büchner and Zelenyi (1989), which results in compression timescale-dependent corrections to the spectrum of the resulting betatron orbit flux. We note however that such non-adiabatic cross-separatrix orbital variations are less significant for the compressing Z pinch than the reconnecting current sheet because the infinite-period separatrix only occurs for purely radial/axial bounce orbits. Non-zero angular momentum of axis-encircling orbits erases the infinite-period separatrix, although all trajectories passing directly through $r=0$ may undergo these chaotic non-adiabatic changes. Non-adiabatic transitions across the betatron-cyclotron separatrix of axis-crossing orbits in a compressing Z pinch deserves deeper investigation in a future study.

Sections III.1 and III.2 established that the CGL model applies to plasma composed of either cyclotron or betatron trajectories, in the usual way for the cyclotron plasma and a modified way for the betatron plasma. Specifically, the models predict the same magnitude of anisotropic heating from flux compression of a pinched plasma displaying either flavor of motion, but predict differing principal anisotropy axes; either the magnetic field lines for fluid consisting of cyclotron trajectories, or the electrical current lines for a fluid of betatron trajectories.

We now combine Sections III.1 and III.2 to propose an anisotropy model for a compressing Z pinch plasma, similar to the Ochs-Fisch model Ochs and Fisch (2018) but valid for all radii. The model consists of estimating a radially uniform magnitude of anisotropy due to flux compression using the CGL model and calculating the principal axis of the anisotropic Maxwellian thus generated in terms of the local partial pressures of cyclotron and betatron trajectories. Based on numerical experiments (Section X), this model appears accurate under the assumptions:

1. 1.

   The flux function varies slowly compared to the Alfvén transit time, and thus changes uniformly throughout the plasma;

2. 2.

   The change in particle magnetization during compression (due to separatrix crossing) is small enough that the relative densities of cyclotron and betatron trajectories are invariant.

Assumption 1 means that Eq. 6b may be used to estimate the magnitude of the generated anisotropy for all radii, and Assumption 2 that the anisotropy vector components in the final state may be parametrized using the initial state’s linear density.

Given Assumptions 1 and 2, the anisotropy model consists of the equations,

	$\displaystyle\alpha$	$\displaystyle\equiv$	$\displaystyle 1-\Big{(}\frac{I_{2}}{I_{1}}\Big{)}\Big{(}\frac{r_{2}}{r_{1}}\Big{)}$		(16a)
	$\displaystyle\alpha_{rz}(r)$	$\displaystyle=$	$\displaystyle\Big{(}\frac{n_{\beta}(r)}{n}\Big{)}\alpha,$		(16b)
	$\displaystyle\alpha_{r\theta}(r)$	$\displaystyle=$	$\displaystyle\Big{(}\frac{n_{c}(r)}{n}\Big{)}\alpha$		(16c)

where $n_{\beta}/n$ and $n_{c}/n$ are the radially varying betatron and cyclotron orbital density fractions (cf. the prequel for the analytic fractions computed for the Bennett pinch Crews, Meier, and Shumlak (2024)). The anisotropies appearing in Eqs. 16 are defined as

	$\displaystyle\alpha_{rz}$	$\displaystyle\equiv$	$\displaystyle 1-\frac{T_{z}}{T_{r}},$		(17a)
	$\displaystyle\alpha_{r\theta}$	$\displaystyle\equiv$	$\displaystyle 1-\frac{T_{r}}{T_{\theta}},$		(17b)

and are referred to as the agyrotropic and gyrotropic anisotropies in the following.

Figure 2 shows the gyrotropic and agyrotropic ion anisotropy profiles predicted by Eqs. 16 for various ion Budker parameters. The model predicts that agyrotropic anisotropy (induced by collisionless compression of ion betatron orbits) is dominant for small Budker parameters, while gyrotropic anisotropy (from collisionless compression of ion cyclotron trajectories) is dominant in the opposite limit. In the intermediate regime, $\mathfrak{B}_{i}=\mathcal{O}(10)$, ion agyrotropy terminates in the vicinity of the pinch radius, which has important consequences for the transport of axial momentum and diamagnetic response of the compressing plasma discussed in the remainder of this work.

In addition to adiabatic compression, agyrotropic anisotropy may be induced non-adiabatically due to collisionless shock heating Gedalin (2022) during fast Z pinch compression. Agyrotropy is induced in collisionless shock heating due to electric field gradients on the scale of the ion inertial length. In addition, radial electric field gradients induced by entropy mode fluctuations Ricci, Rogers, and Dorland (2006) and ensuing kinetic turbulence could also contribute to agyrotropy in ion distributions. Agyrotropic heating is frequently observed in simulations of kinetic plasma turbulence on ion inertial length scales Servidio et al. (2015); Yang et al. (2023). In this connection, note that species inertial length $\delta_{s}=c/\omega_{ps}$ is directly proportional to pinch radius through the species Budker parameter Mahajan (1989); $\mathfrak{B}_{s}=\frac{r_{p}^{2}}{4\lambda_{s}^{2}}$. Finally, strong electric fields can significantly modify particle orbits and drifts, particularly around magnetic O-points Joseph (2021).

Figure 3 presents a visual aid to the discussion of Section III.2, presenting an example of betatron heating and acceleration in a compressing azimuthal magnetic flux. An ion with $P_{z}>0$ is evolved in the model time-dependent $\hat{z}$-directed vector potential $A_{z}=-A_{0}(t)\frac{r^{2}}{2a^{2}}$, with $a=1$ and $A_{0}(t)$ a specified flux compression function. The ion is released from $r=1$ with $v_{r}=0$ and $v_{z}=1$, and the flux factor $A_{0}(t)$ is linearly and slowly varied from an initial value to a final value.

The inductive electric field $E_{z}=-\partial_{t}A_{z}$ does work on the ion, increasing both the radial bounce energy and the axial drift/bounce energy in step (minus $P_{z}^{2}/2m$). Most of the axial oscillation energy is channeled into the axial drift motion, increasing the axial bounce energy in the drifting frame much less than the radial bounce energy, because most of the axial energy in Fig. 3 is composed of the Galillean energy boost $\langle v_{z}\rangle\cdot\delta v_{z}$ from the oscillation-center frame to the lab frame. Thus, betatron heating displays radial/axial temperature anisotropy because the change in radial bounce energy is distinct from the change in axial bounce energy in the drifting frame.

During collisionless adiabatic flux compression of a large Larmor radius Z pinch, mean ion betatron drift changes self-consistently with plasma current and density, and mean radial bounce energy varies so that radial pressure balances the self-magnetic energy $\frac{\mu_{0}}{8\pi}I^{2}$ (Bennett relation). Axial temperature is unconstrained by radial force balance, but is constrained by kinetic equilibrium considerations (which was not appreciated in early Z-pinch studies, see for example the claim in Bennett 1934 that axial temperature is irrelevant to equilibrium Bennett (1934)).

To begin, we present the aforementioned method Suzuki and Shigeyama (2008) to generate a self-consistent kinetic equilibrium using the properties of Hermite polynomials. Consider the ansatz distribution

$$f(A_{z},H)=n_{0}\Big{(}\frac{\beta}{2\pi}\Big{)}^{3/2}\Big{(}\sum_{n=0}^{\infty}g_{n}(A_{z})H_{n}(v_{z})\Big{)}e^{-\beta H}$$

(18)

with $H=\frac{m}{2}(v_{x}^{2}+v_{y}^{2}+v_{z}^{2})+q\varphi$ the energy, $\beta=(k_{B}T)^{-1}$ the inverse temperature, $H_{n}$ the Hermite polynomials, $A_{z}$ the $z$-component of magnetic potential, and $\varphi$ the electric potential. Equation 18 is a kinetic equilibrium with $f=f(P_{z},H)$ provided that

$$\sum_{n=0}^{\infty}(g_{n}^{\prime}(a_{z})H_{n}(v_{z})-g_{n}(a_{z})H_{n}^{\prime}(v_{z}))=0$$

(19)

as found by substitution of Eq. 18 into the kinetic equation. Here the vector potential is written in lower case when it takes on the units of velocity, $a_{z}\equiv qA_{z}/m$, and it is normalized to the thermal speed $v_{t}$. Equation 19 is satisfied if each term of the series balances such that

$$\frac{dg_{n}}{da_{z}}H_{n}=g_{n}\frac{dH_{n}}{dv_{z}}.$$

(20)

Now, the Hermite polynomials are the orthogonal family with the property $H_{n}^{\prime}(x)=2nH_{n-1}(x)$. Substitution of the Hermite property into Eq. 20 shows that the coefficient functions $g_{n}$ must satisfy the recurrence relation

$$g_{n+1}=\frac{1}{2(n+1)}\frac{dg_{n}}{da_{z}}.$$

(21)

Iteratively substituting Eq. 21 into Eq. 18 gives the series

$$f=n_{0}\Big{(}\frac{\beta}{2\pi}\Big{)}^{3/2}\Big{(}\sum_{n=0}^{\infty}\frac{1}{2^{n}n!}\frac{d^{n}g_{0}}{da_{z}^{n}}H_{n}(v_{z})\Big{)}e^{-\beta H}$$

(22)

with $g_{0}=g_{0}(a_{z})$ the first Hermite coefficient. The particle density and current density follow as

	$\displaystyle n_{\alpha}(\varphi,a_{z})$	$\displaystyle=$	$\displaystyle n_{0}g_{0\alpha}(a_{z})e^{-\beta_{\alpha}q_{\alpha}\varphi},$		(23a)
	$\displaystyle j_{z,\alpha}$	$\displaystyle=$	$\displaystyle q_{\alpha}n_{0}v_{t\alpha}\frac{dg_{0,\alpha}}{da_{z}}e^{-\beta_{\alpha}q_{\alpha}\varphi}$		(23b)

because $H_{0,1}(z)=1,2z$ are orthogonal to all other $H_{n}$ over the Gaussian weighting ($e^{-\beta H}$). Equation 22 solves the kinetic equation, and its moments consistently solve the Vlasov-Maxwell system provided it satisfies the nonlinear Poisson equations

	$\displaystyle\nabla^{2}\varphi$	$\displaystyle=$	$\displaystyle\varepsilon_{0}^{-1}\sum_{\alpha}q_{\alpha}n_{0\alpha}g_{0\alpha}(a_{z})e^{-\beta_{\alpha}q_{\alpha}\varphi},$		(24a)
	$\displaystyle\nabla^{2}A_{z}$	$\displaystyle=$	$\displaystyle\mu_{0}\sum_{\alpha}q_{\alpha}n_{0\alpha}v_{t\alpha}\frac{dg_{0\alpha}}{da_{z}}e^{-\beta_{\alpha}q_{\alpha}\varphi}.$		(24b)

Thus, the function $g_{0}=g_{0}(a_{z})$ is a generating function for the equilibrium, determining density, current density, and distribution function in the laboratory frame through its derivatives. Recalling that axial magnetic vector potential is magnetic flux per unit length ($\psi_{\theta}^{\prime}\equiv\int_{0}^{r}B_{\theta}(r^{\prime})dr^{\prime}=-A_{z}$), this means that the generating function and all of its derivatives (the fluid moments) are flux functions. Notably, the toroidal flux function in kinetic tokamak equilibrium can be treated analogously through an azimuthal momentum ($P_{\theta}$) expansion Kaltsas et al. (2024). Intuitively, the generating function of a distribution canonical in the energy is expressed in Hermite series due to the properties of the Weierstrass transform, for which $e^{-\beta H}$ is the kernel Allanson et al. (2016). We now explore a few paradigmatic equilibria in terms of simple generating functions.

The kinetic equilibrium of radially uniform axial velocity is the isothermal Bennett pinch, whose distribution is semi-canonical (namely, a shifted Maxwellian). The function $g_{0}(a_{z})$ is exponential in $a_{z}$,

$$g_{0s}(a_{z})=\exp(2m_{s}\beta_{s}u_{0s}a_{z}).$$

(25)

Repeated differentiation of Eq. 25 and application of the Hermite generating function,

	$\displaystyle\frac{d^{n}g_{0s}}{da_{z}^{n}}$	$\displaystyle=$	$\displaystyle\Big{(}\frac{2u_{0}}{v_{ts}}\Big{)}^{n}g_{0}(a_{z}),$		(26a)
	$\displaystyle\sum_{n=0}^{\infty}H_{n}(v_{z})\frac{u_{0}^{n}}{n!}$	$\displaystyle=$	$\displaystyle\exp(2u_{0}v_{z}-u_{0}^{2}),$		(26b)

arrives at the shifted Maxwellian distribution,

$$f_{s}=n_{0s}\Big{(}\frac{\beta_{s}}{2\pi}\Big{)}^{3/2}e^{-\beta_{s}q_{s}(\varphi-2u_{0s}A_{z})}e^{-\beta_{s}m(v_{x}^{2}+v_{y}^{2}+(v_{z}-u_{0})^{2})/2}.$$

(27)

Both the Bennett pinch (cylindrical) and Harris sheet (Cartesian) are generated by Eq. 25 as they satisfy (under quasineutrality in the center-of-charge frame) the Liouville equation

$$\nabla^{2}A_{z}=-Ke^{2A_{z}},$$

(28)

with $K$ the curvature constant Ceccherini et al. (2005). The cylindrical solution to Eq. 28 is the Bennett pinch’s vector potential, given by

$$A_{z}=-A_{0}\ln(1+(r/r_{p})^{2})$$

(29)

where $r_{p}$ is the characteristic pinch radius and $A_{0}=\frac{\mu_{0}}{4\pi}I_{\infty}$ is the self-magnetic flux per unit length ($I_{\infty}$ is the total current enclosed to infinity). Regarding the generating function, it is useful that the following combination is unity,

$$2q_{s}\beta_{s}u_{0}A_{0}=\frac{\mu_{0}}{4\pi}\frac{N_{s}q_{s}^{2}}{m_{s}}\Big{(}\frac{u_{0}}{v_{ts}}\Big{)}^{2}=1$$

(30)

for each species $s$, where $N_{s}$ is the linear density (particles per unit length) of species $s$. The first factor is the species Budker parameter $\mathfrak{B}_{s}$. For the Bennett pinch, the inverse Budker parameter is the drift Mach number for that species Crews, Meier, and Shumlak (2024),

	$\displaystyle\mathfrak{B}_{s}$	$\displaystyle\equiv$	$\displaystyle\frac{q_{s}^{2}}{4\pi\epsilon_{0}}\frac{N_{s}}{m_{s}c^{2}},$		(31a)
	$\displaystyle\mathfrak{B}_{s}^{-1}$	$\displaystyle=$	$\displaystyle\Big{(}\frac{u_{0}}{v_{ts}}\Big{)}^{2}.$		(31b)

From Eq. 30, it follows from the generating function that the Bennett pinch particle and current densities are given by Meier and Shumlak (2021)

	$\displaystyle n_{s}(r)$	$\displaystyle=$	$\displaystyle\frac{N_{s}}{\pi r_{p}^{2}}(1+(r/r_{p})^{2})^{-2},$		(32a)
	$\displaystyle j_{z}(r)$	$\displaystyle=$	$\displaystyle\frac{I_{\infty}}{\pi r_{p}^{2}}(1+(r/r_{p})^{2})^{-2}.$		(32b)

In the isothermal Bennett kinetic equilibrium, both ions and electrons have uniform flow velocities, so their relative velocity is uniform too. As long as the relative velocity of two species is radially uniform then the magnetic potential will satisfy Eq. 28 and will be of the Bennett form, Eq. 29. The equilibrium has the important property that there is no radial electric field in the center-of-charge frame Gratreau (1978) where $u_{i}=-u_{e}$ (assuming $Z_{i}=1$).

The generator of the next simplest equilibrium is quadratic in the flux, namely $g_{0}=\exp(c_{0}a_{z}+c_{1}a_{z}^{2})$ for constants $c_{0}$, $c_{1}$. This equilibrium consists of a two-temperature bi-Maxwellian distribution with an agyrotropic pressure anisotropy supporting a radially sheared axial flow Del Sarto, Pegoraro, and Califano (2016). Anisotropic collisionless equilibria with sheared flows were explored by Mahajan and Hazeltine Mahajan and Hazeltine (2000), whose discussion we summarize here in the context of the Z pinch.

This sheared anisotropic equilibrium is the simplest kinetic equilibrium displaying the relationship between shear and anisotropy, enabling an understanding of how these two phenomena interact in the collisionless regime, which is relevant to Z pinches approaching fusion conditions. All sheared-flow collisionless kinetic equilibria are “inviscid” because they represent a stationary flow profile, which is a notable property. However, as will be discussed later, the collisionless magnetized dynamics should still be considered highly viscous in the sense that initially out-of-equilibrium flow profiles will quickly find a kinetic equilibrium supporting a steady flow.

Proceeding by direct calculation from an ansatz distribution is the simplest way to show existence of the equilibrium and determine its most fundamental properties. Suppose the distribution function is a two-temperature (bi-)Maxwellian whose first velocity moment $\langle v_{z}\rangle=u_{z}(r)$ is radially sheared (with species subscripts suppressed),

$$f(r,\bm{v})=\frac{\beta_{\perp}}{2\pi}\sqrt{\frac{\beta_{z}}{2\pi}}n(r)e^{-\beta_{\perp}\frac{m}{2}(v_{r}^{2}+v_{\theta}^{2})}e^{-\beta_{z}\frac{m}{2}(v_{z}-u_{z}(r))^{2}}.$$

(33)

Substituting Eq. 33 into the Vlasov equation leads to

$$\partial_{t}f+v_{r}\Big{(}\frac{n^{\prime}}{n}+\beta_{\perp}q\varphi^{\prime}+m\beta_{\perp}u_{z}(r)\omega_{c}\Big{)}f+m\beta_{z}v_{r}(v_{z}-u_{z})A_{rz}f=0$$

(34)

where $\omega_{c}=qB/m$, a shear stress drive term $A_{rz}\equiv u_{z}^{\prime}-\alpha_{rz}\omega_{c}$, and the distribution anisotropy is $\alpha_{rz}\equiv 1-T_{z}/T_{r}$ (for $T_{r}=(k\beta_{\perp})^{-1}$). Taking the first radial moment of Eq. 34 yields

$$\partial_{t}(n\langle v_{r}\rangle)=-n^{\prime}kT_{r}-qn\varphi^{\prime}-qnu_{z}B_{\theta}$$

(35)

giving the equilibrium radial force balance $\nabla p=qn(\bm{E}+\bm{v}\times\bm{B})$ because $\langle v_{r}\rangle=0$. The distribution function supports this force balance in equilibrium (that is, with $\partial_{t}f=0$) if the driver $A_{rz}=0$, namely when the anisotropy and shear profile are related by

$$\frac{du_{z}}{dr}=\alpha_{rz}\omega_{c}.$$

(36)

Having assumed a radially constant anisotropy $\alpha_{rz}$, the axial species velocity is then given by $u_{z}=u_{0}-\alpha_{rz}a_{z}$, meaning that the equilibrium velocity profile is a linear flux function.

The bi-Maxwellian distribution is arrived at from the generating function using the completeness relation for Hermite polynomials Wiener (1988), as the derivatives of the generating function are themselves Hermite polynomials. Applying the completeness relation for the sum $\sum_{n}\frac{1}{n!4^{n}}H_{n}(u_{0}+\alpha_{rz}a_{z})H_{n}(v_{z})$ results in a bi-Maxwellian. The algebra is simpler to proceed as above, and as conducted in Mahajan and Hazeltine Mahajan and Hazeltine (2000), who ultimately find the bi-Maxwellian equilibrium for the sheared Bennett equilibrium (following section) to be

$$f(P_{z},H)=\exp(-\beta(H+u_{0}P_{z}+\gamma P_{z}^{2}/2m))$$

(37)

where $\gamma$ is a constant related to the anisotropy. Equation 37 is a nice closed form result for the distribution function whose sheared flow is a linear flux function. Higher-order dependence of the velocity on the flux function is explored in Section IV.4, but the series appear not to sum cleanly into exponentials as in Eq. 37.

Collisionless kinetic equilibrium of a sheared-flow self-pinched plasma is not necessarily a bi-Maxwellian distribution, which occurs for an isothermal pinch of radially constant $\alpha_{rz}$ temperature anisotropy in which flow velocity is a linear function of magnetic flux. Experimental and numerical simulation dynamics likely display both gyrotropic and agyrotropic anisotropies, due both to flow shear and to the asymmetric responses of cyclotron and betatron trajectories to changes in flux; hence general kinetic equilibria should display both axial and azimuthal temperature gradients. Non-isothermal kinetic equilibria are reserved to a future work.

But even considering the isothermal solutions, velocity need not be a linear flux function. Equilibria can be found as any polynomial function of flux. This section studies more general kinetic equilibria such that the axial velocity is some polynomial of the flux function, and studies the characteristic features of these solutions. First observe that the axial velocity is given by a logarithmic derivative of the generating function,

$$\langle v_{z}\rangle=\frac{1}{g_{0}}\frac{dg_{0}}{da_{z}}=\frac{d\ln g_{0}}{da_{z}}.$$

(39)

Thus, each exponential polynomial $g_{0}(a_{z})$ is associated with a velocity profile as a polynomial in the magnetic potential,

$$g_{0}(a_{z})=\exp\Big{(}\sum_{n=0}^{\infty}\frac{c_{n}}{n+1}a_{z}^{n+1}\Big{)}\implies\langle v_{z}\rangle=\sum_{n=0}^{\infty}c_{n}a_{z}^{n}$$

(40)

where $c_{n}$ are some coefficients. Equation 40 is a formal power series expansion of axial velocity about the pinch axis. With the velocity expressed as such a flux function (Eq. 40), we then seek the Hermite series coefficients to determine the equilibrium distribution function.

Substitution of Eq. 40 into Eq. 22 produces the distribution $f=\widetilde{f}(P_{z})e^{-\beta H}$ by the associated function of $P_{z}$, whose Hermite series coefficients are the derivatives of $g_{0}(a_{z})$. Let

$$h(a_{z})\equiv\sum_{n=0}^{\infty}\frac{c_{n}}{n+1}a_{z}^{n+1}$$

(41)

so that $g_{0}=\exp(h(a_{z}))$. Application of Faà di Bruno’s formula to Eq. 40 yields the Hermite coefficients as

$$\frac{d^{n}g_{0}}{da_{z}^{n}}=e^{h(a_{z})}B_{n}(h^{\prime},\cdots,h^{(n)})$$

(42)

where $B_{n}$ is the n’th complete Bell polynomial Bell (1934). Thus, the semi-canonical distribution function associated with the postulated velocity profile of Eq. 40 is

$$f=e^{-\beta H+h(a_{z})}\sum_{n=0}^{\infty}\frac{1}{2^{n}n!}B_{n}(h^{\prime},\cdots,h^{(n)})H_{n}(v_{z}).$$

(43)

It is intuitive to first consider the forward process, namely viscous damping of sheared flows. To begin, we extend the equilibrium analysis of Mahajan and Hazeltine by considering the evolution of an agryotropic distribution initially out-of-equilibrium, that is, when the equilibrium anisotropy condition $A_{rz}=0$ is not met. Assuming initial radial force balance and introducing a collision operator $C[f]$, Eq. 34 can be written as

$$\frac{\partial f}{\partial t}\Big{|}_{t=t_{0}}=C[f]-m\beta_{z}v_{r}(v_{z}-u_{d})A_{rz}f$$

(48)

instantaneously at the initial time $t=t_{0}$. Given the bi-Maxwellian form of $f$, the collisionless term on the RHS of Eq. 48 may be rewritten as derivatives of $f$, giving

$$\frac{\partial f}{\partial t}\Big{|}_{t=t_{0}}=C[f]+\frac{\partial^{2}}{\partial v_{r}\partial v_{z}}(v_{th,r}^{2}A_{rz}f)$$

(49)

and revealing it to be a diffusive operator, where $v_{th,r}^{2}=kT_{r}/m$ is the radial thermal velocity. The mixed-derivatives term of Eq. 49 satisfies many of the features of a collision operator, such as conservation of particles, energy, and momentum, but drives shear stresses by not conserving the off-diagonal pressure tensor components $\Pi_{rz}=\int v_{r}v_{z}fd\bm{v}$.

Adopting a simple collision model such as the Lenard-Bernstein form shows the mixed derivatives to introduce off-diagonal terms to the velocity space diffusion matrix. The Fokker-Planck equation then has the paradigmatic form

	$\displaystyle\frac{df}{dt}$	$\displaystyle=$	$\displaystyle-\nu\nabla_{v}\cdot(\bm{v}f)+v_{t}^{2}\nabla_{v}\cdot(\mathcal{D}\nabla_{v}f),$		(50a)
	$\displaystyle\mathcal{D}$	$\displaystyle=$	$\displaystyle\begin{bmatrix}\nu&A_{rz}/2\\ A_{rz}/2&\nu\end{bmatrix}$		(50b)

where $\nu$ is the relaxation rate and $A_{rz}$ the shear stress drive frequency. The off-diagonal terms produce shearing in the velocity space (demonstrated in Fig. 5) as they are symmetric.