Full-scale ab initio simulations of laser-driven atomistic dynamics

The coupling of excited states and ionic dynamics is the basic and challenging point for the materials response at extreme conditions. In laboratory, the intense laser produces transient nature and complexity with highly nonequilibrium states, making it extremely difficult and interesting for both experimental measurements and theoretical methods. With the inclusion of laser-excited states, we extended ab initio method into the direct simulations of whole laser-driven microscopic dynamics from solid to liquid. We constructed the framework of combining the electron-temperaturedependent deep neural network potential energy surface with hybrid atomistic-continuum approach, controlling non-adiabatic energy exchange and atomistic dynamics, which enables consistent interpretation of experimental data. By large scale ab inito simulations, we demonstrate that the nonthermal effects introduced by hot electrons play a dominant role in modulating the lattice dynamics, thermodynamic pathway, and structural transformation. We highlight that the present work provides a path to realistic computational studies of laser-driven processes, thus bridging the gap between experiments and simulations.


INTRODUCTION
Intense laser-matter interaction plays an important role in many applications including inertial confinement fusion [1], laser micromachining [2], and material synthesis [3].Ultrafast laser excitation can drive matter into extremely non-equilibrium states, in which the hot electron and cold lattice coexist.The subsequent atomistic dynamics is therefore a long-standing challenge, because it is governed by the interplay between excited-electronmodulated potential energy surface (PES) [4], electronion coupling [5], and geometric characteristics of irradiated samples [6].
Tremendous efforts based on time-resolved probing techniques and simulations have provided valuable insights into the nonthermal behaviors [7][8][9][10][11], kinetics of laser-driven melting [12][13][14][15], and electron-phonon coupling [16][17][18].The related processes from cold solid to hot liquid and plasma are the typical multiscale dynamics due to the cascade of interrelated processes triggered by the laser excitation, both in time scale and size scale.Therefore, it is of great difficulty and importance to construct a well-coordinated picture between experimental and theoretical efforts.For example, the dynamics of laser-excited Au is still under debates [4,7,19], regarding the phonon behaviors driven by laser excitation.In these cases, different priori assumptions on material response were usually made [20,21].
The above obscure stems from the technical limitations that the present methods can not capture both the nonthermal and intrinsic scale of laser-induced process at the same time.For ab initio methods such as time-dependent density functional theory (DFT), the sizes are limited to 10 1 ∼ 10 3 atoms and 10 2 ∼ 10 4 fs, unable to access real-istic representation of structural transformations of irradiated samples.While for classical molecular dynamics simulations coupled with two-temperature model (TTM-MD) [22,23], the implementation of empirical potential like embedded-atom-method (EAM) is limited in prior knowledge and model complexity, thus can hardly capture the high dimensional dependence of PES on both atomic local environments and electron occupations for a wide range of temperature and density [19,24], leading to the inadequate description of nonthermal nature of laser-driven processes.Therefore, bringing the advantage of ab initio and large-scale molecular dynamics including nonthermal effect becomes the route one must take.
In this paper, we developed ab initio atomisticcontinuum model by combining two-temperature-model (TTM) with an extended deep potential molecular dynamics (DPMD), as illustrated in Fig. 1.When the ultrafast laser interacts with solids, the electrons quickly thermalized at timescaes of femtoseconds, producing highly non-equilibrium states (electron temperature T e ≫ ion temperature T i ).The hot T e will result in the redistributed charge density firstly and then modifies the PES of ions.To capture this physics, we introduced the laser-excited PES by constructing electron-temperaturedependent deep neural network potential firstly, and coupled the PES into additional electron continuum subsystem via TTM-MD framework.By this way, we can directly simulate the whole electron-ion coupled dynamics during the laser-driven processes with large-scale simulations within ab initio accuracy.We take tungsten as an example and systematically validate the accuracy of our model in describing lattice dynamics, thermophysical properties, and laser heating process in both equilibrium and laser-excited states, by comparing with the related experimental results recently [13].(b) iterative concurrent learning scheme is used to efficiently sample atomic configurations for a wide range of equilibrium and non-equilibrium conditions.(c) hybrid atomistic continuum approach.The evolution of electron subsystem allows atomistic system transits between different PES, and the Langevin thermostat is introduced to mimic non-adiabatic energy exchange between electron and lattice.

RESULTS
Construction of laser-excited PES.Recent efforts have demonstrated the success of machine learning model towards large-scale simulations of ab initio quality at extreme conditions [25][26][27][28][29][30][31][32], but most of the studies focus on the equilibrium-state and ground-state applications.
Here the electron-temperature-dependent deep potential (ETD-DP) model is implemented in the framework of deep potential method [25,33,34] to model the laserdriven dynamics.
To avoid constructing hand-crafted features or kernels for different types of bulk systems, a general end-to-end symmetry preserving scheme is adopted [35].As illustrated in Fig. 1(a), the ETD-DP model consists of an embedding network and a fitting network.The embedding network is designed to transform the coordinate matrices R I to symmetry preserving features, encoded in the descriptor D I .And the fitting network is a standard fully connected feedforward neural network, mapping the descriptor to the atomic contribution of total energies.
The newly introduced parameter, electron temperature T e , is used to characterize the laser modulation on PES, in which electron occupation distribution is far away from electron-ion equilibrium states.This ETD-DP is defined as where A(R, T e ) is the potential energy depends on the local atomic environment (R) and T e , N αi denotes the neural network of specified chemical species of α i of atom i, and the descriptors D αi describes the symmetry preserved local environment of atom i with its neighbor list n(i) = {j|r ji < r cut }, respectively.
To generate an ETD-DP, the new degree of freedom, T e , will dramatically expand the sampling space in the data labelling process, introducing expensive computational costs.Therefore, an iterative concurrent learning scheme [36] is highly required to efficiently sample atomic configurations under both equilibrium (T e = T i ) and non-equilibrium conditions (T e ̸ = T i ).As shown in Fig. 1(b), to explore the density-temperature space with different electron occupations (ρ, T i , T e ), a variety of crystal structures are used as the initial configurations to run multiple DPMD simulations.And an ensemble of ETD-DP is trained with the same dataset but with different parameter initializations.The model deviation, denoted as the maximum standard deviation of the predicted atomic forces by the ensemble of ETD-DP, is used to evaluate whether the explored atomic configurations should be send to generate referenced ab initio energies, forces, and virial tensors.
Two-temperature model coupled DPMD (TTM-DPMD).To model the whole ultrafast laserdriven processes from cold solid to plasma, we should couple quantum electron subsystem and strongly coupled ionic subsystem.Here, we implemented our laser-excited PES into the TTM-MD framework [22,23,37], going beyond traditional ground-state EAM and neuralnetwork-driven PES descriptions.As shown in Fig. 1(c), the heat conduction equation of electron continuum characterizes the temporal evolution of electron occupations, thus governing the transition of ionic system between different T e -dependent PES.Langevin dynamics is incorporated to mimic the dynamic electron-ion collisions [23,[37][38][39].The TTM-DPMD is defined as follows, where C e is the electronic heat capacity, κ e the electronic thermal conductivity, g ei the electron-phonon coupling constant, S(r, t) the laser source.The ions evolves on the T e -dependent PES A(R, T e ), and suffers fluctuationdissipation forces −γ i v i + Fi (t) from electron sea.Here γ i is the friction parameter that characterizes the electronion equilibration rate, relating to the electron-phonon coupling constant through γ = g ei m i /2n i k B , where n i the ion number density.The Fi (t) term is a stochastic force term with a Gaussian distribution, whose mean and variance is given by ⟨ Fi (t)⟩ = 0 and ⟨ Fi (t) In TTM-DPMD, by practically choosing the electron temperature or ionic temperature in the meshgrid as the additional parameter in ETD-DP model, the ions can evolve under laser-excited PES (T e ≫ T i ) or ground-state PES (T e = T i ), so that we can separate the nonthermal effects defined by the electronic excitation from thermally driven atomic dynamics and phase transformation.
Validating neural network model for laserexcited tungsten.To validate the effectiveness of extended DP model, we chose tungsten as our target system.Tungsten is a typical transition metal, with half-filled d bands that is sensitive to T e .Upon laser excitation, tungsten is expected to go through a complicated dynamic process including possible nonthermal solid-solid phase transition [8,20,21,24], attracting much attention but remains ambiguous.Here we generate a T e -dependent deep-neural-network tungsten model by learning from DFT data calculated with the generalized gradient approximation (GGA) of the exchangecorrelation functional [40] using VASP package [41].The atomic configurations used in the training set are collected from a wide range of (ρ, T i , T e ) condition, covering the phase space of the body-centered-cubic (BCC), closepacked structure, uniaxially distorted crystalline, and the liquid structures.More details about DP training can be found in the supplemental materials [42].
Here we pay special attention to thermodynamic properties of equilibrium tungsten that are closely related to laser heating process.The melting temperature predicted by ground-state DPMD (3550 K [42]) is in consistence to the previous DFT-MD (3450 ± 100 K [46]) and Gaussian approximation potentials simulations (3540 K [47]), which confirms that the present PES can reproduce melting with DFT accuracy.Furthermore, the dependence of DPMD-predicted enthalpy on temperature along isobaric heating condition is shown in Fig. 2(a), and the experimental data agree very well with our DPMD predictions, especially in the liquid regime [44].The estimated enthalpy of fusion at the melting point (∆H m = 237 ± 20 kJ/kg) is also close to the DFT-MD values (250 kJ/kg) [43] and other experimental values (see Table .S1 in [42]).
Based on calculated thermophysical properties, we can determine the complete melting threshold ϵ m , which is the laser energy that is sufficient to drive the complete melt of the samples.We found ϵ m = 0.92 ± 0.04 MJ/kg, corresponding absorbed pump fluence is 53.0 ± 2.2 mJ/cm 2 for 30-nm-thick tungsten film [42].Such values are in agreement with the estimated values from experimental results [13,44,[48][49][50][51], in which energy density is approximately 0.94 MJ/kg (pump fluence of 53.8 mJ/cm 2 ).The density decrease at elevated tem- perature as shown in Fig. S3 [42], is also consistent to the experimental measurements [50][51][52].
The lattice dynamics, that requires high-order derivatives of PES, were further investigated.As shown in Fig. 2(b), the phonon dispersion curves of BCC tungsten under both equilibrium and non-equilibrium states are well-reproduced compared with the DFT results.In particular, comparing with the phonon dispersion at T e = 300 K, the directional phonon softening is observed along the H − N and N − Γ path in the first Brillouin zone at elevated electron temperature (T e = 10000 K), which can be attributed to the delocalization half-filled d bands [8,24].The depopulation of such a strong directional component in electronic bonding weakens the directional forces and may drive the crystalline structures towards to close-packed forms.These results indicate that the neural network PES can provide faithful prediction related properties in consistency to experiments or ab initio method.
Direct ab initio simulations of laser-driven dynamics.
It is stressed that our explicit electrontemperature-dependence PES can well capture the nonthermal nature of laser-excited metals.When implemented in TTM-DPMD framework, it allows us to establish a comprehensive understanding on laser-induced non-equilibrium states within ab initio accuracy.Recent time-resolved ultrafast electron and X-ray diffraction experiments collects direct quantitative structural information of laser-driven processes [13], providing a benchmark for the validation of the present newly-developed methods.Here, we apply TTM-DPMD to directly simulate the dynamic response of tungsten nanofilm under different absorbed laser energy densities.
In TTM-DPMD simulations, full-scale ab initio description in one-dimension of polycrystalline (PC) 30nm-thick tungsten nanofilm is considered, according to relevant UED experiment [13].For PC systems, large size included 752,650 atoms is used to describe crystal grains with random shapes, orientations, and different types of boundaries.The size of each grain ranges from ∼ 5 nm to ∼ 7 nm and each grain contains more than 10 4 atoms (totally reaching to millions of atoms), which cannot be achieved by the traditional time-dependent DFT simulations.Moreover, extra 30 nm vacuum space perpendicular to laser incident direction is set to allow free surface response to the internal stress relaxation, and extra spring forces are introduced for atoms in the bottom regime relating to their initial lattice site, to present bonding to the substrate (see supplementary materials [42]).
Considering the ballistic transportation of excited electron in tungsten (the mean free path ∼ 33 nm), we assumed the uniform deposition of laser energy with relatively low energy density of 0.08 MJ/kg (corresponding to absorbed laser fluence of 4.8 mJ/cm 2 ).In this case, a moderate two-temperature state is created at the initial stage, where maximum electron temperature can reach to 4400 K. Through electron-ion energy exchange, the system quickly reaches thermal equilibrium (T e = T i ∼ 920 K) at t = 5 ps.
The structure factor is calculated to extract the decay dynamics of Laue diffraction peak (LDP) [42,53], which is an important quantity to diagnose the structural dynamics in experiments [13].As shown in Fig. 3(a), based on TTM-DPMD simulations with the inclusion of laserdriven excited states, the temporal evolution of normalized intensity of (211) LDP agrees well with UED measurements.Conversely, the results from simulations by ground-state PES deviate from experimentally measured values significantly.It is interesting to say that the thermal process (T e = T i ) exhibits remarkably slower decay dynamics than the process with excited states (T e ≥ T i ) upon such laser fluence.By further investigating the lattice vibration of bulk tungsten, we note that even under moderate non-equilibrium state (T e = 5000 K), a relative increase of over 10% in mean square displacement (MSD) can be observed under isobaric heating condition, as shown in Fig. 3(b).The enhancement of lattice vibration can be attributed to the hot-electroninduced phonon softening (Fig. 3(c)).Such nonequilibrium and nonthermal effects therefore modify the dynamics of diffraction signals according to Debye-Waller formula [42], in which the decay of LDP is relating to temporal evolution of lattice temperature and temperature dependence of MSD.The quantitative consistency between our simulations and experiments validates our model further, and then provide a chance to further elucidate laser excitation effects.
By increasing laser energy density up to 0.80 MJ/kg, the irradiated tungsten nanofilm starts with more severe nonequilibrium states (T e = 11200 K, T i = 300 K).As presented in Fig. 4(a)(b), the evolution of tungsten nanofilm predicted by ground-state PES (T e = T i ) is a purely thermal process governed by electron-ion coupling.With increased lattice temperature, the system firstly evolves along the equilibrium isochore in the first 4 ps, where the ionic kinetic pressure accumulates to ∼ 10 GPa.Then the thermal pressure is gradually released due to existence of free surface.Although the thermal expansion process leads to temperature and density decrease, the gradient in thermodynamic profile is slight and the whole system can be considered as homogeneous.
When laser-induced changes in the PES is included (T e ≥ T i ), the thermodynamic pathway and thermodynamic profile is totally different.As shown in the Fig. 4(c)(d), the ultrafast excitation of electrons results in the buildup of extra pressure on a sub-picosecond timescale (more details in Fig. S7).Such hot-electroncontributed pressure increases monotonically with increased laser energy density, from ∼ 1 GPa with T e,0 = 4400 K (ϵ = 0.08 MJ/kg) to ∼ 17 GPa with T e,0 = 11200 K (ϵ = 0.80 MJ/kg).The tungsten nanofilm then quickly responds to this nonthermal internal stress, triggering anisotropic volume relaxation dynamics.As a result, a significant inhomogeneity is demonstrated in the thermodynamic profiles.In Fig. 4(d), the propagation and reflection of stress waves can be identified with velocity of ∼ 4.3 km/s, accompanied with the density decrease of ∼ 1 g/cm 3 .With existence of free surface, the build up and uniaxial relaxation of nonthermal stress can strongly influence the thermodynamic pathway especially under high laser fluence, which cannot simply be assumed to be isochoric or isotropically isobaric.We highlight that such real-time material response captured by TTM-DPMD simulation provides unique insights into previous controversial issues on nonthermal behavior of laser-excited matter [19,20].

DISCUSSION
In this work, we developed the deep learning model to perform large-scale ab inito simulations on the laserinduced atomistic dynamics, with quantum accuracy on the non-thermal effects.To validate the accuracy, special attention is paid to recent experiments.We successfully reproduce the experimental data with our model.It is therefore verified that the laser-excited states have profound effects on the thermodynamic evolution and structural transformation dynamics.More importantly, the combination of deep learning techniques with hybrid continuum-atomistic approach bridges the theoretical method and experimental observations, providing a new path to establish accurate and complete understanding of the atomistic dynamics under ultrafast laser interactions.

METHOD
DP training.The ETD-DP models for tungsten are generated with DeePMD-kit packages [54] by considering T e as atomic parameter.Deep Potential Generator (DP-GEN) [36], has been adopted to sample the most compact and adequate data set that guarantees the uniform accuracy of ETD-DP in the explored configuration space.We consider BCC structure (54 atoms) and liquid structure (54 atoms) as the initial configurations and run DPMD under NVT and NPT ensemble (both isotropic and uniaxial constrains are considered), where temperatures ranges from 100 K to 6000 K, pressure ranges from -15 to 60 GPa, and corresponding electronic temperature ranges from 100 K to 25000 K.The training sets consist of 6366 configurations under equilibrium condition (T e = T i ) and 6820 configurations sampled under two-temperature state (T e > T i ).
For DP training, the embedding network is composed of three layers (25, 50, and 100 nodes) while the fitting network has three hidden layers with 240 nodes in each layer.The total number of training steps is set to 400 000.The radius cutoff r c is chosen to be 6.0 Å.The weight parameters in loss function for energies p e , forces p f , and virials p V are set to (0.02, 1000, 0.02) at the beginning of training and gradually change to (1.0, 1.0, 1.0).
The self-consistency calculations are all performed with the VASP package [55].
TTM-DPMD simulation setting.We perform TTM-DPMD simulations with LAMMPS package [59] through modified EXTRA-FIX packages [23].The elec-tronic heat capacity is calculated by individual DFT calculations C e = T e ∂Se ∂Te , which is consistent with previous calculations [60].The electron-phonon coupling factor is set to constant (G 0 = 2.0 × 10 17 W m −2 K −1 ) according to relevant ultrafast electron diffraction experiments [13].The electron thermal conductivity is described by the Drude model relationship, κ e (T e , T i ) = F C e (T e )τ e (T e , T i ), where v F is Fermi velocity and τ e (T e , T i ) is the total electron scattering time defined by the electron-electron and electron-phonon scattering rates, 1/τ e = 1/τ e−e + 1/τ e−ph = AT 2 e + BT i .The coefficients A = 2.11 × 10 −4 K −2 ps −1 , B = 8.4 × 10 −2 K −1 ps −1 , v F = 9710 Å ps −1 are adopted [61].The duration of laser pulse is set to 130 fs.Since the mean free path of laser-excited electrons is ∼ 33 nm in tungsten [62], the electrons are heated uniformly due to the ballistic transport.Therefore, optical penetration of laser energy can be neglected for simplicity.
For atomic system, the simulation size of polycrystalline sample is set to 30 nm×20 nm×20 nm, containing 752650 atoms, with extra 30 nm vacuum space along the x direction to minic the free boundary condition.Extra spring forces are introduced for atoms in the bottom 5 Å relating to their initial lattice site to present bonding to the substrate.
phonon spectra calculation.To validate the accuracy of ETP-DP model, we investigate the lattice dynamics that need high order derivatives of PES.We use finite displacement method to calculate the phonon dispersion with ALAMODE package [63] as a postprocessing code.The forces are calculated in 5 × 5 × 5 supercell with cell lattice parameter a 0 = 3.17104 Å.The atomic displacement is set to 0.01 Å, and the interatomic force constants are extracted from KS-DFT and DPMD calculation respectively.The dynamical matrices are derived from these force displacement data to obtain phonon dispersion spectra.
Ultrafast electron diffraction pattern.To extract the decay of Laue diffraction peak (LDP) intensities as in the UED experiments, we performed the ultrafast electron diffraction simulations with DIFFRACTION package [53] to obtain the structure factor S(Q x , Q y ) defined as follows, where Q = (Q x , Q y , Q z ) the wave vector, f i the atomic scattering factor, λ the wavelength of incident electron, r i the coordinates of atom i.Here, simulated 3.2MeV electron radiation (λ ∼ 0.34 pm) is used to create selected area electron diffraction (SAED) patterns according to relevant experiments [13], and the SAED patterns aligned on the [100] axis (Q z = 0) are constructed by selecting reciprocal lattice points intersecting a 0.01 Å−1 thick Ewald

Figure 1 .
Figure 1.Schematic diagram of workflow for efficient and accurate simulation of laser-driven atomistic dynamics.(a) ETD-DP model.Te is the electron temperature, regarding to the electron occupation distribution.The free energy A, force F, virial Ξ, electronic entropy Se, and electronic heat capacity Ce can be inferred through backpropagation algorithm.(b)iterative concurrent learning scheme is used to efficiently sample atomic configurations for a wide range of equilibrium and non-equilibrium conditions.(c) hybrid atomistic continuum approach.The evolution of electron subsystem allows atomistic system transits between different PES, and the Langevin thermostat is introduced to mimic non-adiabatic energy exchange between electron and lattice.

Figure 2 .
Figure 2. Validating accuracy of ETD-DP model.(a) Temperature dependence of enthalpy under isobaric heating (p = 1 bar) with the reference temperature of 300 K.The blue line, the black cross, and grey square denotes the DPMD results, previous DFT-MD prediction [43], and isobaric expansion experimental data [44] respectively.(b) phonon dispersion of laser-excited tungsten (ρ0 = 19.15g cm −3 ).The black cross and white squares represent the individual KS-DFT calculation and experimental measurements [45].

Figure 3 .
Figure 3.Capturing nonthermal effect with TTM-DPMD approach.Comparison of (a) temporal evolution of (211) diffraction peak intensity in structure factor under absorbed laser energy density of 0.08 MJ/kg, compared with experimental data[13].(b) Temperature dependence of mean square displacement with isobaric constrains and (c) phonon density of states (PDOS), obtained under equilibrium condition (blue) and nonequilibrium condition (orange).

Figure 4 .
Figure 4. Hot electron modifies the thermodynamic pathway.Comparison of (a)(c) thermodynamic pathway (b)(d) temporal evolution of thermodynamic profile of nanofilm, predicted by ground-state PES and laser-excited PES, respectively.In (a)(c), the red arrows indicate the evolution path of selected regime (z = 14.0 nm) in the tungsten nanofilm, and the thermodynamic state is highlighted by colored stars every 1 ps.In (d), the black dashed lines are used to highlight the propagation of stress waves, whose slope represents a constant propagation speed of ∼ 4.3 km/s.