Measurement of the sound speed in dense fluid deuterium along the cryogenic liquid Hugoniot

As the lightest and most abundant element in the universe and the dominant constituent in stars and gas giants, the behavior of hydrogen at high densities and temperatures is of fundamental importance for condensed matter physics,¹ planetary sciences,^2,3 helioseismology,⁴ Jovian Seismology,⁵ and inertial confinement fusion (ICF) research.^6,7 Because the equation of state (EOS) and transport properties of hydrogen rich fluid are needed to connect astronomical observations and planetary interiors models, improved knowledge on the equation of state (EOS) of hydrogen is of critical importance to shed new light on the formation, structure, and evolution of gas giant planets and exoplanets.^2,8 This is much needed as the on-going Juno mission is currently measuring the gravity and magnetic fields of Jupiter with an unprecedented accuracy and has already returned puzzling data.³ 

A detailed understanding of the EOS of dense hydrogen isotopes is also critical for inertial confinement fusion (ICF) research which aims at studying thermonuclear fusion of heavy hydrogen isotopes: deuterium (D) and tritium (T) in a laboratory setting. The ICF research effort at the National Ignition Facility (NIF) uses a laser-driven, quasi-isentropic spherical compression^6,9 of an ∼1 mm diameter capsule consisting of an ∼100 μm thick carbon, polymer, or metallic outer shell ablator and an inner ∼70 μm thick DT ice fuel layer. During the initial phase of compression, a series of shock waves of increasing strength (1–100 Mbar) are carefully timed such that shock coalescence does not occur within the fuel. This multi-shock approach aims to approximate isentropic compression of the DT. The compressed state produced by the first of those shocks plays a key role in determining the speed at which subsequent shock waves travel and more importantly the total entropy of the final compression state. Predictive capability of shock timing campaign can be improved by providing benchmarking experimental data about the first shock state.

Furthermore, achieving high-fusion yields on the NIF requires an accurate understanding of hydrodynamic instabilities to achieve as close as possible to an ideal spherical implosion. This challenging goal is being tackled with different philosophies, leading to various experimental platforms that differ in particular by the strength of the first shock that is launched into the DT ice layer. We show in Fig. 1 a summary of experimental data on the pressure-density compressibility along the cryogenic Hugoniot for D₂, along with several EOS model predictions.^10–16 The strength of the first shock state that would be launched in a D₂ fuel surrogate is shown for several platforms and varies from 50 to 450 GPa, illustrating the need for an accurate knowledge of the EOS over a broad range of pressures.

By leveraging a recently developed analysis method^21–23 for high-resolution Doppler interferometric velocity data, we present Eulerian sound speed measurements in compressed deuterium to shock pressures between 50 and 200 GPa along the locus of shock states starting from cryogenic liquid deuterium. We recall that the Eulerian sound speed is the pressure derivative with respect to density (ρ) at constant entropy (s): $cs=(dP/dρ)|s$⁠. In addition to the sound speed, we also determined the Grüneisen parameter defined as the pressure derivative with respect to the internal energy E at constant volume V, multiplied by the volume: Γ = V (dP/dE)$|V$⁠. Because the pressure P is itself the first derivative of the free-energy F with respect to density at constant temperature: P = ρ²(dF/dρ)$|T$⁠, both the sound speed and the Grüneisen parameters are second derivatives of the free energy. They are consequently crucial variables to benchmark EOS models which define the variation of the free energy as a function of density and temperature or equivalently the variation of pressure and internal energy as functions of density and temperature.

These sound speed results, extracted from previously published laser-driven shockwave experiments,²⁰ are found to agree with several of the most accurate equation of state models for deuterium at those conditions up to ∼150 GPa. However, the data indicate that these models fail to reproduce the experimentally observed sound speed at higher pressures approaching 200 GPa. In particular, we unveil a discrepancy between the experimental results and the equation of state model that is most commonly used in inertial confinement fusion at the National Ignition Facility (LEOS 1014).^6,7

In Sec. II, we summarize the experimental results from Hicks et al.²⁰ and show how those results can be utilized to extract additional EOS information. For clarity, we briefly repeat the description of the experiments performed on the Omega laser at the University of Rochester to measure the pressure-density compressibility of deuterium along the cryogenic Hugoniot. In Sec. III, we present out sound speed and Grüneisen parameter measurements and compare those results with recent EOS models, followed by concluding remarks in Sec. IV.

In this work, we revisit the experimental measurements of Hicks et al.²⁰ in order to extract additional thermodynamic information. The previous work of Hicks et al.²⁰ was performed on the Omega laser at the University of Rochester, a frequency-tripled 0.35 μm light neodymium-doped phosphate glass laser system. The focus of that work was to measure the single shock compression of fluid deuterium from 45 to 220 GPa using laser-driven shock compression. The targets used in those experiments (see the inset of Fig. 2) consisted of planar assemblies with a 20 μm plastic polymer (CH) ablator, a 50 or 90 μm Al baseplate, and a cryogenic liquid D₂ sample side-by-side with an α–quartz witness plate.

Flat-top, 3.7 ns long pulse shapes were used and the laser energy and the number of beams were varied to probe different D₂ pressure states. The laser focal region was smoothed using distributed phase plates to produce an 800 μm uniformly irradiated spot ensuring that the drive pressure was uniform over the two section targets during the experiment. In our analysis, we require planar regions of ∼150 μm in each target region (equivalent to ∼5 fringes). On average, the edge rarefactions arrive at the edge of the ∼150 μm planar region ∼4 ns after the shock enters the sample. At this time, the drive is no longer planar/equivalent for both sections of the target and the velocity interferometer system for any reflector (VISAR) data taken after 4 ns are not used in this analysis.

The principal diagnostic was a line-imaging velocity interferometer system for any reflector (VISAR)²⁴ using Mach-Zehnder interferometers and an image relay of the target onto the entrance slit of streak-camera detectors. In the shock pressure regime explored in these experiments, both the deuterium and the quartz are transformed into optically reflecting states with sub-micrometer optical depth at the wavelength (λ = 532 nm) of the probe laser, so that the VISAR tracked the shock velocity histories both in the D₂ sample and in the quartz witness. Due to the short optical depth of these shock compressed states, only the leading shock front is tracked at each time during the transit through the package. Two interferometers with different velocity sensitivities were used to resolve the 2π phase-shift ambiguities. The VISAR fringe positions were determined to be within <3% of a fringe which resulted in a shock velocity precision of ∼1%. The temporal resolution of the streak cameras was ∼40 ps. All inferred quantities presented here were based on a complete reanalysis of the original experiments starting from the raw experimental VISAR images.²⁰ 

In the original report,²⁰ the shock velocity history in the quartz plate was used to infer the shock velocity in the aluminum pusher right prior to the transmission of the shock from the pusher to the deuterium sample. The inferred aluminum shock velocity and the measured deuterium shock velocity were then used in an impedance matching analysis²⁵ to obtain a single point on the pressure-density shock compressibility Hugoniot curve. One difficulty of this analysis is to account for the slight variations of the shock velocity—and pressure—as it travels through the multi-layer target packages, arising primarily from small variations of the laser intensity during the laser pulse. Interestingly, these perturbations can be used to determine the sample sound speed.^21,26

One can extract useful additional information from these records by analyzing these unintended shock velocity perturbations, and in particular, unfolding the differences between the time it takes them to catch up with the leading shock front traveling in deuterium, compared to the time it takes in the reference quartz witness material. This analysis assumes that the pressure drive is planar and uniform across both sections of the target, so that the correlation of the timing and amplitude of acoustic perturbations reaching the leading shock front in the sample and in the witness can be used to determine the sample sound speed, relative to the sound speed in the quartz witness.²¹ 

As described elsewhere^21,26 and demonstrated on B₄C²² and SiO_2,²³ this first-order perturbation analysis requires prior knowledge of the compressibility and sound speed—at the relevant shock conditions—in both the witness plate and baseplate materials as well as the shock compressibility of the sample. We relied on previous studies on aluminum²⁵ and quartz²⁷ to model the EOS properties of those materials. For the compressibility of the deuterium cryogenic Hugoniot, we fit all available $Us−Up$ high-pressure (P > 50 GPa) Hugoniot data^18,20,28,29

Note that this techniques requires a nearly steady shock wave (<10% deviation about the mean shock velocity), a spatially uniform drive across both sections of the target and drive perturbations at the ablator-aluminum interface are identical for both sections of the target (the experiment is terminated when the aluminum/quartz or aluminum/deuterium rarefaction waves arrive at the CH/aluminum interface).

In more detail, one has to model how the acoustic perturbations that are generated at the ablation front refract and change amplitude as they propagate through the shock compressed material from the drive surface to the leading shock front. As derived elsewhere,²¹ the correlation between pressure perturbations into the aluminum baseplate and the D₂ shock front is defined as

where $ΔtD2Δt0$ is the relative change in arrival of perturbations at the D₂ shock front to that of the source, and $MAlD2,d$ and $MAlD2,u$ are downstream and upstream Mach numbers. The terms “upstream” and “downstream” are used in the context of wave fronts to denote the material regions ahead of and behind the wave front, respectively, when viewed from a reference frame fixed to the front. The Mach numbers (⁠$MAlD2,d$ and $MAlD2,u$⁠) are defined relative to the aluminum centered rarefaction fan generated from impedance matching with the deuterium. These parameters are determined from the known aluminum EOS.²⁵ $MD2,d$ is the D₂ downstream shock front Mach number defined as

where $PD2,cs,D2,up,D2,ρD2$ are the pressure, Eulerian sound speed, particle velocity, and density along the D₂ Hugoniot, respectively.

Similarly, the correlation between pressure perturbations in the aluminum baseplate and in the quartz shock front is defined as

where $ΔtQz2Δt0$ is the relative change in arrival of perturbations at the quartz shock front relative to the source and $MAlQz,d$ and $MAlQz,u$ are downstream and upstream Mach numbers defined relative to the aluminum centered rarefaction fan (or reshock depending upon impedance matching) determined from the impedance matching with quartz.²⁷ 

Provided that the EOS of aluminum, quartz, and the D₂ pressure-density shock compressibility is known, Eqs. (2), (3), and (4) can be combined into a single equation with two unknowns: $ΔtQzΔtD2$ and $cs,D2$⁠.³⁰ The parameter $ΔtQzΔtD2$ represents the relative time relationship between the quartz witness and deuterium sample shock velocities which is measured in these experiments (as shown in Fig. 2). By experimentally measuring this parameter, we are then able to determine the sample sound speed.

To determine the Doppler shift between the quartz and deuterium shock front velocities $(ΔtQzΔtD2)$⁠, we use a Levenberg-Marquardt nonlinear least squares optimization routine to linearly transform (in time and amplitude) the perturbations about the average deuterium shock velocity (⁠$u¯D2$⁠) onto the perturbations about the average quartz shock velocity (⁠$u¯Qz$⁠) as shown in Fig. 2. We aim to determine the temporal $(ΔtQzΔtD2)$ and amplitude scaling factors $(ΔuQzΔuD2)$⁠, such that the measured quartz shock front velocity perturbations [$δuQz(t)$] are mapped onto the measured deuterium shock velocity perturbations [$δuD2(t)$].

The measured quartz and deuterium shock front velocity perturbations are determined by subtracting the average shock state from the measured shock front velocity profiles or

and

It is important to ensure that the range of fluctuations are equivalent for $δuQz(t)$ and $δuD2(t)$⁠. We define the average quartz shock velocity to be

where t₀ is the time the shock first enters the quartz. To ensure that the range of fluctuations are equivalent, $u¯D2$ is defined as

where

In this routine, three free parameters are optimized: a time scaling parameter $(ΔtQzΔtD2)$⁠, an amplitude scaling parameter $(ΔuQzΔuD2)$⁠, and a DC offset parameter which allows the velocity average to be shifted by up to ±0.1 km/s $(δu)$⁠. The δu parameter was added to the fitting routine, to account for the random uncertainties associated with the average shock velocity which improved the goodness of fit. We look to minimize

The implementation of this routine is illustrated within Fig. 2. Using VISAR, we measure the quartz and deuterium shock front velocities shown as the blue and red lines, respectively, within Fig. 2(a). A time fitting range is defined, and a scaling factor is estimated a priori in order to determine the average shock state. The perturbation velocities are determined by subtracting the average shock state from the measured shock front velocities, as shown in Fig. 2(b) as the solid blue and red dashed line. Next, the parameters $ΔtQzΔtD2, ΔuQzΔuD2$⁠, and δu are determined by minimizing Eq. (10). The result of this minimization is shown in Fig. 2(b) as the solid red line.

Table I provides the experimental observables from this work. We report the average shock velocity in the quartz and D₂ as well as $ΔtD2ΔtQz$ for the average shock velocity in the deuterium and the quartz sample. These values are experimental observations and are not dependent upon an equation of state for their determination. Once these are determined, we then solve for the deuterium Lagrangian sound speed using Eqs. (2), (3), and (4) and the EOS of aluminum, quartz, and the D₂ pressure-density shock compressibility.

For small isentropic fluctuations about a shock state, the Grüneisen coefficient Γ can be determined by isentropic expansion of a state on the principal Hugoniot³¹ 

where $dP/dρ|H$ is the pressure derivative with respect to density along the Hugoniot. One can therefore use the sound speed for the shock compressed states to infer the Grüneisen coefficient using a fit to previous experimental data along the cryogenic Hugoniot data and the measured sound speed obtained from the perturbation analysis.

Various deuterium EOS models have been crafted over the years, either focused in planetary relevant pressure-temperature ranges or even expanding to describe extreme states of matter at which thermonuclear fusion is possible. We summarize basic information about the most accurate current wide-range EOS available to date in Table II. Briefly, one can distinguish between models that consider a collection of different species in chemical equilibrium, from models that are more or less directly fit to results from numerical methods such as density-functional-theory based molecular dynamics or path integral Monte-Carlo simulations. Considering the variety of approaches to describe the EOS of hydrogen over a broad range—orders of magnitudes—of pressure and temperatures, it is remarkable to see the overall global agreement between these models in Figs. 1, 3, and 4.

The evolution with increasing shock pressure of our Eulerian sound speed measurements along the Deuterium cryogenic Hugoniot is shown in Fig. 3 and reported in Table III. This illustrates that our data expand the regime previously explored experimentally^32,33 by a factor of four in shock pressure.

Figure 3 reveals a fluctuation in sound speed near 30–50 GPa. This fluctuation is a signature of molecular dissociation³⁵ which is not captured by all of the models. All of our experimental measurements occur above this transition. In addition, while all models essentially converge in their prediction for the cryogenic Hugoniot in the pressure-density plane (Fig. 1), the predictions for the sound speed at high shock pressure span a range of several km/s—larger than 10% in absolute value—near 200 GPa.

While our data do exhibit some scatter—consistent with the estimated experimental uncertainties—it is in good agreement with the majority of the EOS models up to ∼100 GPa as shown in Fig. 3. At higher pressure, the data tend to indicate that the models predicting the highest sound velocities (LEOS 1014)—that is a stiffer response for the deuterium—are less likely to be correct than the other models predicting a softer response such as LEOS 1017. Our two highest pressure measurements may suggest a slower increase in the sound speed with increasing shock pressure that is not currently captured in any of these models as shown in Fig. 3. Additional high-pressure sound speed measurements with more modern diagnostics and improved capabilities for laser pulse shaping, beam smoothing, and target fabrication would be worth performing to resolve this discrepancy.

Within Fig. 4, we compare our Grüneisen coefficient measurements with model predictions along the principal Hugoniot. The onset of the molecular dissociation near 30–50 GPa results in a drastic reduction in the Grüneisen coefficient (from values greater to 1 to less than 0.4 for all models) which then recovers upon completion of the transition to ∼0.6–0.8. In the high-pressure limits (>1000 GPa), several EOS models asymptote to the ideal gas limit $limP→∞Γ=2/3$⁠.

Interestingly, the predictions from the different models vary more dramatically—more than 20%—than the predictions for the shock compressibility or for the sound speed, making it easier to distinguish the different models of Fig. 4. For example, Saumon-2012 follows LEOS 1014 quite closely even though they are based on different assumptions and differs from Saumon-2007 which predicted an oscillatory behavior in this range. Similarly, while several of the EOS models (LEOS 1017, Caillabet, Hu, and HREOS) are based on similar density functional theory (DFT) based molecular dynamics simulations in this regime, they do not predict similar trends for the Grüneisen coefficient. This reveals the sensitivity of EOS derivative quantities to the assumptions or modeling choices made to build a broad range EOS table from a finite set of density functional theory-molecular dynamics (DFT-MD) data and should be seen as a cautionary tale regarding the validity of first-principles equation of states. Our Grüneisen coefficient data validate the increasing trend with increasing pressure but remain on the lower side of the prediction span. In particular, the highest pressure data above 100 GPa are found to follow the trends of LEOS 1014, Saumon-2012, and Hu, while they are inconsistent with several models such as LEOS 1017, Caillabet, and Saumon-2007. Our sound speed and Grüneisen parameter measurements suggest inaccuracies in the D₂ EOS that is currently used to model and help interpret most of the ICF implosion experiments at the NIF (LEOS 1014). These inaccuracies lay in the description of the second derivative of the free-energy surface, which is the derivative of the pressure-density compression curve along the principal Hugoniot.

Resolving these issues might prove important in improving the predictive capability of radiation-hydrodynamic simulations and in particular their ability to correctly capture the evolution of hydrodynamic instabilities³⁶ seeded at the interface between the ablator and the deuterium-tritium fuel (the sound speed of the deuterium-tritium fuel is a key parameter to govern the growth and evolution of these instabilities). In addition, improved EOS models might help mitigate the current difficulty to accurately reproduce the observed relative timing between the various shockwaves that are launched inside the capsule to achieve the quasi-isentropic spherical compression needed to achieve high-fusion yield implosions.^37,38 Our experimental data might prove instrumental in achieving a better understanding of thermonuclear fusion experiments at the NIF.

Overall, we find that the preferred deuterium EOS used to design and interpret many integrated ICF experiments at the NIF, LEOS 1014, shows significant disagreements with the experimental data for shock compressibility, sound speed, and Grüneisen coefficient. The present experimental data provide clear benchmarks for future improved EOS models which might contribute to better understanding of spherical implosion experiments producing higher thermonuclear fusion yield within the laboratory.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344.