<figure><img src="image_1.png" /></figure>

HAL
open science

# Unveiling Ultrafast Dynamics: Atomistic-Continuum Modeling of Laser-Induced Phase Transitions in Silicon

Dmitry S Ivanov, Tatiana E Itina

► To cite this version:

Dmitry S Ivanov, Tatiana E Itina. Unveiling Ultrafast Dynamics: Atomistic-Continuum Modeling of Laser-Induced Phase Transitions in Silicon. Physical Chemistry Chemical Physics, 2025, 10.1039/D5CP00964B. ujm-04784471v3

HAL Id: ujm-04784471

https://ujm.hal.science/ujm-04784471v3

Submitted on 18 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unveiling Ultrafast Dynamics: Atomistic-Continuum Modeling of Laser-Induced Phase Transitions in Silicon

Dmitry S. Ivanov¹ and Tatiana E. Itina¹*

¹ Université Jean Monnet Saint-Etienne, CNRS, Institut d'Optique Graduate School, Laboratoire Hubert Curien UMR 5516, F-42023, Saint-Etienne, France.

*tatiana.itina@univ-st-etienne.fr

## Abstract

Modern semiconductor applications demand precise laser processing at the nanometer scale, requiring a detailed understanding of phase transition and structural modification mechanisms. Accurate control over laser-induced processes in semiconductors is essential for generating predesigned surface structures and modifying the surface properties. In this study, we present a numerical investigation of non-equilibrium laser-induced phase transitions in silicon (Si) using a hybrid atomistic-continuum model. The model combines the strengths of Molecular Dynamics (MD) simulations for atomistic-scale descriptions of non-equilibrium phase transitions with a continuum approach to account for the effect of laser-generated free carriers. As compared to the ordinary continuum or MD approaches, this advanced framework, therefore, captures the kinetics of melting and ablation phenomena on one hand, and generation and diffusion of the electron-hole pairs, thermal diffusion, and the electron-phonon coupling processes during laser energy deposition on the other hand. We applied the model to determine the melting depth as a function of fluence for a 100 fs laser pulse at 800 nm. The results show that the stand-alone continuum approach underestimates the melting threshold as compared to the hybrid atomistic-continuum model by 46% originating from the detailed description of the melting kinetics, inherently given by MD. Additionally, we explored the effect of crystal orientation on melting dynamics and compared the results with the corresponding experimental measurement. Finally, the MD model is used to identify the conditions leading to the amorphization of the Si surface. The determined corresponding cooling rates are referred to the conditions realized in the experiments on femtosecond laser-silicon interaction. These findings provide valuable insights into experimental observations of Si surface structuring induced by ultrashort laser pulses.

**Key words:** laser-material interactions, laser ablation, laser-induced melting, surface amorphization, atomistic-continuum model, ultrashort laser pulses, free carrier dynamics

## 1. Introduction

Silicon, the cornerstone of modern microelectronics, continues to push the boundaries of technological advancement, especially in optoelectronics, and solar cells industry. As the demand for faster, smaller, and more efficient devices continues to grow, being able to control material properties at the nanoscale is crucial. One of the keyways to achieve this is laser processing, which allows us to modify materials with high accuracy and on a nanometer scale [1,2]. In
microelectronics, lasers are used for many purposes, including doping, annealing, and etching of semiconductors [3]. These processes are essential for creating the tiny, complex patterns needed in integrated circuits. Ultrafast laser pulses, with their short duration, allow for extreme precision, minimizing damage to surrounding material. The ability to create localized phase changes, such as surface melting or amorphization, allows fine control of material properties, such as electrical conductivity or optical properties [4,5]. Laser-induced amorphization of silicon is particularly useful for phase-change memory, a technology used in data storage, where switching between the crystalline and amorphous states can store information [6]. Understanding how these laser-induced phase transitions happen is of extreme importance for improving current and future microelectronics technologies [7].

Many experiments have been done to study how short and ultrashort laser pulses affect silicon. These studies have shown that ultrafast pulses can cause rapid melting of silicon, followed by the formation of an amorphous layer if the material's cooling rate is too high for recrystallization. Using ultrafast time-resolved techniques, scientists can observe how silicon melts and re-solidifies in real-time. The melting process depends on factors like the energy of the laser pulse (fluence), laser wavelength, pulse duration, and crystal orientation. For example, at fluences high enough, *homogeneous melting* can occur, where the entire structure melts almost at once, driven by significant overheating exceeding crystal stability [8]. In contrast, *heterogeneous melting* starts at surfaces, defects, dislocations, or imperfections in the material, requiring lower temperatures and occurring unevenly across the sample. But even more, pressure and temperature dynamics play the key roles in these mechanisms—higher temperatures tend to favor homogeneous melting, while pressure variations can reinforce the crystal stability against of melting or weaken it. Also, heterogeneous melting, depending on the presence of defects, is always limited in the speed of solid-liquid interface propagation, and can be overwhelmed by homogeneous mechanism of melting [9]. The experimental measurements are also limited in what they can reveal about the atomic-scale events, so we rely on simulations for complementary research of these processes. To advance the understanding of the fundamental processes occurring in laser-irradiated silicon and to benefit from insights into how laser parameters influence the underlying restructuring mechanisms, we utilize a robust computational model in this work. The model captures laser-induced phenomena in the solid state and can be readily employed to interpret experimental measurements. Its application can thus support the design of experiments aimed at generating structures with desired properties.

To better understand the non-equilibrium phase transition mechanisms, several models were developed that simulated the complex laser-induced processes in solids, including semiconductors [3, 13, 10, 12]. One of the common and most used methods is the Two-Temperature Model (TTM), which shows how the deposited energy is transferred between electrons and the atomic lattice after the laser pulse. However, this stand-alone model can't capture important atomic-level details, such as how melting and amorphization occurs. In other words, the kinetics of melting is not included *a priory* in the continuum description of the above process, and number of add-on models and approximations are frequently utilized. Also, under extreme conditions, realized during the ultrashort laser pulse interaction with semiconductors, in the presence of the resulting string pressure and temperature gradients, and when the size of new phase nuclei getting as small as just several interatomic distances, the applicability of the continuum approaches is hampered. On the other hand, the atomistic approach with Molecular Dynamics technique seems to be more appropriate for the detailed description of the kinetics of the laser-
induced phase transition process, including the transition to amorphous state [15], but this method
does not account on the effect of free laser-generated free carriers (the electron-hole pairs), the
presence of which can play a determinant role during the ultrashort pulse laser processing of
semiconductors. To solve this, researchers combine TTM with Molecular Dynamics (MD)
simulations. In particular, because of the additional descriptor of the free carrier density *n*, for
semiconductors this hybrid approach, known as the nTTM-MD model [9], combined the
advantages of both methods into the frame of a single computational approach. It accounts for
multi-photon and free carriers' absorption of the laser energy, the energy transfer from hot carriers
to the lattice, fast electron thermal diffusion, the carriers' diffusion process, and in the same time
it reveals small-scale atomistic description of the non-equilibrium phase transition. The nTTM-
MD model is especially useful when accurate description of the laser light absorption, the energy
transport, melting, and amorphization processes become essential.

One of the biggest challenges in modeling laser-induced phase transitions in silicon is an
accurate prediction of the melting threshold and depth, especially at lower fluences. Traditional
continuum models often underestimate the extent of melting because they don't fully account for
atomic-scale processes. By using the nTTM-MD model [9], a much better prediction of melting
dynamics can be obtained, particularly for ultra-short laser pulses at low laser fluences. This model
correctly reproduces non-equilibrium phase transitions, so that we explore the difference between
heterogeneous and homogeneous melting and show how each mechanism leads to different
outcomes in melting depth. To highlight the advantages of the combined MD-nTTM approach in
describing the non-equilibrium laser melting process, we compare the melting depths obtained
using the combined model with those from the stand-alone nTTM approach. We also study the
effects of crystal orientation on melting. Our results show that crystal orientation significantly
affects how silicon melts, with different orientations producing different melting patterns and
rates.worthy, the effect of crystal orientation can be barely accounted for in the alone nTTM model.

Another important challenge is the precise control over the depth of amorphization. While
laser-induced amorphization can be useful in some cases, such as memory storage, depending on
its depth it can be either an advantage or a disadvantage in microelectronics. In fact, amorphization
might reduce the electron mobility or change the material's electrical properties. Reducing the
depth of amorphization without sacrificing the ability to make fine changes to the material's
structure is a key goal of the current experimental and theoretical research. Thus, our study focuses
on refining the hybrid nTTM-MD model to better control melting and recrystallization processes
while minimizing unwanted amorphization.

In this paper, we use a hybrid atomistic-continuum model (nTTM-MD) to study how free
carrier temperature, lattice temperature, and free carrier density evolve during ultrafast laser
excitation at a wavelength of 800 nm. We also investigate the relationship between melting depth
and fluence for a 100 fs laser pulse, identifying the roles of heterogeneous and homogeneous
melting. By comparing atomic configurations of silicon samples with <001> and <111> crystal
orientations, we show how crystal structure affects the melting process. Finally, we discuss
strategies for reducing amorphization depth.

The paper is organized as follows: Section 2 explains modeling details for both nTTM and
the hybrid nTTM-MD model, parameters and numerical methods, the possible obstacles and the
solution methods. Section 3 demonstrates the evolution of free carrier dynamics and lattice
temperature, the melting depth variation due to laser pulse fluence, revealing also the distinct
melting behavior for different crystal orientations, and possible conditions for amorphization and
its depth minimization. In the final Section 4, we discuss how the obtained results relate to
experimental results and applications in microelectronics. In particular, we discuss the key findings, conclusions and future directions of the presented work.

## 2. Modeling details

The continuum part of the developed combined atomistic-continuum model is based on the approach proposed by van Driel [16], and further developed by Chen et al in Ref. [17]. The corresponding description of the effect of laser-generated free carriers is similar to a so-called Two Temperature model (TTM), which was formulated for metals by Anisimov et al in Ref. [18]. Unlike the TTM, however, where the conservation of energy is due to the balance between free carriers and phonons an addition equation on the carrier's density $N_e$ is required to account for multiphoton and free carrier's absorption process, impact ionization, Auger recombination and the carrier's flux (carrier's diffusion). In this way, the modified TTM for Si material describes the dynamics of free carrier's density as well is referred as nTTM. This model is based on the following system of equations [17]:

$$ 
\left\{
\begin{aligned}
\frac{\partial N_e}{\partial t} &= \frac{\alpha I(z,t)}{hv} + \frac{\beta I^2(z,t)}{2hv} + \delta N_e - \gamma N_e^3 - \nabla \cdot \vec{J}(N_e, T_e, E_g) \\
\frac{\partial U_e(N_e T_e, E_g)}{\partial t} &= \left(\alpha + \alpha_{FCA}(\varepsilon(N_e))\right) I(z,t) + \beta I^2(z,t) \\
-\nabla \cdot \vec{W}(N_e, T_e, E_g) &= \frac{c_{e-h}(N_e, T_e, E_g)}{\tau_e} (T_e - T_l) \\
C_l \frac{\partial T_l}{\partial t} &= \vec{\nabla} \cdot (\kappa \vec{\nabla} T_l) + \frac{c_{e-h}(N_e, T_e, E_g)}{\tau_e} (T_e - T_l)
\end{aligned}
\right.
\quad (1) $$

where $\alpha$, $\beta$, $\delta$ and $\gamma$ and are coefficients of correspondingly 1-, 2-photon absorption, impact ionization and Auger recombination, and $\alpha_{FCA}$ is free carriers absorption coefficient, depending on the dielectric function $\varepsilon$, which in turns depends on the carrier's density $N_e$. The free carriers (electron-hole pairs) and phonons interact via the carrier-phonon coupling, which is defined here as a ratio of free carrier's heat capacity $C_{e-h}$ and the carrier-phonon energy relaxation time $\tau_e$. Noteworthy, the energy flux $W$ in Eq. (1) includes both the carrier's flux $J$ and heat diffusion defined via the conductivities of both electrons and holes:

$$ 
\left\{
\begin{aligned}
\vec{W} &= \{E_g(N_e, T_e) + 2k_b T_e [H_0^1(\eta_e) + H_0^1(\eta_h)]\} \vec{J} - (k_e + k_h) \vec{\nabla} T_e \\
\vec{J} &= -D \left\{
\begin{aligned}
\vec{\nabla} N_e &+ \frac{N_e}{k_b T_e} \left[ H_{-\frac{1}{2}}^1(\eta_e) + H_{-\frac{1}{2}}^1(\eta_h) \right]^{-1} \vec{\nabla} E_g(N_e, T_e) \\
&+ \frac{N_e}{T_e} \left[ \frac{H_0^1(\eta_e) + H_0^1(\eta_h)}{H_{-1/2}^{1/2}(\eta_e) + H_{-1/2}^{1/2}(\eta_h)} - \frac{3}{2} \right] \vec{\nabla} T_e
\end{aligned}
\right\} \\
D &= \frac{k_b T_e}{q_e} \frac{\mu_e \mu_h H_{1/2}^0(\eta_e) H_{1/2}^0(\eta_h)}{\mu_e H_{1/2}^0(\eta_e) + \mu_h H_{1/2}^0(\eta_h)} \left[ H_{-1/2}^{1/2}(\eta_e) + H_{-1/2}^{1/2}(\eta_h) \right]
\end{aligned}
\right.
\quad (2) $$

and $\mu_e$ and $\mu_h$ are the mobilities of electrons and holes respectively. The total energy of free carriers $U_e$ in Eq. (1) is, therefore, given by the sum of both the kinetic and potential energies, which comprise the energy gap $E_g$ (that varies during the laser-irradiated target evolution due to dynamical change of $N_e$ and $T_e$) as well [16] as follows:

$$ U_e(N_e, T_e, E_g) = N_e E_g(N_e, T_e) + \frac{3}{2} N_e k_b T_e \left[ H_{1/2}^{3/2}(\eta_e) + H_{1/2}^{3/2}(\eta_h) \right] \quad (3) $$
where $H_j^i(\eta_c) = F^i(\eta_c)/F^j(\eta_c)$, and $F^j(\eta_c)$ is Fermi-Dirac integrals of order “i” of the carriers with the reduced chemical potential $\eta_e = (E_F-E_c)/k_bT_e$ and $\eta_h = (E_V-E_F)/k_bT_e$ of correspondingly the electrons “e” and holes “h” with respect of their quasi-Fermi level to the levels of conduction and valence bands respectively, and $k$ is Bolzman constant. The parameters of the continuum nTTM model are given in Table 1.

<table><caption>**Table 1:** The parameters of the continuum nTTM model given by the Eq. (1) – (3)</caption><thead><tr><th>Quantity</th><th>value/approximation</th><th>Ref</th></tr></thead><tbody><tr><td>Initial carrier density</td><td>$N_0 = 10^{16} \text{ m}^{-3}$</td><td>[19]</td></tr><tr><td>Initial lattice and carrier temperature</td><td>$T_0 = 300 \text{ K}$</td><td>[20]</td></tr><tr><td>Lattice specific heat</td><td>$C_l = 1.978 \times 10^6 + 3.54 \times 10^2 T_l - 3.68 \times 10^6 / T_l^2$, J/(m<sup>3</sup>K) (T<sub>l</sub> in K)</td><td></td></tr><tr><td>Lattice thermal conductivity</td><td>$k_l = 1.585 \times 10^5 \times T_l^{-1.23}$, W/(m·K) (T<sub>l</sub> in K)</td><td>[21]</td></tr><tr><td>Carrier thermal conductivity</td><td>$k_e = k_h = -3.47 \times 10^{18} + 4.45 \times 10^{16} T_e$, eV/(s m K)</td><td>[22]</td></tr><tr><td>Band gap</td><td>$E_g = 1.170 - 4.73 \times 10^{-4} \times T_l^2 / (T_l + 636 \text{K}) - 1.5 \times 10^{-10} \times N_e^{1/3}$ if $1.170 - 4.73 \times 10^{-4} \times T_l^2 / (T_l + 636 \text{K}) - 1.5 \times 10^{-10} \times N_e^{1/3} \ge 0$ and $0$ otherwise, eV (T<sub>l</sub> in K, $N_2$ in m<sup>-3</sup>)</td><td>[23, 24]</td></tr><tr><td>One-photon absorption</td><td>$\alpha = 1.34 \times 10^5 \exp(T_l/427)$, m<sup>-1</sup> (T<sub>l</sub> in K)</td><td>[25]</td></tr><tr><td>Two-photon absorption</td><td>$\beta = 15 \text{ cm/GW}$</td><td></td></tr><tr><td>Reflectivity</td><td>$R = 0.329 + 5 \times 10^{-5} (T_l - 300 \text{K})$ (T<sub>l</sub> in K)</td><td>[26]</td></tr><tr><td>Auger recombination coefficient</td><td>$\gamma = 3.8 \times 10^{-43}$, m<sup>6</sup>/s</td><td>[27]</td></tr><tr><td>Impact ionization coefficient</td><td>$\delta = 3.6 \times 10^{10} \exp(-1.5 E_g / k_B T_e)$, s<sup>-1</sup></td><td>[28]</td></tr><tr><td>Free-carrier absorption cross section</td><td>$\Theta = 2.91 \times 10^{-22} T_l/300$, m<sup>2</sup> (T<sub>l</sub> in K)</td><td>[29]</td></tr><tr><td>Electron-phonon relaxation time</td><td>$\tau_{ep} = 0.5 \times 10^{-12} (1 + N_e / (2 \times 10^{27}))$, s ($N_e$ in m<sup>-3</sup>)</td><td>[30]</td></tr><tr><td>Electron effective mass</td><td>$m_e^* = 0.36 m_e$</td><td>[31]</td></tr><tr><td>Hole effective mass</td><td>$m_h^* = 0.81 m_e$</td><td>[31]</td></tr><tr><td>Mobility of electrons (taken at 1000K)</td><td>$\mu_e = 0.0085 \text{ m}^2/(\text{V} \cdot \text{s})$</td><td>[32]</td></tr><tr><td>Mobility of holes (taken at 1000K)</td><td>$\mu_h = 0.0019 \text{ m}^2/(\text{V} \cdot \text{s})$</td><td>[32]</td></tr></tbody></table>
While the parameters entering the continuum part of the combined model can be varied and justified either by experimental measurements [16] or theoretical calculations [33,34], the strategy used in the development of the combined model for Si remains the same as the one used for MD-TTM for metals [35]. Namely, we replace the equation. in the system (1) for lattice temperature with MD integrations. As fot the interatomic potential $U(r_{ij}, T_e)$ for Si, there are several classical and modified interatomic potentials that can be used. The classical Stillinger-Webber parametrization [36] allows for an accurate reproduction of the equilibrium melting temperature of 1687K, which is an important point in our research. The classical Tersoff potential [37], although deviating in the melting point and underestimating the heat capacity, is readily used when description of elastic properties becomes important. The more advanced modern interatomic potentials for Si, but a way more complicated from the point parametrization and from the point of computational efficiency, can even include the dependences on the electronic temperature [38] and possibly on the electronic density [39]. Based on the reliable reproduction of thermophysical properties (density and melting temperature) [9] and the simplicity of implementation, we chose the Stillinger-Webber USW interatomic potential function for the calculation of interatomic forces in the MD part of the model as follows:

$$m_i \frac{d\vec{r}_i}{dt^2} = -\sum_{j=1,N} \vec{\nabla} U(\vec{r}_{ij}, T_e) + \xi m_i \vec{V}_i^T \quad (4)$$

where $r_i$ is the position of “i” atom with the mass $m_i$, $N$ is the number of atoms in the discretized cell, $V_i^T$ is thermal velocity of the atom “i”, and the coefficient $\xi$ is determined from the energy conservation law [35]. Because of the additional equation for the carrier’s density, we refer to this combined method as MD-nTTM [9] as schematically presented in Fig. 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the combined MD-nTTM model for silicon (Si). The system of free carriers is described in the continuum (upper part), while the atomistic behavior is modeled using Molecular Dynamics (MD) (bottom part). This integrated approach enables an accurate description of non-equilibrium phase transitions. The free carrier system, modeled by nTTM, is divided into a number of cells, with each cell corresponding to the same volume in the MD portion. The values of $T_i$, $N_i$, $C_i$, and $E_{gi}$ are defined at the centers of the cells whereas the first derivatives and the fluxes $J_{i+1/2}$ and $W_{i+1/2}$ are determined at the cells’ interfaces. The nTTM framework captures the processes of laser light absorption, free carrier generation, and their subsequent diffusion and thermal diffusion. By synchronizing the time steps of nTTM and MD, where $\Delta t_{MD} = n\Delta t_{nTTM}$, the energy transfer from hot carriers to the corresponding MD cell is equally distributed over the atoms constituting that MD cell. For thick targets, the application of Non-Reflective Boundary (NRB) conditions [40] allows one to limit the MD application to the top of the irradiated target only, leaving its bottom part of the computational cell described completely via TTM model</figcaption></figure>
[35,41]. The utilization of NRB allows to avoid pressure wave reflection artifacts from the MD boundaries.

It is worth mentioning that unlike a similar model for Si published in Ref. [42], where non-degenerate model based on the Boltzmann distribution of the electron-hole pairs was assumed (indicating negligible excitation), the complete description of the laser-generated free carriers dynamics via Fermi-Dirac functions is applied in the presented model. This important modification results in the better applicability of the MD-nTTM model for investigation of laser melting mechanism of Si at the fluences up to the ablation threshold and above.

To solve the proposed system (Eq. 1), the finite difference grid mesh is used, as sketched in Fig. 1. The sample is divided into cells according to the scheme, and the thermodynamic parameters are calculated in each cell. The spatial derivatives of $n$, $T_e$, $T_a$, $J$, $W$, and $E_g$ at the interior points are approximated with the central difference, and those at the boundary are evaluated with the first-order approximation. The explicit algorithm for integration of the finite difference method, where all new time step values are explicitly defined via the current time values is straightforward. In order to advance each time step, any cell only needs the characteristics of its closest neighbors. So, the explicit algorithm for the finite difference method is well parallelized and can be computationally efficient on the one hand. However, in the case of Si, we deal with non-constant carriers' density. Given by the Neuman stability criterium, the choice of the time step $\Delta t$ in the solution of nTTM approach with the explicit finite-differences scheme must satisfy the condition $\Delta t \le \Delta z^2/(2D_N)$, where $\Delta z$ is spatial discretization of the system along the laser pulse propagation, Fig. 1, and $D_N$ is the diffusion coefficient of free carriers and $D_N \sim k_e/C_{e-h}$ and $C_{e-h}$ is heat capacity of free carriers taken as $C_{e-h} = dU_{e-h}/dT_e$. Because of the number of free carriers in the major part of the irradiated target is negligible (except the depth with the laser light absorption), the value of $D_N$ is numerically uncertain. Thus, the application of the explicit scheme of the finite differences method in the solution of the system given by Eq. 1 results in an extremely shallow time step (~$10^{-25} \div 10^{-22}$ s), depending on spatial discretization (1nm – 10 nm) dictated by the modern technological requirements. As the MD time step is normally as small as ~$10^{-15}$ s, the combined MD-nTTM model in this case becomes either very time-consuming or poor on resolution [43]. To overcome this obstacle, a more advanced implicit algorithm based on the Crank-Nicolson scheme [44] was utilized that allowed us to raise the nTTM time step by few orders of magnitude reaching the range of $10^{-17} \div 10^{-16}$ s [45] making the suggested MD-nTTM approach much more efficient to study strongly non-equilibrium laser-induced phase transitions processes in Si. In the implicit algorithm, however, the new time step values are defined via the same new time step values, which leads to the sequential solution of each cell via three diagonal matrix equations (tensor of second order). And because of the major sequential part in this algorithm, the implicit scheme of the Finite Differences (FD) method can be barely parallelized for multiprocessing regime. Thus, despite the efficient integration from the point of advancement in time, the application of standard Crank-Nicolson scheme to solve the hstem given by Eq. 1 frequently leads to the limited material volume for numerical processing and in one dimension (1D) only in the direction along the laser pulse propagation. While this problem does not concern us in the present study, the necessity of implementation of Crank-Nicolson scheme in three dimensions (3D) and in multiprocessing regime for investigation of the Si surface functionalization/nanostructuring mechanisms becomes obvious. Likely, an elegant modification of Crank-Nicolson scheme in 3D was recently suggested in Ref. [46] and will be utilized in our future studies.
## 3. Results and Discussions

### 3.1 Free Carrier and Lattice Dynamics

To clearly demonstrate the applicability and the accuracy of the nTTM model implementing Crank-Nicolson algorithm for description of free carriers' temperature and density dynamics, we perform the initial modeling for the ultra-short laser pulse with temporal width of 100fs at the incident fluence of 0.2 J/cm² interacting with 2000nm free standing Si film. To describe the reflectivity function and the absorption by free carriers we used Drude model as it was suggested in Ref. [33] assuming the wavelength of 800nm. In this case we use laser intensity is given by Gaussian distribution as follows:

$$I_0(t) = F(1 - R(\varepsilon))^{\frac{1}{\tau}} \sqrt{\frac{\sigma}{\pi}} \exp\left(-\sigma \frac{(t-t_0)^2}{\tau^2}\right) \quad (5)$$

where F is the average incident fluence, τ is pulse duration, the pulse is centered at $t_0 = 2.5\tau$ and $\sigma = 4\ln2$. Assuming that the laser spot size is much larger than the size of the computational domain, we can handle with the combined MD-nTTM approach. Here and below we apply 1D heat conduction model in the solution of Eq. 1 and 1D (along the direction of beam propagation in Z) laser light absorption model as well:

$$\frac{dI(z,t)}{dz} = -\left(\alpha + \alpha_{FCA}(\varepsilon(N_e))\right) I(z,t) - \beta I^2(z,t) \quad (6)$$

To have a possibility of referring our simulation results to the experimental data, we need furthermore an accurate description of the reflectivity function $R(\varepsilon)$, which depends on the level of excitation via the electronic temperature, affecting therefore the complex dielectric function $\varepsilon$:

$$\varepsilon(\omega_L) = \varepsilon_r(\omega_L) - \frac{N_e e^2}{\varepsilon_0 \omega_L^2} \left[ \frac{1}{m_{e,cond}^*(1+i\frac{\nu_e}{\omega_L})} + \frac{1}{m_{h,cond}^*(1+i\frac{\nu_h}{\omega_L})} \right] \quad (7)$$

where $\varepsilon_r$ is the intrinsic dielectric constant, $\omega_L$ is the laser frequency, and $m_{e,cond}^*$ is the conductivity effective mass defined to reproduce the electrical conductivity and the susceptibility [47]. The Drude collision frequency of electrons and holes, respectively, is denoted as $\nu_c$. By supposing that the collision frequencies of electrons and holes are equal and by introducing a common optical effective mass $1/m_{opt}^* = 1/m_e^* + 1/m_h^*$, Eq. 5 gives an ordinary Drude model for dielectric function [48], where complex refractive index $n^* = \varepsilon^{1/2}$ and the reflectivity function R and free carriers' absorption $\alpha_{FCA}$ are given correspondingly:

$$R = \frac{[n^*-1]^2}{[n^*+1]^2} \quad \text{and} \quad \alpha_{FCA} = \frac{Im(n^*)\omega_L}{c_0} \quad (8)$$

where $c_0$ is the speed of light in vacuum. The utilization of Drude model is justified here due to strong variation of the complex dielectric function during the pulse evolution. This makes the optical penetration depth ~ 1µm in the pre-pulse regime, whereas it shrinks to just few tens of nm when the pulse reaches its maximum and free carriers' density rises by several orders of magnitude.

For the above laser pulse parameters, Fig. 2 shows the calculated electron temperature, lattice temperature, and free carriers' density dynamics obtained by using the nTTM model Eq. (1-3) together with Drude model given by Eq. (5-8).
<figure><img src="image_4.png" /></figure>

**Figure 2.** The dynamics of free carrier temperature, lattice temperature, and free carrier density are illustrated for the excitation of a 2000 nm free-standing silicon (Si) film using a 100 fs laser pulse with a wavelength of 800 nm at an incident fluence of 0.2 J/cm². Please note that the total energy conservation and the intensity of the laser irradiation are not in scale.

With the realized Crank-Nicolson algorithm for the solution of the system given by Eq. 1, the computational cell was spatially discretized into 400 cells (cell size is 5nm), which gives a high enough precision on one hand and the efficiency of calculation (low computational cost) on the other hand. The accuracy of the algorithm implementation is guaranteed by the total energy conservation (not in scale), Fig. 2. The initial jump of the electronic temperature is explained by a negligible density of free carriers at the beginning of the pulse resulting in their extremely low heat capacity. The following deviation of the electronic temperature from what is normally observed in traditional TTM is due to the dynamically changing chemical potential of free carriers and the energy gap as well. Then, the calculated lattice temperature and free carriers' dynamics, referred to the shape of the pulse intensity (not in scale) are reasonable. The realized Crank-Nicolson scheme for the solution of nTTM model is therefore acceptable to be combined with the classical MD method

**3.2 Melting Depth Dependency on Laser Fluence**

The relevant simplicity and efficiency of the nTTM numerical method looks frequently
attractive for the interpretation of many experimental findings [49]. The standing alone nTTM
model, however, does not include the kinetics of melting mechanism of Si, and does not account
on the peculiarities of Si crystal structure under extreme conditions realized in the material upon
the ultrashort laser pulse irradiation. For this purpose, we perform MD-nTTM study of the melting
depth variations in Si due to 100fs pulse at 800nm for a range of fluences. Because laser sport size
is assumed to be large, we limit the MD simulation to its center. By accounting for the absence of
the temperature and pressure gradients in lateral directions until microseconds, we apply periodic
boundary conditions in directions X and Y.

The initial sample with <001> crystal orientation and dimensions of 5.5 x 5.5 x 1000 nm was prepared and equilibrated at 300K. The sample thickness of 1000 nm was chosen based on the efficiency of calculations on one hand, and to exclude the effect of sample thickness on the obtained result on the other one. For instance, we obtain the same melting threshold for 1000nm
sample as it was detected in large modeling with 2000 nm sample well. Finally, we run the MD-nTTM simulations for the range of fluences beginning with the melting threshold up to the ablation and above: 0.3 – 1.5 J/cm². The results obtained with MD-nTTM model are compared with simulations given by standalone nTTM for the same conditions and shown in Fig.3 (a).

The melting depth in our simulations as determined with the help of Central Symmetry Parameters (SCP) [50], that was successfully utilized in our previous works as well as. The CSP was shown as an accurate method for identification of rapid laser-induced melting and recrystallization processes in solids. There are also other methods of liquid phase identification. Lindeman criterium [51] however, suffers from low precision in the case of ultrafast non-equilibrium phase transitions. Alternatively, the Local Order Parameter [35] is basis orientation dependent, which makes it useless for the identification of recrystallized solid via homogeneous mechanism, when the nuclei of crystal phase can be generated in bulk with random crystal orientation. In this sense the choice of CSP is justified by its high sensitivity and not being dependent on the basis vectors' choice.

Interestingly, the melting depth given by the MD-nTTM model reveals quite a complex behavior as a function of the incident laser fluence. First, the calculated melting threshold is equal to 0.3 J/cm² and it agrees with number of the previous experimental measurements obtained under similar conditions [52,53,54]. Although our calculation results are closer to the experimental measurements obtained by the group of Bonse [53], some of the results such as melting depth and melting threshold vs. crystal orientation obtained by the group of Grojo [52] are also reproduced.

<figure><img src="image_5.png" /></figure>

**Figure 3.** (a) Si melting depth calculated for a 100 fs laser pulse as a function of incident fluence at a wavelength of 800 nm is presented, as determined by the MD-nTTM and nTTM models. The presented results distinguish between the depth of pre-surface melting occurring due to the heterogeneous mechanism and the depth obtained due to the homogeneous melting mechanism. The latter one is calculated by summation of the melted volume in bulk. (b) A zoomed-in view of the melting depth, measured as the surface melting and the depth associated with the total melted volume as given by the MD-nTTM model at low fluences is also shown.

We also note that the resulting melting depth at low fluence regime is clearly attributed by both heterogeneous and homogeneous melting mechanisms. This can be seen in Fig. 4, where the final atomic configuration reveals several melted volumes in bulk, including the surface melting.
On one hand, material melting occurring due to the classical heterogeneous mechanisms starts from free surface and propagates deeper with a well-distinguished solid-liquid interface. In the other hand, material melting tacking place due to homogeneous mechanism is observed in the bulk due to high degree of overheating of the crystal. In the latter case, the material losses its stability and spontaneous solid-liquid transition takes place. We therefore calculate two depths of melting provided by these two different mechanisms. One of those values is the position of the solid layer closest to the surface, which can be detected in X-ray experiments, whereas the other one is accumulated melted volume recalculated into the corresponding melting depth value under the established thermodynamical conditions. Therefore, Fig. 3 (a) shows both the maximum melting depth due to heterogeneous phase transition mechanism and the green squares indicate the accumulated depth of melting (based on the summed melted volume) due to homogeneous mechanism as well, Fig 3 (b). However, as laser fluence increases, the difference between both measurements disappears. This can be explained by the fact that under laser irradiation strong enough, the relatively slow heterogeneous mechanism contributes only slightly to the final state as compared to homogeneous melting mechanism that takes a massive character.

Finally, at the fluence above 0.7 J/cm² the melting depth starts scaling rather linearly, due to the onset of the spallation mechanism of target damage, with removal of the material's surface volume, which possesses high thermal energy. This closely corresponds to the experimental measurements of Bonse group [55]. Here we refer to the spallation mechanism of the ablation process, when as compared to the common thermal damage the material undergoes to rather mechanical damage due to extraction of large chunks and even presurface layers.

Additionally, one of the frequently utilized methods in the estimation of the melting depth when interpreting the experimental results is given by a simple expression as follows

$$d_{melt} = A\ln(F/F_{th}) + B, \qquad (9)$$

where F and Fth are the incident and thresholds fluences correspondingly. After the analysis of the obtained melting depth in our simulations, it was found that Eq. 9 is applicable in out MD-nTTM modeling results as well. The corresponding fitting coefficients were determined to be A = 600nm and B = -220nm.

On the other hand, the melting depth obtained with a standalone nTTM model under the same conditions reveals much larger values (100 – 200 nm) at low fluences and even lower at high energies. The melting threshold given by nTTM was found at the fluence of 0.16 J/cm², which completely reproduces the result obtained in work by Ramer [33]. We note, however, that these results are almost twice as low as the one obtained in MD-nTTM simulations. Such a large difference can be attributed to the absence of melting kinetics in the nTTM model that is, however, naturally included in the MD method. Also, the effect of pressure and the latent heat of melting are missing in the standard nTTM. This results in underestimation of the melting threshold provided by the continuum methods such as nTTM. At high fluence regime, however, the melting depth given by MD-nTTM approach becomes even higher. This is the result of the onset of the spallation regime of the material removal, which is also absent from the description by nTTM alone.

### 3.3 Melting Dynamics and the Role of Crystal Orientation

The next part of our study is devoted to the effect of crystal orientation on the melting mechanism and the depth of melting. There are many experimental measurements indicating the smaller thresholds of melting, ablation, and amorphization for Si crystals with <111> orientation [56]. To clarify this discrepancy and to investigate the crystal orientation dependent melting mechanism, we prepared two samples of the same dimensions, but with two different crystal
orientations <001> and <111>. Then, both samples were irradiated with 100 fs laser pulse at 800 nm at the same fluence of 0.3 J/cm², that was above identified as the melting threshold $F_{th}$. The modeling results can be seen in Fig. 4 revealing the melting dynamics of both samples as a sequence of the atomic configurations taken at the same times. The atoms are color-coded based on the Central Symmetry Parameter (CSP), which allows us to quantitatively distinguish between atoms in a solid environment and those surrounded by the liquid phase. A sharp jump is observed at the empirically determined threshold value of CSP = 0.97, marking the phase transition. Additionally, the rarefaction of the atomic configuration due to the negative volume change during melting results in some atoms having fewer neighbors, leading to lower local symmetry. This is reflected by a blue color.

For both <001> and <111> samples, the relative contribution of the melting due to heterogeneous and homogeneous mechanism results in accordingly roughly 5 nm of surface melting and the melted bulk volume corresponding to 25 nm. For the complex analysis of the melting mechanism, the temperature-pressure dynamics are also plotted for the same time delays as a function of depth. Thus, we can identify the drop of temperature upon the phase transition due to part of the thermal energy transferring to the latent heat of melting. Also, the phase transition of Si can be correlated to the drop of pressure values.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Two sets of sequences of atomic configurations illustrating the differences in the laser-induced melting process for Si samples with: (a) <001> and (b) <111> crystal orientations. Both samples were irradiated with a 100 fs pulse at 800 nm, using the same incident fluence of 0.3 J/cm². Snapshots are presented at identical times (30ps, 50ps, and 100ps) for both cases, with atoms color-coded according to the Central Symmetry Parameter (CSP), enabling the distinction between solid (CSP > 0.97) and liquid (CSP < 0.97) phases. The empirically found threshold value of 0.97 for distinguishing solid and liquid phase of Si was confirmed by a careful analysis. The corresponding profiles of the lattice temperature and pressure are added for a better quantitative analysis accounting on the temperature-pressure interplay.</figcaption></figure>

Despite a similar melting dynamic of both samples, the result obtained with MD-nTTM model reveals that the sample with <001> orientation contained about 25% more material in liquid
state than the one with <111> orientation. This simulation result, although partially controversial with some experimental data [52], can be theoretically supported. Namely, as it is well-known, fast heating of the irradiated solid target due to ultrashort laser pulse induces a fast increase in the internal pressure as well [57], as if the material would be heated under the constant volume. This process occurs when the rate of mechanical sample relaxation is slower as compared to the rate of its heating. The subsequent relaxation of the laser-induced stresses inside the material, however, can vary depending on the irradiation condition and, especially, on the crystal structure. It was shown in Ref. [9] that the effect of pressure on melting dynamics of Si is quite different from that of metals. For metals, mostly exhibiting fcc crystal structures, the increase of pressure leads to the reinforcement of crystal stability against the homogeneous melting. The phase transition itself results in an increase of material volume by 3-5% / [35]. However, in case of semiconductor materials, exhibiting an open diamond crystal structure, we have an opposite effect. The increase in pressure in the solid phase of Si weakens its stability, and the material loses its volume upon the phase transition by ~7% [9]. Therefore, upon reaching some critical values of the initially rising internal thermal stresses, a rapid solid-liquid phase transition takes place and an abrupt drop in the internal pressure due to the negative volume of melting [58]. These values depend on the local temperature.

We note that the interplay between temperature and pressure due to laser heating leads to conditions favorable for the onset of homogeneous melting at a certain distance from the surface, within a spatial range of 10 nm to 40 nm (Fig. 4). Within this range, both temperature and pressure are sufficiently high to initiate the formation of liquid nuclei. Nucleation occurs randomly within this region, and due to the negative volume change upon melting, the solid-liquid phase transition causes immediate mechanical relaxation in the neighboring crystalline regions. This is reflected in the rapid growth of liquid nuclei, with a relative stabilization of the remaining crystalline parts. Continuous heat flux from the surface further drives the gradual expansion of the liquid phase through a slow heterogeneous melting mechanism. However, due to the high computational cost of the simulations and the limited lateral size of the target, periodic boundary conditions lead to the merging of liquid nuclei after 100 ps, resulting in a stripe-like liquid-solid final configuration. For a larger sample, a distribution of liquid volumes with rough spherical shapes and partial coalescence is expected.

In our simulations, as can be seen in Fig.4, the final configuration of the <111> crystal is expanded a bit more than the one with <001> orientation. This result can be due to a higher value of the speed of unloading pressure wave <111> Si. To verify this, we performed two identical MD simulations, in which the densification of atomic system by 5% was created at the edge of both samples. Following the pressure wave propagation, it was found that the speed of the pressure wave is 8.1 nm/ps in <001> direction, whereas it is 8.74 nm/ps in <111>, that is 8% greater. This can be also seen in Fig. 4, where at the time of 50ps the pressure wave tail is still visible within the frame scale for the crystal <001>, whereas it is gone for the crystal <111>. As a result, a faster mechanical relaxation of the <111> crystal leaves it in a preferable condition for supporting it in the solid phase. The temperature drop at the surface is associated with the part of thermal energy going to the latent heat of melting upon the phase transition at the surface.

For a quantitative analysis of the melting mechanisms, the temperature and pressure dynamics are plotted as a function of depth for both samples at the same time delays in Fig. 4. This allows us to observe the temperature drop during the phase transition, which occurs as part of the thermal energy is transferred to the latent heat of melting. Additionally, the phase transition of Si can be correlated with the drop in pressure values. However, we note that the temperature drop at
the surface is more pronounced, as part of the thermal energy loss to the latent heat of melting occurs under nearly zero pressure due to the proximity to the surface, resulting in an immediate mechanical relaxation of the melted volume. In contrast, for deeper melted regions, the solid-liquid phase transition occurs under constant volume conditions, which, along with the negative volume change upon melting, compensates for the thermal energy loss through an input of potential energy.

When comparing our modeling results for both samples with the experimental findings, we keep in mind that during the simulation process we used the same Drude-Lorenz model with the same parameterization as in Ref. [33]. We note, that there are studies considering the optical properties of Si to be dependent on the crystal orientation [59]. This fact is actually not accounted for in our model due to lack of the published data. Nevertheless, the difference in the absorption process, depending on the crystal orientation, can potentially play a more important role leading to a smaller melting threshold for the Si crystal with <111> direction, as was detected in the experiments [52,53,54].

To provide a more quantitative view of the melting process difference between different Si crystal orientation, we present the melting dynamics using CSP variations in the contour plots shown in Fig. 5. This allows us to distinguish and quantify the contributions of two melting mechanisms: heterogeneous and homogeneous. Heterogeneous melting, a classical process, begins at the free surface and progresses into the bulk, though it is limited by the propagation speed. It was also shown in our previous theoretical simulations with MD-TTM approach that the maximum speed of the solid-liquid front propagation due to heterogeneous mechanism is significantly overestimated in experimental data interpretation, when assumed to be limited by the speed of sound in the material [60]. In contrast, the homogeneous melting occurs when the material is overheated (~25%) above the melting point, making the surface-driven melting impossible and causing rapid bulk melting. In Fig. 5, the yellow and orange regions readily represent the melting dynamics and the resulting melting depth, consistent with Figure 4. However, fluctuations are observed before 20 ps, along with greenish areas that persist until the end of the simulation, and thin yellow-green regions appear around 50 ps below the main melting depth.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Calculated dynamics of the CSP parameter for the irradiated Si target with <001> (a) and <111> (b) crystal orientations are shown as a contour plot for a 100fs pulse at a wavelength of 800 nm and an incident fluence corresponding to the melting threshold interacting with 1000 nm samples. The CSP parameter is represented by color with the empirically chosen value of 0.97 as a threshold value for distinguishing the solid (CSP > 0.97) and liquid (CSP < 0.97). The greater</figcaption></figure>
expansion of the irradiated <111> sample due to faster speed of sound is also visible in (b) as a
ledge below zero position (depth) in green.

From the above discussions we can draw an important conclusion regarding the kinetics of laser-induced non-equilibrium melting in Si. Namely, the noticeable contribution of both heterogeneous and homogeneous mechanism into the melting process is observed at the fluences of just above the melting threshold. In this case, the relatively slow heating rate leads to a non-significant rise in internal stresses. Thus, the homogeneous mechanism of liquid phase nucleation due to weakening of the crystal stability is limited and clearly attributed with the classical heterogeneous mechanism of the surface melting, Fig. 4. This fact explains the different values of the melting depth in Fig. 3 (b) defined as the surface melting and the accumulated melting depth due to homogeneous mechanism of the phase transition in bulk (black dots and green squares correspondingly). However, at higher fluences, the rate of heating becomes faster, leading to stronger accumulated internal stresses, that causes Si crystal losing its stability over a shorter time period and the contribution of the homogeneous melting mechanism is increasing. Eventually, already at the fluence of 0.5 J/cm² the contribution of the heterogeneous melting becomes negligible, whereas the homogeneous melting takes a massive. This situation is reflected in Fig. 3 (a-b) when surface melting depth and accumulated melting depth given by black dots and green squares correspondingly become indistinguishable. Thus, the homogeneous melting mechanism in Si due to an ultrashort laser pulse is basically the only melting mechanism for laser fluence of 0.5 J/cm² and above. It should be emphasized that such a non-obvious result is very different from that obtained in laser melting of metals with fcc crystal structure [60,61,62,63].

The emergence of localized melted regions in the simulations of ultrashort laser irradiation of silicon raises important questions about the underlying mechanisms governing non-uniform phase transitions. The above applied model, however, could not consider explicitly all the factors that may play an important role in a real experiment during laser-driven melting process, so that the obtained here results may have a semi-quantitative character. Thus, in addition to free carrier dynamics and ultra-fast energy deposition, few more factors may contribute to this phenomenon particularly in multi-pulse irradiation experiments, including pre-existing defects and dislocations, as well as localized variations in electronic and thermal transport properties. Importantly, ultra-fast energy deposition definitely plays an crucial role. The non-equilibrium conditions promote rapid heating in localized regions where electron-lattice energy exchange occurs. This effect is particularly pronounced in materials with complex band structures, such as silicon, where free carrier dynamics and carrier localization can influence the spatial distribution of absorbed energy. Additionally, structural defects and dislocations within the crystalline lattice can act as preferential sites for melting initiation. These defects locally alter the bonding environment, reducing the effective melting threshold compared to pristine crystalline regions. Dislocations, in particular, serve as stress concentration points where atomic vibrations are enhanced, potentially leading to local temperature spikes upon laser irradiation. The presence of impurity atoms or vacancies can, furthermore, modify the local electronic structure, affecting carrier mobility and electron-phonon coupling efficiency. The interaction between laser-induced free carriers and defect states could further enhance localized heating, triggering spatially inhomogeneous melting patterns. Finally, the interplay between free carrier transport and thermal diffusion also plays a critical role in the formation of melted spots. Silicon exhibits a strong anisotropy in both electronic and phononic transport properties, leading to direction-dependent energy redistribution. Given the rapid heating timescales, non-uniform carrier diffusion and recombination processes may contribute to spatial variations in energy density. Moreover, if carrier recombination occurs at different rates across the
irradiated region, some areas may reach the melting threshold faster than others, reinforcing heterogeneous melting.

Nevertheless, the applied above model accounted on the effect of laser-generated free carriers' dynamics via their density and temperature, that in the case of ultrashort laser pulses plays a determinant role. The model also captured the main driving mechanisms of the melting process, and it reproduced the melting threshold, matching its value to many experimental measurements. This is an important step forward towards the detailed understanding of the laser-induced processes in the solid and generation of Si surfaces with predesigned electrical and optical properties.

### 3.4 Laser-induced amorphization of Si

The amorphous phase of Si in the presurface region is frequently registered upon the ultrashort laser excitation [52] in the experiments and technologically attracts a high interest due to its unique properties [6]. The electrical and optical properties of the amorphous Si layer depend on its thickness, which in turn is a function of the laser irradiation parameters. So, a better understanding of the mechanism of the amorphous layer formation is of particular interest nowadays, especially in microelectronics. The detailed understanding of the mechanism of the amorphization process due to an ultrashort laser pulse also helps in generating the well reproducible structures with specific parameters under controllable condition.

In this study, we seek to identify the necessary cooling conditions for achieving amorphization and establish their correspondence with optimal laser-Si interaction regimes in experiments [36]. To achieve this, we prepared an initial sample at a temperature 25% above the equilibrium melting point, allowing it to undergo heterogeneous melting from the surface. We then applied a range of cooling rates from 1 to 50 K/ps, systematically observing the recrystallization process. The amorphous state was identified when the propagation of the solid-liquid interface ceased for a prolonged period (250 ps and more), even if the sample is cooled down to approximately 1000 K, which is well below the melting point of 1687 K Given that CSP does not inherently distinguish between liquid and amorphous states, the amorphous configuration was determined based on the stalling of the solid-liquid interface motion and the final atomic structure, as depicted in Fig. 6.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The obtained in simulations presurface amorphous state of Si was achieved through a cooling process of the melted region with a certain rate. Atoms are colored by the CSP, indicating that amorphization occurred when prolonged immobility of the solid-liquid interface takes place, even after the sample temperature dropped to 1000 K, well below the melting point of 1687 K.</figcaption></figure>

The process of recrystallization is governed by the mobility of atoms and is related to the diffusion coefficient, and, hence, on a local temperature. Importantly, Gibbs free energy of the solid state should be lower than it is for the liquid state under the current conditions. As a result, the crystalline phase would be more preferable from the energetical point of view. Meanwhile, the local temperature should be high enough, so that the mobility of the atoms would let them take the equilibrium positions corresponding to the crystal structure. If the cooling rate is, however, high enough, the mobility of atoms can rapidly decrease, leading to the situation when the atoms have not enough time for the reconstruction of the crystal structure. If this happens due to laser cooling, the formation of the amorphous phase is actually observed. When considering the amorphization process, however, the CSP used above for quantitative description of the melting dynamics becomes useless for the identification of current phase. Namely, because any order is missing in bonding construction, the atoms in amorphous state, highlighted with CSP, are indistinguishable from those in the ordinary liquid phase. We could quantitatively detect the liquid-amorphous phase transition by considering the kinetics of the solid-liquid interface propagation. When a transition to the amorphous phase occurs, the interface becomes immovable, even though the local temperature is well below the equilibrium melting point. In addition, we could also compare the density of the crystalline phase and the phase that we assume amorphous. The straightforward calculations under the conditions depicted in Fig. 6 yield 4.856*10²² at/cm³ for the solid and 5.179*10²² at/cm³ for the state we suspect to be amorphous correspondingly. As one can see, the obtained densities differ by -6.25%, which is a bit less than the volume of melting (-7%) due to a lower temperature of the obtained state. So, the immovable interface together with the characteristic density difference helps us to identify this phase as amorphous.
Molecular dynamics (MD) simulations indicate that a cooling rate of approximately 30 K/ps is required to achieve an amorphous Si layer. However, when performing nTTM simulations to replicate these cooling conditions, the maximum cooling rate found near the melting threshold was only about 2.5 K/ps, which is an order of magnitude lower than the one required by MD results. This significant discrepancy, on one hand, highlights the limitations of continuum models, such as nTTM, that are missing the kinetics of recrystallization. But on the other hand, the choice of interatomic potential plays a crucial role in the modeling of Si amorphization dynamics.

Herein we have employed the standard parametrization of the Stillinger-Weber (SW) potential [64], which reliably reproduces the thermodynamic properties of crystalline Si. This choice enables the direct comparison of simulation results with experimental measurements. An alternative parametrization of the SW potential [65], specifically optimized for describing amorphous phases, has demonstrated appropriate accuracy in modeling amorphous Si. However, it fails to reproduce the correct melting and point, thereby hampering the direct interpretation of the simulation results and their comparison with the experimental data. Furthermore, at higher fluences when non-thermal effects become dominant, traditional interatomic potentials may no longer be used. In such cases, advanced interatomic potentials incorporating the electronic temperature and density effects [66] or the potentials based on machine-learning technique [67] could offer a more accurate description of phase transitions under extreme conditions. These newly developed potentials may enable better predictions of ultrafast laser-induced structural transformations, making them a promising tool for future studies on Si amorphization. For the considered here conditions (moderate excitations above the melting threshold), a good compromise would be probably the modified Tersoff potential [68]. In fact, it was claimed to be more accurate for description of the solid-liquid phase transition as compared to the standard form of Tersoff potential parametrization [69]. Although they both underestimate the latent heat of melting, a significant improvement of the modified Tersoff interatomic potential in description of the laser-induced phase transition was recently demonstrated [70]. Furthermore, the modified Tersoff was also recently shown as an appropriate interatomic potential for investigation of the amorphous phase of Si [71]. Unlike the modified SW potential, the modified Tesoff interatomic function accurately reproduces the melting point and, therefore, can be applied in the simulations targeting the direct comparison of the simulation results with the experimental data. Therefore, we are also working on the implementation of such modified Tersoff interatomic potential in our simulation code for future detailed investigations of the amorphization kinetics of Si due to ultrashort laser pulses. These results, as well as the comparison between the results obtained by using different potentials will be presented elsewhere.

Once the appropriate interatomic potential is selected, no further assumptions are made in MD regarding the nature of the laser-induced processes under study. However, when combining the atomic description of non-equilibrium phase transitions with the continuum description of free carrier dynamics, the parameters involved in such models—such as those in TTM [18] and nTTM [17]—are derived from different models. Depending on the specific experimental geometry, irradiation parameters, and excitation levels, it is important to assess whether the chosen parametrization of the nTTM approach provides the required accuracy. For example, when comparing the modeling results with experimental measurements, the model for reflectivity becomes crucial. The selected Drude-Lorentz model is appropriate for optical wavelengths, but modifications are necessary for infrared (IR) and longer wavelengths. Finally, while progress is being made to overcome these limitations, the challenge of running simulations in a multiprocessing regime for modeling laser-induced processes on scales directly comparable to
experimental conditions remains the most difficult task, primarily due to the complexity of
implementing the Crank-Nicolson scheme in 3D space.

Based on the above results on melting threshold, melting dynamics, and on the achievement
of the fastest cooling rate upon the laser pulse irradiation of Si target, we can now state that the
required minimization of thickness of the surface amorphous layer of Si can be obtained if one
works at the fluences just above the melting threshold. This is an important conclusion of our
research, which should be verified experimentally. At present, the predicted melting threshold
value from the MD-nTTM simulations aligns with the measurements by Bonse [53] and partially
reproduces the data obtained by Grojo [52]. Further verification of the MD-nTTM model's
applicability to ultrashort laser-silicon interaction experiments can be conducted by comparing the
results for the obtained melting depth (Fig. 3). However, the corresponding TEM experiments
would require a complex experimental setup and significant effort to achieve the required
resolution. From this perspective, the detection of the amorphous phase in experiments and the
investigation of amorphization conditions using MD-nTTM simulations appear promising and will
be the focus of our future research.

The findings of this study underscore the critical role of cooling dynamics in laser-induced
amorphization and highlight the necessity of employing appropriate computational models. The
discrepancy between MD and nTTM predictions emphasizes the need for refined approaches that
integrate atomistic and electronic effects in laser-matter interaction simulations. These insights
pave the way for improved control over amorphous layer formation in microelectronic and
optoelectronic applications, ultimately enhancing the precision and reproducibility of laser-based
material processing techniques.

4. Summary and Conclusions

This study presents a comprehensive investigation of ultrafast laser-induced phase
transitions in silicon (Si) using a hybrid atomistic-continuum model (nTTM-MD). By
incorporating a refined Drude-Lorentz approach for the description of free carrier absorption
process, the model accurately captures the interplay between electronic excitation, energy transfer,
and structural transformations under femtosecond laser irradiation. The simulations employ a laser
pulse of 100 fs duration and 800 nm wavelength, with fluences ranging from 0.3 to 1.5 J/cm²,
reflecting experimentally relevant conditions for ultrafast laser processing.

Our results demonstrate that precise modeling of free carrier dynamics is essential for
predicting phase transitions in Si. The hybrid nTTM-MD model determines the melting threshold
to be approximately 0.3 J/cm², in agreement with experimental data, whereas conventional
continuum models such as the standalone nTTM significantly underestimate this value at 0.16
J/cm². This finding underscores the importance of atomic-scale simulations in achieving
quantitative agreement with experiments, particularly for applications in microelectronics and
optoelectronics, where precise control over phase transitions is critical for device performance and
reliability.

The comparison between the MD-nTTM and nTTM results emphasize the benefit of
atomistic modeling in capturing the fundamental mechanisms governing ultrafast laser-induced
phase transitions in silicon. The substantial discrepancy in melting thresholds (0.3 J/cm² vs. 0.16
J/cm²) and melting depths at low fluences highlights the limitations of conventional continuum-
based approaches, which fail to account for key physical effects such as melting kinetics, pressure
relaxation, and latent heat absorption. The absence of these factors in nTTM leads to a significant
underestimation of the melting threshold and an overestimation of melting depth at lower fluences.
Conversely, at high fluences, the MD-nTTM model reveals a transition to a spallation-driven regime, where the onset of material ejection results in an accelerated scaling of melting depth. This effect, which is absent in nTTM, aligns closely with experimental observations and further validates the necessity of hybrid atomistic-continuum modeling for accurately predicting material response under extreme conditions. The extracted logarithmic fit, provides a robust parameterization of the melting behavior, demonstrating that MD-nTTM not only enhances predictive accuracy but also enables systematic optimization of laser processing conditions. These insights are crucial for applications requiring precise control over phase transitions, such as laser annealing, surface structuring, and semiconductor fabrication, where the balance between thermal and mechanical effects determines the final material state.

Our analysis elucidates the competition between heterogeneous and homogeneous melting mechanisms. At fluences near the melting threshold, heterogeneous melting has a noticeable contribution along with the homogeneous mechanism, resulting in presurface melting depths of ~6 nm, with an additional ~24 nm contribution from homogeneous nucleation. As fluence increases above 0.5 J/cm², homogeneous melting becomes the primary and the only mechanism, leading to a linear scaling of melting depth with fluence. This transition is crucial for laser-based material processing, including nanostructuring, thin-film fabrication, and selective ablation, where depth control is paramount for industrial applications.

The study further reveals a strong dependence of the melting dynamics on crystal orientation. Under an incident fluence of 0.3 J/cm², Si samples oriented along <001> exhibit ~25% greater resulting melted volume than <111>-oriented samples, an effect attributed to faster pressure relaxation along the <111> direction. This anisotropy has direct implications for laser-based surface modifications, where controlled energy deposition is essential for achieving uniform material processing in applications such as photonic circuits and advanced semiconductor manufacturing.

Moreover, our simulations provide insight into the conditions required for laser-induced amorphization of Si, a technologically significant phenomenon in phase-change memory and semiconductor device engineering. The results indicate that achieving a stable amorphous layer necessitates rapid cooling rates on the order of 30 K/ps—significantly exceeding the ~2.5 K/ps predicted by conventional nTTM models. This discrepancy highlights the limitations of continuum approaches in capturing epy atomic-scale quenching dynamics and emphasizes the necessity of hybrid modeling techniques to achieve accurate predictions. Based on the results on depth of melting, the melting dynamics, and the fastest possible cooling rate that could take place due to ultrashort laser pulse energy deposition, we can draw a conclusion that minimization of the amorphous layer on the Si surface can be obtained if working at the influences just above the melting threshold. Understanding these cooling constraints enables precise control over amorphization depth, enhancing reproducibility and functionality in next-generation semiconductor devices.

Finally, the femtosecond laser parameters explored in this work not only align with experimental conditions but also provide a framework for optimizing laser-based manufacturing techniques. By ensuring minimal thermal damage to surrounding material, these findings support the advancement of high-precision applications such as phase-change memory fabrication, integrated circuit design, and optoelectronic device production. Future research directions include implementation of the modified Tersoff interatomic potential into the MD part of the combined MD-nTTM model for an accurate description of the amorphization process, extending the model to account for multi-pulse irradiation effects, investigating alternative semiconductor materials,
and integrating real-time experimental validation. The incorporation of advanced 3D computational techniques for Crank-Nicolson algorithm of the nTTM model solution, and time-resolved pump-probe measurements will further enhance the accuracy and applicability of the combined model to industrial-scale ultrafast laser processing.

In summary, this study establishes the hybrid nTTM-MD approach as a powerful tool for simulating laser-induced phase transitions in Si, providing fundamental insights into the role of electronic excitation, pulse parameters, and crystal orientation. The results not only advance our understanding of ultrafast laser-matter interactions but also offer a solid foundation for optimizing femtosecond laser technologies in microelectronics, photonics, and beyond.

## Author Contributions

DSI has developed the model, performed the calculations and analysis and has written the manuscript. TEI has conceived the ideas, helped in the modelling and in the analysis, edited the manuscript, and obtained funds for the project.

## Conflicts of interest

The author declares no conflicts of interests.

## Acknowledgements

The authors gratefully acknowledge the project ANR-23-CE08-0029 "Laser-forming of ultrathin amorphized Si layer for microelectronics applications" (LAMORSIM) for financial support. Computer support was provided by CINES of France under the project AD010814604R1.

## References

[1] Th. J.-Y. Derrien, J. Kruger, T.E. Itina, S. Hohm, A. Rosenfeld, J. Bonse, Rippled area formed by surface plasmon polaritons upon femtosecond laser double-pulse irradiation of silicon, *Optics Express* **21**, 29643-29655 (2013)

[2] E. Granados, M. Martinez-Calderon, B. Groussin, J.-Ph. Colombier and I. Santiago, "Highly uniform silicon nanopatterning with deep-ultraviolet femtosecond pulses", *Nanophotonics* **13**, 4079-4089 (2024).

[3] Dieter Bäuerle, "Chemical Processing with Lasers", Springer Berlin, Heidelberg, 2013. https://doi.org/10.1007/978-3-662-02505-5

[4] Daniel Puerto, Javier Solis, Jan Siegel, Dynamics of phase transformations in Si and Ge upon strong excitation with UV femtosecond laser pulses, *Applied Surface Science* **666**, 160372 (2024).

[5] A. Dostovalov, K. Bronnikov, V. Korolkov, S. Babin, E. Mitsai, A. Mironenko, M. Tutov, D. Zhang, K. Sugioka, J. Maksimovic, T. Katkus, S. Juodkazis, A. Zhizhchenko, A. Kuchmizhak, "Hierarchical anti-reflective laser-induced periodic surface structures (LIPSSs) on amorphous Si films for sensing applications", *Nanoscale* **12**, 13431-13441 (2020).

[6] T. Lipka, L. Moldenhauer, J. Müller, and H. K. Trieu, "Photonic integrated circuit components based on amorphous silicon-on-insulator technology", *Photonics Research* **4**, pp. 126-134 (2016)

[7] Y. Siegal, E. N. Glezer, L. Huang, and E. Mazur, "Laser-Induced Phase Transition in Semiconductors", *Annu. Rev. Mater. Sci.* **25**, 223-47(1995).
[8] P. Stampfli and K. H. Bennemann, “Dynamical theory of the laser-induced lattice instability of silicon”, *Phys Rev B* **46**, 10686 (1992).

[9] V.P. Lipp, B. Rethfeld, M.E. Garcia, D.S. Ivanov, “Atomistic-Continuum Modeling of Short Laser Pulse Melting of Si Targets”, *Phys. Rev. B* **90**, 245306 (2014).

[10] A Mittiga, L Fornarini, R Carluccio, “Numerical modeling of laser induced phase transitions in silicon”, *Applied Surface Science* **154–155**, 112-117 (2000).

[11] T. E. Itina, O. Uteza, N. Sanner, and M. Sentis, “Interaction of femtosecond laser pulses with dielectric materials: insights from numerical modelling”, *Journal of Optoelectronics and Advances Materials*, **12**, 470-473 (2010).

[12] I.B. Bogatyrev, D. Grojo, T. Itina, Ph. Delaporte, S. Leyder, M. Sentis, W. Marine, “Non-linear absorption of 1.3µm wavelength femtosecond laser pulses focused inside semiconductors: FDTD-TTM combined computational study”, *J. Appl. Phys.* **110**, 103106 (2011).

[13] G. Calogero, D. Raciti, P. Acosta-Alba, F. Cristiano, I. Deretzis, G. Fisicaro, K. Huet, S. Kerdilès, A. Sciuto and A. La Magna, “Multiscale modeling of ultrafast melting phenomena”, *npj Computational Materials* **8**, 36 (2022).

[14] G. Calogero, D. Raciti, P. Acosta-Alba, F. Cristiano, I. Deretzis, G. Fisicaro, K. Huet, S. Kerdilès, A. Sciuto and A. La Magna, “Atomistic Insights into Ultrafast SiGe Nanoprocessing”, *J. Phys. Chem. C* **127**, 19867–19877 (2023).

[15] L. Pelaz, L.A. Marques, M. Aboy, and J. Barbolla, and George H. Gilmer, “Atomistic modeling of amorphization and recrystallization in silicon”, *Appl. Phys. Lett.* **82**, 2038 (2003).

[16] H.M. van Driel, “Kinetics of high-density plasmas generated in Si by 1.06- and 0.53-µm picosecond laser pulses”, *Phys. Rev. B* **35**, 8166-8176 (1987).

[17] J.K. Chen, D.Y. Tzou, and J.E. Beraun, “Numerical Investigation of Ultrashort Laser Damage in Semiconductors”, *Int. J. Heat and Mass Transfer* **48**, 501 (2005).

[18] S.I. Anisimov, B.L. Kapeliovich, T.L. Perel'man, Electron emission from metal surfaces exposed to ultrashort laser pulses, Sov. Phys. JETP **39**, 375 (1974).

[19] A.B. Sproul and M.A. Green, Improved value for the silicon intrinsic carrier concentration from 275 to 375 K, *J. Appl. Phys.* **70**, 846 (1991).

[20] R.F. Wood and G. E. Giles, Macroscopic theory of pulsed-laser annealing. I. Thermal transport and melting, *Phys. Rev. B* **23**, 2923 (1981).

[21] R.F. Wood and G. E. Giles, Macroscopic theory of pulsed-laser annealing. I. Thermal transport and melting, *Phys. Rev. B* **23**, 2923 (1981).

[22] D. Agassi, Phenomenological model for picosecond pulse laser annealing of semiconductors, *J. Appl. Phys.* **55**, 4376 (1984).

[23] C.D. Thurmond, The Standard Thermodynamic Functions for the Formation of Electrons and Holes in Ge, Si, GaAs, and GaP, *J. Electrochem. Soc.* **122**, 1133 (1975).

[24] R. Vankemmel, W. Schoenmaker, and K. De Meyer, A unified wide temperature range model for the energy gap, the effective carrier mass, and intrinsic concentration in silicon, *Solid State Electron.* **36**, 1379 (1993).

[25] G. E. Jellison and F. A. Modine, Optical absorption of silicon between 1.6 and 4.7 eV at elevated temperatures, *Appl. Phys. Lett.* **41**, 180 (1982).
[26] G.E. Jellison and F. A. Modine, Optical functions of silicon between 1.7 and 4.7 eV at elevated temperatures, *Phys. Rev. B* **27**, 7466 (1983).

[27] J. Dwiezor and W. Schmid, Auger coefficients for highly doped and highly excited silicon, *Appl. Phys. Lett.* **31**, 346 (1977).

[28] J. Geist and W.K. Gladden, Transition rate for impact ionization in the approximation of a parabolic band structure, *Phys. Rev. B* **27**, 4833 (1983).

[29] J.R. Meyer, M.R. Kruer, and F.J. Bartoli, Optical heating in semiconductors: Laser damage in Ge, Si, InSb, and GaAs, *J. Appl. Phys.* **51**, 5513 (1980).

[30] D. Agassi, Phenomenological model for picosecond pulse laser annealing of semiconductors, *J. Appl. Phys.* **55**, 4376 (1984).

[31] Ioffe Physical Technical Institute, Electronic Archive "New Semiconductor Materials. Characteristics and Properties", http://www.ioffe.rssi.ru/SVA/NSM/Semicond/Si/index.html.

[32] J.R. Meyer, M.R. Kruer, and F.J. Bartoli, Optical heating in semiconductors: Laser damage in Ge, Si, InSb, and GaAs, *J. Appl. Phys.* **51**, 5513 (1980).

[33] A. Rämer, O. Osmani, and B. Rethfeld, "Laser damage in silicon: Energy absorption, relaxation, and transport", *Journal of Applied Physics* **116**, 053508 (2014).

[34] T.E. Itina and N.S. Shcheblanov, "Effect of Laser field on Collision Frequencies and Absorption during Ultra-Short Laser Interactions with Dielectric Materials", *Fundamentals of Laser-Assisted Micro- and Nanotechnologies* 2013, Proc. SPIE **9065**, 906506 (2013).

[35] D.S. Ivanov and L.V. Zhigilei, "Combined Atomistic-Continuum Modelling of Short-Pulse Laser Melting and Disintegration of Metal Films", *Phys. Rev. B* **68**, 064114 (2003).

[36] F. Stillinger and T.A. Weber, "Computer simulation of local order in condensed phases of silicon", *Phys. Rev. B* **31**, 5262 (1985).

[37] J. Tersoff, "Empirical interatomic potential for silicon with improved elastic properties," *Phys. Rev. B* **38**, 9902-9905 (1988).

[38] L. Shokeen and P. K. Schelling, "Thermodynamics and kinetics of silicon under conditions of strong electronic excitation", *J. Appl. Phys.* **109**, 073503 (2011).

[39] B. Bauerhenne, V.P. Lipp, T. Zier, E.S. Zijlstra, M.E. Garcia, "Self-Learning Method for Construction of Analytical Interatomic Potentials to Describe Laser-Excited Materials", *Phys. Rev. Lett.* **124**, 085501 (2020).

[40] L. V. Zhigilei and B. J. Garrison, "Pressure Waves in Microscopic Simulations of Laser Ablation", *Mater. Res. Soc. Symp. Proc.* **538**, 491 (1999).

[41] D.S. Ivanov, P.N. Terekhin, S.I. Kudryashov, S.M. Klementov, A.V. Kabashin, M.E. Garcia, B. Rethfeld, "The atomistic perspective of nanoscale laser ablation", chapter in "Ultrafast Laser Nanostructuring - The Pursuit of Extreme Scales" Springer Nature (2023). DOI:10.1007/978-3-031-14752-4

[42] Y. Gan and J. K. Chen, "A hybrid method for integrated atomistic-continuum simulation of ultrashort-pulse laser interaction with semiconductors", *Computer Physics Communications* **183**, 278 (2012).

[43] M.S. Grigoryeva, I.A. Kutlubulatova, S.Yu. Lukashenko, A.A. Fronya, D.S. Ivanov, A.P. Kanavin, V.Yu. Timoshenko and I.N. Zavestovskaya, "Modeling of Short-Pulse Laser Interactions with
Monolithic and Porous Silicon Targets with an Atomistic-Continuum Approach", *Nanomaterials* **13**, 2809 (2023).

[44] J. Crank and P. Nicolson, "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type", *Proc. Cambridge. Philos. Soc.* **43**, 50-67 (1947).

[45] V.P. Lipp, B. Rethfeld, M.E. Garcia, D.S. Ivanov, "Solving a System of Differential Equations Containing a Diffusion Equation with Nonlinear Terms on the Example of Laser Heating in Silicon", *Appl. Sci.* **10**, 1853 (2020).

[46] Wang, T.-Y.; Lee, Y.-M.; Chen, and C.C.-P., "3D Thermal-ADI: An Efficient Chip-Level Transient Thermal Simulator", In Proceedings of the 2003 International Symposium on Physical Design, Monterey, CA, USA, 6-9 April 2003; pp. 10-17.

[47] W. G. Spitzer and H. Y. Fan, "Determination of Optical Constants and Carrier Effective Mass of Semiconductors", *Phys. Rev.* **106**, 882 (1957).

[48] K. Sokolowski-Tinten and D. von der Linde, "K. Sokolowski-Tinten and D. von der Linde, *Phys. Rev. B* **61**, 2643 (2000).", *Phys. Rev. B* **61**, 2643 (2000).

[49] J. Bonse and S. Gräf, "Maxwell Meets Marangoni—A Review of Theories on Laser-Induced Periodic Surface Structures", *Laser and Photonics Review* **14**, 2000215 (2020).

[50] Cynthia L. Kelchner, S. J. Plimpton, and J. C. Hamilton, "Dislocation nucleation and defect structure during surface indentation", *Phys. Rev. B* **58**, 11085 (1998).

[51] G. Blatter, M. V. Feigel'man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, "Vortices in high-temperature superconductors", *Rev. Mod. Phys.* **66**, 1125 (1994).

[52] M. Garcia-Lechuga, N. Casquero, J. Siegel, J. Solis, R. Clady, A. Wang, O. Utéza, D. Grojo, "Amorphization and Ablation of Crystalline Silicon Using Ultrafast Lasers: Dependencies on the Pulse Duration and Irradiation Wavelength", *Laser Photonics Rev.* **18**, 2301327 (2024).

[53] J. Bonse, K. W. Brzezinka, A. J. Meixner, "All-optical characterization of single femtosecond laser-pulse-induced amorphization in silicon", *Appl. Surf. Sci.* **221**, 215 (2004).

[54] S. Besner, J. Y. Degorce, A. V. Kabashin, M. Meunier, "Influence of ambient medium on femtosecond laser processing of silicon", *Appl. Surf. Sci.* **247**, 163 (2005).

[55] J. Bonse, S. Baudach, J. Krüger, W. Kautek, M. Lenzner, Appl Phys A **74**, 19 (2002).

[56] W. Yang, Y. Li, "The Influence of Crystal Orientation on Subsurface Damage of Mono-Crystalline Silicon by Bound-Abrasive Grinding". *Micromachines* **12**, 365 (2021).

[57] C.C. Yang, J.C. Li, Q. Jiang, "Temperature-pressure phase diagram of silicon determined by Clapeyron equation", *Solid State Communications* **129**, 437-441 (2004).

[58] M.S. Grigoryeva, I.A. Kutlubulatova, S.Y. Lukashenko, A.A Fronya, D.A. Ivanov, A.P. Kanavin, V.Y. Timoshenko, I.N. Zavestovskaya, "Modeling of Short-Pulse Laser Interactions with Monolithic and Porous Silicon Targets with an Atomistic-Continuum Approach", *Nanomaterials* **13**, 2809 (2023).

[59] M.M. El-Nahass, H.A.M. Ali, Estimation of optical parameters of silicon single crystals with different orientations, Materials Science-Poland, DOI: 10.2478/msp-2019-0003

[60] D.S. Ivanov and L.V. Zhigilei, "Kinetic Limit of Heterogeneous Melting in Metals", *Phys. Rev. Lett.* **98**, 195701 (2007).
[61] M. Z. Mo, M. Z. Mo, Z. Chen, R. K. Li, M. Dunning, B. B. L. Witte, J. K. Baldwin, L. B. Fletcher, J.B. Kim, A. Ng, R. Redmer, A.H. Reid, P. Shekhar, X.Z. Shen, M. Shen, K. Sokolow ski-Tinten, Y. Y. Tsui, Y. Q. Wang, Q. Zheng, X. J. Wang, and S. H. Glenzer, "Heterogeneous to homogeneous melting transition visualized with ultrafast electron diffraction", *Science* **360**, 1451-1455 (2018).

[62] M.I. Arefev, M.V. Shugaev, and L.V. Zhigilei, "Kinetics of laser-induced melting of thin gold film: How slow can it get?", *Sci. Adv.* **8**, eab02621 (2022).

[63] V. I. Mazhukin, A. V. Shapranov, O. N. Koroleva, A. V. Mazhukin, "Modification of the Wilson-Frankel kinetic model and atomistic simulation of the rate of melting/crystallization of metals", *Mat. Model.* **35**, 103–121 (2023).

[64] R.L.C. Vink, G.T. Barkema, W.F. vaan Weg, and N. Mousseau, "Fitting the Stiiger-Weber potential to amorphous silicon", *Journal of Non-Crystalline Solids* **282**, 248-255 (2000).

[65] F. Stillinger and T.A. Weber, "Computer simulation of local order in condensed phases of silicon", *Phys. Rev. B* **31**, 5262 (1985).

[66] Bauerhenne, B.; Lipp, V.P.; Zier, T.; Zijlstra, E.S.; Garcia, M.E. Self-Learning Method for Construction of Analytical Interatomic Potentials to Describe Laser-Excited Materials. *Phys. Rev. Lett.* **124**, 085501 (2020).

[67] V.L. Deringer, N. Bernstein, A.P. Bartók, M.J. Cliffe, R.N. Kerber, L.E. Marbella, C.P. Grey, S.R. Elliott and Gábor Csányi, "Realistic Atomistic Structure of Amorphous Silicon from Machine-Learning-Driven Molecular Dynamics", *J. Phys. Chem. Lett.* **9**, 2879–2885 (2018).

[68] T. Kumagai, S. Izumi, S. Hara, and S. Sakai, "Development of bond-order potentials that can reproduce the elastic constants and melting point of silicon for classical molecular dynamics simulation", *Comput. Mater. Sci.* **39**, 457–464 (2007).

[69] J. Tersoff, "Empirical interatomic potential for silicon with improved elastic properties," *Phys. Rev. B* **38**, 9902–9905 (1988).

[70] P.K. Schelling, "Phase behavior and kinetics of a new bond-order potential for silicon", *Computational Materials Science* **44**, 274–279 (2008).

[71] M. He and L.V. Zhigilei, "Multiscale modeling of short pulse laser induced amorphization of silicon", *J. Appl. Phys.* **136**, 213103 (2024).