Ultrafast laser excitation of a metal is able to bring the material into a state far from equilibrium.1,2 The preferential and rapid heating of one subsystem over the other leads to a system of highly coupled cold ions immersed in a partially degenerate electron sea.3 These transient states commonly occur during the formation of high energy density plasmas, including warm dense matter (WDM), with particular relevance to laser micromachining4–7 and inertial confinement fusion experiments.8 In the laboratory, these transient states serve as a testbed where quantum mechanical theories of electron–ion interactions, nuclear dynamics, and phase transitions can be validated.9–12

Historically, the response of the electron subsystem has been measured in optical pump–probe experiments.13–17 However, these model-dependent techniques provided only a surface measurement of the electron properties and limited information on the ionic response. More recently, x-ray scattering experiments have measured the bulk electron temperature by observing the electron plasmon feature.18 By contrast, the bulk ion temperature has only been inferred from the atomic structure, measured through ultrafast electron10 or x-ray diffraction.19,20 When the ionic system is still crystalline, the Debye–Waller factor (i.e., the reduction in the intensity of the Laue diffraction peaks) is a practical, albeit model-dependent, method.19,20 However, modeling the decay of the Laue diffraction peak, which is ultimately dependent on the root-mean-square (rms) deviation of the atoms from their lattice positions, depends on several physical parameters:•the energy density of the sample ε,•the electron–ion equilibration rate g_ei(T_e, T_i),•the electronic heat capacity C_e(T_e),•the ionic heat capacity C_i(T_i),•the Debye temperature T_D(T_e, T_i),

where the generally assumed dependences on electron temperature T_e and ion temperature T_i are shown. While the electronic and ionic heat capacities are well constrained for gold,21,22 the values found in the literature for the remaining three parameters exhibit large uncertainties. In particular, theoretical predictions of the equilibration rate vary by up to an order of magnitude.22–29

During the last decade, three different experiments have used ultrafast electron diffraction to measure the decay of the Laue diffraction peaks in laser-irradiated gold. They each attempted to elucidate the effects of nonequilibrium species on the interatomic potential. i.e., to measure the nonequilibrium Debye temperature and answer the long-standing question surrounding the existence of bond hardening/softening in warm dense gold.10–12,30 For comparable fluences, each experiment measured similar decays in the Laue diffraction peaks. However, they reached opposing conclusions; this was primarily due to different assumptions regarding the behavior of the electron–ion equilibration rate and the initial energy density deposited into the sample. For example, Daraszewicz and co-workers found that the assumed initial energy density E, calculated by taking into consideration purely the reflected and transmitted light, was reduced by a factor of η:11ε=ηE,(1)where ε is the corrected energy density. They obtained a value of η ∼ 0.5, measured through independent optical absorption methods at lower laser fluences.31 At higher fluences, they found it necessary to treat η as a free parameter. This extra energy-loss pathway, specific to the target and experimental geometry, was attributed to ballistic/fast electron transport away from the system—despite contemporary claims that the value of η was negligible.32

One method of modeling such experiments is classical molecular dynamics (MD). When coupled to an electronic system through an appropriate thermostat, such simulations can, in one dimension, simulate the full target thickness, thus capturing the necessary strongly coupled nature of the ions.33 However, recent efforts using the highly optimized and validated interatomic potential developed by Sheng et al.34,35 were unable to reproduce the experimental results for all energy densities in tandem.36

In contrast to previous simulation work, we treat both the electron–ion equilibration rate g_ei and the initial energy density ε as free parameters to provide a consistent description of melting in warm dense gold. Initially, we make use of a constant g_ei throughout each simulation, as well as exploring three different interatomic potentials. In all cases, we are able to obtain excellent agreement between the decays of Laue diffraction peaks obtained experimentally and those calculated from synthetic diffraction patterns. We find a strongly energy-dependent electron–ion equilibration rate, in addition to considerable differences between the assumed energy density E and the energy density in the simulations ε. In fact, for each interatomic potential that we have tested, we have found that the success of the model is contingent on allowing energy to escape from the target region. This additional energy loss, which we also characterize by η, scales strongly with laser intensity. Additionally, we have implemented several published electron-temperature-dependent (ETD) and ion-temperature-dependent g_ei models.21–23,25,26,29 For all forms of g_ei, consistency between the simulated and measured diffraction patterns is only achieved by reducing the absorbed energy density.

By introducing the free parameter η, we present a consistent description of the melting process in warm dense gold, resolving long-standing discrepancies between MD simulations and experiments. Our finding of a non-negligible η parameter is in good agreement with the work of Daraszewicz and co-workers; however, we additionally observe a strong dependence on laser intensity.

We perform large-scale MD simulations in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package37 in the canonical ensemble. The ions, treated explicitly, are coupled to an electron subsystem through a Langevin thermostat within the two-temperature model (TTM).38,39 Justified by the rapid electron thermalization time,40 we heat the electrons instantaneously to simulate the laser–matter interaction; we verify the accuracy of this procedure in the supplementary material by comparing the results of the forthcoming analysis with simulations taking the temporal width of the laser pulse into account. Synthetic diffraction patterns are produced for comparison with the experimental electron diffraction data of Mo et al.10

In this work, we focus on benchmarking the data obtained by Mo et al., which is expected to be below the threshold for considerable changes in interatomic bond strength.10,11 In keeping with the experimental conditions, we model free-standing 35 nm single-crystalline gold films and compare our Laue peak decay with measured results for the three assumed energy densities E=0.18, 0.36, and 1.17 MJ/kg, referred to here as the low, intermediate, and high energy density cases.

We model the nonequilibrium conditions of the laser-irradiated system via a traditional TTM that describes the evolution of the electron temperature throughCe(Te)∂Te∂t=−gei(Te−Ti),(2)where C_e(T_e) is the temperature-dependent electron heat capacity. All simulations utilize an electronic heat capacity derived from ab initio simulations21 and found to be in good agreement with previous thermal conductivity experiments.17 Within each simulation, the electron–ion equilibration rate is treated as a constant, and we employ a 1 fs time step throughout. In addition, before coupling the two subsystems, we equilibrate just the ions for 6.5 ps using a velocity-scaling thermostat and Berendsen barostat41 in order to bring the system to a temperature of 300 K at atmospheric pressure. This introduces an initial energy density into the ionic subsystem of ∼0.05 MJ/kg.

The ballistic electrons, accelerated by the optical laser, have a mean free path of ∼100 nm, approximately three times the thickness of the foil, and therefore we approximate the heating as isochoric.13,42–44 The high thermal conductivity45 also allows us to retain a single spatially invariant electron temperature and thus avoid complications due to lattice expansion.36 This is particularly important to avoid spurious edge effects, which may initiate heterogeneous melting. For the duration of the simulation, to model the electron–ion energy exchange, we couple the ions to this single-temperature electron bath via a Langevin thermostat given by46,47m∂vi∂t=Fi−γvi+fL(Tt),(3)where v_i is the velocity of ion i, F_i is the force acting on ion i due to its interaction with the surrounding atoms at time t, m is the ion mass, and γ is the friction parameter that characterizes the electron–ion equilibration rate. The friction parameter is related to the equilibration rate through γ = g_eim/3nk_B, where n is the ion number density and k_B is Boltzmann’s constant.38,39 The f_L(T_t) term is a stochastic force term with a Gaussian distribution and a mean and variance given by38,48〈fL(Tt)〉=0,(4)〈fL(Tt)⋅fL(Tt′)〉=2γkBTeδ(t−t′).(5)

We utilize the highly optimized embedded atom method (EAM) interatomic potential developed by Sheng et al.34 Under equilibrium conditions, the properties of this potential are closely matched with experimental values, as shown in Table I. Additionally, as demonstrated in Fig. 1, the atomic rms deviation and corresponding Debye temperature, obtained in x-ray diffraction experiments, are well matched at equilibrium temperatures up to melting. In nonequilibrium matter, the Sheng et al. potential has been able to successfully match experimental lattice disassembly times after laser irradiation.34,35

We utilize two simulation geometries: either a parallelepiped of (86a₀ × 15a₀ × 15a₀) or one of (86a₀ × 70a₀ × 70a₀), where a₀ = 4.078 Å is the lattice constant of gold predicted by the Sheng et al. potential. The smaller volume contains 77 400 atoms and the larger volume 1 685 600 atoms. In the larger direction, which corresponds to the 35 nm thickness of the gold foils used in the experiment, the simulation box is much larger than the size of the target geometry, creating front and rear surfaces capable of expansion. In the two smaller directions, we adopt periodic boundary conditions. The smaller geometry is employed when calculating the decay of the Laue peaks. The larger volume, which produces the same Laue decay curve as the smaller volume, allows for a higher resolution in reciprocal space and is used to create the synthetic spatially dependent diffraction patterns out to 50 ps.

We use the decay of the (220) and (420) Laue diffraction peaks as a metric to compare the simulation and experiment. We obtain the intensity of the diffraction peaks from the static structure factor of the system, which is easily calculated from the Fourier transform of the atomic positions.55,56 For 3.2 MeV electrons, whose wave vector (∼2 × 10¹³ rad/m) is orders of magnitude higher than their scattering vector (∼10¹⁰ rad/m), the Ewald sphere can be considered flat over a finite region of reciprocal space. As such, we take a Q_x–Q_y slice of reciprocal space, where the x and y coordinates run parallel to the surface of the foil. The scattering intensity is calculated by integrating over a specific reflection.

In finding the best fits for all three experimental datasets, we have performed hundreds of simulations utilizing the small simulation geometry. We have considered initial electron temperatures Te0 that, when converted into energy densities via the electron heat capacity,21 correspond to a range in η of 0 ≤ η ≤ 1, with a step size ΔTe0=100 K. We have considered values of g_ei in the range of 0 ≤ g_ei ≤ 20 × 10¹⁶ W m⁻³ K⁻¹, with a step size of approximately Δg_ei ∼ 0.3 × 10¹⁶ W m⁻³ K⁻¹. We quantify the best fit by identifying which (Te0,gei) pair minimizes the rms difference between the experimental and calculated decays in the intensity of the (220) and (420) Laue diffraction peaks. For each of the three energy densities investigated, this difference is shown in Figs. 2(a)–2(c), with the corresponding best decay curve and temperature evolution for the subsystems given in Figs. 2(d)–2(i). The black dashed line in Figs. 2(a)–2(c) encloses the region in which the simulation results lie within the experimental error bars. In each case, we are able to find a (Te0,gei) pair that gives excellent agreement with the experimentally measured decay curves.

For the low energy density case, we find that values of Te0=7800±300 K and g_ei = 2.2 ± 0.5 × 10¹⁶ W m⁻³ K⁻¹ produce the best fit. The obtained value of g_ei is in good agreement with experimental measurements of g_ei in room temperature gold.13,57 The initial electron temperature corresponds to a corrected energy density in the electron subsystem of ε ∼ 0.17 ± 0.01 MJ/kg, corresponding to η ∼ 0.92. In the intermediate energy density case, we find that Te0=8500±300 K and g_ei = 5.0 ± 0.5 × 10¹⁶ W m⁻³ K⁻¹ produce the best fit, corresponding to a corrected initial energy density for the electron subsystem of ε ∼ 0.21 ± 0.01 MJ/kg, with η ∼ 0.57. Finally, in the high energy density case, we find that the best fit to experimental data comes from values of Te0=9900±300 K and g_ei = 15.0 ± 2.0 × 10¹⁶ W m⁻³ K⁻¹, corresponding to a corrected initial energy density for the electron subsystem of ε ∼ 0.3 ± 0.01 MJ/kg, and thus η ∼ 0.26. Our analysis demonstrates a strong dependence of the parameter η and the electron–ion equilibration rate on the assumed energy density, i.e., the energy density calculated by taking into consideration purely the reflected and transmitted light (cf. the dashed black lines in Fig. 3).

The bond strength in nonequilibrium gold is predicted to change when the electron temperature is significantly higher than the ion temperature. We check the effect of this on our conclusions by repeating the analysis with two additional interatomic potentials designed to take into account higher temperature electrons,49 both of which are validated against finite temperature Kohn–Sham density functional theory.58 The first, shown by the dashed blue lines in Fig. 3, is calculated for electrons at a temperature of 0.1 eV (∼1200 K) and found to be in close agreement with those produced by the interatomic potential developed by Sheng et al. (a comparison of the properties of these two interatomic potentials can be found in Table I). Here, we find there to be a ∼11% difference between the best-fit results produced by the potential developed by Sheng et al. and that developed by Norman et al.34,49 For the second potential, calculated for an electron temperature of 1.5 eV (∼17 400 K), we only consider results in the high energy density case. Here, we find η = 0.15, which is less than that in the case of the Sheng potential, corresponding to an even more substantial reduction in absorbed energy (cf. the square red points in Fig. 3). Figures S1 and S2 in the supplementary material depict the full results of this analysis.

We find, across multiple potentials, an η parameter between 0.2 and 1 that scales strongly with the assumed energy density. Our result agrees qualitatively with previous work, which measured a value of 0.5 for similar conditions.11 In that work, it was suggested that the origin of the lost energy is energy dissipation into the supporting grid by ballistic electrons or electron ejection from the rear of the target.11 Other studies have suggested that the temporary formation of an electron sheath around the surface of the target59 accounts for a decrease in the energy density on a timescale commensurate with the decay of the Laue peaks.60,61 Our scaling with the assumed energy density, which is proportional to the laser fluence for constant target thickness, corroborates this conclusion. We note that this interpretation makes the η parameter specific to the target and experimental geometry and thus complicates the interpretation of such experiments.

For each of the three energy density cases, we take the best simulation achieved with a constant g_ei, and predict spatially resolved synthetic diffraction patterns. For the first 50 ps, these are created by Fourier transforming the atomic positions taken from simulations using the larger geometry and taking a Q_x–Q_y slice of reciprocal space. In each case, we confirm that these larger simulations reproduce the Laue peak decay seen in the smaller simulations. The diffraction patterns, shown in Figs. 4–6, correspond to the three energy density cases. The timescales of the experimental diffraction pattern measurements necessitate the use of the smaller simulation geometry in producing Figs. 4–6 for times above 50 ps. Also shown, for comparison, are the experimentally obtained diffraction patterns.

For a direct comparison with the measured spatially resolved electron scattering data of Mo et al., the simulated structure factor is multiplied by the electron scattering form factor, which is calculated using the Mott–Bethe formula62 utilizing tabulated x-ray elastic scattering cross-sections.63 Above ∼2 Å⁻¹, we find a negligible difference between the use of a cold or singularly ionized form factor. In addition, a background term of the form64,65dσineldΩ=AQ41−1(1+BQ2)2+mQ+c(6)is introduced to account for inelastic scattering. The coefficients A = 350 and B = 1.3 are kept constant for all comparisons, whereas small variations are required in the linear component between shots at different laser intensities. In each case, the background term is plotted as a dashed line in the figures and represents a small component of the overall signal. Finally, the entire image is convoluted with a pseudo-Voigt profile that represents the spatial profile of the electron beam66 (FWHM_Gauss = 17 meV and FWHM_Lorentz = 6.5 meV).

Initially, we find that the diffraction peaks from the solid gold are an order of magnitude larger than the liquid scattering signal, an effect not seen in the experimental data, where the difference is closer to a factor of two. We attribute this discrepancy to misalignment and beam divergence in the experiment that is not present in the simulations. Therefore, we average over a small ΔQ_z, reducing the intensity of the solid diffraction peaks while leaving the liquid diffraction peaks, which are broad in reciprocal space, unchanged. We choose ΔQ_z, which remains constant in our modeling, by matching the spatially resolved scattering signal at both early (t = 0 ps) and late (t = 17 ps) times in the high energy density case. With these changes, we are able to match, in each of the three energy density cases, the angularly resolved line-outs obtained in the experiment. It is important to note that the applied modifications to the raw data, used to obtain the angularly resolved data, affect only the qualitative comparisons provided in Figs. 4–6, and not the intensity decay curves in Fig. 2.

In our analysis of the melting behavior of the sample, we utilize the smaller simulation geometry. This is because of a lack of available computational resources that would allow us to perform simulations in the large simulation geometry that capture the full duration of the low and intermediate energy density cases (3000 and 800 ps, respectively). We are however, able to compare the small and large geometries, in all three energy density cases, out to 50 ps. In all energy density cases, we find the melting behavior and the expansion rates of the large and small geometry samples to be equivalent out to 50 ps.

From the time for the total disappearance of the (200) peak, Mo et al.10 concluded that the three energy density cases correspond to the incomplete, heterogeneous, and homogeneous melting processes, respectively. On this basis, they obtained a threshold for complete melting (assuming total absorption of the laser energy, η = 1) of ∼0.25 MJ/kg and a transition between heterogeneous and homogeneous melting of ∼0.38 MJ/kg. Of note, both of these values are higher than those predicted by molecular dynamics simulations employing the interatomic potential of Zhou et al.67,68

In our analysis, utilizing the interatomic potential developed by Sheng et al.,34 we are able to match the entire decay of the experimentally measured Laue peak. Using the OVITO visualization tool69 together with polyhedral template matching,70 we can directly assess the structure and melting process in each of our cases (cf. Fig. 7). In contrast to the original work, we find that the low energy density case exhibits heterogeneous melting, with a clear propagation of the melt front [see Figs. 7(a)–7(c)] and no nucleation of melt inside the target. The two higher energy densities both exhibit homogeneous melting, with the intermediate case appearing to be on the cusp of the two regimes, as determined from nonuniform melt regions inside the target. Taking this into consideration, and the lower predicted energy density from the η parameter, we find a total energy density threshold for complete melting below 0.22 MJ/kg (which includes the 0.17 MJ/kg present in the electron subsystem and the 0.05 MJ/kg already present in our 300 K gold sample). This value is in agreement with the expected value of 0.22 MJ/kg.71 For the transition between homogeneous and heterogeneous melting, we find a value that lies between 0.22 and 0.26 MJ/kg. With the introduction of the η parameter, we find that the complete melting threshold, and heterogeneous-to-homogeneous melting boundaries, are both lower than previously thought.10

Finally, we note that despite achieving excellent agreement with the experimental decay curves, we are not able to reproduce the low intensity, long lived, diffraction features present in the experiment [cf. Figs. 5(b) and 5(f)]. We find that the intermediate energy density case melts within 100 ps, whereas the original work concludes that it does not melt until 800 ps. To investigate if this is a fault of the equilibration rate, we manually adjust the value of g_ei throughout the simulation to best match the experimental results; we are unable to match both the early-time Laue peak decay and late time long lived diffraction patterns. In the simulation, this may highlight a deficiency in either the TTM model or the interatomic potential. In the experiment, this could be attributed to inhomogeneous heating.

The electron–ion equilibration rate g_ei is predicted to scale strongly with both electron and ion temperature. Thus, to supplement our study, we employ seven models of temperature-dependent equilibration rates.21–23,25,26,29 This set of simulations is performed under the same specifications outlined in Sec. II and use the smaller simulation geometry. As in Sec. III, we use the decay of the (220) and (420) Laue diffraction peaks to draw comparisons between simulation and the experimental results obtained by Mo et al.10 In each energy density case, we consider only a temperature range corresponding to η in the range ∼0 ≤ η ≤ 1, so as to consider only physically possible values. The performance of each model can be seen in Fig. 8. In the low energy density case, shown in Fig. 8(a), we find that the Medvedev and Milov23 model, calculated for a constant ion temperature of 300 K, and the Migdal et al.26 model perform best. For the intermediate energy density case, shown in Fig. 8(b), we find that the model developed by Smirnov25 has the lowest rms difference. Finally, in the high energy density case, shown in Fig. 8(c), we find that the models developed by Lin et al.,22 Holst et al.,21 and Petrov et al.29 perform best. Considering all three energy density cases simultaneously, the model developed by Smirnov25 performs best.

The success of these theoretical models is contingent on a similar inclusion of energy loss as was found in the constant-g_ei description. Using the best-fit value of Te0, we can calculate the η parameter for each theoretical model and each energy density case. The evolution of the η factor for both the constant-g_ei case and the temperature-dependent g_ei models is shown in Fig. 3(a). The electron-temperature-dependent models display a similar trend in η as the constant case, except for the solely electron-temperature-dependent model developed by Medvedev and Milov,23 calculated for a constant ion temperature of 300 K, and the model developed by Migdal et al.26 Excluding these outliers on the basis that they predict either energy gain or no energy loss in the system, we calculate the average η across all of our tested models. We find values of η̄∼0.83, 0.57, and 0.40 for the low, intermediate, and high energy density cases respectively, which are very close to that for the constant-g_ei case.

In Fig. 3(b), as a method of comparing models of g_ei with different functional dependences, we plot the average value of g_ei in the best-fit simulation for each theoretical model, as well as the constant-g_ei cases discussed in Sec. III. In each simulation, we average the value of g_ei out until the time at which the decrease in the Laue decay curves stops; in the low, intermediate, and high energy density cases, these times correspond to 50, 25, and 10 ps, respectively. In Figs. 3(a) and 3(b), the best-fit results for the constant and temperature-dependent g_ei models are in close agreement with one another for the low and intermediate energy density cases. However, they diverge in the high energy density case. This discrepancy could be due to bond softening, which is predicted to occur in gold, at its equilibrium volume, for electron temperatures T_e > 9000 K (or an energy density of ∼0.24 MJ/kg).11 Our results suggest that in the high energy density case, the initial electron temperature of the sample is expected to be above Te0∼9000 K.

We have benchmarked the molecular dynamics simulations of the ultrafast excitation of thin gold films against experimental MeV electron diffraction results published by Mo et al.10 We have carried out this analysis using a constant electron–ion coupling rate, seven prominent temperature-dependent electron–ion coupling rate models, and three different interatomic potentials. Using the Laue peak decay as a metric, and treating the initial energy density and the electron–ion equilibration rate as free parameters, we have found a strong dependence of the η parameter on initial energy density, with increased energy loss at higher energy densities (and higher laser fluence), which is in good agreement with the work of Daraszewicz et al.11 We believe the need for an additional energy loss pathway, and the presented scaling with increased fluence, to be reliably established. Ultimately, this ambiguity needs to be resolved via direct measurement of the subsystem temperatures.

Acknowledgment. The authors acknowledge M. Z. Mo for providing the experimental data that were used throughout this work. This work was funded in part by the U.S. Department of Energy, National Nuclear Security Administration (NNSA) (Award No. DE-NA0004039). This material is based upon work supported by the National Science Foundation under Grant No. 2045718. The authors acknowledge the support of Research & Innovation and the Office of Information Technology at the University of Nevada, Reno for computing time on the Pronghorn High Performance Computing Cluster. J.M.M. acknowledges funding and support from the Nevada NSF-EPSCoR NEXUS Award and the McNair Scholars Program at the University of Nevada, Reno.