Thermodynamical Analysis of the Formation of α-Si Ring Structures on Silicon Surface

Superficial modifications on silicon wafers produced by single-shot focused femtosecond laser irradiation having a 1030 nm wavelength and 300 fs pulse duration were experimentally and theoretically analyzed. The laser fluence window when the amorphous silicon phase develops, resulting in a ring-like modification shape, was experimentally estimated to be between 0.26 J/cm2 and 0.40 J/cm2 and was independent of the silicon dopant type and laser focusing conditions; however, the window was narrower when compared to results reported for shorter pulse durations. In addition, we present a simplified numerical model that can explain and predict the formation of these patterns based on the caloric coefficients of silicon and the energy distribution of the deposited material.


Introduction
Direct laser writing (DLW) and surface patterning with ultrafast laser pulses remain one of the most versatile methods for tailoring material properties [1,2]. A wide range of light-matter interaction phenomena could lead to different technological applications such as precise laser ablation, enabling laser micromachining with minimal thermal impact [3], and the self-organized generation of nanostructures on the surface [4] helps to achieve unique surface properties while a selective modification of materials' optical properties can open new ways for integrated photonics and similar applications [5]. The universality of light absorption mechanisms for ultrafast laser pulses makes the DLW method applicable for virtually all materials: dielectrics, metals, and semiconductors; however, the reaction of each material to the absorbed radiation may vary.
For example, in semiconductors, such as crystalline silicon (c-Si), it is very common to observe amorphized regions (α-Si) at the boundaries of laser-irradiated areas. It is generally accepted that amorphization in silicon happens when the molten silicon layer cools down at a rate such that the resolidification front moves faster than the critical speed, which is in the range of 1-25 m/s [6]. At this cooling speed, the nucleation process cannot develop, resulting in amorphous silicon. The high thermal conductivity of silicon can efficiently cool down the solid-liquid interface and reach ultracritical speed in thin liquid layers. As the refractive index of (α-Si) is slightly higher than that of c-Si [7], resulting in a higher reflectivity under visible light, this region can be clearly distinguished using an optical microscope. The typical ultrafast laser fluence window for amorphization can be found in the literature and is in the range of 0.26 J/cm 2 -0.52 J/cm 2 [8]. At higher fluences, the amorphized regions abruptly disappear, indicating silicon recrystallization. Therefore, a bright ring-shaped pattern can be observed if the laser's fluence is higher than the indicated window. Such rings were first observed by Bloembergen and Liu in the 1980s using the first picosecond lasers, while the α-Si phase was confirmed using electron diffraction patterns [9]. Utilizing the sharp visibility of the α-Si phase, Liu suggested a simple method for measuring the focused laser waist size, assuming that threshold values were constant and depended only on the fluence [10]. In later studies, it was shown that more complex structures could form on silicon after local melting induced by a single ultrafast laser pulse, including amorphization, recrystallization, native oxide removal or melting, oxidation, and ablation [8,11], each of which has its own threshold fluence.
The selective superficial amorphization of silicon in tandem with ultrafast DLW has several emerging applications. α-Si has a lower selectivity for wet chemical etching that can be utilized for maskless lithography [12], and DLW flexibility and selective surface patterning [13] can be used for custom MEMS devices. In addition, laser patterning can generate features with subwavelength resolution, as demonstrated by the production of highly periodic LIPPS with alternating α-Si and c-Si regions [14][15][16][17]. α-Si surrounded by c-Si can support waveguiding [18] and the research on the DLW integration of such devices was carried out in several studies [19]; however, the successful integration has yet to be demonstrated-the laser-produced α-Si layer is too thin and prevents waveguiding for practical wavelengths. Using fluences well over the ablation threshold [20] or selecting the laser wavelength inside the silicon transparency window [19] still produces an amorphous layer of tens of nanometers in thickness.
The extensive real-time dynamic studies of semiconductor behavior under a femtosecond laser pulse were carried out by the Siegel and Bonse groups, which helped to establish time scales as well as experimentally characterize the damage/modification morphology [17,21,22]. However, there is still a lack of theoretical analysis that could explain the amorphization effects from light-matter interaction effects and the thermodynamic properties of silicon. In this study, we investigated the formation of amorphous ring structures on the silicon surface methodically and examined the ring diameter dependence on irradiation conditions and the silicon material itself (different doping concentrations and crystal orientations). The analysis of the experimental results leads to a theoretical investigation. The derived simple numerical simulations give a very good match to the experimental data when laser-deposited energy is connected to the caloric coefficients.

Experiment
An experiment was conducted to investigate the formation of ring-shaped structures on the silicon surface after single-shot laser irradiation. An Yb:KGW laser Pharos (Light Conversion Ltd., Vilnius, Lithuania) with a beam quality of M 2 equal to 1.1 that operated at a 1030 nm wavelength and produced 300 fs duration pulses was used as a light source. A silicon sample was placed on the XY translation stage ANT 180 (Aerotech, Pittsburgh, PA, USA), while a laser beam was focused with a 100 mm lens mounted on the Z axis translation stage ANT 130 (Aerotech, Pittsburgh, PA, USA). The beam waist at the lens focus was 15.4 µm (measured by Liu's method [10]). Several different beam diameters were also tested in the experiment, which were produced by shifting the sample slightly above the geometrical focus as shown in Figure 1. In this way, we were able to analyze the ring formation dependence on the beam diameter. The silicon surface was patterned in an X-Y grid with different pulse energies, which were attenuated from 0.2 µJ (no visible damage was observed on a sample at such energy) up to 8 µJ. All experimental procedures were automated using the laser microfabrication software application SCA (Workshop of Photonics, Vilnius, Lithuania)). Samples were investigated using an optical microscope profilometer, PLµ2300 (Sensofar Ltd., Barcelona, Spain) with an EPI 100X, NA = 0.9 objective, and no morphological surface change was observed till the damage of the material. Therefore, an optical microscope was used at the standard regime, and images were processed to extract the inner and outer ring diameters of amorphized regions. Scanning electron imaging was also performed; however, due to the poor contrast of the images of the ring structures, it was not investigated. Selected samples were imaged with an Olympus BX51 microscope using an oil immersion objective with NA 1.25. These are depicted in Figure 2.  In total, six polished silicon samples with different dopant concentrations and crystal orientations were used in the experiment (see Table 1). A typical result of a single laser pulse effect on the silicon surface is shown in Figure 2. When laser fluence is relatively modest, the irradiated spot appears brighter in microscope images because the amorphized silicon region, which has a higher refractive index, reflects more light. At a certain fluence, however, the reflectance reduces sharply, resulting in a well-defined ring-shaped pattern that is clearly seen in the pictures. With increasing fluence, the ring diameter expands while the ring width decreases, indicating the existence of a specific fluence window for silicon amorphization. The rings are only slightly irregular due to the quality of the laser beam, which was measured to be M 2 = 1.1. The ring's inner and outer diameters were measured by fitting two coaxial circles that best fit the data (see Figure 2 bottom right image). We did not identify any variations in ring formation thresholds between different dopant levels in silicon samples. A slight reduction of damage threshold for doped silicon might be visible in Figure 2 when the highest fluence was used. There was a striking difference in crystal orientation: amorphized rings were only distinguishable in the 111 crystal orientation cut wafers, whereas the 100 samples were scarcely discernible-only with digitally enhanced images was the modified region evident, as seen in the bottom section of Figure 2. Despite the difference in visibility, the modified region size for 111 samples and 100 samples under the same focusing condition was nearly identical, indicating the same modification threshold. In this paper's Supplementary Material, we provide raw optical microscope images comparing the modifications on all samples under identical focusing conditions (see Animation S1), as well as confocal microscope images made with a monochromatic (460 nm) LED illumination with enhanced contrast (see Animation S2). The later regime helped to identify and measure the boundaries of the inner and outer modification regions with greater clarity; therefore, these data were used in the analysis.
It was reported that high-contrast amorphized regions in 100 silicon developed only after multipulse irradiation [13]; nevertheless, Florian et al. showed that clear rings could also be produced with single-pulse irradiation [11]; however, in that study, the laser source had a pulse duration of up to an order of magnitude shorter than used in our study, demonstrating that pulse duration may be critical even at the ultrafast time scales.
Due to the greater visibility of the rings, the 111 silicon wafers were extensively investigated and used for systematic measurements. The Si 100 orientation laser-generated amorphization rings were faint, therefore difficult to accurately measure. Only the outer ring radius was measured for the smallest beam diameter and is depicted by squares in Figure 3b,e,h. The Si 100 and Si 111 orientations' outer ring radius were similar, the Si 100 results being slightly scattered due to the poor contrast of the rings, leading to errors in measurement.
For each type of silicon doping and beam diameter, we analyzed the dependence of the inner and outer ring radii on the beam's maximum fluence. The results for silicon without impurities are shown in Figure 3a-c, for n-type silicon in Figure 3d-f and for (p-type) silicon in Figure 3g-i. Figures 3a,d,g show that for all types of silicon and all focusing conditions, there was a fluence threshold below which the inner ring was not observed. The same conclusion followed for the outer ring (see Figure 3b,e,h. To determine the threshold thr. for the inner and outer radii, we used the empirical formula where w is the laser beam radius on the surface of a crystal, r i is the inner radius of the ring, r o is the outer radius of the ring, and F is the fluence of a laser pulse at the beam center. The solid lines in Figure 3

Theory
The formation of the ring-shaped structure on the crystal surface is directly related to the amount of absorbed energy of the incident pulse and how this energy is distributed inside the crystal. Our hypothesis is that the absorbed energy density value ρ E (r, z), where z is the distance from the crystal surface and r is the distance from the laser beam center (cylindrical symmetry is assumed), determines the outer and inner ring radii essentially. As a result, the energy density must have two critical values, ρ E1 and ρ E2 , which correspond to the outer and inner radii of the ring on the surface (z = 0): Here, ρ E1 is the energy density needed to raise the local temperature of the silicon from room temperature up to the melting temperature T m , and ρ E2 = ρ E1 + ρ EL , where ρ EL is the latent energy density needed to transform solid silicon to the liquid phase. It is well-known from elementary thermodynamics that when a quantity of heat dQ is applied to a body of mass m, its temperature changes according to the law mCdT = dQ, where C is the specific heat of that material. Taking into account that m = ρV, we rewrite our equation in the form ρCdT = dQ/V ≡ dρ E . From this expression, we can find the critical energy density corresponding to the melting temperature where ρ is the density, and T 0 is the initial temperature of the silicon. Given that the isobaric and isochoric specific heats of solids differ only slightly, we assume that C = c p (T)/M, where M is the molar mass and c p (T) is the molar isobaric heat of silicon (see Table 2). Evaluating the integral, we find ρ E1 = 2.987 nJ/µm 3 . If the silicon's molar latent heat is denoted as L (J/mol), then the corresponding energy density is ρ EL = ρ L/M = 4.171 nJ/µm 3 , and consequently, ρ E2 = 7.158 nJ/µm 3 . In order to estimate the amount of absorbed energy density in the crystal, we considered a simplified model of the interaction between laser radiation and matter [25]. Given the short pulse duration, we accounted only for single and two-photon absorption, free carrier absorption (FCA), carrier relaxation, and Auger recombination but ignored heat transport. The characteristic heat transport distance per pulse duration can be estimated using the expression l = (Kτ 1/2 /(ρC)) 1/2 which is derived from the fundamental solution of heat diffusion equation. For silicon, the typical value of the thermal conductivity is K = 1.3 W/(cm·K) and therefore, l ≈ 4 nm which is an order 10 4 times smaller than the laser beam diameter on the crystal surface or pulse propagation distance inside the crystal. In such circumstances, the effect of heat diffusion is negligible and, to a first approximation, can be ignored. The system of equations describing the silicon-radiation interaction is ∂N ∂t Here, I(r, z, t) denotes the intensity of the laser pulse, N(r, z, t) is the free electron concentration, β is the two-photon absorption coefficient,hω is the photon energy, τ R is the relaxation time constant of the electron concentration, and γ is the coefficient for Auger recombination. Absorption coefficient α accounts for the intrinsic silicon absorption and the absorption due to free electrons: where m * e,opt. is the joint optical effective mass of an electron in a crystal, 1/m * e,opt. = 1/m * e,c + 1/m * h,c , τ c = 1/ν c is the collision period of electrons and holes with the lattice. For silicon, the exact value of the collision frequency ν c is not known, but it could be estimated from expression ν c = AT, where proportional constant A = 3.9 × 10 11 s −1 K −1 [25]. From this expression, it follows that the collision period varies from τ c = 33 fs for 300 K down to 8.3 fs for T = T m . In our model, we used τ c as a free parameter. The coefficient α SPA = α(0) accounts only for single-photon absorption due to bound electrons. The boundary conditions are where N 0 = 10 10 cm −3 is the free electron concentration at room temperature, w is the beam radius on the crystal surface, τ p = τ 1/2 / √ 2 ln 2, I 0 = ( √ 2/π/τ p )F 0 , F 0 is the laser pulse fluence, and R is the reflectivity under normal incidence With the solutions of Equations (4) and (5), we can calculate the energy density absorbed in silicon using the expression During the simulation, we varied three parameters: the electron and lattice collision period τ c , the electron concentration relaxation time τ R , and the two-photon absorption coefficient β. The best agreement with the experiment was obtained when τ c = 64.6 fs, τ R = 161 fs, and β = 18.0 cm/GW. In particular, on the crystal surface for this set of parameters, we obtained early calculated absorbed energy densities ρ E1 and ρ E2 for fluences F 1 = 0.260 J/cm 2 and F 2 = 0.408 J/cm 2 , respectively, which were in very good agreement with the experimental values of threshold fluences (see Figure 4a-c). An interesting feature of the model was that it predicted a relatively sharp outer edge of the ring-shaped structure (see Figure 4d). It was directly related to critical electron concentration N c which could be determined from equation Re(ñ 2 ) = 0. In our case, N c = 0.0409 × N at , where N at = 4.99 × 10 22 cm −3 was the atom concentration. Our simulations showed that as soon as the electron concentration approached the critical value N c (Figure 4g), and there was enough pulse energy, the effective absorption coefficient α started to grow rapidly ( Figure 4e) and the absorption rate increased significantly (Figure 4f). The dependence of the absorbed energy density on the crystal depth over the interaction time (about 1 ps) is shown in Figure 5. The process of crystal amorphization (ρ E > ρ E1 ) began when the energy fluence exceeded the value F 1 = 0.26 J/cm 2 . It is important to note that the amorphization depth found from the graph represented only the lower value of the real depth. For example, the model predicted that the amorphization depth was 23 nm for the fluence F = 0.32 J/cm 2 (see the line marked by the symbol), whereas the experiment showed it to be around 45 nm [19]. One reason for this discrepancy is the fact that we did not take into account the heat transport after the interaction. Knowing that the laser beam size on the surface is of the order of µm, and the absorbed energy depth is of the order of nm, due to the gradient of heat, the heat transport is mostly in the depth direction. The generation of shock waves might also help to deepen the amorphization layer [26]. Another and potentially more important reason is that the plasma gets very reflective and changes the group velocity of the pulse when it reaches close to the critical plasma density, which also may change the depth of the amorphization.
Despite a lot of postpulse effects emitting drawbacks in our numerical analysis, we showed that this simplified numerical model which did not demand a high computation power worked surprisingly well on the analysis of the surface modification of silicon. For example, we found experimentally that the surface ablation threshold for undoped silicon was about 0.62 J/cm 2 . For such fluence, the absorbed energy density on the surface was very close to the limiting value ρ E3 = 11.334 nJ/µm 3 required to reach the boiling point of silicon (see Figure 5, line marked by symbol).
The simplified numerical model suggested that a ring-shaped structure on the silicon surface was formed due to a change in the aggregate state of the material. From a thermodynamic point, it is safe to say that when r > r o , the crystal surface remains unaffected. In the region r o < r < r i , one has a transient solid-liquid state where T = T m . In that region, bonds between atoms begin to break and the formation of randomly oriented and partially oxidized microcrystals is possible. The random orientation of microcrystals can explain the strong light scattering from the surface. In the region r < r i , all silicon on the surface should be liquid for some time. It is very difficult to expect that silicon in this region can recrystallize into the former lattice without any consequences. The formation of the oxide layer or modification of the native oxide layer is not excluded, and it can form on the silicon surface, especially when it is heated and in a liquid state [11]. However, other experiments in inert gases or liquids also show the same ring formation [27]. When laser energy absorption is high enough, i.e., ρ E > 11.3 nJ/µm 3 , the boiling of silicon can be reached, and this can be associated with the ablation threshold of the material. ρ E1 = 2.987 nJ/µm 3 (requires to reach the melting point), ρ E2 = 7.158 nJ/µm 3 (requires to transform solid silicon into liquid phase), ρ E3 = 11.334 µm 3 (requires to reach the boiling point).

Summary and Conclusions
We showed the formation of rings on the surface of silicon with various beam diameters and provided a simplified numerical tool that accurately predicted the ring formation when the latent heat was included. We experimentally demonstrated that the laser fluence threshold of the outer ring and inner ring does not depend on the beam diameter, material doping, or crystal orientation, when it is irradiated with a 1030 nm wavelength and 300 fs duration laser pulses. Surprisingly, the numerically estimated depth of the amorphous layer was similar to what other research groups measured. This simplified numerical tool can be used to further analyze the amorphization of the silicon surface with altered beam shapes and/or pulse duration. It was shown that the depth of the amorphization was clearly dependent on the laser energy absorption and heat; we think the use of the laser burst mode (a couple of pulses that are separated in time with delays much shorter than the repetition rate of the laser) could increase the depth of the heat transfer and increase the depth of the amorphous layer.  Data Availability Statement: Any further details relevant to this study may be obtained from the authors upon a reasonable request.