A modified two temperature molecular dynamics (2T-MD) model for cascades

Two-Temperature molecular dynamics (2T-MD) is a common approach for describing how electrons contribute to the evolution of a damage cascade by addressing their role in the redistribution of energy in the system. However, inaccuracies in 2T-MD’s treatment of the high-energy particles have limited its utilisation. Here, we propose a reformulation of the traditional 2T-MD scheme to overcome this limitation by addressing the spurious double-interaction of high-energy atoms with electrons. We conduct a series of radiation damage cascades for 30, 50, and 100 keV primary knock-on atoms in increasingly large cubic W cells. In the simulations, we employ our modified 2T-MD scheme along with other treatments of electron–phonon coupling to explore their impact on the cascade evolution and the number of remnant defects. The results suggest that with the proposed modification, 2T-MD simulations account for the temperature time evolution during the ballistic phase and remove arbitrary choices, thus providing a better description of the underlying physics of the damage process.


Introduction
Exposure to radiation dramatically alters the physical properties of materials.Such alterations are often detrimental to the performance of the material in its given use case, for example, those that weaken the integrity of structural materials in nuclear fission and fusion reactors [1][2][3].In the case of nuclear reactors, the most relevant process involves high-energy neutrons impinging on crystalline materials.Neutron-atom interactions ultimately lead to the production of host atoms with a relatively high (non-thermal) kinetic energy.The newly energized host atom (or 'primary knock-on atom' -PKA) then goes on to displace other atoms in the host matrix, and the process repeats in a so-called 'cascade' while the energy transferred to a host atom is sufficient to overcome the threshold displacement energy [4,5].The difference in structure pre-and post-event is, generally, radiation damage and is what ultimately leads to a modification of the host's properties.
Given the short temporal and spatial scales of the radiation damage process, the current experimental technologies are incapable of providing direct observation of cascade events.Due to their ability to describe phenomena occurring on the order of nanometres and picoseconds, atomistic simulations, particularly molecular dynamics (MD), become a suitable technique to predict the atomic dynamics following irradiation damage.However, these simulations employ empirical interatomic potentials that, by construction do not account for the role of electrons, limiting the accuracy of MD results [6][7][8], as the process of energy transfer involves both atomic and electronic degrees of freedom and their interaction.At low projectile velocities, energy is transferred via elastic collisions with the target's nuclei.However, at higher velocities, the cross-section for interaction with target nuclei decreases and electronic stopping dominates, while the energy subsequently transferred to the electronic system may be spatially redistributed before returning to the lattice.This redistribution may act to considerably alter the evolution of a cascade and hence final damage outcomes.Only recently have simulations capable of fully incorporating this full flow of energy become possible [5].
The general problem of including electronic excitation effects in classical MDs simulations has previously been reviewed [9].One modeling approach which aims to treat atomic and electronic degrees of freedom in tandem is twotemperature MDs (2T-MD) [10], in which, a traditional MDs supercell is coupled to an underlying electronic continuum mimicking the electronic excitation in the system.In this model, energy is able to move between the ions and the electrons at a rate dictated by the difference in the temperature between the two systems and the material-dependent electronphonon coupling strength.In addition to being applied to the simulation of radiation damage cascades [11], 2T-MD has been applied to processes such as swift heavy ion irradiation [12,13] and ultrafast laser irradiation [14].Ultrafast laser irradiation in particular has enabled the validation of the technique against experiments conducted over the same time and length scales [15].
Tungsten is a candidate material for the first-wall and divertor for future nuclear fusion reactors [16][17][18].Industrial challenges facing such plasma-facing materials relate to neutroninduced radiation damage and its consequences for thermal and structural performance.The distribution of vacancy clusters (voids), and networks of dislocations are direct consequences of many individual neutron-induced radiation damage events and their repeated influence.In this context, experiments are now able to map void size distributions and image dislocation loops in situ [19,20].Furthermore, numerical simulations can also provide insights into the accumulation of damage over repeated radiation damage events (see for example [21][22][23][24][25]).Nonetheless, these simulations often neglect electronic effects entirely or represent electronic effects using a simple friction model.A predictive approach for evaluating the degree of damage accumulated in single events, which considers the dynamics of energy transfer between ions and electrons in a way which is insensitive to particular model assumptions, is therefore desirable, and key to the development of the quantitative approaches aimed at accurately characterizing damage.
Previously, 2T-MD simulations have shown that electronic effects in the radiation damage process in tungsten can be significant [26].Related 2T-MD studies on tungsten have also probed other relevant physical phenomena, including surface erosion [27], ultra-fast (fs) laser-induced melting [28,29], and laser-induced structural phase changes [30].Nonetheless, in traditional 2T-MD, with finite electron-phonon coupling and electronic stopping, the physical mechanisms by which highenergy atoms are stopped are often unphysically compounded.Fast-moving atoms are not only coupled to electrons directly via the stopping term but also via the electron-phonon coupling through the MD temperature estimator.Hence, as electronic stopping may include effects of electron-phonon coupling, there is a double-counting of energy transfer to the electronic system from fast-moving atoms.In the real system, these two mechanisms are not artificially separate, but a continuation of the same sets of guiding physics, i.e. that fast-moving atoms lose energy to other atoms and to electronic excitation.Modern approaches such as that of Tamm et al [31] aim to put these effects on a proper level basis and to allow for a less ad-hoc ab initio treatment of the crossover between electron-phonon coupling and high-energy stopping physics.Although they are non-perturbative and apply out of thermal equilibrium, these approaches require detailed parameterisation using time-dependent DFT, creating a significant barrier to adoption.Therefore, investigations into alternatives for improving the modeling of electronic excitation effects in MDs simulations are warranted.A middle-ground remedy to the problem posed by double-counting in 2T-MD is simply to remove one of the effects which is double-counted.This maintains the temporal and spatial scalability of the approach while removing an accuracy-limiting factor.
In the present work, we propose reformulating the traditional 2T-MD approach to overcome the inaccurate treatment of high-energy particles.In this modification, the ionic temperature entering the electron-phonon coupling will be calculated solely from the atoms close to equilibrium as these atoms are in a thermodynamically well-defined state.Additionally, fastmoving atoms which exceed a threshold velocity will be subject to electronic drag only ensuring equivalence with experimental observation and electronic stopping in other codes, including stopping range of ions in matter (SRIM) [32,33].To test our modified 2T-MD model, a series of cascade simulations are performed in tungsten to study the damage process by employing different schemes for the treatment of electronphonon coupling.
Note that this study is focused on the methodology for the treatment of electron-phonon coupling that addresses the spurious double-interaction of high-energy atoms with electrons.Direct comparison of the results of the present work with experiment is complicated by other factors in the model, including the choice of potential and the parameterisation of the specific heat and thermal conductivity.It would be interesting to examine whether the observed defect production could be related to differences in modelled spectroscopic signals from experiments, as was done in [34], however, this is beyond the scope of the present work.
The remainder of the manuscript is organized as follows.Section 2 describes the 2T-MD methodology, explains our modification to the 2T-MD approach, and provides details about our simulations.Section 3 presents our findings on radiation damage cascades employing different schemes for the treatment of electron-phonon coupling.Finally, section 4 summarizes the results and draws conclusions.

Two-temperature MDs (2T-MD)
Here we briefly outline the modern 2T-MD theory, which is introduced fully in [10] and the subject of various other studies in addition to those already referenced [35][36][37][38][39].In 2T-MD the MDs equation of motion for each atom i is modified by the presence of electronic stopping and electron-phonon coupling via the introduction of frictional and stochastic force, i.e.: where, γ i is the friction coefficient on a given atom, and F is the stochastic force.The value of γ i depends on the velocity of the considered atom, where, γ ep,s are the friction coefficient due to the electronphonon coupling interaction and the electronic stopping, respectively.The parameter γ s is normally taken to be a constant for a given material, in order to match tabulated electronic stopping power data, while γ ep varies depending on the electronic-temperature-dependent electron-phonon coupling parameter.The velocity parameter, v 0 , is a cutoff, beyond which electronic stopping applies in addition to electronphonon coupling.
A value of twice the cohesive energy is often taken to define v 0 [9,40], although there is no strong physical basis for this choice.Linear stopping manifests from velocities, where atom motion is uncorrelated with that of host atoms, up to the electronic Fermi velocity v F [32].As v F is relatively high, this defines an extensive range.Below v 0 , we are supposing that all electronic stopping physics is captured by the (screened) electron-phonon interaction.At slightly higher energies, beyond v 0 , purely electronic linear stopping is considered to match the SRIM table data, whose physical origin is in the screened atom-electron interaction between atoms and (predominantly) valence host electrons [32,33].However, note that the construction of this data may include remnant lowenergy electron-phonon effects, owing to the largely experimental sources of SRIM stopping powers.This intermediate range has been the subject of a recent study on stopping from the point of view of first-principles time-dependent DFT calculations [41].Ultimately, at very high atom velocities (≳ v F ), stopping peaks and begins to recede.There, moving atoms are partially or fully ionized, and stopping from this point on is accurately described by the Bethe-Bloch formula [42,43].In summary, for extremely high-energy PKA events, there ought to be adjustments made to the electronic stopping as currently implemented.
The energy lost from atomic stopping processes enters the electronic subsystem as a spatially dependent source, where, C e (T e ) is the electronic specific heat, κ the electronic thermal conductivity, G ep (T e ) the electron-phonon coupling, G s an effective coupling which balances energy lost from atoms due to electronic stopping, and T e and T a the electronic and ionic temperatures respectively.The parameters (C e , G ep , and κ) depend on the electronic temperature, atomic temperature, structure, and density.Nonetheless, the role of electronic temperature on C e and G ep is predominant for cascade simulations [38,40].Equation ( 3) is solved discretely, with sources and diffusion occurring in and between cells (or voxels) which are commensurate with corresponding cells in the atomic system.Let the volume of these cells be defined as V, while considering a uniform array of cubic cells each having the same volume, with each containing enough atoms to define statistically meaningful values of T a .In order to balance energy on average3 , the parameter T ′ a is required to be equal to the 'average temperature' of all atoms moving with v > v 0 .This has no physical interpretation as a temperature other than as a numerical device for energy balance.Indeed, the fast-moving particles to which this stopping applies may well be out of thermal equilibrium with the remaining atoms in the system.The same energy balance requirement leads to the relations between γ ep,s and G ep,s as, with N ′ the number of atoms in a cell affected by electronic stopping in a given step, and k B the Boltzmann's constant.In fact, G ep (T e ) will be used in our simulations to determine γ ep on a per-step, per-cell basis (i.e. as a function of the local electronic temperature), while a fixed γ s will determine G s .All of the associated energy balance restrictions still apply on a per-cell basis.Energy transfer from the electronic subsystem to the atomic subsystem occurs via the random force, which is white-noise-correlated with a magnitude determined from the fluctuation-dissipation relation, which may be integrated into the dynamics via the use of, for example, a Langevin impulse method [45].
The current 2T-MD model conducts an artificial distinction between 'fast' and 'slow' moving projectiles based on an arbitrary threshold.Further, the existing scheme for the treatment of high-energy particles turns-off electron-phonon coupling entirely after the initiation of a cascade until some 'thermalization time' [46].The rationale for this is that electronphonon coupling is a close-to-equilibrium process that can be applied only in situations where the atomic temperature is well defined.In this thermal offset scheme, during the ballistic phase, the mechanism for energy loss from the lattice directly to the electronic system is limited.Only electronic stopping can facilitate such energy loss, and only for the fastest particles in the simulation.This implies neglecting all of the important short-time dynamics involving the electronic and ionic systems, which partly motivate the use of 2T-MD in the first place.Moreover, by dismissing electron-phonon coupling, it results in an artificially hot ionic system whose particles below the stopping threshold cannot lose energy to the electrons.After the thermalization time, the switching on of electron-phonon coupling acts as an instantaneous quench of lattice temperature.This occurs at a time when elevated lattice temperatures may have only recently been established, forming the thermal spike, and offering the system a means of removing energy, potentially freezing in an artificially enhanced amount of damage.Furthermore, from a theoretical standpoint, the stage at which one defines a system as 'thermalized' holds ambiguity.Therefore, the arbitrary definition of fast and slow particles and the choice of initialisation time for electron-phonon coupling have a significant impact on the cascade evolution and the number and morphology of the remnant defects.

Stopping modifications to 2T-MD
We propose instead that, for each atom, electronic stopping should apply above a cutoff v 0 , but at that point, electronphonon coupling should not, i.e.: This has the consequence, via energy balance requirements (see appendix), that the atomic 'temperature' entering the electronic heat equation is limited to exclude high-energy atoms, which do not then couple directly to electrons by means of an equilibrium-parameterized electronic-temperature-dependent electron-phonon coupling parameter.Therefore, in the evaluation of T a , we exclude particles moving with velocity v i > v 0 .This is in direct analogy with the energy-balancing scheme for electronic stopping.In practice, our modified scheme assumes that electronic stopping is dominant beyond the cutoff, effectively preventing the disproportionate impact of a few shortlived high-energy atoms on heat exchange with the electronic system.Within this framework, it is not necessary to define an initiation time for electron-phonon coupling.As opposed to the alternative of disabling electron-phonon coupling altogether for some time period after the initiation of a damage event [46].Moreover, the resulting dependence on cutoff parameter v 0 is expected to be weak, owing to the large separation of scales between thermal and radiation energy scales (a few meV, to tens or hundreds of keV), and the short amount of time that fast atoms are expected to spend with intermediate velocities on the order of the cutoff.Electron-phonon interactions and electronic stopping do not have a clear separation in reality.In that sense, the cutoff v 0 is just a device to avoid any double counting in electron-atoms interactions.
As an example motivating the removal of athermal effects from the electron-phonon coupling, consider the instantaneous introduction of a PKA having additional kinetic energy ∆E4 into a cell of N atoms.The cell, assuming a naive temperature estimator based on numerical averaging of v 2 (as in traditional 2T-MD), suffers a temperature change associated with the introduction of the PKA (∆T) of, For a cubic cell of W of side 150 unit cells long (6.75 million atoms), a PKA gaining 50 keV in energy results in a temperature jump of around 57.3 K.For a 300 keV PKA, the jump is 6 times as large, around 343.8 K.This effect is amplified in 2T-MD simulations, where the effective cell size implies N eff ≪ N.Moreover, in 2T-MD simulations, it is localized to the ionic cell containing the PKA, resulting in ionic temperatures of thousands of Kelvin.This spurious rise in temperature is the reason the electron-phonon coupling is turned off in the first steps of the cascade.Here, the assumed thermal equilibrium is the problem, and in this case, the treatment of electronphonon coupling as two thermal baths (ionic and electronic) in separate thermal equilibriums is not suitable.In reality, electronic stopping processes and kinetic energy loss from collisions slow fast-moving particles over a short timescale of O(100) fs (the 'ballistic phase').The ionic temperature profile does not simply experience a sudden rise immediately after the inception of a PKA.Instead, it is associated, at first, with the formation of a thermal spike in the region surrounding the PKA.As the energy dissipates more slowly due to stopping collisions, the subsequent thermalization of this region occurs over a timescale of several tens of ps (the 'thermal phase').

2T-MD simulations
All of our 2T-MD simulations were performed using the DL_POLY code [47].We have simulated radiation damage cascades of 30, 50 and 100 keV PKAs in increasingly large cubic W cells, at a target temperature of 300 K. We have chosen these PKAs, as both nuclear and electronic stopping are active in this range.The cells were created by taking 94, 150, and 190 repetitions of the W bcc unit cell in each of the cardinal directions, resulting in simulation boxes of 1661168, 6750000 and 13718000 atoms.To describe the W interactions, we employed the extended Finnis-Sinclair potential of Dai et al [48] joined to the short-range repulsive Ziegler-Biersack-Littmark [49] to screen nuclear repulsion for describing highenergy collisions between atoms.This potential [48] reproduces many of the experimentally observed properties of W, including a lattice constant of 3.160 Å, while also having a relatively simple functional form.We used variable time steps to describe dynamics on short time scales after the introduction of fast particles, such that minimum and maximum move sizes are recognized and time steps altered to ensure dynamics always falls within these limits.The maximum allowed time step was overridden and has a value of 1 fs, while the maximum displacement was 0.1 Å.In 2T-MD, both the ionic and the electronic subsystem are divided into finite element cells, i.e. voxels that form two grids in order to determine coarse-grained temperature estimators, as discussed in section 2.1.In our simulations, ionic grids of 20 × 20 × 20, 33 × 33 × 33, and 41 × 41 × 41 (for PKAs of 30, 50, and 100 keV) in each of the cardinal directions correspond to an average of ca.200 atoms per cell ensuring that the electronic and atomic temperatures are not subject to harsh oscillations.Whereas larger electronic grids of 30 × 30 × 30, 43 × 43 × 43, and 51 × 51 × 51 extending five more grid points than the ionic grids in each of the cardinal directions allow a better representation of the electronic thermal heat spread away from the system.Moreover, the Dirichlet boundary conditions were applied to the outer edges of the electronic grids.The parameters of the simulations are summarized in table 1.With respect to the parameters of the two-temperature model, we have taken electronic stopping friction terms to match SRIM data, with a value of γ s = 1.1 ps −1 .We have chosen to make clear our parametric dependence on T e in the parameters C e (T e ) and G ep (T e ) and assumed no lattice or electronic temperature dependence on κ to guarantee a more consistent energy exchange through electron-phonon coupling in our simulations, in agreement with previous works [26,38,40,50].We have used C e (T e ) and G ep (T e ) as parameterized by DFT calculations in [51], and the value of κ was defined at the target lattice temperature of 174 W m −1 K −1 .
The simulation box was initially equilibrated for 100 ps, at 300 K under constant pressure, and temperature (NPT) conditions using the Berendesen et al [52] barostat and thermostat.Followed by a simulation under constant volume and temperature (NVT) conditions using the Nosé-Hoover [53] thermostat with a relaxation time of 1 fs for 100 ps.Once equilibration was complete, the ionic system was connected to the electronic (continuum) system, and the initial electronic temperature was set to 300 K during 200 ps.
To conduct a cascade event, one atom at the center of the simulation box was given an initial velocity to create a PKA that shocks other atoms and create a cascade collision.For each of the schemes for the treatment of electron-phonon coupling considered in this work, we simulated 12, 16, and 14 cascade events of 30, 50, and 100 keV PKAs, with randomly chosen directions.In all cases, the data we present are ensemble averages over several PKA directions.We have assumed that individual samples from the distribution of all PKA directions are normally distributed in forming means and standard deviations.Furthermore, in this work, we label a cascade as 'ended' when the defect counts are approximately static in time.This is given that long-term defect dynamics acts over a time scale upwards of ns and beyond into years or decades, which are prohibitive for the type and number of simulations required.

Electron-phonon coupling methods.
We considered four schemes for the treatment of electron-phonon coupling, as follows: (i) Cutoff (ceph): electron-phonon coupling evaluated excluding fast-moving particles from the temperature (Ta) estimator, while these fast-moving particles transfer (lose) energy with the electronic system via electronic stopping only, section 2.2.(ii) Thermal offset (teph): existing approach of disabling electron-phonon coupling for a 'thermalization time' [46].Unless otherwise stated, this time fixed at 0.2 ps.(iii) Bare (bare): all atoms (regardless of velocity, or simulation time) enter the temperature estimators and are coupled to the electronic system.(iv) Friction-only (fric): disable electron-phonon coupling entirely at all times (note the electronic subsystem is present to enable evaluation of electronic temperature).
Scheme (i) is the focus of the adjusted method we propose and is outlined in section 2.2.Scheme (ii) is the existing approach described in section 2.1.Scheme (iii) is generally agreed to have a significant drawback, given that the interaction between fast-moving PKAs and secondaries coupled to the electronic system through both stopping processes leads to a dubious energy exchange between electrons and ions.Scheme (iv) removes electron-phonon coupling entirely and the ionic system only exchanges (loses) energy with the electronic system via electronic stopping, effectively removing the two-temperature model from the simulations.This scheme is employed to evaluate the electronic temperature and its behaviour compare to the others.
On the one hand, in scheme (iv), all energy which is lost to electronic stopping (friction) is instantaneously removed and has no impact on the progression of primary damage events.On the other hand, in schemes (i), (ii), and (iii) to a lesser or greater degree, the energy lost to electrons by the ions can be redistributed in the electronic system.The physical mechanism for heat loss to a bulk electronic system is allowed (via an extended grid of electronic voxels in 2T-MD, and the boundary conditions on the outside layer), but before energy is lost, it may re-heat, or support a diminished heat loss from, the ionic system.Whether or not such re-heating is relevant in any given material system depends on the timescale for electron-phonon heating (determined by electron-phonon coupling parameter), the timescales for heat loss to the electronic bulk (electronic diffusivity-a function of heat capacity and electronic thermal conductivity), and the timescale for atom-atom heating (the extent of the ballistic damage phase).
It is worth mentioning, that recently developed approaches such as that of Tamm et al [31,54] which make no distinction between electron-phonon coupling and electronic stopping processes would not require decoupling effects, as they incur no such double-counting and do not rely on the ability to define a lattice temperature.However, the salient point remains: whether or not electronic effects matter will depend on what heat lost to the electronic system does once there.In the few 2T-MD calculations employing this method in the literature, in the LAMMPS code [55], there is no temperature variation in the remaining key model parameters for the electronic system, and so a full description of heat redistribution effects from more sophisticated calculations is still absent.

Temperature evolutions
Figures 1 and 2 show the maximal and average temperature in the ionic and electronic grids as a function of time for the 50 keV PKA.The most remarkable observation from these figures is the dramatically different evolution of the lattice temperatures between schemes during the early part of the cascade.For the teph, bare, and fric schemes, the maximum (of the average and maximal) ionic temperatures are present at the start of the cascade.However, for our modified 2T-MD scheme (ceph), fast-moving atoms which exceed threshold velocity are not included in the estimator (T a ), therefore the temperature increases more slowly as energy is transferred from the PKA.After 0.1 ps, the cutoff scheme reaches its maximum ionic temperatures, and after 0.3 ps, it reaches temperatures at similar scales displayed in teph, bare, and fric schemes.This indicates that the thermal state of the ionic cells is similar for the different schemes during this phase, although the fast-moving ions are dismissed for the kinetic energy calculation in ceph.It is important to understand that the transfer of energy, due to electron-phonon coupling is dictated by the temperature difference between the ionic and electronic systems, therefore, the differences in lattice and electronic temperatures plotted in figures 1 and 2 suggest a significant change in energy exchange and distribution between the different techniques.In agreement with previous authors [26], we find that for tungsten, the electronic grid temperature moves back toward values close to its target before the recombination phase is fully developed.For all the schemes, the ionic temperatures are higher than the electronic temperatures after 0.2 ps.Therefore, in our cascades, the electronic system mainly acts to remove energy from the ions.In other words, the electronic heat diffuses away from the cascade region more quickly than it can be fed back to the ionic cells.While there is an overall trend for the electronic system to remove energy from the ions, the exact magnitude depends on the coupling scheme employed.
From figures 1 and 2, we find that our modified 2T-MD scheme yields the lowest electronic temperatures during the ballistic phase.Notably, the bare scheme displays the highest Maximal value for the ionic and electronic grid temperatures in the system for a 50 keV PKA, averaged over multiple cascades.The continuum and dashed lines represent the coarse-grained ionic and electronic temperatures, respectively.The inset displays the temperatures after 0.2 ps has elapsed, i.e.where electron-phonon coupling is activated on the thermal offset scheme.Average value for the ionic and electronic grid temperatures in the system for a 50 keV PKA, averaged over multiple cascades.The continuum and dashed lines represent the coarse-grained ionic and electronic temperatures, respectively.The shaded regions denote the standard deviation in the spread of data taken over multiple cascades.
average electronic temperatures at the start of the cascade, as there is transfer due to both electronic stopping on fast-moving atoms and electron-phonon coupling on all atoms.Further, we find that the temperatures for the cutoff and bare models are similar beyond ca.0.3 ps.For the thermal offset and frictiononly schemes, we see identical results until 0.2 ps (where the electron-phonon coupling is activated in the former).From the temperatures of the thermal offset scheme, we observe that the effect during the (excessive) heat spike region holds importance during the latter stages of the cascade, as the ionic and electronic temperatures are larger than those from ceph and bare schemes during the latter stages of the cascade, although the electron-phonon coupling is present in all of them.From the temperatures of the friction-only scheme, it is clear that the ionic temperature for the model is significantly different from the other methods at the end of the cascade, given that all particles have velocities below the cutoff threshold, therefore, energy cannot move to the electronic system where it can diffuse away.With no other avenue for energy to leave the simulation cell, the ionic temperatures remain elevated.Moreover, as the energy that is transferred early in the cascade rapidly diffuses away and, there is no mechanism for more energy to move into the electronic system, we find that the electronic temperature returns to the target value most rapidly on the friction-only approach.All other approaches include some form of electron-phonon coupling so as the lattice temperature is greater than the electronic temperature at the latter stages of the cascades energy is transferred to the electronic system, resulting in electronic temperatures (slightly) above the target temperature.

Defect counts & final distributions
During the course of a cascade, we employ an on-thefly defect-counting methodology to classify interstitials and vacancies on a spherical site-based approach.Particles which exist within some cutoff distance (1 Å) from a singly-occupied reference position are not counted, the unoccupied sites are listed as vacancies, and interstitials are counted as multiple occupancies of the reference sphere.In figure 3 we show the defects a typical 50 keV PKA cascade for the different schemes.Here, we observe that schemes that neglect electron-phonon coupling at the early stages of the cascade display a higher number of defects and a wider spread at the maximum defect count compared to schemes that include it.Moreover, we observe a larger creation of clusters (consisting of two or more interstitials) taking place in the fric and teph cascades followed by ceph and bare cascades wherein the sum of clusters is decreased.We will perform a more detailed analysis of clusters in section 3.3.Dynamic defect counts for the different PKAs for all schemes are given in figure 4. From this figure, we see a clear distinction between phases of the cascades in terms of defect creation and healing.Furthermore, in the case of W, the treatment of electronphonon coupling via either the cutoff or bare schemes results in a similar impact on the course of cascades and the final defect counts.However, their impact differs from schemes that neglect electron-phonon coupling at any stage of the cascade (teph, fric), which yield a much higher count of defects.The explanation for the similarity between the results for ceph and bare is related to the relative magnitudes of electronic and lattice thermal conductivities in W (where around 90% of the total thermal conductivity is attributed to electrons [56][57][58]) and in similar metals where electronic thermal conductivity dominates [26].Here, the electronic system acts predominantly as a homogeneous Langevin bath at the system's target temperature, as explained in section 3.1.This goes in line with the average ionic and electronic temperatures of the cutoff and bare schemes beyond 0.3 ps, cf figure 2. On the other hand, concerning the higher defect counts on schemes that neglect electron-phonon coupling during the cascade.We note that, by the initial absence of electron-phonon coupling as a loss mechanism, the time offset and friction-only schemes will lead to a hotter heat spike region at the beginning of the cascade.Notably, in these schemes, particles below the stopping threshold do not lose energy to electrons at all during the ballistic phase, and only after the 'thermalization' time on the time offset scheme, energy will be transferred as a consequence of the lack of thermal equilibrium in the collective ensemble.Physically speaking, these particles should be coupled to electrons at every time.Given that their excessive energy could have a great influence on the physics of the system during and after the time offset has elapsed.
To perform a more detailed analysis of the final configurations after the cascade has ended, the distributions of remnant defects are given in figure 5. First, as expected, we observe a monotonic increase in the number of remnant defects as a function of the PKA energy.From this figure, one can note that the friction-only simulations display more defects at the end of the cascade, followed by the thermal offset scheme, while the cutoff and bare schemes display similar values across all PKA energies.This shows that the atoms decoupled from the electronic system during the thermalization time produce an excess of defects that cannot be recovered during the subsequent mitigation of the damage.Further, the contribution of the electron-phonon coupling shows to be enough to cause the thermal offset and friction-only schemes to differ, despite the electronic system acting (predominantly) as a Langevin heat bath.From figure 5, we observe that for the friction-only and thermal offset schemes, the data present a greater spread depicted by the extent of the T-shaped whiskers.Our results suggest that the arbitrary definition of fast and slow particles, along with the selection of the initiation time for electron-phonon coupling, have a substantial impact on the remaining defects.

Cluster analysis
To gain further insight into the type of remnant defects, we perform a more detailed analysis of the final configurations after  the cascade has ended.In this analysis, we employ the Wigner-Seitz defect analysis to identify the formation of clusters.
Where clusters consist of 2 or more interstitial defects close to each other, identified using the first nearest neighbor distance as implemented on Ovito [59].Table 2 shows the average values for the largest, mean, median, and number of clusters across different 100 keV cascades for all schemes.From these results, we observe that the friction-only cascades display the highest values for the largest and mean cluster size, followed by thermal offset.Further, the number of clusters is similar between schemes.Moreover, we classified the clusters into 5 ranges depending on the number of interstitials within them.The ranges are 2, 2-15, 15-30, 30-45, and more than 45 interstitials in the same cluster.Figure 6 shows (a) the averaged number of atoms (interstitials), and (b) the averaged number of clusters per range for all cascades.From figure 6(a), we observe that (on average) schemes that turn off electron-phonon coupling at the start of the cascade display very large clusters wherein more than half of their interstitials exist in these types of clusters.Whereas for schemes that consider some form of electronphonon coupling during the ballistic phase, the share of interstitials contained within large clusters is significantly lower.Indeed, from figure 6(b), we observe that very large clusters are absent in the bare scheme, and only represent a small fraction on ceph.This is depicted when defect counts are measured via the sphere criterion, cf figures 3(e)-(h).These results indicate that ceph and bare hinder the existence of large clusters at the end of the cascade, whereas they are ubiquitous on teph and fric.This occurs as a consequence of a hotter spike region caused by turning off electron-phonon coupling during the early stages of teph and fric cascades, enabling higher temperatures associated with increased kinetic energies and the formation of larger clusters.

Parametric dependences
We have simulated a separate set of simulations at PKA energies of 50 keV to examine the particular dependence of the time offset and cutoff velocity on the physics of the cascades on the teph and ceph schemes, respectively.We explore cutoff values of 1, 2, and 4 times the cohesive energy for the cutoff cascades, and electron-phonon coupling activation times of 0.1, 0.2, and 0.3 ps for the thermal offset cascades.The distributions of remnant defects of these cascades are given in figure 7. From this figure, we see that the individual parametric variation for a 50 keV cascade in W is lesser, compared with the variation between schemes themselves.Furthermore, we observe that the thermal offset scheme presents a greater spread with a more random distribution, depicted by the extent of the T-shaped whiskers and the means and medians within the boxes.

Conclusions
We have examined the effects of a proposed modification of the 2T-MD method, which removes the spurious doublecounting of atom-electron interactions for relatively fast atoms in radiation damage cascade simulations.This will enable a better description of the radiation damage cascades within the intermediate range for the projectile energy wherein both electronic and nuclear stopping are important.In this modification, fast-moving atoms exceeding the threshold velocity are excluded from the temperature estimator and, consequently, from the energy exchange via electron-phonon coupling.Thus, these atoms only exchange energy with the electrons via electronic stopping.We find the general conclusion that 2T-MD schemes where electron-phonon coupling terms are disabled during the ballistic phase produce markedly higher estimates of primary damage quantities, with a larger number and spread of defects along with larger types of clusters, than schemes that either include the spurious ionic temperature during thermalization or which address its limitation via another means (our modification).We give a physically motivated explanation for why the former approaches may overestimate damage, in terms of the dynamical behaviour of the thermal heat spike region formed at short times after the initiation of a cascade.Our results suggest that with the proposed modification, 2T-MD simulation produces a better and unambiguous description of the underlying physics of a cascade by establishing a more physics-based treatment of the fast-moving ions.This reformulated 2T-MD scheme offers a new approach that accounts for the temperature time evolution during the ballistic phase while mitigating the problem of double-counting.Besides, it avoids the need to select a specific thermalization time, which could have a significant impact on the remaining defects.

Figure 1 .
Figure 1.Maximal value for the ionic and electronic grid temperatures in the system for a 50 keV PKA, averaged over multiple cascades.The continuum and dashed lines represent the coarse-grained ionic and electronic temperatures, respectively.The inset displays the temperatures after 0.2 ps has elapsed, i.e.where electron-phonon coupling is activated on the thermal offset scheme.

Figure 2 .
Figure 2.Average value for the ionic and electronic grid temperatures in the system for a 50 keV PKA, averaged over multiple cascades.The continuum and dashed lines represent the coarse-grained ionic and electronic temperatures, respectively.The shaded regions denote the standard deviation in the spread of data taken over multiple cascades.

Figure 3 .
Figure 3. Representative snapshots of a 50 keV PKA with an initial direction of 0.95559, 0.17764, and 0.34629 in W, for all schemes.(Top) Defects at the maximum extent of defect counts attained at (a) 1.43878, (b) 1.80263, (c) 1.37672, and (d) 1.85213 ps after the initialisation of the PKA.The number of defects is 72688, 156100, 72020, and 166540 for figures (a)-(d) respectively.(Bottom) Final position of the defects, the number of defects is 184, 522, 182, and 430 for figures (e)-(h) respectively.The red sphere and blue squares represent the interstitials and vacancies, respectively.The length of each size of the simulation box represented by the cyan lines is 473.96Å.

Figure 4 .
Figure 4. Dynamic defect counts as per the sphere criterion with a cutoff distance of 1 Å.(a) 30 keV PKA.(b) 50 keV PKA.(c) 100 keVThe orange, blue, red, and green lines represent the ceph, teph, bare, and fric schemes.The shaded regions denote the standard deviation in the spread of data taken over multiple cascades.

Figure 5 .
Figure 5. Boxplot showing the distribution of defects remaining at the end of the cascades.The lower and upper bounds of the boxes are the first and third quartiles, the horizontal black lines the median, and the black dots the mean values of defects.The T-shaped whiskers indicate the last point, which is still within 1.5 times the interquartile range.Note, that the y-axis is in logarithmic scale.

Figure 6 .
Figure 6.Stacked-bar charts representing the average values of (a) atoms per clusters and (b) number of clusters per cluster's size range.The cluster sizes are classified into 5 ranges of 2, 2-15,15-30, 30-45, and more than 45 interstitials within the same cluster, as described in the label of (a).Note that interstitials classified as single, i.e. not part of any cluster, are disregarded from this analysis.

Figure 7 .
Figure 7. Boxplot showing the distribution of defects remaining at the end of the cascades.The boxes represent the cutoff (orange), and time offset (blue) schemes at different parameters, as detailed in the x-label.The lower and upper bounds of the boxes are the first and third quartiles, the horizontal black lines the median, and the black dots the mean values of defects.The T-shaped whiskers indicate the last point, which is still within 1.5 times the interquartile range.The cascades correspond to the 50 keV PKA.

Table 1 .
Simulations parameters.CIT and CET are the number of coarse-grained ionic and electronic temperature voxels (grid points).

Table 2 .
Cluster size analysis.The table presents the average size of the largest, mean, median, and number of clusters across different cascades.