Two-Temperature Model

Introduction

Traditionally, the modelling of radiation damage has been confined to cases in which the radiation (i.e. a projectile) collides elastically with the nuclei of the target material. Collisions of this type are said to be dominated by nuclear stopping (energy transfer from the projectile to the nuclei). However, there are a broad range of scenarios (high energy collision cascades, swift heavy ion irradiation, and laser excitation) in which a significant portion (high energy collision cascades) or all (lasers and swift heavy ions) of the energy is transferred to the electrons in the target material. These cases, in which electronic stopping cannot be neglected, are impossible to account for using traditional molecular dynamics simulations (which assume equilibrium conditions between the target nuclei and electrons). This chapter describes the implementation of a hybrid continuum-atomistic implementation within DL_POLY_5 that incorporates these electronic excitations.

The model is based on the traditional two-temperature model (TTM) ⁶⁷. This model splits the nuclei and electrons into two separate but interacting subsystems, each evolving according to a modified version of Fourier’s heat equation. The continuum implementation of the nuclei (assumed to be in a lattice) in this form of the model is unable to track individual atomistic trajectories, thus information on superheating, recrystallisation, and pressure waves is lost. These limitations were overcome by Duffy and Rutherford ^30,31, by replacing the continuum representation of the nuclei with an MD cell (the details of which will be described in this chapter). This hybrid continuum-atomistic model, a combination of TTM and MD known as the two-temperature molecular dynamics (2T-MD) model, has been successfully used to model high energy cascades in iron ¹⁵⁰, ultrafast laser irradiation of gold nanofilms ²⁴, and swift heavy ion irradiation of silicon ⁵⁹.

Methodology

Electronic subsystem

The electronic temperature is treated as a continuum system, and evolves according to Fourier’s law of heat conduction:

(144)\[C_e(T_e) \frac{\partial T_e}{\partial t} - \nabla . [\kappa_e \nabla T_e ] = -G_{ep}(T_e - T_a) + G_s T^{\prime}_a + A(r,t),\]

where \(C_e(T_e)\) is the electronic volumetric heat capacity (equal to the product of specific heat capacity and density), \(\kappa_e\) is the electronic thermal conductivity, \(G_{ep}(T_e)\) the electron-phonon coupling, \(G_s\) the electronic stopping term (which is significant for atomistic/ionic temperatures \(T^{\prime}_a\) greater than a velocity cutoff, as defined later), and \(A(r,t)\) the temporal and spatial dependent source term. This equation describes how energy evolves in the electronic system as follows: energy is dumped into the electronic system via the source term (for swift heavy ions and laser excitation), \(A(r,t)\), the electronic volume-specific heat determines the electronic temperature (\(T_e\)) rise due to this deposition, the thermal conductivity describes how energy dissipates throughout the electronic subsystem, and the electron-phonon coupling determines energy transfer from the electrons to the MD cell (and is proportional to the temperature difference between \(T_e\) and the atomistic temperature, \(T_a\)). Equation (144) can be equated to the more general heat diffusion equation,

(145)\[\frac{\partial T}{\partial t} - \alpha \nabla^2 T = \frac{\dot{q}}{C},\]

\(T\) is temperature, \(t\) is time, \(\dot{q}\) is a heat source or sink, \(C\) is heat capacity, and \(\alpha\) is thermal diffusivity. This partial differential equation is solved using Euler’s method, which is a space-centred, forward-in-time integration algorithm. For a temperature \(T_{i}^{n}\) (at time step \(n\) and grid point \(i\)), the forward-in-time implementation can be Taylor-expanded and rearranged to:

\[\bigg(\frac{\partial T}{\partial t}\bigg)_i = \frac{T_{i}^{n+1} - T_{i}^{n}}{\Delta t} - \frac{\Delta t}{2}\bigg(\frac{\partial^2 T}{\partial t^2}\bigg)_i - \frac{\Delta t^2}{6}\bigg(\frac{\partial^3 T}{\partial t^3}\bigg)_i - \dots \approx \frac{T_{i}^{n+1} - T_{i}^{n}}{\Delta t} ,\]

for equally-spaced timesteps \(\Delta t\), and leads to a truncation error \(O(\Delta t)\). The one-dimensional space-centred integration is

\[\begin{aligned} \bigg(\frac{\partial T}{\partial x}\bigg)_n &= \frac{T_{i+1}^{n} - T_{i-1}^{n}}{\Delta x} - \frac{\Delta x^2}{6}\bigg(\frac{\partial^3 T}{\partial x^3}\bigg)_i - ... \approx \frac{T_{i+1}^{n} - T_{i-1}^{n}}{\Delta x}\end{aligned}\]

for equally spaced grid lengths \(\Delta x\), and leads to a truncation error of \(O(\Delta x^2)\). The second derivative can be calculated as follows:

\[\begin{split}\begin{aligned} \bigg(\frac{\partial^2 T}{\partial x^2}\bigg)_n =& \bigg[ \frac{\partial}{\partial x} \frac{\partial T}{\partial x}\bigg]_n \nonumber \\ =& \lim\limits_{\Delta x \to 0} \left(\frac{\text{forward difference - backwards difference}}{\Delta x}\right) \nonumber \\ \approx& \frac{\frac{T_{i+1}^{n} - T_{i}^n}{\Delta x} - \frac{T_{i}^{n} - T_{i-1}^{n}}{\Delta x} }{\Delta x} \nonumber \\ =& \frac{T_{i+1}^{n} - 2T_{i}^n + T_{i-1}^{n}}{\Delta x^2}\end{aligned}\end{split}\]

Inserting these numerical solutions into Equation (145), the one-dimensional heat diffusion equation can be expressed via a finite-difference scheme as

\[\frac{T_{i}^{n+1} - T_{i}^n}{\Delta t} - \alpha \left(\frac{T_{i+1}^{n} - 2T_{i}^n + T_{i-1}^{n}}{\Delta x^2}\right) = \frac{\dot{q}}{C}.\]

Rearranging for \(T_{i}^{n+1}\) gives

\[T_{i}^{n+1} = T_{i}^n + \Delta t \bigg[ \alpha \left(\frac{T_{i+1}^{n} - 2T_{i}^n + T_{i-1}^n}{\Delta x^2}\right) + \frac{\dot{q}}{C} \bigg],\]

which is also known as the one-dimensional explicit finite-difference solution to Fourier’s law of heat conduction. This scheme is illustrated in Figure (8): it is explicit as the temperature at time \(n+1\) explicitly depends on the temperature at time \(n\), and is evidently forward-in-time and space-centred.

Fig. 8 One-dimensional finite-difference schematic (boundary nodes indicated by dark vertical lines)

The timestep and lattice spacing, \(\Delta t\) and \(\Delta x\) respectively, must be chosen carefully to ensure the stability of this algorithm, and this is provided by defining the Fourier mesh number, \(F\), as

\[F = \alpha \frac{\Delta t}{\Delta x^2}\]

This can be thought of as the ratio of timestep to the time it takes to equilibrate a region of length \(\Delta x\). In this one-dimensional case, the value of \(F\) must satisfy \(0 < F < \frac{1}{2}\), or else the algorithm becomes unstable and oscillates wildly.

In three-dimensions, if \(\Delta x = \Delta y = \Delta z\), the finite-difference solution becomes

(146)\[\begin{split}\begin{aligned} T_{i,j,k}^{n+1} =& T_{i,j,k}^n + \Delta t \bigg[ \alpha \left(\frac{T_{i+1,j,k}^n + T_{i-1,j,k}^{n} + T_{i,j+1,k}^n + T_{i,j-1,k}^{n} + T_{i,j,k+1}^n + T_{i-1,j,k-1}^{n} - 6T_{i,j,k}^n}{\Delta x^2}\right) + \frac{\dot{q}}{C} \bigg] \nonumber \\ =& T_{i,j,k}^n + F [T_{i+1,j,k}^n + T_{i-1,j,k}^{n} + T_{i,j+1,k}^n + T_{i,j-1,k}^{n} + T_{i,j,k+1}^n + T_{i-1,j,k-1}^{n} - 6T_{i,j,k}^n] + \Delta t \frac{\dot{q}}{C} \end{aligned}\end{split}\]

with a new stability criteria of \(0 < F < \frac{1}{6}\). Thus, the size of the timestep must satisfy \(\Delta t < \frac{\Delta x^2}{6 \alpha}\). Equation (146) applies under the assumption that the thermal diffusivity \(\alpha\) is a constant value, i.e. \(\nabla . [\alpha \nabla T ] = \alpha \nabla^2 T\), but the more general (and hence more complicated) case, where \(\alpha\) can vary spatially, takes the form

(147)\[\begin{split}\begin{aligned} T_{i,j,k}^{n+1} =& \frac{\Delta t}{\Delta x^2} \left( \frac{\kappa \big[ \frac{1}{2} (T^n_{i+1,j,k} + T^n_{i,j,k}) \big]}{C(T^n_{i,j,k})} (T^n_{i+1,j,k} - T^n_{i,j,k}) + \frac{\kappa \big[ \frac{1}{2} (T^n_{i-1,j,k} + T^n_{i,j,k}) \big]}{C(T^n_{i,j,k})} (T^n_{i-1,j,k} - T^n_{i,j,k}) \right) \nonumber \\ +& \frac{\Delta t}{\Delta y^2} \left( \frac{\kappa \big[ \frac{1}{2} (T^n_{i,j+1,k} + T^n_{i,j,k}) \big]}{C(T^n_{i,j,k})} (T^n_{i,j+1,k} - T^n_{i,j,k}) + \frac{\kappa \big[ \frac{1}{2} (T^n_{i,j-1,k} + T^n_{i,j,k}) \big]}{C(T^n_{i,j,k})} (T^n_{i,j-1,k} - T^n_{i,j,k}) \right) \nonumber \\ +& \frac{\Delta t}{\Delta z^2} \left( \frac{\kappa \big[ \frac{1}{2} (T^n_{i,j,k+1} + T^n_{i,j,k}) \big]}{C(T^n_{i,j,k})} (T^n_{i,j,k+1} - T^n_{i,j,k}) + \frac{\kappa \big[ \frac{1}{2} (T^n_{i,j,k-1} + T^n_{i,j,k}) \big]}{C(T^n_{i,j,k})} (T^n_{i,j,k-1} - T^n_{i,j,k}) \right) \nonumber \\ +& \Delta t \frac{\dot{q^n_{i,j,k}}}{C^n_{i,j,k}} \end{aligned}\end{split}\]

Here the electronic thermal conductivity has an explicit spatial dependence. To simplify this relationship, \(\kappa\) can be assumed to be constant locally, and is taken to be the average value between the current and neighbouring cells. An adaptive timestep is also utilised, so at each timestep a fraction of the ’worst case scenario’ for the Fourier mesh number, \(F\), is chosen, ensuring the stability of the electronic subsystem.

Various boundary condition choices are available for the edge cells in Figure (8), which surround the simulation in all three dimensions. These are:

• Dirichlet Boundary Conditions: Also known as infinite-flux boundary conditions, the edge cell is fixed at a finite temperature, \(T = T_0\), where \(T_0\) is the target (system) emperature. Dirichlet BCs are usually chosen for cascade simulations.

• Neumann Boundary Conditions: Also known as zero-flux boundary conditions, the temperature of the edge cell is taken to be the value of its corresponding neighbour, thus \(\frac{dT}{dx} = 0\) in this region.

• Robin Boundary Conditions: Also known as partial or variable flux boundary conditions, the temperature of the edge cell is taken to be a fixed proportion of the neighbouring cell’s temperature. Thus \(\frac{dT}{dx} = -k (T-T_0)\), where \(k\) is the fraction of the neighbouring temperature that is ’targeted’.

The electronic energy contained in a voxel representing an electronic temperature grid point can be calculated by integrating the volumetric heat capacity between a datum temperature (e.g. system temperature) and the local electronic temperature, i.e.

\[E_{e}^{j} = \int_{T_{0}}^{T_{e}^{j}} C_{e} dT.\]

A quantity of electronic energy can be added to a voxel by setting the local electronic temperature to a new value, such that the integral of volumetric heat capacity between the original and new temperatures gives the required energy.

The atomic temperatures \(T_a\) and \(T_a^{\prime}\) can be calculated from kinetic energies of atoms in each voxel, i.e. for cell \(j\):

\[\begin{split}\begin{aligned} T_{a}^{j} &= \frac{\sum_{p \in \vec{x}^j} m_p v_p^2}{3 k_B N} \\ T_{a}^{\prime j} &= \frac{\sum_{p \in \vec{x}^j, |v_p | > v_{cut}} m_p v_p^2}{3 k_B N^{\prime}}\end{aligned}\end{split}\]

where \(v_{cut}\) is the cut-off velocity beyond which electronic stopping is significant, \(N\) is the total number of atoms in the voxel and \(N^{\prime}\) is the number of atoms in the voxel with velocities greater than \(v_{cut}\). To account for centre-of-mass drift, the atomic velocities \(\underline{v}_p\) should be corrected by the voxel’s velocity, calculated by \(\underline{v}^j_{com} = \frac{\sum_{p \in \vec{x}^j} m_p \underline{v}_p}{\sum_{p \in \vec{x}^j} m_p}\).

If there are too few (or no) atoms in a voxel, it is considered to be inactive as no definable ionic temperatures can be calculated. Equation (147) does not have to be applied to inactive voxels (setting electronic temperatures to the datum value \(T_0\) and source terms to zero), while temperature gradient terms involving inactive voxels in the same equation can be omitted for their neighbours.

MD system

The principal idea is to modify the MD equations of motion according to Langevin dynamics, which describes the movement of particles in a viscous medium. This viscous medium will represent the electronic subsystem, and the modified equation of motion takes the form

(148)\[{d \underline{v}_{p}(t) \over d t} = {\underline{{f}_{p}(t)+\underline{R}_{p}(t)} \over m_{p}} - \chi_{p} \; \underline{v}_{p}(t),\]

where \(m_p\) and \(\vec{v}_p\) are the mass and velocity of atom \(p\) at time \(t\), \(\underline{f}_{p}\) is the deterministic force on atom \(p\) due to the interatomic potential, \(\underline{R}_{p}(t)\) is a stochastic force with random magnitude and orientation and \(m_{p} \chi_{p} \underline{v}_{p}(t)\) is the frictional force due to the electrons. These last two terms in Equation (148) are the Langevin modifications to Newton’s second law, which allow energy to be lost and gained by the MD system.

The stochastic force \(\underline{R}_{p}(t)\) returns energy from the electrons to the ions and is formulated as \(\underline{R}_{p}(t) = \sqrt{\Gamma} \vec{A}_{p}(t)\), where \(\vec{A}_{p}(t)\) is a three-dimensional vector with components randomly distributed in \([-1,1]\), and \(\Gamma\) is the stochastic friction coefficient. \(\underline{R}_{p}(t)\) must satisfy two important time-averaged conditions:

(149)\[\begin{split}\begin{aligned} \langle \underline{R}_{p}(t) \rangle =& 0, \\ \langle \underline{R}_{p}(t) \cdot \underline{R}_{q}(t^{\prime}) \rangle &\propto& \delta_{pq} \delta (t-t^{\prime}) \end{aligned}\end{split}\]

The first condition states that over a significant period of time, \(\underline{R}_{p}(t)\) must not behave as a net source or sink. Equation (149) is known as the fluctuation-dissipation theorem, which describes how the drag felt by a particle as it moves through a viscous medium can give rise to Brownian motion. In the standard homogeneous Langevin thermostat, every atom in the MD simulation is thermostatted to a target temperature. The inhomogeneous case allows for each atom to be thermostatted to the electronic temperature of the corresponding continuum electronic cell. This leads to the stochastic friction term

\[\Gamma = \frac{6 m_{p} \chi_{ep}^j k_B T_e^j}{\Delta t},\]

where \(\chi_{ep}^j\) is the electron-phonon friction of the \(j^{th}\) electronic finite-element cell, \(T_e^j\) is the electronic temperature of the corresponding cell, \(k_B\) is the Boltzmann constant, and \(\Delta t\) is the timestep. The electron-phonon friction term is thus calculated at each point in the finite electronic temperature grid:

(150)\[\chi_{ep} = \frac{G_{ep} \Delta V}{3 k_B N},\]

where \(\Delta V\) is the volume of the cell (\(\Delta x \Delta y \Delta z\)), \(G_{ep}\) is the electron-phonon coupling constant of the material, and \(N\) is the number of atoms in the cell.

The friction term in Equation (148) is made up of two forms of energy loss: the previous discussed electron-phonon friction, and electronic stopping, which is inelastic electron scattering of ballistic atoms. The total electron friction coefficient \(\chi_p\) for atom \(p\) is given by

\[\begin{split}\chi_{p} = \left\{ \begin{array} {l@{\quad:\quad}l} \chi_{ep} + \chi_{es} & | v_p | > v_{cut} \\ \chi_{ep} & | v_p | \le v_{cut} \end{array} \right. ,\end{split}\]

where \(\chi_{es}\) is the electronic stopping friction, \(\underline{v}_{p}\) is the velocity of atom \(p\), and \(v_{cut}\) is the cut-off velocity for which electronic stopping becomes significant. The electronic stopping friction term can be calculated in a similar fashion to the electron-phonon term:

(151)\[\chi_{es} = \frac{G_{s} \Delta V}{3 k_B N^{\prime}}\]

where \(N^{\prime}\) is the number of atoms in the cell with velocities greater than \(v_{cut}\). Note that this term is set to zero when \(N^{\prime} = 0\).

From Equations (148), (150) and (151), the differences between the contributions from electron-phonon coupling and electronic stopping are evident. Electron-phonon coupling allows energy to flow to and from the lattice (depending on the temperature gradient between ions and electrons), whereas electronic stopping acts solely as an energy loss mechanism for the lattice.

Fig. 9 Schematic of thermodynamic coupling and processes in 2T-MD model

Figure (9) illustrates these processes for swift heavy ion simulations, and highlights how the MD cell is now thermostatted to a heat bath. The lattice will reach local equilibrium with the electron gas, which is thermostatted to the heat bath, thus eventually driving both subsystems to the chosen ambient temperature. Energy can only be removed from the system via the electron gas; this is justified due to how slow lattice heat diffusion is in comparison to electronic heat diffusion.

It is possible to use the inhomogeneous Langevin thermostat (Equation (148)) on its own without coupling it to the electronic temperature grid, but still enhancing the total Langevin friction term for atoms with velocities greater than a cut-off value:cite:zarkadoula-13a. In this case, the stochastic friction coefficient \(\Gamma\) is modified to use the system temperature \(T_0\) instead of a local electronic temperature and to take advantage of the enhanced friction coefficient when electronic stopping applies, i.e.

\[\Gamma = \frac{6 m_{p} \chi_{p} k_B T_0}{\Delta t}.\]

Simulation setup

There are three distinct types of irradiation that can be simulated using the TTM (2T-MD) implementation in : swift heavy ions, laser excitation, and high-energy cascades. These are conducted by splitting the MD cell into discrete coarse-grained lattice ionic temperature (CIT) voxels, and discretising the electronic system into coarse-grained electronic temperature (CET) voxels. Energy can thus be exchanged between the voxels and subsequently passed to or from the atoms within each respective CIT. The volume of each CIT voxel must contain a sufficient number of atoms so that thermal fluctuations of ions are negligible and an ionic temperature can be defined: a good general choice is a cube of length 10 Å in each direction. There is more flexibility in choosing the number of CET voxels, as long as an integer number of these overlap with the CIT grid: to simplify the connections between the CET and CIT grids, equal-sized voxels for both systems will be assumed from here on.

Cascades

Fig. 10 Schematic of cascade simulation setup

High-energy cascades require no initial energy deposition into the electronic system (i.e. \(\frac{dE}{dx} = 0\)): instead, an ion is initialised with a very high velocity. The electronic temperature (CET) voxels extend further than the ionic temperature (CIT) voxels in all directions, with open (Dirichlet) or semi-open (Robin) boundary conditions in all dimensions to represent thermal electronic conduction into the bulk, allowing the electronic temperature to converge towards the initial system temperature \(T_0\). (Figure (10) gives an example schematic of this simulation setup.) Stochastic boundary conditions can be applied in the ionic system to dampen the shock wave generated by the displacement spike.

Swift heavy ions

Fig. 11 Simulation setup for swift heavy ion impact.

Swift heavy ion systems can be modelled using an initial Gaussian spatial energy deposition into the electronic system (i.e. \(\frac{dE}{dx} > 0\)) with either Gaussian or exponentially decaying temporal distribution in electronic temperature. The size of the electronic temperature (CET) grid in the z-direction is set equal to the size of the ionic temperature (CIT) grid in the same dimension, while the CET voxels are extended over the corresponding CIT voxels in the x- and y-directions. Boundary conditions can be set with no energy flux in the z-direction and open or semi-open boundary conditions in x- and y-directions. (Figure (11) gives a schematic of this simulation setup.) Stochastic boundary conditions can be applied to the lattice system in lateral directions only to represent non-negligible phononic thermal conductivity in semiconductors into the builk. Similarly, while electronic thermal conduction in the lateral directions is allowed, conduction parallel to impact is not. This reflects the fact that the simulation represents a small cross-section of the evolution of a micron-sized track.

Laser excitation

Fig. 12 Simulation setup for laser irradiation.

Laser excitation systems can be modelled with an initial homogeneous spatial energy deposition into the electronic system (either in all three directions or in x- and y-directions with exponential decay in the z-direction) with either Gaussian or exponentially decaying temporal distribution in electronic temperature. (The energy deposition can be specified for the fully homogeneous case either by setting \(\frac{dE}{dx} > 0\) or by giving values for the absorbed fluence and penetration depth from the laser.) The size of the electronic temperature (CET) grid is set to the same size as the ionic temperature (CIT) grid, with zero-flux (Neumann) boundary conditions in all directions. This setup (shown in Figure Fig. 12) represents a homogeneous laser excitation with the simulated part as a small section of a larger photoexcited sample.

It is possible in such simulations for voxels to become empty due to displacement of atoms from the laser source. These voxels are omitted from electronic heat diffusion calculations, setting their electronic temperatures to the background value \(T_0\) and their source terms to zero. Their associated spatial gradients in Equation (147) are also omitted for neighbouring voxels.

Implementation

TTM with MD (2T-MD) has been implemented in DL_POLY_5 to take advantage of the domain decomposition scheme used by the code, by dividing up the coarse-grained ionic (CIT) and electronic (CET) temperature voxels as evenly as possible among the processors based on location. This avoids the need for each processor to hold copies of the entire CIT and CET grids and provides good to excellent parallel scalability for larger scale problems.

Coarse-grained ionic temperature (CIT) voxels are divided among processors with overlapping voxels between two or more processors assigned to the first processor in each direction. A boundary halo of voxels is also included to allow communication of contributions to voxel momenta, kinetic energies and atom counters between processors for calculations of ionic temperatures. Since ionic temperatures are only needed for finite-difference calculations of Equation (144), some of these communications only need to be applied in one direction for each dimension.

The coarse-grained electronic temperature (CET) grid is considered as integer multiples of the ionic temperature grid, with equal numbers in both directions of each dimension. While this may provide more CET voxels than requested by the user, the application of boundary conditions in the correct places means that the finite-difference calculations can be carried out in superfluous voxels without affecting the result. The centre of the CET grid is located at the same place as the CIT grid, matching the two up precisely: the electron-phonon, electronic stopping and energy deposition source terms are only applied in these CET voxels. Communications of electronic temperature can be carried out both within each ‘copy’ of the ionic temperature grid and between them: these need to be applied for each iteration (timestep) of the finite-difference solver.

Communications to and from boundary halos for both CIT and CET grids make use of MPI derived data types, which allow for single MPI send and receive calls for grid values without needing to pack and unpack data. This is the same communication technique used in DL_MESO for its lattice Boltzmann equation code:cite:seaton-13a and has been shown to give near-perfect parallel scaling to thousands of processors.

Functionality and directives

The CONTROL file direction ttm_calculate can be used to switch on the two-temperature model (2T-MD) as described above. If no other information is provided, DL_POLY_5 will use default values for certain required properties, but some information must be provided: if this information is unavailable, DL_POLY_5 will terminate. The list of TTM directives is given as part of Section The CONTROL File Directives: more details about these directives are given below.

The inhomogeneous Langevin thermostat can be activated using the directive ensemble_method with the ttm option in the CONTROL file, specifying the electron-phonon friction term (\(\chi_{ep}\), in ps\(^{-1}\) using ttm_e-phonon_friction), electronic stopping friction term (\(\chi_{es}\), in ps\(^{-1}\) using ttm_e-stopping_friction) and the cutoff atomic velocity for electronic stopping (\(v_{cut}\), in Å ps\(^{-1}\) using ttm_e-stopping_velocity). This thermostat is required for 2T-MD calculations but can also be used independently: these CONTROL file directives therefore do not automatically switch on the two-temperature model, but DL_POLY_5 will quit with an error message if TTM is selected and the ensemble is not NVT or the ensemble method is not ttm.

By default, the inhomogeneous Langevin thermostat will be applied only to particle thermal velocities, i.e. velocities that have been corrected to remove the centre-of-mass flow calculated from its coarse-grained ionic temperature (CIT) voxel. The directive ttm_com_correction directive can be set to full for this default behaviour, zdir to only correct the z-direction velocity component and use total velocity components in the x- and y-directions, or off to omit all velocity corrections and use total particle velocities. A warning message will also be printed if the fixed_com option is switched on, as removal of total centre-of-mass motion may affect the dynamics of systems with electronic stopping effects.

The number of coarse-grained electronic temperature (CET) voxels is specified in the CONTROL file with the directive ttm_num_elec_cellst: by default, a grid of \(50 \times 50 \times 50\) will be used if this information is not supplied by the user. The number of coarse-grained ionic temperature (CIT) voxels in the z-direction is specified in CONTROL using the directive ttm_num_ion_cells: the default number is 10, and the number of voxels in x- and y-directions will be determined automatically based on system size. The number of CET voxels must be at least the same as CIT or larger: too few CET voxels will cause DL_POLY_5 to terminate.

The volumetric heat capacity \(C_e\) can be obtained for TTM calculations in four different forms with their corresponding options for the CONTROL file directive ttm_heat_cap_model:

• A temperature-independent constant value (constant)

  \[C_e = C_0 \rho\]

• A linear function of temperature up to a maximum at the Fermi temperature \(T_{f}\) (linear)

  \[C_e (T) = C_0 \rho \max \left(T/T_{f}, 1\right)\]

• A hyperbolic tangent function of temperature (tanh)

  \[C_e (T) = A \rho \tanh \left(10^{-4} B T\right)\]

• A tabulated function of temperature supplied in a Ce.dat file (tabulated).

Note that the values given for the first three options are specific heat capacities (\(C_0\) and \(A\) given in \(k_B\) per atom by the ttm_heat_cap directive), which are converted to volumetric heat capacities by multiplying by the atomic density \(\rho = \frac{N}{\Delta V}\). The atomic density is assumed to be constant throughout the system, and this value can be set in one of three ways: (1) a value can be calculated from the provided configuration at the start (assuming all ionic temperature cells are active), (2) the user can specifiy a value using the ttm_dens directive in the CONTROL file and apply it using the constant option for the ttm_dens_model directive, or (3) after energy deposition to the electronic temperature grid, the value can be calculated dynamically from active CITs when using the dynamic option for the ttm_dens_model directive. In the latter case, the system or user-specified atomic density is used during energy deposition. Tabulated volumetric heat capacities are given in the Ce.dat file as J m\(^{-3}\) K:math:^{-1}, which are converted to \(k_B\) Å:math:^{-3}. If no heat capacity information is supplied, DL_POLY_5 will assume a constant volumetric heat capacity of 1 \(k_B\) per atom by default. In all cases, the electronic energy of a given voxel can be determined from the product of cell volume and the integral of the volumetric heat capacity between the system temperature \(T_0\) and its current electronic temperature \(T_e\).

If the system is metallic (specified by the CONTROL directive ttm_metal as on), a thermal conductivity needs to be supplied: no default value is provided by . Four options are available with their corresponding keywords for the CONTROL file directive ttm_elec_cond_model:

• An infinitely large value (infinite)

  \[K_e = \infty\]

• A temperature-independent constant value (constant)

  \[K_e = K_0\]

• A Drude-model (linear) function of temperature (drude)

  \[K_e (T) = K_0 \frac{T}{T_0}\]

• A tabulated function of temperature supplied in a Ke.dat file (tabulated).

All values (constants or tabulated) are supplied in W m\(^{-1}\) K:math:^{-1}. In the case of infinitely large conductivity, all heat diffusion in the electronic subsystem is instantaneous and all active CET voxels without source terms will be at the same mean electronic temperature. The Drude model uses the electronic temperature for a given voxel, while the tabulated function will use either the ionic temperature for overlapping cells or the system temperature for CET voxels beyond the CIT system.

If the system is non-metallic (specified by ttm_metal as off in CONTROL), a thermal diffusivity needs to be supplied: no default value is provided by . Three options are available with their corresponding options for the CONTROL file directive ttm_diff_model:

• A temperature-independent constant value (constant)

  \[D_e = D_0\]

• A reciprocal function of temperature up to the Fermi temperature \(T_{f}\) (recip or reciprocal)

  \[D_e (T) = D_0 \frac{T_0}{\min \left(T, T_{f}\right)}\]

• A tabulated function of temperature supplied in a De.dat file (ttm detab).

All values (constant and tabulated) are supplied in m\(^2\) s:math:^{-1} and subsequently converted to Å\(^2\) ps:math:^{-1}.

The electron-phonon coupling friction term \(\chi_{ep}\) can either be held at the constant value given in ensemble nvt ttm, or it can be dynamically varied according to electronic temperature. The CONTROL file directive ttm_variable_ep can be used to specify that the electron-phonon coupling constant \(G_{ep}\) is supplied in tabulated form from a g.dat file, with values of \(G_{ep}\) given in W m\(^{-3}\) K:math:^{-1} and converted to \(\chi_{ep}\) values (using Equation ([eq:chiep]) with the mean value for the number of atoms per voxel, \(N = \rho \Delta V\)) in ps\(^{-1}\). Two variants of dynamic coupling are available: homogeneous coupling uses the mean electronic temperature to calculate a system-wide value of \(chi_{ep}\), while heterogeneous coupling uses local values of electronic temperature (based on values in CET voxels) to calculate \(chi_{ep}\) for each atom.

Boundary conditions to the electronic temperature system are applied using the ttm_boundary_condition directive in the CONTROL file:

• dirichlet

• neumann

• robin

• periodic

with the ttm_boundary_xy switch available to apply the Dirichlet or Robin boundary conditions in x- and y-directions with Neumann boundary conditions in the z-direction. These boundary conditions have no direct connection to any boundary conditions used for the MD system: conditions for the latter should be chosen carefully by the user to give the desired effects. In the case of Robin boundary conditions, the ttm_boundary_heat_flux directive can be used to indicate the proportion of energy flux (given as a percentage) leaving the electronic system at the boundaries.

An energy deposition (\(A(r,t) = f(r,z) g(t) \Delta V\)) can be applied using the CONTROL file directive ttm_spatial_dist, which can use the following types:

• gaussian: specifies a Gaussian distribution in x- and y-directions from the centre of the system (applied homogeneously in the z-direction), i.e.

  \[f(r) = \frac{dE}{dx} \frac{1}{2 \pi \sigma^2} e^{-\frac{r^2}{2\sigma^2}}\]

  with standard deviation \(\sigma\) (ttm_spatial_sigma), extending up to a cut-off distance given in multiples of \(\sigma\) (ttm_spatial_cutoff);

• flat: specifies a homogeneous distribution in x-, y- and z-directions

  \[f(r) = \frac{dE}{dV} = \frac{dE}{dx} \frac{1}{L_{x} L_{y}}\]

  where \(L_{x}\) and \(L_{y}\) are the MD system dimensions in x- and y-directions;

• laser: specifies a homogeneous distribution in x- and y-directions with either a homogeneous distribution or exponential decay in the z-direction (from the centre of the system), originating from a laser.

In the case of the Gaussian or flat distributions, the stopping power for the depositions are given as \(\frac{dE}{dx}\) (ttm_stopping_power) in eV/nm. For lasers this quantity is obtained from the fluence \(F\) (ttm_fluence) in mJ cm\(^2\) and penetration depth \(l_p\) (ttm_penetration_depth) in nm. The directive ttm_laser_type enables the type of laser (flat, exponential) to be selected: the flat type is similar to the equivalent non-laser option with \(\frac{F}{l_p} = \frac{dE}{dV}\), while the exponential function

\[f(r,z) = \frac{F}{l_p} \exp\left(-\frac{|z|}{l_p}\right)\]

uses the penetration depth \(l_p\) as the exponential scaling factor. The depositions can be applied temporally in one of four ways with the ttm_temporal_dist option:

• delta: specifies a Dirac delta function in time

  \[g(t) = \delta (t - t_0)\]

  which is approximated as an energy injection during a single diffusion timestep;

• square: specifies a square pulse function in time

  \[\begin{split}g(t) = \left\{\begin{array}{lr} \frac{1}{\tau}, & t - t_0 < \tau\\ 0, & t<t_0 \text{ and } t - t_0 \ge \tau \end{array} \right.\end{split}\]

  over a period \(\tau\);

• gaussian: specifies a Gaussian function in time

  \[g(t) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\left(\frac{t-t_0}{\sigma}\right)^2}\]

  over a period of a multiple of standard deviations \(\sigma\);

• exponential: specifies a decaying exponential function in time

  \[g(t) = e^{-\frac{t-t_0}{\tau}}\]

  over a period of a multiple of \(\tau\).

with the ttm_temporal_duration giving the time duration of the pulse in ps and ttm_temporal_cutoff gives the deposition duration as multiples of the duration. The exponential function is scaled to ensure the correct total energy (within a tolerance of \(\pm 1\)%) is deposited to the electronic system. Energy depositions are achieved by determining the required increases in electronic temperature to add the required energy to the CET voxel.

To calculate ionic temperatures in CIT voxels, a minimum number of atoms is required: this can be specified by the CONTROL file directive ttm_min_atoms, although a default value of 1 is used if this is not supplied by the user. Any CIT voxel that includes this number of atoms will be considered as active and an ionic temperature can be calculated from the mean kinetic energy of atoms in the voxel (after centre-of-mass motion is removed). If the voxel does not contain enough atoms, it is considered inactive: no electron-phonon coupling is applied to the corresponding electronic temperature cell, although thermal diffusion is still applied as normal. All voxels within the CIT grid are checked at each MD timestep and can be reactivated once the minimum number of atoms can be found. If the dynamic option for the ttm_dens_model is also specified, the average cell density (number of atoms per unit volume) will be recalculated after energy deposition by only considering active CIT voxels: note that the system-based or user-specified atomic density (the latter given by ttm_dens) will be used before and during energy deposition.

If the ttm_redistribute option is in use, any recently deactivated voxels have their electronic energy transferred to neighbouring active voxels by increasing their electronic temperatures: these are used for localised Dirichlet (constant temperature) boundary conditions during thermal diffusion in the single MD timestep. The corresponding electronic temperature cell is also deactivated: its electronic temperature is set to the system temperature and the voxel is excluded from thermal diffusion calculations by setting the electronic temperature gradient between itself and any neighbours to zero. This option requires at least one electronic temperature cell beyond the ionic system in each direction to ensure electronic energy from deactivated cells is not lost: if this condition does not hold, the energy transfer option is switched off.

To equilibrate the electronic system prior to any energy deposition, the CONTROL file directive ttm_time_offset can be used to allow the electronic and ionic systems to settle: this switches off the electron-phonon terms in both electronic heat diffusion calculations and the inhomogeneous Langevin thermostat, only applying electronic stopping effects. The ttm_oneway directive restricts electron-phonon coupling to coarse-grained ionic temperature (CIT) voxels and their atoms when the electronic temperature is higher than the ionic temperature: both the electron-phonon coupling term in the thermal diffusion equation and the inhomogeneous Langevin thermostat are switched off if this condition is not met.

To overcome the limitations on the treatment of the high-energy particles with electrons, the CONTROL file directive ttm_e-phonon_cutoff_velocity can be used to address spurious double-interactions: this switches off the electron-phonon terms in fast-moving particles that exceed a threshold velocity, only applying electronic stopping effects (removing the need to specify a thermalization time). Additionally, particles subject to electronic stopping are excluded from the CIT voxel to prevent inaccurate temperature elevation ¹⁰⁰.

The CONTROL file directives ttm_stats_frequency and ttm_traj_frequency switches on outputs of statistical information and limited temperature trajectory snapshots of the two-temperature model simulation. Ionic temperatures (minimum, maximum, mean and sum) are supplied at user-defined intervals in the file PEAK_I, while electronic temperatures and energies are supplied in the PEAK_E file. Ionic and electronic temperature profiles along the y-direction at the centre of the xz plane are given in the LATS_I and LATS_E respectively at user-defined intervals.