Three-temperature modeling of laser-induced damage process in silicon

In the recent years, we have been able to use intense laser pulses of duration ranging from femto- to pico-second time scales to investigate laser-matter interaction. In particular, laser processing of semiconductors using ultrashort pulses attracts great interest, due to its application in nano-structuring [1, 2]. Ultrafast laser pulse enables precise processing of the material, without causing damage in the surrounding area [3].

Laser processing involves various aspects of physics such as optics, quantum mechanics, material science and thermodynamics, to name a few. Thus, numerical modeling is crucial in order to understand and optimize the process. Numerical modeling of laser excitation has seen considerable interest over the years [4, 5, 6, 7, 8, 9]. In particular, the processing of silicon has important application, not only in understanding the fundamental science but also industrially [10, 11].

When the laser pulse is shorter than the time-scale of phonon excitation, electronic excitation by photo-absorption is the dominant process. The non-equilibrium phase with hot electron-hole pairs and cool lattice relaxes into thermal equilibrium via electron-phonon interaction.

In order to study non-equilibrium dynamics of electrons and lattice, we have to develop a numerical model including dynamics of electro-magnetic fields, electrons, lattice, and the interaction among them. The Two-Temperature Model (TTM) has been proposed by Anisimov et al. to describe the electron-lattice coupling in metals [12]. The TTM has been extended to study excitation in semiconductors [5] and different models have been developed, one of them being the self-consistent density-dependent TTM (nTTM) [6]. One of the critical approximations in TTM for semiconductors is assuming thermal equilibrium in electrons. There are various differences in the TTMs for semiconductors found in existing literature, such as the treatment of excitation process and the values of optical parameters of silicon [4, 8, 9, 13].

It has been observed that electron and hole relaxation in conduction and valence bands is faster than relaxation of the system as a whole [14] . As the electron temperature increases, electron-hole scattering frequency decreases significantly [15]. The calculations presented using first-principles numerical simulation show that electrons and holes have different energy, depending on the band structure. The energy difference between electrons and holes suggests that quasi-temperatures are different in the conduction band (CB) and valence band (VB) [16]. Also, the decreased electron-hole scattering frequency indicates that it is important to take quasi-temperatures of electrons and holes into consideration.

In this work, we would like to propose a modified nTTM treating the evolution of quasi-temperatures in CB, VB and lattice (3TM). A similar study has been presented by Silaeva et al. [17], in order to study laser assisted atom probe tomography.

The 3TM is modeled to consider three sub-systems, electrons, holes and lattice, and their temperatures are calculated separately. The carrier densities are also calculated distinctively for electrons and holes. We treat the dynamics of laser field employing the one-dimensional Maxwell’s equations, which are solved using the Finite Difference Time Domain (FDTD) approach. In order to understand and define damage in silicon, we assume several processes that may be responsible for structural changes and compare their calculated thresholds with the experimental damage thresholds.

The 3TM is similar to the previous nTTM formalisms [6] in its formulation. The major difference lies in the separate evolution of electron and hole temperatures and densities. 3TM also considers the changes in optical parameters of silicon during the excitation processes which in turn, affect the evolution of carrier densities and temperatures. The single and double photon excitation processes have been modified from previous studies to include the effect of band re-normalization. Drude model is used to calculate dielectric function, including the effect of electron-hole-phonon collisions.

The time-evolution of electron and hole densities, $n_{e}$ and $n_{h}$ is described as:

where $\omega_{0}$ is the laser frequency and $\alpha$ is the single photon absorption coefficient for transition from VB to CB [18]. $\beta$ is the two-photon absorption coefficient for which we use the DFT calculation when $2\hbar\omega_{0}>E_{g}$ [19] and around the band gap energy we employ interpolation to the model described in Ref.[20]. $\gamma_{e(h)}$ is the Auger re-combination coefficient [17] and $\theta_{e(h)}$ is the impact ionization coefficient [6]. Equation 1 also includes the effect of spatial charge distribution and the associated electric field, and $J_{e(h)},D_{e(h)}$ and $\vec{F}$ are the charge current, diffusion coefficient and the electric field induced by the electron–hole separation, respectively [21].

The total dielectric function along with the effect of band structure re-normalization [22] is expressed by

where $n_{0}$ is the density of valence electrons. It should be noted that $n_{e}$ and $n_{h}$ are nearly the same due to the effect of $\vec{F}$ and can be approximated at $n_{e}$. $\epsilon_{L}(\omega)$ is the innate dielectric function, $\delta E_{g}$ represents the band re-normalization by carrier density, and $\epsilon_{D}$ is the complex dielectric function calculated from Drude model. The temperature dependent optical parameters of silicon are referred to from Ref.[18]. $\epsilon_{D}$ accounts for the effect of plasma in the excited system. Considering the electron and hole sub-systems,

where $m^{*}_{e(h)}$ is the effective mass and $\nu_{e}$ and $\nu_{h}$ are the collision frequencies, describing the electron-hole ($e-h$), electron-phonon ($e-ph$) and hole-phonon ($h-ph$) collisions. The $e-ph$ and $h-ph$ collisions are assumed to have the same frequency which is dependent on lattice temperature [13]. Effect of electron-ion core collisions is also considered [23, 21]. The collision frequencies for $e-h$ interactions are calculated as per the model presented in Ref.[15]. The total one-photon absorption coefficient including free-carrier absorption is

where $\lambda$ is the laser wavelength.

Since 3TM considers three sub-systems, viz. electron, hole and lattice, their temperatures also evolve separately. The temperature evolution is expressed as:

The third and fourth terms in Eq. (5) account for the loss of energy due to electron-lattice interaction and energy current. The last two terms on right hand side include the changes in carrier density and band gap energy. Here, $H_{\xi}^{\zeta}(\eta)=F_{\zeta}(\eta)/F_{\xi}(\eta)$ and $F_{\xi}(\eta)$ is the Fermi integral [21]. The heat capacities $C_{e(h)}$ are calculated from the carrier densities and temperatures [21]. $T_{l}$ is calculated following the empirical model, where the term for carrier temperature is replaced by terms for electron and hole temperatures, as described in Eq.  (6). Here, $\kappa_{l}$ is the thermal conductivity [21].

The propagation of the laser pulse is described by solving the Maxwell’s equations using FDTD method. Mur’s absorbing boundary condition is employed to prevent reflection from the boundary [24]. Another salient feature is that the electric field is considered to be complex for the calculation of laser intensity, to ensure a non-zero field at all points in time and space. Assuming a one-dimensional system, the electric field is:

Here $E_{0}$ includes Gaussian envelope defining the shape of the pulse as $\exp[-(t-4T)^{2}/T^{2}]$, where $T=t_{p}/(4\sqrt{\ln 2})$, $t_{p}$ being the FWHM pulse duration.

where $I(x,t)$ is the laser intensity. Evaluation of charge current induced by the laser field is a crucial part of the module. We calculate the current with and without excitation i.e., for photo-absorption and dielectric response. For dielectric response, $j_{0}(x,t)$ is calculated as:

where $A(x,t)$ is the vector potential, $P$ is the polarization and $\chi_{r}$ is the real part of susceptibility ($\chi=(\epsilon-1)/4\pi$). The Maxwell’s equation thus, becomes

where

is the current associated with photo-absorption.

The time evolution of 3TM is calculated using $4^{th}$ order Runge-Kutta method. Euler’s method is used in the FDTD method to solve Maxwell’s equations. In our study, the 3TM time step ($dt$) is of the order $10^{-1}$ fs. For FDTD calculation, we define the another time step , $dt_{m}\sim 10^{-3}$ fs, because the time-scale of electromagnetic field dynamics is faster than that of electron dynamics. For spatial parameters, the simulation grid size is taken to be 150 Å  and the silicon film thickness is 60.5 $\mu$m, which is thicker than the penetration depth, which is $\sim$10 $\mu$m for 775 nm laser.

Figure 1 shows the time evolution of electrons, holes and the lattice dynamics in silicon surface for interaction with a 775 nm, 350 fs pulse having incident fluence 0.18 J/cm². We assume the one-dimensional system for 3TM. Figure 1(a) shows the electron density, along with the pulse profile in the background. The steep rise in $n_{e}$ is due to photo-absorption, followed by fall in density due to Auger re-combination. The electron, hole and lattice temperature evolution in Fig. 1(b) shows that the three sub-systems eventually reach equilibrium at $\sim 1.85\times 10^{3}$ K. Initially the carrier temperatures rise steeply due to photo-absorption and subsequent plasma heating. This is followed by re-combination and carrier-phonon dynamics which re-distribute the energy. There is transfer of energy from electrons and holes to the lattice, causing the lattice temperature to rise. Figure 1(c) shows the evolution of electron energy and the total energy per atom. Since the penetration depth for the pulse wavelength of 775 nm pulse is high, the gradients of $n_{e(h)}$ and $T_{e(h)}$ are small. Small $\nabla n_{e}(r)$ and $\nabla T_{e}(r)$ suggest low charge and energy current. Therefore, the total energy is nearly constant with time, which also demonstrates the stability of the 3TM simulation.

The damage in a semiconductor is defined in different ways throughout the published literature [25, 26, 27, 28, 29]. Broadly speaking, the structural and phase changes observed in the lattice of laser excited material are defined as ’damage’. More specifically, there are different processes leading to damage, depending on the conditions of laser excitation and material properties.

One of the most prominent processes for damage is thermal melting. The carrier-phonon interactions lead to re-distribution of energy among the carriers and lattice, causing the lattice temperature to rise up to melting point, which can manifest as damage in the material [15]. Thermal melting threshold can be studied from two perspectives, initial commencement of melting which leads to what is known as the ’slush’ phase (partial melting) [4], and the completion of melting in crystalline silicon (c-Si) when bond configuration changes to tetrahedral after absorption of latent heat from the carriers (complete melting). The two different approaches to study melting in silicon [30] lead to different understanding of the onset of damage.

Critical density has also been considered as a crucial parameter to study damage [31] and the threshold is considered when the electron density reaches the critical density and plasma state exists in the material. The Drude-model indicates an abrupt increase in absorbed energy above the critical density threshold, which leads to sudden increase in melting thresholds.

The first-principles calculation by the time-dependent density functional theory (TDDF) indicates that the ablation threshold of $\alpha$-Quartz lies between melting and cohesive energy [32]. The TDDFT results suggests the bond-breaking energy is a reasonable candidate for the damage threshold. Threshold for breaking of bonds is calculated as the threshold fluence when energy per atom reaches 2.3 eV, which is half the cohesive energy of a silicon crystal [33].

Emission of electrons from the surface into vacuum is also an important candidate for understanding damage. When the electron temperature is high, electron emission (e-emission) from the surface induces an impulsive force, Coulomb explosion being the extreme case [34]. We estimate the e-emission threshold as the threshold when the electron energy exceeds the work function for silicon (4.65 eV).

Figure 2 shows the numerical results for damage threshold by 3TM with a 775 nm laser. Experimental data is also included in the figure. Data from Allenspacher et al. shows results for damage threshold in laser processed silicon using 775 nm pulse of varying duration [25], Bonse et al. data shows ablation damage threshold for a 800 nm pulse of duration 130 fs [26] and Izawa et al. data shows the crystallization threshold range for a $800$ nm pulse of duration 100 fs [27]. The damage in these experimental studies are defined differently, which is why there is a considerable difference in the threshold data. Allenspacher et. al. [25] observed morphology of the surface and threshold was estimated as the damage probability becomes zero. Bonse et. al. [26] defined damage threshold as the fluence at which area of damaged region becomes zero by extrapolating the fitted line. Izawa et al. [27] defined a range of incident fluence for crystallization of amorphous silicon (a-Si) on c-Si.

In Fig. 2, the experimental data suggests weak $t_{p}$-dependence in the short pulse regime (SPR, $t_{p}<1$ ps), while it is sensitive to $t_{p}$ in the long pulse regime (LPR, $t_{p}>1$ ps). In SPR, experimental results for threshold lie between partial melting and complete melting. At $\sim 1$ ps, there is a reversal in the partial melting and e-emission threshold, and e-emission coincides with the experimental results in LPR. The $t_{p}$-dependence for LPR indicates a power law of $t_{p}^{\eta}$, which has also been reported with fused silica and CaF2 ($\eta=1/2$) [35]. In the LPR, the experimental data fits well with a $t_{p}^{0.15}$ dependence, and the e-emission threshold fits with a similar power law of $t_{p}^{0.14}$. The agreement of e-emission threshold with experimental data indicates that the impulsive Coulomb force induces damage above 1 ps. The $t_{p}$ dependence of threshold suggests that thermal melting and electron emission are necessary and sufficient conditions for causing damage in silicon.

A possible explanation to the damage process is a cooperative relationship between bond softening by excitation and Coulomb force due to electron emission. The melting threshold indicates Si-Si bond-softening, while electron emission indicates the presence of impulsive Coulomb force. Although the melting process takes time, the impulsive Coulomb force can enhance the lattice dynamics. On the other hand, Coulomb force may not be strong enough to break the cold Si-Si bonds. The thermal melting and e-emission thresholds cross at $\sim 1.4$ ps, indicating that both these processes are together responsible for causing damage in silicon.

The amorphization threshold of Bonse et al. coincides with the calculated threshold for breaking of bonds. The Izawa et al. results give a range of crystallization for amorphous silicon. As it can be seen in Fig. 2, the threshold for crystallization lies between the melting and bond-breaking threshold.

Although considered as an important measure of studying damage, we find that in our model critical density does not qualify as the dominant reason for damage in silicon. When the plasma frequency increases beyond laser frequency, photo-absorption by free carriers is assumed to be significant. As it can be seen in Fig. 3, the laser fluence dependence of the absorbed energy shows monotonic increase for various pulse duration. However, we found that the fluence dependence can be fitted by the quadratic function (black solid lines in Fig. 3). The quadratic increase of the energy above the critical density thresholds indicates that the critical density cannot be a favorable criterion for laser damage.

The e-emission is an important process in LPR and its corresponding threshold indicates a significant pulse duration dependence, as compared to thresholds related to melting. Since $T_{e}$ defines the e-emission threshold, it should depend on the $e-ph$ relaxation time. Figure 4 shows the threshold for thermal melting and e-emission for different $e-ph$ relaxation time $\tau=\tau_{0}[1+\{n_{e}/8\times 10^{20}(1/$cm${}^{3})\}^{2}]$, by taking the constant $\tau_{0}$ as $\tau_{0}=240$ fs and 500 fs. The e-emission threshold shows a qualitative shift towards higher $t_{p}$ dependence, which also depends on $\tau_{0}$. The threshold can be studied in two parts, where higher dependence is seen for longer $t_{p}$. We present the fitted line with $t_{p}^{\eta}$ for e-emission threshold in LPR, where $\eta=0.14$ and $0.13$ for $\tau_{0}=240$ fs and $500$ fs, respectively. The deviation from the fitted line occurs at $1.5\tau_{0}$ for $\tau_{0}=240$ fs and $1.4\tau_{0}$ for $\tau_{0}=500$ fs, respectively. The crossing points of thermal melting and e-emission threshold occur at $5\tau_{0}$ and $6\tau_{0}$ for $\tau_{0}=240$ fs and 500 fs. Due to slower energy transition in case of $\tau_{0}=500$ fs, the crossing point shifts higher. Also, since the $t_{p}$ dependence follows almost the same power law for both cases of $\tau_{0}$, it shows that the physical process is almost the same when $t_{p}$ is sufficiently long. The threshold for thermal melting is not affected by the change in $\tau_{0}$, suggesting that the temperature dependence of photo-absorption by electron-hole system does not affect the lattice temperature significantly.

In this work, we studied the laser excitation and possible damage processes in silicon numerically using 3TM combined with Maxwell’s equations. The model allowed us to study in detail the excitation dynamics and energy re-distribution in silicon for electrons, holes and the lattice. Pulse duration dependence of the damage threshold was also studied in detail. The 3TM simulations reproduced the experimental damage thresholds and the qualitative dependence on incident pulse duration. The comparison of experimental threshold alongside the calculated thresholds indicated that thermal melting and e-emission are both important processes, contributing to damage in silicon. However, critical density threshold is not an equitable criterion for damage. Also, variation in the $e-ph$ relaxation time showed that $e-ph$ interaction dynamics plays an important role in damage process with longer pulse duration.