Research ArticlePHYSICS

Ultrafast optical field–ionized gases—A laboratory platform for studying kinetic plasma instabilities

See allHide authors and affiliations

Science Advances  06 Sep 2019:
Vol. 5, no. 9, eaax4545
DOI: 10.1126/sciadv.aax4545


Kinetic instabilities arising from anisotropic electron velocity distributions are ubiquitous in ionospheric, cosmic, and terrestrial plasmas, yet there are only a handful of experiments that purport to validate their theory. It is known that optical field ionization of atoms using ultrashort laser pulses can generate plasmas with known anisotropic electron velocity distributions. Here, we show that following the ionization but before collisions thermalize the electrons, the plasma undergoes two-stream, filamentation, and Weibel instabilities that isotropize the electron distributions. The polarization-dependent frequency and growth rates of these kinetic instabilities, measured using Thomson scattering of a probe laser, agree well with the kinetic theory and simulations. Thus, we have demonstrated an easily deployable laboratory platform for studying kinetic instabilities in plasmas.


It is well known that more than 99% of the observable matter in our universe is in the plasma state. When the velocity distribution functions of the plasma electrons, ions, or both are nonthermal, plasmas are susceptible to kinetic instabilities (13). Experimental validation of the theory of these instabilities and computer codes is predicated upon having a direct knowledge of the initial velocity distribution functions of the plasma species. This has proven to be illusive in the laboratory. With the advent of intense ultrashort-pulse, near-infrared lasers, it has become possible to ionize atoms and/or molecules of a gas in a few laser cycles (a few femtoseconds—a time duration faster than inverse of the electron plasma frequency) and to generate anisotropic and/or nonthermal electron velocity distribution (EVD) functions by a process known as optical field–induced or tunnel ionization (OFI) (4). This ability to initialize the velocity distribution functions affords an opportunity to quantitatively test the kinetic theory of plasmas (5, 6) on ultrafast time scales before electron-electron (e-e) collisions thermalize the electrons. The mechanisms and time scale by which plasma electrons evolve from an anisotropic to a thermal state are important but unsolved experimental problems in basic plasma science. The kinetic plasma instabilities due to anisotropic plasmas are thought to occur in recombination x-ray lasers based on OFI plasmas (7), gamma ray flashes (8) and bursts (9), electron-positron plasmas (10), collisionless shocks (11), growth and filamentation of magnetic fields in plasmas (12), electron cloud effects in storage rings (13) and proton synchrotrons (14), solar corona and interplanetary medium (15), solar wind trapped in the ionosphere by Earth’s magnetic field (16, 17), neutrino-plasma interactions that may occur in supernova explosions (18), quark-gluon plasmas produced in heavy-ion beam collisions (19, 20), and hot-electron transport in fast ignition-inertial confinement fusion (21, 22).

Given the extremely broad range of situations that give rise to kinetic instabilities, a voluminous theoretical body of work exists on the kinetic theory of plasmas. We first give a brief description of three most frequently studied kinetic instabilities that OFI plasmas have enabled us to quantitatively study in the laboratory. It is known that when plasma electrons composed of two or more co- or counterpropagating streams (beams), they can be unstable to both the streaming and filamentation instabilities (23). The streaming instability is electrostatic (k × E = 0), which generates longitudinal waves (charge density modulations) on the beams via inverse Landau damping at the expense of the directional energy of the streaming electrons (6). If the streams are symmetric, then the two-stream instability is aperiodic, which means that it grows without oscillation; otherwise, the instability has a nonzero oscillation frequency. The filamentation instability grows simultaneously with the streaming instability. It begins with azimuthal magnetic fields surrounding noise currents in the plasma. Attraction (repulsion) between co- (counter-) propagating currents causes the currents to coalesce and thereby to amplify the magnetic fields that surround the filaments (2). In general, the filamentation instability is not purely transverse (k · E ≠ 0) unless the counterpropagating currents are symmetric so that they can pinch at the same rate (23, 24). Therefore, it is possible to probe the electrostatic component of the filamentation instability using Thomson scattering (TS). The anisotropic EVD of the plasma electrons may also cause the Weibel instability to grow via a similar mechanism as the filamentation instability (2, 3, 25), but the growth rates and the spectral evolution of the unstable modes could be different. One distinct feature of the filamentation and Weibel instabilities is their aperiodicity, which means that they have a near-zero frequency.

Although a great deal of theoretical work exists on the kinetic instabilities in plasmas (26), there have been relatively few direct laboratory verifications (2731) of these instabilities because a suitable (and relatively simple) plasma platform that would allow initialization of known anisotropic EVD functions was not available. In the past, these instabilities have been studied by passing relativistic electron beams through plasmas (27, 31) or by creating two interpenetrating plasmas (28, 29). We note that solar wind measurements have been used to test the kinetic theory of the proton firehose and mirror instabilities (32). Here, we show that an ultrafast OFI helium plasma with a known polarization-dependent anisotropic EVD (33) is susceptible to kinetic streaming, filamentation, and Weibel-like filamentation instabilities and measures the growth rates and frequencies of these instabilities using time-resolved TS (34). The measurements are then compared against self-consistent (exact) particle-in-cell (PIC) computer simulations and with theory with good agreements.


EVD function initialization

In the experiment and simulations, we initialize the anisotropic EVD functions by ionizing the first (that produces He1+ ion; IP, 24.6 eV; Ith, ~3 × 1015 W/cm2) and the second (that produces He2+ ion; IP, 54.4 eV; Ith, ~3 × 1016 W/cm2) helium electron using either circularly polarized (CP) or linearly polarized (LP) Ti-sapphire (λ0 = 800 nm) laser pulses. Here, IP is the ionization potential of the electron, and Ith is defined as the required laser intensity to ionize >90% of the He atoms via a tunneling mechanism (35). The details are given in a separate publication (33), but we summarize them here, as they are central to the rest of the paper. The bi-Gaussian pump laser pulse has a full width at half maximum (FWHM) duration of ~50 fs and a spot size of w0 ≈ 8 μm (at which radius the intensity drops to 1/e2 of the peak intensity), and the peak intensity of the laser pulse is ~1 × 1017 W/cm2, sufficient to ionize the first He electrons early during the rising intensity of the laser pulse and the second He electrons approximately 10 fs later. In this process, electrons are either ejected in the (x-y) plane transverse to the direction of propagation of the laser pulse (z) in the CP or along y, the direction of polarization in the LP case. In both cases, the EVD function of the second He electron is “hotter” than that of the first He electron, as shown in Fig. 1 (A and B). These results, taken right after the passage of the laser pulse, are from a three-dimensional (3D) PIC simulation using OSIRIS (36) that uses the Ammosov-Delone-Krainov (ADK) (4) theory to determine the times and the phases of the laser electric field when ionization occurs and then self-consistently determines the energies of the electrons after the passage of the pulse. As explained in (33), the electron momentum distribution resembles a “double donut” shape in the transverse (x, y) plane for the CP case (Fig. 1A) and a two-temperature distribution in the laser polarization direction (y) for the LP case (Fig. 1B). The projected momentum distributions in y and x directions (integrated over px and py, respectively) are shown by the blue lines.

Fig. 1 Initial EVD of OFI helium plasma.

EVDs (A) for CP and (B) for LP laser pulse from 3D OSIRIS simulations. The solid blue lines in (A) and (B) show the projected distributions. In the CP case (A), the projected distribution deviates significantly from a Maxwellian distribution having the same root-mean-sqaure (rms) temperature of 470 eV, as shown by the red dashed line. In the LP case (B), the projected distribution can be well approximated by a two-temperature (1D Maxwellian) distribution with THe1+=60 eV and THe2+=60 eV = 214 eV. The blue lines in (C) and (D) show the measured TS spectrum for CP (C) and LP (D) for an initially fairly low plasma density of 6.6 × 1017 cm−3. The red dashed lines in (C) and (D) are fits to the measured spectrum (see the main text).

We have confirmed that OFI plasmas produced by CP and LP laser pulses do indeed have EVD functions, as depicted in Fig. 1 (A and B) by TS of a probe laser pulse. The blue curves in Fig. 1 (C and D) show the measured TS light spectrum (33). In the CP case (Fig. 1C), because the EVD is grossly nonthermal, we fit the measured TS spectrum by assuming that the frequency of the probe laser light is simply Doppler-shifted by an amount kv by the electrons having a velocity distribution shown in Fig. 1A. This procedure produces an excellent fit (red dashed line) to the measured scattered photon spectrum of Fig. 1C. In the LP case (Fig. 1D), the scattered light can be fitted with a two-temperature but 1D Maxwellian distribution with T1 = 22 eV and T2 = 180 eV (red dashed line). These measured values agree reasonably with the expected values from the PIC simulation, T1 = 45 eV and T2 = 160 eV, where a cos2(π/6) factor has been included to account for the projection due to the scattering angle.

Origin of the kinetic instabilities in OFI plasmas

Figure 2 shows the results of 2D PIC code simulations of the kinetic instabilities triggered by OFI in helium plasma, which is considered here. The plasma is ionized by a CP laser and therefore expected to have an initial distribution function similar to that shown in Fig. 1A. Because such a distribution consists of radially counterpropagating streams, the plasma is unstable to both the streaming and filamentation instabilities, as previously explained. These instabilities grow rapidly following the passage of the laser pulse, and they are self-consistently seeded by the noise fluctuations in PIC simulations. The characteristic electrostatic (Ey) and magnetic (Bx) fields and the density fluctuations associated with these instabilities are shown in Fig. 2 (A, B, and C, respectively). In Fig. 2 (D and E), we zoom in on the two regions marked by the boxes in Fig. 2C to show the details of the density fluctuations at two different places in the simulation and therefore at two different time delays with respect to the pump laser. Panels F and G of Fig. 2 show the 2D fast Fourier transform (FFT) of the density fluctuations (the k-space) shown in panels D and E (Fig. 2), respectively.

Fig. 2 2D simulations show OFI-triggered kinetic streaming and filamentation instabilities in a helium plasma.

The plasma (ne = 5 × 1018 cm−3) is ionized by a CP laser (τ = 50 fs, w0 = 8 μm, I = 1.6 × 1017 W/cm2). The Ey field, Bx field, and density fluctuations associated with the instability are shown in (A), (B), and (C), respectively. (D) and (E) are zoom in of the regions marked by the boxes in (C). The corresponding k-space of these density fluctuations is shown in (F) and (G), where the two dots mark the k of the waves being measured in experiments and where the 400-nm (800 nm) probe is used for CP (LP) pump pulses. (H and I) and (J and K) show the transverse phase space of He1+ and He2+ electrons ionized by CP and LP lasers, respectively. These results are from simulations with higher resolutions. The color bars represent the density of the electrons [in arbitrary units (a.u.)]. The simulation box is 35 μm wide in y. Because the laser only ionizes the central 20 μm of He, a 30-μm window is shown in these plots. In all cases, the electrons inside a Δz = 2-μm slab at z = 20 μm are used to show the phase space. (H) and (I) are taken 0.14 ps while (J) and (K) are taken 1.9 ps after the laser has passed the slab. The gray dashed lines mark the locations of the thin sheaths. The direction of the arrows indicates the shift of the momentum distributions.

One can see from Fig. 2 (A and B) that both the streaming (Ey field) and the filamentation (Bx field) instabilities begin to grow soon after the plasma is created. Initially, the wave vector of the streaming (filamentation) instability is primarily along (perpendicular to) the streaming direction of the electrons (i.e., kStreamky and kFkz) as expected. In the linear regime that lasts for a very short time (≪1 ps), there is no coupling between these instabilities; therefore, they grow independently. In the nonlinear regime, these instabilities start to couple with each other. Accordingly, the density fluctuations initially consist of strips along the z axis (Fig. 2D) but quickly break up at various angles in the y-z plane and eventually become small-scale and randomly distributed speckles (Fig. 2E). There is also a small wakefield induced by the laser pulse, as seen in Fig. 2A, but it is irrelevant to the topic of this work as it simply adds an additional oscillating energy component to the longitudinal or transverse momentum of the electrons. This is so because the ks being probed in this experiment are far off from that of the wake feature, as illustrated in Fig. 2F. In addition, one can see that the plasma is bounded by two sheaths—the outer one separates the plasma from the surrounding neutral gas, and the inner one separates the core of the plasma where helium is double-ionized from the surrounding single-ionized region. The sheath fields have a large electric field (Ey ~ 2 GV/m), and electrons that have energies less than the sheath potential are reflected by the sheaths.

It can be seen from the k-space evolution of the density fluctuations (Fig. 2, F and G) that the wave vector of the instability is primarily parallel to the y direction when the laser has just passed because the instability has the largest growth rate in the y direction. Nevertheless, after only 0.4 ps, the wave vector of the instability has spread along kz, as shown in Fig. 2G. As will be discussed later, we measured waves at oblique angles with different ky and kz at locations marked by the circles in Fig. 2 (F and G). In these measurements, although the wave vector being probed, k, is fixed by choosing a particular scattering angle, it is possible to simultaneously probe the streaming and the filamentation instabilities because they have different frequencies.

The streaming instability saturates and damps very quickly [~1 ps or in ~300 μm in space, as seen in Fig. 2 (A and B)] because of the isotropization of the counterpropagating streams caused by collisionless transverse phase space diffusion, as shall be shown later. We expect the filamentation instability, which is also driven by these streams, to have a similar temporal behavior. However, the magnetic field Bx continues growing past 1 ps, as shown in Fig. 2B, suggesting that, at late times, a Weibel-like filamentation instability driven by a reduced but finite temperature anisotropy of the electrons begins to dominate in the plasma.

Onset of the kinetic instabilities in simulations

Snapshots of the phase space distribution of the He1+ and He2+ electrons for both the CP and LP cases are shown in Fig. 2 (H to K). Note that these results are from 2D simulations with higher resolutions (see Materials and Methods). We can see from Fig. 2 (H and I) that, in the CP case, the streaming motion of electrons is initiated directly by the ionization process and therefore includes streaming not only between the counterpropagating He2+ electrons but also between the co- and counterpropagating He2+ and He1+ electrons. One can see in Fig. 2 (H and I) the modulation of the electron density of the streams. The wave amplitude n1/n0 is just few percent, yet the waves trap background electrons (seen here as closed loops in the py versus y phase space). The He1+ electrons are heated in this process at the expense of He2+ electrons. In the LP case (Fig. 2, J and K), the streaming is mainly between the highest energy He2+ electrons as they bounce off the sheath fields and the remaining quasi-stationary electrons. Because the electrons produced by an LP laser are cooler than those produced by a CP laser and in any case the streams are not established until the He2+ electrons bounce off the sheath, we expect the streaming instability to grow more slowly in the LP case compared to the CP case.

In our case, there exist asymmetric streams (e.g., the copropagating He1+ and He2+ electrons in the CP case and the reflected He2+ electrons streaming against most quasi-stationary electrons in the LP case). Here, asymmetric means the streams have nonzero average drift velocities, different densities, and different transverse temperatures. Kinetic theory shows (see Materials and Methods) that, in this case, the fastest growing streaming mode has a nonzero oscillation frequency (hereafter referred to as the electron feature) proportional to the kvd, where vd is the relative drift velocity (this may be the drift velocity between the two different species at a particular angle that changes with time) that depends on the ionization process and not on the plasma density. As can be seen in Fig. 2 (H and I), this fastest growing streaming mode is initially driven primarily by the copropagating He1+ and He2+ electrons: The instability grows only in locations where copropagating He1+ and He2+ electrons overlap with one another. In the LP case, the two-stream mode is primarily driven by the reflected He2+ electrons and the slow-moving remaining electrons. Therefore, in the CP case, the streaming instability grows simultaneously everywhere in the plasma, whereas in the LP case, the instability starts to grow from the outer edge of the plasma and propagates inward (see the Supplementary Materials).

The nonoscillating filamentation and Weibel instabilities have near-zero frequencies (23) (hereafter referred to as the zero-frequency feature). As mentioned earlier, the filamentation and Weibel mode will also have an electrostatic component (23, 24) that can be probed by the TS diagnostic. Kinetic theory enables us to predict the frequency and the growth rates of the unstable modes being probed as a function of time at different plasma densities so that a comparison can be made between the theory, simulations, and the experiment.

Observation of kinetic instabilities in optical field–ionized plasmas

In the experiment, we choose to probe a wave vector that will give information about the kinetic instabilities but will not be affected by the wakefield effects. Figure 3A shows the sketch of the k-matching diagram for probing the instabilities used in our experiments. We use either a 400-nm laser (probe 1, blue) or an 800-nm laser (probe 2, red) with a 5-nm bandwidth and a nominally 100-fs pulse width to probe the electrostatic components of the instabilities by TS. The scattered light is collected at an angle of either 60° (scattered beam 1) or 150° (scattered beam 2) with respect to the incident probe. The wave vector of the instability being probed, k, is determined by the k-matching condition, ks = ki ± k, as shown by the dashed magenta lines. The locations of these measured ks are marked by the circles in Fig. 2 (F and G).

Fig. 3 TS diagram and examples of measured TS spectra.

(A) k-matching diagram where a helium plasma produced by a 50-fs, 800-nm CP (LP) pump laser is diagnosed by a 400-nm, probe 1 (800-nm, probe 2) laser traversing through the plasma with a variable delay. The measured time-resolved TS spectra are shown in (B) and (C) for the CP and LP pump, respectively. Note that the time scales for the two polarizations are different. The dashed lines mark the position of the expected plasma frequency corresponding to the plasma density. The entire dataset is obtained by scanning the timing in 50- to 200-fs steps, and each step is the average of 20 individual scattering events. Time t = 0 is defined as the time when pump and probe overlap with one another (determined by locating the position of the ionization front seen in a shadowgram formed by the probe at the same location as the probe beam).

The spectra shown in Fig. 3B (C) are measured with a CP (LP) pump and a 400-nm (800-nm) probe. There are two major notable features of these spectra. First, the electron feature grows to a saturated level and then damps in a time duration that is much shorter than the e-e collision time as we shall see. Second, the spectral shift of the electron feature does not follow the familiar ωp shift as would be the case if the probe photons were being scattered from the usual Langmuir wave (Fig. 3, B and C, black dashed lines). In the CP case, the frequency spectrum is very broad extending from or even above ωp all the way to the zero-frequency feature. One may also notice that the blue satellite of the electron feature is somewhat stronger than the red satellite. In other runs, the asymmetry is far less pronounced and shows similar trends as the blue satellite. The peak frequency of the electron feature and the existence of the zero-frequency feature are essential evidence for the streaming and filamentation instabilities, respectively, while the longer-time behavior of the zero-frequency feature is consistent with the Weibel-like filamentation mode.

Streaming instability induced by circular polarization

We will first discuss the case of the CP-induced instabilities (Fig. 3B). The amplitude of the blue-shifted electron satellite is plotted in Fig. 4A as a function of the probe delay in log scale (blue circles). The signal (data in arbitrary units of counts) is raised to the power of 0.5 to obtain a quantity that is proportional to the wave amplitude (density fluctuation) because the measured signal is proportional to the scattered power, which, in turn, is proportional to the square of the density fluctuations (34). The red line shows the wave amplitude obtained from a 2D PIC simulation performed using the experimental parameters (see Materials and Methods). The black dashed line shows the initial (exponential) growth rate predicted by the kinetic theory (see Materials and Methods). We note that, apart from an ~100-fs offset, the temporal variation of the measured signal shows a very good agreement with the kinetic theory and simulations. Some of this discrepancy in the time of onset is due to an uncertainty in knowing the shape of the rising intensity of the laser pulse, the uncertainty in measuring the relative delay between the plasma pump and the probe (~100 fs), and the differences in the noise level between the experiment and the simulations. The instability grows at the expense of the directional energy of the electrons; therefore, the growth rate of the instability decreases with time until it saturates at ~0.5 ps (~70 ωp−1) and then begins to damp. The nonlinear phase of the instability follows, when the change of the distribution function due to the interaction of the waves with electrons is substantial, and lasts for another ~0.5 ps. During this phase, some of the electrons’ directional energy transferred to the electrical field energy of the unstable waves is returned back to the electrons via electron trapping. As a result of these wave-particle interactions, the transverse phase space diffuses and the process of isotropization begins.

Fig. 4 Streaming, filamentation, and Weibel instabilities induced by circular polarization.

(A) The blue circles show the measured magnitude of the electron feature as a function of time (derived from Fig. 3B) for a plasma density of 6 × 1018 cm−3. The red line shows temporal evolution of the density fluctuation magnitude (instability wave amplitude) for the same probed value of k obtained from a 2D PIC simulation. The black dashed line shows the exponential growth of a linear wave from kinetic theory. The horizontal error bars mark the uncertainty in determining the t0 when ionization is completed. The vertical error bars represent the ±2σ confidence level. (B) Measured (magenta circles), predicted (blue line), and simulated (green squares) growth rates of the instability. (C) Measured (magenta circles) and predicted (blue line) frequencies of the streaming instability for the same range of densities as in (B). (D) Measured magnitude of the zero-frequency feature (green) and the electron feature (blue) within the first 2 ps. (E) The measured (circles) and calculated (green line) initial growth rates of the filamentation instability. The horizontal error bars show the uncertainty of density measurement, and the vertical error bars represent the ±σ confidence interval of the deduced growth rate. The blue dashed line shows the growth rate of the zero-frequency mode of the streaming instability. (F) Measured magnitude of the zero-frequency feature as a function of time (green squares, derived from Fig. 3B). The solid (dashed) purple line shows the evolution of the amplitude of the electron (ion) density fluctuation at the same k that is being probed in the experiment. The red dotted-dashed line shows the maximum growth rate of the Weibel instability calculated using the simulated EVD at t = 1 ps.

We repeated the measurements at four different plasma densities, and from each of these measurements, we calculated the initial growth rate and the frequency of the fastest growing mode at the measured k of the streaming instability. The results are shown as magenta circles in Fig. 4 (B and C, respectively). The blue lines show predictions of the kinetic theory (see Materials and Methods). The green squares in Fig. 4B show the 2D simulation results of the growth rate of the instability having the same k (i.e., ky = −k0 and kz = 1.7 k0 for the blue satellite) as that measured in the experiment. Here, k0 = 2π/λ0 is the wave vector of the 800-nm laser. The frequency of the instability shown in Fig. 4C was calculated using the spectral shift of the peak of the blue satellite (see Fig. 3B) before the instability saturates. The measured initial growth rates show reasonable agreement with the theory. One of the possible reasons the measured growth rate is smaller than the analytical and computational predictions might be that the experiment is in 3D, whereas the calculations and simulations are in 2D. Collisions may also reduce the growth rate of the instability (37) (also see the Supplementary Materials).

In contrast, the measured and predicted frequencies, using the kinetic theory, of the streaming instability show an excellent agreement. This is because the frequency of the streaming instability is given by ω ≈ kvd = kyvd, where vd ≈ 0.04c is the drift velocity between the copropagating He1+ and He2+ electrons, which is predominantly determined by the ionization process itself. Because of the large spread in vd, a wide spread of unstable streaming modes is seen to grow (at a given k) in the TS electron spectrum (Fig. 3B). The frequency of the streaming modes driven by the counterstreaming He1+ (He2+) species is zero because the two beams are symmetric; therefore, they do not contribute to the measured electron feature.

The phase velocity of the streaming instability is vϕ ≡ ω/kvdcosθ, which is independent of the plasma density. Here, θ is the angle between the wave vector k and the drift velocity vd. This is confirmed in the experiment as shown in Fig. 4C, where the measured oscillation frequency of the instability follows the analytical prediction, which is almost linear with a slope of vϕ. As described earlier, initially, the wave vector of the streaming instability is primarily along the streaming direction so that cosθ ~ 1 and therefore vϕvd; thus, strong resonance occurs between the copropagating He1+ and He2+ electrons, as evidenced by the energy gained by He1+ electrons and energy lost by He2+ electrons seen in Fig. 2 (H and I). This resonant interaction hastens the collisionless phase space diffusion from initially double donut–shaped distribution to a single quasi-Maxwellian distribution in approximately 1 ps, leaving a bitemperature distribution (hot in the transverse plane but cold in the longitudinal direction), which is unstable to the Weibel-like filamentation instability.

Filamentation and Weibel-like filamentation instabilities with circular polarization

The radially counterpropagating streams initiated by the tunnel ionization using a CP laser readily also drive the filamentation instability. Because of the anisotropy of the EVD, Weibel-like filamentation instability may also grow. In Fig. 4D, we plot the measured growth of the zero-frequency feature (green squares) and the electron feature (blue circles) together to illustrate the very similar behavior of these two signals. For the k being probed, the zero-frequency filamentation mode is actually first to appear above the measurement threshold, followed by the streaming mode (see Figs. 3B and 4D), but both have similar growth rates. Following a rapid growth, both decay within ~1 ps, a time duration shorter than the ion plasma period (2πωpi−1 ≈ 3 ps), which suggests that the zero-frequency mode corresponds to the filamentation instability instead of the usual ion acoustic waves.

Figure 4E shows the initial growth rate of the filamentation mode, where the magenta circles show the measurements and the green line shows the solution of the dispersion relation (see Materials and Methods). Once again, a very reasonable agreement between the measurements and the kinetic theory is seen. We note that kinetic theory predicts that there may also be a nonoscillating branch of the streaming instability. The growth rates of the two are very comparable, and therefore, this nonoscillating branch of the streaming instability also likely contributes to the measured zero-frequency feature. However, it is the recurrence of the measured zero-frequency feature (see Fig. 4F) that is the strongest indicator of the filamentation and Weibel-like filamentation instability because recurrence of the streaming mode is not expected once the streams no longer exist, as evidenced by the decay of the electron feature.

In Fig. 4F, the green line shows the temporal evolution of the zero-frequency feature in our experiment (400-nm peak in the spectra shown in Fig. 3B), which shows two distinct peaks. As discussed earlier, the first peak (0 and 1 ps) is likely due to the filamentation instability, which rapidly reduces the anisotropy of the EVD. Next, we will show using simulation that the second peak (1 to 7 ps) is consistent with a Weibel instability driven by the residual anisotropy in the distribution function. In our simulations, we cannot model both ionization of He and ion motion at the same time. The latter is important for accurately and directly determining the electrostatic component of the Weibel instability. Therefore, we performed a simulation where the plasma was preionized and had an EVD similar to that of the OFI plasma (see Fig. 5A). In this simulation, the motion of He ions was self-consistently included. The amplitudes of the density fluctuations of electrons (δne) and ions (δni) at the same k as that being probed in the experiment as a function of time are shown in Fig. 4F by the solid and dashed purple lines, respectively. The evolution of electron density fluctuation δne(kexpt) shows a peak at t ≈ 0.5 ps, which fairly well tracks the first peak of the measured zero-frequency feature. Because there has not been enough time for the heavy ions to follow the motion of electrons, the filamentation instability contributes to the electron density fluctuations through the different pinching rates of the electron streams. However, as time progresses, the scale length of these density fluctuations evolves to larger k, allowing the Weibel instability to grow to detectable levels through the ion density fluctuations. In the experiment, this manifests as the recurrence of the zero-frequency feature, while in the simulations, this is seen as the growth of the ion density fluctuations after the first ~2 ps, as shown by the dashed purple line. After that, δni(kexpt) reaches a saturated level with small-amplitude oscillations. The oscillation period is 3.7 ps, which is slightly longer than the ion plasma period (3 ps) but agrees with the ion acoustic wave frequency that corresponds to Te = 150 eV [ω = kcs, cs = (γZκTe/mi)1/2 is the ion sound velocity, γ is the adiabatic index, k = 2k0, and k0 is the wave vector that is being probed]. This late time oscillation of the zero-frequency feature was not unambiguously seen in the experiment in the CP case, although an oscillatory behavior was observed in the LP case seen in Figs. 3C and 6B. In Fig. 4F, we also show the theoretical growth rate of Weibel instability by the red dotted-dashed line. This line represents the maximum Weibel growth rate (γ ≈ 0.02ωp) calculated using the simulated EVD at t = 1 ps, which shows a good agreement with both the simulated and the measured ion density fluctuations. We should note that kinetic theory predicts that the unstable range of Weibel instability covers 0 < k < A1/2ωp/c (38), which is much smaller than the k that is being probed in the experiment. However, it has been shown that the unstable k of Weibel instability can extend to large k through turbulence cascade (39, 40).

Fig. 5 Evolution of the temperature anisotropy of the OFI plasma.

The upper (lower) row in (A) shows the py (pz) distribution function of electrons at t = 0, 1, and 6 ps. The dashed gray line is a Gaussian fit to the distribution. The initial distribution can be approximated by four drifting Maxwellian beams in the transverse plane as indicated by the red line and the arrows. The red dashed line is a Gaussian fit to the pz distribution. (B) The blue line shows the anisotropy from the same simulation as in (A), which does not include collisions. The red line shows the simulation of anisotropy evolution of a preionized plasma with only Coulomb collisions included. (C) The average magnetic field energy as a function of time shows two distinct growth phases corresponding to filamentation and Weibel regimes, respectively.

Fig. 6 Instabilities in a plasma ionized by an LP laser.

(A) Measured (blue) and simulated (red) evolutions of the magnitude of the electron density fluctuations of the streaming instability. (B) The measured magnitude of the zero-frequency mode as a function of time, displaying an oscillatory behavior with a roughly ion acoustic period.

Isotropization of the plasma electrons

We have tracked the evolution of the EVDs and the temperature anisotropy A ≡ (T/T − 1) of the OFI plasma in a 2D simulation, and the results are shown in Fig. 5. In this simulation, the ionization and evolution of the plasma were self-consistently modeled, but Coulomb collisions were not included to isolate the effect of the instabilities on the temperature anisotropy. The upper (lower) row of Fig. 5A shows the py (pz) distribution of electrons within a thin slab (Δz = 2 μm and Δy = 2 μm) at z = 17 μm in the 2D simulation at t = 0, 1, and 6 ps. As can be seen in Fig. 5A, the initial py distribution is highly nonthermal (deviates significantly from the Maxwellian distribution shown by the gray dashed line). The initial distribution can be approximated by four drifting Maxwellian beams (the red line; the streams are indicated by the black arrows) in the transverse (perpendicular) plane. The initial distribution has an extremely low temperature (~1 eV) in the longitudinal (parallel) direction (the red dashed line). However, the multiple beam structure disappears because of the collisionless phase space diffusion (Landau damping and particle trapping by the waves excited on/by the streams and the v × B motion of the electrons due to filamentation and Weibel modes). The distribution approaches a quasi-Maxwellian state within a few picoseconds. Consequently, the root-mean-square (rms) temperature T drops, while the rms temperature T increases significantly.

As a result of these kinetic instabilities, the anisotropy of the plasma drops rapidly from an initial value of a few hundreds to ≈10 within 1 ps, as shown by the blue line in Fig. 5B. Thereafter, the Weibel-like filamentation instability continues to isotropize the plasma but with a smaller rate to further reduce the anisotropy to <1 in about 7 ps. The magnetic fields (Fig. 5C) also show two distinct growth phases that correspond to the filamentation and the Weibel instability, respectively. Simulation shows that, as the Weibel instability saturates, the magnetic fields self-organize to a quasi-static helical structure as predicted in (38).

To ensure that such a rapid drop of the anisotropy is due to instabilities, we performed a simulation using a preionized plasma with similar initial EVD (represented by the red lines in Fig. 5A) and included only the Coulomb collisions [both e-e and electron-ion (e-i)]. Without the presence of the kinetic instabilities, electrons efficiently exchange their energy (momentum) only through the e-e (e-i) collisions so that they isotropize and thermalize after tens of picoseconds (see the Supplementary Materials). The result is shown by the red line in Fig. 5B. It can be seen that, except for the first few hundred femtoseconds, the anisotropy in this collision-only simulation drops much slower than that in the simulation that only includes instabilities. Thus, we confirm that the collisions do not play a significant role in the first ten picoseconds after the formation of the plasma over the range of plasma densities used and that, during this time, the isotropization of the plasma is dominated by the kinetic instabilities. Nevertheless, collisions will eventually thermalize the plasma (see the Supplementary Materials).

Kinetic instabilities induced by linear polarization

In contrast to the CP case where electron streams are directly initiated by the ionization process, in the LP case, the initial EVD consists of two 1D Maxwellian distributions—not a clear “bump on tail” distribution that is unstable to kinetic instability (6). The electron streams are established only after the hotter He2+ electrons have bounced off the sheath electric field at the plasma-vacuum boundary (see Fig. 2K). In this case, the instability is driven by the reflected electrons, as they propagate through most slower moving electrons, and therefore has a nonzero density-independent frequency of ω ≈ kvd = kyvd, where vd ≈ 0.04c is the drift velocity of the reflected electrons, which was confirmed in the experiment. The frequency spectrum of this mode is narrower than that in the CP case because there are only two streams in the LP case.

The growth rate of the streaming instability in the linear polarization case is smaller because only a small fraction of the electrons is reflected to drive the streaming instability. Therefore, it takes a longer time for the streaming instability to grow and saturate. This prediction was also confirmed in the experiment, as shown in Fig. 6A, where the blue circles show the measured magnitude of the streaming instability (blue satellite in Fig. 3C) whereas the red line shows the corresponding simulation result. A remarkable agreement is found between the measurements and the simulation. When making this comparison, the fully kinetic PIC simulation can also be seen as giving the prediction of kinetic theory. Figure 6B shows the magnitude of the measured zero-frequency feature. Once again, the filamentation mode at the k being probed (ky = −1.86 k0 and kz = 0.5 k0) appears earlier than the streaming mode. These two modes grow and damp with a fixed phase relationship with one another; however, both the streaming mode and the zero-frequency mode show similar recurrent peaks as were seen in the mobile ion PIC simulations done with CP. In other datasets, we observed the recurrence of the streaming mode and the same phase relationship between the streaming and the zero-frequency modes (see the Supplementary Materials).


In addition to the examples that we have shown, it is possible to generate other “designer” EVDs by using combinations of different laser polarizations, wavelengths, intensity profile, and ionizing medium. For instance, by using an elliptically polarized laser to ionize the He atoms, one can initialize counterpropagating streams along one direction (33). The drift velocity and the transverse temperature of these streams can be controlled by changing the polarization ellipticity and therefore can be used to suppress one of the (streaming or filamentation) instabilities. Because the EVDs are predominantly determined by the OFI process, one may expect similar instabilities even in solid density plasmas. The growth of the instabilities will be faster (γ ∝ ωp), the wavelength of the predominant mode will be shorter (λcωp1), and the lifetime may be shorter because of an extremely high collision rate. Therefore, one may require the use of femtosecond soft x-ray laser pulses to probe these instabilities in solid-state plasmas. Simulations show that the most unstable k of the instability changes as the plasma evolves, which can be confirmed by simultaneously collecting the scattered light from multiple angles. One can simultaneously measure the time evolution of the electron and ion features at several ks, which will enable the hierarchy of kinetic instabilities to be established (41, 42) and track the evolution of the plasma from a quiescent to turbulent state (15).

In summary, we have shown that ultrafast OFI plasmas are nonthermal and have a large velocity anisotropy and that these plasmas undergo the streaming and filamentation instabilities, followed by a Weibel-like filamentation instability that act to isotropize the plasma. The polarization-dependent frequency and growth rates of these kinetic instabilities, measured using TS of a probe laser, agree well with the kinetic theory and simulations. Therefore, we have demonstrated an easily deployable laboratory platform for studying kinetic instabilities in plasmas.


Experimental setup

The laser used in the experiment had a pulse duration of 50 ± 5 fs and a central wavelength of 800 nm. It was split into two pulses: the pump and the probe. The polarization of the pump laser was changed to circular by sending the LP laser pulse through a quarter-wave plate before being focused down to a spot size w0 of 8 μm by an F/6 off-axis parabolic mirror. The laser ionized neutral helium gas ejected from a supersonic nozzle with a diameter of 1 mm. The density of the plasma can thus be adjusted by changing the backing pressure of the jet. The probe beam contained ~50% of the total energy (20 ± 4 mJ) and was p-polarized with respect to the scattering plane. The probe was sent through a 1.5-mm-thick potassium dihydrogen phosphate (KDP) crystal to double its frequency. The pulse duration of the probe beam was measured to be 50 ± 5 fs. The probe beam was colinear with the pump but delayed by ~0.3 ps (estimated) by the group delay induced by the KDP crystal and the wave plate. The scattered light was collected using an F/6 achromatic lens at 60° with respect to the probe propagation direction, and its spectrum was recorded using a spectrometer (Horiba HR-320 equipped with a PI-MAX 4 camera). The slit width at the entrance of the spectrometer was set to 200 μm. The probe beam was changed to be perpendicular to the pump for measuring the instabilities. The delay of the probe was controlled using a trombone arm with a motorized translation stage in the probe beamline. A 0.8-μm probe (duration of 200 ± 10 fs, calculated) was used when the polarization of the pump beam was switched to linear, and the scattering angle was correspondingly changed to 150° to fulfill the k-matching condition to probe the most unstable mode of the instability.


The PIC simulations were performed using the OSIRIS code (36) on the Edison cluster at National Energy Research Scientific Computing Center (NERSC). To obtain the initial velocity distribution shown in Fig. 1 (A and B), 3D simulations (simulation 1 and simulation 2) were performed. The dimensions of the simulation box were 500 × 280 × 280 c0, where ω0 is the frequency of the ionization laser. The box was divided into 2500 × 560 × 560 cells. Eight particles with a quadratic shape were initialized in each cell. A 50-fs (FWHM), 0.8-μm laser with a sin2 intensity profile was injected from the left wall of the simulation box and propagated to the +z direction. The laser pulse had a Gaussian radial intensity profile with a spot size w0 of 8 μm. The peak intensity of the laser was 1.6 × 1017 W/cm2, which gives a normalized vector potential a0 = 0.19 (CP; a0 = 0.27 for the LP case). Transversely uniform helium gas was initialized inside the box with a longitudinal profile that increases linearly from zero at z = 10 c0 to 2.5 × 1018 cm−3 at z = 20 c0 and then remains flat before decreasing linearly from the peak density at z = 490 c0 to zero at z = 500 c0. The ionization of the helium gas in OSIRIS was modeled using the ADK theory (4), and for the above simulations, the ions were immobile. The initial distribution functions of the plasma were calculated using electrons within a Δz = 2-μm-wide slab (we confirmed that the size of the slab did not affect the result) around the focal position of the laser at the time when the laser had just passed. Limited by available computational resources, 2D simulations with higher resolutions were performed to model the evolution of the instabilities. In simulation 3, which generates Fig. 2 (A to G), the simulation box had dimensions of 4800 × 480 c0, which was divided into 24,000 × 1440 cells. Sixty-four particles with a quadratic shape were initialized in each cell. All other parameters were the same as those in simulation 1 except that the uniform region was larger to match the enlarged box. Two simulations (simulation 4 and simulation 5) with higher resolutions were performed to track the evolution of the instabilities for CP and LP cases, respectively. The simulation box had dimensions of 500 × 280 c0 and was divided into 5000 × 2800 cells. Sixty-four particles with a cubic shape were initialized in each cell. All other parameters were the same as those in the 3D simulation. The red line in Fig. 4A was obtained as follows. First, we did a 2D FFT of the density perturbation δne inside a box (90 ≤ z ≤ 410 c0, r ≤ 68 c0) to get the amplitude spectral density; then, we took the average magnitude of the signal at (ky = −1 ± 0.05 k0 and kz = 1.73 ± 0.05 k0), where k0 is the wave number of the 800-nm laser. The green squares in Fig. 4B were obtained by changing the plasma density from 5 × 1017 cm−3 to 8 × 1018 cm−3. In simulation 5, for the LP case, the instability signal was taken at (ky = −1.86 ± 0.05 k0 and kz = 0.5 ± 0.05 k0). The result corresponds to the red line in Fig. 6A. A 2D simulation (simulation 6) with motion of He ions self-consistently included was performed to show the evolution of both the electron and the ion density fluctuations (the solid and dashed purple lines in Fig. 4F, respectively). The simulation box had dimensions of 25 × 12 cp and was divided into 1000 × 480 cells. Periodic boundary conditions were used in both directions. Two hundred fifty-six particles with a quadratic shape were initialized in each cell. Preionized plasma with EVDs shown by the red lines in Fig. 5A was used. The neutralizing He2+ ions were mobile. Collisions were not included in these simulations to isolate the effect of the kinetic instabilities on isotropization of the electrons and also to save simulation time. The results in Fig. 5 are from 2D simulations (simulation 7 and simulation 8), where one included ionization but no collisions (simulation 7, blue lines) and the other used a preionized plasma with a similar initial distribution and included collisions (simulation 8, red lines). For the 2D simulation that included ionization, all the parameters were the same as those in simulation 4 except that the number of cells was reduced to 2500 × 1400 to save time, and 256 particles were used in each cell to suppress artificial numerical collisions to ensure the accuracy within the longer (7 ps) time window. In simulation 8, we used a 1 × 1 cp box with periodic boundary conditions in both directions, and the box was divided to 48 × 48 cells. Preionized plasma with distributions shown by the red lines in Fig. 5A was used. Two hundred fifty-six particles were initialized in each cell. The He2+ ions were mobile. The particle pusher was turned off so that the particles can only change their velocities through e-e and e-i collisions.

Kinetic theory of the instabilities

The frequency and the growth rate of the streaming instability were calculated by solving the dispersion relation 1Σjωpj22(ky2vthjy2+kz2vthjz2)Z(ωkyvjykzvjz2ky2vthjy2+kz2vthjz2)=0 using an open source code (43), where ωpj, vj, and vthj are the plasma frequency, the drift velocity, and the thermal velocity of the jth beam, respectively. The subscripts y and z represent the velocity component in the corresponding direction. For beams drifting along the y direction and that are cold in the z direction, the dispersion relation reduces to the standard form given by 1Σjωpj22ky2vthjy2Z(ωkyvjy2kyvthjy)=0. The growth rate of the filamentation mode was calculated by solving the dispersion relation 1ωp2+c2kz2ω2Σjωpj2ω2vjy2+vthjy22vthjz2Z(ω2kzvthjz)=0. In the above dispersion relations, Z′(ξ) = − 2[1 + ξZ(ξ)] is the derivative of Z(ξ)=π1et2tξdt, which is the plasma dispersion function introduced by Fried. The parameters of the beams were obtained by fitting the initial velocity distribution in the 2D simulation with four drifting Maxwellian distributions (see the top left subplot in Fig. 5A). The parameters of the streams are as follows: n1 = n2 = 0.385n0 and n3 = n4 = 0.115n0, where n0 is the total plasma density; v1y = −v2y = 0.0204c, v3y = −v4y = 0.0616c, and vjz = 0; and vth1y = vth2y = 0.0138c and vth3y = vth4y = 0.0093c. The thermal velocities of these streams in z direction are vth1z = vth2z = 4.2 × 10−4c and vth3z = vth4z = 9.7 × 10−4c. The dotted-dashed line in Fig. 4F corresponds to the maximum growth rate at the most unstable k ≈ 2ωp/c and is calculated using equation 1 in (38) with the anisotropy A = 10.6 and T = 48 eV.


Supplementary material for this article is available at

Supplementary Text

Fig. S1. A time sequence of snapshots of electron density fluctuations.

Fig. S2. Comparison of PIC simulations without and with e-e collisions.

Fig. S3. Relaxation of nonthermal plasma due to Coulomb collisions.

Fig. S4. Recurrence of the streaming mode and the fixed phase relationship between the streaming and the filamentation modes.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported by the Office of Naval Research (ONR) MURI (4-442521-JC-22891), AFOSR grant FA9550-16-1-0139, U.S. Department of Energy grant DE-SC0010064, and NSF grant 1734315. The simulations were performed on the Edison cluster at the NERSC. Author contributions: C.-K.H., C.Z., and K.A.M. performed the experiments. C.-K.H., C.Z., K.A.M., C.E.C., and C.J. analyzed the data. C.Z. developed the theory and performed the numerical simulations with the help of W.B.M. All authors discussed the results and contributed to the preparation of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data and the codes for analyzing the data may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article