Molecular dynamics simulations of the laser ablation of silicon with the thermal spike model

The purpose of this work is to model laser ablation of silicon on an atomistic scale in combination with a mesoscale model for the description of the electron-phonon interaction and an electron-temperature dependent interaction potential. The well-known continuum two-temperature model (TTM) for solids with highly excited electrons is extended from metals to silicon by explicitly taking charge carrier transport effects into account (nTTM). This is accomplished by the drift-diffusion limit of the Boltzmann-transport equation leading to the so called thermal-spike model (TSM). The model is further enhanced by extending the static modified Tersoff potential to a dynamical carrier excitation dependent interaction potential. We compare the TSM and nTTM with regard to physical correctness, numerical stability and applicability in the context of large-scale massive parallel high performance computing.


Introduction
The non-equilibrium phenomena in highly excited crystalline solids induced by strong laser radiation fields have received much attention in recent years. Despite of many theoretical and computational investigations these ultra-fast processes are still not well understood. This is especially the case for band gap materials, where laser induced ultra-fast structural changes were observed [1][2][3][4][5].
In the context of laser ablation simulations the main properties of interest are the spatially dependent absorption of laser light, the thermal conduction of the excited charge carrier sub-system and the characteristics of energy exchange of the Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. carrier subsystem and the ion-lattice. In the case of metals, these effects are most commonly described by the well established two-temperature model which leads to predictions of ablation thresholds and depths which are in good agreement with experimental data. In contrast to this, modeling of transport effects in highly excited semiconductors and dielectrics have been recognized as being more challenging. This is due to the complex nature of excited charge carrier dynamics in band gap materials, which in turn influence transport properties. While in metals transport properties are dominated by charge carrier temperature, transport effects in band gap materials show an additional dependence on charge carrier distribution. The charge carrier distribution does not only influence the local heat capacity of the carrier system, leading towards inhomogeneous thermal conduction but also influences the dynamic of the charge carrier system itself. For example these effects include charge carrier confinement [6] due to spatially inhomogeneous band gap narrowing. Band gap narrowing in turn influences charge carrier distribution which again, influences the spatial profile of laser light absorption which again influences local charge carrier generation rates. Furthermore spatially inhomogeneous charge carrier distributions will lead to charge carrier diffusion. Charge carrier diffusion is driven by charge carrier density and charge carrier temperature gradients and leads to energy transport by convection, ultimately determining the spatial dependence of the local energy transfer towards the ion lattice [7].
To tackle the problem of modeling this electronic subsystem the density dependent two-temperature model (nTTM) was introduced [8] and incorporates carrier-lattice interaction and electron-hole recombination processes but lacks energy convection terms. In this work we want to extend the nTTM towards the so called thermal-spike model (TSM) [9] in the context of laser ablation, leading to a combined self-consistent continuum-atomistic model, which explicitly includes energy convection transport effects as well as band gap effects. In addition, the excitation dependence of inter-atomic bonding strength of strongly excited covalently bonded semiconductors are taken into account by applying a carrier temperature dependent inter-atomic interaction potential.
In this work we use multi-million particle molecular dynamics (MD) simulations to study the laser ablation in covalently bonded materials by applying the nTTM as well as the TSM under explicit and semi-implicit solution schemes. The excitation dependence of inter-atomic bonding strength is modeled using the novel electron temperature dependent modified Tersoff-potential (MOD * ) [10], which enables us to incorporate laser induced bond-weakening into the model. This in turn allows for investigation of ultra fast phenomena like non-thermal melting on an atomic scale. In addition, the optical properties i.e. excitation dependent absorption mechanisms are modeled with inclusion of non-linear effects, namely two-photon absorption (TPA).
This work is structured into two main parts. First we will discuss modeling of laser ablation simulations by giving an introduction to combined continuum-atomistic simulations, as well as a comparison of the applied models. In the second part we will compare the predictions of the nTTM and TSM of the spatial temperature and carrier density distributions, as well as the numerical stability of both models.

Continuum-atomistic modeling
In this work the heat transport equations for charge carriers in semiconductors, as well as the particle transport equation have been solved by a continuum Finite-Difference (FD) method coupled to MD simulations. This section is divided into the modeling applied to laser field absorption 2.1, the modeling of the continuum transport equations (2.2), the coupling of these transport equations towards MD simulations 2.3 and a brief summary of the applied inter-atomic potentials 2.4.

Modeling of laser absorption
We begin with the modeling of laser pulse propagation. Assuming that the x-axis is the laser beam direction, the power density I of a laser pulse at the surface of the sample can be given as a Gaussian function in space and time: where R is the reflectivity, ϕ is the laser fluence, b is the laser beam diameter, t p is the laser pulse duration, r = y 2 + z 2 is the orthogonal distance from the optical axis and t 0 is the time at the maximal laser intensity. For spatially homogeneous distributed laser power density (b → 0) this equation simplifies to The absorbed energy density per propagation length is given by the rate equation [8] dI(x, t) dx Here α and β are the one-and two-photon absorption coefficients and Θ the free-carrier absorption cross section with the carrier density n c . This rate equation extends the well-known Lambert-Beer law and takes into account a variable number density n c of free carriers in covalently bonded semiconductors. Note that this model explicitly takes band-gap effects into account as the one-and two-photon absorption coefficient account for inter-band transitions. In contrast, the free carrier absorption cross section models the direct excitation of conduction band electrons.
For modeling the reflectivity and laser field absorption process, we use the Drude formula for the dielectric function [8,11] which leads to the reflectivity and the one-photon absorption coefficient of where ω L is the laser frequency, ν the collision frequency parameter, c the speed of light, ω p the plasma frequency,ñ the complex refractive index and ε r the relative dielectric constant. From non-linear optics we know the TPA coefficient of band gap materials as [12] with the band-gap scaling function K is a free parameter and was varied to reproduce a TPA coefficient of of β = 2 cm/GW at equilibrium conditions at 300 K [13]. Since the complex refractive index and the dielectric function are connected byñ = √ ε and both, the plasma frequency and the collision frequency ν are electron density dependent, the reflectivity R and one-photon absorption coefficient α are charge-carrier density dependent as well. Here e is the elementary charge and m e the electron mass. For the the collision frequency, the dominant processes of electron-hole (e-h) and electron-phonon (e-p) collisions are taken into account. The resulting collision frequency is given according to the Matthiessen rule by the linear combination of both collision processes. The expression applied for the electron-hole collision frequency was derived in [14]. The result in the non-degenerate formulation reads The used electron-phonon collision frequency was derived by [15] and reads where k s is a fitting parameter. In both equations T c and T ph respectively represent the electron-and phonon temperature, k b is the Boltzmann-and ℏ the reduced Plank-constant. m * e,DOS and m * h,DOS are the density of states effective masses for electrons (e) and holes (h). In addition laser field absorption leads to electron-hole pair generation, which is determined by the local laser field and the local charge carrier concentration. This dynamic is given by the rate equation where the last two terms correspond to Auger-recombination with the recombination coefficient γ and impact ionization with the ionization coefficient δ.

Modeling of the electronic subsystem
Modeling of the energy transport of the charge carrier subsystem for metals is typically done in the framework of the well-established two-temperature model (TTM). The assumption of a non-degenerate carrier gas leads to an internal energy U ∝ T c proportional to the carrier temperature T c . In the case of semiconductors the theory includes energy elevation by the band gap E g and reads Based on this, the TTM was extended to semiconductors taking into account the density of excited carriers to [8,[16][17][18] with κ c representing the carrier heat conduction coefficient. The term 1 τ e−ph 3n c k B (T c − T ph ) represents the carrier-phonon energy exchange rate in Fokker-Plank relaxation time approximation [18] with the energy relaxation time τ e − ph . The differential equation (22) for energy transport of the phononic subsystem is derived analogously and will be discussed shortly in Section 2.3. Both transport equations (15) and (22) combined with the rate equation for carrier excitation (13) lead to the carrier density dependent two-temperature model (nTTM).
The downside of this approach is that only heat transport effects of the carrier system are taken into account but no particle transport effects. For a quick estimation: At a room temperature of 300 K undoped mono-crystalline silicon shows a equilibrium charge carrier density is about n c ≈ 10 16 m −3 which is very small compared to a typical equilibrium carrier density of n c ≈ 10 28 m −3 for metals. In comparison this leads to a vanishing electron heat capacity under equilibrium conditions for semiconductors and therefore to high carrier temperatures under laser field absorption. Further more, according to (3) the carrier system excited under laser irradiation gets spatially non-uniform, leading to a spatially nonhomogeneous carrier density, specific heat capacity and band gap energy. Effects like carrier diffusivity and particle fluxes due to temperature and band gap gradients lead to energy advection, which on the other hand influences total temperature and energy distribution after laser irradiation. Due to the smaller equillibrium charge carrier density, these gradient driven effects are expected to be way more pronounced in semicondurctors, giving rise to the necessity to include them in the used transport model.
Starting with the drift-diffusion limit of the Boltzmann equation and ambipolar approximation, Driel [18] showed that the total carrier flux J can be written as with the ambipolar diffusion coefficient D 0 . The energy flux W can be written as a combination of the heat conduction by the well known Fourier law and the advection by the particle flux, leading to The continuity equations for the carrier density and energy read and Inserting the respective fluxes into these continuity equations leads to the so called thermal-spike model (TSM). Further insertion of the laser ablation specific absorption and the electron-phonon coupling terms leads to the applied system of equations: Note that for vanishing particle fluxes the TSM reduces to the nTTM.

Modeling of the phononic subsystem
The energy transport equation for the phononic subsystem corresponding to (15) can be derived analogously and reads with the phononic specific heat capacity C ph and phononic heat conduction coefficient κ ph . These two coupled partial differential equations, well known as the Two-Temperature Model (TTM) have been established for continuum modeling of electronic and lattice sub-systems after ultra-short laser irradiation. For a more detailed description of the non-equilibrium lattice dynamics and the inclusion of atomistic resolution we replace equation (22) by a set of MD simulation equations of motion for each atom i with mass m i at position r i namely These are Newton's equations of motion with a coupling term for the energy exchange with the carrier subsystem. The additional friction term models energy dissipation from the phononic subsystem via a virtual friction term against the thermal velocities v T i with ξ as the friction coefficient. The force F i = −∇ ri V acting on particle i is given by the interatomic potential V({r}) and the spatial configuration of the sample {r} and will be further discussed in section 2.4. From the requirement of energy conservation, the assumption (14), and the carrier-phonon energy exchange rate in Fokker-Plank relaxation time approximation, ξ can be directly derived as [18][19][20] Here m k and v T k are the mass and the thermal velocity of the k-th atom, N FD is the number of electronic iterations within a single MD step and N V is the number of atoms in a volume V FD . Within the framework of the Drude model, the electronphonon relaxation time is given as the inverse of the Drude collision frequency presented in equation (10).

Inter-atomic potential
For the simulations an electron-dependent interaction potential was developed. The starting point was the Modified Potential (MOD) [21].
In the original Tersoff interaction the total potential energy V was modeled as a sum of pair-like repulsive V R and attractive V A interactions with an environment-dependent coefficient b. The potential reads The modified angular-dependent term and cutoff function were introduced in the MOD potential to improve the melting temperature value. Under strong laser irradiation the anti-bonding states of covalently bonded semiconductors are occupied. As a consequence the potential energy surface and thus the interatomic interactions change almost instantaneously. The resulting inter-atomic forces can induce non-thermal processes in the lattice. To take these effects into account, the MOD potential for silicon was extended to a new interaction potential denoted MOD * where the parameters of MOD become carrier temperature dependent. A polynomial fit of the form was used to achieve a smooth dependency on the carrier temperature T c for each potential parameter P listed in table 2 [22].
Here the a 0 coefficients of MOD * correspond to the MOD parameters for silicon at a carrier temperature of T c = 0 K. The carrier temperature dependent interactions were derived from finite-temperature density functional theory (FTDFT) calculations [10]. For the fitting procedure of the effective MOD * interaction to the FTDFT results, the freely available software package potfit [23] was used. The cohesive energy, the lattice constants and the elastic properties of silicon depending on the electronic temperature have been calculated and can be found in [10].

Simulation results
The combined continuum-atomistic simulations of laser ablation were performed on two samples. The first sample represents an ultra-thin silicon film with a sample size of 82 nm × 11 nm × 11 nm containing 4.8 · 10 5 silicon atoms. The second sample represents bulk silicon with a sample size of 5500 nm × 5.5 nm × 5.5 nm containing 8 · 10 7 silicon atoms. All simulations were performed with a MD time step of 1.01806 fs. The simulation domain is sketched in figure 1. The simulation domain was divided into FD cells along the x-axis, which corresponds to the [1 0 0] crystallographic direction. Each cell possesses a depth of 11 nm in the x-direction. In the case of the bulk sample each of these FD cells contains approximately 10000 atoms. Periodic boundary conditions were applied in y-and z-direction and free boundariy condistions along the x-direction. The electronic-temperature dependent MOD * potential was used throughout this work.
The nTTM and TSM for the electronic system was solved on a regular finite difference grid with a time step of 0.001806 fs i.e. 1000 electronic iterations within a single MD step. The material parameters for silicon are are listed in table 3. The laser pulses were modeled with a Gaussian temporal profile according to (2) with the full width at half-maximum of 100 fs and a wavelength of 800 nm. This wavelength corresponds to a photon energy of 1.54 eV thus exceeding the the band gap. First we performed MD simulations at room temperature (300 K) and a static pressure of p = 0 Pa for a few thousand steps in order to reach an equilibrium atomic configuration. During equilibration, the carrier temperature was always kept at T ph = T c .
All simulations were implemented and carried out using our in-house laser-ablation-MD simulation software IMD (ITAP Molecular Dynamics of the former Institute for Theoretical and Applied Physics at the University of Stuttgart) [24,25].

Model comparison
3.1.1. Ultra-thin films. First we want to compare the nTTM and TSM with respect to excitation of the electron subsystem in the case of free-standing ultra-thin silicon film. For this, we use the simplest optical parameterization and set the TPA coefficient to zero. Further, we want to exclude ablation effects. Therefore, simulations of the ultra-thin silicon-film sample under homogeneous laser irradiation with a fluence below the damage threshold were performed. The nTTM and TSM was solved with the semi-implicit Crank-Nichelson algorithm. The sample exhibits a depth of 82 nm which is roughly a tenth of the characteristic absorption length of unexcited silicon at the applied wavelength of λ = 800 nm. In this simplified case, the sample gets homogeneously excited, thus transport phenomena as well as spatial carrier density gradients can be neglected. The temporal evolution of spatial averaged observables of the nTTM and TSM under sub-damage threshold laser irradiation is displayed in figure 2. It is visible that the nTTM and TSM show qualitatively and quantitatively identical excitation mechanisms when applied to an ultra-thin sample. All models show an electron temperature plateau at low laser intensities followed by rapid conduction band electron generation and excitation when the laser intensity reaches its peak value. Postpulse mechanisms are also equivalently featuring electron and   (7) band gap energy Eg lattice temperature equilibration due to electron-phonon coupling, as well as Auger-heating. Auger-heating originates from conduction band electron recombinations leading to a decrease of the specific heat capacity and therefore heating the electron subsystem. Because of the identical parameterization the carrier density dynamics is similar for both models. This still holds for β ̸ = 0.

Bulk.
In the case of laser irradiation of bulk silicon transport effects can no longer be ignored. Linear and nonlinear laser absorption processes lead to a spatially inhomogeneous carrier excitation, thus inducing carrier temperature and carrier-density gradients followed by inhomogeneous lattice heating. Laser ablation simulations on bulk silicon under homogeneous laser irradiation were performed. For the optical parameterization we chose the full Drude parameterization including contributions of electron-hole and electron-phonon collisions and an initial TPA of β = 2.0 cm GW −1 . For the laser parameters a fluence of σ = 320 mJ cm −2 and a FWHM of τ p = 100 fs was chosen. Figure 3 shows the time evolution of the nTTM and TSM observables within the first 10 ps of simulation time, averaged over the 11 nm-thick layer of surface FD-cells. It is visible that the nTTM and TSM no longer show identical carrier excitation in terms of locally generated carrier-density and temperature. The nTTM generates a higher carrier-density at the surface due to the lack of carrier transport. Since the heat capacity of the carrier subsystem is proportional to the carrier-density, this changes first, the local energy exchange rate between electronand lattice sub-system and second, the role of temperature in this context. In addition, the greater surface carrier-density leads to a spatially confined recombination process resulting in a vast carrier temperature difference due to Auger heating.  Table 4 shows the peak positions of the carrier temperature for both models. Of course, the peak temperatures reached differ drastically with varying laser parameters, but the nTTM yielding a ≈ 28% higher surface carrier temperature for this parameter-set demonstrates a significant difference. This is also reflected in a vast difference in the final lattice temperature T ph . While the TSM results in T ph = 1850 K, the nTTM predicts T ph = 2250 K.
The observed differences between the models may all be attribute to transport phenomena into the bulk. Despite identical carrier excitation at the time of laser peak intensity, a distinct difference in magnitude of transport effects within the models can be observed, leading to a difference in surface properties. Figure 4 shows the spatial distribution of carriertemperature and carrier-density along the x-axis at t = 3.5 ps simulation time. This corresponds to a delay of ∆t = 3.4 ps with respect to the time of peak laser intensity. The main transport contributions of each model are obtained by comparing the PDEs. For the nTTM and TSM the impact of the the definition of energy flux is directly visible in the particle flux definition The nTTM transport phenomena are purely driven by temperature gradients while for the TSM additional gradients in carrier density contribute to energy transport. On top of that carrier temperature gradients enhance the particle flux J which in turn enhances the energy flux further. This carrier temperature driven carrier flux leads to a more pronounced carrier wave propagating into the sample, again enhancing energy transport into the sample. This again leads to both, reduced surface carrier density and therefore reduced electron-phonon coupling and reduced surface temperature but a spatially extended laser affected region. Both models vary drastically in their predictions of the dynamic and final result of the spatial energy distribution within the ion-lattice. For a long time, the nTTM was the accepted standard transport model in the context of modeling laser ablation of covalently bonded semiconductors. This status was achieved because the observables of interest like damage thresholds and crater depths were predicted by the nTTM within the experimental error bounds. However, all those simulations used classic inter-atomic interaction potentials which where independent of electronic excitation. This approach should be scrutinized. The key mechanisms in modeling laser ablation phenomena are • parameterization of laser light absorption and electronlattice relaxation, • parameterization of transport mechanisms and • parameterization of inter-atomic interactions.
All these items greatly influence the resulting ablation depth in a mutually dependent way. As a result optimizing one of these parameterizations while keeping the remaining fixed can easily lead to the wrong conclusions. This is especially misleading in the case of the TPA coefficient, where experiments report values ranging from β = 0 to β = 55 GW/cm 2 [5,13,18,[26][27][28] for a wavelength of λ = 800 nm. Free choice of β within this range implicitly leads to a free choice of absorption length and therefore ablation depth, which in turn can be matched with experimental data. Intrinsically, the resulting parameter set reproduces the experimental data it was fitted to but does not necessarily include all relevant physical phenomena, thus showing false ablation dynamics. At this point we want to emphasize that our laser light absorption parameterization is carrier density dependent and included scattering effects are directly derived from theory rather than being experimental fitting-laws. The MOD * is constructed by FTDFT abinitio calculations rather than experimental data. This is a step towards a parameterization based on mutual independent theories, rendering the model predictions more trustworthy. FTDFT simulations suggest that all the silicon bonds become purely repulsive at an excitation of 2.14 eV, which translates into a carrier temperature of 24.834 K. This mechanism results in the prediction of lower ablation thresholds when the MOD * is used. The TSM, however, predicts lower excitation temperatures but a broader energy distribution, increasing the predicted ablation threshold. By performing laser ablation simulations under various fluences, models, and parameterizations, we found that the combinations nTTM+MOD, with a specificly chosen static TPA for each fluence and TSM+MOD * both predict the ablation depths and thresholds correctly. The combination nTTM+MOD * fails at this task, due to a combination of underestimating the thermal heat conduction but inclusion of the effects of non-thermal melting. The result is instantaneous ablation of the laser affected zone, without heat conduction into the sample.This underestimates the ablation depth. Only the combination of TSM+MOD * includes effects like non-thermal melting, ultra-fast phase transitions and sudden evaporation by excitation without being tuned to the desired results. The examination of resulting ablation dynamics of the combination TSM+MOD * will be part of our future work.

Simulation limits
The main challenge in performing simulations of ultra-fast laser ablation simulations of covalently bonded semiconductors is located in the equilibrium conditions. For metals, the assumption of an infinite charge carrier density contributing to energy conduction within the carrier subsystem holds, while for semi-conductors only excited conduction-band electronhole pairs contribute to energy transport and enter into the nTTM or TSM. This gives rise to a set of numerical problems.
In the case of explicit numerical solution schemes, the stability of the applied scheme can be determined with the well-known von-Neumann stability analysis. The von-Neumann stability analysis results in the condition CFL ≤ 1 with CFL being the Courant-number obtained by the stability analysis. The CFLs for the homogeneous part of both, the nTTM and TSM transport equations, are given in table 5. By this analysis it is evident, that the contribution of carrier density to the specific heat capacity within the conduction band leads to the inverse proportionality between the CFL number and carrier density. In our applied parameterization with a FD-Grid spacing of ∆x = 11 nm the CFL-number calls for a FD-time step of ∆t = 4.78 · 10 −27 seconds, which translates into the requirement of roughly 2.2 · 10 12 FD-steps per MD-step to achieve numerical stability. This is not suitable for a combined continuum-atomistic approach on today's supercomputers.
Other groups suggested to use semi-implicit solution techniques like the Crank-Nichelson algorithm, which is unconditionally stable for varying n c and arbitrary ∆t. The downside of semi-implicit solution schemes is that implementation is considerably more difficult task than the implementation of explicit solution schemes accompanied with increased computational cost. In an explicit forward-in-time-centered-in-space (FTCS) solution scheme the status of an FD cell at time t i + 1 can directly be calculated by the status of the surrounding FDcells at time t i . In case of a semi-implicit solution scheme the status of all FD-cells at both times t i + 1 and t i has to enter the solution scheme. This is done by creating and solving a linear equation system for the complete FD-lattice within each FDtime step. The semi-implicit approach was shown to be fruitful to predict charge-carrier dynamics in laser excited silicon for small timescales [8,27].
However the semi-implicit solution scheme shows an extreme downside for longer simulation times, especially when combined with MD. The MD part of the simulation is necessary to model the occuring phase transitions throu nonequillibrium pathways thus making predictions of the state of sample after laser ablation took place. At fluences greater than the ablation threshold we want to investigate effects like evaporation, phase explosions, melting and general material removal. This results in the necessity of FD-cell deactivation when the corresponding material within the FD-cell is removed. In the context of a semi-implicit solution scheme this directly leads to a under-determined system of equations where convergence of the respective solutions is not given anymore.

Critical carrier-density.
When silicon is excited by laser light with a suitable wavelength, the carrier density increases suddenly within the laser effected zone by many orders of magnitude. The result is that explicit solution schemes get applicable again for already excited silicon.
We performed laser ablation simulations under varying starting carrier densities n eq , with both the nTTM and TSM model, and with both the semi-implicit Crank-Nicholson algorithm and an explicit FTCS scheme. The usage of the semi-implicit Crank-Nicholson solution algorithm allows to surpass the stability limit of an initial carrier density n crit eq given by von-Neumann stability analysis. All simulations were carried out with sub-ablation threshold laser fluences and a simulation time of 500 ps.
We observed that for values of n eq smaller than the expected maximum excitation n max no relevant changes in electronic excitations given by (T c , n c ) or lattice dynamics given by (T ph ) are detectable in our simulations. All other observables are functions of (T c , T ph , n c ), and intrinsically give the same result for identical arguments.
The results of the simulations performed for varying initial carrier densities n eq are shown in figure 5. Here T l denotes the final lattice temperature T ph (t → ∞) averaged spatially over the whole sample and averaged temporally over the last 10 ps of simulation time. RMS denotes the root mean square deviation of the lattice temperature against the averaged lattice end temperature and Diff denotes the difference between lattice end temperature, obtained with the correct initial carrier density of n lit eq = 10 10 cm −3 , to the resulting lattice end temperature obtained with the corresponding initial carrier density n eq . In our parameterization n eq was determined to n eq = 10 19 cm −3 . For initial carrier-densities in excess of n eq = 10 19 cm −3 , an artificial electron temperature surplus is generated by Auger-heating even for pre-pulse timescales. In addition the artificially elevated initial carrier density greatly affects the free carrier absorption coefficient α leading to a greater electron temperature increase as well as free carrier generation, which in turn again leads to greater Auger-heating. It is noteworthy that it can be assumed that the influence of the surplus charge carrier generation on the free carrier absorption coefficient will greatly affect the absorption length on simulations of bulk systems and therefore manipulate spatial energy deposition. For n eq ≤ 10 19 /cm 3 the difference between the resulting lattice end temperature and the lattice end temperature obtained with n eq = n lit eq is several orders of magnitude smaller than the root mean square deviation of the temperature fluctuations of the MD cells. Thus showing identical results within the error bars of statistical fluctuation.
In conclusion this means that an explicit FTCS scheme with an artificially high n eq will give the same results in the context of TSM and nTTM as an semi-implicit Crank-Nicholson scheme with the correct n eq . In addition, explicit schemes are usually easier to implement and are also computationally cheaper, especially in the context of large-scale high performance computing when inter-CPU communication is crucial for performance. This renders an explicit FTCS scheme with an artificially high n eq as the way to go. We want to emphasize that our code IMD is capable of solving the TSM on a three-dimensional FD-lattice using this approach. This opens the new opportunity to perform large scale laser ablation simulations on covalently bonded semiconductors with spatially non-uniform laser profiles and investigation of the results of the MOD * for the ablation crater. To our knowledge this was not possible for covalently bonded semiconductors before.

Stability observations
Besides the physical advantages of the TSM including convection terms, a purely numerical advantage is noteworthy. In conclusion of section 3.2.1 it is legitimate and recommendable to use explicit simulation schemes in the context of laser ablation simulations. We applied the nTTM and TSM, both solved with an explicit FTCS scheme, to a FD-MD simulation run on the bulk sample. In figure 6, the resulting carrier density and the carrier temperature after laser irradiation with a fluence of σ = 750 mJ/cm 2 are shown at different simulation times. For various parameter sets, we observed that in case of the nTTM, numerical instabilities occur at the rear of the sample, where carrier excitation is comparable low. These instabilities propagate towards the surface, causing the simulation to crash. In contrast, the TSM behaves numerically stable for the same parameter set. The coupling of the transport equations to MD-fluctuations makes it appropriate to view the problem as a system of stochastic non-linear partial differential equations. Since a rigorous and precise mathematical stability analysis for such a system is a very complex and protracted task, we want to give a phenomenological reasoning of the increase in stability. First, we ignore all non-diffusive parts of the respective homogeneous transport equations of the nTTM and the TSM, allowing us to define a Courant-number (CFL) [29] for stability analysis. The resulting CFL-numbers are shown in Table 5. First we observe that the CFL-number for the nTTM is inversely proportional to the carrier density which is vanishing in the case of equilibrium conditions of silicon. This explains why the instabilities occur at the back of the sample where laser excitation and therefore charge carrier excitation shows the smallest impact, resulting in a more demanding CFL. Furthermore, we observe that in the TSM the drift-terms contribute to the CFL-number in a carrier density independent way, greatly increasing stability in regions of small carrier densities.
Secondly, in both models the band gap E g is dependent on the MD-temperature and therefore shows statistical fluctuations. Laser field absorption yield a MD-temperature dependence, causing the carrier-as well as the temperature distribution after laser irradiation to reflect these fluctuations. The nTTM only uses the rate equation (13) for carrier dynamics, allowing no diffusion of charge carriers and leading to an accumulation of noise. This noise enters the computation of heat propagation favoring numerical instabilities. The TSM on the other hand, includes diffusion terms for carrier dynamics as well as drift terms for heat propagation, which lead to a smoothing of temperature and carrier distributions, averaging numerical noise.

Summary and conclusions
In this work we compared the charge carrier concentration dependent two-temperature model (nTTM) with the thermalspike model (TSM) in the context of laser ablation simulation of covalently bonded semiconductors. The models were implemented into our in-house simulation software package IMD. We carried out laser ablation simulations on samples of size 5500 nm × 5.5 nm × 5.5 nm and 82 nm × 11 nm × 11 nm, under laser irradiation of a 100 fs FWHM pulse with varying fluence over a time-span of up to 100 ps. In all cases the nTTM and TSM were solved using two methods, an explicit FTCS scheme and a semi-implicit Crank-Nicholson scheme. In the case of ultra-thin films where transport effects can be neglected, the nTTM and TSM yield identical results. However when applied to bulk samples, transport effects lead to different spatial energy distributions. As a result the nTTM predicts higher surface temperatures than the TSM model, while the TSM model shows a spatially extended energy distribution because of the incorporation of carrier-diffusivity. We observed that the combinations of the nTTM with the MOD or MOD * are not suitable for laser ablation simulations on bulk samples since it either reproduces the wrong ablation mechanisms or ablation depths. In contrast, the TSM with the electron excitation dependent modified MOD * yields the correct ablation depths with incorporation of the experimentally observed ablation mechanics, especially the mechanic of ultra fast non-thermal melting [1][2][3][4][5]. We showed that the combination TSM+MOD * is the more physical parameterization due to the incorporation of excitation effects. On top of that we demonstrated the superiority of the TSM in terms of numerical stability. In the past the instability of the nTTM was recognized [17], which lead to the implementation of the nTTM with semi-implicit solution schemes [8,27]. The usage of the semi-implicit algorithms in the context of continuum-atomistic modeling leads to two problems. First semi-implicit solution schemes are more complicated to implement and show a higher computational cost. Secondly when coupled to a MD simulation it is necessary to deactivate FD-cells when the atoms contained in these cells are ablated, leading to under-determined equation systems. By comparing the applied models under variation of initial carrier density and different solution schemes we demonstrated the existence of an artificial high initial carrier density, at which an implicit FTCS solution scheme becomes numerically stable without deviation from the solution obtained by a more complicated and computational more expensive semi-implicit algorithm. This exposes the main advantage of the TSM in the context of laser ablation on band-gap materials: For circumstances where the already established nTTM is applicable, it yields identical results while being computational more efficient. In circumstances where the nTTM fails [27] because of the need for semi-implicit solution algorithms, namely at above ablation threshold fluences, the TSM is stable and solvable with simpler explicit solution schemes. This makes the TSM more suitable for high performance computing in the context of laser irradiated covalently bonded semiconductors, and opens the possibility of simulation of three dimensional heat propagation and investigation of consequent phenomena like ablation mechanisms, behavior of the ejected material or shock wave propagation in novel laser generated geometries like ablation craters.