Using 3D Finite-Difference-Time-Domain simulations, we study the morphology of the laser-created damage inside fused silica. Among the competing effects limiting the intensity in the dielectric, we find the most important is the pulse defocusing by the plasma lens, partially balanced by the Kerr effect. Less important are collisional energy dissipation and laser depletion by multi-photon absorption. We also found that the profile of generated plasma is asymmetrical in the transverse cross-section, with the plasma extended along the direction perpendicular to the laser polarization.
© 2011 Optical Society of America
Recent advances in laser science have enabled the precise control of the modification of dielectric materials via intense femtosecond pulses. An ultrafast pulse can be focused within the bulk of a dielectric, which is transparent in the linear regime, but the intensity can be made high enough that nonlinear interactions such as multiphoton ionization cause permanent modification of the target region near the laser focus. The damage can be used to tailor the structural properties of dielectrics in three dimensions on a sub-wavelength scale. Depending on the energy deposited in the region, different types of damage have been observed, including changes to the index of refraction at intensities just above threshold , periodic nano-cracks inside the material above threshold , and the production of void-like structures at even higher intensities . The numerous potential applications of this technology include fabrication of three-dimensional photonic crystals, microfluidic devices, volume gratings, optical memory devices and optical waveguides [4, 5].
Precise control of the laser damage requires a detailed analysis of the dynamics of the laser-material interaction. There have been a number of numerical studies of the multi-photon ionization process and the resulting nonlinear electromagnetic wave propagation in the large-bandgap dielectrics such as water or fused silica  – . Of particular interest is identification of the particular mechanisms responsible for the morphology of the damaged region. Burakov et al.  use a two-dimensional nonlinear optical Schrodinger equation to analyze the ionization patterns produced by femtosecond and picosecond lasers. They attribute the elongated modified region to a balance between the nonlinear defocusing due to the plasma dispersion and the self-focusing of the Kerr effect. Rayner et al.  and Grojo et al.  propose that the Kerr effect is not important and the shape arises from the depletion of the incident laser beam. They argue that since there are orders of magnitude more ionizable molecules in solids than photons to excite these molecules, the depletion of the beam becomes an important factor in the nonlinear process. They hypothesize that the intensity will not exceed the threshold intensity because, once the intensity is high enough, ionization will occur and absorb photons, depleting the beam. In a focused laser, as the beam width narrows and the peak intensity increases, a larger portion of the intensity curve will be cut off. This causes the narrowing of the ionizing region. Rayner et al. model this effect by cutting the intensity at the threshold and assuming that the energy from the photons was used to ionize the material. Which of these models is correct, and in which parameter regimes, is still up for debate.
In this paper, we introduce a 3D dynamic model for simulating nonlinear ionization in dielectrics that enables a systematic analysis of the underlying processes. In particular, we investigate the time evolution of plasma creation to understand the microscopic mechanisms responsible for the morphology of the damage region in silica. The region is considered “damaged” in the text below wherever there is an appreciable level of ionization happened. We extended a standard Finite-Difference Time-Domain (FDTD) algorithm  to include multiphoton ionization (MPI), multiphoton absorption (MPA), the Kerr effect, and plasma dispersion, where the generation and optical response of the plasma is described via a plasma fluid model. The model introduced is intrinsically capable of modelling propagation of a tightly focused laser pulse in a dielectric medium with the possible over-dense plasma generation, i.e., it does not possess limitations of the uni-directional nonlinear models presented in the majority of the previous work. In this study we will concentrate on a relatively mild focusing with a numerical aperture (NA) equal to 0.65, as in Ref. .
In addition to the identification of the contribution of various mechanisms to the morphology of the damaged region, our 3D model enabled us to study effects of laser polarization. In particular, we have found that generated plasma pattern possesses an asymmetry in the transverse cross-section that can be explained by the refraction index mismatch between the plasma and the surrounding dielectric.
2.1. Basic equations
In our 3D numerical model, we solve Maxwell equations,20]. In Eq. (1), E⃗ and B⃗ are the components of electromagnetic field, D⃗ electric field displacement vector, H⃗ the magnetic field auxiliary vector, t is the time, and c is the speed of light. We extend the Yee discretization algorithm  by introducing the consitutive relations,
The free particle current J⃗p is calculated from the fluid equations for the electron component of the generated plasma:7] 7, 17] and ns is the saturation particle density. The quantity ns can be estimated as the electron density when every molecule of silica is singly-ionized: ns = ρmNA/Mμ, where ρm is the mass density of fused silica, Mμ its molar mass and NA Avogadro number. Substitution of the appropriate constants into this expression results in ns ≈ 10 ncr, where ncr is the critical particle density for electrons, i.e., the density for which the electron plasma frequency equals the laser frequency; ncr = 1.75 × 1021 cm−3 for the wavelength considered. In the actual calculations, we approximate the laser intensity by the square of instantaneous electric field: Eq. (5) assumes that the ionization occurs on a time scale slower than a wave period. The approximation given by Eq. (6) introduces an inaccuracy of order , where the angle brackets denote a cycle-average, which results in an under-estimation of the threshold intensities by approximately 35%. Performing an accurate cycle averaging in 3D is numerically unfeasible, given the large scale of our calculations as mentioned below, since each time point over a wave period at each grid cell location would need to be stored in memory.
The quantity Γ in Eq. (4) is the damping factor representing collisional energy dissipation. In the present work, we consider it to be approximately constant. This constant is estimated as follows. If we assume that the kinetic energy of an electron in the laser pulse is dissipated after a few collisions, then18]
E is the electric field in the laser pulse, ω0 is laser frequency and e, m are electron charge and mass, respectively. The laser intensity in the dielectric is limited by the ionization threshold, equal to (1 ÷ 1.5) × 1013 W/cm2 . Substituting this into (9), one obtains vq ≈ 3 × 10−3c, where c is the speed of light in vacuum. This velocity corresponds to an electron kinetic energy in eV range. Substituting this value into Eq. (8), one can obtain a pessimistic estimation for the damping factor
To solve Eq. (4), we make further approximations. To achieve hydrodynamic closure we assume a cold plasma model, for which p = 0. We also assume that the plasma fluid velocity u⃗ in the collision-dominated regime is small and thus (u⃗∇)u⃗ ≈ 0. Finally, we assume that the dominant mechanism causing changes to the particle density of electrons is ionization rather than the hydrodynamic advance of the plasma, i.e., ∇(nu⃗) ≪ ∂nMPI/∂t. Under these approximations, Eq. (4) simplifies to19]. One may calculate the current of free plasma electrons
The effective current J⃗MPA is obtained by setting the laser energy density depletion rate equal to the product (J⃗MPA · E⃗) . Assuming that J⃗MPA is parallel to E⃗ one obtains
A similar model was developed in Ref.  for a 2D plane geometry case.
2.2. Model benchmark
Now we demonstrate that the approximations made in Sec. 2.1 are reasonable and reproduce the basic physics of interaction.
Our first test is propagation of an electromagnetic field in an infinite cold plasma of predefined particle density with no ionization, multi-photon absorption or Kerr effect. The dispersion relation for a plane wave in a cold plasma isEq. (14).
A comparison, not presented here, of Eq. (16) with our actual FDTD simulation results revealed perfect matching. A similar test measuring the amount of reflection from underdense and overdense semi-infinite plasmas has also shown perfect agreement with the appropriate analytical formula.
Next we demonstrate the multi-photon ionization and absorption processes. For this test, we assume Jp ≡ 0, χk = 0, and propagate a plane electromagnetic wave through a 40-micron-thick slab of fused silica. Figure 1a shows the absolute value of Poynting vector of the pulse just before it enters the slab (at t = 0) and right after it exits the slab (t = 350 fs). The transmitted pulse is cut at intensity ∼ 1013 W/cm2. The dependence of the peak transmitted intensity on peak incident pulse intensity is given in Fig. 1b. It is seen that as soon as the incident laser exceeds the threshold intensity ∼ 1013 W/cm2, the pulses become depleted via ionization, and the transmitted intensity is capped by the threshold value. Due to the reflections from the vacuum-dielectric and dielectric-vacuum interfaces, the actual intensity threshold inside the dielectric should be slightly larger. The found value is consistent with the experimental measurements . An additional energy conservation test (not shown) demonstrated that the total energy in the system ℰem + ℰabs is conserved throughout the simulation with Jp ≡ 0, where ℰem is the total electromagnetic energy in the domain and ℰabs the absorbed energy.
Given the results of the demonstrated tests, we conclude that the code generates a reliable plasma response and correctly accounts for the multi-photon ionization and absorption processes.
2.3. Laser source excitation
In our 3D numerical model, the laser pulse is excited at the boundary of the box using the Total-Field-Scattered-Field approach [?]. This method generally requires an exact knowledge of the electromagnetic field being solution to the Maxwell equations at the boundary. To evaluate this field we use the model of the focused laser pulse, previously employed in , which allows us to use high NA-optics, including NA > 1. We give only a short description of this model here.
Let the focusing optics be represented by a paraboloidal mirror (Fig. 2), with a wide Gaussian beam incident onto the mirror. Right at the mirror surface, the field is evaluated using the proper boundary conditions at the ideally reflecting surface. In any point in the space, the field is evaluated using Stratton-Chu integrals :
In the present paper we use NA = 0.65, for a Gaussian beam incident of radius Rg, with Rg = 0.5Rm. The actual mirror radius used in the simulations was equal to 1 mm. The resulting focused pulse had full width at intensity half-maximum (in silica, in the absence of nonlinear effects) equal to 1 μm in the transverse direction and 7.5 μm in the longitudinal direction.
3. Morphology of the generated plasma
In our analysis of the generated plasma morphology we will consider the following laser parameters: laser wavelength λ = 800 nm, peak laser intensity at the focus (in the absence of nonlinear effects and dispersion) I = 1014 W/cm2 and the laser pulse length (defined as the full width at the half-maximum within a Gaussian envelope) equal to 50 fs. The corresponding laser pulse energy is 4 × 10−7 J. In all the examples shown below we assume the medium possesses a linear refraction index n0 ≈ 1.45. The Kerr susceptibility was assumed to be χ3 = 1.9 × 10−15 esu.
In our simulations we used a 50μm × 32μm × 32μm simulation domain, with 453 grid points per μm3. These simulations were performed on a 300-cpu computer cluster, with a single simulation requiring 650 Gb of RAM and approximately 24 hours of runtime. For all the results shown below, the laser pulse is propagating from left to right, with the laser axis passing trough the center of the simulation domain.
As the laser propagates in the dielectric, it ionizes it and creates a plasma. As was estimated above, the energy of quiver motion of electrons in the laser pulse is on the ∼eV scale. The corresponding sound velocity of plasma,
Figure 3a shows the profile of the generated plasma in the polarization plane passing through the laser focus for the fused silica occupying the entire computational domain. The laser pulse propagates in the +x̂ direction. The profile is characterized by a narrowing shape in the direction of laser propagation (cf. , ), with the maximum particle density being before the actual geometrical focus of the laser pulse.
To consider the effect of laser depletion, we have performed simulations with J⃗MPA = 0, i.e., with the effect of multi-photon absorption turned off, which is equivalent to assuming an infinite supply of photons. The resulting plasma profile is shown in Fig. 3b. It is seen that in the absence of multi-photon absorption, the maximum density of the generated plasma increases by ∼ 50%. However, the shape itself is essentially unchanged. Further, contrary to our expectations, no catastrophic ionization occurs. On the other hand, if all the nonlinear effects except for multi-photon absorption are switched off, the maximum density of plasma approaches saturation density ns = 10 ncr (Fig. 3c). The shape of the damaged region also changes considerably. We conclude that although multi-photon absorption is an important effect, it is not the main limiting mechanism of the laser intensity inside the plasma, at least in our considered parameter regime.
Figure 3d shows the plasma pattern for all the effects accounted for but with Γ = 0, i.e., for collisionless plasma. It will be shown below that collisional energy dissipation is responsible for a large part of absorbed energy by the medium. However, as follows from Fig. 3a,d, this energy dissipation has a little effect on the maximum plasma density, though it does increase the total number of electrons.
If the Kerr effect is turned off (not shown), the plasma shape remains close to the one shown in Fig. 3a, however with the maximum density decreased by ∼ 50%. This is consistent with what we would expect by removing self-focusing.
Figure 4 shows the details of interaction process. The left column of Fig. 4 shows the absolute value of the Poynting vector, |S⃗|, and the right column the plasma density at t = 300 fs (top), t =320 fs (middle), and t =340 fs, where t is the time after the start of the simulation. As the laser focuses, its intensity grows to an above-threshold value already at the pulse leading edge, corresponding to |S⃗|max ≈ (2 ÷3) × 1013 W/cm2 (t = 300 fs in Fig. 4). As a result, plasma is generated. As the plasma density grows, the laser field appears to be expelled from the damaged region, and thus the intensity profile forms a ring-like structure in the transverse cross-section (t = 320, 340 fs). Due to this modification the laser intensity decreases to approximately the threshold value and as a result the ionization process stops. The ring-like intensity merges back to a spot behind the generated plasma (t = 340 fs), but intensity there is still below the ionization threshold, and almost no new plasma is generated.
There are two possible reasons of the laser pulse being expelled from the region occupied by plasma: 1) as plasma forms, it acts as a diverging lens for the remaining part of the pulse, and thus the pulse gets defocused; 2) due to the relatively large damping factor Γ the plasma quiver motion is converted to heat at a fs time scale, and thus the laser energy dissipates inside the underdense plasma.
To find out which mechanism is dominant, let us refer to the energy diagnostic shown in Fig. 5. The actual simulations for the energy diagnostic were done in a larger simulation domain than that of Figs. 3, 4 to fit the entire laser pulse into it. The diagnostic in Fig. 5 shows the electromagnetic energy over the whole domain ℰem, total ionization-absorbed energy ℰabs, instantaneous kinetic energy of all the free particles in the domain ℰk and their sum ℰem + ℰabs + ℰk.
At t = 0, there is no electromagnetic field in the domain, the total energy is zero. At t > 0 the pulse starts entering the domain from its left boundary. By t ≈ 260 fs, most of the pulse is inside the domain. At t ≈ 290 fs, the ionization process starts. Eventually the total energy starts to decrease, indicating conversion of the particle kinetic energy into heat. At the end of simulation (t = 425 fs), the total electromagnetic energy in the domain has decreased approximately two times, with ∼ 1/3 of this change caused by the multi-photon absorption and ∼ 2/3 of the change by the thermal energy dissipation.
Taking into consideration that for a collisionless plasma (Fig. 3d) the morphology of the ionized medium is almost unchanged, we conclude that among the three competing effects the most important effect responsible for the shape of the damaged region is the defocusing of the laser pulse by the generated plasma, partially balanced by Kerr effect. An important effect responsible for dissipation of a large portion of the laser energy is the collisional damping of the quiver motion of plasma particles, and the least important effect is the laser energy depletion due to multi-photon absorption.
The interaction scenario is, in this way, temporally asymmetric. Whereas the leading edge of the pulse propagates in the dielectric medium almost unchanged, the peak and the trailing edge of the pulse experience a strong defocusing by the generated plasma. The intensity of the defocused peak of the laser pulse is still high enough to cause the broad low density plasma at back end, as seen in the the left sides of the plasma density profiles in Figure 4 at t = 320 and 340 fs. The rest of the defocused pulse does not contribute to additional plasma generation. A portion of the first half of the pulse that has created the plasma and does not experience defocusing continues propagating undisturbed in the forward direction. This results in an asymmetric plasma profile, elongated and sharpened in the forward direction, as shown in Figs. 3a and 4 at t = 340 fs.
Figure 6a shows the plasma density after the passage of the laser along its axis. The plasma density spans a distance ∼ 15μm, a figure on the order of the laser Rayleigh length, with the maximum value before the geometrical focus of the focusing mirror. The longitudinal profile is characterized by a relatively sharp growth (x ≈ −7 μm in Fig. 6a), a pedestal (−6 μm ≲ x ≲ −3 μm), an extremum (x = 0) and a relatively slow decrease (x > 0). Of a practical interest is the transverse size of the damaged region. This size vs. longitudinal coordinate is given in Fig. 6b, where each plot gives the width defined for constant electron density values, as indicated in the figure caption. This provides an effective contour plot of the plasma.
4. Laser polarization effect
A consequence of importance of the interaction with generated plasma is asymmetry that we have observed in the transverse cross-section. Figure 7a shows the plasma density versus distance, for the directions along the laser polarization and perpendicular to it, in the plane parallel to the laser focal plane and passing through the maximum plasma density point.
The plasma density along the laser polarization decreases faster with the radius than that in the direction perpendicular to the laser polarization. The plasma patterns appear deformed, with more plasma along the direction perpendicular to the laser polarization. There is an asymmetry of the focused laser spot intensity in the transverse direction due to the longitudinal field components near the laser focus , . However, this does not explain our observations, since this asymmetry is in the opposite direction, and for NA = 0.65 is nevertheless very small (see Fig. 7b). We conclude that the plasma pattern asymmetry is caused entirely by the field enhancement due to the plasma refraction index smaller than unity. We describe this below with a simple model.
Consider a dielectric cylinder of radius r, having its axis along x̂, located in space filled with another dielectric. Let the cylinder be composed of dielectric with permittivity ɛ2, and the background dielectric have permittivity ɛ1. Let also the field far from the cylinder be parallel to ŷ and the origin be located at a point at the cylinder axis. Then it can be shown thatEq. (20), if ɛ2 < ɛ1, as is the case for a plasma in a dielectric background, E1 < E2. Thus there is a field enhancement along direction perpendicular to the external field. This is what we observe in our simulations. Although the interaction problem is electromagnetic rather than electrostatic, this reasoning should be approximately valid for a cylinder radius smaller than the laser wavelength.
In this way, the laser propagating in the dielectric creates plasma, and this plasma in turn diverts the laser field. Due to Eq. (20), the diverted field is asymmetric, with enhancement along the direction perpendicular to laser polarization (i.e., ẑ for a laser field polarized along ŷ). This field enhancement produces more plasma in the ẑ direction. This process results in the asymmetry of the generated plasma pattern shown in Fig. 7a.
We have generally found that the extent of asymmetry increases with higher NA. We will quantitatively characterize this extent in a future study.
In the present study we have discussed the importance of particular effects for formation of a damaged region during propagation of a short intense laser in a bulk dielectric, for material parameters corresponding to fused silica. We have developed a 3D numerical model that accounts for multi-photon ionization, electromagnetic response of the resulting plasma, energy dissipation due to the plasma electron collisions, laser energy depletion due to the photon absorption, and Kerr effect. The analysis given in the paper has demonstrated that, unlike 1D case, the most important effect limiting the laser intensity in the 3D geometry is the refraction by the created plasma lens. The energy absorbing processes (heating and photon absorption), although are responsible for absorption of a considerable part of laser energy, do not have a drastic effect on the morphology of the damaged region. We have also observed an asymmetry of the plasma pattern in the transverse direction. This effect was found to be caused by plasma dispersion.
We would like to acknowledge fruitful discussions with P. Corkum and M. Gertsvolf. This work was supported in parts by Natural Sciences and Engineering Research Council of Canada, Ontario Ministry of Research and Innovation, Canada Research Chairs program and Canada Foundation for Innovation.
References and links
1. D. Homoelle, S. Wielandy, A. L. Gaeta, N. F. Borrelli, and C. Smith, “Infrared photosensitivity in silica glasses exposed to femtosecond laser pulses,” Opt. Lett. 24, 1311–1313 (1999). [CrossRef]
2. V. R. Bhardwaj, E. Simova, P. P. Rajeev, C. Hnatovsky, R. S. Taylor, D. M. Rayner, and P. B. Corkum, “Optically produced arrays of planar nanostructures inside fused silica,” Phys. Rev. Lett. 96, 057404 (2006). [CrossRef] [PubMed]
3. E. N. Glezer and E. Mazur, “Ultrafast-laser driven micro-explosions in transparent materials,” Appl. Phys. Lett. 71, 882–884 (1997). [CrossRef]
4. R. R. Gattass and E. Mazur, “Femtosecond laser micromachining in transparent materials,” Nat. Photonics 2, 219–225 (2008). [CrossRef]
5. G. D. Valle, R. Osellame, and P. Laporta, “Micromachining of photonic devices by femtosecond laser pulses,” J. Opt. A: Pure Appl. Opt. 11, 013001 (2009). [CrossRef]
6. L. Sudrie, A. Couairon, M. Franco, B. Lamouroux, B. Prade, S. Tzortzakis, and A. Mysyrowicz, “Femtosecond laser-induced damage and filamentary propagation in fused silica,” Phys. Rev. Lett. 89, 186601 (2002).
7. J. R. Peñano, P. Sprangle, B. Hafizi, W. Manheimer, and A. Zigler, “Transmission of intense femtosecond laser pulses into dielectrics,” Phys. Rev. E 72, 036412 (2005). [CrossRef]
9. A. Q. Wu, I. H. Chowdhury, and X. Xu, “Femtosecond laser absorption in fused silica: numerical and experimental investigation,” Phys. Rev. B 72, 085128 (2005). [CrossRef]
10. C. L. Arnold, A. Heisterkamp, W. Ertmer, and H. Lubatschowski, “Computational model for nonlinear plasma formation in high NA micromachining of transparent materials and biological cells,” Opt. Express 15, 10303–10317 (2007). [CrossRef] [PubMed]
11. I. M. Burakov, N. M. Bulgakova, R. Stoian, A. Mermillod-Blondin, E. Audouard, A. Rosenfeld, A. Husakou, and I. V. Hertel, “Spatial distribution of refractive index variations induced in bulk fused silica by single ultrashort and short laser pulses,” J. App. Phys. 101, 043506 (2007). [CrossRef]
12. P. P. Rajeev, M. Gertsvolf, C. Hnatovsky, E. Simova, R. S. Taylor, P. B. Corkum, D. M. Rayner, and V. R. Bhardwaj, “Transient nanoplasmonics inside dielectrics,” J. Phys. B 40, S273–S282 (2007). [CrossRef]
13. 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). [CrossRef]
14. D. Grojo, M. Gertsvolf, H. Jean-Ruel, S. Lei, L. Ramunno, D. M. Rayner, and P. B. Corkum, “Self-controlled formation of microlenses by optical breakdown inside wide-band-gap materials,” App. Phys. Lett. 93, 243118 (2008). [CrossRef]
15. A. Taflove and S. C. Hagness, Computational Electrodynamics, 3rd. ed. (Artech House, 2005), pp. 58–79.
16. K. S. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. AP-14, 302–307 (1966).
17. L. V. Keldysh, “Ionization in the field of a strong electromagnetic wave,” Sov. Phys. JETP 20, 1307–1314 (1965).
18. NRL Plasma Formulary (2002), p. 28.
19. C. A. Brau, Modern Problems in Classical Electrodynamics (Oxford Univ. Press, 2004), pp. 342–347.
20. A. Taflove and S. C. Hagness, Computational Electrodynamics, 3rd. ed. (Artech House, 2005), pp. 186–213.
21. K. I. Popov, V. Yu. Bychenkov, W. Rozmus, R. D. Sydora, and S. S. Bulanov, “Vacuum electron acceleration by tightly focused laser pulses with nanoscale targets,” Phys. Plasmas 16, 053106 (2009). [CrossRef]
22. J. A. Stratton and L. J. Chu, “Diffraction theory of electromagnetic waves,” Phys. Rev. 56, 99–107 (1939). [CrossRef]
23. S. Quabis, R. Dorn, M. Eberler, O. Glöckl, and G. Leuchs, “Focusing light to a tighter spot,” Opt. Commun. 179, 1–7 (2000). [CrossRef]
24. K. I. Popov, V. Yu. Bychenkov, W. Rozmus, and R. D. Sydora, “Electron vacuum acceleration by a tightly focused laser pulse,” Phys. Plasmas 15, 013108 (2008). [CrossRef]