Structural dynamics of laser-irradiated gold nanofilms

Szymon L. Daraszewicz,1 Yvelin Giret,1,2 Nobuyasu Naruse,2 Yoshie Murooka,2 Jinfeng Yang,2 Dorothy M. Duffy,1 Alexander L. Shluger,1 and Katsumi Tanimura2 1Department of Physics and Astronomy, University College London, Gower Street, WC1E 6BT London, United Kingdom 2The Institute of Scientific and Industrial Research (ISIR), Osaka University, Mihogaoka 8-1, Ibaraki, Osaka 567-0047, Japan (Received 20 March 2013; revised manuscript received 30 August 2013; published 5 November 2013)


I. INTRODUCTION
2][3][4][5][6] In this context, metals are often used as a playground to study the dynamics of an excited electron gas.The fast and complex nuclear dynamics following electronic excitation can now be probed more directly by ultrafast electron diffraction (UED) and time-resolved x-ray diffraction (tr-XD). 5,7][10][11][12][13][14][15][16][17][18][19][20][21] Different regimes of excitation have been identified, depending on the energy deposited by the laser pulse in a gold nanofilm: a low-fluence regime in which sample does not melt, where coherent acoustic phonon generation is observed; 13 a medium-fluence regime in which the sample undergoes melting, where a competition between homogeneous and heterogeneous melting is identified; 9 a high-fluence regime in which electronic effects are expected to affect the melting process. 10,11,22The initial stages of the dynamics are extremely difficult to unravel and recent results still stir a lot of controversy 10,11,14,15 even for relatively simple gold nanofilms.In particular, it remains unclear whether strong photoexcitation leads to bond softening or hardening and results in thermal or nonthermal melting and at which absorbed fluences these processes should occur.In this paper, we attempt to shed more light on these issues by developing a theoretical model that allows us to make a direct quantitative comparison to UED signals for a large range of fluences.
The UED experiments provide global information on the evolution of the crystalline order inside the sample but do not deliver detailed atomistic resolution directly.Complete description of the photoinduced structural evolution requires an atomistic model capable of quantitatively reproducing the time evolution of diffraction intensities.A two-temperature (2T) model 23,24 is often used for describing the dynamics of electronically excited solids.The method exploits the notion that electronic excitation quickly (within several femtoseconds) creates a very high electronic temperature, while the nuclear subsystem remains cold.6][27] It is, therefore, referred to as the two-temperature MD (2T-MD) method.This method has already been applied to simulate the atomistic dynamics of gold nanofilms excited by different laser fluences. 9,28However, the results have not been directly compared with experimental data.
In this work, we measure the time evolution of Bragg peak intensities by relativistic UED and use the 2T-MD method to model the behavior of gold nanofilms under different absorbed fluences.Because the film thickness and the temporal scales of UED measurements are the same for both experiment and theory we can compare the Bragg peak evolutions directly.At moderate laser fluences our results are in excellent quantitative agreement with experiment 22 and with the theoretical results of Refs. 9 and 29.However, at high fluences, the time evolution of Bragg peaks calculated using the conventional 2T-MD method disagrees with experiment.We demonstrate that a much better agreement with UED data can be obtained by using an interatomic potential which directly depends on the electron temperature.Our calculations show that the volume of the gold lattice depends strongly on the electron temperature and therefore the laser heating in a freestanding thin-film setup cannot be treated as an isochoric process.The quantitative agreement between the temporal evolutions of the experimental and theoretical Bragg peaks at all fluences suggests that the 2T-MD method provides a faithful atomistic representation of the structural evolution of photoexcited gold films.
Below we provide the experimental data and details of calculations, emphasizing the differences with the previous work.Comparison with experiment requires that UED signals are obtained in the kinematic regime, i.e., where only single scattering events occur free from any multiple diffraction effects, 22 which would induce a decrease of the zero-order peak. 30We focus mainly on the detailed analysis of the structural information regarding the mechanism of photoinduced melting of the gold nanofilms, especially in the high-fluence regime.The results for low fluences will be published in more detail elsewhere. 22

II. EXPERIMENTAL SETUP
We used high-resolution UED experiments with relativistic 3.0-MeV electrons to probe the time evolution of diffraction peaks after the laser excitation.Ultrashort electron probe pulses were generated with a custom designed 1.6-cell S-band rf gun with a magnetic solenoid.Laser pulses generated by a Ti:Sapphire laser were time-synchronized with rf by adjusting the oscillator cavity length to phase-lock the laser output with the 79.3 MHz rf generated as the 36th subharmonic of the 2856 MHz accelerating rf.The copper photocathode was illuminated by the third harmonic of a fundamental wave (770 nm) of the laser output, with a 90-fs pulse full width at half maximum (FWHM), and the photoelectrons were accelerated in a high rf field to 3.0 MeV with a 10-Hz repetition rate, and collimated by a solenoid with a nearly Gaussian transverse shape (2 mm FWHM).The electron pulses were collimated to a 200-μm diameter by an aperture constructed from graphite before entering the diffraction chamber.
The 10 ± 2 or 35 ± 5-nm-thick single crystal gold films were placed on a gold mesh in the diffraction chamber.The samples were excited with 90-fs pulses of 3.1-eV photons.The diameter of the excitation beam was 800 μm, i.e. much larger than that of the probe electron beam.The incident angle of the pump-laser light was 14 • from the surface normal, and transmission electron diffraction was measured along the (001) direction of the specimens.To achieve high sensitivity to MeV electrons and a high damage threshold, a CsI(Tl) scintillator equipped with fiber optic plates (Hamamatsu photonics FOP11) was used to convert the diffraction pattern into an optical image with a spatial resolution of 50 μm.The optical image was then reflected at 45 • using a 5-μm-thick optical mirror onto a CCD camera.The system temporal resolution was determined to be 180 fs, including timing jitter between rf and fs laser.In order to detect diffraction patterns with fine line widths, a condenser lens (CL) of the diffraction chamber precisely collimates 200-μm-diameter beams on the sample, a diffraction lens (DL) provides a back-focal plane for expanded diffraction images, and a projection lens (PL) displays the diffraction patterns with desired fashion onto the detector.
The relativistic electron energy of the probe beam gives two crucial advantages over conventional UED systems. 11irstly, this allows us to minimize space-charge effects and hence to perform high-quality single shot measurements while maintaining the pulse width less than 100 fs.This feature is crucial for studying irreversible phase transformations, such as laser-induced solid/liquid transitions.Secondly, our diffractometer provides structural information almost free from any multiple diffraction and possible inelastic effects inducing transient (000)-order attenuation. 11,30At low probe electron energies, extinction distance is smaller than the sample thickness and hence quantitative interpretation of diffraction requires detailed analysis in terms of dynamical theory of electron diffraction.For 3.0-MeV electrons, the extinction distance for (200)-order in Au is 186 nm, much larger than the sample thickness ( 35 nm); hence multiple diffraction effects are negligible. 14,31In fact, for 3.0-MeV energy of probe electrons, the (000)-order peak intensity remains constant in our measurements, hence the kinematic theory assuming single scattering events can be applied.This is in contrast with the previous results obtained by conventional UED, where transient (000)-order attenuation, characteristic of multiple scattering processes, is unavoidable. 11,30

III. DETAILS OF CALCULATIONS
The 2T-MD method has been described in detail in Refs.9, 23, and 25-29, and therefore we describe it only briefly here.We calculate real-space atomistic correlations to obtain the structure factor 28,29 required for describing structural and thermal contributions in the temporal evolution of Bragg peaks.In order to characterize the effects of electronic excitations on the interatomic interaction, we carried out ab initio calculations of the phonon spectra at the free energy minimum volume for different values of electronic temperature and found the conditions where the electronic excitation significantly modifies the interatomic interactions.

A. Hybrid continuum-atomistic 2T-MD method
2T-MD solves the diffusion equation for the electronic temperature simultaneously with the modified MD equations of motion, which incorporate an electron-ion energy exchange term.Below, we briefly describe the method and focus on the computational setup used in our simulations.

Description of electronic subsystem
The rapid thermalization of electrons after absorption of the laser energy in thin gold films 14 (∼100 fs) allows us to assume a well-defined electronic temperature T e (z,t) throughout the sample, where z is the distance from the surface and t is the time.The laser spot diameter is typically much larger than the probed area, therefore lateral energy redistribution can be safely neglected.The 2T-MD model assumes that the T e evolution follows the heat diffusion equation: 23,24 C e (T e ) ∂T e ∂t = ∇ • (κ e ∇T e ) − G(T e ) (T e − T l ) + S(z,t), (1)   where C e (T e ), G(T e ), κ e , and T l are the electronic specific heat, the electron-ion coupling, the electronic heat conductivity, and the lattice temperature, respectively.T e -dependent parameters of Eq. ( 1) are obtained from ab initio calculations (see Sec. IIB).
The laser source term S(z,t) in Eq. ( 1) is described by a pulse with a Gaussian shape in time t and an exponential decreasing amplitude with respect to z: 19

S(z,t)
where F is the absorbed fluence, l p the optical penetration depth of the sample at the wavelength of the pump pulse, t p the duration of the pulse taken at the FWHM, and t 0 the time zero defined as the arrival of the maximum of the laser pulse on the sample surface.We note that one can easily consider different pulse shapes or multiple-pulse excitation through modification of Eq. ( 2).
As it has been experimentally confirmed that gold films up to 100 nm are homogeneously excited, 24 initial homogenous excitation is a fair assumption for our 10 and 35 nm films.Thus κ e and the z dependencies disappear from Eqs. ( 1) and (2).Finally, we neglect the blast force 16 resulting from the gradient of electronic temperature as we assume homogeneous excitations.

Description of ionic subsystem
Concurrently with the description of the electronic subsystem through Eqs. ( 1) and ( 2), the ion subsystem is modeled with classical MD, according to modified equations of motion: 26,27 where v i denotes velocity of atom i of mass m.Here, F i is the classical force on an atom resulting from the gradient of the interatomic potential, while Fi is the additional driving term, based on a modified Langevin thermostat formulation: where γ represents the frictional drag and f L is the stochastic force.In order to satisfy the fluctuation-dissipation theorem, it is assumed that the stochastic force has a Gaussian distribution: where k B is the Boltzmann's constant and T t the target temperature of the thermostat.
In order to represent a mechanism for energy transfer from electrons to ions, we consider an "out-of-equilibrium" Langevin thermostat, where T t in Eq. ( 6) is equal to the electronic temperature T e , giving an expression for the frictional drag term: 26,27,32,33 where V and N are the volume and number of atoms in a given coarse-grained temperature cell (voxel).The local lattice temperature is computed from the atomistic velocities: and takes into account the random (thermal) motion only, discounting for the velocity of the center of mass of the temperature voxel due to possible expansion of the film 29 (v i ): We simultaneously resolve Eqs. ( 1)-( 3), where γ is calculated at each time step from Eq. ( 7), T l from Eq. ( 8), and the stochastic force from: 26,27,32,33 where t is the MD time step and R a vector generated from three random numbers uniformly distributed between −1 and +1.
Our approach differs slightly from the one proposed in Ref. 9 in which the electron-ion coupling term that appears in Eq. ( 3) is expressed as an external force proportional to the thermal velocity of a particular atom and the electronphonon coupling strength.However, the formulation based on an inhomogeneous Langevin thermostat reflects the statistical nature of the energy-exchange process better and, in principle, allows for selective phonon excitations through modifications of the random force spectrum.Nevertheless, we are able to reproduce the results of Ref. 9 using the same parameters, but employing the modified Langevin thermostat, showing the similarities of these two formulations in the case of laserexcited metal targets.

Simulation setup
To calculate the classical forces F i , we have used a recent highly optimized embedded atom model (EAM) Au potential developed by Sheng et al., 34 which correctly reproduces the thermal and structural properties, such as the melting temperature and phonon spectrum.This potential was developed by fitting the potential energy surface derived from first-principles calculations and scaled to match the experimental reference data.Overall, it performs much better than the previously developed EAMs for Au, such as the commonly used ones developed by Johnson, 35 Foiles et al., 36 Lee et al., 37 and Gronchola et al. 38 Ensuring correct thermal parameters in the model is crucial for correctly describing the melting dynamics as a function of the energy delivered by the electron-phonon coupling.
We used an MD cell containing 250 000 atoms with a size of 20.4 × 20.4 × 10.2 nm to represent a 10 nm 001 -orientated Au film, and a cell containing 860k atoms with a size of 20.5 × 20.5 × 35.9 nm to represent a 35 nm 001 -orientated Au film.The overlaying T e and T l voxels are cubes with a side length of ∼1.4 nm.The MD cell boundary conditions are periodic in the lateral directions with two free (001) surfaces (see Fig. 1) to allow for uniaxial expansion.We have implemented this 2T-MD model in a local version of the DL-POLY (4.01) code. 39e used a constant time step of 1 fs for the MD part of the calculations and checked that the total energy (electrons and ions) is conserved.Since we consider a uniform excitation at each time step, the finite difference solver time step for Eqs. ( 1) and (2) (i.e., the 2T continuum part) was equal to the MD time step.The atomistic configurations were pre-equilibrated in a constant pressure and temperature (NPT) ensemble (1 atm, 300 K) and subsequently in a constant volume and temperature (NVT) ensemble at 300 K.
Finally, to distinguish solidlike from liquidlike atomistic surroundings, we used a nearest-neighbor averaged centrosymmetry parameter ( i ) for each atom i, computed The energy of the laser pulse is initially given to the electronic temperature grid points, which can exchange this energy (e-p coupling process) with the coarse-grained cells (voxels) of ionic temperature.The conductivity in the electronic system is assumed to be infinite and the initial electronic energy distribution is homogeneous (see text).When the MD system expands along the free surfaces, the ionic temperature voxels become activated once a sufficient number of atoms occupies them (after Ref. 9). according to 28,40 where d j and d −j are pairs of vectors that connect an atom i to the opposite nearest neighbors j and −j .The value of the centrosymmetry parameter is zero for atoms in perfect crystalline surroundings and sharply increases as the local atomistic environment becomes disordered.The values of the centrosymmetry parameter are dimensionless through normalization by the square of the lattice parameter.

Bragg peaks calculation
The atomistic information obtained from 2T-MD enables us to calculate the structure factor S(Q) from real-space correlations. 29The time evolution of Bragg peaks can be obtained from S(Q) through a one-dimensional sine Fourier transform of the pair density function: where r ij is the distance between atoms i and j , on the positive radial distance (r) axis: 29 where W (r) is a damping modification function to suppress artificial ripples resulting from a cutoff of r max applied in the evaluation of ρ(r).This approach offers significant advantages over Debye-Waller factor calculations, as it describes both structural and thermal contributions to the Bragg peak changes.Nonetheless, this scheme assumes a scattering equation for a random sample orientation (i.e., in the case of a crystalline powder). 41Therefore care should be taken when applying it to the analysis of diffraction patterns of monocrystalline films taken from a particular beam axis (also referred to as zone axis), as it would allow for forbidden reflections to appear in the computed pattern.Since the experimental electron beam axis was placed normal to the (001) film surface (i.e., on [001] zone axis) any peak splitting or shifts, which result from uniaxial expansion in the (001) directions, will not be detectable.For this reason, when computing the areas under the selected Bragg peaks any visible splitting in the computed S(Q) pattern is neglected.In particular, in the integration procedure we neglect the peak that splits of the main one, when it is clearly separated.Figure 2 shows an example of a computed temporal evolution of a pair density function and its corresponding S(Q) spectrum.

B. Ab initio calculations
The values of C e (T e ) and G(T e ) that appear in Eq. (1) were obtained through ab initio calculations performed using the ABINIT code. 42,43Calculations are based on the local density approximation (LDA) 44 together with the norm-conserving pseudopotential method, 45 where the 5d and 6s electrons are retained as valence electrons.The valence pseudo-wavefunctions are expanded into plane waves up to a 60-Ha cutoff and a 16×16×16 Monkhorst-Pack k-points mesh is used.In this framework, basic structural and electronic properties are in good agreement with literature values.We tried different pseudopotentials, different flavors of the exchange-correlation energy, with and without spin-orbit coupling, without noticing significant deviations from the LDA results.
Calculations at high electronic temperature implement the generalization of the density functional theory (DFT) 47,48 given by Mermin. 49The electronic specific heat is evaluated from the internal energy E e calculated for different electronic temperatures: C e (T e ) = ∂E e /∂T e .It has been shown recently that using the T e dependence of C e obtained from the free electron gas model (C e (T e ) ∝ T e ) can strongly overestimate or underestimate the electronic temperature for a given absorbed fluence. 19,50Moreover, by using results obtained from ab initio electronic energy for different T e instead of using ground state electronic density of states (DOS), 50 we implicitly take into account both the effect of the electronic DOS and its modification at elevated electronic temperatures. 10he effective e-p coupling function G(T e ) is also estimated from ab initio calculations where we used the experimental value at room temperature 24 G 0 = 2.1 × 10 16 W m −3 K −1 and calculated the dependence on electronic temperature using ab initio densities of states as proposed in Refs.50 and 51: where ε F is the Fermi energy, g(ε) the electronic DOS, and f the Fermi-Dirac distribution.It has recently been shown that this form of G(T e ) is appropriate for metals. 52The values obtained for both C e (T e ) and G(T e ) are similar to those by Lin et al. 50We neglected the effects of lattice temperature and structural changes on the electron-phonon coupling as they are unknown for gold. 53

Phonon spectra calculations
In order to characterize the effects of electronic excitations on the interatomic interaction, independently from the 2T-MD simulations, we carried out ab initio calculations of the phonon spectra at the free energy minimum volume for different values of T e .The free energy minimum volume is denoted as V eq (T e ).Lattice stability of a crystal can be calculated using the density functional perturbation theory (DFPT), 54,55 as implemented in the ABINIT code. 42,43Dynamical matrices were computed on a 8 × 8 × 8 q-points grid in the Brillouin zone and used for interpolation to obtain phonon spectra for different T e along the [100], [110], and [111] high-symmetry directions.The phonon spectrum calculated for V eq (300 K) is in good agreement with the experimental measurements, 46 as shown in Fig. 3(a).The phonon spectrum for high T e calculated for the 300 K equilibrium volume V eq (300 K), shown in Fig. 3(b), agrees well with the previous theoretical calculations 10 and demonstrates a strong hardening of the phonon branches when the unit cell volume is kept constant.Our calculations show that there is no hardening effect at constant volume for electronic temperatures lower than T e ∼ 25 000 K.
Recently, UED measurements on gold nanofilms have been interpreted in terms of electronic bond hardening, 11 in order to explain a slower decay of the Debye-Waller factor than expected from an isochoric continuum 2T model.This effect has been first discussed theoretically by Mazevet et al., 10 where the evolution of the phonon spectrum of gold as a function of T e calculated at the equilibrium volume V eq (300 K) demonstrated a strong hardening of the dispersion branches at elevated T e .This has been interpreted as a reduction of the screening by the excitation of the 5d electrons making the electron-ion potential more attractive.To investigate this effect further, we calculated the electron density of gold at equilibrium volume for different T e .Figure 4 shows the difference in the total electron density n(r) in the primitive cell between the excited (T e = 50 000 K) and the room-temperature gold (T e = 300 K).It demonstrates that, overall, there is a net migration  of electrons from the volume around the nuclei to the center of the unit cell.This redistribution of the electron density reduces the screening between the nuclei and results in an increased repulsive interaction.Our calculations also show that 5d electrons are slightly more localized around the nuclei at high T e and that the delocalization is dominated by 6s electrons.Importantly, the results agree with those of Mazevet et al., 10 as long as the volume of the unit cell does not change.
Figure 5 shows the ab initio calculations of the Au cell parameter [corresponding to the minimum of the free energy: a 0 (T e ) = 3 √ 4 V eq (T e )] as a function of T e (black circles).One can then expect that above T e ∼ 9000 K ( a 0 > 1%) electronic effects may have a significant impact on lattice dynamics, as the system will expand in response to the extra force coming from the shift in the interatomic potential minimum.Furthermore, the phonon spectra for T e = 300 and 10 000 K presented in Fig. 3(c) demonstrate that, if one is using V eq 0 (300 K) to calculate the phonon spectrum at T e = 10 000 K, almost no modifications occur compared to the room temperature phonon spectrum.However, when using V eq (10 000 K), we observe a noticeable softening of all branches.Crucially, this softening trend persists for all electronic temperatures when calculations are made at V eq (T e ).

T e -dependent potential
The significant shift in the interatomic potential minima at T e 9000 K predicted by our ab initio calculations ( a 0 >1%) implies that in the high absorbed fluence regime, the use of a potential that does not capture the modified interactions resulting from the redistribution of the electron density is no longer a good approximation.We have investigated the accuracy of the only available electronic temperature dependent (ETD) potential for Au recently developed by Norman et al., 56 which was parameterized with respect to several T e points by the force-matching technique. 57We interpolated it using cubic splines in 0.05-eV increments and examined whether it reproduces the increase in the ab initio equilibrium lattice spacing in the electronic temperature range considered here.
Figure 5 shows that the mismatch between the ETD and ab initio lattice parameter increases as a function of T e .We obtain a linear evolution of the lattice parameter from our interpolated ETD potentials (blue diamonds), which can be attributed to the very limited number of T e points in the ETD potential (only k B T e = 0.01, 3, and 6 eV have been calculated so far) and possibly to the inherent limitations of the force-matching technique.To ensure that the interatomic potential used in the model accurately reproduces the ab initio lattice parameter at different electronic temperatures in the range considered here (<1 eV), we matched the T e values obtained from the 2T-MD simulations with the ab initio lattice parameters, and subsequently map these onto the interpolated ETD potential.For the F = 25 mJ cm −2 simulation (described in Sec.III), we choose an ETD potential at a particular increment corresponding to the T e averaged over the first 3 ps of the simulation.For simplicity, the ETD potential is kept at constant T e during the simulation required to reproduce the Bragg peak intensities, which is run only for 8 ps.Such setup is a reasonable representation of physical processes as during that time period T e does not change appreciably [see Fig. 7(a)].When analyzing the long-term behavior (> 8 ps) the ETD potential was changed in 0.05 eV increments.

A. Bragg peak evolution
The 2T-MD method described above is free from adjustable parameters and allows one to calculate Bragg-peak intensities from atomistic dynamics for the given energy absorbed by the Au film.Hence, to compare directly theory and experiment, it is important to have an accurate estimate of the adsorbed fluence, F .The absorbed fluence is determined from the incident fluence (F inc ) as well as reflectivity (R), transmission (T ) coefficients of the film and a coefficient dependent of the experimental setup (η exp ): 22 where L is the film thickness.The coefficient η exp describes losses due to the energy dissipation into the supporting grid (via ballistic electrons 16 ) and/or electron ejection. 12Our independent time-resolved optical absorption measurements 22 gave the value of η exp = 0.5.In this work, we present the results for three incident fluences F inc = 27, 41 (10-nm films) and 108 mJ cm −2 (35-nm film), which we thereby relate to the corresponding absorbed fluences F = 3.0, 4.5, and 25.0 mJ cm −2 .Some of the results for low fluence excitation of the 10-nm gold film are discussed in Ref. 22.We present a comparison between the time evolution of the experimental and theoretical (200) Bragg peaks at F = 3.0 and 4.5 mJ cm −2 as well as for the high fluence 25.0 mJ cm −2 in Fig. 6(a).For this fluence, we show the measurements for the (200) and (220) peaks, where the dashed and solid lines represent the results obtained with the ground-state and ETD EAM potentials, respectively.The agreement is excellent for all considered fluences within the whole measured time domains with no fitting parameters used.This demonstrates that the 2T-MD model captures the essential physics of the film behavior and that the time-resolved optical absorption measurements 22 provided accurate values of the absorbed fluence.We now discuss the structural changes observed in our 2T-MD simulations in more detail.

B. Atomic structure evolution
Depending on the fluence, we observed three different types of melting dynamics: a slow heterogeneous melting (F < F m ), a rapid homogeneous melting (F m < F < F e ), and an ultrafast nonthermal/thermal melting (F > F e ), where F m denotes the threshold fluence of thermal melting and F e a threshold above which a significant shift in the interatomic potential minima occurs.The values of these thresholds (F m = 3.1 mJ cm −2 and F e = 5.0 mJ cm −2 for a 10 nm film) are calculated in Appendix.Figure 7 shows the time evolution of the averaged temperatures, self-diffusivity, and atomic density in each of these three fluence regimes, which highlight the differences between the types of the melting processes, which we discuss in detail below.

Low-fluence regime
At low fluence (F < F m ), we observe premelting of the free surfaces and heterogeneous thermal melting by melt front propagation (see Fig. 8).The averaged lattice temperature (T l ) reaches the melting temperature already at 12.5 ps, however, the middle part of the film remains crystalline, while the melt fronts slowly propagate towards the center.The local temperature never exceeds the limit of crystal stability 9 (1.25 T m ) and therefore the solid remains superheated until it is overrun by the propagation of a relatively cooler melt front.The melting is accompanied by the temperature and density decrease and the speed of the melt front propagation decreases after ∼150 ps never exceeding a few percent of Heterogeneous melting leaves a superheated sheet of solid in the middle of the sample until the surface start to join up at 400 ps.The liquid at the surfaces is typically undercooled.At 1.5 ns a small crystalline pocket coexists with the liquid.The time frames are chosen in order to include the beginning and the end of the melting process.
the sound velocity.The same order of magnitude has recently been reported for the velocity of the melt front propagation in photoexcited Ni films. 25uring the cooling process, the sample thickness, and the ionic temperature in the middle of the film oscillate with a frequency of 2L/v c , where v c is the sound velocity at equilibrium conditions.This is in an agreement with very recent experimental observations of coherent acoustic phonon generation in a laser-excited thin gold films. 13We observe that the two melt fronts join up by a thin filament at 400 ps, which subsequently grows and at 1.5 ns the film melts almost entirely with only a small pocket of crystalline gold remaining.The almost complete melting of the sample is consistent with the fact that the nanofilm cannot be found on the UED sample holder mesh after the measurements are taken.

Medium-fluence regime
At medium fluence (F m < F < F e ), the sample expands more rapidly and the pre-melting is more pronounced [see Fig. 7(c)].This is accompanied, as in the previous case, by oscillations in the temperature and film thickness.The average temperature reaches T m at 6 ps and subsequently (6-12 ps) homogenously distributed seeds of low-density molten phase are created and destroyed in an superheated state (T m < T l < 1.25 T m ).When the sample reaches the limit of crystal stability locally (after 12 ps), these molten seeds, which served as nucleation sites, grow and coalesce until the sample melts entirely at around 20 ps (see Fig. 9).Notably, the average sample temperature does not exceed the crystal stability limit and therefore the homogeneous melting proceeds at a slower rate than in Ref. 9, where a collapse of the entire lattice is observed within 3 ps.Melting by solid-liquid interface The sample is superheated at 6 ps and the melting has already started from the surfaces.However, the middle of the sample remains largely crystalline until it locally approaches the crystal stability temperature (after 12 ps) when the homogenously distributed molten sites rapidly grow and coalesce completing the melting process at 20 ps.The time frames are chosen in order to include the beginning and the end of the melting process.
propagation is also observed, but it is less significant than in the heterogeneous case.As in the low fluence scenario, melting is accompanied by a decrease of density and temperature.Such complex melting dynamics at intermediate fluence is consistent with recent theoretical predictions, 28,29 however, alternative theories exist. 58urthermore, inside the molten pockets and when the sample is entirely melted, the calculated self-diffusivity is the same as equilibrium liquid gold at equivalent conditions (see Fig. 7).This shows that the laser-induced disordered state of thin gold films corresponds to an equilibrium liquid gold, in contrast with some semiconductors, where a highly coordinated liquid has been speculated. 59

High-fluence regime
The agreement of the 2T-MD model with UED experiment suggests that the character of interaction between gold atoms remains relatively unchanged at the fluences described so far.However, at higher fluence (F > F e ), the effect of the electronic excitations on the interatomic interactions becomes noticeable, as indicated by the results from ab initio calculations shown in Figs. 2, 3, and 5.In particular, the rapid drop of the Bragg peak intensities (to ∼ 50% in 1 ps) cannot be explained by a thermal model which assumes that the e-p coupling is the only energy transfer channel (see Fig. 6).
To include the effect of modified interactions at high T e , we employ the ETD potential that takes into account the reduced nuclear screening caused by the excited electronic distribution 56 (see Sec. IIC).The reduced screening results in a net repulsive force between the atoms, which causes the surfaces to rapidly expand and melt (see Fig. 10).The core of the sample becomes compressed due to the build-up of mean pressure (calculated from the stresses in the nanofilm) from the modified forces.The mean global pressure change, which results from the modified interatomic interactions, is 3 GPa, and the mean pressure reaches a maximum of 5.5 GPa at 3 ps due to the e-p coupling process.For comparison, a maximum of 1.5 GPa is reached in the F = 4.5 mJ cm −2 case.
Subsequently, after ∼2 ps, the remaining crystalline parts of the sample melt by the growth of nucleation sites concentrated around the melt fronts propagating from the surfaces.The expansion is very violent in this regime, and the fronts of decreased density coming from the surfaces precede the melt fronts propagation, as can be seen in Fig. 10.After melting, other fronts of decreased density propagate, and the film continues to expand until voids are created in the middle (see Fig. 11), which coalesce almost breaking the film into two parts at later stages (100 ps).As in the previous regimes, we did not see appreciable differences in self-diffusion coefficient and pair density functions between the molten film and equilibrium liquid gold.However, due to the change in the equilibrium lattice spacing at high T e , the density of the nonequilibrium liquid state of matter at 9 ps (1800 K, 2 GPa) of 0.051 atoms/ Å3 is approximately that of an equilibrium (i.e., T e = T l ) liquid at 2000 K and 2 GPa.The corresponding diffusion coefficients of 3.5 × 10 −9 m 2 /s (calculated with ETD potential) and 4.1 × 10 −9 m 2 /s (ground-state potential) are also very similar.
Melting at this fluence involves both nonthermal (T edependent interatomic forces) and thermal (e-p coupling) effects, with the former process dominating within the first ps.Therefore the decay of Bragg peak intensity at high laser fluence is a direct manifestation of the evolution of the energy landscape experienced by the atoms under photoexcitation.In this regime, we cannot qualify the melting as solely nonthermal, as it occurs above the melting temperature, and hence perhaps this process should be referred to as "nonthermally accelerated melting."

V. CONCLUSIONS
We have employed ultrafast electron diffraction measurements and a 2T-MD model to obtain detailed atomistic information on the structural dynamics of thin Au films following laser irradiation at a range of fluences.The model quantitatively reproduces the time evolution of the Bragg peak intensities measured by UED without any fitting parameters, and therefore it provides with a faithful representation of the real atomistic dynamics on a picosecond time scale.
At moderate fluences we observed heterogeneous melting, where the melt front propagates from the free surfaces and homogeneous melting, where molten pockets are formed and grow in the center of the film in agreement with previous calculations. 9,28In these two cases, density and temperature oscillations were identified.
In the high-fluence regime, the high electronic temperatures reached took the simulations into a regime where the groundstate interatomic potential did not provide a good description of the interatomic interactions.In this regime, we employed an interatomic potential that accounts for the expansion induced by the reduced screening due to the spatial redistribution of the conduction electrons.This nonthermal mechanism resulted in rapid expansion, and rapid melting, of the film and no density oscillations were observed.
Our ab initio calculations of phonon spectra suggest bond softening, if gold samples are allowed to expand freely under electronic pressure, and bond hardening, if they are constrained in all three dimensions.This is in contradiction to earlier theoretical research, which imposed isochoric constraints and, consequently, identified bond hardening in Au thin films under laser irradiation.The results of our work, however, do not provide unambiguous evidence regarding bond softening or hardening as the sample was constrained in two dimensions in UED experiments.This issue therefore requires further investigation.
We show that the volume of the gold lattice depends strongly on T e and therefore the laser heating in a free standing thin film setup cannot be treated as an isochoric process (see also Ref. 8 for discussion).Our conclusions regarding the shift of the equilibrium lattice spacing leading to rapid surface expansion can be verified experimentally by measuring the speed of surfaces expansion at a picosecond resolution.Such experiments, based on the Fourier Domain Interferometry technique, are currently available. 15The ab initio parameterized 2T-MD methodology with a T e -dependent potential described here can be readily applied to describe and understand photoinduced phase transitions in different metal films.
To conclude, the results demonstrate that the 2T-MD method complements UED with detailed atomistic dynamics and provides unprecedented insight into the melting behavior of metals following laser irradiation.

APPENDIX: THRESHOLD FLUENCE CALCULATIONS
We estimated the thermal melting fluence threshold from thermodynamic considerations by calculating the deposited energy required to increase the ionic and electronic temperatures up to the melting point (T m ) and to overcome the enthalpy of melting ( H m ): where C l and C e are the electronic and lattice specific heats, respectively, L the depth of the film, and T 0 = 300 K. From the characteristics of the EAM potential employed, we estimate that above the threshold fluence of F m ∼ 3.1 mJ cm −2 (10-nm film) enough energy is provided to cause complete melting of a film.
Based on the calculated phonon spectra presented in Sec.IIC, we have estimated an electronic temperature threshold (T th e ∼ 25 000 K) for lattice hardening at constant volume in gold.To make a direct link with experiments, we converted this threshold to an absorbed fluence.We calculated the maximum T e as a function of F by combining Eq. ( 1) (with G(T e ) = 0 and ∇T e = 0) with Eq. ( 2  where the duration t p of the laser pulse, unlike the thickness L of the sample, does not influence the reached T e .To solve Eq. (A2), we used a time step of 0.1 fs and a laser duration of 100 fs.
Figure 12(a) shows the maximum electronic temperature calculated from Eq. (A2) as a function of F for gold films of different thicknesses.We can see that a high absorbed fluence is required to reach the hardening regime for Au even for a 10-nm film (∼45 mJ cm −2 ).This result is in contradiction with the analysis of Ernstorfer et al. 11 where bond hardening for a 20-nm-thick gold film is claimed for an absorbed fluence of ∼47 mJ cm −2 .We were unable to reproduce their theoretical analysis based on the continuum 2T model using the same parameters, even when employing the same (unjustified) rescaling of the Debye temperature estimated from the phonon spectra calculations. 10Despite the high quality of the experimental study, 11 we do not believe that bond hardening has yet been observed in thin gold films.
Figure 12(b) shows the calculated threshold absorbed fluence as a function of the film thickness, where we found a perfect linear relationship with a slope of 4.37 × 10 7 mJ cm −3 .This result could be used in order to determine the range of fluences to use for a given film thickness in order to possibly measure electronic bond hardening in an experimental configuration that does not allow sample expansion.
Nevertheless, to define a threshold above which electronic effects on laser-induced gold become significant, we used the T e -dependent unit cell volume (shown in Fig. 5).Above T e ∼ 9000 K, the equilibrium lattice parameter changes by more than 1%, and the electronic effects on the interatomic potential become important.This corresponds to an absorbed fluence of F e ∼ 5 mJ cm −2 for a 10-nm film.

FIG. 1 .
FIG. 1. (Color online) Schematic representation of the 2T-MD simulation setup.The energy of the laser pulse is initially given to the electronic temperature grid points, which can exchange this energy (e-p coupling process) with the coarse-grained cells (voxels) of ionic temperature.The conductivity in the electronic system is assumed to be infinite and the initial electronic energy distribution is homogeneous (see text).When the MD system expands along the free surfaces, the ionic temperature voxels become activated once a sufficient number of atoms occupies them (after Ref. 9).

2 )FIG. 2 .
FIG. 2. Temporal evolution of the pair density function (left panel) and the corresponding structure factor (right) computed from the F = 4.5 mJ cm −2 (10 nm) 2T-MD simulation.Peak splitting is due to the uniaxial expansion of the sample, however, it is not detected in the experimental setup because of the beam axis direction.To compute the area under the peaks, a Lorenzian fit with linear background subtraction was used.

FIG. 4 .
FIG. 4. (Color online) Total electron density difference n(r) at equilibrium volume V eq 0 (300 K) between T e = 50 000 and 300 K. Blue, red, and white colors represent positive, negative, and zero difference, respectively.One can see the net migration of electrons from the volume around the nuclei to the center of the unit cell.

FIG. 5 .
FIG. 5. (Color online) Interpolated (see text) electronic temperature dependent (ETD) interatomic potential (blue) and ab initio (black) equilibrium lattice spacing increase as a function of electronic temperature.The maximum electronic temperatures reached in the three different simulations are represented by vertical dashed lines.

FIG. 7 .
FIG. 7. (Color online) Time evolution of global physical properties for F = 3.0, 4.5, and 25 mJ cm −2 : (a) electronic and lattice temperatures, where the light-blue solid line represents the melting point (T m ) of the ground-state EAM potential at 1281 K, and the light-blue dashed line represent the crystal stability limit of 1.25 T m .(b) Self-diffusion coefficient (D) computed from the mean-square displacements in the lateral directions for the whole sample (circles) and molted pockets (squares) as identified by the centrosymmetry parameter.The EAM potential gives D ∼ 1.8 × 10 −9 m 2 /s for an equilibrium bulk liquid gold at 1300 K (grey line).(c) Global density shows periodic oscillations in the case of 3.0 and 4.5 mJ cm −2 fluences, whilst at 25 mJ cm −2 the sample becomes continuously less dense.The figure inset shows the surface expansion velocities.

FIG. 8 .
FIG. 8. (Color online) A simulation cross-section showing centrosymmetry (top), density (middle) and lattice temperature (bottom) evolutions for F = 3.0 mJ cm −2 .The red atoms in the top panel correspond to nearest-neighbor averaged > 0.45 and blue atoms to 0.45 [see Eq. (10)].The local atomic densities and temperatures are averaged over neighboring atoms within 19 and 12 Å, respectively.Heterogeneous melting leaves a superheated sheet of solid in the middle of the sample until the surface start to join up at 400 ps.The liquid at the surfaces is typically undercooled.At 1.5 ns a small crystalline pocket coexists with the liquid.The time frames are chosen in order to include the beginning and the end of the melting process.

FIG. 9 .
FIG. 9. (Color online) Centrosymmetry (top), density (middle), and lattice temperature (bottom) evolutions for F = 4.5 mJ cm −2 .The sample is superheated at 6 ps and the melting has already started from the surfaces.However, the middle of the sample remains largely crystalline until it locally approaches the crystal stability temperature (after 12 ps) when the homogenously distributed molten sites rapidly grow and coalesce completing the melting process at 20 ps.The time frames are chosen in order to include the beginning and the end of the melting process.

FIG. 11
FIG. 11. (Color online) The F = 25 mJ cm −2 simulation at 15-100 ps with variable ETD potential, color-coded according to the local density.At longer times, the 35-nm sample breaks into two parts held by a thin liquid filament after formation, and subsequent growth and coalescence of voids.

FIG. 12
FIG. 12. (Color online) (a) Maximum electronic temperature as a function of absorbed fluence for different film thicknesses [see Eq. (A2)], the horizontal black dashed line represents the threshold above which hardening can occur at constant volume.(b) Absorbed fluence threshold required to reach bond hardening regime as a function of the film thickness. ):