This study presents a novel numerical model for laser ablation and laser damage in glass including beam propagation and nonlinear absorption of multiple incident ultrashort laser pulses. The laser ablation and damage in the glass cutting process with a picosecond pulsed laser was studied. The numerical results were in good agreement with our experimental observations, thereby revealing the damage mechanism induced by laser ablation. Beam propagation effects such as interference, diffraction and refraction, play a major role in the evolution of the crater structure and the damage region. There are three different damage regions, a thin layer and two different kinds of spikes. Moreover, the electronic damage mechanism was verified and distinguished from heat modification using the experimental results with different pulse spatial overlaps.
©2013 Optical Society of America
The display industry is currently trending towards thinner and lighter devices such as smart phones and tablet computers. Thin glass sheets (typical thickness of < 1 mm), which are widely used in almost all types of flat panel displays, are mechanically delicate and require precise and controllable cutting. Laser ablation with ultrashort pulses is a promising solution for precise glass cutting because of the following advantages. Micromachining of transparent dielectrics [1, 2] based on ultrashort-pulsed laser ablation can provide a localized material removal in the micro/nanometer scale through nonlinear absorption of photons. The separation of energy deposition and material removal in time reduces the destructive thermal and thermo-mechanical effects.
There are two characteristic aspects in ultrashort pulsed laser glass cutting, namely, ablation and damage beneath the ablated crater. The observable qualities of laser ablation are the crater profile and the ablation rate. Laser damage, such as refractive index modification [3–5], molten region  and micro-cracks , is crucial for the quality of the glass product, such as the edge strength. Therefore, the laser-glass interaction in glass cutting with picosecond pulsed lasers must be thoroughly investigated to control laser ablation and avoid laser damage.
Several numerical models [8–10] for laser ablation with ultrashort pulses have been presented and the critical free-electron density for laser breakdown has been widely used as an ablation criterion. For laser ablation with a single pulse, time-resolved ultrafast radiative transport equation (TE) was used in  to simulate the energy deposition and predict the crater shape. The TE model, however, cannot describe the diffraction and refraction of the following laser pulse induced by the ablated crater, and is thus limited to the single-pulse ablation. For the multi-pulse case, Vázquez de Aldana et al. [9, 10] introduced a numerical model for the propagation of the ablation channel in dielectrics with femtosecond laser pulses based on the integration of the scalar wave equation. However, this model is preferred when the laser pulse is shorter than 50 fs. Enormous computational power and time are required for longer pulses because of the high resolution in space and time. Therefore, the time-resolved wave equation is not favorable for laser ablation with multiple picosecond pulses.
Meanwhile, laser-induced damage in bulk dielectrics has been widely investigated both experimentally and theoretically [3, 11–14]. In the formation of bulk damage, the nonlinear ionization and the nonlinear beam propagation effects, such as filamentation, are the most important. However, the surface damage mechanism and the underlying physical phenomena in the laser ablation process have not yet been studied extensively.
The current study presents a numerical model to study both laser ablation and laser damage in transparent dielectrics with multiple picosecond laser pulses. The code combines a beam propagation model (BPM) and a plasma model. The former is based on the paraxial Helmholtz equation, whereas the latter describes the nonlinear absorption of photons in dielectrics. Simulations were performed to study the evolution of the ablated crater and the modified/damaged region induced by multiple incident laser pulses in the glass.
2. Numerical model
In the numerical ablation model, the laser energy deposition in the sample is determined in BPM coupled to a plasma model. In BPM, the paraxial Helmholtz equation is used. The plasma generation is described by a rate equation for free-electron density. The ablated crater and the damage region are defined by the resulting free-electron density distribution from the model.
2.1 Beam propagation model (BPM)
The laser beam propagates into the glass sample in the z direction. The x–y plane is perpendicular to the z direction. The electric field is expressed asEq. (1) into the Helmholtz equation and assume that the envelope slowly varies in the z direction, a paraxial Helmholtz equation can be derivedEq. (2), the dielectric function ε  of the material is given by9]. The damping constant χ is the reciprocal of electron relaxation time, which varies from 0.1 fs to 10 fs in transparent dielectrics [15, 16]; thus the maximum value of the damping constant χmax is 10 fs−1. χ saturates to χmax for greater free-electron densities. From the calculation of the dispersion data of the sample in the visible spectrum  with Eq. (3), ρB = 1.46 × 1023 cm−3 and ω0 = 19 fs−1.
The response of the material to the incident radiation is determined by the complex dielectric function in Eq. (3), the real and imaginary parts of which represent refraction and loss terms, respectively. Thus, the diffraction and refraction of the electric field caused by the crater wall, crater entrance and laser-induced solid-state-plasma are taken into account in the BPM model, as shown in Eq. (2).
2.2 Plasma model
In transparent dielectrics, the electrons in the valence band can be excited as free electrons by strong electric fields via photo-ionization. Then it is heated to even higher energy levels by the inverse bremsstrahlung absorption to initialize the cascade ionization. Meanwhile, there are some losses for the free-electron density due to electron-hole recombination and electron diffusion. In the picosecond scale, the recombination dominates and the diffusion can be ignored. The temporal evolution of free-electron density ρ is described by a rate equation of the generic form [18–20]
Based on the Keldysh theory of electronic ionization  in condensed media, multi-photon ionization (MPI) and tunneling ionization take place at high frequency and low frequency limits, respectively. The Keldysh parameter γ = ω(cε0n0mΔ/I)1/2/e is a quantitative indicator of the regime in which the ionization processes may occur. In this equation, m is the reduced mass of the electron and the hole (m = me/2), I is the laser intensity, and c is the speed of light in vacuum. Figure 1 shows the comparison of the photo-ionization rates obtained with different models. When γ ≥ 1, the Kennedy’s approximation in the MPI limit is consistent with the complete Keldysh model. For our study, the laser intensity was up to 1013 W⋅cm−2 and γ is larger than 1.8, thus the MPI coefficient is chosen in the Kennedy’s approximation 
The free electrons can absorb energy from the electric field by inverse bremsstrahlung absorption after being generated. When the kinetic energy of the free electron reaches the critical energy-level Ecrit = Δ, it will shortly produce a new free electron during collisions with atoms or molecules. The collisional losses and the excitation time of free-electrons from the bottom of the conduction band to the critical energy level are neglected. The cascade ionization coefficient  takes the reduced form as
The lifetime of laser-induced plasma in Soda Lime and K9 glasses  depends on the free-electron density, and the average lifetime is 100 ps when the plasma density is about 1020 cm−3. These glasses have similar characteristics ; thus, the alkaline earth boro-aluminosilicate glass has the same quality for laser-induced plasma as that of K9 glass. The electron recombination cannot be described by a characteristic constant lifetime. Therefore a two-body recombination mechanism is assumed in the plasma model. The effective lifetime tr of free electrons depends on the electron density, tr = 1/(ηrecρ). The recombination coefficient ηrec is taken as 1 × 10−10 cm3/s to obtain the plasma lifetime of 100 ps when the plasma density is 1020 cm−3.
Although the laser pulse shape along the pulse-propagating direction in the sample is modified due to the time-dependent nonlinear absorption [8, 14], for simplicity, the pulse shape is assumed as a Gaussian profile I(t) = I0 exp[-4 ln2(t/tp)2], where tp is the pulse duration and I0 is the peak intensity given by I0 = cε0n0|E|2. Based on this assumption and Eq. (4), the maximum free-electron density within the pulse duration only depends on the peak intensity and an interpolated function ρmax = F(I0) can be obtained in a wide range of laser intensities. With this function, the temporal evolution of the electric field can immediately provide the maximum electron density. Thus, the required computational power for pulse propagation with BPM is greatly reduced.
2.3 Ablation and damage criteria
In ultrashort pulsed laser ablation, energy deposition (heating electrons during the pulse) and material removal (heating lattice after the pulse) are separated in time. The electric field excites the electrons via MPI or tunneling ionization, followed by avalanche ionization. Free-electron density increases until reaching the critical electron density  ρcrit = (ω2meε0)/e, which leads to strong absorption, then reflection and ultimately a breakdown. Laser ablation induced by electron energy release takes place after the laser pulse passed, allowing the use of the critical free-electron density as the ablation criteria in the model. That is, the material will be removed in the region where the free-electron density exceeds the critical value, ρcrit = 3.95 × 1021 cm−3.
The material near the ablated-crater wall will be modified due to the energy released by high-density free-electrons. In the bulk of fused silica, the free-electron density corresponding to the damage threshold is close to 1020 cm−3 [3, 11–13]. The damage tracks are formed in the zone where the electron density is above this threshold electron density. For boro-aluminosilicate glass sample in this study, the threshold electron density for laser damage is assumed in the same order of magnitude as that for fused silica. Thus, the 1020 cm−3 value of the free-electron density, which is 0.025 ρcrit, is used as the damage criteria in the model.
Figure 2 shows the structure of the numerical model for laser ablation and laser damage; the model has three steps. In step 1, the maximum density of free-electron was pre-calculated and interpolated as a function of laser intensity, ρmax = F(I0), using Eq. (4). In step 2, coupled with the plasma model and the electron-density dependent dielectric function Eq. (3), the beam propagation and energy deposition were resolved within Eq. (2). The distribution of free-electron density inside the material was then calculated by the interpolated function ρmax = F(I0). In step 3, the ablation and damage criteria were used to determine the ablated crater and the damage region, respectively. The resulting crater surface serves as a new air-solid interface for the next pulse. For each laser pulse, the second and the third steps were repeated, and the damage region in the bulk was accumulated until it was ablated.
The simulation parameters were obtained from our experiments, where a frequency-doubled mode-locked laser system was used. The laser system delivers 40 μJ, 10 ps laser pulses at 532 nm with a repetition rate of 400 kHz. With an objective of 63 mm focal length, the laser beam was focused on the surface of the glass sample with a beam size of 6.5 μm in radius. The sample is an alkaline earth boro-aluminosilicate glass (Corning Eagle XG®) with a thickness of 0.3 mm. Its band-gap Δ is 3.45 eV, measured by the Tauc plot method . The refractive index n0 of the glass was 1.51 at the wavelength of 532 nm.
Several different spatial overlaps between the adjacent pulses (0% to 90%) can be obtained in the experiments by controlling the scan speed of the laser beam. For the simplest case, the scan speed of the laser beam is chosen to be large enough to ensure that no overlap between the adjacent pulses (if not specified) exists. Several passes on top of each other, with a time interval of several milliseconds, were selected to provide the multi-pulse behavior of laser ablation.
3. Results and discussion
3.1 Ablation and damage dynamics
Using the numerical model, the ablation processes for 10 incident pulses were obtained. Figures 3(a) -3(e) show the cross sections for laser intensity in air (a), laser intensity in the glass (b), free-electron density (c), crater profile (d), and damage region (e), respectively. The laser pulse propagates into the target dielectrics from the top in the normal direction and the initial planar surface of the glass located at z = 18 µm. The intensity distribution in Figs. 3(a) and 3(b) and the electron density distribution in Fig. 3(c) were generated in the second step of the model. The crater profile in Fig. 3(d) and the damage region in Fig. 3(e) were determined by the corresponding criteria in the third step, as described in Section 2.3.
The characteristics of the ablated crater are the most intuitive and of primary interest in the laser ablation process. The profiles and the fine structure of the ablated craters are shown in Fig. 3(d). The crater shape after the first pulse is nearly parabolic because of energy deposition through the initial planar interface and the Gaussian beam shape. This non-planar crater will reflect, refract and diffract the incident electric field when the second laser pulse arrives. The interference pattern of the laser field in the crater cavity, as shown in Fig. 3(a), induces the fine structure of the crater wall. The crater structure becomes more complex and pronounced for the subsequent laser pulses. Therefore, the interference effect of the laser beam is the main cause for the fine structure of the crater wall.
Figures 4(a) and 4(b) show the cross sections of the ablated craters and the damaged regions after 5 and 10 pulses in the experiment, which agree well with the numerical results shown in Figs. 4(c) and 4(d). The experimental results show that the crater width after 10 pulses was slightly larger than that after 5 pulses, whereas it was constant in the numerical results. The increasing of the crater width is due to the backward reflection on the crater wall, which induces a part of the laser beam irradiating the edge of the crater cavity and causes its ablation. In our model, however, the paraxial Helmholtz equation, which is restricted to the propagation of the forward propagating waves, cannot take into account the backward reflection.
In the experimental results shown in Figs. 4(a) and 4(b), we can observe a thin dark layer near the crater wall and some needle-like structures beneath the crater wall, namely, the damage region. Due to the laser-induced free-electron density and its energy release, the refractive index of the sample was permanently changed and some other defects were generated, such as color centers , influencing the interaction of the material with the next incident pulse. In the numerical simulation, using the free-electron density of 1020 cm−3 as the damage criteria, the damage regions shown in Figs. 4(c) and 4(d) were well reproduced compared with the experimental results. As shown in Fig. 3(e), a regular 2 μm-thick layer of damage region is generated near the crater surface after the first pulse. For the subsequent pulses, refraction and diffraction on the crater wall play a significant role in the propagation of the electric field and drastically change the energy deposition in the material. As a result, the damage regions become different from that after the first pulse, and several needle-like damage regions in a length of several microns were formed beneath the crater wall and surrounding the whole crater. In particular, two extremely long spikes of about 10-µm depth inside the sample body were generated by the crater edges after the second pulse, which can also be observed in the experimental results shown in Fig. 4(a).
Comparing the experimental and numerical results, three kinds of damage regions exist in the sample: a thin damage layer, some spikes beneath the ablation surface, and two long spikes from the crater edges. The thin layer is the penetration zone formed by the nonlinear absorption of laser energy. This zone will be modified after the energy of free-electrons is transferred to the lattice. The needle-like spike regions are formed by the diffraction of the electric field in the crater. The two long spikes, induced by the refraction of laser beam at the crater edges, are at a small angle with respect to the forward direction. The latter two damage regions are possible sources of crack formation in the later application of glass piece, ultimately leading to a decreased fracturing strength of the glass product.
In order to check the impact of beam propagation effects such as interference and diffraction on the evolution of the crater and the damage region, we ignore the transverse Laplacian part Δ⊥E in Eq. (2) and get the transport equation (TE)Equation (7) was also implemented with the same plasma model to simulate the ablation process with multiple picosecond laser pulses. The results are shown in Fig. 5 . The crater profiles obtained with the TE model after 1, 5, and 10 pulses are quite smooth as shown in Fig. 5(a). Only a thin and smooth penetration layer of the damage region is obtained as shown in Figs. 5(b)-5(d) while the crater profiles and the damage regions obtained with the BPM model are bumpy and noisy. Because the TE model cannot describe the diffraction and refraction of laser beam on the crater, the beam propagation effect is crucial in the ablation and damage dynamics.
3.2 Electronic damage region and heat affected region
In Section 3.1, the good agreement of the damage region in the experimental results and the numerical ones reveals that the damage beneath the crater wall is mainly induced by the high electron density and the subsequent atomic/molecular processes. In order to verify this electronic damage mechanism and distinguish it from heat modification, the effect of energy accumulation on the damage region was investigated experimentally by taking different spatial overlaps of the adjacent pulses.
With a repetition rate of 400 kHz and pulse energy of 40 μJ (50 μJ for 54% case), the overlap was adjusted from 0% to 90%, and the spatial distance between the adjacent pulses changed from 18 μm to 1.25 μm with the same time interval of 2.5 µs. The experimental cross sections of the damage regions are shown in Fig. 6 . With the overlaps of 0% and 54%, the contours of the damage region are irregular, with some similar spikes as that of the numerical results shown in Figs. 4(c) and 4(d). This characteristic of the damage region indicates that the electronic damage was dominant in the damage mechanism, and is represented by the dark black region. A description of the damage regions with the overlaps of 85% and 90% are given in Figs. 6(c) and 6(d). A broader affected zone and a smoother contour are observed, indicating that the large-area modified zone is due to a diffusion process. The energy accumulation becomes significant for higher spatial overlaps, and heat diffusion leads to a wide heat-modified zone, whereas the electronic damage (dark black region) is still present and just covered by the diffusion region (light black region). Therefore, the different pulse overlaps and the energy accumulation will change the relative role of electronic damage and heat modification in the evolution of the damage zone.
This study has focused on the dynamics of laser ablation and laser damage in glass with multiple incident picosecond pulses. The mechanism of laser ablation and laser damage was investigated numerically using the presented ablation model in transparent dielectrics. The ablated crater and the damage morphology agree well with our experimental results. The laser-induced damage in glass cutting depends essentially on the characteristics of beam propagation in the material bulk, especially the mechanism of refraction and diffraction. Three damage mechanisms were revealed in the material: regular energy deposition leads to a penetration layer beneath the crater wall, whereas the spikes were identified as bulk damage induced by peaks of intensity patterns, which themselves create peaks of electron density. Finally, the electronic damage was verified and compared with the heat-affected zone using the experimental results at different spatial overlaps of the adjacent pulses. The electronic damage region was irregular due to the beam propagation effect, while the heat-affect zone was the characteristic thermal diffusion region and much smooth. The electronic damage is covered by heat modification at high pulse overlaps due to heat accumulation effect.
References and links
1. R. R. Gattass and E. Mazur, “Femtosecond laser micromachining in transparent materials,” Nat. Photonics 2(4), 219–225 (2008). [CrossRef]
2. G. Della Valle, R. Osellame, and P. Laporta, “Micromachining of photonic devices by femtosecond laser pulses,” J. Opt. A, Pure Appl. Opt. 11(1), 013001 (2009). [CrossRef]
3. 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(18), 186601 (2002). [CrossRef] [PubMed]
5. L. Shah, A. Y. Arai, S. M. Eaton, and P. R. Herman, “Waveguide writing in fused silica with a femtosecond fiber laser at 522 nm and 1 MHz repetition rate,” Opt. Express 13(6), 1999–2006 (2005). [CrossRef] [PubMed]
6. A. Ben-Yakar, A. Harkin, J. Ashmore, R. L. Byer, and H. A. Stone, “Thermal and fluid processes of a thin melt zone during femtosecond laser ablation of glass: the formation of rims by single laser pulses,” J. Phys. D Appl. Phys. 40(5), 1447–1459 (2007). [CrossRef]
7. H. Varel, D. Ashkenasi, A. Rosenfeld, M. Wähmer, and E. E. B. Campbell, “Micromachining of quartz with ultrashort laser pulses,” Appl. Phys., A Mater. Sci. Process. 65(4-5), 367–373 (1997). [CrossRef]
8. L. Jiang and H. L. Tsai, “Prediction of crater shape in femtosecond laser ablation of dielectrics,” J. Phys. D Appl. Phys. 37(10), 1492–1496 (2004). [CrossRef]
9. J. R. Vázquez de Aldana, C. Méndez, L. Roso, and P. Moreno, “Propagation of ablation channels with multiple femtosecond laser pulses in dielectrics: numerical simulations and experiments,” J. Phys. D Appl. Phys. 38, 3764 (2005).
10. J. R. Vázquez de Aldana, C. Méndez, and L. Roso, “Saturation of ablation channels micro-machined in fused silica with many femtosecond laser pulses,” Opt. Express 14(3), 1329–1338 (2006). [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. Appl. Phys. 101(4), 043506 (2007). [CrossRef]
12. D. G. Papazoglou and S. Tzortzakis, “Physical mechanisms of fused silica restructuring and densification after femtosecond laser excitation,” Opt. Mater. Express 1(4), 625–632 (2011). [CrossRef]
13. S. S. Mao, F. Quéré, S. Guizard, X. Mao, R. E. Russo, G. Petite, and P. Martin, “Dynamics of femtosecond laser interaction with dielectrics,” Appl. Phys., A Mater. Sci. Process. 79, 1695–1709 (2004). [CrossRef]
15. D. Arnold, E. Cartier, and D. J. DiMaria, “Acoustic-phonon runaway and impact ionization by hot electrons in silicon dioxide,” Phys. Rev. B Condens. Matter 45(3), 1477–1480 (1992). [CrossRef] [PubMed]
16. J. R. Gulley, S. W. Winkler, W. M. Dennis, C. M. Liebig, and R. Stoian, “Interaction of ultrashort-laser pulses with induced undercritical plasmas in fused silica,” Phys. Rev. A 85(1), 013808 (2012). [CrossRef]
17. EAGLE XG® Material Information Sheet. http://www.corning.com/displaytechnologies/en/products/eaglexg/index.aspx.
18. A. Vogel, J. Noack, G. Hüttman, and G. Paltauf, “Mechanisms of femtosecond laser nanosurgery of cells and tissues,” Appl. Phys. B 81(8), 1015–1047 (2005). [CrossRef]
19. P. K. Kennedy, “A first-order model for computation of laser-induced breakdown thresholds in ocular and aqueous media: Part I—Theory,” IEEE J. Quantum Electron. 31(12), 2241–2249 (1995). [CrossRef]
20. Y. Wang, Y. Zhao, J. Shao, and Z. Fan, “Effect of native effect and laser-induced defects on multi-shot laser-induced damage in multilayer mirrors,” Chin. Opt. Lett. 9(9), 093102–093105 (2011). [CrossRef]
21. L. V. Keldysh, “Ionization in the field of a strong electromagnetic wave,” Sov. Phys. JETP 20, 1307–1314 (1965).
22. S. Quan, J. Hong-Bing, L. Yi, Z. Yong-Heng, Y. Hong, and G. Qi-Huang, “Relaxation of dense electron plasma induced by femtosecond laser in dielectric materials,” Chin. Phys. Lett. 23(1), 189–192 (2006). [CrossRef]
23. Types of Glass.http://www.cmog.org/article/types-glass.
24. J. Tauc, R. Grigorovici, and A. Vancu, “Optical properties and electronic structure of amorphous germanium,” Phys. Status Solidi 15(2), 627–637 (1966). [CrossRef]
25. D. Ashkenasi and A. Lemke, “Picosecond laser-induced color centers in glass optics,” J. Laser Appl. 23(1), 012007 (2011). [CrossRef]