1. INTRODUCTION

Recent advances in ultrafast laser material processing enable nanosized structures to be directly fabricated in various materials [1 –4]. The basic processes during femtosecond laser ablation are absorption by electrons, energy transfer to the lattice, and subsequent material removal. These processes can in a first approximation be treated as temporally separated [5]. The first step, the absorption of light, is a complex and interesting issue. In the field of laser nanoprocessing, subwavelength features are created using diffraction-limited laser spots. As the laser pulses are focused into a spot size comparable to the wavelength of the laser light, the interaction between the light and the material takes place in an equally confined volume. Furthermore, due to the high peak intensities involved, nonlinear optical effects play a dominant role. Typically, when a laser pulse propagates through a semiconductor or an insulator, multiphoton absorption takes place during the leading part of the pulse. This leads to the generation of a high concentration of charge carriers. When a laser pulse impinges on a metal, the existing free electrons in the metal will be strongly heated by the leading part of the pulse. In both cases the leading part of the pulse dramatically alters the optical properties of the material, which implies that the trailing part of the laser pulse interacts with a material whose optical properties are locally significantly different from those of the unexcited material. Due to strong focusing of the laser and the inherent nonlinearity of the processes, the changes in the optical properties actually occur on length scales below the diffraction limit. A detailed numerical modeling of this process may serve to provide insight into the complex mechanism of energy deposition under these conditions and is therefore crucial to describe laser nanoprocessing.

In earlier one-dimensional (1D) models of the absorption of femtosecond laser pulses in silicon, the dynamically changing optical properties due to the generation of free carriers were taken into account. The results of these models are in good agreement with reflectivity measurements performed with weakly focused laser beams [6,–8]. However, these 1D models are expected to be inadequate to describe the laser–matter interaction in subwavelength volumes and with the large focusing angles obtained with high-numerical-aperture (NA) objectives. In other studies, the propagation of femtosecond laser pulses in nonlinear media is simulated by solving the nonlinear Schrödinger equation (NLSE) [9,10]. However, in experiments involving subwavelength structures, the NLSE does not hold. For instance, as shown in one of our previous publications [11], the subwavelength plasma formed by the laser is actually a Mie-like scatterer induced inside the silicon-on-insulator (SOI) device layer, which couples incident light in a slab mode. This implies the breakdown of the slowly varying amplitude approximation on which the NLSE is based. Additionally, the high-NA objective used to focus the beam implies the breakdown of the paraxial approximation, another approximation used in the derivation of the NLSE. Due to these limitations of the NLSE, there is an increasing need for the development of numerical models that resolve the full set of Maxwell’s equations coupled to the equations describing the changes in the materials induced by the pulse under tight focusing conditions [12 –15]. However, these studies do not find quantitative agreement with experiments or lack a comparison to experiments. Furthermore, no existing studies address subwavelength laser–matter interactions in metals. Finally, all mentioned models neglect the fact that the effective mass of the charge carriers actually changes due to the increasing temperature during the pulse.

In this paper, we numerically model the first phase of femtosecond laser processing, which is the laser energy deposition. Furthermore, we present self-reflectivity measurements of a strongly focused femtosecond laser beam in single-shot ablation experiments on silicon and gold. The model simulates the propagation of light using a two-dimensional (2D) finite-difference time-domain (FDTD) method, coupled to a set of differential equations that describe the changes in the material driven by the laser light. The model results are compared with the self-reflectivity measurements. We show that the model excellently describes the self-reflectivity measurements on the four types of specimens we investigated, namely, two SOI samples with different device layer thickness, bulk silicon, and gold. We further show that, in the case of strong focusing, a 1D model fails to reproduce the experimental results and that a higher dimensional model is required. Another important result is that it is crucial to take into account that the effective mass of the charge carriers changes during the pulse. As our model excellently agrees with the experimental results, with no open parameters, it can be applied to study and optimize the energy deposition in femtosecond laser nanoprocessing of materials.

The paper is organized as follows. In Section 2 we discuss the theoretical model to describe the propagation of an intense laser pulse in a medium in which the optical properties are changing during the pulse. In Section 3, we explain in detail the numerical implementation of that model. In Section 4, we compare the results of the model to experimental results. A summary and conclusion are presented in Section 5.

2. THEORY

The propagation of electromagnetic waves is governed by Maxwell’s equations,

The above relations are specific to the medium and can be evaluated using microscopic theories. For instance, in a linear and homogeneous medium, we would write $\mathbf{P}={ϵ}_0\chi \mathbf{E}$ and $\mathbf{M}={\chi }_m\mathbf{B}/{\mu }_0$, where $\chi$ is the electric susceptibility and ${\chi }_m$ is the magnetic susceptibility. For electromagnetic waves at optical frequencies, the magnetic susceptibility is in practice negligibly small. We will therefore neglect it in the remainder of this paper and consider only the electric susceptibility. In the simple linear and homogeneous case, the above set of equations can be cast into a wave equation, which can be solved analytically. However, in many cases, the susceptibility is not homogeneous. Then one generally needs numerical methods to solve Maxwell’s equations. Furthermore, when sufficiently strong fields are applied, nonlinear effects can come into play. For example, absorption of light by a semiconductor or dielectric medium can result in the generation of free charge carriers. These free carriers will contribute to the susceptibility. As this will cause the susceptibility to depend on the intensity of light inside the medium, free-carrier generation inherently leads to inhomogeneity even in the case of an initially homogeneous medium. Local heating of the material has a similar effect, as it locally changes the refractive index and thus makes the medium inhomogeneous.

Generally the changes in the local susceptibility are small and slow enough to be negligible when considering the propagation of light. However, when using focused picosecond or femtosecond pulses, dramatic changes in local susceptibility may take place already during the incident pulse. As illustrated in Fig. 1, this requires one to solve the equations governing the dynamics of the properties of the material, in the presence of a driving optical field. The intensity of the optical field is obtained by solving Maxwell’s equations, with the material properties as an input.

 

Fig. 1. Diagrammatic view of the model. The iteration starts at Maxwell’s equations, which solve the propagation of the laser pulse in an initially unexcited material. Once every optical half-cycle, the optical intensity inside the material can be derived using the electric field obtained with Maxwell’s equations. This optical intensity is used as an input to the material equations that describe the response of the material to the intensity. By solving the material equations, parameters of the excited material, such as electron density and electron temperature, are obtained. Subsequently, a new susceptibility can be deduced from the parameters of the excited material. This new susceptibility and thus the dielectric function are further inserted into Maxwell’s equations to complete a single iteration. The iterations are conducted once every optical half-cycle until the end of the pulse.

Download Full Size | PDF

If we assume that the susceptibility changes slowly with respect to the oscillation period of the light, we can neglect time derivatives of the susceptibility and thus write the time derivative of the displacement field as $\partial \mathbf{D}/\partial t=(1+\chi )\partial \mathbf{E}/\partial t$. This reduces Maxwell’s equations to

To determine the susceptibility, we need to know certain position- and time-dependent properties of the material. For instance, we require the carrier density $N(\mathbf{r},t)$, which we can obtain by integrating the diffusion equation

The carrier density distribution obtained from the diffusion equation is subsequently used to calculate the position- and time-dependent susceptibility (and thus the dielectric function):

In the next section, we will follow the outline indicated above to derive a model to compute the self-reflectivity of silicon and gold samples subjected to intense femtosecond laser illumination. We will explicitly discuss the physical processes that should be taken into account and give details of the numerical implementation of the model.

3. THE MODEL

A. Laser–Matter Interaction for Silicon

The initial process during the absorption of a femtosecond laser pulse by a semiconductor is the excitation of electrons from the valence band to the conduction band. Depending on the bandgap of the material and the incident photon energy, the excitation can be either one-photon absorption (OPA) or multiphoton absorption, or both. Subsequently, electrons already excited to the conduction band can gain energy in the laser field via free-carrier absorption. If the excess energy of a conduction electron is sufficiently high, it may excite another electron in the valence band to the conduction band, by a process known as impact ionization. When this process occurs multiple times, it is referred to as avalanche or cascade ionization. For silicon at 800 nm excitation wavelength, free carriers are created by OPA, two-photon absorption (TPA) [7,16], and impact ionization [17]. Generally, impact ionization is a very complicated process that involves the energy distribution of the electrons [18]. Here we use a simplified and convenient expression for the impact ionization term [19,20]. We rewrite Eq. (3) as

Although the carrier temperature does not explicitly appear in the above equations for the dielectric function, we nevertheless have to calculate it. This is because the carrier diffusivity $D_0$ and the optical effective mass $m^{*}$ depend on the carrier temperature. In Ref. [21], a complete two-temperature model is given. In our case, however, we can safely neglect the change of lattice temperature, since our pulse width is much shorter than the electron–phonon coupling time, which is of the order of a picosecond [16]. Thus, in order to model the change of the carrier temperature $T_c$ during the pulse, we use a simplified carrier heat equation [21,7,16]:

One should distinguish ${\alpha }_{\mathrm{ex}}$ in Eq. (13) from ${\alpha }_0$ in Eq. (5) since ${\alpha }_{\mathrm{ex}}$ includes both the contribution of intraband and interband transition. The total energy density $U_c$ of the excited carriers is given by

Rigorously speaking, the carrier–carrier collision time ${\tau }_d$ also dynamically changes due to the change of carrier temperature during the pulse. However, as we will discuss in Section 4, in silicon the effect of the changes in the carrier–carrier collision time can be neglected with respect to those of the carrier density.

B. Laser–Matter Interaction for Gold

In the case of femtosecond laser ablation of gold, the leading optical absorption process is different. As there is already a high density of free electrons present in the unexcited material, the number of free electrons is not increased by orders of magnitude during the laser pulse as is the case in silicon. There is an increase of free electron density in gold due to the transition from the $d$- to the $s$-band, as discussed by many groups [24,25]. However, its effect on reflectivity is minor compared with the effect of electron heating due to intraband absorption. This will be discussed in detail later. In analogy with Eq. (13), we use

The dielectric function of excited gold can then be written as

C. FDTD Model

To solve Eqs. (1) and (2), we use a FDTD method. As dictated by the Courant condition of the FDTD method (see, for instance, [28]) and the fine spatial grid used in our simulation, the $E$ and $H$ fields are updated hundreds of times per optical cycle. Once every optical half-cycle, we extract the electric-field amplitude $E_0(\mathbf{r},t)$ from that optical half-cycle. Details of how we extract the amplitude can be found in Appendix A. From the electric-field amplitude, we determine the intensity $I(\mathbf{r},t)$ using Eq. (6). We use this intensity to march Eqs. (5) and (13) (in the case of silicon) or Eq. (18) (in the case of gold) forward in time by a single step using an implicit Euler method. As implicit Euler methods are unconditionally stable, we can choose the time step in the Euler method as half an optical cycle. The other quantities, such as diffusivity, heat capacity, and heat conductivity, are subsequently calculated. After this, we use the carrier density $N(\mathbf{r},t)$ and/or the electron temperature $T_e(\mathbf{r},t)$ to obtain a new dielectric function ${ϵ}_{\mathrm{ex}}(\mathbf{r},t)$. This updated dielectric function is used in the next optical half-cycle of the FDTD simulation.

The self-reflectivity of a strongly focused laser pulse is in principle a three-dimensional (3D) problem. However, 3D-FDTD simulations are notoriously time and memory consuming. We therefore instead run 2D simulations for the TE and TM cases and use those to approximate the 3D reflectivity. This requires an extra step that we discuss at the end of this subsection. A schematic representation of our 2D-FDTD simulation box is shown in Fig. 2. The grid sizes for the simulations are chosen to be around 15 nm such that the errors in the absolute reflectivity are smaller than 0.01 (see Appendix B) and the device layer thicknesses can be written as integers times the grid size. The width of the simulation box is 2 μm, which is 2 times the size of the focused laser spot ($1\mathrm{\mu m}@1/{\mathrm{e}}^2$ of intensity). The incident pulse duration used in the simulation is 126 fs, which is the value measured experimentally using a single-shot autocorrelator. The source plane is located at 200 nm above the sample surface, while the near-field detector plane is located one cell above the sample surface. There are 30 grid points in each of the four perfectly matched layers (PMLs) to ensure negligible reflections at the boundaries. We tested the accuracy of our method by comparing to several benchmarks, as discussed in Appendix C. The FDTD simulation of dispersive and lossy material as the excited silicon and gold we study in this paper is implemented using the auxiliary differential equation (ADE) method [28]. Details of this implementation are given in Appendix D.

 

Fig. 2. Layout of the 2D-FDTD simulation box. The FDTD grid is excited by a soft source that is located 200 nm above the silicon–air interface. Scattered $\mathbf{E}$ and $\mathbf{H}$ near-field values are recorded at the detector plane to extract the reflectivity.

Download Full Size | PDF

To determine the field that will be scattered/reflected back by the sample, we run the simulation with and without the sample and take the difference in the electric field at the detector plane as the scattered near field. We therefore filter the high spatial frequency components from the scattered near fields, as described in Appendix E, in order to obtain the scattered far field. From the resulting fields, we calculate the reflected pulse fluence $F_{\mathrm{refl}}^{\mathrm{TE},\mathrm{TM}}(y)$ and the incident pulse fluence $F_{\mathrm{inc}}^{\mathrm{TE},\mathrm{TM}}(y)$. Here, the superscripts TE and TM denote the results for the TE and the TM case. To obtain the total reflected/incident pulse energy, we add the TE and TM contributions. To approximate the 3D results, we treat the $y$ coordinate in our simulation as the radial coordinate in a polar coordinate system and integrate the reflected/incident fluence over the area of the incident focal spot:

This equation yields the value for $R$ that we will compare with experimental measurements in the next section.

4. RESULTS AND COMPARISON WITH EXPERIMENTS

In the case of silicon we have carried out simulations on thin-film SOI samples as well as bulk silicon. Experimental results for these samples can be found in Ref. [11]. For the simulations, a number of input parameters need to be specified. The values we used are listed in Table 1. In addition to these parameters, values for the Si and ${\mathrm{SiO}}_2$ layer thicknesses are required for the SOI samples. Table 2 lists the parameters used in our simulations. As can be seen, we inserted values for the layer thicknesses slightly deviating from their measured values. These adjusted values were chosen in order to yield the correct self-reflectivity at vanishing fluence. It should be noted that the reflectivity in this regime depends only on the thicknesses of the layers and the refractive indices of the unperturbed media.

Table 1. Material Parameters Used in the Simulation as Obtained from Literature^a

View Table | View all tables in this article

Table 2. Sample Parameters^a

View Table | View all tables in this article

Figure 3 shows the self-reflectivity calculated using the model and the experimental data of the bulk silicon and the SOI samples. As is shown for the ${\mathrm{SOI}}_1$ sample, the experimental data lie between the calculated self-reflectivity for the TM and the TE modes and agree well with the calculation considering both modes. For clarity, we do not show the TM and TE modes separately for the ${\mathrm{SOI}}_2$ and bulk samples. The experimental data scatters at low fluences, due to electronic noise and quantization errors in the data acquisition. This is because we designed our setup to work at fluences at and above the ablation threshold. A more detailed description of the setup can be found in Ref. [11]. At fluences above $0.05\mathrm{J}/{\mathrm{cm}}^2$ the scatter of the experimental data quickly reduces. With small adjustments in the values of $\theta$ and ${\tau }_d$, the agreement with the experimental data is even better [11]. We find that the reflectivities for the bulk and the ${\mathrm{SOI}}_2$ samples are very similar, whereas the reflectivity of the ${\mathrm{SOI}}_1$ sample deviates strongly from the other two samples. This is because the device layer of the ${\mathrm{SOI}}_1$ sample is thick enough to allow for constructive interference in the layer at 800 nm, impossible for the bulk sample and the 100 nm thin device layer of the ${\mathrm{SOI}}_2$ sample. As the incident fluence increases, the reflectivity drops to a minimum and then increases again. This behavior suggests a typical free-carrier (Drude) response. If one applies Eq. (7) with the Fresnel equation to a simple homogeneous plasma, one would find that after the critical density (where the plasma frequency equals the incident light frequency), as the carrier density further increases, the reflectivity drops to a minimum and then increases. The dashed lines in Fig. 3 show the results of 1D-FDTD calculations using the same parameters as the 2D calculations. The disagreement with the experimental data of the 1D calculations directly shows that the substantial range of incident angles must be taken into account when a high-NA objective is used, as is the case in nanoablation experiments. Note that, in the above calculations, we used a constant carrier collision time ${\tau }_d=1\mathrm{fs}$. As mentioned earlier, the carrier collision time is actually both density and temperature dependent. A standard result from Landau-liquid theory is that the mean collision time of a thermalized electron bath of temperature $T_e$ is given by [30]

 

Fig. 3. 2D-FDTD calculations and experimental measurements of self-reflectivity for bulk silicon and SOI samples. Open and closed symbols indicate results of two independent experimental runs for bulk silicon (circles), ${\mathrm{SOI}}_1$ (squares), and ${\mathrm{SOI}}_2$ (diamonds). The solid lines show the reflectivity calculated by 2D-FDTD simulations with the TM and TE modes combined for bulk (green), ${\mathrm{SOI}}_1$ (blue), and ${\mathrm{SOI}}_2$ (red). The dashed lines show the results obtained from a 1D-FDTD simulation. The blue dashed–dotted and dotted lines show the reflectivity calculated for ${\mathrm{SOI}}_1$ using either the TE or the TM mode, respectively. The left inset shows the reflectivity calculated using a dynamically changing collision rate with a $\sqrt{T}$ dependence, and the right inset shows the reflectivity calculated using a collision rate with a $T^2$ dependence.

Download Full Size | PDF

We now turn to the importance of impact ionization in the carrier creation process. In Fig. 4, we show the calculated results with the impact ionization coefficient $\theta$ set to zero. The disagreement with the experimental data in Fig. 4 and the good agreement with the experimental data in Fig. 3 demonstrate directly that impact ionization plays a significant role in the development of the dense electron-hole plasma in silicon induced by a single femtosecond laser pulse. It should be pointed out that in earlier work [7] the role of impact ionization was ignored, which resulted in an underestimation of the carrier density and a fitted (static) optical effective mass of $m^{*}=0.18m_e$. However, as shown by Riffe’s theoretical work, at a carrier temperature of 3000 K the optical effective mass of silicon already exceeds $0.24m_e$ [23]. Considering that carrier temperatures above ${10}^4\mathrm{K}$ are reached in the experiments [6,7], the value $m^{*}=0.18m_e$ reported in Ref. [7] is far too low. To see whether these conditions also occur in our model, we inspect the carrier density and the carrier temperature as calculated in our model. In Fig. 5, we plot the results of those quantities. Specifically, we plot the values obtained at the surface of the sample, directly after the pulse. In Fig. 5(a) we find that the carrier density is clearly beyond ${10}^{22}{\mathrm{cm}}^{-3}$ at the excitation level relevant for this work.

 

Fig. 4. FDTD calculations and experimental measurements of self-reflectivity on the ${\mathrm{SOI}}_1$, ${\mathrm{SOI}}_2$, and bulk samples. For the calculations, an impact ionization coefficient $\theta =0{\mathrm{cm}}^2/\mathrm{J}$ is used. The lines and symbols have the same meaning as in Fig. 3.

Download Full Size | PDF

 

Fig. 5. Calculated (a) carrier density, (b) carrier temperature, and (c) optical effective mass on the surface of the sample after the end of the pulse. In each plot, the green line shows the data for the bulk silicon sample, the blue line for the ${\mathrm{SOI}}_1$ sample, and the red line for the ${\mathrm{SOI}}_2$ sample.

Download Full Size | PDF

In Fig. 5(b), we see that the calculated carrier temperature is indeed larger than ${10}^4\mathrm{K}$ for all but the lowest fluences. Finally, in Fig. 5(c), we plot the optical effective mass obtained from our model. It is also clearly beyond its unperturbed value ($0.15m_e$). This is caused by the high carrier temperature reached in the laser-induced plasma.

To test the versatility of our method, we also carried out experiments and simulations on a 400 nm thick gold film grown on a glass substrate by vapor deposition. As the available electron heat-capacity data of gold obtained by ab initio calculations is limited to temperatures below 4.3 eV [26], the simulation results in Fig. 6 have a reduced accuracy for fluences $>1\mathrm{J}/{\mathrm{cm}}^2$. As mentioned in Section 3.B, in contrast to the case of semiconductors, metals already have such a high carrier concentration that the heating of those carriers has the dominant effect on the optical properties of the samples. In Fig. 6(a) we can clearly see this simplicity in the experimental results: the self-reflectivity starts out high for low fluences and decreases slowly with increasing fluence. The lines in the plot are the results from our model, where no free parameters where used. All fixed parameters, as listed in Table 3, were taken from literature. We can see that the model predicts the self-reflectivity very well from the unperturbed reflectivity without free parameters. The calculation also yields the electron temperature and the Drude damping time, as shown in Figs. 6(b) and 6(c), respectively. We find also, in the case of gold, final temperatures beyond ${10}^4\mathrm{K}$ ($\sim 1\mathrm{eV}$) for all but the lowest fluences. In order to separate out the effect of the increased conduction-electron density on the reflectivity from the decreased electron collision time, we further calculated the reflectivity using the constant electron density $N=5.9\times {10}^{22}/{\mathrm{cm}}^3$ of gold at room temperature [24]. The result is plotted as the dashed line in Fig. 6(a). We see that the increasing conduction-electron density indeed has a very small effect compared to the decreasing collision time, but certainly neglecting the increased conduction-electron density does lead to a slight overestimation of the change in reflectivity, especially in the high fluence regime.

 

Fig. 6. (a) FDTD calculation and experimental measurements of self-reflectivity on a 400 nm gold film on glass. Open and closed squares are the experimental data from two independent runs. The solid line shows the calculated self-reflectivity taking both the TM and the TE mode into account. The dotted and dashed–dotted lines show the self-reflectivity taking only the TE mode or only the TM mode into account, respectively. The dashed line shows the self-reflectivity calculated with the conduction electron density at room temperature. (b) Calculated electron temperature on the surface of the sample. (c) The corresponding Drude damping time.

Download Full Size | PDF

Table 3. Physical Parameters Used in the Simulation for Gold^a

View Table | View all tables in this article

5. SUMMARY AND CONCLUSIONS

We presented a model describing the propagation and absorption of a strongly focused femtosecond laser pulse used for single-shot laser ablation in semiconductor and metal samples, based on an extended 2D-FDTD method. The model is compared with self-reflectivity measurements of strongly focused femtosecond laser pulses used for single-shot femtosecond laser ablation on two SOI samples, bulk silicon, and gold. We obtain good agreement between simulation and experiments, using the unperturbed reflectivity to fix material and sample specific constants; the self-reflectivity at high fluences follows without open parameters. This confirms the accuracy and robustness of the model. The model shows the important role of impact ionization for the carrier generation in silicon induced by a femtosecond laser pulse of 800 nm. The model also indicates a marked increase of the optical effective mass in silicon due to the elevation of the carrier temperature during the pulse. Furthermore, the simulations demonstrate that the dominant factors determining the self-reflectivity are the carrier density in the case of silicon and the electron–electron collision time in the case of gold.

These results prove that FDTD simulations incorporating production and heating of free carriers and the free-carrier Drude response faithfully describe the behavior of the self-reflectivity of strongly focused femtosecond laser beams under the ablation conditions. As the simulation accurately predicts the self-reflectivity, it also gives a detailed understanding of the energy deposition of femtosecond laser pulses in metals and semiconductors. The method can easily be further extended to study the effects of heat and charge transport after the excitation, which could be compared to transient reflectivity measurements. In conclusion, this extended FDTD method is a flexible and indispensable tool in the study of femtosecond laser nano-structuring of materials.

APPENDIX A: EXTRACTING THE COMPLEX AMPLITUDE FROM FDTD SIMULATION

The FDTD method calculates real-value, time-varying electric and magnetic fields. However, some relevant physical quantities, such as the intensity [Eq. (6)], are more conveniently expressed in the amplitude of the oscillation. We extract the amplitudes of the fields from the simulation as follows. If we assume the amplitude is slowly varying with respect to the optical cycle, we can write the electric field on time $t$ as

(A1)$$E(t)=E_0\cos ({\omega }_{0}t+\varphi ),$$

where ${\omega }_0$ is the frequency of the light, and $E_0$ and $\varphi$ are the amplitude and phase, respectively. If we integrate $E^2(t)$ over half an optical cycle $T/2$, we find

(A2)$$\begin{gathered} {\int }_0^{T/2}{E}^2(t)\mathrm{d}t={\int }_0^{T/2}{E}_0^2{\cos }^2({\omega }_{0}t+\varphi )\mathrm{d}t \\ =\frac{1}{4}T{E}_0^2, \end{gathered}$$

and thus

(A3)$$E_0=\sqrt{\frac{4{\int }_0^{T/2}{E}^2(t)\mathrm{d}t}{T}}.$$

In the FDTD simulation, the integration is approximated as a summation:

(A4)$${\int }_0^{T/2}{E}^2(t)\mathrm{d}t\cong \sum _{n=n_1}^{n_1+m}{E}^2(n\mathrm{\Delta }t)·\mathrm{\Delta }t,$$

where $n_1$ is the starting time step of the summation and $m$ is the total time steps contained in half an optical cycle.

To extract the phase $\varphi$ we consider the integral

(A5)$${\int }_0^{T/2}{E}_0\cos ({\omega }_{0}t+\varphi )e^{-i{\omega }_{0}t}\mathrm{d}t=\frac{E_0\pi }{2{\omega }_0}(\cos \varphi +i\sin \varphi ).$$

Thus, the phase of the electric field is the phase of above integral.

APPENDIX B: FDTD ACCURACY

The FDTD method has intrinsic second-order accuracy because it uses central difference for both the time and space derivative. A second source of error of a FDTD code is the small residual reflection at the PML boundary. In this appendix, we analyze the numerical error due to the FDTD algorithm. Based on this result, we deduce the optimal grid size and PML thickness for the simulations performed in this paper. To analyze the error, we calculate the reflectivity of a plane wave under normal incidence using 1D-FDTD method and compare the result to the exact Fresnel result [33]:

(B1)$$R={\left(\frac{n-1}{n+1}\right)}^2,$$

which for silicon at 800 nm incident wavelength yields a reflectivity of 0.3287. The 1D-FDTD simulation space we use consists 50 nm of vacuum and 1 μm silicon, sandwiched between two PML layers. In the simulations we record the time-varying $E$ and $H$ field at the detector that is located one grid space above the silicon/vacuum interface. From the recorded time-varying fields, one can calculate the reflected Poynting vector and thus deduce the reflectivity $R^{\prime }$.

Figure 7 shows the resulting absolute error $\mathrm{err}=R^{\prime }-R$ as a function of grid spacing $\mathrm{\Delta }x$ and PML layer thicknesses $n_{\mathrm{pml}}$ measured in the number of cells in the PML layer. From the figure, we conclude that in order to keep the error in the absolute reflectivity smaller than 0.01, one needs a grid spacing $\le 15\mathrm{nm}$ and a PML layer not less than 20 points.

 

Fig. 7. Error of the FDTD method. Results with a range of PML thicknesses are presented.

Download Full Size | PDF

APPENDIX C: BENCHMARK

As a test for the validity of our self-developed 2D-FDTD code, we calculate the reflectivity of both unperturbed SOI with a 230 nm thick device layer and unperturbed bulk silicon with our FDTD code and compare the results with a commercial FDTD solver [34]. The oxide layer of the SOI wafer is 1 μm. The grid size in the simulation is set as 13 nm. We extract the reflectivity from the field values calculated by the simulations. The results are summarized in Table 4. We further tested the code by calculating the reflectivity of a micrometer-sized scatterer (dispersive and lossy) embodied in the device layer of the SOI wafer in order to test our implementation of the ADE method. As can been seen, the results obtained by the FDTD Solutions and our FDTD code are very close to one another. The small differences are most likely due to differences in PML layer thicknesses and/or the grid spacing.

Table 4. Comparison between the Results of Our 2D-FDTD Code with the Results from the Commercial Software FDTD Solutions^a

View Table | View all tables in this article

APPENDIX D: DISPERSIVE AND LOSSY MEDIA

In two dimensions, Maxwell’s equations reduce to independent equations for the TM and TE modes [33]. For the TM mode, these become [28]

(D1)$$\frac{\partial {\tilde{D}}_z}{\partial t}=\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}(\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}),$$

(D2)$${\tilde{D}}_z(\omega )={ϵ}_r^{*}(\omega ){\tilde{E}}_z(\omega ),$$

(D3)$$\frac{\partial H_x}{\partial t}=-\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}\frac{\partial {\tilde{E}}_z}{\partial y},$$

(D4)$$\frac{\partial H_y}{\partial t}=\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}\frac{\partial {\tilde{E}}_z}{\partial x}.$$

For the TE mode they become

(D5)$$\frac{\partial {\tilde{D}}_x}{\partial t}=\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}\frac{\partial H_z}{\partial y},$$

(D6)$$\frac{\partial {\tilde{D}}_y}{\partial t}=-\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}\frac{\partial H_z}{\partial x},$$

(D7)$${\tilde{D}}_x(\omega )={ϵ}_r^{*}(\omega ){\tilde{E}}_x(\omega ),$$

(D8)$${\tilde{D}}_y(\omega )={ϵ}_r^{*}(\omega ){\tilde{E}}_y(\omega ),$$

(D9)$$\frac{\partial H_z}{\partial t}=-\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}(\frac{\partial {\tilde{E}}_y}{\partial x}-\frac{\partial {\tilde{E}}_x}{\partial y}).$$

In the above equations, we normalize the electric and displacement fields as

(D10)$$\tilde{\mathbf{E}}=\sqrt{\frac{{ϵ}_0}{{\mu }_0}}\mathbf{E},$$

(D11)$$\tilde{\mathbf{D}}=\sqrt{\frac{{ϵ}_0}{{\mu }_0}}\mathbf{D},$$

to make the electric and magnetic fields the same order of magnitude, and thus increase the numerical precision.

Note that Eqs. (D2), (D7), and (D8) are expressed in the frequency domain, whereas the FDTD method works in the time domain. To bring the frequency domain equations into the time domain, we need to assume they are of a known analytical form. Here, we use a Drude model:

(D12)$$\begin{gathered} {ϵ}_r^{*}(\omega )={ϵ}_{\mathrm{ex}} \\ ={ϵ}_{\mathrm{Si}}-{\left(\frac{{\omega }_p}{\omega }\right)}^2\frac{1}{1+i\frac{1}{\omega \tau }}. \end{gathered}$$

Using partial fraction expansion and switching the imaginary unit from $i$ to $j$ ($j=-i$, as is conventional in engineering), Eq. (D12) can be written as

(D13)$$\begin{gathered} {ϵ}_{\mathrm{ex}}={ϵ}_{\mathrm{Si}}+\frac{{\omega }_p^2\tau }{j\omega }-\frac{{\omega }_p^2{\tau }^2}{1+j\omega \tau } \\ ={ϵ}_{\mathrm{Si}}+\frac{{\sigma }_{\text{Drude}}}{j\omega {ϵ}_0}+\frac{\chi }{1+j\omega \tau }, \end{gathered}$$

where we introduced the parameters

(D14)$${\sigma }_{\text{Drude}}={ϵ}_0{\omega }_p^2\tau ,$$

(D15)$$\chi =-{\omega }_p^2{\tau }^2,$$

where ${\sigma }_{\text{Drude}}$ is the conductivity of the plasma. The $\chi$ term causes additional dispersion. This is referred to as the Debye formulation [28] of the Drude model. Inserting this dielectric function, the electric displacement reads

(D16)$$\begin{gathered} D(\omega )={ϵ}_{\mathrm{ex}}(\omega )E(\omega ) \\ ={ϵ}_{\mathrm{Si}}E(\omega )+\frac{{\sigma }_{\text{Drude}}}{j\omega {ϵ}_0}E(\omega )+\frac{\chi }{1+j\omega \tau }E(\omega ) \\ =D^{\prime }(\omega )+S(\omega ), \end{gathered}$$

where the first two terms of the right-hand side of the equation are summarized as $D^{\prime }(\omega )$ and the last term is written as $S(\omega )$.

As the FDTD operates in the time domain, Eq. (D16) must be transformed into the time domain. We first transform $D^{\prime }(\omega )$ into the time domain. Recall that $1/j\omega$ in the frequency domain corresponds to integration in the time domain, so $D^{\prime }(\omega )$ becomes

(D17)$$D^{\prime }(t)={ϵ}_{\mathrm{Si}}E(t)+\frac{{\sigma }_{\text{Drude}}}{{ϵ}_0}{\int }_0^{t}E(t^{\prime })\mathrm{d}{t}^{\prime }.$$

This integral is approximated as a summation over the time steps $\mathrm{\Delta }t$:

(D18)$$D^{\prime n}={ϵ}_{Si}{E}^n+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}\sum _{i=0}^n{E}^i,$$

where $n$ indicates the time step at $t=n\mathrm{\Delta }t$. In the FDTD algorithm, we use this equation to determine the current $E^n$ from the current $D^{\prime n}$ and the previous values of $E$. We do this by first separating the $E^n$ term from the rest of the summation,

(D19)$$D^{\prime n}={ϵ}_{Si}{E}^n+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}{E}^n+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}\sum _{i=0}^{n-1}{E}^i,$$

and solving that equation for $E^n$:

(D20)$$E^n=\frac{D^{\prime n}-\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}\sum _{i=0}^{n-1}{E}^i}{{ϵ}_{Si}+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}}.$$

We now treat the $S(\omega )$ term in Eq. (D16) in a similar manner. We convert the term

(D21)$$S(\omega )=\frac{\chi }{1+j\omega \tau }E(\omega )$$

to the time domain to find

(D22)$$S(t)=\frac{\chi }{\tau }{\int }_0^t{e}^{-\frac{t-t^{\prime }}{\tau }}E(t^{\prime })\mathrm{d}{t}^{\prime }.$$

We approximate the integral as a summation:

(D23)$$\begin{gathered} S^n=\chi \frac{\mathrm{\Delta }t}{\tau }\sum _{i=0}^n{e}^{-\frac{\mathrm{\Delta }t(n-i)}{\tau }}{E}^i \\ =\chi \frac{\mathrm{\Delta }t}{\tau }[E^n+\sum _{i=0}^{n-1}{e}^{-\frac{\mathrm{\Delta }t(n-i)}{\tau }}{E}^i]. \end{gathered}$$

We can add Eq. (D23) to Eq. (D19) to find the electric displacement for all three terms,

(D24)$$\begin{gathered} D^n={ϵ}_{\mathrm{Si}}{E}^n+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}{E}^n+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}\sum _{i=0}^{n-1}{E}^i \\ +\chi \frac{\mathrm{\Delta }t}{\tau }[E^n+\sum _{i=0}^{n-1}{e}^{-\frac{\mathrm{\Delta }t(n-i)}{\tau }}{E}^i], \end{gathered}$$

which, when we solve for $E^n$, yields

(D25)$$E^n=\frac{D^n-\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}\sum _{i=0}^{n-1}{E}^i-\chi \frac{\mathrm{\Delta }t}{\tau }\sum _{i=0}^{n-1}{e}^{-\frac{\mathrm{\Delta }t(n-i)}{\tau }}{E}^i}{{ϵ}_{Si}+\frac{{\sigma }_{\text{Drude}}\mathrm{\Delta }t}{{ϵ}_0}+\chi \frac{\mathrm{\Delta }t}{\tau }}.$$

As in the case of pure damping, we calculate $E^n$ (the current value of $E$) from the current value of $D$ and the summation of all previous values of $E$.

We can easily add the contributions of two-photon absorption (TPA) and the optical Kerr effect into Eq. (D25) by introducing an extra conductivity term due to the TPA and an extra change in the real part of the dielectric constant due to the optical Kerr effect. Together with the OPA, Eq. (D25) can finally be written as

(D26)$$E^n=\frac{D^n-\frac{\sigma \mathrm{\Delta }t}{{ϵ}_0}\sum _{i=0}^{n-1}{E}^i-\chi \frac{\mathrm{\Delta }t}{\tau }\sum _{i=0}^{n-1}{e}^{-\frac{\mathrm{\Delta }t(n-i)}{\tau }}{E}^i}{R\mathrm{e}[{ϵ}_{\mathrm{Si}}]+\mathrm{\Delta }ϵ+\frac{\sigma \mathrm{\Delta }t}{{ϵ}_0}+\chi \frac{\mathrm{\Delta }t}{\tau }},$$

where $\sigma ={\sigma }_{\mathrm{OPA}}+{\sigma }_{\mathrm{TPA}}+{\sigma }_{\text{Drude}}$ and $\mathrm{\Delta }ϵ$ is the change in the real part of the dielectric constant due to the optical Kerr effect. The term ${\sigma }_{\mathrm{OPA}}$ can be linked to the imaginary part of the dielectric constant of unexcited silicon. The terms $\mathrm{\Delta }ϵ$ and ${\sigma }_{\mathrm{TPA}}$ are linked to the real and imaginary parts of ${\chi }^{(3)}$, respectively.

APPENDIX E: SPATIAL FILTERING

In the simulation, interaction between the light and the localized laser-induced plasma gives rise to field components with high spatial frequencies. When we decompose the wave vector into the transverse and the longitudinal component, they are related by

(E1)$$k_x=\sqrt{k^2-k_y^2},$$

where $k$ is the wavenumber in air and $k_y$ is the transverse wavenumber. We observe from the above equation that for $k_y^2>k^2$ the longitudinal wavenumber $k_x$ is purely imaginary; therefore, the wave associated with it is evanescent and will not propagate into the far field. Prior to calculating the reflectivity, these components should therefore be filtered out. The transverse wavenumber can be conveniently written as $k_y=k\sin \theta$, where $\theta$ is the angle of reflection. As the numerical aperture of the objective used in our experiments is $\mathrm{NA}=\sin \theta =0.8$, this means that scattered waves associated with a transverse wavenumber $k_y>0.8k$ will not be collected by the objective. This means that what is measured in Figs. 3 and 6 is in fact the scattered field associated with a transverse wavenumber $k_y<0.8k$. So in the FDTD simulation, spatial-frequency components with $k_y>0.8k$ should be filtered out. This is accomplished by the spectral decomposition fields [35],

(E2)$$\widehat{A}(k_y)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }A(y)e^{-i{k}_{y}y}\mathrm{d}y,$$

and by transforming this composition back using a truncated inverse Fourier transform:

(E3)$$A^{\prime }(y)=\frac{1}{\sqrt{2\pi }}{\int }_{-0.8k}^{0.8k}\widehat{A}(k_y)e^{i{k}_{y}y}\mathrm{d}{k}_y.$$

ACKNOWLEDGMENTS

We thank Cees de Kok and Paul Jurrius for discussions and technical assistance. We thank Sandy Pratama for proofreading the manuscript. HZ acknowledges financial support from the China Scholarship Council. DMK acknowledges support from NSF grant DMR 1206979.

REFERENCES

1. E. G. Gamaly, S. Juodkazis, K. Nishimura, H. Misawa, B. Luther-Davies, L. Hallo, Ph. Nicolai, and V. T. Tikhonchuk, “Laser-matter interaction in the bulk of a transparent solid: Confined microexplosion and void formation,” Phys. Rev. B 73, 214101 (2006).

2. Y. Liao, Y. Shen, L. Qiao, D. Chen, Y. Cheng, K. Sugioka, and K. Midorikawa, “Femtosecond laser nanostructuring in porous glass,” Opt. Lett. 38, 187–189 (2013). [CrossRef]  

3. H. Zhang, D. van Oosten, D. M. Krol, and J. I. Dijkhuis, “Saturation effects in femtosecond laser ablation of silicon-on-insulator,” Appl. Phys. Lett. 99, 231108 (2011). [CrossRef]  

4. M. Li, K. Mori, M. Ishizuka, X. Liu, Y. Sugimoto, N. Ikeda, and K. Asakawa, “Photonic bandpass filter for 1550 nm fabricated by femtosecond direct laser ablation,” Appl. Phys. Lett. 83, 216–218 (2003). [CrossRef]  

5. B. Rethfeld, K. Sokolowski-Tinten, D. Von Der Linde, and S. I. Anisimov, “Timescales in the response of materials to femtosecond laser excitation,” Appl. Phys. A 79, 767–773 (2004). [CrossRef]  

6. D. Hulin, M. Combescot, J. Bok, A. Migus, J. Y. Vinet, and A. Antonetti, “Energy transfer during silicon irradiation by femtosecond laser pulse,” Phys. Rev. Lett. 52, 1998–2001 (1984). [CrossRef]  

7. K. Sokolowski-Tinten and D. von der Linde, “Generation of dense electron-hole plasmas in silicon,” Phys. Rev. B 61, 2643–2650 (2000).

8. K. Sokolowski-Tinten, J. Bialkowski, and D. von der Linde, “Ultrafast laser-induced order-disorder transitions in semiconductors,” Phys. Rev. B 51, 14186–14198 (1995).

9. R. W. Boyd, Nonlinear Optics3rd ed. (Academic, 2008).

10. Th. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. 78, 3282–3285 (1997). [CrossRef]  

11. H. Zhang, D. M. Krol, J. I. Dijkhuis, and D. van Oosten, “Self-scattering effects in femtosecond laser nanoablation,” Opt. Lett. 38, 5032–5035 (2013). [CrossRef]  

12. L. Hallo, A. Bourgeade, V. T. Tikhonchuk, C. Mezel, and J. Breil, “Model and numerical simulations of the propagation and absorption of a short laser pulse in a transparent dielectric material: Blast-wave launch and cavity formation,” Phys. Rev. B 76, 024101 (2007).

13. C. Mézel, L. Hallo, A. Bourgeade, D. Hébert, V. T. Tikhonchuk, B. Chimier, B. Nkonga, G. Schurtz, and G. Travaillé, “Formation of nanocavities in dielectrics: a self-consistent modeling,” Phys. Plasmas 15, 093504 (2008). [CrossRef]  

14. I. B. Bogatyrev, D. Grojo, P. Delaporte, S. Leyder, M. Sentis, W. Marine, and T. E. Itina, “Non-linear absorption of 1.3-μm wavelength femtosecond laser pulses focused inside semiconductors: finite difference time domain-two temperature model combined computational study,” J. Appl. Phys. 110, 103106 (2011). [CrossRef]  

15. H. Schmitz and V. Mezentsev, “Full-vectorial modeling of femtosecond pulses for laser inscription of photonic structures,” J. Opt. Soc. Am. B 29, 1208–1217 (2012). [CrossRef]  

16. T. Y. Choi and C. P. Grigoropoulos, “Plasma and ablation dynamics in ultrafast laser processing of crystalline silicon,” J. Appl. Phys. 92, 4918–4925 (2002). [CrossRef]  

17. P. P. Pronko, P. A. VanRompay, C. Horvath, F. Loesel, T. Juhasz, X. Liu, and G. Mourou, “Avalanche ionization and dielectric breakdown in silicon with ultrafast laser pulses,” Phys. Rev. B 58, 2387–2390 (1998).

18. N. Medvedev and B. Rethfeld, “A comprehensive model for the ultrashort visible light irradiation of semiconductors,” J. Appl. Phys. 108, 103112 (2010). [CrossRef]  

19. E. Yablonovitch and N. Bloembergen, “Avalanche ionization and the limiting diameter of filaments induced by light pulses in transparent media,” Phys. Rev. Lett. 29, 907–910 (1972). [CrossRef]  

20. B. C. Stuart, M. D. Feit, S. Herman, A. M. Rubenchik, B. W. Shore, and M. D. Perry, “Nanosecond-to-femtosecond laser-induced breakdown in dielectrics,” Phys. Rev. B 53, 1749–1761 (1996).

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

22. N. W. Ashcroft and N. D. Mermin, Solid State Physics (Brooks/Cole, 1976).

23. D. M. Riffe, “Temperature dependence of silicon carrier effective masses with application to femtosecond reflectivity measurements,” J. Opt. Soc. Am. B 19, 1092–1100 (2002). [CrossRef]  

24. C. Fourment, F. Deneuville, D. Descamps, F. Dorchies, S. Petit, O. Peyrusse, B. Holst, and V. Recoules, “Experimental determination of temperature-dependent electron-electron collision frequency in isochorically heated warm dense gold,” Phys. Rev. B 89, 161110 (2014).

25. E. Bevillon, J. P. Colombier, V. Recoules, and R. Stoian, “Free-electron properties of metals under ultrafast laser-induced electron-phonon nonequilibrium: a first-principles study,” Phys. Rev. B 89, 115117 (2014).

26. Z. Lin, L. V. Zhigilei, and V. Celli, “Electron-phonon coupling and electron heat capacity of metals under conditions of strong electron-phonon nonequilibrium,” Phys. Rev. B 77, 075133 (2008). Updated calculations at http://www.faculty.virginia.edu/CompMat/electron-phonon-coupling/.

27. K. Eidmann, J. Meyer-ter-Vehn, T. Schlegel, and S. Huller, “Hydrodynamic simulation of subpicosecond laser interaction with solid-density matter,” Phys. Rev. E 62, 1202–1214 (2000). [CrossRef]  

28. D. M. Sullivan, “Electromagnetic simulation using the FDTD method,” IEEE, 2000.

29. A. D. Bristow, N. Rotenberg, and H. M. Van Driel, “Two-photon absorption and Kerr coefficients of silicon for 850–2200 nm,” Appl. Phys. Lett. 90, 191104 (2007). [CrossRef]  

30. B. Y. Mueller and B. Rethfeld, “Relaxation dynamics in laser-excited metals under nonequilibrium conditions,” Phys. Rev. B 87, 035139 (2013).

31. S. J. Youn, T. H. Rho, B. I. Min, and K. S. Kim, “Extended Drude model analysis of noble metals,” Phys. Status Solidi B 244, 1354–1362 (2007). [CrossRef]  

32. K. Vestentoft and P. Balling, “Formation of an extended nanostructured metal surface by ultra-short laser pulses: single-pulse ablation in the high-fluence limit,” Appl. Phys. A 84, 207–213 (2006). [CrossRef]  

33. M. Born and E. Wolf, Principles of Optics (Cambridge University, 1999).

34. Lumerical Solutions, Inc.http://www.lumerical.com/tcad-products/fdtd/.

35. L. Novotny and B. Hecht, Nano Optics (Cambridge University, 2006).