Finite element simulation of the liquid-liquid transition to metallic hydrogen

TL;DR

This paper presents a finite element simulation of the liquid-liquid transition to metallic hydrogen, confirming experimental signatures and suggesting complex transition dynamics, while proposing experiments to distinguish optical property changes.

Contribution

The study introduces a detailed finite element model that accurately simulates laser-heated hydrogen and analyzes the phase transition signatures, highlighting the need for large latent heat values and complex transition dynamics.

Findings

Simulated heating curves match experimental plateaus indicating phase transition.

Large latent heat values are required to explain observed signatures.

Proposed experiments can differentiate between metallic transition and bandgap closure.

Abstract

Hydrogen at high temperature and pressure undergoes a phase transition from a liquid molecular phase to a conductive atomic state, or liquid metallic hydrogen, sometimes referred to as the plasma phase transition (PPT). The PPT phase line was observed in a recent experiment studying laser-pulse heated hydrogen in a diamond anvil cell in the pressure range $\sim 100 - 170 \text{GPa}$ for temperatures up to $\sim 2000 \text{K}$. The experimental signatures of the transition are (i) a negative pressure-temperature slope, (ii) a plateau in the heating curve, assumed to be related to the latent heat of transformation, and (iii) an abrupt increase in the reflectance of the sample. We present a finite element simulation that accurately takes into account the position and time dependence of the heat deposited by the laser pulse. We calculate the heating curves and the sample reflectance and transmittance. This simulation confirms that the observed plateaus are related to the phase transition, however we find that large values of latent heat are needed and may indicate that dynamics at the transition are more complex than considered in current models. Finally, experiments are proposed that can distinguish between a change in optical properties due to a transition to a metallic state or due to closure of the bandgap in molecular hydrogen.

Tables

Table 1. Table 1: Reflectance coefficients R 𝑅 R of the non-active boundaries. Transmittance can be found using T ​ r = 1 − R 𝑇 𝑟 1 𝑅 Tr=1-R . The indices refer to the interfaces as listed in the text.

Boundaries	$R_1$	$R_2$	$R_4$
$\lambda =500$nm	0.174	0.011	0.174
$\lambda =1064$nm	0.168	0.003	0.168

Table 2. Table 2: The values for the material parameters used in the simulation. A dash means the parameter is not used in the simulation. Stressed diamond is diamond with a tenfold lower thermal conductivity.

	H2	MH	Tungsten	Al2O3	Diamond	Stressed diamond
$\rho$ $\left(\frac{kg}{m^3}\right)$	850	850	25100	5000	4400	4400
$k$ $\left(\frac{W}{mK}\right)$	20	20	-	20	2200	220
$\rho c$ $\left(\frac{MJ}{m^{3}K}\right)$	10	20	3.5	1.2	1.7	1.7
$\tilde{N}\left(\lambda =1064\text{nm}\right)$	2.31	5.22+6.48i	3.12+3.67i	2.16	2.39	2.39
$\tilde{N}\left(\lambda =500\text{nm}\right)$	2.46	2.94+4.58i	3.31+2.74i	2.20	2.43	2.43
$\rho L$ $\left(\frac{MJ}{m^3}\right)$	1700		-	-	-	-
$T_m\left(K\right)$	1500		-	-	-	-
Sources	Rillo2019 ; Hanfland1994 ; Evans1998 ; Morales2010 ; Moroe2011		Haynes2015 ; Raju1997	Malitson1963 ; Duan1999	Haynes2015 ; Olson1993	-

Equations

$$\frac{\partial q}{\partial t}=-\frac{\partial S}{\partial z}+\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right).$$

$$\Delta q_i^{(n)}=\Delta t\left[S^{(n)}\left(z_i-\frac{\Delta z_i}{2}\right)-S^{(n)}\left(z_i+\frac{\Delta z_i}{2}\right)+\frac{T_{i+1}^{(n)}-T_i^{(n)}}{r_{i,i+1}^{(n)}}-\frac{T_i^{(n)}-T_{i-1}^{(n)}}{r_{i-1,i}^{(n)}}\right].$$

$$r_{i,i+1}^{(n)}=\frac{z_{i+1}-z_i}{k_{i,i+1}^{(n)}},$$

$${\tilde{N}}_{\text{eff}}^{\frac{2}{3}}=f{\tilde{N}}_{\text{MH}}^{\frac{2}{3}}+(1-f){\tilde{N}}_{{\text{H}}_2}^{\frac{2}{3}}.$$

$$E(z)$$

$$E_j(z)$$

$$\left\{\begin{array}[]{rcl}a_{j}e^{i\kappa_{j}d_{j}}+b_{j}e^{-i\kappa_{j}d_{j}}&=&a_{j+1}+b_{j+1}\\ a_{j}\kappa_{j}e^{i\kappa_{j}d_{j}}-b_{j}\kappa_{j}e^{-i\kappa_{j}d_{j}}&=&a_{j+1}\kappa_{j+1}+b_{j+1}\kappa_{j+1}\end{array}\right.\hskip 5.0pt.$$

$$\left\{\begin{array}[]{l}a_{1}=1\\ b_{N}=0\end{array}\right.\hskip 10.0pt\text{or}\hskip 10.0pt\left\{\begin{array}[]{l}a_{1}=0\\ b_{N}=1\end{array}\right.\hskip 5.0pt,$$

$$S_{j}^{(n)}(z)=-I_{0}(t_{n})\times\left[\begin{array}[]{cc}&\frac{n_{j}}{n_{1}}\left(|a_{j}|^{2}e^{-\frac{4\pi}{\lambda}K_{j}(z-Z_{j})}-|b_{j}|^{2}e^{\frac{4\pi}{\lambda}K_{j}(z-Z_{j})}\right)\\ &-\frac{K_{j}}{n_{1}}2|b_{j}^{*}a_{j}|\sin\left(\frac{4\pi}{\lambda}n_{j}(z-Z_{j})+\text{arg}(b_{j}^{*}a_{j})\right)\\ \end{array}\right].$$

$$I_0(t)=\frac{2At}{{\tau }^2}{e}^{-{\left(\frac{t}{\tau }\right)}^2},$$

$$Tr=\frac{n_N}{n_1}|a_N{|}^2,R=|b_1{|}^2.$$

$$k_{i,i+1}=k_{{\text{H}}_2}^{1-<f>}{k}_{\text{MH}}^{<f>},<f>=\frac{f_i+f_{i+1}}{2}.$$

$$r_{{\text{Al}}_2{\text{O}}_3}=\frac{d}{k_{{\text{Al}}_2{\text{O}}_3}}+r_{\text{bound},1}+r_{\text{bound},2}.$$

$$\Delta T_i^{(n)}:=T_i^{(n+1)}-T_i^{(n)}=\frac{\Delta q_i^{(n)}}{\rho c\Delta z_i}.$$

$${\rho }_{\text{eff}}{c}_{\text{eff}}=f{\rho }_{\text{MH}}{c}_{\text{MH}}+(1-f){\rho }_{{\text{H}}_2}{c}_{{\text{H}}_2}.$$

$$\Delta f_i^{(n)}:=f_i^{(n+1)}-f_i^{(n)}=\frac{\Delta q_i^{(n)}}{\rho L\Delta z_i},$$

Full text

I Introduction

Metallic Hydrogen (MH) was studied by Wigner and Huntington Wigner1935 in 1935; they predicted that compressing solid molecular hydrogen (H2) to a pressure (P) of 25 GPa or higher would destabilize the molecular bond and cause a phase transition from an insulating solid molecular state to a solid atomic metallic state. In its metallic form, hydrogen is expected to have many interesting properties, amongst which high-temperature superconductivity stands out Ashcroft1968 . After Wigner and Huntington’s prediction it became clear that the transition pressure lies at much higher pressure McMahon2012 ; McMinis2015 . The Wigner-Huntington transition has recently been observed at 495 GPa Dias2017 . The subject of this paper is a second pathway (Pathway II) to metallic hydrogen shown in Fig. 1. This transition is at lower pressures and is achieved at high temperatures due to thermal dissociation to form liquid atomic metallic hydrogen (LMH).

The great outer planets such as Jupiter are mainly composed of hydrogen. Descending through Jupiter’s atmosphere, the pressure and temperature (T) increase until a transition takes place to liquid metallic hydrogen. For planetary modeling the location of the transition temperature and pressure is very important Guillot2005 . An early theory based on a chemical model predicted high temperatures for the transition Saumon1992 . Density functional theories predicted much lower temperatures, and pressures in the 100-200 GPa range for the phase transition line Scandolo2003 ; Lorenzen2010 ; Morales2010 . The common denominator of these theories is a first-order phase transition from liquid molecular to liquid atomic metallic hydrogen with latent heat of transformation, a phase line with a negative P,T slope with increasing pressure, and a critical point (Fig. 1), such that the transition is continuous for lower pressures and higher temperatures. However, DFT-based predictions had a large scatter of values because the results depend strongly on the choice of the density functional. An alternative route using quantum Monte-Carlo simulation has also been giving converging results recently Pierleoni2016 ; Mazzola2018 .

The ideal experiment would have a thermally isolated sample in a pressure cell and add heat to the sample to warm it through the transition, enabling a determination of the transition P and T and the latent heat. In this case a plot of the temperature vs. heat added would have a rising slope with a plateau at the temperature of the transition to enable a direct measurement of the latent heat. There are two ways of achieving the high pressures and temperatures required to observe the PPT: dynamic (shock) compression and static compression in a diamond anvil cell (DAC). Neither of these techniques provide thermal isolation of the sample. In dynamic experiments P and T conditions last for tens of nanoseconds; in DACs high pressures can be maintained continuously. DACs cannot maintain their integrity if continuously heated to the high temperatures of the transition ($\sim$1000-2000 K), so pulsed heating is utilized to heat the sample.

Liquid metallic hydrogen was first observed in dynamic experiments using reverberating shock waves Weir1996 , however the measurement was not sensitive to determining the phase line. In a single-shock Hugoniot such as in Ref. Loubeyre2012 the temperature rise can be $\sim$10000-20000 K. A recent dynamic experiment on deuterium at the National Ignition Facility Celliers2018 investigated the transition at temperatures under 2000K. These dynamic experiments either did not probe the phase transition or did not determine the transition temperature. In DACs one can measure both pressure and temperature (by fitting blackbody radiation to a Planck curve). By using pulsed lasers to heat the sample for short periods of time (few hundred nanoseconds) a standard DAC can be used in which mainly the sample is heated and not the body of the DAC itself, enabling high static pressures and temperatures of thousands of degrees K.

The technique of pulsed laser heating led to the recent experimental determination of the phase line of the PPT for pressures in the range of 100-170 GPa and temperatures up to  2000 K by Zaghoo, Salamat, and Silvera (ZSS) Zaghoo2016 ; these studies were further refined to determine the conductivity, degree of dissociation, and Drude parameters that lead to the dielectric function Zaghoo2017 . In the present paper we focus on this pathway to LMH. We note that the DAC is not a thermally isolated system and in the geometry used in the experiment, the energy from the pulsed laser heated sample can both go into heating the sample as well as diffuse into the diamonds (Fig. 2a). Nevertheless a plateau was found in the heating curves; below we shall show that this can be simulated. The metallization transition changes the optical properties as the sample becomes absorbing and reflective. Detecting this change and fitting it to metallic optical conductivity has been the method used to detect the transition to the metallic phase. In particular, in Zaghoo2016 , the PPT is seen to lead to an abrupt increase in the reflectance, with frequency dependence in qualitative agreement with conductivity predicted by a Drude free-electron model of a metal. Additionally, the heating curve relating temperature to applied laser power is accompanied by the appearance of a plateau, indicating the presence of latent heat of transition. Several authors Geballe2012 ; Knudson2015 ; Goncharov2017 ; Celliers2018 interpreted the plateau differently, instead relating it to a continuous change in optical parameters due to temperature dependent bandgap closure.

In the experiment the properties change rapidly in time. The pressurized molecular hydrogen in the DAC is heated by a pulsed laser beam (280 ns pulse width), partially absorbed by a thin tungsten layer deposited on one of the diamonds (Fig. 2). The heated tungsten film heats the hydrogen pressed to its surface, creating a layer of LMH when sufficiently hot. This layer appears, thickens, and then shrinks as the pulse cycle is traversed. Increased laser power leads to thicker films of LMH growing at the LMH-molecular H2 interface. Continuous wave laser light shines on the sample to measure transmission and reflection during the heating cycle using time-resolved spectroscopy. The optical properties are strongly time-dependent and moreover feed back into the position dependent absorption of the heating laser pulse, as the LMH itself also absorbs the laser power.

In this paper we simulate these dynamics with a finite element analysis (FEA) of the pulsed-laser heated DAC. Previous simulations Geballe2012 ; Montoya2012 investigated the case of a thicker ($\sim\mu$m) laser absorber embedded in an insulating material, with double sided heating; the absorber itself developing a steep thermal gradient. In these cases the absorber was also the sample. In the experiment of ZSS the absorber is only several nanometers thick, heats uniformly and heats the sample. A recent FEA by Goncharov and Geballe Goncharov2017 attempting to simulate the experiment of ZSS did not find a plateau. Their simulation did not include thermodynamic properties of LMH, only molecular hydrogen so that it was not designed to simulate the experiment and reported non-physical results Silvera2017 .

The goal of the present paper is to provide a finite element simulation that better corresponds to the experimental configuration of ZSS Zaghoo2016 . We show that it is important to take into account the spatial and temporal dependence of the heat delivered to the cell by the laser pulse. When this improvement is taken into account, we find that plateaus can appear in the heating curve, due to the latent heat of transformation. The temperature at which the heating curve develops the plateau is essentially the metallization temperature and can be used as an estimate for this temperature.

The structure of this paper is as follows. In Sec.II, we outline the numerical simulation method, consisting of three computation steps per time step: the laser-pulse heating, the heat conduction, and the temperature change. The results are presented in Sec.III for the heating curves and for the reflectance and transmittance. In Sec.IV these results are discussed and compared to other simulations and to the experiment. The importance for further experiments is discussed in Sec.V. We conclude in Sec.VI.

II Numerical methods

The experimental setup of Ref. Zaghoo2016 is shown in Fig. 2a. The sample chamber is about 20 $\mu$m in diameter and has cylindrical symmetry. It consists of multiple thin layers, as modeled and illustrated schematically in Fig. 2b. The largest portion of the sample chamber consists of a 2 $\mu$m thick layer of molecular hydrogen, transparent to the heating laser. The diamond culets are covered with a 50 nm layer of amorphous alumina (Al2O3); an $\sim 8$ nm layer of tungsten (W) is deposited on one of the diamonds, which is covered by a 5 nm Al2O3 layer. The transparent alumina coatings inhibit the diffusion of hydrogen into the diamonds and the tungsten. The partially transmitting tungsten layer acts as an absorber of laser light. A pulsed heating laser (with pulse duration $\tau=280$ ns and wavelength $\lambda=1064$ nm) is incident on the top (hydrogen side) to heat the sample. Upon heating to the metallization temperature, part of this molecular hydrogen layer converts into metallic hydrogen.

Probe laser light (wavelengths of 514 nm, 633 nm, 808 nm, and 980 nm) is continuously shone into the cell from the top (hydrogen side) to measure the transmittance and reflectance. Both the heating and probe lasers are slightly defocused, so that the light intensity can be assumed roughly uniform along the surface of the cell. Finally, the temperature is determined by collecting the blackbody radiation from the top of the DAC and fitting it to a Planck curve to enable determination of the peak temperature in a pulse Rekhi2003 .

In our numerical simulation, we use a one-dimensional (1D) model that incorporates the different layers with dimensions shown in Fig. 2b. We also include 100 $\mu$m of the bottom diamond in our simulation and assumed, as our boundary conditions, that the rest of the (2 mm thick) diamond, as well as the top diamond, remains at room temperature during the heating pulse. Since our model is 1D, neither the conductive energy loss due to the gasket nor the Gaussian profile of the heating laser can be taken into account. However, both the sample chamber and the laser profile are wide in the radial direction, meaning the 1D model acts as a good first approximation able to capture the essential physics.

To describe the thermodynamics of the system, the one-dimensional form of the diffusion equation was used. We define the z-axis to be normal to the plane of the sample. Call $q(z,t)$ the local energy density ($J/m^{3}$), $T(z,t)$ the temperature, $S(z,t)$ the Poynting vector of the heating laser in the positive z-direction, and $k$ the thermal conductivity of the material. Then, the continuity equation implies:

[TABLE]

To proceed, this equation was discretized in the z- and t-domains. A uniform temporal grid $t_{n}=n\Delta t$ and a non-uniform spatial grid $z_{i}$ with element thicknesses $\Delta z_{i}$ were chosen. Replacing the energy density per unit volume $q(z,t)$ by the energy density per unit surface $q_{i}^{(n)}=q(z_{i},t_{n})\Delta z_{i}$, and defining $\Delta q_{i}^{(n)}:=q_{i}^{(n+1)}-q_{i}^{(n)}$, equation (1) becomes:

[TABLE]

Here a superscript (n) was used to indicate a variable is evaluated at time $t_{n}$, and $r^{(n)}_{i,i+1}$ denotes the thermal resistance between elements $i$ and $i+1$, given by:

[TABLE]

with $k^{(n)}_{i,i+1}$ the thermal conductivity of the material between elements $i$ and $i+1$.

As an amount of heat $\Delta q_{i}^{(n)}$ is added to the element $i$, one of two things can happen. First, the temperature in the element can increase. Second, if element $i$ is a hydrogen element, the hydrogen molecules can dissociate, creating metallic hydrogen. To model this, a metallization fraction $f_{i}^{(n)}$ between [math] and $1$ is assigned to the hydrogen elements: $f_{i}^{(n)}=0$ means the element is in the molecular state; when $f_{i}^{(n)}=1$ it is metallic.

For our grid, we chose elements of thickness $\Delta z_{i}=1$nm in the hydrogen and alumina layers, and a non-uniform grid of thicknesses $\Delta z_{i}>5$nm in the lower diamond. The tungsten is considered to be its own element with uniform temperature. Both doubling and halving the number of elements did not change the results.

We start with a uniform temperature of $T_{j}^{(0)}=300$K and no MH: $f_{j}^{(0)}=0$. Every time step in the simulation is divided in three parts. First, the Poynting vector $S^{(n)}(z)$ is determined in every element from the optical properties of each element. Second, the thermal conductivity between adjacent elements and the contribution of heat conduction to $\Delta q_{i}^{(n)}$ is calculated. Third, $\Delta q_{i}^{(n)}$ is calculated from (2) and is used to compute the resulting new temperature and metallization fraction for each element. The simulation proceeds in time by repeating this process with the new values $T_{i}^{(n+1)}$ and $f_{i}^{(n+1)}$, which allows us to calculate $\Delta q_{i}^{(n+1)}$, and so on. In the following subsections, each of these parts is described in detail. The values for the material parameters are listed in the final subsection.

II.1 Heat absorption from the laser pulse: calculation of the Poynting vector

The Fresnel formalism is suitable to describe the optical properties of any number $N$ of thin layers (thinner than the wavelength $\lambda$). The layers in the chosen geometry are shown in Fig. 2b. The complex index of refraction for layer $j$ is denoted by $\tilde{N}_{j}=n_{j}+\textrm{i}K_{j}$, and $\kappa_{j}=\frac{2\pi}{\lambda}\tilde{N}_{j}$ is the corresponding complex wavenumber. Furthermore, $d_{j}$ are the thicknesses of the layers, and $Z_{j}$ are the positions of the boundaries between the layers. The thickness of the layer of metallic hydrogen is found from the metallization fractions $f_{i}^{(n)}$ mentioned above. In our simulation, we always find a sharply bounded region where $f_{i}^{(n)}=1$, so the layer thickness of metallic hydrogen can be clearly determined. For a mixture of molecular and metallic hydrogen, the following complex index of refraction was used LandauLifshitz :

[TABLE]

The laser is assumed to be incident along the normal, so the electric field is tangential to the interfaces between the layers. The magnitude of the electric field throughout the sample is given by the piecewise function

[TABLE]

with unknown coefficients $a_{j}$ and $b_{j}$. The boundary conditions are such that $E$ and $\partial E/\partial z$ are continuous across the boundaries. Applied to (6), this gives the following $2N-2$ equations:

[TABLE]

Imposing these, along with the boundary conditions

[TABLE]

depending on the incident direction of the laser, gives a system with a unique solution for all $a_{j}$ and $b_{j}$, which is found numerically. Calculating the time-averaged Poynting vector and normalizing it to the time-dependent intensity of the laser pulse $I_{0}(t)$ then gives the laser power profile as a piecewise function in space:

[TABLE]

Note that $S_{j}(z)$ is constant if $K_{j}=0$, meaning that no power is lost in non-absorbing materials as expected. The Poynting vector is a continuous function of $z$ because of the imposed electromagnetic boundary conditions. For the time-dependent laser pulse intensity we use the following function, that simulates the pulses used in the experiment:

[TABLE]

where $A$ is the pulse energy (per unit area) and $\tau$ the full width at half maximum (FWHM) time. The transmittance $Tr$ and reflectance $R$ of this multilayered system are calculated by comparing the transmitted and reflected laser power to the incident laser power assuming the outermost layers are not absorbing ($K_{1}=K_{N}=0$). We get

[TABLE]

The transmission and reflection are calculated at every point in time for a wavelength of 500 nm.

For thickness variations comparable to the wavelength of light, the Fresnel equations lead to Fabry-Pérot type interferences. However, this is an artifact related to the assumption of coherent, collimated illumination. Since the heating laser is not collimated, the cell chamber should not be treated as a Fabry-Pérot interferometer. To correct for this, the sample chamber is divided into the following subsystems:

The air-diamond interface at the top 2. 2.

The alumina coated boundary between the top diamond and the hydrogen 3. 3.

The “active region”, starting from the H2-MH interface down to the interface between alumina and the lower diamond, and 4. 4.

The diamond-air interface at the bottom.

Absorption takes place only in the active region. This region is treated within the Fresnel formalism. The remaining boundaries are treated in the incoherent approximation: only the intensity changes from reflection and transmission are taken into account. The reflection coefficients used for the boundaries outside the active region are fixed, and are listed in table 1. The (multiple) reflections between the boundaries outside the active region lead to power influx on the active region from both sides. Using the previously outlined Fresnel method on the active region with this power influx as input, the heat deposited in each element can be calculated.

II.2 Heat redistribution from thermal diffusion: calculation of the thermal resistance

If elements $i$ and $i+1$ are associated with the same material, the thermal resistance can be calculated with formula (3) where $k_{i,i+1}:=k$ is the constant thermal conductivity of that material. For the hydrogen elements, we need to be more careful, since in general, the thermal conductivity of hydrogen depends on its state. Therefore, the following interpolation is used Woodside1961 :

[TABLE]

The above method works for all pairs of neighboring elements except those where either $i$ or $i+1$ is the tungsten element. The thermal resistance to the tungsten element from either side is taken as the thermal resistance of the respective alumina layer. This includes any thermal boundary resistance. In general, the thermal resistance of an alumina layer with thickness $d$ can be written as follows:

[TABLE]

It is possible that the thermal boundary resistances $r_{\text{bound},1}$ and $r_{\text{bound},2}$ cannot be neglected, since the thermal boundary resistance of amorphous alumina to platinum Cappella2013 and to aluminum Hopkins2007 is of the same order as the bulk thermal resistance of the alumina layers.

II.3 Change in temperature and metallization fraction

Let $\Delta q_{i}^{(n)}$ be the total heat generated in element $i$. This is the sum of the heat delivered by the laser pulse described in subsection II.1 and the heat conducted into the element as discussed in subsection II.2. From the mass density $\rho$ and the specific heat per unit mass $c$ of the material, the specific heat per unit volume $\rho c$ is calculated. From its definition, it follows that the change in temperature is given by:

[TABLE]

For elements with a mixture of metallic and molecular hydrogen, the following relation based on the Dulong-Petit form is assumed:

[TABLE]

For the hydrogen elements the temperature can be increased up to the metallization temperature $T_{m}$. In this paper we present a simulation at one pressure condition ($P=140$GPa) corresponding to a transition temperature of $T_{m}=1500K$. Then, we assume a first-order phase transition occurs, which can be modeled by changing the metallization fraction as follows:

[TABLE]

where $L$ is the latent heat per unit mass and $\rho L$ is the latent heat per unit volume. The metallization fraction can increase until it reaches $f_{i}=1$, after which any additional heat again leads to a temperature increase. If $\Delta q_{i}^{(n)}$ is negative, the temperature/metallization fraction decreases, so the process occurs in reverse order.

II.4 Values for the material parameters

The material parameters used are given in Table 2. These values correspond to estimates at 140 GPa. For the thermal conductivities of molecular hydrogen and amorphous alumina no reliable values were found. However, from studies of similar materials in the literature we could make an estimate of these values. For H2, an estimate is made using values measured for a pressurized gas of normal H2 Moroe2011 at pressures of up to $\sim$100 MPa. Extrapolation of these results leads to a thermal conductivity of the same order of magnitude as that of MH. The specific heat and the optical properties were obtained by extrapolating values measured at ambient conditions to high pressure, using the Vinet equation of state Vinet1987 . The complex index of refraction was extrapolated using the Drude model for metals, and the Clausius-Mossotti relation (in both cases assuming that the dielectric function $\epsilon$ satisfies $\epsilon=\tilde{N}^{2}$). For $H_{2}$, the index of refraction measured at high pressure Evans1998 was extrapolated to the desired wavelengths using the Lorentz model. The FWHM time duration of the laser pulse was chosen as $\tau$ = 280 ns. The laser pulse energy per unit area is varied in the range $47$ kJ/m{}^{2}\leq A\leq 54\kJ/m2. A thermal boundary resistance of $10^{-7}$ m2 K/W was added between amorphous alumina and diamond. The other thermal boundary resistances were assumed small enough to be neglected.

III Simulation results

III.1 Heating curves

In order to compare with the experimental results Zaghoo2016 , heating curves were simulated. These show the temperature in the sample as a function of the input pulse energy. Experimentally, the peak temperature in a heating cycle is determined from a fit of the blackbody radiation, using a procedure discussed by Rekhi et al. Rekhi2003 . Numerically, this temperature is approximated by the maximum temperature at any point in time in any element $T_{max}$. Such a procedure is necessary since no steady state is reached. The heating curves are found by plotting $T_{max}$ as a function of the pulse energy per unit area. The result is shown in Fig. 3a); the heating curves measured by ZSS Zaghoo2016 are shown in Fig. 3b).

In Fig. 3a), the open circles connected by a (blue) dashed line show the results obtained with the parameters from Table 2. For pulse energies above $49.2$ kJ/m2 metallization is reached; the temperature scales linearly for lower pulse energy. When the pulse energy is high enough to reach metallization, the absorption strongly increases, leading to a quick growth of the layer of metallic hydrogen and an abrupt rise in temperature, disallowing a plateau. For higher pulse energies, the absorption mainly takes place at the boundary between the metallized layer and the molecular hydrogen as the tungsten heater is increasingly obscured by the MH layer. Additionally, the pulse energy goes into thickening of the layer rather than to further increase the temperature. The object of a simulation is to reproduce the experimental observation. The absence of a plateau in the heating curve does not agree with the experimental observation Zaghoo2016 .

In order to obtain heating curves with a plateau, it proved necessary to substantially increase the latent heat of transformation. This removes the rapid layer growth, and allows the tungsten layer to remain active as a heater. By increasing the latent heat by factors of 25, 50 and 100, the curves with asterisks, triangles and squares (Fig. 3a) were obtained, respectively. The sudden jump in temperature disappears, and is gradually replaced by a plateau if the latent heat is larger than $85000$ MJ/m3. As the latent heat is increased, the width of the plateau increases and the thickness of the MH layer decreases. It should be noted that the MH layer remains thin ($\leq 15$ nm) for any point on the plateau, so that both the tungsten and the MH layer absorb the laser energy and act as heating sources in this case. The resulting slope of the heating curve after the plateau is significantly smaller than the slope before the plateau. At the onset of the transition (at pulse energies of $49.2$ kJ/m2), the maximum temperature in the DAC is located at the tungsten and is slightly higher than the transition temperature $T_{m}=1500$ K, as the thin (5 nm) alumina layer on the tungsten acts as a small thermal impedance.

III.2 Reflectance and transmittance

In the experiment Zaghoo2016 , the presence of the plateau is taken as a sign of a first-order phase transition. On the other hand, several authors Celliers2018 ; Rillo2019 have suggested that the plateau corresponds to the onset of absorption rather than metallization. A recent Quantum Monte Carlo calculation Rillo2019 indicates that the absorption starts to increase near the experimental plateaus, whereas the increase in reflection occurs at higher temperatures/pressures. This highlights the difficulty of correctly identifying the transition pressure and temperature. Here, we will follow Ref. Zaghoo2016 and interpret the onset of the plateau as indicative of the phase transition.

To determine whether the phase that results after the transition contains free electrons, ZSS measured the transmittance and reflectance and compared the results to a Drude free-electron model. It has also been pointed out that the free-electron model is not valid at the transition, but applies beyond the phase line when the electron density is higher Zaghoo2017 ; Morales2010 . In the plateau region close to the transition, molecules dissociate and reform on a short time scale, so the charge density exists but is small Morales2010 . Experimentally Zaghoo2016 the reflectance change is very small in this region, whereas the absorbance (transmittance change) is larger and grows with thickening of the MH. Density functional theory Morales2010 predicts a reasonable but small optical conductivity, yet a small but finite DC conductivity in the plateau region. The DC electrical conductivity seems to be a good order parameter and should have a discontinuity for this first order phase transition. The onset of the plateau correlates with the electron density and seems to be a proper choice for the transition line, not the abrupt rise of the reflectance that starts deeper in the plateau region.

Using expression (10) the transmittance and reflectance can also be calculated in our simulation. The minimal transmittance and maximal reflectance as a function of time are plotted as a function of input pulse energy, shown in Fig. 4 for $\lambda=500$ nm and $\rho L=170000$ MJ/m3. Both the transmittance and reflectance remain constant before the transition starts, since there is no change in optical parameters in our simulation. When the transition starts, the transmittance and reflectance change due to the layer of MH that is forming. As the pulse energy is increased, the MH layer grows thicker and thicker. The transmittance approaches zero and the reflectance approaches its bulk value. Near the final data points of Fig. 4, which correspond to a MH layer thickness of $\sim 20$ nm, the transmittance and reflectance are almost converged to the bulk value.

IV Discussion

In this paper we have used the proper geometry and thermodynamic properties for the experimental demonstration of LMH Zaghoo2016 ; Zaghoo2017 , whereas earlier simulations are not relevant to this experiment Geballe2012 ; Montoya2012 ; Goncharov2017 . We have found that the results of the simulation are quite sensitive to material properties such as thermal conductivity (TC) and latent heat, but also the optical properties of all materials in the DAC, as these play a large role in determining the absorbed laser power. We explored a large range of parameters. Of particular importance is the thermal conductivity of various materials and we have attempted to consider more realistic values. Thus we have added thermal boundary impedance, as well as considering that the highly stressed part of the diamond (culet region) will have a reduced conductivity due to deformation and defects, etc. Another interesting aspect was the thermal conductivity of the liquid molecular hydrogen that originally was thought to be much smaller than that of the LMH. To estimate this we used the results of Moroe et al. Moroe2011 who studied the TC of a pressurized gas of normal H2 (75% ortho, 25% para). Since we expect the TC of a high-density gas to be similar to that of a liquid we examined various extrapolations to higher T and P, and surprisingly found values similar to LMH, and thus chose similar values for the FEA.

In the FEA the observation of plateaus in heating curves was a challenge and we were only able to create these using very large values of the latent heat of dissociation. Morales et al. Morales2010 have estimated the latent heat or enthalpy change/atom in their DFT and quantum Monte Carlo analysis of the PPT. At $T=1000$ K and a pressure of 140 GPa, they found $\rho L=1500-1900$ J/cm3. This corresponds to $\sim$0.04-0.05 eV/atom, whereas the dissociation energy of a free molecule is 4.47 eV/atom. To obtain plateaus in our simulations, the latent heat needs to be a factor of about 50 larger. This corresponds to a dissociation energy of around $2$ eV, close to the value of 1-1.5 eV reported in shock-compression experiments Holmes1995 for corresponding densities.

It is possible that the dissociation of molecular hydrogen in the liquid phase is a multistep process. Norman and Saitov Norman2017_1 ; Norman2017_2 ; Norman2018 have studied the pair correlation function using DFT theory. They propose that the transition involves a two step process: molecules are ionized to form H${}_{2}^{+}$; this is followed by the H${}_{2}^{+}$ interacting with H2 to form H${}_{3}^{+}$ + H. The enthalpy of formation has not been calculated in the high density liquid phase; however we note that for a free molecule the ionization energy is $\sim 15.4$ eV Liu2009 . We also note H${}_{3}^{+}$ has been observed as a stable species in ionized solid hydrogen Momose2001 . Subsequently, Pierleoni, Holzmann, and Ceperley Pierleoni2017 , who also find H${}_{3}^{+}$ in their theory, have shown that it is unstable. Thus, the formation of atomic metallic hydrogen by a multistep process may be involved in the transition.

Comparing our results from Fig. 3a) to the heating curves measured by ZSS from Fig. 3b), it can be seen that there is a difference in the shape of the heating curves. ZSS measure heating curves with similar slopes before and after the plateau, whereas in our simulations the slope after the plateau is smaller than that before the plateau. This depends on the material parameters: for example, decreasing $k_{\text{MH}}$ in the simulation should mean higher temperatures since the heat is more localized, while increasing $k_{\text{H}_{2}}$ should decrease the initial slope since the heat would diffuse out into the molecular liquid more readily. Variation of these parameters may lead to results comparable to experiment Zaghoo2016 .

In the simulation, the maximum temperature in the DAC is determined as a function of the laser power entering the DAC, whereas in the experimental results, the heating curves were reported as a scaled function of the laser power measured by a power meter. Moreover, in the experiment the laser is defocused, so that only a fraction is absorbed in the tungsten film and part shines on the gasket. Hence, an undetermined transfer function is needed for a quantitative comparison of the heating curves.

Finally, a comparison is made between our results and the results measured by ZSS for the transmittance and reflectance. In their paper ZSS used a multilayer analysis to isolate the reflectance and transmittance of hydrogen from that of the tungsten layer. Here, in Fig. 4 we show the absolute reflectance and transmittance of the system including the tungsten layer. There are two important experimental observations: in the plateau region the reflectance change is very small while the transmittance already has a measureable decrease, and as the LMH thickens the reflectance, $R$, saturates to a bulk value, which is about 0.5 as observed in the experiment. Our simulation supports this behavior. When saturated so that there is no transmittance, the absorption of the LMH is $1-R$, or about 0.5.

V Determination of the Phase Line P-T Points.

In a dynamic experiment on liquid metallic deuterium, with rising pressure and temperature, Knudson et al. Knudson2015 observed very weak transmission of visible light (they measure at only one wavelength in the visible) and then a steep rise in reflectance attributed to metallization. Their sample was microns thick, compared to the DAC sample of MH which was nanometers thick; this can explain the differences in observed transmission (absorption). Knudson et al. ascribed the absorption to conduction to valence band transitions as the bandgap closes down. Their experiment evidently was not sensitive to observation of a plateau, possibly due to a paucity of data points. They used the onset of the rapid rise of the reflectance as the signature of the pressure of the transition; the temperature could not be measured and was calculated. They determined a semi-empirical phase line for the PPT. An issue with this phase line was that it had a very large isotope shift when compared to the hydrogen PPT line and by extrapolation was in strong disagreement with earlier dynamic measurements on deuterium.

Recently, there have been two new experiments on the PPT of deuterium. Measurements at static pressures in a DAC were used to determine a phase line using two criteria: onset of plateaus and onset of reflectance Zaghoo2018 . The line using the onset of reflectance criterion was used to compare to dynamic measurements that also use this criterion. There was a large disagreement with that of Knudson et al., while there was no conflict with earlier measurements, within the uncertainties. Comparison to the hydrogen phase line using the same criteria provided an experimental determination of the isotope shift. Celliers et al. Celliers2018 made dynamic measurements on deuterium at the National Ignition Facility. They found a phase line that differed from that of Knudson et al. Knudson2015 , and was in much better agreement with earlier measurements, the phase line of ref. Zaghoo2018 , and theory Pierleoni2016 ; Mazzola2018 . Just as in Knudson2015 they presented a semi-empirical phase line, as temperatures were calculated, not measured. In their measurements they first observed absorption, then reflection, and used the latter as the criterion for the phase line. They too attributed absorption to bandgap transitions.

At this point in the study of the PPT there is now reasonable agreement between recently published experiments. There is a smaller disagreement for the isotope effect measured in static experiments Zaghoo2016 ; Zaghoo2017 ; Zaghoo2018 ; Ohta2015 and theory Pierleoni2016 . Yet there is disagreement as to the method of determining the P-T points of the transition: either the onset of the plateau or the abrupt rise in reflectance. Although the plateaus seem to be most appropriate from the consideration of latent heat and our FEA (see also Fig. 4), it would be useful to have an experimental determination. An ideal order parameter for metallization is the DC electrical conductivity. If electrical wires were inserted into the sample one could measure the conductivity and heating curves simultaneously to make the determination.

It is also possible to make a determination with only optical measurements. In dynamic experiments the pressure is ramped up in time. If in the dynamic experiments optical transmission was measured at two wavelengths (say blue and red) a determination could be made. In the case of metallization creating free electron absorption, one would see simultaneous absorption of both wavelengths. Alternately, in the case of bandgap closure one would first see absorption in the blue and then the red with increasing pressure. In DAC experiments the sample is isochoric (fixed density) and the pressure is essentially constant; the temperature is pulsed up for each pressure or density. In the case of creation of free electrons due to the PPT all wavelengths would be absorbed simultaneously at the transition. In the case of bandgap closure, at a pressure sufficient to close the bandgap into the visible, one would observe absorption first in the blue, and then as pressure was increased, absorption would occur in both the blue and the red. This absorption should be approximately temperature independent. In the DAC experiments visible/IR absorption was not observed until the pressure reached the plateau region and then all wavelengths were absorbed. Evidently, the transition takes place before the bandgap is in the visible region. We conclude that the onset of the plateaus is the proper criterion for the transition.

VI Conclusion

We have presented a complete FEA of the liquid-liquid phase transition to LMH, also known as the PPT, using a simulation program that we have developed to include the many subtle conditions in the experimental study at static pressures and high temperatures in a DAC. We calculated heating curves for a number of conditions and were able to find plateaus. To simulate reasonable plateaus comparable to those observed in experiment it was necessary to use a substantially larger value of the latent heat than has been calculated for this many-body system. We discussed a more complicated multi-step process that may possibly explain this discrepancy and may lead to new, as of yet unexplored physics.

Optical transmittance and reflectance were also simulated, and are in reasonable agreement with the experimental results. We found that the onset of the plateau is a good indicator for the insulator-to-metal transition. We discussed experiments to determine whether a change in optical parameters is due to bandgap closure or a transition to a metallic state. Possible further steps in this investigation could be to expand the simulation to two/three dimensions, to include a more accurate optical response in the transition regime, and to include better estimates of the thermal (boundary) conductivities of the various materials. The current simulation scheme, which makes use of the Fresnel equations for the pulse absorption, will be useful to analyse other high-pressure experiments with laser-pulse heated DACs.

Acknowledgements.

We would like to thank Mohamed Zaghoo, Rachel Husband, Kaan Yay and Ori Noked for suggestions, useful advice and discussions of the results. M.H. acknowledges support from the University Research Fund (BOF) of the Antwerp University. I.F.S. thanks the NSF, grant DMR-1308641 and the DoE Stockpile Stewardship Academic Alliance Program, grant DE-NA0003346 for support of this research.