## Abstract

We present a numerical model that describes the propagation of a single femtosecond laser pulse in a medium of which the optical properties dynamically change within the duration of the pulse. We use a finite-difference time-domain method to solve the Maxwell’s equations coupled to equations describing the changes in the material properties. We use the model to simulate the self-reflectivity of strongly focused femtosecond laser pulses on silicon and gold under laser ablation conditions. We compare the simulations to experimental results and find excellent agreement.

© 2015 Optical Society of America

## 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}={\u03f5}_{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.

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):

where we have now explicitly written the susceptibility as a function of the carrier density. The dots in Eq. (4) indicate that the susceptibility is also a function of other properties of the medium, such as electron temperature and lattice temperature. If these are expected to vary on the timescale of the pulse, additional equations need to be included to describe their dynamics. For instance, the susceptibility in gold also depends strongly on the electron temperature. This will be discussed further in the following sections.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

*et al.*, who measured the impact ionization rate (${\mathrm{s}}^{-1}$) for the dielectric breakdown in silicon with 786 nm fs laser pulses [17]. Due to the creation of a high density of carriers, the optical response of silicon under ablation conditions is dominated by the free-carrier response, which can be calculated using the Drude model [7,16]. In our model, we also take into account the changes in the dielectric constant due to the optical Kerr effect and TPA. This is accomplished by an extra ${\u03f5}_{\mathrm{NL}}$ term in the total dielectric function ${\u03f5}_{\mathrm{ex}}$, as shown in Eqs. (7) and (9). Thus the dielectric function of strongly excited silicon ${\u03f5}_{\mathrm{ex}}$ can be written 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

with ${C}_{c}$ the carrier heat capacity, ${E}_{g}$ the bandgap energy, and ${\kappa}_{c}$ the thermal conductivity of the carriers. The carrier density and temperature we find through the simulation in silicon around the ablation threshold are $N$ on the order of ${10}^{22}/{\mathrm{cm}}^{2}$ and ${T}_{c}\sim {10}^{4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}$ (Fig. 5). These numbers are in general agreement with the estimation made by Sokolowski-Tinten and von der Linde from time-resolved pump–probe experiments [7]. The Fermi temperature of a free-electron gas of this density is ${T}_{F}\sim 2\times {10}^{4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}$. As this is on the same order as the temperature we obtained from the simulation, the carrier gas is neither in the degenerate nor in the nondegenerate limit. Thus we determine the carrier heat capacity by taking the temperature gradient of the total energy density of the free carriers, assuming a free electron density of states [22]. We subsequently determine the carrier heat conductivity using where ${k}_{B}$ is the Boltzmann constant, ${v}_{c}$ the thermal velocity of the free carrier, and ${\tau}_{d}$ the carrier–carrier collision time. The carrier diffusivity is then found using the Einstein relation [22]: with an effective mass ${m}^{*}$. In Ref. [7] the authors determine a static optical effective mass of the femtosecond-laser-induced plasma in silicon using time-resolved pump–probe experiments. However, as the optical effective mass is both temperature and density dependent, it will actually change during the pulse. When the electrons are heated up far beyond the conduction band edge, the curvature of the band decreases, giving rise to a larger effective mass. As experimental measurements of the optical effective mass are based on the Drude model, they depend only on the ratio $N/{m}^{*}$ [see Eq. (10)]. So to measure ${m}^{*}$ one needs to know the carrier density $N$, which in our case is*a priori*unknown. Therefore, we use the relation suggested by Riffe’s theoretical calculation, which takes the detailed band structure of silicon into account. His calculation shows that the optical effective mass can in the nondegenerated limit be approximated by [23] where the optical effective mass of unperturbed silicon ${m}_{0}^{*}=0.15{m}_{e}$ [7,23]. We extract the slope ${m}_{k}=3.1\times {10}^{-5}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{K}}^{-1}$ from Riffe’s calculations [23]. Although the above equation is valid for the nondegenerate limit, the elevation of carrier effective mass in the process is clear.

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

*et al.*, obtained with

*ab initio*calculations [26]. For very high electron temperatures, impact ionization or ballistic transport of electrons could actually change the electron density on the local scale. These effects are expected to have a small influence on the self-reflectivity, but could play an important role in the subsequent transfer of heat to the lattice. Very recently, Fourment

*et al.*[24] demonstrated the importance of $d$-band electrons for the electron–electron collision rate in noble metals. They experimentally determined the conduction electron density as well as the electron–electron collision rate in gold as a function of electron temperature in the range 0.6–5 eV. We incorporate the increase of conduction electrons into our FDTD model by fitting a simple spline function to their experimental results. The authors find that the electron–electron collision rate due to $d$-band electrons is given by ${\nu}_{\mathrm{ee}}^{d}=A(11-{Z}_{\mathrm{eff}})({Z}_{\mathrm{eff}}-1)$, where ${Z}_{\mathrm{eff}}$ is the effective number of conduction electrons per atom. They further find that for electron temperatures in the range of 0.6–5 eV, the collisions due to $d$-band electrons are dominant. We therefore use this expression and adjust the coefficient $A$ to best fit their experimental data. Finally, we obtain the following expression for the Drude damping term ${\tau}_{d}$ in gold: where ${\nu}_{\mathrm{ph}}=0.129\times {10}^{15}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ is the contribution of electron–phonon collisions in gold [27]. Because, in the temperature range and timescale considered, the electron–electron collision dominates [24], we fix the electron–phonon collision rate in our model.

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 ${\u03f5}_{\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.

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.

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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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]

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.18{m}_{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.24{m}_{e}$ [23]. Considering that carrier temperatures above ${10}^{4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}$ are reached in the experiments [6,7], the value ${m}^{*}=0.18{m}_{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}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-3}$ at the excitation level relevant for this work.

In Fig. 5(b), we see that the calculated carrier temperature is indeed larger than ${10}^{4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.15{m}_{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{K}$ ($\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.

## 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

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 findTo extract the phase $\varphi $ we consider the 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]:

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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and a PML layer not less than 20 points.

## 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.

## 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]

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:

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

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

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

to the time domain to findWe 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

## 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

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\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\theta $, where $\theta $ is the angle of reflection. As the numerical aperture of the objective used in our experiments is $\mathrm{NA}=\mathrm{sin}\text{\hspace{0.17em}}\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], and by transforming this composition back using a truncated inverse Fourier transform:## 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 Optics*3rd 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).