# A local two-temperature model for electronic heat conduction in molecular dynamics simulations

###### Abstract

We present a local two-temperature molecular dynamics (MD) model for coupled simulations of electronic and phonon heat conduction in nanoscale systems. The proposed implementation uses a master equation to perform local heat conduction of the electronic temperature eschewing the need to use a basis set to evaluate operators. This characteristic allows us to seamlessly couple the electronic heat conduction model with molecular dynamics codes without the need to introduce an auxiliary mesh. We implemented the methodology in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code and through multiple examples, we validated the methodology. We then study the effect of electron-phonon interaction in high energy irradiation simulations and the effect of laser pulse on metallic materials. We show that the electrons and phonons exchange energy properly, and the energy is conserved. The parallel performance and some aspects of the implementation are presented.

###### keywords:

Nanoscale heat transport, electronic heat conduction, two temperature model, kinematic mean field theory, irradiation damage, laser ablation.## 1 Introduction

Classical Molecular Dynamics (MD) simulations have been used widely in materials physics to study the lattice and defect dynamics at the atomic scale in both equilibrium and non-equilibrium conditions. In classical MD, the thermal conductivity due to lattice contribution is considered, since the phonons are explicitly resolved. This makes MD simulations one of the most used methods to measure thermal properties of materials where heat is mainly transported by phonons. Unfortunately, the electronic contribution is neglected in MD, and thus the study of materials where electronic heat transport is important cannot be performed with accuracy in MD. This render MD simulations well suited but also, essentially limited to the simulation of low level of electronic excitation and ionization events, semiconductors and insulators. For example, phenomenon like swift heavy ion irradiation, short-pulse laser interaction, and metal-semiconductor joints where electronic energy loss/gain is significant cannot be described accurately using MD.

A preferred method to address this issue is two-temperature molecular dynamics (2T-MD) model. In this approach, a continuum model of the electronic energy transport is combined with classical MD in a concurrent multiscale scheme. In this model, the temperature evolution in electronic subsystem is described by a heat diffusion equation Kaganov and Tanatarov (1957); Ivanov and Zhigilei (2003),

(1) |

In the above equation, , , and are electronic heat capacity, thermal conductivity, local electronic and ionic temperature, respectively. is a electron-phonon coupling constant between electron and phonons.

One appealing aspect of the 2T-MD is that it can concurrently simulate energy exchange between phonons and electrons during a simulation. This has been achieved by a number of means, including finite differences (FD) Ivanov and Zhigilei (2003); Duffy and Rutherford (2007); Schäfer et al. (2002); Rutherford and Duffy (2007) and finite elements method (FEM) Wagner et al. (2008). In implementations using FD, the atomic system is divided into small cells (voxels) that generate a mesh and the FD method is used to evaluate the Laplacian of the temperature and solve Eq. 1 for each voxel. Similarly, in implementations using FEM, the MD system is embedded into a FEM mesh, where some of the atoms generate the nodes of this mesh. Thus, the coupling is done by passing the MD simulation data into the FEM mesh, which solves Eq. 1. Other methodologies have coupled both phonons and electrons using a heat conduction equation as Eq. 1 while completing removing the lattice dynamics, leading to free-energy type methods to describe thermo-mechanical response of materials Ariza et al. (2012); Ponga et al. (2015, 2017).

Regardless the choice of the method to evaluate the Fourier law for the electronic temperature, previous implementations have been using finite size sub-divisions of the system and linked the behavior of a cluster of atoms into the FD grid or FEM nodes. However, in applications where large amounts of energy are driving systems out of equilibrium, it is unclear the validity of Fourier law which is only valid under weak non-equilibrium conditions, e.g., small temperature gradients. Secondly, due to the large temperature gradients that happen in short length and time scales, the spatial discretization of the mesh might not be adequate to properly describe the temperature variation. This indicates that heavily excited atoms and electrons can quickly exchange energy in local regions of the system and thus, small resolution is required to resolve this issue. Hence, recent studies have been indicated the need for local 2T-MD methods that can accurately account for this Zarkadoula et al. (2016, 2017). The 2T-MD approach has been applied to a number of technological problems of interest, most notably in laser ablation and high energy ion irradiation simulations Leino et al. (2018); Shugaev et al. (2017), shock induced melting Koči et al. (2006), and energy transport across semiconductor-metal interfaces Wang et al. (2012).

In this work, we develop a local version of the two-temperature molecular dynamics (2T-MD) method by replacing the classical Fourier law by a Fokker-Planck equation to simulate the exchange of energy between electrons. The Fokker-Planck equation, used in the form of a master equation, computes the probability that two adjacent electronic sites exchange energy through electronic heat flow. This eschews the need of using an auxiliary mesh in order to evaluate the Laplace operator and retains, at the same time, the smallest resolution that one can achieve in MD systems, e.g., the atomic resolution. Another remarkable aspect of this approach is that can be applied to situations where the system is far from the thermodynamic equilibrium, e.g., large temperature gradients. Additionally, since the method does not need an auxiliary mesh, it can be seamlessly coupled to MD codes. Hence, we implemented the 2T-MD method into the LAMMPS code Plimpton (1995) and make it available through its online site l2T (????).

The manuscript is organized as follows. First, we start by introducing the formulation and performing the energy balance to determine the coupling between electrons and phonons. Next, we demonstrate the methodology by using a simple one-dimensional example using the master equation approach to simulate the electronic and phonons heat conduction. We compare this analytical example with the newly implemented methodology in the LAMMPS code. We investigate the conservation of the energy in the system and the energy exchange between electrons and phonons. The parallel performance of the method is also studied. Finally, we study two problems of technological interest, an irradiation damage problem in Ni using a 50 KeV cascade simulations; and a laser spallation simulation in Ni. We conclude the work with the main outcomes of this work.

## 2 Methodology

We now describe the 2T-MD method. Let us now consider a system of particles at finite temperature . This temperature is the contribution of two players, e.g., phonons and electrons. MD simulations explicitly model the phonons as a result of atomic interactions accounted for by the interatomic potential. Unfortunately, the electronic contribution to the temperature is neglected and difficult to model. The two-temperature model assumes that the temperature field can be split into two interacting subsystems contribution, e.g., phonons and electrons. Thus, one has an electronic temperature, , and a lattice temperature given by the classical definition, i.e., , where is the total kinetic energy of the atoms, and are the mass and the velocity of the atom, respectively. is the Boltzmann constant. We assume that associated with each atom, there is a local electronic temperature, i.e., that can fluctuate from atom to atom. Introduce the normalized temperature field and ; where and are the maximum and minimum allowed electronic temperatures in the system, respectively.

Now, we are interested in simulating heat transport in molecular dynamics simulations and include the effect of both phonons and electrons. In order to simulate heat conduction of the electronic temperature, we follow the approach proposed by Ponga and Sun Ponga and Sun (2018) and modified it to consider coupling between electrons and phonons. Thus, we use the following modified master equation

(2) |

where is a pair-wise exchange rate thermal coefficient for the electronic temperature , and is the normalized difference in the electronic temperature between the and atoms, and is the electronic heat capacity . The first term on the right hand side of Eq. 2 measures the rate of electronic energy exchanged between two adjacent particles and it is arbitrarily chosen to be around the nearest neighbors of the particle, which makes the model fully local. According to Ponga and Sun Ponga and Sun (2018), the pair-wise exchange rate thermal coefficient can be computed from an asymptotic expansion of the master equation into the continuum limit leading to a relation with the (macroscopic) electronic thermal conductivity () as

(3) |

where is the dimension of the system, is the coordination number, is the Burgers vector of the material.

The second term on the right hand side of Eq. 2 is a traditional coupling term between the electrons and phonons. is a constant that defines the strength of the coupling between the atoms and the electrons and has units of and is the temperature difference between the electronic and the lattice temperature of the particle. denotes a local lattice temperature of the atom defined as

(4) |

where denotes the velocity of the atom between a cut-off radius, i.e., . The cut-off radius is arbitrarily selected to be the cut-off of the potential function.

The amount of energy exchanged per unit of time between the particle and the lattice is . This flux of energy needs to be introduced (removed) in the lattice by modifying its kinetic energy. In order to do so, we introduce this energy flux by modifying the dynamics of the lattice Ivanov and Zhigilei (2003), i.e.,

(5) |

where is acceleration of the atom, is the force acting on the atom and is the positive damping force that appears due to the lattice heating (cooling) coming from the electronic temperature. is a local coupling coefficient of the atom whose strength needs to be computed at each time step. In order to do so, let us now analyze the amount of energy exchange by the electronic temperature per unit of time. Letting be the atomic volume, the energy exchange between the lattice and the atom is

(6) |

The rate of work (power) done by the damping force is

(7) |

It is easy to see that has to be computed at each time step in order to conserve energy in the system due to the exchange of energy between the electrons and the phonons. The rate work done by the damping force has to be equal to the amount of energy exchange by the electrons in Eq. 6. leading to

(8) |

In order to avoid numerical instabilities, Eq. 8 is modified to take into account the kinetic energy of a set of particles between the cut-off radius of the potential. In this way, we avoid situations where the kinetic energy goes to small values close to zero.

Before closing this section, we make some remarks. First, we notice that the formulation uses local electronic and lattice temperature, and therefore, the energy exchange happens locally in the atoms. The strength of the coupling between electrons and phonons, given by changes atom by atom and it is a function of the time. Hence, in the current form, the electronic temperature serves as a local thermostat that tends to equilibrate the electronic and lattice temperatures.

### 2.1 Implementation

We have implemented the proposed 2T-MD model in the LAMMPS code. We have done so by including an electronic temperature variable per atom. Then, the master equation Eq. 2 was implemented using a fix operation, which is called every integration step to evaluate the electronic temperature. The electronic temperature for each atom is calculated within a fixed cut-off radius that can be given by the user in the script. Once that the electronic temperature is updated, the coupling strength between phonons and electrons (Eq. 8) is computed and the equation of motion (Eq. 5) is modified to include the damping force at the atom. The fix was implemented to work for parallel calculations using the Message Passing Interface (MPI) library. Since the electrons thermal diffusivity is very high, the time step for the integration of Eq. 2 has to be a fraction of the integration time step in MD simulations ( 1 fs). Thus, during a lattice integration step, we actually allow to perform multiple integration steps of the heat equation, and this can be selected by the user from the simulation script. Additionally, in large energy cascade simulations ( 10 KeV), the velocity of the ions could be extremely high, and thus sometimes it is desirable to have a variable integration time step for the lattice. Our implementation is capable of handling this if the user specifies it. Cascade simulations described in Section 4 were carried out using a variable time step.

## 3 Heat Conduction in a Ni bar

Our first test consists of the classical 2T model, where heat diffusion is modeled for both phonons and electrons. In the classical approach, a Fourier law is used to model electron and phonon transport, i.e,

(9a) | |||

(9b) |

In our test, we first replace Eq. 9a by the master equation in Eq. 2 in order to show the validity of the approach. We solved the lattice heat conduction on a 1D-Ni bar using Eq. 9b and electronic heat conduction using the master equation (Eq. 2).

The equations were solved using a custom code in MATLAB. The mesh points coincided with the atomic positions, which are separated a distance where is the lattice parameter. The evolution of the lattice temperature was carried out using the FD method. We used a three-point stencil, the temperature was integrated using a Euler forward algorithm and periodic boundary conditions were applied along the -axis for both electrons and phonons. The evolution of the electronic temperature was carried out with the master equation (Eq. 2) using a nearest neighbor interaction.

Both lattice and electronic subsystems were kept at 300 K at the beginning of the simulation. A hot zone was created by applying 6000 K to the electronic subsystem to a span of length 2 at the center of the bar. The dimensions of the Ni-bar was = 70 = 25 nm, = = 2 = 0.71 nm, where = 0.3562 nm is the lattice constant of fcc Ni. The material constants used for the description of the electronic subsystem are, electron-phonon coupling () = 3.6 10 J m K s, electronic heat capacity () = 3.6886 10 J mK and thermal conductivity () = 91 J s m K Ivanov and Zhigilei (2003); Hohlfeld et al. (2000). Although, all these quantities are in fact temperature-dependent, we used a constant value through our work for simplicity. However, we note that using temperature dependent coefficients can be easily handled in our method. In literature both temperature independent and dependent quantities are used in 2T-MD simulations Ivanov and Zhigilei (2003); Zarkadoula et al. (2016). Hence, we choose constants that are temperature independent and calculated at 300K. The pair-wise exchange rate thermal coefficient = 2019 ps was calculated using Eq. 3, where = 3, = 12 and = 0.247 nm. The heat capacity () = 400 J Kg K and density () = 8890 Kg m were used for the lattice heat diffusion. To compare the classical approach with MD, a similar study was also done using the same parameters in the newly implemented 2T-MD method in LAMMPS. A constant time step of 0.1 fs was used for both simulations.

The temperature profiles of the classical approach and 2T-MD models are shown in Fig. 3. Let us now focus our attention on the classical approach, shown in Fig. (a)a. In the figure, we observe how the heat is being diffused following a Gaussian profile for the lattice and electrons, as predicted by the analytical solution. We see that the master equation predicts the same response than the Fourier law (difference is less than 0.1 %) when the temperature gradients are small (Fourier results have been omitted in the plot for simplicity). We also see that some of the energy injected into the electrons is transferred to the lattice, and then the lattice temperature is increased with time. After 3 ps, both electronic and lattice temperature are in equilibrium.

Let us now describe the 2T-MD model for the same example. Fig. (b)b shows the time evolution of the electronic temperature in the MD sample. We see a similar time scale for heat diffusion and relaxation after approximately 3 ps. Remarkably, once the system is in equilibrium, both electronic and lattice temperatures reach the same value, which is very close to the classical model.

We make the following remarks regarding the MD model. First, the example is fully three-dimensional besides that the heat is diffusing in one direction. Second, since the lattice temperature is difficult to obtain for individual atoms, we could not plot the local lattice temperature of the sample. Having said that, the final mean temperature of both electronic and lattice subsystem reached equilibrium at 310 K. In metals, the phonon contribution to the thermal conductance is usually negligible and difficult to measure experimentally compared to electronic contribution. In general, thermal conductivity for metals means electronic thermal conductivity. For the example, the thermal conductivity for both subsystems was taken equal to the electronic one, . We have also tested with several lower values for lattice thermal conductivity; this increases the time for both systems to equilibrate and changes the transient solution making it more similar to the 2T-MD. Thus, the similarities of the temperature profile, relaxation time and final equilibrium temperature confirm the good qualitative agreement between the classical approach and 2T-MD models.

The energy balance of the 2T-MD model is shown in Figure 4. The total energy (TE) is conserved in the simulation. Initially, the electronic energy was high as a high temperature was given to the electronic subsystem. The energy was gradually transferred to the atomic subsystem and increasing both kinetic and potential energies. We see that after 3 ps, the energy exchange reduces significantly and the three quantities remain approximately the same, as expected. Figure (8(a-c)) also show the snapshots of electronic temperature distribution of 2T-MD simulation at different time steps. We see the ability of the model to retain an atomic resolution for the electronic temperature.

Figure 9 shows the performance of the implementation in the LAMMPS code for the uniaxial bar problem with a larger number of atoms, i.e., . As we can see, the implementation has good scalability up to 32 processors (80% efficiency). However, after this point, the communication becomes the bottleneck and the efficiency gets reduced ( 60% for 96 processors). This is due to two main reasons, the number of atoms per processor reduces to less than 1000 atoms/processors, and therefore, it is difficult to scale since the ratio between computing over communication time is reduced significantly. Second, in order to the electronic heat diffusion, we need to communicate information to the ghost atoms in the simulation, which is an overhead with respect to MD. Thus, this reduces the performance. Reduction of communication can be achieved by using hybrid implementations with threads in OpenMP. This will be addressed in the future.

## 4 Irradiation damage modeling

Now, we turn our attention to simulate far from equilibrium phenomenon like ion irradiation, where electronic excitations and ionizations are significant. Irradiation in metallic systems has been studied in several works including only the effect of the phonons Aidhy et al. (2015); Ullah et al. (2016). Here, we include the effect of the electrons and phonons using the proposed approach. We performed 50 KeV Ni cascade simulations to model irradiation induced damage of elemental Ni using two different models (classical MD and 2T-MD). The cascade simulations were done using the newly implemented 2T-MD method in LAMMPS Plimpton (1995) code. We used the embedded atom method (EAM) based potential of Bonny et al. Bonny et al. (2013) to model elemental Ni, however, the 2T-MD method works with all potentials available in LAMMPS. This potential was developed to study the production and evolution of radiation defects, particularly, point defects. The simulation cell size was 343434 nm containing 3,500,000 atoms. The lattice and electronic subsystems were initially at equilibrium having a temperature of 300 K. Periodic boundary conditions were applied to all directions to replicate an infinitely large bulk system. A fixed time step of 1 femtosecond (fs) was used to relax the system and a variable time step was used having the maximum time step restricted to 0.01 fs for cascade simulation. At the beginning of the simulation, a recoil energy of 50 KeV was given to a central atom at a random direction. The simulation run at 300 K for 16 ps. A total of 12 such simulations were performed for each simulation model with the same recoil direction to collect damage statistics.

Figure 13(a) shows the time evolution of mean lattice and electronic temperatures of the simulation system. Both subsystems are at equilibrium for sometime (3 ps) before the recoil starts. The lattice temperature rises immediately to the peak with the recoil then decays very fast ( 5 ps) to equilibrium temperature. Whereas, the mean electronic temperature rises slightly slower than the lattice temperature. The electronic and lattice temperatures are equilibrated after 6 ps ( 3 ps from the moment recoil starts). The maximum electronic temperature is shown in Fig. 13(b) is the highest atomic value calculated at the corresponding time step. We saw that while the mean electronic temperature remains low (330 K), there are heavily excited electrons in small regions that can reach temperatures in excess of 6500 K. The ability of the model to predict such large temperature gradients is remarkable. The maximum lattice kinetic energy is shown in Fig. 13(c) indicates two important things. First, the rate of change of the kinetic energy is much faster for the 2T-MD model, since the effect of electrons are taking into account. As a result, the maximum kinetic energy in the sample for the regular MD model remains large (35 KeV) for approximately 0.25 ps, while in the 2T-MD model, the maximum kinetic energy is reduced immediately. For instance, the kinetic energy in the 2T-MD model is reduced from 35 KeV to 1 eV in around 0.3 ps while in the regular MD model this takes 0.5 ps. This impacts the generation of defects, as we will explain later on. Second, the plot shows that local atoms in the sample are heavily excited as indicated by their kinetic energy and thus, are far from equilibrium in both MD and the 2T-MD models. While the atomic interactions are handled locally through the interatomic potential it was not possible until now to have a local treatment of the electrons. Thus, the current implementation has the ability to model problems that are very far from the thermodynamic equilibrium accounting for electronic contributions in a local sense to diffuse heat in the sample.

Figure 16, shows the electronic temperature distribution of the 50 KeV recoil event after 0.03 ps together with atoms showing exceeding amounts of kinetic energy. We see that only a small fraction of atoms have kinetic energy in excess of 0.4 eV. However, the electronic temperature can reach exceedingly large values as shown in the color map. The rainbow color map shows the radial distribution of the electronic temperature on the plane, characteristic of a Gaussian type distribution. To have a better understanding, 16(b) shows the temperature distribution along the narrow slit (brown line) on the plane. We observed a large temperature gradient ( 150 Knm) from the central hot zone to the surrounding. This is a clear indication of a large temperature gradient happening in a very small length scale, denoting the local character of the energy exchange in the system. Our analysis illustrates the need for models that are able to have a fine resolution in order to properly predict the electronic effects in these situations.

The Voronoi cell method was used to identify the final damage production Nordlund et al. (1998). In this method, Voronoi polyhedra were centered on each atom position in the ideal fcc crystal and compared with atoms in the defective crystal. Polyhedra with no atoms were labeled as vacancies and polyhedra with 2 or more atoms were labeled as interstitials. The total number of Frenkel pairs (FPs) are 102 7 and 144 14.5 for 2T-MD and classical MD cascade simulations, respectively. The number of FPs and the standard error of the mean is calculated form 12 simulations for each simulation condition. The low number of FPs in 2T-MD simulation is expected due to energy transfer to the electronic subsystem as shown in the kinetic energy plot (Fig : (c)c). Figure 19 shows the final point defects of a representative classical MD simulation and the 2T-MD simulation having same recoil direction. The features of cascade damage is well captured in the simulation. It is visible from the snapshot that the interstitials are more aggregated due to high mobility and diffusivity compared to vacancies.

## 5 Laser ablation

Laser simulations were performed by giving a Neumann boundary conditions to the simulation cell. This boundary condition specifies the rate of change of the electronic temperature by giving a heat flow that follows the following expression Ivanov and Zhigilei (2003)

(10) |

Here, is the maximum absorbed intensity [temperaturetime], reflects the amount of energy absorbed, and is the reflection parameter, and are the first components of the position vector of the -atom and the free-surface end, where the laser pulse is applied, respectively. is the penetration length used to simulate ballistic effects of the heat conduction in the electrons, is the time at which the pulse is applied and is the standard deviation of the Gaussian profile. The heat flow intensity can be related to the laser fluency by taking and , and , then

(11) |

where is a characteristic time and is a factor that accounts for the atomic positions.

In our simulations, we selected Jm ( Kps), the half-width duration and the duration of the pulse were ps, and ps, respectively. The penetration length was selected to be nm. The simulation cell was containing 160,000 atoms and it was carried out with a constant time step of ps. The material’s constants and potential were specified in Section 4.

Fig. 23 shows the time evolution of the change of energy between electrons, phonons, and potential energy; the temperature evolution with and without including the center of motion; and the fraction of the atoms in a particular atomic structure. Fig. 23(a) shows that initially, a large amount of energy is deposited on the electrons due to the laser pulse. However, since heat conduction is allowed, this energy is quickly distributed among the electrons and phonons. This is indicated by the raise in the kinetic energy. After 2 ps, the phonons reach a constant level of energy, and the potential energy increases significantly. The constant kinetic energy, with the increase in the potential energy, is an indication that part of the material is absorbing energy at constant temperature and changing its physical state from solid to liquid. This is usually known as latent heat. Around ps, another large change in the kinetic energy with a decrease of the potential energy is observed, which indicates a large portion of the atoms in the system being melt and increasing their temperature. Then, after ps, the sample stabilizes and the three quantities remain approximately the same, an indication of a stable equilibrium in the system.

Fig. 23(b) shows the evolution of the atoms according to their local structure using the common neighbor analysis Stukowski (2012). Initially, the sample contains mainly fcc atoms (99.9%) and few atoms in other types due to surface effects (0.1%). As the laser pulse gives energy to the sample, the fcc fraction quickly decreases to 50% in about 5 ps, while the fraction of atoms in hcp and bcc structures changes from 0 to 20% in about 4 ps, to get reduced again to 0 after 8 ps. This is a clear indication that a wave front is traveling in the sample. Interestingly, the fraction of atoms in other types (amorphous) structure increases significantly after 5 ps. This is due to the fact that they have absorbed sufficient energy to change phase, from solid to liquid. This is also evident in the amount of kinetic energy in the sample (Fig. 23(a)), that increases during this period. Fig. 23(b) indicates that after 10 ps, the fraction of atoms in the amorphous phase remains constant and no further changes in the sample are seeing.

Fig. 23(c) shows the time evolution of the lattice temperature for the whole sample, with or without the velocity of the center of mass. As the thin film is free to move, the temperature of the center of mass is more representative of the system’s temperature. We observed that the temperature of the sample is significantly increased during the melting, and then part of this energy is inserted back in the system in the form of potential energy.

Figure 24 shows the hydrostatic pressure evolution computed using the virial tensor for the sample as a function of the time. We observed that a large compressive shock wave is propagated through the sample as can be seen in Fig. 24. After ps, the hydrostatic pressure changes direction and reflected at the free-surface as a tensile shock wave, nucleating and growing nanovoids in the sample (see Fig. 31(b) ). This, ultimately, leads to the fragmentation of the sample, as can be clearly see in Fig. 31(c). Figs. 31(d-f) show the atoms colored according to their local atomic structure for the previous snapshots shown in Figs. 31(a-c). We also observed that the energy is conserved, although there is some small drift as expected.

## 6 Conclusions

We have developed, implemented and validated a local version of the two temperature model coupled with MD simulations (2T-MD). The methodology uses a master equation to compute effective rates of energy exchange in the electronic subsystem and it is coupled to the phonon vibrations of the atoms through classical MD. We have demonstrated that the energy in the system is preserved, that the electrons and phonons exchange energy properly, and the system reaches equilibrium upon sufficient time. In the cascade simulations, we found that while the mean average electronic temperature remains relatively small, the maximum temperature in the system can reach exceedingly large values. For instance, in the 50 KeV simulation, the mean electronic temperature increases up to 450 K, while the maximum reaches about 7,000 K. We also showed that the temperature gradients can be very large of the order of 150 Knm, as demonstrated in the electronic temperature across a section of the simulation. Such a large gradient, illustrates the ability of the master equation and the current approach to transport energy at the atomic scale including situations where the system is far from the thermodynamic equilibrium. Interestingly, the introduction of electrons in the irradiation damage simulations leads to a significant reduction of the numbers of Frenkel pairs with respect to MD, illustrating the need for using the 2T-MD model.

We also illustrated a laser ablation simulation, where a large amount of energy was introduced in the sample, generating melting of the sample and shock waves that generate nanovoid nucleation, growth and coalescence of voids, leading to an ablation of the material. The current approach provides an alternative and mesh free way for coupling electrons with phonons in MD software. This is important because the implementation can be seamlessly implemented in multiple MD codes. We have done so using the LAMMPS code, providing capabilities for simulating atomic systems with massively parallel clusters. Future directions include the coupling of phonons with magnons to simulate ferroelectric materials.

## 7 Acknowledgments

We gratefully acknowledge the support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grant under Award Application Number RGPIN-2016-06114 and the support of Compute Canada through the Westgrid consortium.

## 8 Data availability statement

The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations.

## References

## References

- Kaganov and Tanatarov (1957) M. I. Kaganov, I. M. L. L. V. Tanatarov, Relaxation between elecrons and crystalline lattice, Sov. Phys. JEPT 4 (1957) 173.
- Ivanov and Zhigilei (2003) D. S. Ivanov, L. V. Zhigilei, Combined atomistic-continuum modeling of short-pulse laser melting and disintegration of metal films, Physical Review B 68 (2003) 064114.
- Duffy and Rutherford (2007) D. M. Duffy, A. M. Rutherford, Including the effects of electronic stopping and electron-ion interactions in radiation damage simulations, Journal of Physics: Condensed Matter 19 (2007) 016207.
- Schäfer et al. (2002) C. Schäfer, H. M. Urbassek, L. V. Zhigilei, Metal ablation by picosecond laser pulses: A hybrid simulation, Physical Review B 66 (2002) 115404.
- Rutherford and Duffy (2007) A. M. Rutherford, D. M. Duffy, The effect of electron?ion interactions on radiation damage simulations, Journal of Physics: Condensed Matter 19 (2007) 496201.
- Wagner et al. (2008) G. Wagner, R. Jones, J. Templeton, M. Parks, An atomistic-to-continuum coupling method for heat transfer in solids, Computer Methods in Applied Mechanics and Engineering 197 (2008) 3351 – 3365. Recent Advances in Computational Study of Nanostructures.
- Ariza et al. (2012) M. P. Ariza, I. Romero, M. Ponga, M. Ortiz, Hotqc simulation of nanovoid growth under tension in copper, International Journal of Fracture 174 (2012) 75–85.
- Ponga et al. (2015) M. Ponga, M. Ortiz, M. Ariza, Finite-temperature non-equilibrium quasi-continuum analysis of nanovoid growth in copper at low and high strain rates, Mechanics of Materials 90 (2015) 253 – 267. Proceedings of the IUTAM Symposium on Micromechanics of Defects in Solids.
- Ponga et al. (2017) M. Ponga, M. Ortiz, M. P. Ariza, A comparative study of nanovoid growth in fcc metals, Philosophical Magazine 97 (2017) 2985–3007.
- Zarkadoula et al. (2016) E. Zarkadoula, G. Samolyuk, H. Xue, H. Bei, W. J. Weber, Effects of two-temperature model on cascade evolution in ni and nife, Scripta Materialia 124 (2016) 6 – 10.
- Zarkadoula et al. (2017) E. Zarkadoula, G. Samolyuk, W. J. Weber, Effects of the electron-phonon coupling activation in collision cascades, Journal of Nuclear Materials 490 (2017) 317 – 322.
- Leino et al. (2018) A. A. Leino, G. D. Samolyuk, R. Sachan, F. Granberg, W. J. Weber, H. Bei, J. Liu, P. Zhai, Y. Zhang, Gev ion irradiation of nife and nico: Insights from md simulations and experiments, Acta Materialia 151 (2018) 191.
- Shugaev et al. (2017) M. V. Shugaev, I. Gnilitskyi, N. M. Bulgakova, L. V. Zhigilei, Mechanism of single-pulse ablative generation of laser-induced periodic surface structures, Physical Review B 96 (2017) 205429.
- Koči et al. (2006) L. Koči, E. M. Bringa, D. S. Ivanov, J. Hawreliak, J. McNaney, A. Higginbotham, L. V. Zhigilei, A. B. Belonoshko, B. A. Remington, R. Ahuja, Simulation of shock-induced melting of ni using molecular dynamics coupled to a two-temperature model, Phys. Rev. B 74 (2006) 012101.
- Wang et al. (2012) Y. Wang, X. Ruan, A. K. Roy, Two-temperature nonequilibrium molecular dynamics simulation of thermal transport across metal-nonmetal interfaces, Phys. Rev. B 85 (2012) 205311.
- Plimpton (1995) S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117 (1995) 1–19.
- l2T (????) l2t-md web site, https://mech-modsim.sites.olt.ubc.ca/software/l2t-md-a-local-two-temperature-model-for-electron-phonon-coupling-in-md/, ???? Accessed: 2018-10-25.
- Ponga and Sun (2018) M. Ponga, D. Sun, A unified framework for heat and mass transport at the atomic scale, Modelling and Simulation in Materials Science and Engineering 26 (2018) 035014.
- Hohlfeld et al. (2000) J. Hohlfeld, S.-S. Wellershoff, J. Güdde, U. Conrad, V. Jähnke, E. Matthias, Electron and lattice dynamics following optical excitation of metals, Chemical Physics 251 (2000) 237–258.
- Aidhy et al. (2015) D. S. Aidhy, C. Lu, K. Jin, H. Bei, Y. Zhang, L. Wang, W. J. Weber, Point defect evolution in ni, nife and nicr alloys from atomistic simulations and irradiation experiments, Acta Materialia 99 (2015) 69 – 76.
- Ullah et al. (2016) M. W. Ullah, D. S. Aidhy, Y. Zhang, W. J. Weber, Damage accumulation in ion-irradiated ni-based concentrated solid-solution alloys, Acta Materialia 109 (2016) 17 – 22.
- Bonny et al. (2013) G. Bonny, N. Castin, D. Terentyev, Interatomic potential for studying ageing under irradiation in stainless steels: the fenicr model alloy, Modelling and Simulation in Materials Science and Engineering 21 (2013) 85004.
- Nordlund et al. (1998) K. Nordlund, M. Ghaly, R. S. Averback, M. Caturla, T. Diaz de la Rubia, J. Tarus, Defect production in collision cascades in elemental semiconductors and fcc metals, Phys. Rev. B 57 (1998) 7556–7570.
- Stukowski (2012) A. Stukowski, Structure identification methods for atomistic simulations of crystalline materials, Modelling and Simulation in Materials Science and Engineering 20 (2012) 045021.