Sub-surface modifications in silicon with ultra-short pulsed lasers above 2 µm

Nonlinear optical phenomena in silicon such as self-focusing and multi-photon absorption are strongly dependent on the wavelength, energy, and duration of the exciting pulse, especially for wavelengths > 2 µ m. We investigate the sub-surface modiﬁcation of silicon using ultra-short pulsed lasers at wavelengths in the range of 1950–2400 nm, at a pulse duration between 2 and 10 ps and pulse energy varying from 1 µ J to 1 mJ. We perform numerical simulations and experiments using ﬁber-based lasers built in-house that operate in this wavelength range for the surface and sub-surface processing of Si-wafers. The results are compared to the literature data at 1550 nm. Due to a dip in the nonlinear absorption spectrum and a peak in the spectrum of the third-order nonlinearity, the wavelengths between 2000 and 2200 nm prove to be more favorable for creating sub-surface modiﬁcations in silicon. This is the case even though those wavelengths do not allow as tight focusing as those at 1550 nm. This is compensated for by an increased self-focusing due to the nonlinear Kerr-effect around 2100 nm at high light intensities, characteristic forultra-shortpulses.


INTRODUCTION AND STATE OF THE ART
Silicon wafers are currently largely separated from the bulk mono-crystalline Si-block using thin diamond saws, which introduce a loss of material of up to 50% [1]. For extremely thin wafers with a thickness ≤ 100 µm, the loss increases to 70% [1]. Alternative technologies of separating the chiplets from the wafer aimed at lowering the cost and losses, as for example, stealth dicing [2,3], were introduced several years ago, but on one hand, the technology is not perfect, and on the other hand, to the best of our knowledge, there is no such technique available for separating the wafer plates from the bulk crystal itself (similar to stealth dicing, when separating the chiplets from each other).
Thus, alternative methods for wafer separation have been in development, such as epitaxial Si lift-off, stress-induced spalling, and smart-cut [4,5]. This latter technique employs the fact that by the introduction of defects in a target layer below the Si wafer surface, this layer will be weakened, allowing the wafer itself to be removed. There are multiple techniques for introducing such defects into bulk Si. The laser-based technique is only one of them. The expansion coefficient difference for two materials can be used [6] without adding modifications inside the crystal. However, the crystal itself can be modified internally by the use of a layer of porous Si [7], the use of high-energy protons [5], or the use of two-photon absorption of laser pulses in the transparency wavelength range of silicon, e.g., at 1550 nm [8]. Compared to those results, based on the study we already published [9], we suggest to move to longer wavelengths in the three-photon-absorption range.
Even though much work has been conducted at common wavelengths such as 1064 nm, 1200-1300 nm, and 1550 nm, both experimentally and numerically [2,[9][10][11][12][13][14][15][16][17], only a small amount of work has been done at longer wavelengths, including 1965 nm [9], 1970 nm [18], and 2300 nm [9,19], where threephoton absorption comes into play, promising to produce even better results than at 1550-1960 nm [9]. Chambonneau et al. [18] used a Tm-fiber laser for generating pulses with a duration Research Article of 2 ps, 9 ps, and 400 ps. Furthermore, they investigated the damage threshold necessary to modify silicon, depending on the pulse duration. They showed an increased threshold for femtosecond pulses compared to picosecond pulses, with the lowest damage threshold of 0.2 µJ being around 9 ps. In the latter work from Nejadmalayeri et al. [19], a femtosecond optical parametric oscillator (OPO) was used to produce modifications below a SiO 2 layer, but not within bulk Si. In our initial study [9], we have noticed that around 2000-2200 nm due to the wavelength dependence of several nonlinear optical phenomena in silicon an ultra-short pulse propagating through the material should be more efficient in modifying the material compared to longer and shorter wavelengths.
The aim of the present work was to carry out numerical investigations that study this phenomenon in more detail by taking into account different pulse energies and additional wavelengths for the numerical simulations. Furthermore, we conducted experiments to support the numerical simulations and to compare the results. For that, we applied lasers operating at different wavelengths around 2100 nm (an optimum wavelength according to Richter et al. [9]) for generating sub-surface modifications in silicon, and compared the results with numerical simulations.
The paper is structured as follows. In the first two sections, we numerically investigate the influence of the wavelength, pulse energy, and pulse duration on the defect formation of silicon using both numerical simulations and experiments for wavelengths between 1950 and 2400 nm. The results obtained are compared with the existing data for 1550 nm, 1970 nm, and 2300 nm. Thereby we can confirm the preliminary results we presented in Ref. [9], that an optimal wavelength range for silicon processing exists.
On the experimental side, this paper is separated into two parts. Several different lasers have been used to modify a Sicrystal, both on the surface and within the bulk material by using wavelengths longer than 1950 nm. In order to optimize the process in both quality and speed, several parameters such as pulse energy, pulse duration, and wavelength had to be optimized. It has been already demonstrated that shorter pulses provide better cut quality and precision compared to longer pulses [20,21], while lowering the ablation threshold. Due to the shorter pulse duration (and resulting higher peak intensities), the material around the target area has less time to react to the incoming pulse and to warm up, which in turn results in more localized modifications only at the targeted spot. We have previously shown that longer wavelengths in the mid-IR range between 2000 and 2200 nm are beneficial for processing of bulk Si [9,15]. This is due to the combination of lower nonlinear absorption [22] and higher self-focusing [23]. The lower nonlinear absorption leads to less absorption of the pulse energy while propagating through the material, except at the focal spot. At this point, the higher self-focusing results in a very localized intensity peak, which then in turn triggers nonlinear absorption and creates localized modifications. Therefore, we present results showing the creation of sub-surface modifications within silicon by using pulses at nano-joule levels at two different wavelengths (1965 nm and 2350 nm) in the first sub-part of the experimental section of this paper. Afterwards, we present the generation of sub-surface modifications in silicon using micro-joule level pulses at a wavelength of 2090 nm in the second sub-part.
In summary, we present here for the first time to our knowledge the results of both numerical and experimental studies above 1550 nm using fiber based lasers with wavelengths of 1965 nm, 2090 nm, and 2350 nm at different pulse energies. We use a simple numerical model based on the nonlinear Schrödinger equation. The results show not only that silicon processing is favorable at longer wavelengths (as was shown in previous studies [9,18]), but also that an optimum wavelength range for modifying silicon exists. We explain this by the pronounced wavelength dependencies of the nonlinear refractive index n 2 and the multi-photon absorption coefficients β and γ for two-and three-photon absorption, respectively.

METHODOLOGY
The spectrum of the nonlinear multi-photon absorption calculated (in our case) as a sum of two-photon and three-photon absorption has a "dip" within the wavelength range we target (between 2000 and 2200 nm), as shown in Ref. [9,15]. The values for two-photon absorption decrease above ≈ 1800 nm, while the three-photon absorption is noticeably rising above 2300 nm, thus leaving a gap in between [9,15]. These values are shown in Fig. 1. Here we show the values for the twophoton absorption coefficient β in Fig. 1(a), the three-photon absorption coefficient γ in Fig. 1(b), and two different sets of measurements for the nonlinear Kerr coefficient in Figs. 1(c) and 1(d), depending on the wavelength λ. While a continuous function can be given for β and γ , only discrete values can be given for n 2 , with errors up to 30% within each measurement [22,23]. The difference of the measured Kerr coefficients can be up to five times in absolute values. The exact values of the Kerr coefficient depend strongly on material properties (e.g., singlecrystal versus polycrystalline material, n-doped versus p-doped). Therefore, we decided to include two different values for the Kerr coefficient, to cover a larger range of possible material configurations.
The nonlinear figure of merit (NFoM) characterizing the relative merit of the Kerr nonlinearity coefficient versus the two-photon nonlinearity is usually defined as the ratio between both values, divided by the optical wavelength in vacuum [26]. Extending this expression to multi-photon absorption results in the NFoM written in the following form [23]: with λ the wavelength in vacuum, n 2 (λ) the wavelengthdependent refractive index, I the field intensity, and β (K ) the coefficient for K -photon absorption. Based on Eq. (1), an optimal wavelength range for generating sub-surface modifications could be estimated to be in the range of 2000-2300 nm [9]. The exact wavelength depends on the exact value of the material parameters involved, such as the value for the wavelength-dependent Kerr nonlinearity n 2 (λ), which again depends strongly on the material properties, as mentioned above. Thus, one can only predict a certain wavelength range, Fig. 1. Values for the nonlinear parameters for (a) two-photon absorption β(λ) [24], (b) three-photon absorption γ (λ) [25], and (c), (d) two different Kerr-value levels n 2 (λ) [22,23], plotted against the wavelength. The Kerr coefficients are shown together with their respective measurement errors.
not an exact value of the optimal wavelength for generating sub-surface modifications.
A peak of the NFoM similar to the main peak in the range of 2000-2200 nm can be calculated around 3200 nm, at the transition between three-photon absorption and four-photon absorption. Still, according to Zavedeev et al. [15], the deposited energy is ≈ 6 times higher for a pulse in the wavelength range of 2000-2200 nm compared to the longer wavelengths, as shown in Fig. 2. This confirms the higher efficiency of the wavelength range around 2100 nm compared to longer and shorter wavelengths.
To model for the propagation of the pulse through the material, we use the nonlinear Schrödinger equation [15,27,28] Equation (2) describes the propagation of a time (t)dependent and transversely inhomogeneous laser field E through bulk Si along the z axis. The diffraction is taken into account by the transverse Laplacian, where k 0 is the wave number. The nonlinear effects of self-phase modulation and self-focusing are described by the term proportional to n 2 , and τ c = 3.5 fs is the free-carrier momentum scattering time [15]. The nonlinear multi-photon absorption for K -photon absorption (β (K ) ), the absorption caused by free carriers with the concentration N, and the inverse bremsstrahlung absorption σ are taken into account, as well. Since the investigated wavelengths are in the range where the single-photon absorption and higher-photon absorption can be neglected [29], we have Research Article considered the action of only two-and three-photon absorption. The numerical simulation was done while using support functions from the library deal.II [30].
The paraxial approximation introduces several simplifications into the beam propagation equation, but is valid only until a certain beam diameter [31] and a maximum NA of the lens of 0.2 [32]. Thus, if the beam is focused too strongly because of the optical Kerr effect, Eq. (2), which uses that approximation, loses its validity. This effect can be seen for Figs. 4(a), 4(c), 5(a), and 5(c). A high value for n 2 will result in a collapse of the beam, which in turn results in oscillations of the intensity and a beam diameter below the accuracy of the paraxial approximation. Deviations from paraxiality also have to be taken into account when simulating lenses with a NA of ≥ 0.2 [32].
There are different approaches for calculating the response of the material if irradiated with ultra-short pulses, which are extensively discussed in Ref. [33]. Examples include the two-temperature model (TTM) and its extension, the threetemperature model [34], a hydrodynamic approach, and molecular dynamics (MD) simulations. Furthermore, hybrid models such as a combination of a TTM and a MD are suggested [33]. In order to simplify calculations, in this study we take into account only the generated carrier density. This value can be calculated by using [15,27,28] where |E | 2 ≡ I . Recombination of the free carriers was not taken into account due to the long lifetimes in Si, up to ≈ 10 ns [10]. Furthermore, we decided to neglect the impact ionization coefficient. The coefficient cannot exceed 3.6 × 10 10 s −1 [15], and thereby it is significant only for pulse durations ≥ 10 ps.

NUMERICAL SIMULATIONS
An initial rough comparison of the action produced on silicon by ultra-short pulses at different wavelengths can be done using Eq. (2), by comparing the transported energy towards the focal spot. The results of such simulations for the wavelengths 1550 nm, 1950 nm, 2150 nm, and 2350 nm with pulse energies of 1 µJ, 10 µJ, and 100 µJ for two different values of n 2 [22,23] are shown in Fig. 3. The results denoted with "high n 2 -values" and "low n 2 -values" use n 2 -values taken from Refs. [23] [shown in Fig. 1(c)] and [22] [shown in Fig. 1(d)], respectively. The pulse duration is set to 5 ps, to correlate as close as possible with the experiments presented in Section 4. In these experiments, a large range of different lenses and focal depths was used. To limit the presented data to a reasonable amount, while still being able to demonstrate representative results, in the simulations, we selected a lens with a focal length of 8 mm. The beam radius was adjusted to have a comparable focal spot area for linear focusing, with a beam radius of 1 mm for the wavelength of 1550 nm. To achieve a similar fluence for the other wavelengths used in the simulations, the beam radius r λ 1 for a wavelength λ 1 was adjusted using This resulted in a similar fluence level at the focal spot for each wavelength, regardless of the used wavelength. On the other hand, this beam radius modification results in NAs, which in turn influence the level of the self-focusing. A smaller NA will result in a higher intensity closer towards the lens, thereby triggering the self-focusing and nonlinear absorption earlier. Therefore, we have to investigate the influence of changing the NA in a future investigation. While the pulse with the shortest wavelength is absorbed quite significantly already on the way to the focal spot, the other wavelengths can retain more energy, before depositing it closer to the focal spot. In addition, the absorption close to reaching the focal spot depends on the value of the third-order nonlinearity, as shown, e.g., when comparing Figs. 3(a) and 3(b). Here one also can see the ripples originating from the strong self-focusing within the material, followed by local absorption and collapse, which show up in Figs. 4 and 5, too.
In Figs. 3(a)-3(f ), most of the pulse energy is absorbed already in the first millimeter of the material, leading to a loss of 50%-99% of the pulse energy, depending on the initial pulse energy and wavelength. This process is independent of the value of the Kerr nonlinearity. This high loss also means that the pulse has ≈ 0.5 µJ left for 1 µJ initial pulse energy, ≈ 0.9 µJ is left for 10 µJ initial pulse energy, and ≈ 0.9 µJ is left for 100 µJ initial pulse energy, which in turn shows that an increase in pulse energy does not necessarily increase the energy transported to the focal spot, but instead increases the chance of damaging material on the way to the focal spot. Due to the similarity of the transported energy for the pulse with 10 µJ and with 100 µJ, the profile for the generated intensity and carriers are very similar. Thus, those figures are omitted in Figs. 4 and 5.
When considering the achieved intensities, one can see in Fig. 4 both the effect of the nonlinear refractive index and the collapse of the pulse.
In Fig. 4, the self-focusing of the pulse at wavelengths of 1950 nm and 2150 nm starts at a significantly smaller depth than for the other two wavelengths. Self-focusing at a smaller depth leads to a modification zone that is spread out much more, compared to a pulse focused at a larger depth. Nevertheless, it also shows the collapse of the beam in Figs. 4(a) and 4(c), which breaks the paraxial approximation and leads to non-physical results. This breakdown results in oscillations of the beam. For a lower Kerr coefficient [as shown in Figs. 4(b) and 4(d)], the created plasma and the nonlinear absorption stop that breakdown before it can occur. Therefore, for larger Kerr coefficients, a non-paraxial model is required.
The achieved intensity for the shortest wavelength (1550 nm) is higher than for other wavelengths when assuming the lowest reported values for the Kerr effect in silicon, even though-as shown in Fig. 3-a much lower portion of the pulse energy is transported to the focal spot for this wavelength. This can be explained by the interplay of stronger focusing at shorter wavelengths, something compensated for at longer wavelengths with an increasing Kerr effect around 2000-2200 nm and thus even tighter focusing.
Finally, the effect already manifested in Fig. 3 also shows up in the intensity curves: since most of the energy of the pulse is absorbed in the first small part of the beam path, the resulting energy at the focal spot is approximately the same for 10 µJ and 100 µJ, being only slightly higher than that for 1 µJ. This results in the same intensity levels for the former two, and only slightly lower intensity levels for the latter one.
Higher intensity also means higher absorption, which can be calculated using [15] Here the absorption depends not only on the nonlinear absorption, but also on the free-carrier density, which in turn is increased by the pulse itself [see Eq. (3)]. This generated free-carrier density is shown in Fig. 5, for both low and high n 2 -values. Again, the pulses at 1950 nm and 2150 nm can generate more free carriers, compared to the other two wavelengths, which in turn increases absorption and diffraction. This also holds true when considering the lower Kerr nonlinearity. Even though the shortest wavelength could generate the highest intensity, the amount of the generated free carriers is approximately the same for all wavelengths, because of the higher free-carrier-generation efficiency at longer wavelengths. As mentioned above for Figs. 4(a) and 4(c), we still have to consider the breakdown of the paraxial approximation for Figs. 5(a) and 5(c), and thereby have to consider the results as not physically relevant. Nevertheless, we still can use the results for the lower Kerr coefficients. Here we can see an elongated structure within the material, which closely corresponds to the structures observed in Fig. 12.
Moreover, because of the similar intensities for all pulse energies, the generated free-carrier densities are also at approximately the same level, thereby confirming again that increasing the pulse energy from 10 to 100 µJ does not increase modifications at the focal spot, but rather hinders them by modifying the material on the way there. This will result in modified material at places along the beam axis that are not targeted, thereby decreasing the value of the generated cut.
When extending the equation responsible for the reaction of the material [here Eq. (3)] towards models such as TTMs (as for example in Refs. [14,16]), the free carrier density also plays an important role by coupling the energy of the generated free carriers to the lattice itself by the coupling constant [14,35] with τ c = 240 fs the carrier independent relaxation time [35], N the carrier density, and N crit = 6 × 10 26 m −3 the critical density. The values are valid for N ≤ 2 × 10 26 m −3 , a value we do not reach in our calculations. Thus, a higher free-carrier density makes the energy transfer from the excited carriers to the lattice more efficient. This assumption must still be investigated further by coupling the pulse propagation to TTM. Finally, two different processes can lead to melting [36,37]. On one hand, the carrier density can be increased up to a critical carrier density N c ,melt [16,38]. At that density, enough free carriers are excited that the material can show non-thermal melting [21,37,39,40]. This critical density (n c ≈ 1 × 10 22 cm −3 ≡ 1 × 10 28 m −3 [38]) is reached when approximately 10%-20% of the valence electrons are excited. This state shows a behavior similar to melting, while the lattice temperature is still below the melting temperature. The other process is conventional thermal melting, when the lattice reaches a temperature above the melting temperature. This critical density is wavelength independent, and is not to be confused with the critical density in a plasma, which we can define as This density describes the density of a plasma at which it becomes reflective for the incoming laser light at a wavelength ≥ λ 0 , thereby blocking it from propagating further. Contrary to that, N c ,melt is not depending on the wavelength. Still, there are several different possibilities for determining the modification threshold inside bulk Si, and this critical density is only one of them, as shown in Ref. [41]. Due to the long upper lifetime of free carriers (up to 10 ns [15]), these models become of interest when going to simulation times in the nanosecond range, after there will not be a significant heat transfer before that. For shorter simulations, especially in the lower pico-and femtosecond ranges, Eq. (3) can already give a good estimation of the response of the material. Based on the numerical simulation data shown in Figs. 3-5, we can conclude that the created modifications in Si using the lens with f = 8 mm and f i = 4 mm have a thermal nature, if we use the criterion for non-thermal melting as defined above, since the generated free-carrier density did not exceed the critical density N c ,text . Our numerical simulations also confirm the initial assumption, that the wavelengths in the range of 2000-2200 nm are more efficient for modifying Si compared to wavelengths above 2200 nm and below 1900 nm.
We also noted that the plasma contributed a significantly lower amount to the pulse propagation compared to the twoand three-photon absorption. Even though the plasma intensity could reach values at ≈ 1 × 10 26 m −3 , the absorption of the light due to multi-photon absorption was several orders of magnitude higher. Therefore, the influence of plasma absorption and de-focusing contributed significantly less to the beam propagation compared to the nonlinear absorption processes and can be neglected when considering the arrest of the collapse of the beam. Nevertheless, for high values of the Kerr coefficient, this did not stop the collapse of the beam, as shown in Figs. 4(a), 4(c), 5(a), and 5(c) compared to lower values, as displayed in Figs. 4(b), 4(d), 5(b), and 5(d).

A. Summary
We investigated the wavelength dependence of different nonlinear effects acting on a pulse propagating through silicon using a simple numerical model. Furthermore, we observed the correlation between the reached peak intensity, the amount of generated free carriers, and the initial pulse energy. In our modeling, we used two values of Kerr coefficients reported in the literature, "high values," shown in Fig. 1(c) based on the results from Wang et al. [23], and "low values," displayed in Fig. 1(d) based on the results from Lin et al. [22]. For the lower Kerr coefficient, we obtained a prolonged shape for both the intensity values and the generated carrier density, similar to the structures shown in Fig. 12, due to a balance between the self-focusing and nonlinear absorption of the pulse. For the higher Kerr coefficient, we could show a collapse of the pulse and oscillations of both the intensity values and the generated carrier density, Research Article especially for 1950 nm and 2150 nm, indicating a breakdown of the paraxial approximation and therefore non-valid results. Nevertheless, we can say that significantly lower pulse energies at 1950 nm and 2150 nm are necessary for reaching a similar effect compared to the longer and shorter wavelengths, when considering the higher Kerr coefficients. For the lower Kerr coefficients, we could not detect significant differences (≥ 10%) between the different wavelengths.
We also showed that the increased laser pulse energy does not necessarily result in the higher pulse energy reaching the focal spot inside the silicon, particularly for pulse energies above 10 µJ. The excess energy will be absorbed while propagating through the material towards the focal spot, increasing the probability to induce damage both on the surface and in the region between the surface and the focal spot. Finally, we have shown the importance of balancing the high nonlinearity of silicon in the region between 2000 and 2200 nm by careful choice of the laser pulse energy, since the excess energy at the conditions of high nonlinearity might cause the self-focusing to happen too far away from the focal point, thus leading to pulse collapse. In this case, the area with peak light intensity ("modification zone") will be prolonged along the pulse propagation axis [Figs. 4(a) and 4(c)], which is not beneficial for microprocessing.
The results from Zavedeev et al. [15] suggest a critical role of plasma defocusing, which we could not confirm in our simulations. We cannot draw an exact comparison, because we were using a relatively low NA of 0.125 to stay below the limitations of the paraxial approximation limited to an NA of <0.2 [32], while Zavedeev et al. used an NA of 0.3 [15]. The effects at high NA require careful non-paraxial simulations, which will be a topic of a separate study.

A. Laser Sources
In order to get experimental data that can be used as verification of the numerical results, modifications on and below the surface of Si-wafers were made using three different home-built lasers based on specialty fibers: an all-fiber large mode area (LMA)-based master oscillator power amplifier (MOPA) system based on Tm operating at a central wavelength of λ = 1965 nm, τ p = 2 ps, and W max = 980 nJ (see Section 4.A.1), a Ho-fiberbased CPA laser operating at λ = 2090 nm, τ p = 5 ps, and W max = 760 µJ (see Section 4.A.2), and an Er-fiber-based Cr:ZnS laser operating at λ = 2350 nm, τ p = 1.7 ps, and W max ≈ 20 nJ (see Section 4.A.3). The fiber MOPA (labeled laser "A" in this paper) consists of an oscillator, mode-locked by a semiconductor saturable absorber (SESAM), and two amplification stages. Here the oscillator generates laser pulses at a wavelength of λ = 1965 nm with a repetition rate of f = 7.6 MHz, pulse duration of τ p = 2 ps, and pulse energy of W = 0.1 nJ. By using a pre-amplifier, this pulse energy is amplified to W = 0.5 nJ. Finally, the pulse is amplified in a second amplifier to a pulse energy of W = 980 nJ [42]. For experiments at 2090 nm, we used a compact ultra-short pulsed-fiber-laser-based Ho:YAG amplifier system from ATLA Lasers AS (labeled laser "B" in this paper). At 10 kHz repetition rate, the system delivered up to 760 µJ pulses with 5 ps duration from the amplifier stage.
3. Er-Fiber-Based Cr:ZnS Laser at λ = 2350 nm The Er-based fiber laser (labeled laser "C" in this paper) was a mode-locked Cr:ZnS oscillator from ATLA Lasers AS working in the normal dispersion regime [43], pumped by the Er:fiber laser. The combination of dispersion-managed mirrors and material dispersion compensation in the cavity ensured the normal-dispersion operation regime of the laser. The laser generated linearly chirped pulses centered at 2350 nm with a duration of about τ p = 1.7 ps and pulse energy up to 48 nJ at a pulse repetition rate of 3.8 MHz. Around 50% of the pulse energy was lost in the isolator, beam delivery optics, and focusing objective, resulting in slightly above 20 nJ delivered to a surface of a silicon sample.

B. Experimental Setup
All lasers were coupled to a processing stage. The light was focused using several lenses with focal distances between 0.9 and 11 mm either onto the surface of the sample (as in Fig. 7) or inside the sample [as, e.g., in Fig. 10(b)]. The typical beam diameter before the focusing optics was ≈ 2 − 4 mm. The sample itself was mounted onto a four-axis translation stage, with two axes controlled by an electronic controller. Electronic control allowed precise movement with a fixed speed, variable between several µm s −1 and several mm s −1 . The samples were mono-crystalline silicon wafers with a thickness of 0.6-2 mm, a surface orientation of (111), and resistivity of 10 cm −1 . For processing, the sample was moved either parallel to the incoming beam [as shown in Fig. 6(b)] or perpendicular to it [as shown in Fig. 6(a)].

C. Investigation of the Processed Samples
The processed samples were investigated using three different techniques, including optical microscopy, transmission IR microscopy, and scanning electron microscopy (SEM). The first two techniques are non-destructive, while the cross-sectional imaging using the SEM required splitting of the sample. Optical microscopy alone is applicable to study of the surface, while transmission IR microscopy was used to look through the sample and investigate sub-surface modifications without damaging the sample. As mentioned in Section 1, both surface and subsurface silicon modifications were created using pulses at both nJ and µJ pulse energy levels. In Section 4.C.1, the modifications created by pulses at sub-µJ pulse energy levels are investigated. In Section 4.C.2, we investigate the modifications created by pulses with energies above 1 µJ.

Silicon Modifications Created at nJ Pulse Energy Levels
We used the Tm-fiber-MOPA (laser A), operating at f = 7.6 MHz and W = 28 nJ, and the Er-based fiber laser, operating at λ = 2350 nm, f = 3.7 MHz, and W = 17 nJ (laser C), to create modifications both on the surface and below the surface at nJ pulse energy levels. As an initial test, we created lines on the surface of the silicon sample, which could be observed directly with an optical microscope. Those lines are shown in Fig. 7.
As one can see in Fig. 7, the cut done at the longer wavelength [i.e., Fig. 7(b)] is significantly smaller in diameter (<2 µm compared to 3 µm) and cleaner with less residual around the cut, thus suggesting that using a longer wavelength might be beneficial for cutting Si.
This size reduction could originate from a smaller spot size. Assuming that the beam was focused perfectly onto the surface, the focal spot size for laser A was 4.25 µm, while for laser C, the spot size was 1.6 µm. Furthermore, we assumed that the absorption for laser A was based purely on two-photon absorption, and the absorption from laser C was based purely on three-photon absorption. This means that the effective spot size (if the threshold is reached for the whole Gaussian pulse) is 3 µm for laser A and 0.92 µm for laser C. This corresponds roughly to the experiment.
After we confirmed that we could create modifications at the surface of silicon, we also conducted sub-surface modifications in silicon. For this task, the same lasers as mentioned above were used, with varying pulse energy levels. Those modifications were observed using both a transmission Fourier transform infrared spectroscopy (FTIR) and an IR microscope. The FTIR microscopy was carried out in transmission mode on a Bruker Hyperion 3000 microscope with a 64 × 64 pixel mercurydoped CdTe focal plane array detector interfaced to a Bruker Vertex 70v FTIR spectrometer. A 10× Cassegrain objective was used for imaging. Background spectral images were recorded on an empty sample stage and used to calculate extinction spectra at each pixel after measuring the laser processed Si sample. For spectral acquisition, 100 scans with a spectral resolution of 4 cm −1 were averaged. Spectral images were obtained from two spectral regions. Extinction maps from the region 3.965 − 3.027 µm (2522 − 3304 cm −1 ) were obtained as false-colored maps from integrating extinction without baseline correction. As the system does not show any absorption in this spectral range, extinction here is caused by scattering. Maps were also produced in the spectral region 10.320-7.553 µm (969 − 1324 cm −1 ), containing the absorption peaks of Si-O stretching modes. In this region, a linear baseline was subtracted from the spectra between the left and right ends of the spectral range. In that manner, maps show the distribution of absorption from Si-O modes. Results of those measurements for the laser  Fig. 9, the continuation of the sample is shown as a visible image outside the IR maps. The movement of the sample is perpendicular to the laser beam, as shown in Fig. 6(a). (a) Mapped integrated extinction in the range 3.965-3.027 µm, (b) mapped integrated absorbance in the region 10.320-7.553 µm, and (c) example spectra at the spots marked in the respective color in Fig. 8(a).

Silicon Modifications Created at µJ Pulse Energy Levels
By using laser B operating at λ = 2090 nm, we increased the pulse energy to several µJ. We used this increased pulse energy for creating larger sub-surface modifications, and investigated their shape afterwards. In addition, for faster non-destructive detection of sub-surface modifications in bulk silicon (compared to the FTIR used in Section 4.C.1), we designed a special transmission IR microscope operating at λ = 1300 nm. To increase the energy exposure for every point on the beam axis without changing the writing speed, and also increase visibility of the written structure during the analysis at the IR transmission microscope, the sample was moved parallel to the laser beam (so-called longitudinal writing approach, in contrast to the transverse writing approach demonstrated in Figs. 8 and 9).  Fig. 10(b)]. After having successfully demonstrated volume modification of silicon using the longitudinal writing approach, we applied the transverse approach [as shown in Fig. 6(a)]. The movement speed and laser parameters were similar to those described in Fig. 10.
The results of the experiment are shown in Figs. 11-13. In addition to the surface modification lines visible on the standard microscope image, there are several buried modification lines visible only on the transmission microscope image (Fig. 11).
Each of the detected buried lines was labeled with a number from 1 to 5, and the numbering is kept throughout Figs. 11-13. One can notice that lines 2 to 5 are not continuous, indicating that the processing parameters are rather close to the modification threshold.
The profile of the buried modifications can be seen on the cross-section transmission microscopy view of the sample (Fig. 12). It must be noted that lines 1, 3, and 5 are stacked lines, i.e., several lines drawn on top of each other; however, lines 2 and 4 are single lines. Each "layer" of the stacked line is labeled with an index letter. From the cross-section view, one can clearly see that each single modification has a triangular (conical in volume) elongated shape with the sharp corner pointing out the focal point of the focusing objective. The fact that modifications started significantly before the focal spot indicates that processing parameters must be optimized. The shapes of the modifications will be further discussed in Section 5.
After the induced buried modifications were confirmed by the IR transmission microscopy, a more traditional SEM technique was implemented to visualize the shape of the modified regions. The thin cross-section slice was cut from the original sample, polished from both sides, and etched for 4.6 h in KOH solution. The SEM image of the cross section is shown in Fig. 13, both front view [Figs. 13(a) and 13(b)] and back view [Figs. 13(c) and 13(d)] of the sample. Similar elongated conical shapes are observed. One can note that some layers and lines are missing, which is the consequence of the dashed (noncontinuous) geometry of the buried structures. The observed lines were identified in correspondence with Fig. 12 and labeled with corresponding line numbers and index letters.

D. Summary
We showed three different lasers operating at central wavelengths of 1965-2350 nm and two different levels of pulse energies (nJ and µJ) and the resulting modifications of silicon both on the surface and below. The modifications were investigated in more detail for λ = 2090 nm, which we identified as a "sweet spot" for silicon modifications. The sub-surface modifications we generated with this laser could be visualized  The buried modification lines are numbered 1 to 5, corresponding to Fig. 11. Layers of the stacked lines are shown by indices, e.g., 1a, 1b, and 1c.
both by a custom-built transmission IR microscope and a SEM. The structures visible with both methods correlated in size. Sub-wavelength feature size was evidenced in Fig. 7 but could not be observed yet for sub-surface structures as in Figs. 10-13.
A possible reason for that is the nonlinear absorption, which is starting already before the focal spot, thereby creating modifications with sizes larger than the laser wavelength as shown in the numerical simulation in Fig. 4, and in the experiment in Figs. 10-13.

DISCUSSION
As can be seen in Fig. 12 and other figures, the created modifications within Si are not perfectly shaped around the focal spot, but rather elongated and start significantly before the focal spot. One possible explanation for that is that nonlinear absorption happens already before reaching the targeted spot within the material. This means that the critical intensity for material modification is reached already in the prefocal region, thereby starting the modification process before reaching the focal point.
A similar behavior is shown in Fig. 4. One would expect the intensity to rise sharply to one single peak exactly at the focal spot, but instead it stays at the level close to the peak value for a certain distance before and around the focal spot. In that period, the pulse energy is absorbed within the material, and depending on the pulse energy, leads to modifications at points that are not targeted. This modified area length depends on the pulse energy, but also on the value for the nonlinear Kerr value [as can be seen in Fig. 4(a) compared to Fig. 4(b) for 1 µJ, and similarly in Figs. 4(c) and 4(d) for 10 µJ].
In addition to the elongated shapes, one can notice certain asymmetry of the modification traces visible in Figs. 12 and 13. This effect was also observed by Kämmer et al. [44], where the processing of silicon was done at a wavelength of 1550 nm. Several factors might have contributed to this, including purely experimental issues such as laser beam quality, pulse-to-pulse stability, and asymmetry of the focused beam due to imperfect adjustment of the focusing lens, as well as more fundamental aspects such as aberrations caused by the focusing optics (such as spherical aberration and coma), local impurities, and defects both on the surface and in the volume of the silicon wafer, and, finally, strong nonlinearity of the self-focusing process resulting in highly distorted intensity distribution around the focal region. A more detailed analysis is required to identify the most contributing factor.
This also means that a higher Kerr coefficient allows processing of Si at lower pulse energies. Lower energies will result in lower absorption on the way to the focal spot, while the selffocusing will occur later, thereby resulting in a smaller modified area, which in turn explains the sub-surface modifications shown in Figs. 8 and 9. Both modifications were made at energy levels significantly below 1 µJ and are thereby unable to damage the surface, while still modifying material below the surface.
Similar results are also reported by Chambonneau et al. [18], who reported a lower threshold of 300 nJ for creating subsurface modifications in silicon when applying multiple pulses at the same spot. Finally, their reported fluence distribution matches the shape of the lines shown in Fig. 12, which thereby provides experimental confirmation for their simulations with our experimental data.

CONCLUSION
In this paper, we present a numerical model that describes the laser pulse propagation in silicon and gives insights into the dependence of the intensity and distribution of the light field in the volume of silicon on the laser pulse intensity and wavelength. The model is based on the nonlinear Schrödinger equation and accounts for the wavelength dependence of the nonlinear refractive index and multi-photon absorption. Furthermore, to initially confirm the validity of the model, we experimentally demonstrated volume modifications in silicon by ultrashort laser pulses emitting at wavelengths of 1965 nm, 2090 nm, and 2350 nm at two different pulse energy levels (nJ and µJ).
As the main conclusion, we can say that the model predicts that there is an optimal wavelength range around 2000-2200 nm, where subsurface processing of Si occurs at comparatively low pulse energies. By applying the picosecond laser emitting within this wavelength range, we demonstrated the writing of buried modifications inside a silicon wafer using the transverse approach, i.e., when the sample is moved perpendicular to the laser beam direction during the processing. To the best of our knowledge, this is the first experimental demonstration of the transverse-written buried structures in silicon achieved with a picosecond laser in this wavelength range.
Our results show that the wavelength, the value of the Kerr coefficient, the pulse energy, and the pulse duration have a strong impact on the volume structuring of silicon. For nonoptimal pulse parameters, the modified area will significantly extend closer towards the focusing lens. Further research is required to optimize the processing parameters and obtain controllable localized modifications in the focal spot.