Shock Hugoniot

Recent Nova laser shock wave experiments on pre-compressed liquid deuterium (Da Silva, 1997; Collins et al., 1998) provided the first direct measurements of the high temperature equation of state of deuterium for pressures up to 330 GPa. It was found that deuterium has a significantly higher compressibility than predicted by the semi-empirical equation of state based on plasma many-body theory and lower pressure shock data (see SESAME model by Kerley (1983)). In an earlier series of experiments using the two-stage gas gun (Nellis et al., 1983; Holmes et al., 1995), pressures of up to 23 GPa were reached. The laser experiments are of particular importance to this work because they represent the only experimental data our PIMC simulation results can be directly compared to. The temperatures reached in gas gun experiments did not exceed 4400 K. Therefore a direct comparison with PIMC simulations is currently not feasible.

Figure: Density profile in an idealized shock experiment, in which a driving force moves the pusher at constant velocity $u_p$. The resulting shock wave travels at velocity $u_s$.

Shock wave experiments are an established technique (Zeldovich and Raizer, 1966) to determine the equation of state at high pressures and temperature, which has been applied to a wide range of materials including aluminum, iron, and water. In the experiment, a driving force is utilized to propel a pusher at constant velocity $u_p$ into a material at predetermined initial conditions ( $\varrho_0, p_0, T_0$) as shown in Fig. . The impact generates a planar shock wave, which travels at the constant velocity $u_s$ with $u_s>u_p$. The shock compression drives the material to a point on the principle Hugoniot, which is the locus of all final states characterized by ($\varrho,p,T$) that can be achieved by a single shock wave passing through. Under the assumption of an ideal shock (see below), the conservation laws of mass, momentum, and energy require only the measurement of the velocities of the pusher $u_p$ and the shock front $u_s$ in order to obtain an absolute equation of state data point. Pressure and density of the shock material are related to $u_s$ and $u_p$ by,

$$p-p_0 = \varrho_0 u_s u_p$$ (203)

$$\frac{\varrho}{\varrho_0} = \frac{u_s}{u_s-u_p} \quad.$$ (204)

The internal energy in the final state follows from the conservation laws,

$$H = E-E_0+\frac{1}{2}({V^{\!\!\!\!\!\!\:^\diamond}}-{V^{\!\!\!\!\!\!\:^\diamond}}_0)(p+p_0)=0 \quad,$$ (205)

which is called the Hugoniot relation. It can also be used to determine the Hugoniot curve analytically from a given equation of state.

By an ideal shock, one means that a planar pusher is driven at constant velocity into the sample. The resulting shock wave is characterized by a planar shock front that travels at constant velocity during the measurement. Furthermore, one assumes that unshocked material remains at known initial conditions and not preheated as for example by x-rays created at the laser target interaction. Under these assumptions, one can determine the equation of state from the measured velocities $u_s$ and $u_p$ using Eqs. -.

Figure: Comparison Hugoniot function calculated with PIMC simulations of different accuracy using $n_{\rm A}=1$, and $n_{\rm E}=2$: FP nodes with $N_P$=32 ($\triangle$ for $\tau^{-1}=10^6\rm K$ reported by Militzer et al. (1998), $\rhd$ for $\tau^{-1}=2 \cdot 10^6 \rm K$, $\bigtriangledown$ for $\tau_F^{-1}=8 \cdot 10^6 \rm{K}$ and $\tau_B^{-1}=2 \cdot 10^6 \rm{K}$) and $N_P$=64 ($\Box$ for $\tau^{-1}=2 \cdot 10^6 \rm{K}$) as well as with VDM nodes and $N_P$=32 ($\circ$ for $\tau^{-1}=10^6\rm{K}$ and $\bullet$ for $\tau^{-1}=2 \cdot 10^6 \rm{K}$). Beginning at high pressures, the points on each Hugoniot correspond to the following temperatures $125\,000, 62\,500, 31\,250, 15\,625,$ and $10\,000\,\rm{K}$. The dashed line corresponds to a calculation using the VDM alone.

In the recent laser shock experiments, a shock wave is propagating through a sample of pre-compressed liquid deuterium characterized by an initial state, ($E_0$,  ${V^{\!\!\!\!\!\!\:^\diamond}}_0$, $p_0$) with $T=19.6\,\rm {K}$ and $\varrho_0=0.171\,\rm {g/cm^3}$. In our calculations, we set $E_0$ to its exact value of $-15.886\rm {eV}$ per atom (Kolos and Wolniewicz, 1964) and $p_0 = 0$ because $p \gg p_0$. Using the PIMC simulation results for $p$ and $E$, we calculate $H(T,\varrho)$ from Eq.  and then interpolate $H$ linearly at constant $T$ between the two densities corresponding to $r_s=1.86$ and $2$ to obtain a point on the Hugoniot in the $(p,\varrho)$ plane. Results at $r_s=1.93$ confirm that the function is linear within the statistical errors. The PIMC data for $p$, $E$, and the Hugoniot are given in Tab. .

Table: Pressure $p$ and internal energy per atom $E$ from PIMC simulations with $32$ pairs of electrons and deuterons. For $T \ge 250\,000\,\rm{K}$, we list results from simulations with FP nodes and $\tau_F^{-1}=8 \cdot 10^6\,\rm{K}$ and $\tau_B^{-1}=2 \cdot 10^6\,\rm{K}$, otherwise with VDM nodes and $\tau^{-1}=2 \cdot 10^6\,\rm{K}$.

$T (\rm {K})$	$p (\rm {Mbar})$	$E (\rm {eV})$	$p (\rm {Mbar})$	$E (\rm {eV})$	$\varrho^{\rm Hug} (\rm {gcm}^{-3})$	$p^{\rm Hug} (\rm {Mbar})$
	$r_s=2~~$	$r_s=2~~$	$r_s=1.86$	$r_s=1.86$
1000000	53.79 (5)	245.7 (3)	66.85 (8)	245.3 (4)	0.7019 (1)	56.08 (5)
500000	25.98 (4)	113.2 (2)	32.13 (5)	111.9 (2)	0.7130 (1)	27.48 (4)
250000	12.12 (3)	45.7 (2)	14.91 (3)	44.3 (2)	0.7242 (1)	12.99 (2)
125000	5.29 (4)	11.5 (2)	6.66 (2)	11.0 (1)	0.7300 (3)	5.76 (2)
62500	2.28 (4)	-3.8 (2)	2.99 (4)	-3.8 (2)	0.733 (1)$\:\:$	2.54 (3)
31250	1.11 (6)	-9.9 (3)	1.58 (7)	-9.7 (3)	0.733 (3)$\:\:$	1.28 (5)
15625	0.54 (5)	-12.9 (3)	1.01 (5)	-12.0 (2)	0.721 (4)$\:\:$	0.68 (4)
10000	0.47 (5)	-13.6 (3)	0.80 (8)	-13.2 (4)	0.690 (7)$\:\:$	0.51 (5)

In Fig. , we compare the effects of different approximations made in the PIMC simulations such as time step $\tau$, number of pairs $N_P$ and the type of nodal restriction. For pressures above 3 Mbar, all these approximations have a very small effect. The reason is that PIMC simulation become increasingly accurate as temperature increases. The first noticeable difference occurs at $p \approx 2.7 \rm {Mbar}$, which corresponds to $T=62\,500\,\rm {K}$. At lower pressures, the differences become more and more pronounced. We have performed simulations with free particle nodes and $N_P=32$ for three different values of $\tau$. Using a smaller time step makes the simulations computationally more demanding and it shifts the Hugoniot curves to lower densities. These differences come mainly from enforcing the nodal surfaces more accurately, which seems to be more relevant than the simultaneous improvements in the accuracy of the action $S$, that is the time step is more constrained by the Fermi statistics than it is by the potential energy. We improved the efficiency of the algorithm by using a smaller time step $\tau_F$ for evaluating the Fermi action than the time step $\tau_B$ used for the potential action. Unless specified otherwise, we used $\tau_F=\tau_B=\tau$. At even lower pressures not shown in Fig. , all of the Hugoniot curves with FP nodes turn around and go to low densities as expected.

As a next step, we replaced the FP nodes by VDM nodes. Those results show that the form of the nodes has a significant effect for $p$ below 2 Mbar. Using a smaller $\tau$ also shifts the curve to slightly lower densities. In the region where atoms and molecules are forming, it is plausible that VDM nodes are more accurate than free nodes because they can describe those states (see chapter 3). We also show a Hugoniot derived on the basis of the VDM alone (dashed line). These results are quite reasonable considering the approximations (Hartree-Fock) made in that calculation. Therefore, we consider the PIMC simulation with the smallest time step using VDM nodes ($\bullet$) to be our most reliable Hugoniot. Going to bigger system sizes $N_P=64$ and using FP nodes also shows a shift towards lower densities.

Figure: Comparison of experimental and several theoretical Hugoniot functions. All three PIMC curves were calculated with $\tau^{-1}=2 \cdot 10^6\,\rm{K}$, $n_{\rm A}=1$, $n_{\rm E}=2$, and $32$ pairs of electrons and deuterons. They were obtained with free particle nodes, VDM nodes, and VDM nodes plus the corrections from Eq. .

Figure: Logarithmic Hugoniot graph as in Fig. . including the gas gun experiments by Holmes et al. (1995) and Nellis et al. (1983).

Fig. compares the Hugoniot from laser shock wave experiments (Da Silva, 1997; Collins et al., 1998) with PIMC simulations (VDM nodes, $\tau^{-1}=2 \cdot 10^6\,\rm {K}$) and several theoretical approaches: SESAME model by Kerley (1983) (thin solid line), linear mixing model (dashed line) by Ross (1998), DFT-MD by Lenosky et al. (2000) (dash-dotted line), Padé approximation in the chemical picture (PACH) by Ebeling and Richert (1985a) (dotted line), and the work by Saumon and Chabrier (1992) (thin dash-dotted line).

Figure: Shock velocity vs. pusher velocity as directly measured in the shock wave experiment and the comparison with estimates from PIMC simulation with VDM nodes.

The differences of the various PIMC curves in Fig.  as well as in Fig.  are small compared to the deviation from the experimental results by Da Silva (1997) and Collins et al. (1998). One finds that the corrections from Eq.  have only a small effect on the Hugoniot. In the experiments, an increased compressibility with a maximum value of $6 \pm 1$ was found while PIMC predicts $4.3 \pm 0.1$, only slightly higher than that given by the SESAME model. Only for $p>2.5 \rm Mbar$, does our Hugoniot lie within experimental error bars. In this regime, the deviations in the PIMC and PACH Hugoniot are relatively small, less than $0.05 \, \rm gcm^{-3}$ in density. In the high pressure limit, the Hugoniot goes to the FP limit of 4-fold compression. This trend is also present in the experimental findings. For pressures below 1 Mbar, the PIMC Hugoniot goes back to lower densities and shows the expected tendency towards the experimental values from earlier gas gun work Nellis et al. (1983); Holmes et al. (1995) and lowest data points from Da Silva (1997); Collins et al. (1998). This trend can be studied best in the logarithmic graph shown in Fig , where we also included our lowest available pressure point on the Hugoniot, which was obtained from simulations with 32 pairs of electrons and deuterons and the time step $\tau^{-1}_F = 8 \cdot 10^6\,\rm K$ and $\tau^{-1}_B = 2 \cdot 10^6\,\rm K$. Within the statistical error bars, the PIMC Hugoniot curve tends towards the results from gas gun experiments. For these low pressures, one also finds that the differences between PIMC and DFT-MD are also relatively small compared to the deviation from the laser shock data.

Using the PIMC equation of state, one can also determine the shock and pusher velocity on the Hugoniot. From Eqs.  and , one finds,

$$m_D u_s^2 = \frac{\varrho}{\varrho-\varrho_0} \frac{p}{\varrho_0}\quad,$$ (206)

$$m_D u_p^2 = p\left( \frac{1}{\varrho_0} - \frac{1}{\varrho}\right)\quad.$$ (207)

The results are shown in Fig.  and compared to the experimental shock and pusher velocities published in (Collins et al., 1998). First of all, one finds the differences between theory and experiment are not as pronounced as in the $p$-$\varrho$ graph in Fig. . This fact simply follows from Eq.  where one divides by the difference of $u_s$ and $u_p$. It means that one needs a high accuracy in the measurements of both velocities, since error bars will increase substantially when the density is determined. However, it is not obvious why the PIMC and DFT-MD results are in within the corners of the experimental $u_s$-$u_p$ error bars. A simple error propagation of these $u_s$ and $u_p$ error bars leads to much bigger error bars in the density than those reported by Da Silva (1997) and Collins et al. (1998) and shown in Fig. . Possibly the experiments allow a more accurate determination of the difference $u_s-u_p$ than of the individual velocities.

Summarizing, one can say that PIMC simulations predict a slightly increased compressibility of $4.3 \pm 0.1$ compared to the SESAME model but they cannot reproduce the experimental findings of values of about $6 \pm 1$. Further theoretical and experimental work will be needed to resolve this discrepancy.