Dissociation and ionization of warm dense hydrogen

In 1912, Langmuir [39] immersed a hot tungsten wire in a hydrogen atmosphere, and found that above 3000 K, heat was carried away from the wire at a rate much higher than would be expected by convection alone. The abnormally high conductivity appears because hydrogen molecules dissociate into atoms at that temperature, absorbing heat which is later released when the atoms recombine. At higher temperatures ($\sim$10000 K), hydrogen atoms separate into protons and electrons. Heavier atoms ionize at even higher temperatures, and Saha [40] proposed in 1920 that one could infer the temperature of stars from their relative concentration of ions.

For an equilibrium of ideal gases $\mathrm{C \rightleftharpoons A + B}$, we can write the dissociation fraction as a function of the gas temperature and density using the Saha equation (Appendix B, [40]):

$$\frac{f^{2}}{1-f} = 1.667 \times 10^{-4} \cdot \left(\frac{m_{A} ... ...rac{Z_{A} Z_{B}}{Z_{C}} \cdot \frac{T^{3/2}}{r_{s}^{3}} \exp(-\Delta E_{d}/k T)$$ (4.8)

where $f$ is the fraction of dissociated species, $\Delta E_{d}$ is the energy of dissociation of C, $m$ are the masses of the species in amu, $T$ is the temperature in Kelvin, and $Z$ are the vibrational-rotational-electronic partition functions. The density is characterized by the parameter $r_{s}$, so that each atom takes up a volume $4/3 \pi rs^{3}$. For hydrogen gas at one atmosphere and room temperature, $r_{s}$ is 86.

Figure 4.15: Previously proposed plasma phase transition [42,43] where hydrogen dissociates and ionizes simultaneously

Applying the Saha equation to the reactions $\mathrm{H_{2} \rightarrow H + H}$ and $\mathrm{H \rightarrow p^{+} + e^{-}}$, we find that in dilute gases, dissociation is a gradual process, not an abrupt transition (Figure 4.15). These reactions are entropy driven -- the bond dissociation energy of $\mathrm {H_{2}}$ is $\sim$50,000 K and the ionization potential of H is $\sim$150,000 K, yet dissociation and ionization occur at much lower temperatures, driven by the separation of one particle into two. Two particles take up more space than one, which means that by La Chatlier's principle, compressing a gas shifts the equilibrium toward association. Thus the temperature of dissociation and the temperature of ionization increase with increasing density.

This analysis shows that in dilute gases, dissociation and ionization are two separate events. At higher densities, however, it should become easier to ionize hydrogen, since as atoms are squeezed together, the band gap decreases. At extreme compressions ($r_{s} = 1$), hydrogen becomes metallic, and the electrons move freely as if in a uniform sea of background positive charge. This pressure ionization occurs even at absolute zero [44]. There has been speculation that at intermediate densities, the temperature needed to ionize hydrogen decreases with temperature, and at some point matches the temperature required to dissociate hydrogen molecules [43]. At such a plasma phase transition, hydrogen would simultaneously dissociate and ionize, and properties like pressure or conductivity could change abruptly with variations in temperature or pressure (depending on the order of the phase transition).

If a plasma phase transition existed, it could lead to the revision of astrophysical models [42] -- for example, giant planets like Jupiter have dense hydrogen near their core, and an abrupt phase change would change the way helium partitioned itself between molecular and metallic phases of hydrogen. Recently there has been a renewed interest in studying dense hydrogen, stemming from (1) the development of path-integral Monte Carlo methods to calculate hydrogen equations of state ab initio [45,50,46,47,48,49] (2) shock hugoniot experiments with gas guns [51], lasers [54,55], and exploding wires [52,53] able to access densities and temperatures near the postulated PPT ($\sim$15000 K, $r_{s} \approx 2$ bohr, according to a chemical model [43]. We demonstrate that the electron force field gives results consistent with the most recent high-level theory [49] and shock Hugoniot experiments [53]

In our simulations, we placed hydrogen molecules (64 nuclei) in a cubic periodic box, set atom velocities randomly from a Boltzmann distribution, then integrated the dynamics equations of motion with fixed volume and energy and a time step of 0.01 fs. We set the electron mass to be the same as the proton mass, making our simulation a plausible model of deuterium. We calculated the instantaneous temperature as the total kinetic energy of nuclei and electrons divided by $3/2 k T$; we computed electrostatic energies by the Ewald method; and we averaged thermodynamic data over 1 ps following a 200 fs equilibration period.

Holding density fixed ($r_{s}$ = 2 bohr) and performing simulations at a range of temperatures, we observed a thermal transition from a molecular to an atomic fluid (Figure 4.16). Proton-proton pair distribution functions plotted as a function of temperature (Figure 4.17) show a gradual transition between molecular and atomic fluid extremes, with an intermediate point at $T = 15400$ K, which compares well with the phase transition temperature of 15300 K estimated from a chemical model [43]. The pair distribution curves look similar to ones obtained using path-integral Monte Carlo [46], where an intermediate point (roughly estimated by looking at pair distribution curves by eye) occurs at 10000 K.

Figure 4.16: Dynamics snapshots showing deuterium dissociation as temperature is raised.

Figure 4.17: Proton-proton pair distribution function shows gradual dissociation.

We can also calculate pressure using the virial expression

$$P = \frac{2\braket{K_{total}}}{V} - \frac{1}{3 V} \left< \sum_{\textrm{all d.o.f.}} x_{i} \frac{\partial E}{\partial x_{i}} \right>.$$ (4.9)

We note that for a purely electrostatic system, we have the simple virial

$$\sum_{i<j} \vec{F}(\vec{r}_{ij}) \cdot \vec{r}_{ij} = \sum_{i<j} -\frac{q_{i} q_{j}}{r^{2}} \hat{r}_{ij} \cdot \vec{r}_{ij}$$ (4.10)

$$= -\sum_{i<j} \frac{q_{i} q_{j}}{r_{ij}}$$ (4.11)

$$= -U.$$ (4.12)

However, since our system includes Pauli exclusion forces as well, such a calculation gives the wrong result, whereas ours agrees with numerical differentiation of $-dE/dV$. We hold $r_{s}=2$ bohr and plot the equation of state, which we find agrees well with the chemical model at low temperatures and the QMC model at high temperatures, exactly the range of applicability each model is expected to have. We note that early path integral methods [50] predicted a first-order plasma transition with a negative slope $dP/dT$ at 10000 K; however, more accurate calculations with a more accurate nodal surface did not show any evidence of this negative slope [49]. Given that our method handles both low and high temperatures consistently, and we do not see any evidence of a negative $dP/dT$, we conclude that a first-order PPT does not exist in the temperature and density range considered.

Figure 4.18: Equation of state ($r_{s}$ = 2 bohr) shows good agreement with best available theory.

By looking at the distribution of electron sizes over the course of the simulation, we can estimate how many electrons become ionized. At low temperatures, we observe a Maxwell-Boltzmann-like distribution of electron sizes that broadens with increasing temperatures. At higher temperatures, we find that a small fraction of electrons escape and expand to be larger than the size of our periodic box; at that size, they no longer interact strongly with the rest of the system. Taking electrons with $r_{s} > 10$ bohr to be ionized, we find that some ionization occurs for $r_{s} = \mathrm{2.6\ bohr}, L_{box} = \mathrm{16.8\ bohr}$ at $\sim$25000 K and for $r_{s} = \mathrm{2.2\ bohr}, L_{box} = \mathrm{14.2\ bohr}$ at $\sim$30000 K but not for $r_{s} = \mathrm{2\ bohr}, L_{box} = \mathrm{12.9\ bohr}$. The ionization we observe is consistent with thermal ionization of hydrogen atoms in a dilute gas, where the electrons, not having a nucleus to associate with, expand to fill free space; in this regime, it is reasonable to expect the temperature required for ionization to increase with increasing density. However, further work needs to be done to determine whether metallic-like electrons appear at higher density and lower temperature. Metallic electrons in the electron force field would be characterized not by a large size, but by an increased mobility.

Figure 4.19: High densities suppress ionization that occurs at high temperature.

Experimentally, liquid deuterium can be compressed to near-metallic densities using shock waves generated by explosives, exploding wires, or lasers. In these experiments, the deuterium is compressed with a solid pusher; by measuring the position and acceleration of the pusher over a duration of nanoseconds, we deduce a density-pressure relation called a Hugoniot curve that is a characteristic of the material. From conservation of mass, energy, and momentum over the boundary of a shock wave, we know the internal energy, volume, and pressure must satisfy the Hugoniot relation

$$U –- U_{0} + \frac{1}{2}(V-V_{0})(P+P_{0}) = 0.$$ (4.13)

The Hugoniot curve measures how compressible liquid deuterium is to shock. In the last decade, there has been some controversy over compressibility, with laser driven experiments (Nova [54,55]) indicating a maximum compressibility of six times, and gas gun experiments indicating a lower compressibility of four times with a stiffer response. More recent experiments done with exploding wires ( [52,53]) support the stiffer response function. We would like to see what kind of Hugoniot our theory, which has only been parameterized to fit bond lengths of simple alkanes, would produce.

To estimate the Hugoniot curve using eFF, we carry out a series of simulations at fixed volume and different fixed temperatures; measure the pressure; plot the Hugoniot function versus pressure; and find by interpolation the pressure that makes the Hugoniot function zero. As a starting point, we compute a box of liquid hydrogen with $r_{s} = 3.16$ $(\rho = 0.171 \mathrm{g/cm^{3}})$, $T = 19.6 K$; we find $U_{0} = -0.477043$ hartrees/atom and $P = 0$.

Table 4.6: eFF computed Hugoniot curve.

$r_{s}$ (bohr)	$\rho$ (g/ $\mathrm{cm^{3}}$)	P (gPa)	T (K)
2.6	0.31	0.1	375
2.2	0.51	8.1	1768
2	0.68	34.2	6688
1.93	0.75	66.8	12464
1.86	0.84	151.0	23198
1.86	0.84	1216.5	344144

The eFF Hugoniot curve matches the curves obtained by the gas gun and Z machine, but not the ones obtained by the Nova laser. eFF reproduces the nearly vertical curve upward that path integral Monte Carlo shows; the vertical curve in eFF is the result of a Hugoniot function that was zero at two different pressures. It is believed that the true Hugoniot bends to the left slightly at high temperatures; if this is the case, we should be able to run the simulation at a slightly higher temperature, and have the Hugoniot curve be zero at only one point.

While PIMC shows a maximum density of $\sim$0.73 $\mathrm{g/cm^{3}}$ (compressibility of 4.3 times), eFF shows a maximum density of $\sim$0.84 $\mathrm{g/cm^{3}}$ (compressibility of 4.9 times). In contrast, the Nova laser Hugoniot shows a maximum compressibility of 6.0 times. eFF also shows a comparable rise in temperature to PIMC (ours is $\sim$1/3 less) over the course of its Hugoniot.

Another group has performed a WPMD simulation of hydrogen plasma, but used the earlier described Klakow potential, with an additional term added to capture electrostatic energy changes upon antisymmetrization. Their results agreed well with eFF at low temperatures, but deviated significantly at higher temperatures, with a higher compressibility (6.4 times) that better matches the Nova laser data. The same extra repulsion that prevents electrons from inappropriately coalescing in the eFF model may be serving to give a higher compressibility and better agreement with path-integral Monte Carlo over Pauli potentials based on Klakow's expression.

Figure 4.20: Reproduction of the experimental shock Hugoniot curve obtained by gas gun and Z machine; Nova laser remains an outlier.

To summarize, the electron force field shows that at $r_{s}=2$ bohr, dissociation from molecular to atomic fluid is gradual, and the equation of state and proton-proton pair distribution functions are consistent with path-integral Monte Carlo calculations. We observed no evidence of a plasma phase transition at this density. The good agreement with Hugoniot curves obtained from gas gun and Z machine experiments, as well as those obtained from PIMC, further confirms that we are describing the thermal transition from molecules to atoms correctly.

We have examined ionization as well, but only by measuring the concentration of free (large) electrons; we find that at low densities, the temperature required to create large electrons in free space increases with increasing density, as we would expect from a Saha model. Pressure ionization creates metallic electrons that are smaller yet highly mobile. It would be interesting to measure electron mobility at higher densities to determine if eFF can model metallic phases of hydrogen properly, and to better characterize pressure induced transitions from molecular to metallic hydrogen.