Experimental and computational study of phase space dynamics in strongly coupled plasmas with steep density gradients Open Access

Hydrodynamic models accurately predict plasma transport when the mean free path is shorter than the plasma scale length.^1–4 This capability is critical for experiment planning, data analysis, plasma characterization, and parameter optimization.⁵ Radiation-hydrodynamic codes like LILAC,⁶ HYDRA,⁷ and xRAGE⁸ have been shown to successfully capture the temperature, density, and neutron yield of experiments at NIF and OMEGA. They nicely model the influence of target surface defects, voids, and structural imperfections on ICF neutron yield in the hydrodynamic limit.^9,10 Hydrodynamic codes are widely used because they are much faster than kinetic codes and incorporate radiation effects much more easily.¹¹ 

However, many plasmas start far from the hydrodynamic limit.¹² Kinetic effects arise when the velocity distributions are strongly non-Maxwellian, which is the inevitable outcome of nearly every laser–matter interaction, including laser direct-drive experiments. In some cases, when experiment and simulation are directly compared, hydrodynamic codes do not reproduce the experiments even though the hydrodynamic limit is reached,¹³ suggesting the presence of kinetic effects outside the hydrodynamic approximation. In models of laser direct-drive experiments at OMEGA, for example, understanding kinetic effects is critical.⁶ 

In a series of papers, Killian's group at Rice University studied the kinetic-to-hydrodynamic transition using ultracold neutral plasmas (UNPs).^14–16 By placing a series of obstructions in the ionizing laser beam, those authors created a UNP with “sculpted” initial conditions.¹⁷ They found that UNPs with sharp and deep density gradients produced non-thermal ion velocity distributions. They explored the evolution of a UNP with a gap in the middle. Their gap had sharp edges, leading to ion jetting across the gap. They compared their results to a hydrodynamic code and found that shallower gradients, smaller gaps, and lower initial electron temperatures led to more hydrodynamic behavior.

In this paper, we take a closer look at the kinetic-to-hydrodynamic evolution of a UNP with steep density gradients. We create a strongly coupled calcium UNP with steep density gradients. Using precision laser spectroscopy, we determine the plasma ion phase space distribution function. We compare these distributions with predictions from a kinetic code that uses two different collision operators. We find that our simulations generally reproduce plasma behavior, with significant differences in the approach to equilibrium. The kinetic codes rely primarily on low-order moments of the distribution functions and do not capture all of the observed phase-space dynamics. This platform enables precision tests of kinetic theories and opens new possibilities for studying plasma stopping power and flow-induced instabilities in strongly coupled systems.

The Yukawa model does not capture all of the expected physics in HEDPs or UNPs. Electron-ion thermalization,^22–33 for example, is not explicitly included in the Yukawa model. Quantum effects,^34,35 such as bound states, three-body recombination, and quantum degeneracy, are also not included. HEDPs and UNPs share similar physics in these non-Yukawa aspects as well.

The Yukawa model is commonly used to describe ion–ion interactions in both weakly and strongly coupled plasmas.^19,36 Examples of strongly-coupled plasmas with $Γ>1$ include white dwarf stars,³⁸ the cores of Jovian planets,^39,40 some stages of laser-driven plasma experiments^6,41 and X- and Z-pinch plasmas,⁴² dusty plasmas,^43,44 quark-gluon plasmas,^45,46 dusty plasmas consisting of highly charged dust particles,^47–50 non-neutral plasmas,^51–53 ions in UNPs^54–56 and other systems.

Comparing model predictions with experimental data provides an important check of plasma theories in the strongly coupled plasma regime. The range of strong coupling parameter values accessible to UNPs⁵⁷ is in the range of $Γ≈2$⁠. Using sculpted plasmas,¹⁷ it is possible to create kinetic conditions that are far from equilibrium. When $Γ≪1$⁠, collisions are described using small-angle scattering that leads to various prescriptions for Coulomb logarithms.^28,58 The Coulomb logarithm description can be extended up to $Γ∼1$ using a numerically computed binary cross sections.²¹ Ultimately, for $Γ≫1$⁠, one needs to include collective phenomena.⁵⁹ 

Strong coupling also introduces new behaviors in the plasma that arise from spatial order. For example, the scaled viscosity displays a minimum near $Γ∼5$ to 10, increasing toward very large values of the strong coupling parameter.²¹ As $Γ$ increases toward 1, many-body hard collisions become important. The characteristic length scale changes from the Debye length to the inter-particle spacing. As $Γ$ increases beyond $Γ=10$ to greater coupling, viscosity increases again. Isolating transport processes in HEDPs is difficult because many different processes compete.^60,61 Ideally, a precision experiment which enables isolation of one or two simple processes would make it possible to evaluate theoretical treatments. UNPs provide such an opportunity in a range of $Γ$ and $κ$ relevant to HEDP experiments.⁶² 

Three excellent reviews describe how UNPs are generated, diagnosed, and modeled.^54,55,63 UNPs are strongly coupled, non-degenerate, quasi-homogeneous, quasi-steady-state, $Z=1$ plasmas with accurately known initial conditions.⁶⁴ They provide an idealized platform for measuring plasma transport properties,^18,64–70 and can be used to simulate transport in HEDPs.^{16,18,57,71–77} Some studies in our lab include dynamics of strongly coupled neutral Coulomb systems,^78–81 exploring the HEDP/UNP crossover,^18,71 ion transport properties of binary ion mixtures,^64,65,72 and the density evolution of partially magnetized plasmas.^66,82,83

The peak plasma density is in the range of $n0=108–1010 cm−3$⁠, the initial rms size is $σ0∼1 mm$⁠, and the characteristic expansion time is $τe=[mσ02/(kBTe)]1/2$⁠. The resulting plasma is just inside the strongly-coupled plasma regime with $Γ=2$⁠, although a range of $Γ=0.1–11$ is possible.⁵⁷ 

We analyze these plasmas using laser-induced fluorescence (LIF).^64,72,88 A probe laser beam tuned to the Ca⁺ 397 nm resonance [Fig. 1(b)] is shaped to be 1 cm wide in the y-direction and 1 mm thick in the x-direction, and it is centered on the plasma [Fig. 1(c)]. Ion fluorescence in the illuminated portion of the plasma is imaged onto a gated ICCD camera, allowing us to record a cross-Sec. of the ion density distribution near the center of the plasma. The ICCD camera is triggered after a delay time of 0–20  $μ$s and captures fluorescence during a short time window, typically between 25 and 200 ns. Identical plasmas are generated 10 times per second. To improve the signal-to-noise ratio, we repeatedly measure ion fluorescence on 50 sequentially created plasmas with identical measurement conditions.

We adjust the frequency of the probe laser beam in $Δf=20$ MHz steps that span the entire velocity distribution. For ions moving at the appropriate velocity, the probe laser beam is Doppler-shifted into resonance, $v=λΔf$⁠, where $λ$ is the atomic transition wavelength. The ICCD camera then records where ions at that velocity are located in the plasma. By analyzing the images at different probe laser frequencies, we are able to map out a spatially resolved ion distribution function at the chosen delay time. The Doppler shift provides information about the velocity component that is parallel to the probe laser beam propagation direction. The camera gives the location of those ions.

As shown in Fig. 1(d), we integrate the fluorescence image over a range of $±0.5$ mm in y to obtain a portion of the distribution function as a function of z at that ion velocity and time, $f(z,vn,tn)$⁠, where $vn$ is the ion velocity Doppler-shifted into resonance with the probe laser beam and $tn$ is the delay time at which the image was taken. This y-integrated image becomes one row in the phase-space image of the plasma, showing where all the ions of that particular velocity are located in space. We repeat this measurement for a wide range of probe laser beam frequency detunings to generate the image shown in Fig. 1(e).

The image in Fig. 1(e) represents the convolution of the true distribution function with the Lorentzian linewidth of the Ca⁺ 397 nm transition. The convolution appears in both velocity v and position z due to hydrodynamic flow in the z direction. Deconvolving the images is challenging due to the non-trivial hydrodynamic velocity distribution in sculpted plasmas.

We generate sculpted plasmas¹⁷ by spatially structuring the ionizing laser beam. We place a wire in the ionizing laser, directly in front of a lens pair. The wire is 1:1 imaged into the plasma and the ionizing laser beam is nearly collimated at the plasma⁸² when the lenses are identical, spaced by 2f, and the wire-to-plasma distance is 4f. For this work, we use a pair of 200-mm focal length UV-achromatic lenses to image the wire onto the plasma. The plasma forms as the optical “negative” of the wire. A wire is imaged as a sheet at the plasma, creating two hemispherical plasmas separated by a sharp-edged gap. For the data shown in this paper, our wire diameter and the plasma gap are both 0.28 mm.

In our experiment, we collect fluorescence images when the probe laser beam propagates both parallel and perpendicular to the direction of the density gradient. In this way, we generate information about $f(y,vy,t)$ and $f(z,vz,t)$⁠. Due to the symmetry of our experiment, $f(x,vx,t)=f(y,vy,t)$⁠.

Kinetic models track particle distributions in phase space. They accommodate non-equilibrium flow features such as e.g., bump-on-tail interactions, two stream instabilities. The cost is that they are more computationally expensive than hydrodynamic models due to their higher dimensionality.

We use two collision operators in our kinetic code: the Bhatnagar–Gross–Krook (BGK)^89–91 and the Lenard–Bernstein (LB)⁹² models.

For both collision models, the collision times $τ$ are a free parameter that can be tuned to match specific properties, e.g., temperature or momentum relaxation rates. In this work, the Ion–ion $τ$ are chosen to match temperature relaxation using the Stanton–Murillo Transport (SMT) cross sections,²¹ which are valid into moderately coupled plasma regimes. These collision frequencies, as well as the derivations of the BGK and LB forms of the collision operator themselves, implicitly rely on the assumption that the distribution functions are near Maxwellian; however, they are favored for their simplicity relative to the much more computationally expensive Boltzmann or Fokker–Planck models.

Our BGK model is a finite-volume, 1D-3V code capable of modeling multiple kinetic species.⁹³ We use $Nx=496$ spatial grid points, each 5.2  $μ$m wide, and $Nv=1003$ velocity grid points with a resolution of 2.8 m/s. The time step is 1 ns. The ions are the single kinetic species. The simulation is initialized with the experimentally measured $f(y,vy,0)$ and $f(z,vz,0)$⁠.

The ion distribution function is evolved using Eq. (1) with collision operators in either Eq. (6) or Eq. (8). We solve the kinetic equation in conjunction with Eqs. (10), (11), and (12). The nonlinear equation (11) is solved using a standard Newton–Raphson method.⁹³ 

Experimental and simulated distribution function results are shown in Fig. 2. The z-component of the distribution function, $f(z,vz,t)$ is represented using a logarithmic false-color plot. An imperfection in the optical imaging system causes additional fluorescence to appear in the gap between the plasmas. This is most clearly evident in the first and second panels of Fig. 2(a), for 0.1 and 1.1  $μ$s. We have verified that the ion density in the gap is less than 1% of the peak plasma density using an independent measurement. For these experimental data, the peak density is $1.2×109 cm−3$⁠. The initial electron temperature is $Te=96 K$⁠. The ion fluorescence is integrated over a 200 ns window centered on the time displayed in the plots.

The experimental data are compared with the simulation in Figs. 2(b) and 2(c). Both the BGK and LB simulations reproduce the ion jetting and counterpropagating plasma flow at 1.1 and 3.1  $μ$s. The thermalization rate in both simulations is somewhat slower than the experiment. Because the experimental inputs are well known, the simulations are tightly constrained. Referring to Eq. (1), only the terms on the right-hand side are in question. For kinetic codes, the parameters are the electron model, the collision frequency, and the form of the collision operator.

For the electron model, we use a Poisson–Boltzmann equation. The electrons are modeled as a fluid. We tried running the simulation treating the electrons as a second kinetic species. However, the large ion-to-electron mass ratio caused the simulation to run unacceptably slowly.

We vary the form of the collision operator using both BGK and LB collision operators, Eqs. (6) and (8). The primary difference is that the LB operator adds a velocity diffusion term. From the data in Fig. 2, this diffusion term speeds the approach to equilibrium, but in a way that doesn't precisely match the experimental data.

Varying the collision frequency also does not improve the agreement between the simulation and experiment. In both the BGK and LB simulation, decreasing the collision frequency decreases the thermalization rate. On the other hand, increasing the collision frequency dramatically decreases the ion jetting into the gap. Ions accelerate out from the two sides of the gap, collide in the middle without overlapping, and generate a broad Maxwellian velocity distribution. In addition, the increased collision frequency increases the pressure in the gap, restricting the plasma expansion at the gap edge.

The agreement between the experiment and simulations at 0.1 and 1.1  $μ$s in Fig. 2 validates the electric field model described in Sec. IV B. At these times, when the relative ion velocities in the gap are high, we expect the ballistic motion and electric fields to dominate the evolution of the ion distribution making the collision terms negligible.

At even earlier times, the electrons may not be Maxwellian. Nonthermal electrons could race into the gap ahead of the ions, contributing to a greater electric field than predicted by our model.

Both the experiment and the simulation reveal a linear dependence of the ion acceleration on the initial electron temperature. The simulations underpredict the initial acceleration, as expected. However, this underprediction is inconsequential in the overall plasma evolution, because the electrons quickly thermalize. The electron density model in Eq. (10) is a good approximation for most of the plasma evolution.

The data in Fig. 4 provide insights into the details of the kinetic simulations. In the longitudinal direction, the simulations match each other and the experiment for the second moment of the velocity distribution. In the transverse direction, the simulations overpredict the broadening of the second moment. This is likely due to the isotropic nature of the collision model's treatment of post-collisional velocities. The BGK model assumes a spherical Gaussian post-collision velocity distribution. An improved model that used an ellipsoidal distribution would better capture scattering in the direction of particle motion.^94,95 Similarly for LB, which broadens the fastest in Fig. 4, an improved model would allow for non-symmetric diffusion tensor.⁹⁶ 

In the experiment and both simulations, we verify that the transverse distribution functions are nearly Maxwellian at all times. However, the longitudinal distributions are more complicated. In the experiment, the bimodal nature of the longitudinal velocity distribution in the center of the plasma has relaxed by 5  $μ$s [see Fig. 2(a)] and reached a local equilibrium by 7  $μ$s (see Fig. 4). A similar situation is true for the LB simulation. However, the BGK simulation has not fully relaxed by 7  $μ$s.

One of the aims of this study is to probe the validity of kinetic models compared to a simple, well-characterized experiment, one in which the initial ion density, electron temperature, spatial density profile, and charge state are known. In both the BGK and LB simulations, the approach to thermalization is much slower than in the experiment. This is evident in both Figs. 2 and 4. This discrepancy likely stems from fundamental limitations in both kinetic models. The BGK model approximates the collision operator by relaxing the distribution function toward a local Maxwellian at a single characteristic rate. It considers only the first three velocity moments (density, momentum, and energy) and assumes these moments adequately capture the essential physics of thermalization. For example, the second moment of the velocity distribution in the longitudinal direction matches the experiment in Fig. 4. However, this simplified treatment may not properly account for the complex velocity-space structures that develop in strongly-coupled plasmas with steep gradients.

The Lenard–Bernstein operator introduces velocity-space diffusion and friction, potentially capturing some additional physics missing from BGK. However, it uses only the same low-order moments to define these terms. It also relies on a factorized collision frequency that may not accurately represent the true collision processes in our system. The velocity-space derivatives in LB help describe the smoothing of sharp features in the distribution function. In our hands, the LB model over predicts the relaxation of counterpropagating ion flow.

Our electric field model, based on Poisson–Boltzmann in Eq. (12), treats electrons as a fluid rather than a kinetic species. While this assumption significantly improves computational efficiency, it may oversimplify the electron response to density gradients. The fluid treatment cannot capture non-Maxwellian features in the electron distribution function that could affect ion thermalization through modified screening and electric field fluctuations.

Several factors could contribute to observed discrepancies between the simulations and the experiment and also between the two simulations.

• The assumption of a single relaxation time in BGK may be too simplistic, as different moments of the distribution function might relax at different rates in strongly coupled systems.

• Both models account for strong coupling effects that could enhance thermalization through collective modes only implicitly, through the SMT collision cross-sections.

• The use of only three moments in BGK might miss important higher order effects that become relevant when density gradients are steep.

• Both kinetic models derive from a perturbative expansion of the Boltzmann equation. BGK is derived using an expansion around an equilibrium distribution, plus further approximations. LB derives from the Fokker–Planck equation, which approximates the collision physics assuming small angle scattering. These perturbative approaches may not be valid for our highly non-equilibrium plasma in the gap.

• The factorized collision frequencies in both models may not capture the full velocity dependence of the actual collision processes and collision cross sections.

• The fluid treatment of electrons may miss important kinetic effects in the electron response that could modify the electric field and affect ion thermalization.

An extended BGK model (EBGK) that includes additional moments and a more sophisticated treatment of the collision frequency might better match the experimental results. Additionally, treating electrons kinetically, while computationally intensive, could provide insight into whether electron kinetic effects significantly influence the thermalization rate. However, the observed discrepancies suggest that new theoretical approaches may be needed to fully describe thermalization in strongly coupled plasmas with steep gradients.

With a better model, future work could focus on more complex interface structures, such as a sharply pointed triangular gap, or systems in which the plasma density varied significantly on either side of the gap. Adding the possibility of binary ion mixtures further enriches the HEDP applicability of future studies.

This work was supported in part by the National Science Foundation (Grant Nos. NSF-2009999 and NSF-2108505). This work was supported in part by the U.S. Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

The authors have no conflicts to disclose.

Scott Bergeson: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Matthew Schlitters: Data curation (equal); Investigation (equal). Matthew Miller: Data curation (equal); Investigation (equal). Ben Farley: Data curation (equal); Investigation (equal). Devin Sieverts: Data curation (equal); Formal analysis (equal); Investigation (equal). Michael S. Murillo: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Jeffrey R. Haack: Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.