Short pulse laser processing of silicon has been widely adopted in the past decades due to its cost efficiency and suitability for high-throughput manufacturing. A broad range of laser applications in the design and fabrication of semiconductor devices includes surface microtexturing and micromachining [1][2][3], transfer of contacts [4,5], internal structuring for threedimensional silicon photonics [6], and generation of strong solute supersaturation (hyperdoping)
for photovoltaics and infrared photodetectors [7][8][9].
The formation of amorphous phase is commonly observed in short pulse laser processing of silicon and attributed to highly non-equilibrium conditions produced by the laser irradiation [10]. The amorphization can be undesired, e.g., in the manufacturing of photovoltaic devices, where the amorphous inclusions increase charge carrier recombination rate and decrease the solar cell efficiency [ 11]. When the amorphization is done in a well-controlled spatially-resolved manner, however, it opens up new opportunities for device manufacturing and surface functionalization. In particular, the difference in the etch rate of amorphous and crystalline silicon suggests that the direct laser writing of amorphous regions combined with chemical etching can be used as an effective maskless lithography approach [ 12 ]. A higher refractive index of amorphous silicon as compared to the crystalline phase can enable confined light propagation (waveguiding) for silicon photonic devices [13], although the direct laser writing of amorphous tracks that are sufficiently deep for effective waveguiding still presents a challenge [14]. The utility of silicon as a phase-change material for active photonics in the visible spectrum is suggested by a recent demonstration of reversible crystallization and amorphization of nanostructured silicon [15].
At the most basic level, the control over spatially resolved amorphization of silicon surface requires complete understanding of the conditions leading to the generation of the amorphous phase within an individual laser spot. The first studies of the morphology of amorphous regions generated by single picosecond pulse irradiation of silicon (001) and (111) surfaces [16,17] revealed a transition from an amorphous region produced at the center of the laser spot just above the fluence threshold for the onset of melting to an amorphous ring surrounding a crystalline central part of the spot at higher laser fluences. These initial reports were followed by systematic investigations of the dependence of the amorphous phase morphology and depth on the irradiation conditions (laser wavelength, fluence, pulse duration, number of pulses) and crystallographic surface orientation [14,18 -24], with most studies performed with femtosecond laser pulses. It has been demonstrated that the direct writing of partially overlapping amorphous rings can be used as an effective strategy for two-dimensional surface patterning [21], and that the depth of the amorphous region can be increased when the laser processing is done in the presence of an optically transparent overlayer providing an additional cooling channel [14].
The conditions for the laser-induced amorphization can be described in terms of a competition between the rapid cooling of a transiently melted region driven by a steep temperature gradient and the epitaxial regrowth of the crystalline part of the target. The outcome of this competition is largely defined by the non-monotonous temperature dependence of the velocity of crystallization front. According to the results of molecular dynamics (MD) simulations of the solidification kinetics in Si [25 -29], the velocity of the crystal-liquid interface 𝑣𝑣 𝑐𝑐𝑐𝑐 exhibits a nonlinear temperature dependence. As temperature decreases below the equilibrium melting temperature Tm, the initial close-to-linear increase of the interface velocity becomes weaker and, at about 0.85Tm, the velocity reaches its maximum value of 𝑣𝑣 𝑐𝑐𝑐𝑐 𝑚𝑚𝑚𝑚𝑚𝑚 . Further cooling leads to the decrease of the velocity down to the levels that cannot be resolved in MD simulations at T ≤ 0.6Tm.
The slowdown and eventual halt of the crystallization front propagation, combined with the insufficient time for the homogeneous nucleation of the new crystallites, allows the undercooled liquid to transform to the metastable amorphous state.
The distribution of the amorphous phase within an irradiated spot is commonly attributed to the variation of the laser intensity within the laser beam, usually assumed to have a Gaussian spatial profile [14,[17][18][19][20][21][22][23][24]30]. A sequence of local fluence thresholds for material modification have been identified based on the observation of annular structures formed within the laser spot [18,23,30]. The fluence thresholds include those for amorphization, oxidation or modification of native oxide layer, recrystallization, and ablation. The physical processes that define the range of laser fluences producing the amorphous phase (i.e., the lower and upper threshold fluences for amorphization) and control the depth of the amorphous regions are still debated in literature. The parabolic amorphous region thickness profiles evaluated from optical imaging of surface regions modified by femtosecond laser single pulse irradiation [19] suggest a linear relationship between the amorphization depth and the maximum depth of melting. The mapping of the amorphous layer thickness to the local laser fluence deposited by a laser beam with a Gaussian spatial profile reveals that the slope of the linear scaling is about twice smaller for (001) surface orientation as compared to that for (111) silicon surface [23]. This observation cautions against the direct association of the melting depth with the depth of the amorphous region [14], and suggests that the (fluence and surface orientation dependent) epitaxial regrowth of the crystalline part of the target is an important factor affecting the spatial extent of the amorphization. The recently suggested simple criterion associating the lower and upper fluence thresholds for the amorphization with the levels of energy density required for reaching the melting temperature and complete melting of silicon, respectively [24], also requires further justification.
In this paper, we investigate the irradiation conditions leading to the amorphization and defining the morphology of the amorphous regions using a hierarchical multiscale atomisticcontinuum computational approach. Although direct MD simulations of laser ablation, melting, and resolidification of silicon targets have been reported [31 -36], the severe limitations on the time-and length-scales accessible to the atomistic modeling make it difficult to establish quantitative links between the computational predictions and experimental data on the morphology of laser-generated amorphous regions. Therefore, in the present study, we only use MD simulations for obtaining the information required for the design and parametrization of a computationally efficient continuum-level model accounting for the nonequilibrium kinetics of melting and solidification processes. The continuum model is then applied for investigation of the generation of amorphous regions on a (001) silicon substrate irradiated by a circular laser beam with a Gaussian spatial profile. In order to avoid the complications related to the onset of nonthermal melting [22,37,38] and the distortion of the energy deposition profiles by multi-photon absorption [14,24] in the femtosecond laser interactions, the first application of the multiscale model is reported here for the amorphization produced by picosecond laser pulses.
The results of the MD simulations of the order-disorder phase transformations in silicon are reported next, in Section II. The formulation of the continuum model is provided in Section III.
The results of one-dimensional (1D) continuum simulations of laser-induced melting,
The atomistic MD modelling is used in this study to reveal the temperature dependence of the melting/solidification front propagation velocity and establish the conditions for the formation of a stable amorphous phase. The results of MD simulations are then used in the design and parameterization of the continuum model discussed in Section III.
The interactions between silicon atoms are described by a modified Tersoff potential [39] where the angular-dependent term of the original Tersoff potential [40] is modified to improve the description of the melting transition. The modified Tersoff potential reproduces the cohesive energy and lattice parameter of crystalline Si, provides a good description of the experimental elastic constants of Si, and predicts the melting temperature of Tm = 1681 K, which is very close to the experimental value of 1687 K [41]. The description of liquid and amorphous phases by the modified Tersoff potential has also been demonstrated to be improved [42] with respect to the commonly used Stillinger-Weber potential [ 43 ]. The modified Tersoff potential still underestimates the latent heat of crystal-liquid phase transition, predicting 33.7 kJ/mol versus an experimental value of 50.25 kJ/mol [41], although the deviation is smaller compared to the value of 30.9 kJ/mol obtained with the Stillinger-Weber potential [44].
The transition to the amorphous state upon rapid cooling of molten silicon is first investigated in a series of small-scale MD simulations performed for a cubic system of 4096 Si atoms with a size of about 4 nm and periodic boundary conditions applied in all directions. The pressure and temperature are maintained using the Berendsen barostat and thermostat algorithms [45], and the pressure is maintained at a constant zero level in all simulations. Note that pressure fluctuations in MD simulations make the results obtained at zero pressure and 1 atm indistinguishable. Starting from the room-temperature crystalline silicon, the temperature is gradually increased up to the limit of superheating, about 1.3Tm, when a spontaneous homogeneous melting takes place. The liquid system is then thoroughly equilibrated at Tm, and the temperature dependence of the thermodynamic and structural characteristics of the disordered system are obtained by heating and cooling the liquid phase with a rate of 10 12 K/s. As discussed below, this cooling rate is comparable to those realized in the picosecond laser processing of silicon.
The internal energies of the crystalline and disordered phases of silicon are plotted in Figure 1a. The curve for the disordered phase exhibits a pronounced drop of the internal energy in the range of temperatures from 0.83𝑇𝑇 𝑚𝑚 (1400 K) to 0.71𝑇𝑇 𝑚𝑚 (1200 K), marked by grey color in the plot. In the same range of temperatures, the density of undercooled liquid decreases by about 9% (Figure 1b), signifying the transition from the high-density liquid state to a low-density disordered phase. As noted in a previous study [42], this transition is reversible, albeit a small hysteresis of the density curves obtained upon cooling and heating is observed at the rates of temperature variation used in MD simulations. The Abraham factor [46], defined as a ratio of the magnitudes of the pair distribution function at its first minimum and first maximum, 𝑔𝑔 𝑚𝑚𝑚𝑚𝑚𝑚 /𝑔𝑔 𝑚𝑚𝑚𝑚𝑚𝑚 , also decreases rapidly over the same temperature range, and exhibits a discontinuous change in slope at the lower end of this temperature interval. Following the common association of the change in slope of the temperature dependence of the Abraham factor with the transition to an amorphous state [47,48], we define the glass transition temperature as 𝑇𝑇 𝑔𝑔 = 0.71𝑇𝑇 𝑚𝑚 .
We note that, due to the difference in the nature of interatomic bonding of the high-density metallic liquid phase and the low-density covalently-bonded amorphous phase, as well as the distinct thermodynamic and structural characteristics of the two phases, the liquid -amorphous phase transition is usually assumed to be of the first order [49,50]. Thus, the very definition of the glass transition temperature becomes questionable. The results of MD simulations [48,51,52],
however, suggest that the transition to the amorphous state is likely to be preceded by a first-order phase transition of a dense high-temperature liquid phase to a low-density tetrahedrallycoordinated liquid, which then undergoes a glass transition upon further cooling [ 53 ]. The discontinuous (or non-monotonous) temperature dependences of the internal energy and density expected for the first-order liquid-liquid phase transition [51,52] do not appear in the plots shown in Figure 1a,b, due to the high quench rate used in the simulations. A recent analysis of the effect of the quench rate on the evolution of density and Abraham parameter during quenching in MD simulations [48] reveals qualitatively similar trends in the range of the quench rates from 10 12 to 5×10 10 K/s. The quench rates realized in the picosecond laser processing of silicon fall within this range (see Section IV), suggesting that the phase transformations predicted in the MD simulations are representative of those occurring in the laser processing.
Based on the above discussion, we define the disordered phase below 𝑇𝑇 𝑔𝑔 as the amorphous phase. Since the simulations are performed at zero pressure, the distance between the two curves of the internal energy in Figure 1a corresponds to the enthalpy or latent heat of the order-disorder phase transformation. The latent heat drops by about a factor of two upon the amorphization, from ΔHm = Hl -Hc ranging from 0.33 to 0.36 eV/atom (32 to 35 kJ/mol) above 0.83Tm to ΔHa = Ha -Hc ranging from 0.18 to 0.19 eV/atom (17.4 to 18.3 kJ/mol) below 0.71Tm. These values are consistent those reported earlier for the same interatomic potential, ΔHm = 33.7 kJ/mol and ΔHa = 18.6 kJ/mol [42].
The transition from the undercooled liquid to the amorphous state is also signified by a sharp drop, by more than an order of magnitude, in the diffusion coefficient calculated for the disordered phase from the time dependence of the mean square displacement of atoms, Figure 1d.
For the high density liquid phase above Tm, the temperature dependence can be relatively well described by the Arrhenius equation yielding 𝐷𝐷(𝑇𝑇) = 2.03 × 10 -3 𝑒𝑒𝑒𝑒𝑒𝑒( -0.47 eV/𝑘𝑘 𝐵𝐵 𝑇𝑇)
[cm 2 /s], although the undercooled liquid exhibits a significant deviation from the Arrhenius behavior, as also observed in earlier simulations [51]. The data points around and below 𝑇𝑇 𝑔𝑔 can also be fitted the Arrhenius equation, yielding 𝐷𝐷(𝑇𝑇) = 1.45 × 10 -1 𝑒𝑒𝑒𝑒𝑒𝑒( -1.23 eV/𝑘𝑘 𝐵𝐵 𝑇𝑇)
[cm 2 /s], although the parameters of this fit are expected to be sensitive to the quench rate used in the simulations.
The variation of the thermodynamic and kinetic parameters discussed above has direct implications for the temperature dependence of the velocity of the interface between the crystalline and disordered phases. The interface kinetics has been investigated in MD simulations performed with different interatomic potentials and crystallographic orientations of the interface [25][26][27][28][29]. To provide an analytical description of the melting and solidification processes consistent with the discussion of the amorphization provided above, we perform a series of two-phase MD simulations
with the same modified Tersoff potential [39] that is used to produce data shown in Figure 1. The two-phase systems are generated by combining crystalline and disordered parts prepared and equilibrated at the same temperature in separate simulations. Both parts have the same lateral dimensions of about 22 × 22 nm 2 , chosen to match the dimensions of the crystalline part at the target temperature and zero pressure. The two parts are then joined along the axis normal to (001) face of the crystalline part of the system. Periodic boundary conditions are applied in all three directions, leading to the formation of two interfaces separating the disordered and crystalline regions. The total system is composed of 1 228 800 silicon atoms and has a length of about 52 nm along the axis normal to the interfaces. In simulations of solidification performed at T < Tm, the crystalline part in the initial system is four times smaller than the liquid part, while in systems prepared for simulations of melting at T > Tm the proportion is reversed. The simulations of melting and solidification are performed under conditions of zero pressure maintained by the Berendsen barostat algorithm [45] applied only in the direction normal to the interfaces, while the lateral dimensions are fixed to maintain the equilibrium lattice parameter in the crystalline part of the system. To avoid the local increase (decrease) of temperature at the interfaces due to the release (absorption) of the latent heat during the solidification (melting) [54][55][56][57], the system is divided into 4-nm-thick slices parallel to the interfaces, and the local temperature is maintained at a constant level by applying the Berendsen thermostat algorithm [45] to each individual slice. The interface velocity is calculated by tracking the total volume of the two-phase system [58], and the simulations are performed in the regimes of melting (at 1.07Tm, 1.13Tm and 1.19Tm) and crystallization (at 0.95Tm, 0.9Tm, 0.85Tm, 0.8Tm, and 0.75Tm). To obtain statistically reliable data, four simulations are performed at each temperature, and the standard deviations of the mean velocity values are calculated.
The results of the calculations of the melting or crystallization front velocity are shown in Figure 2. Above Tm, the melting front velocity exhibits a super-linear temperature dependence, also observed in simulations of metals [59]. The velocity of the melting front increases up to the limit of crystal superheating, where a rapid homogeneous nucleation of liquid regions takes place inside the superheated crystal [60]. In contrast, the crystallization front velocity below Tm exhibits a sub-linear temperature dependence at moderate levels of undercooling, reaches a maximum value of about 13.8 m/s at about 0.85Tm, and then decreases upon further temperature decrease. In all simulations performed below Tm, the crystallization proceeds by propagation of relatively flat interfaces, with no homogeneous nucleation of new crystallites observed inside the undercooled liquid. The simulations are continued until complete solidification of the system, and the simulation time is varied between 1 and 3 ns depending on the interface velocity. The range of temperatures where the velocity of the interface decreases with increasing level of undercooling largely coincides with that marked by the grey color in Figure 1, where the transformation of undercooled liquid to an amorphous phase is inferred to occur from the analysis of thermodynamic, structural and kinetic parameters of silicon. The lower boundary of this temperature range is identified as the glass transition temperature, Tg. Below Tg, the velocity drops down to the levels that cannot be resolved in nanosecond-scale MD simulations. As shown in continuum-level simulations of picosecond pulse laser processing of Si discussed below, in Sections IV and V, the cooling of the surface region of the irradiated targets also occurs on the timescale of nanoseconds, suggesting that the cooling process can overrun the epitaxial regrowth of the crystalline part of the target, freezing the undercooled liquid into a kinetically stable amorphous phase. Indeed, the drop of the crystallization front velocity down to zero under conditions of deep undercooling appears to be a prerequisite for generation of a stable amorphous phase in rapid quenching experiments. The inability to produce a stable amorphous state for onecomponent face centered cubic (fcc) metals even in experiments where extreme quench rates of about 10 13 K/s are realized [61] can be explained by a weak decay of the crystallization front velocity with increasing undercooling predicted in MD simulations [57,[62][63][64]. On the other hand, the velocity of the crystallization front is predicted to vanish under conditions of deep undercooling for some of the body centered cubic (bcc) metals [29,57,62,63], which explains the formation of stabile amorphous vanadium and tantalum in the same extreme quenching experiments [61].
In order to formulate the predictions of MD simulations in a form suitable for incorporation into a continuum-level model, we describe the MD data points shown in Figure 2 by the Wilson -Frenkel equation. This equation describes the temperature dependence of the velocity of the crystal-liquid interface 𝑣𝑣(𝑇𝑇) based on the consideration of the exchange of atoms between the two phases within the phenomenological transition state theory [65][66][67]:
where 𝐷𝐷(𝑇𝑇) is the atomic diffusivity in the disordered phase, Δ𝐻𝐻(𝑇𝑇) is the latent heat of the order -disorder phase transformation, 𝐶𝐶 0 is a constant pre-factor that depends on the details of atomic rearrangements responsible for the movement of the interface, and 𝑘𝑘 𝐵𝐵 is the Boltzmann constant.
This equation accounts for the temperature dependence of the thermodynamic driving force, i.e., the free energy difference between the two phases approximated here as Δ𝐺𝐺(𝑇𝑇) ≈ Δ𝐻𝐻(𝑇𝑇)(𝑇𝑇 𝑚𝑚 -𝑇𝑇)/𝑇𝑇 𝑚𝑚 , as well as the atomic mobility at the interface described by 𝐷𝐷(𝑇𝑇). As a result, it is capable of providing a unified description of the non-monotonous temperature dependence of the interface velocity in both melting and crystallization regimes. The super-linear increase of the melting front velocity with increasing superheating above Tm is defined by the simultaneous increase of the driving force and atomic mobility. The sub-linear increase turning into a decrease of the crystallization front velocity with increasing undercooling below Tm is defined by the competition of the increasing driving force and the decreasing atomic mobility, with the dominance of the latter at strong undercoolings.
To apply Eq. ( 1) to the description of the MD results, the minimal fitting approach is applied. Namely, the fitting is limited to choosing the value of the constant pre-factor, 𝐶𝐶 0 = 1.12×10 11 m -1 , based on the three data points obtained in MD simulations of melting. The temperature dependences of the latent heat of the phase transformation Δ𝐻𝐻(𝑇𝑇) and the diffusivity in the disordered phase 𝐷𝐷(𝑇𝑇) are taken from the results of independent MD simulations shown in Figures 1a and 1d, respectively. The solidification velocities obtained in MD simulations are not used in fitting of Eq. ( 1), but are still described reasonably well by the red line that shows the prediction of the kinetic equation in Figure 2.
The focus of the present study on the laser-induced amorphization makes the continuumlevel models that assume heat-flow limited interface kinetics formulated within the framework of the Stefan problem [ 68 -70] or describe the deviations from the local equilibrium through a linearized interface response function [71,72] unsuitable. As discussed in Section II, a complete non-equilibrium kinetic description, fully accounting for the non-monotonous temperature dependence of the interface velocity [67], is required for a realistic representation of the amorphization process. In the present study, the temperature dependence of the velocity of the crystal-liquid interface is described by Eq. ( 1
where 𝑘𝑘 𝜕𝜕ℎ , 𝑐𝑐 𝑃𝑃 , and 𝜌𝜌 are the thermal conductivity, heat capacity, and density of silicon, respectively. The source terms correspond to the laser energy deposition, 𝑆𝑆 𝐿𝐿 , the latent heat of melting/solidification, 𝑆𝑆 𝑚𝑚 , and the energy loss due to vaporization from the target surface, 𝑆𝑆 𝑒𝑒 . In the simulations reported in this paper, the laser pulse duration ranges from 30 ps to 1.5 ns, which is sufficient for the thermal equilibration between the excited electrons and phonons [17,73].
Therefore, the electron-phonon energy exchange and two-temperature heat diffusion for electron and lattice subsystems are not considered.
The source term describing the laser energy deposition represents the irradiation of the target surface by a circular laser beam with Gaussian temporal and spatial profiles,
where 𝐼𝐼 0 is the laser peak intensity, 𝑅𝑅 is the reflectivity defined by the temperature and phase state of the surface, 𝐿𝐿 𝑝𝑝 is the optical penetration depth, σ 𝜕𝜕 and σ 𝑟𝑟 are the standard deviations of the Gaussian function describing the temporal and spatial profiles, respectively, and 𝑡𝑡 0 is the time when the laser intensity reaches its maximum. The laser pulse duration 𝜏𝜏 𝐿𝐿 , defined as full width at half maximum of the Gaussian profile, is related to σ 𝜕𝜕 as 𝜏𝜏 𝐿𝐿 = σ 𝜕𝜕 �8ln (2). The laser spot diameter 𝐷𝐷 𝐿𝐿 , defined as full width at 1/e 2 height of the Gaussian spatial profile, is related to σ 𝑟𝑟 as 𝐷𝐷 𝐿𝐿 = 4σ 𝑟𝑟 . The laser fluence 𝐹𝐹 0 is related to the peak intensity 𝐼𝐼 0 as 𝐹𝐹 0 = √2𝜋𝜋σ 𝜕𝜕 𝐼𝐼 0 . To ensure the absorption of the whole laser pulse energy, the center of the Gaussian temporal profile 𝑡𝑡0 is shifted with respect to the start of the simulations by 2.5𝜏𝜏𝐿𝐿, so that the laser energy deposition before the start of the simulations at 𝑡𝑡 = 0 is negligible.
Most simulations reported in this paper are performed at a laser wavelength λ of 532 nm, although several additional 1D simulations are also done at λ = 355 and λ = 1064 nm to evaluate the effect of the laser penetration depth on the amorphization behavior. The laser source term given by Eq. ( 3) assumes the exponential attenuation of the laser intensity with depth under the surface following the Beer-Lambert law, with the temperature dependences of the optical penetration depth Lp taken from experimental studies [74][75][76]. The constant values of reflectivity of liquid silicon and the temperature dependences of the reflectivity of crystalline silicon are also adapted from experimental measurements [77][78][79][80].
The heat diffusion equation, Eq. ( 2) is solved with explicit finite difference method for a 10 μm-thick cylindrical system with a radius Rm chosen depending on the laser spot diameter 𝐷𝐷 𝐿𝐿 .
In particular, 𝑅𝑅 𝑚𝑚 = 30 μm for 𝐷𝐷 𝐿𝐿 = 5 μm, 𝑅𝑅 𝑚𝑚 = 300 μm for 𝐷𝐷 𝐿𝐿 = 10 μm, 30 μm, and 50 μm, and 𝑅𝑅 𝑚𝑚 = 500 μm for 𝐷𝐷 𝐿𝐿 = 70 μm. The thickness and radius of the computational system are chosen to ensure that no significant changes of temperature are observed at the bottom or the sides of the system during the simulations. The temperature of the initial system is set to 300 K, and adiabatic boundary conditions are applied at the bottom and at the sides of the system. The resolutions of the depth and radial spatial discretization are ∆𝑧𝑧 = 2 nm and ∆𝑟𝑟 = 100 nm, respectively. The time step in the finite difference integration is chosen to satisfy the von Neumann stability criterion [81],
The 1D modeling reported in Section IV only accounts for the thermal diffusion along the z-axis, and the contribution of the lateral heat transfer is neglected, i.e., the first term in the righthand part of Eq. ( 2) is omitted. This approximation is suitable for irradiation conditions where the laser spot diameter 𝐷𝐷 𝐿𝐿 is much larger than the characteristic depth of the irradiated target affected by the heat transfer and laser-induced phase transformations.
To provide a realistic description of the spatial propagation of the melting/solidification front in 2D simulations, the crystal-liquid interface is represented by a series of discrete points distributed along the interface. The points are initially placed at the free surface of the target at 𝑧𝑧 = 0, with the distance between points set to around 1 nm. When the laser heating brings the temperature of the target surface above Tm, the points start to move into the bulk of the target, representing the onset of surface melting. The motion of each point proceeds with a velocity of magnitude defined by the local temperature at the location of the point, as specified by Eq. ( 1), and in the direction normal to the liquid-crystal interface. The local temperature at the location of a point is defined by the linear extrapolation between the temperatures of the two closest neighboring cells of the finite-difference discretization along the 𝑧𝑧-axis, i.e., along the direction of the high temperature gradient generated by the laser energy deposition. The direction of the velocity vector is defined by the average angle 𝛼𝛼 between the r-axis and the two segments connecting the point of interest to the two neighboring points along the interface. The origin of the coordinate system is placed at the initial surface of the target in the center of the laser spot, and the 𝑧𝑧and 𝑟𝑟-axes are directed into the bulk of the target and from the center of laser spot to its periphery, respectively. The displacement of an interfacial point in 𝑧𝑧 and 𝑟𝑟 directions during a time step Δ𝑡𝑡 is then expressed as ∆𝑧𝑧 = -𝑣𝑣(𝑇𝑇)∆𝑡𝑡 cos 𝛼𝛼 and ∆𝑟𝑟 = -𝑣𝑣(𝑇𝑇)∆𝑡𝑡 sin 𝛼𝛼, respectively. Within this approach, the negative values of 𝑣𝑣(𝑇𝑇) predicted by Eq. ( 1) correspond to the advancement of the crystal-liquid interface towards the solid bulk of the target, i.e., melting at T > Tm, while the positive values of 𝑣𝑣(𝑇𝑇) correspond to the propagation of the solidification front towards the molten surface, i.e., solidification at T < Tm.
Note that the decrease of the solidification velocity with increasing undercooling, predicted by Eq. ( 1) and illustrated by Figure 2, may prevent the propagation of the solidification front all the way to the surface of the target. When the rate of cooling exceeds the rate of the epitaxial recrystallization, the advancement of the crystallization front can be ceased, and the supercooled liquid can transform into the amorphous phase. Based on the analysis of the results of MD simulations discussed in Section II, the glass transition temperature is defined as Tg = 0.71Tm, and the motion of the points defining the interface is stopped once the condition T < Tg is realized.
The possibility of phase transformations occurring through the homogeneous nucleation of a new phase is also included to account for the formation of liquid regions inside a crystal superheated up to the limit of its stability against melting [82,83] or nucleation and growth of new crystallites in a strongly undercooled liquid region [84,85]. The temperature dependence of the nucleation rate and critical nuclei radius are related to the free energy difference between the solid and liquid phases using the classic nucleation theory [82,86]. The changes of the phase fractions and temperature in each cell of the finite difference discretization due to the nucleation and growth of liquid or crystalline regions are computed based on the Kolmogorov-Johnson-Mehl-Avrami (KJMA) equation [ 87 , 88 ]. For the irradiation conditions considered in this paper, the homogeneous melting is found to play a substantial role at the initial stage of the melting process, while the homogeneous crystallization is not activated. The latter observation can be related to much higher cooling rates produced by picosecond pulse laser irradiation as compared to those where the homogeneous nucleation is predicted to occur [84].
The simulations are performed at laser fluences where the laser-induced phase transformations are limited to melting, crystallization, amorphization, and mild evaporation from the surface, while the ablation proceeding through the collective material ejection through spallation or phase explosion [82] does not occur. The material removal from the target, therefore, is limited to the evaporation described by the Hertz-Knudsen equation [89], with the corresponding energy loss included in Eq. ( 2) through the source term 𝑆𝑆 𝑒𝑒 .
The material parameters in Eq. ( 2) are taken as the experimental values [41,90,91],
including the melting temperature 𝑇𝑇 𝑚𝑚 = 1687 K, enthalpy of melting ∆𝐻𝐻(𝑇𝑇 𝑚𝑚 ) = 50.2 kJ/mol, boiling temperature 𝑇𝑇 𝑏𝑏 = 3505 K, and density 𝜌𝜌 = 2320 kg/m 3 . The thermal conductivity 𝑘𝑘 𝜕𝜕ℎ is implemented as a function of temperature for both crystalline [92] and liquid [93] silicon. The temperature dependence of the heat capacity 𝑐𝑐 𝑃𝑃 is fitted based on experimental measurements for crystal [41], liquid [94,95], and amorphous [96] phases.
Using the 1D version of the continuum-level model described above, the irradiation conditions for the generation of an amorphous surface region are analysed in this section based on the results of simulations performed for different laser fluences, pulse durations, and wavelengths.
The fluence dependence is illustrated in Figure 3, where the time evolution of the temperature and position of the crystal-liquid/amorphous interface are plotted for a pulse duration of 30 ps and a laser wavelength of 532 nm. In the simulation performed at 0.18 J/cm 2 , slightly above the melting threshold of 0.15 J/cm 2 , the temperature increases to about 1.3Tm during the laser pulse, causing the propagation of the melting front down to about 9 nm under the surface.
The steep temperature gradient results in a rapid cooling of the crystal-liquid interface. The velocity of the interface changes its sign when the temperature drops below Tm at about 120 ps, and the melting turns into resolidification. The epitaxial crystallization front reaches the surface by about 0.8 ns (marked by the blue star in Figure 3), thus leading to the recovery of the original crystalline structure of the surface region.
The increase in the laser fluence leads to strong transient superheating of the surface region above 1.3Tm, triggering the rapid homogeneous melting of larger surface regions of the targets.
The larger melt depth allows for more time for the crystal-liquid interface to cool down as it moves towards the surface. As the interface cools down, the velocity reaches its maximum value at about 0.85Tm, and then decreases upon further cooling, as described by Eq. ( 1) and illustrated by Figure 2. The decrease in the velocity of the crystal-liquid interface propagation reduces the rate of the latent heat release at the solidification front and further accelerates the cooling of the interface, as illustrated by the temperature profiles shown for laser fluences of 0.2 J/cm 2 and 0.25 J/cm 2 in Figure 3a. The solidification fronts in these simulations stop when the temperature of the fronts drop down to Tg (marked by hexagons in Figure 3), leading to the formation of kinetically stable amorphous regions with thicknesses of 4.6 nm and 16 nm at 0.2 J/cm 2 and 0.25 J/cm 2 , respectively. Further increase of the laser fluence to 0.3 J/cm 2 , however, does not produce a larger amorphous region but, on the contrary, results in the complete disappearance of the amorphous region. As can be seen from the temperature profile for this fluence in Figure 3a, further decrease in the cooling rate during the resolidification creates the conditions when the release of the latent heat at the solidification front is barely compensated by the heat transfer to the bulk of the target.
As a result, the temperature of the interface stays at an almost constant level, between 0.82Tm and 0.85Tm for about 4 ns, from 1.3 to 5.3 ns, while the corresponding velocity of the interface remains close to its maximum level. The interface temperature and velocity are even increasing from 0.82Tm and 11 m/s at 3.5 ns up to 0.85Tm and 13 m/s at 5.3 ns due to the rapid release of the latent heat of the phase transformation at the rapidly advancing solidification front. At 5.3 ns, the solidification front reaches the surface, leading to the recovery of the original crystalline structure of the target.
Thus, based on the results discussed above, we can define two threshold fluences, 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 and 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ , for the generation and disappearance of an amorphous surface region as the laser fluence increases. Below 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 , the melting depth is too shallow to produce a sufficiently strong undercooling at the solid-liquid interface before the completion of the epitaxial crystal growth.
Above 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ , the temperature gradient established by the time the interface cools down to Tm is smaller, and the corresponding heat flux from the interface is too low for bringing the interface temperature down to Tg before the complete regrowth of the crystal structure. The kinetically stable amorphous phase can, therefore, only be formed in a range of fluence between 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 and 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ .
In addition to the laser fluence, the formation of an amorphous region is also affected by the laser pulse duration. The effect of the laser pulse duration is illustrated in Figure 4 To establish the dependence of the threshold fluences 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 and 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ on the laser pulse duration, the fluence was systematically varied with a step of 0.005 J/cm 2 for each of the four pulse durations, and the conditions for surface amorphization are established in the parameter space of laser fluence and pulse duration, Figure 5a. The fluence range for the amorphization is shrinking with increasing pulse duration, as the lower fluence threshold 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 increases and the upper fluence threshold 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ decreases up to the pulse duration of 92 ps, above which no formation of amorphous structure is observed. Note that the material removal due to the evaporation from the surface is still negligible in the irradiation regime where the surface amorphization is observed, with the evaporation depth of 0.005 nm predicted for F = 0.29 J/cm 2 and 𝜏𝜏 𝐿𝐿 = 30 ps. The threshold for the onset of the phase explosion, defined as a fluence Fpe that brings the surface to Tpe = 0.9Tc [82], where Tc = 6666 K is the thermodynamic critical temperature of Si [97], is substantially higher that the upper fluence threshold 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ for surface amorphization. Similar to the threshold for the surface melting (dashed line in Figure 5), the threshold for the onset of the phase explosion exhibits the expected 𝐹𝐹 𝑝𝑝𝑒𝑒 ~√𝜏𝜏 𝐿𝐿 dependence and increases from 0.37 J/cm 2 at 𝜏𝜏 𝐿𝐿 = 30 ps up to 0.51 J/cm 2 at 𝜏𝜏 𝐿𝐿 = 90 ps. The combined effect of the laser fluence and pulse duration on the thickness of the amorphous region is summarized in Figure 5b. For a fixed pulse duration, the thickness of the amorphous region first increases gradually with increasing fluence and then drops sharply as the fluence approaches 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ , producing the asymmetric bell-shaped curves plotted in Figure 5b. The asymmetric shapes of the curves in Figure 5b can be explained by the competing effects of the increasing maximum melting depth and the decreasing temperature gradient at a time when the crystal-liquid interface cools down to Tm with increasing laser fluence. As discussed above, the increase in the maximum melting depth prolongs the time for cooling of the interface, while the decrease in the temperature gradient reduces the cooling rate and prevents the interface from reaching Tg before the complete crystallization.
An additional parameter that affects the formation of an amorphous region is the laser wavelength λ, which defines the energy deposition to the irradiated material. To investigate the effect of the wavelength, two additional series of simulations are performed for λ = 355 nm and λ = 1064 nm. The temperature dependences of the optical absorption depth for the three values of laser wavelength considered in the simulations are plotted in Figure 6a. While the optical absorption depth of the (metallic) liquid silicon is very low, has weak sensitive to laser wavelength, and ranges from 8 to 12 nm, the absorption depth of crystalline Si varies by orders of magnitudes with wavelength, resulting in the correspondingly large variation in the length-scales of the laser energy deposition. The decrease of λ from 532 nm to 355 nm causes a substantial drop in the laser penetration depth and results in the generation of a of stronger temperature gradient in the surface region of the Si target. As discussed above, the stronger temperature gradient and the correspondingly higher cooling rate at the solidification front are the factors that facilitate the formation of an amorphous region. Indeed, as illustrated in Figure 6b, the range of irradiation parameters leading to the formation of an amorphous surface region is much broader at λ = 355 nm as compared to that at λ = 532 nm, Figure 5a. In particular, the pulse durations that can lead to the amorphization are much longer and extend up to 1.44 ns. The range of fluences producing the amorphization is also wider and corresponds to five-fold variation at 𝜏𝜏 𝐿𝐿 = 250 ps (Figure 6b), as compared to the fluence variation by a factor of 1.5 at 𝜏𝜏 𝐿𝐿 = 30 ps and λ = 532 nm (Figure 5a). The depth of the amorphous region can also be substantially higher at λ = 355 nm, e.g., it is 69.7 nm at 𝜏𝜏 𝐿𝐿 = 250 ps and F = 0.5 J/cm 2 . Similar to the case of λ = 532 nm, the evaporation rate remains low in the regime of surface amorphization at λ = 532 nm, e.g., the maximum evaporation depth is predicted to be about 0.1 nm at F = 0.5 J/cm 2 and 𝜏𝜏 𝐿𝐿 = 250 ps.
In contrast to the results discussed above for wavelengths of 355 and 532 nm, an increase of the wavelength to 1064 nm eliminates the possibility to produce surface amorphization by melting and resolidification with any choice of irradiation parameters. As can be seen from Figure 6a, the room temperature penetration depth at λ = 1064 nm increases up to hundreds of micrometers, which leads to an almost uniform heating of a large part of the irradiated Si target up to a temperature close to Tm. The onset of surface melting, however, leads to a drastic change of the optical properties of the target surface. While the reflectivity more than doubles, from about 0.3 to 0.76 [77,79], the sharp drop, by almost two orders of magnitude, of the optical absorption depth (Figure 6a) leads to a strong localization of the laser energy deposition in the molten surface layer of the target.
The dramatic change of the optical properties upon melting results in a sharp fluence dependence of the maximum surface temperature and the depth of melting just above the threshold for surface melting. This observation is illustrated in Figure 7 by the results of simulations performed at two pulse durations, 30 ps and 1 ns, and two laser fluences, 26 J/cm 2 and 26.2 J/cm 2 .
At both pulse durations, the melting temperature is barely reached by the end of the laser pulse at 26 J/cm 2 , leading to the transient melting of 1.8 nm and 3.8 nm surface layers of the targets for 30 ps and 1 ns pulse irradiation. The resolidification of the molten layers proceeds at very moderate levels of undercooling, with temperature of the crystal-liquid interface decreasing down to 0.988Tm
for 𝜏𝜏 𝐿𝐿 = 30 ps and to 0.986Tm for 𝜏𝜏 𝐿𝐿 = 1 ns by the end of the resolidification. The fluence increase by only 0.8%, from 26 J/cm 2 to 26.2 J/cm 2 , has a large effect on the maximum temperature of the target surface and the depth of the target affected by melting. In particular, the melt depth increases by more than an order of magnitude, from 1.8 nm to 19.6 nm for 𝜏𝜏 𝐿𝐿 = 30 ps and from 3.8 nm to 59.1 nm for 𝜏𝜏 𝐿𝐿 = 1 ns. The following resolidification still proceeds at low levels of undercoolings that do not exceed 0.985Tm in either of the simulations. These low levels of undercooling at λ = 1064 nm are related to the preheating of the deep surface regions of the irradiated targets prior to the onset of the surface melting, which reduces the temperature gradients at the solidification fronts and eliminates the possibility for the generation of an amorphous layer.
The mapping of the predictions of 1D simulations on the morphology of amorphous regions generated by flat-top (top-hat) laser beams is straightforward for beam diameters that are sufficiently large for making the radial (parallel to the surface) heat transfer during the time of the laser-induced melting and resolidification negligible. In these cases, the conclusions on the dependence of the amorphization conditions on laser wavelength, fluence, and pulse duration discussed in the previous section can be directly applied to the whole region of the target irradiated by the laser beam.
In a more common case of the Gaussian spatial distribution of the laser beam intensity, the existence of the upper and lower fluence thresholds for surface amorphization, 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 and 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ , leads to a more complex fluence dependence of the morphology of an amorphous region generated within the laser spot. This is illustrated in Figure 8, where the distribution of the amorphous phase is shown for targets irradiated at a wavelength of 532 nm by 30 ps laser pulses with Gaussian beams characterized by peak fluences of 0.25 J/cm 2 and 0.3 J/cm 2 .
The peak fluence of 0.25 J/cm 2 is above 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 but below 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ . As a result, the amorphous phase is generated within a circular central region of the laser spot, Figure 8a, with the radius of the amorphous region, about 12 μm, defined by the location where the local fluence drops below 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 in the Gaussian laser beam profile, as shown in Figure 8b. Note that the radius of a region where the melting is observed in the simulation is larger, about 14.5 μm. The cooling of the crystalliquid interface during the resolidification, however, is insufficient for the amorphization of a part of the spot between 12 and 14.5 μm, leading to the recovery of the crystal structure of the surface.
In the simulation performed with a peak laser fluence of 0. 5b, where the maximum depth of amorphization is closer to 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ than to 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 . As a result, the increase in the depth of the amorphous ring is steeper on the inner side of the ring in Figure 8c.
The formation of ring-shaped amorphous regions predicted in the simulations is consistent with experimental observations of amorphous rings generated on Si substrates by picosecond [16,17] and femtosecond [18,[20][21][22][23][24] laser irradiation. Moreover, the transition from the ring pattern to the central amorphous spot close to the threshold for the laser induced amorphization has also been reported [16,17]. and 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ , predicted in 1D simulations. The dark gray bands mark the regions where the amorphous structure is formed. The dependences of the upper and lower fluence thresholds for amorphization, 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ and 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 , on laser fluence, pulse duration, and wavelength, established in 1D simulations, can be directly mapped to the morphology of amorphous regions produced by the Gaussian laser beam irradiation. To exemplify such mapping, a series of simulations is performed at a laser wavelength of 532 nm and different values of laser pulse durations and peak fluences. In each simulation, the presence and the shape of the amorphous region is determined, and the data points reflecting this information are used for establishing the boundaries separating the regions where distinct surface structures are produced. The resulting map, plotted in Figure 9a, shows the regions in the parameter space where the formation of a central amorphous spot, an amorphous ring pattern, and the recovery of the crystalline structure are observed. The boundaries between the regions coincide with the corresponding dependences of 𝐹𝐹 𝑚𝑚 𝑐𝑐𝑙𝑙𝑙𝑙 and 𝐹𝐹 𝑚𝑚 ℎ𝑚𝑚𝑔𝑔ℎ on the pulse duration shown in Figure 5a, as expected from the discussion provided above.
The computational prediction on the existence of the maximum value of the pulse duration for surface amorphization (Figure 9a) can be related to the experimental study [20] where the range of local fluences leading to the amorphization is observed to decrease with the increase in the pulse duration from 100 fs to 4.3 ps, and the maximum pulse duration for the amorphization is estimated to be 8 ps. These experiments, however, were performed at a wavelength of 800 nm, where the nonlinear absorption of the laser pulse is essential for creating conditions for surface amorphization (e.g., see discussion of the absence of amorphization for the case of linear absorption at λ = 1064 nm in Section IV). In particular, a substantial increase in both the lower and upper thresholds for amorphization with the pulse duration increasing from 100 fs, to 400 fs, and to 850 fs can only be explained by the sensitivity of the nonlinear absorption to the pulse duration. The nonlinear absorption can produce the concentration of the deposited energy in a nanoscale surface layer of the target necessary for the amorphization [24]. Indeed, time-resolved probing of melting and resolidification induced by a 120 fs laser pulse with a wavelength of 800 nm [21] suggests the formation of an amorphous ring at a timescale of 2-3 ns, followed by a longer solidification recovering the crystal structure in the center of the laser spot. These time scales are in a good agreement with the computational predictions obtained for the picosecond pulses and shorter wavelength, suggesting that the nonlinear absorption of femtosecond pulses may create the energy deposition profiles comparable to those analyzed in the present study. The direct mapping between the predictions of 1D simulations to the morphology of the amorphous regions produced by the Gaussian laser beam irradiation relies on the assumption that the heat transfer in the radial directions can be neglected on the timescale of the laser-induced