Wave steepening and shock formation in ultracold neutral plasmas

A shockwave in a plasma is an abrupt jump in pressure, temperature, and density across a narrow front that propagates through the plasma.¹ Shocks can form when a driver moves through a plasma or a population of particles moves through a background plasma² at relative speeds that exceed the speed of sound. Ordinary waves can develop into shockwaves through wave steepening, in which nonlinear propagation effects lead to a steepening of leading density gradients.³ Plasma shockwaves are observed in many settings, such as laser-induced plasmas,⁴ supernova remnants,⁵ and the solar wind.⁶ 

In this paper, we present experimental observations of wave steepening and signatures of shock formation during the expansion of ultracold neutral plasmas (UNPs). UNPs are created by photoionizing laser cooled atoms^7–9 or cold molecules in a supersonic expansion¹⁰ near the ionization threshold. They have been used to study many different plasma phenomena, such as plasma waves,^11–14 streaming plasmas,¹⁵ three-body recombination,^16,17 instabilities,¹⁸ and equilibration and transport in the strongly coupled regime.^19–21 The expansion into surrounding vacuum of a UNP with an initial Gaussian density distribution has been studied extensively,^11,22–28 including expansion in a uniform magnetic field.^29,30 Recently, expansion of UNPs with an initial density distribution that decays exponentially with distance from the plasma center was characterized experimentally in a field-free configuration³¹ and in a quadrupole magnetic field.³² Molecular dynamics simulations were presented for each case.^33,34

Since UNPs were first created,⁷ the development of shockwaves in this extreme environment has been of interest. Early fluid²² and hybrid molecular dynamics²³ models indicated shockwave formation in the low-density outer regions of expanding plasmas with an initial Gaussian density distribution. Experiments failed to detect such features within the limit of experimental sensitivity.²⁷ It has also been predicted that shocks develop in UNPs with non-Maxwellian electrons.³⁵ More recent molecular dynamics simulations of a Gaussian, quasi-neutral ultracold plasma showed that significant escape of electrons from the outer regions of the plasma and resulting non-neutrality can give rise to shock-like features.³³ 

Early numerical work by Sack and Schamel³⁶ in a very different context than UNPs, showed that expansion into vacuum of a plasma with an initial step-function density distribution yields wave-steepening and formation of a sharp density peak and shock structure. More recent work³⁷ shows that sharp fronts persist during the plasma expansion and points out the role of charge separation. Following this line of reasoning, the possibility of forming shocks in UNPs sculpted to have initial density distributions with steep density features deviating from an ideal Gaussian has been discussed in multiple contexts.^13–15,38 Experimental realizations produced streaming populations of UNPs, displaying a crossover from hydrodynamic^13,14,39 to kinetic behavior¹⁵ with increasing initial density gradients and resulting relative velocities, but signatures of shock formation were not observed. A hydrodynamic model of colliding ultracold neutral molecular plasmas also predicted the formation of shocks.⁴⁰ While not in the UNP context, shockwaves have been seen in the Coulomb explosion of a cold, pure-ion, non-neutral plasma created by photo-excitation well above the ionization threshold such that electrons escape the system.⁴¹ 

We observe that the shape of the initial density distribution of a UNP largely determines the character of the plasma expansion when magnetic fields are negligible.³¹ Plasmas with a Gaussian-shaped initial density distribution, which are well-studied experimentally, expand self-similarly²⁷ and do not display detectable signatures of wave-steepening and shock formation. As shown in this study, however, plasmas formed with a more sharply peaked distribution that decays exponentially in space^31,32 yield more complex dynamics, and after significant expansion, they display many of the traits that characterize the formation of shocks.

After the magnetic trap is loaded with about $108$ atoms, the cooling laser beams and the magnetic field are turned off. A 10 ns pulse of 322 nm light photoionizes the cloud of strontium atoms forming a plasma. The electron temperature is set by how far above the ionization threshold the 322-nm beam is detuned,⁸ as illustrated in Fig. 1. Experiments reported here were performed at electron temperatures (⁠ $Te=160$ K) high enough to avoid the effects of three body recombination, which increases at lower temperatures.¹¹ 

The early stages of UNP equilibration are described in Ref. 8. After photoionization, electrons equilibrate in the first few hundred nanoseconds. On the timescale of the ion plasma oscillation of about 1  $μ$s, the ions undergo disorder-induced heating (DIH)^46,47 to a temperature of $Ti≈e2/(12πϵ0akB)≈1$ K depending on the density, where e is the elementary charge, $ϵ0$ is the vacuum permittivity, $kB$ is the Boltzmann constant, and $a=[3/(4πni)]1/3$ is the Wigner–Seitz radius, and $ni$ is the ion density. Following intraspecies equilibration, the plasma expands into the surrounding vacuum. At an adjustable time after plasma creation, the Sr⁺ ions are imaged using laser-induced fluorescence (LIF) on the $2S1/2$- $2P1/2$ transition at 422 nm.⁴⁸ Ions are excited with a thin (1 mm) sheet of laser light illuminating the central plane of the plasma. LIF emitted perpendicular to the light sheet is imaged onto an intensified charge coupled device (CCD) camera using a 1:1 optical relay with a resolution of 0.1 mm.

Figure 2(a) shows a false-color image of the density profile of a UNP 0.3  $μ$s after plasma formation, which is close to the initial plasma density distribution. The image is of the plane illuminated by the LIF laser, passing through the center of the plasma along the axis of cylindrical symmetry of the initial atomic density distribution.

The steepest plasma gradients are along the x-axis [see Eq. (1)], producing the highest plasma acceleration along this axis and an inversion of the aspect ratio during the expansion as shown in the image taken at 12  $μ$s in Fig. 2(b). Sharp transitions in the density are visible at the edges of the plasma experiencing the fastest expansion. We focus on the evolution of the density along the x-axis to look for signatures of wave steepening and shock formation.

A transect of the initial density at $y≈0$ is shown in Fig. 3. To fit the data, the 2D exponential [Eq. (1)] is averaged over a width in the z-axis corresponding to the thickness of the LIF light sheet. This produces a slight rounding of the peak. The overall match to the data are good, although there is some excess plasma in the regions of approximately $1<|x|<3$ mm. This pedestal may contribute to wave steeping as it creates a region of lower $∇n/n$ and therefore lower initial acceleration velocity [Eq. (2)].

For an exponential density distribution, the hydrodynamic expression for acceleration [Eq. (2)] yields a step function, uniform on either side of the plasma, but with opposite signs. This would initially give rise to a velocity profile that is also a step function. The evolution is more complex than this simple picture because of initial deviations from a perfect exponential and the density distribution quickly evolves away from an exponential.

Figure 4 shows transects of the experimentally measured plasma density (a), bulk velocity (b), and ion temperature (c) along the x-axis at various times after photoionization, compared to fluid simulations produced by SPRUCE, a multipurpose fluid code under development in the Solar Physics Research Group at Rice University.⁴⁹ 

The density [Fig. 4(a)] becomes flat in the center within the first few microseconds of evolution. As time progresses, the size of the central region of flat plasma density increases as the overall plasma expands and the peak density decreases. The leading edge of the plasma expansion accelerates with time as electron thermal energy is converted into directed kinetic energy of the ion expansion.¹¹ 

The expansion for the exponential plasma is clearly not self-similar. It is instructive to contrast it with the expansion of a spherical Gaussian density distribution^27,31 given by $n=n0 exp(−r2/2σ02)$⁠, for which the hydrodynamic force [Eq. (2)] yields an acceleration that increases linearly from zero at the plasma center. This yields a hydrodynamic velocity of $v=γ (t)r$ that increases linearly with distance from the plasma center and initially increases with time, resulting in a self-similar expansion, $σ0→σ(t)$⁠. Both $γ(t)$ and $σ(t)$ can be expressed analytically.^8,27

For the exponential plasma, the acceleration is initially much larger near the center than for a Gaussian plasma of comparable size and electron temperature, but the acceleration of the Gaussian plasma will be greater at large radial distance. An exponential density is therefore susceptible to the development of a velocity distribution that displays extrema away from the plasma center, as shown in Fig. 4(b). For the earliest time point, the sharply peaked center of the plasma may experience significant non-neutrality in a region on the order of the Debye length, which is 23  $μ$ m for the peak density. This may create a local Coulomb explosion that also contributes to the acceleration of ions from the central region.³¹ 

The bulk velocity in the flat central region of the exponential plasma is well-approximated by the linear relationship $v=x/t$ after the first few microseconds, indicating that plasma in this region quickly evolves to a state of expansion with vanishing acceleration, as expected for a flat density [Eq. (2)]. Note that Fig. 4(b) shows the velocity scaled by the ion acoustic wave velocity, which will be discussed in detail as follows and is given for each trace in the caption of the figure. The largest observed peak expansion velocities are $∼200$ m/s at 7  $μ$s.

The peaks of the ion velocity for positive and negative x precisely correspond to the edges of the central density plateau. In the velocity profiles, it is evident that in the range $0.3−7.3 μ$ s, higher-velocity plasma is overtaking lower-velocity plasma in front of it yielding wave steepening. By $∼7 μ$ s, this creates well-defined high-density, high-velocity fronts that propagate through the low-density, low-velocity plasma in front of them, with a large jump in density and velocity across the front.

At $12 μ$ s, the velocity peaks and gradients are extremely sharp, with a length scale for the front (from 75% to 25% of peak density) of 0.4 mm, much smaller than the several-millimeter scale of the overall plasma. The jump in velocity is over 100 m/s. The spatial overlap of the sharp density and velocity features is clearly seen in the expanded views shown in Fig. 5.

The experimental and simulation results agree reasonably well. The simulation treats the plasma as an ideal quasi-neutral fluid in two dimensions with one continuity and momentum equation and separate electron and ion energy equations. We equate the ambipolar electric field force term to the electron pressure gradient through $eniE≈∇Pe$⁠. The fluid equations are solved on a uniform Cartesian grid using Barton's method for monotonic transport to compute transport derivatives, second-order central finite differences for all other derivatives, and a second-order Runge–Kutta time integration scheme.⁴⁹ Barton's method is a simple, yet robust, method of differentiation for discontinuous fluid problems that result in shock conditions without post-shock oscillations present in other methods.⁵⁰ Open boundary conditions are implemented by choosing a large simulation domain $|x|,|y|<1.2$ cm and limiting spatial derivatives at the boundary through the use of two ghost cells, which copy the fluid characteristics of the nearest interior cell.

The plasma density is initialized on the domain by interpolating the density distribution shown in Fig. 2 for distance from the center $|r|<5$ cm and filling all cells exterior to this radius with an exponential fit to the distribution in Fig. 2 superimposed on a small uniform background (⁠ $10−4$ of peak density) to avoid simulation instability at vanishing density. The initial electron temperature is set to 160 K to match the energy expected from the detuning of the ionizing laser above resonance. The ion temperature profile is set to the initial measured value, after disorder-induced heating.

The validity of this quasi-kinetic treatment of the electrons, where the electrons are globally hydrodynamic but locally kinetic, has been demonstrated for UNPs with a variety of initial density distributions such as a spherically symmetric Gaussian density distribution imprinted with planar ion holes and periodic perturbations that excite ion acoustic waves.⁴⁹ The simulations in two dimensions accurately capture the key elements of the plasma expansion dynamics. The pressure gradients that drive plasma acceleration are identical in each case as is the resulting expansion-induced cooling for thermal energy density, $e=P/(γ−1)$⁠, where P is the plasma pressure, which is dominated by the electrons. $γ=2$ is the two-dimensional adiabatic index used in the simulations. Likewise, the plasma sound speed is identical in each case due to its independence of the adiabatic index under locally kinetic conditions where $Te≫Ti$⁠.

One consequence of simulating in two dimensions is the plasma is effectively treated as uniform in the z direction. Thus, the simulations will neglect ion flux out of the central plane, causing the absolute density to be higher in the simulations. The lack of self-similarity in plasma expansion can result in differences in the evolution of the density distribution. However, the absolute density does not play a role in driving the plasma expansion or the sound speed. The agreement between experiment and simulation for the renormalized density profile and absolute velocity distribution gives confidence that the electron temperature predicted by the simulation is a reasonable estimate of the actual electron temperature. We use this electron temperature to calculate the ion acoustic wave speed for Fig. 4(b).

Knowledge of the electron temperature is essential to estimate the ion acoustic wave speed. We are unable to directly measure the electron temperature; however, the initial electron temperature is set by the ionization laser detuning above the threshold^7,8,26 (⁠ $Te(0)=160$ K for the data discussed here). We can also reliably assume global thermal equilibrium for the electrons because the MFP is on the order of the size of the plasma, and the collision time is short compared to the plasma expansion time. Initially, for 160 K and a density of $15×1014$ m⁻³, the MFP for an electron moving with the thermal velocity is 2.4 mm.

We compare three different estimates of the electron temperature and ion acoustic wave speed for each time point during the plasma expansion. An upper bound for both is given by $Te=Te(0)$⁠, which neglects adiabatic cooling of the electrons during the plasma expansion. This yields a lower limit on the Mach number. Alternatively, an approximate value of $Te(t)$ can be obtained from a simple estimate of the volume V of the plasma at each time derived from the plasma images assuming cylindrical symmetry and adiabatic cooling of the electrons according to $TeV2/3=constant$⁠.²³ Finally, from the 2D fluid simulation, we extract the simulated $Te(t)$⁠, which we take as our most reliable value. Figure 4(b) shows the velocity transects scaled by the ion acoustic wave speed using $Te(t)$ from the hydrodynamic simulation. For calculating $vIAW$⁠, the ion temperature is negligible.

Figure 6 shows the Mach number as a function of time for the left side of the plasma (⁠ $x<0$⁠). The relative velocity is taken as the difference between the peak velocity and the mean velocity in the region of relatively constant velocity near the edge of the plasma. Initially, M < 1, but for expansion times greater than $∼7 μ$s, the Mach number slightly exceeds unity (or comes very close to unity when using the lower bound). $7.3 μ$s is also the time where the velocity profile develops a sharp kink at the peak [Fig. 4(b)], further suggesting the formation of a shockwave.

The plasma also displays localized heating in the ion temperature where the density and velocity gradients are highest, as shown by Figs. 4(c) and 5. Initially (e.g., for plasma evolution time $t=0.3 μ$ s), the ion temperature is relatively smooth. Disorder-induced heating (DIH)^46,47 should yield a density-dependent temperature, with the highest value in the plasma center of 1.2 K for density of $14×1014$ m⁻³. The observed temperature is higher (2.6 K in the center). This may arise from kinetic effects and initial non-neutrality in the region of very high density gradient in the plasma center.³¹ The center of the plasma adiabatically cools throughout the plasma expansion, as expected.

As wave steepening develops in the plasma, increases in the ion temperature form and grow at the density and velocity jump. This behavior becomes prominent at 7.3  $μ$s in Fig. 4(c). Ion temperatures in the front exceed 10 K where the density and velocity jumps are greatest. Plotting the density, velocity, and temperature transects on the same plot, as shown in Fig. 5 for 7.3  $μ$s, shows the overlap of the temperature features with the front. The numerical simulation only displays a small temperature rise in the region of the front, most clearly seen in Fig. 5(a). This suggests that the temperature spikes arise from kinetic effects, which are not included in the simulation.

It is important to consider that an overestimation in the ion temperature could arise from artificial spectral broadening introduced by the finite resolution of the imaging system and analysis. This is a concern because of the large velocity gradient at the location of the temperature spike, and because we measure the spectra over a finite region size (0.1 mm).

The magnitude of the largest velocity gradient in the data is $3.6×105s−1$ at 12  $μ$s [Fig. 4(b)]. For a 0.1 mm region, this introduces a velocity spread of $Δv=36$ m/s. A simple upper bound on the increase in the temperature extracted from spectral analysis is $ΔT≈mΔv2/2kB=7$ K, which does not account for the observed temperature in excess of 20 K.

To better estimate the spectral broadening due to the velocity gradient, we developed a simulation that divides each 0.1 mm region into smaller subregions in which the velocity spread can be ignored. The simulation creates spectra for each subregion and then adds them together to create the simulated spectrum for the original 0.1 mm region. Simulated spectra are then analyzed in the same way as the experimental data and yield only about 2 K of increase in the measured temperature due to the finite region size for the largest velocity gradients.

As discussed earlier, the electron MFP is on the order of the plasma size, and the collision time is short compared to the timescale for plasma expansion. This allows us to assume global thermal equilibrium for the electrons. The situation for ions is more complex. The MFP for an ion moving at the thermal velocity for 1 K is 28  $μ$ m for a typical density of $1×1014$ m⁻³. This is small relative to the size of the plasma, which is consistent with the lack of global thermal equilibrium observed in the system [Fig. 4(c)]. In regions away from the expanding fronts, the ions quickly establish local thermal equilibrium because the collision time and the timescale for DIH is $∼1 μ$ s.⁴⁷ 

At the front, however, the large relative velocities of ions in the low density and high density sides of about 100 m/s yield much longer MFPs. For the typical density on the high density side of $1×1014$ m⁻³, the MFP is 6 mm, which is very large compared to the size of the shock front. Therefore, we expect a significant fraction of the population of low density ions in front of the shock to penetrate into the high density region overtaking it. A detailed description of this process would require a calculation of the stopping power of the ions, including possible effects from strong coupling,⁵¹ which is beyond the scope of this paper. We expect that the population of ions in the high density region expanding outward is in local thermal equilibrium with itself due to the high collisionality. Any streaming population is too small to detect relative to the high density population. LIF spectra for plasma in the transition region are fit well with a line shape that assumes local thermal equilibrium.

We have presented evidence of wave steepening and signatures of shock formation in UNPs formed with an initially sharply peaked density distribution. As the plasma expands, sharp transitions in density and velocity form at expanding fronts that also display localized ion heating where plasma from the inner, high-density region overtakes outer, lower-density plasma. Using the electron temperature taken from a fluid simulation, Mach numbers slightly greater than unity are observed. This work establishes that conditions for shock formation, or near to it, can be obtained in UNPs with sculpted density distributions.

These data present a rich challenge for testing hydrodynamic and kinetic theories and for searching for effects of strong coupling on viscosity through damping of the shock features. Future work can explore the variation of the observed features with plasma density and initial electron temperature. Expanding the field of view for plasma imaging would enable exploration of whether wave steepening continues or is eventually balanced by dissipation.

This work was supported by the NSF/DOE Partnership in Basic Plasma Science and Engineering under Award No. 2107709, the National Science Foundation Graduate Research Fellowship Program under Grant No. 1842494, and the NSF-CAREER award AGS-1450230. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

The authors have no conflicts to disclose.

M. K. Warrens: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Supervision (supporting); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). N. P. Inman: Data curation (supporting); Investigation (supporting); Software (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). G. M. Gorman: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). B. T. Husick: Investigation (supporting); Software (supporting); Visualization (supporting); Writing – review & editing (supporting). S. J. Bradshaw: Conceptualization (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (lead); Software (equal); Supervision (equal); Writing – review & editing (equal). T. C. Killian: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

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