Maxwell’s wave equation was solved for fs laser drilling of silicon. The pre-formed hole wall’s influence on the propagation behavior of subsequent laser pulses was investigated. The laser intensity at hole bottom shows distinct profile as compared with that at hole entrance. The multi-peaks and ring structure of the laser intensity were found at hole bottom. The position of maximum laser intensity (MLI) in relation to the wall taper angle was studied. It was found that the position of the MLI point would be closer to the hole entrance with increasing taper angle. This observation provides valuable information in predicting the position of plasma plume which is a key factor influencing laser drilling process. The elliptical entrance hole shape and zonal structure at the hole bottom reported in the literatures have been reasonably explained using the laser intensity distribution obtained in the present model.
© 2015 Optical Society of America
For percussion laser drilling, the material inside the micro hole is removed pulse by pulse. However, in the laser drilling process, the taper , barreling shape  and tail shape  of the hole can hardly be avoided owing to the power intensity distribution of the Gaussian beam and plasma expansion when the hole is going deeper. Because of the electromagnetic nature of the fs laser, the hole geometry formed by the previous laser pulses will interact with the successive pulses. The absorption, reflection, scattering, transmission at every point of the hole wall will significant influence in the following propagation behavior of laser pulse. The various shapes of the hole wall with different angles may act like wave guide  for laser beam or just cause energy loss with absorption. Accordingly, due to the multi-reflection on the hole wall, the energy density distribution of the laser pulse inside the hole exhibits distinctive structures with Gaussian distribution at the orifice of the hole. In the later stage of the laser percussion drilling, the amount of laser energy that can reach the bottom of the hole decides the laser drilling depth. In other words, the aspect ratio and drilling efficiency of the hole is highly related to the hole geometry at the beginning of the laser ablation. However, direct experimental study of the hole wall’s influence on the laser drilling depth and laser energy propagation is rare because of the dynamic nature of the laser drilling process and the lack of in situ characterization tools for the micro features. Thus, it is of great importance to conduct a simulation to investigate the influence of geometrical effects on the laser beam propagation and subsequent laser ablation. In order to fundamentally understand the laser beam propagation behavior in the micro hole, it is essential to solve the Maxwell’s equations which describe the nature of the laser beam.
There have been a number of simulation works [5, 6] using various theories that make contributions to study the hole wall’s influence on the propagation of the laser beam in the micro hole drilling. From the previous works, it can be concluded that it is essential to solve Maxwell’s equations to fully investigate the laser propagation behavior inside the micro-hole. Only by this numerical way, the laser energy distribution inside the micro-hole can be fundamentally understood. However, so far, we have not found any of the modeling works that focus on the investigation of the taper angle’s influence on the laser beam propagation by solving the Maxwell’s equations. This will thus be the main content of this paper.
2. Numerical model
In this study, the steady state analysis of energy distribution of fs pulse laser beam in the micro hole was conducted by numerically solving time harmonic Maxwell’s wave equation (Eq. (1)) . FEM analysis was carried out using COMSOL Multiphysics software . Moreover, the taper angle of hole wall was varied to investigate the influence of multi-reflection and thereby interference of the laser beam on the laser energy distribution at the bottom of the micro hole in silicon. Figure 1 shows the schematic view of the model setup.
There are two regions in computational domain, namely air domain and silicon domain representing the micro-hole in the silicon wafer. In order to reduce the computational time, the thickness of silicon layer (in the both the horizontal and vertical directions) is assumed to be as small as 0.2 µm. The incoming laser beam was applied at the top (boundary A).8, 9]. In this work, a fs laser with the wavelength of 775nm and pulse duration of 250fs is assumed. The laser beam is assumed as linear polarized transverse electric wave which was focused by a lens at the top of computational domain (boundary A). The fs laser spot size on the substrate is assumed as 3.2 µm. The polarization direction is along the x direction in Fig. 1. Boundary B is the interface between air and silicon (Fig. 1). In order to predict the reflection and refraction at this boundary, the tangential component of electric field across the interface is continuous. Meanwhile, tangential component of magnetic field across the interface is also continuous because of the assumption of zero surface charges. The boundary conditions (BCs) at boundary B can be mathematically expressed as Eqs. (2) and (3)  for electric and magnetic field, respectively. is the unit normal vector. The subscripts ‘1’ & ‘2’ denote the respective air domain and silicon domain.11] expressed in Eq. (4) is selected .12, 13]. A finer mesh is required in the silicon domain due to the shortened wavelength when light is coming into the solid material. Thus, in the case of silicon domain, the element size is less than , where is the refractive index of silicon at 775nm. The input average laser intensity used in the model is 1 × 1012 W/cm2. This laser intensity is large enough to induce the nonlinear absorption of laser energy by the silicon. Thus, we consider that the refractive index is changed by varying the intensity of laser. The nonlinear refractive index is written as: , where is the constant part of the refractive index, (3 × 10−14 cm2/W ) is the Kerr coefficient of silicon, is the laser intensity.
3. Result and discussion
In most of the previous studies [16, 17] involving the 3-D or 2-D axial symmetry  modeling of laser ablation, the energy of laser source term is commonly considered as Gaussian distribution. However, what exactly happens for the laser power distribution inside the hole and at the bottom surface is still an unsolved problem. In this study, electric field and magnetic field are obtained by solving the Maxwell’s wave equation inside the micro-hole of silicon wafer. Figure 2 shows the laser intensity by time averaging at the bottom of the micro-hole (z = 0.75 µm) with varying the taper angle. As shown in Fig. 2, the laser intensity distribution varies dramatically when the taper angle increases from 0 degree to 53 degree. The laser beam intensity pattern turns to be very complex when the taper angle is increasing to above 46 degree. The refractive index of the single crystal silicon is 3.71  at a wavelength of 775 nm which indicates that the reflection of optical wave at the air/silicon interface is much stronger than the other transparent material such as glass and polymer . Thus, a considerable portion of the laser beam is reflected into the bottom by the hole wall. This result shows that reflection of the hole wall has a significant influence to the distribution of the laser energy inside the micro-hole. Figure 3(a) shows the maximum laser intensity (MLI) at the bottom of the micro-hole at various taper angles. MLI is defined as the maximum value of calculated laser intensity. As indicated in Fig. 3(a), the MLI firstly increases slowly with increasing the taper angle until it reaches the peak point at 34 degree. Then, the MLI shows a sharp decrease when taper angle is further increased to 53 degree. The critical taper angle in terms of laser beam propagation direction can be calculated using the Eq. (5) that is obtained from the geometrical relation of the entrance hole diameter , hole depth d and taper angle as shown in Fig. 3(c).Fig. 3(a). The deviation between the above two values may be due to increased energy loss on the hole wall when taper angle become larger. Figure 3(c) illustrates the comparison of multi-reflections of the laser beam on the hole wall with taper angles of 25, 33, 40 degree, respectively.
The laser wave that interacts with hole wall for taper angle smaller than 33 degree can propagate to the bottom of the hole and contribute to energy concentration at hole bottom. It is possible that the hole wall acts as waveguide which is helpful to deliver the laser energy to the bottom under these circumstances . However, when taper angle is larger than 33 degree, a large portion of the laser energy can be reflected back to the air or strongly absorbed by the side wall due to increased transmittance resulting from the small incident angle on the hole wall. This can lead to lower MLI as shown in Fig. 3(a). It also can be observed from Fig. 2 that the laser intensity shows a ring shape when taper angles are 46 and 53 degree. In these two cases, laser intensity at center of the bottom is much lower than the surrounding ring shape area that is contrary to the assumptions in most of the previous modeling works where the laser intensity is always assumed Gaussian distribution.
However, in some previous experimental studies , the zonal structure at the bottom of the micro-hole is also evident which suggests the assumption of Gaussian energy distributionat the hole bottom may not be accurate. On the other hand, the current result of the intensity distribution at bottom surface can be related to the zonal structure. The ring shape distributionof the laser intensity resulted from the reflection of the hole wall can lead to a higher material removal rate not at the center of the bottom but at the surrounding annularity area. This high material removal rate probably produces a deeper ring shaper groove at the area surrounding the center part of the bottom surface. Eventually, the zonal structure is formed at the hole bottom.
Additionally, Fig. 2 shows that the energy distribution inside the hole is even not axial symmetry. This is due to fact that the laser beam is a linear polarized plane wave which indicates that the component of electric field in polarized direction is stronger than that in other directions. The laser intensity distribution shows an ellipse shape that can lead to an ellipse shape of the entrance hole drilled by linear polarized laser which is used in this study. Thus, the ellipse shape of entrance hole that is experimentally observed in other studies [22–24] can be explained by the current numerical result. Figure 4 presents the evolution of spatial distribution of laser intensity at the x-z plane and y-z plane inside the micro-hole at various taper angles. The definitions of the x-z plane and y-z plane are illustrated in the Fig. 1.
As shown in Fig. 4, the laser intensity of the x-z plane and y-z plane at same taper angle are distinctive in the distribution pattern. This asymmetric laser intensity distribution in two orthogonal planes is resulted from the linear polarization of the laser beam and the following multi-refection of the plane wave on the hole wall. The maximum values of theenergy intensity in these two planes are however close to each other which indicate the point with maximum intensity in the 3-D space is located at the central axis of the 3-D model. The MLI locations obtained from the above 2-D observation fit well with the data that we extract from the 3-D calculation. In Fig. 3(b), the MLIs at y-z plane and x-z plane of the micro-hole at various taper angles are plotted. The two curves in Fig. 3(b) show a similar trend with the Fig. 3(a) that presents the MLI in the bottom surface. The peak points of the MLI at x-z and y-z planes appear at around 35 degree of the taper angle. This result suggests that there is strong correlation between the energy intensity inside the hole and outflow at the hole bottom surface. The reflection of laser beam by the hole wall at smaller taper angle is helpful to focus the laser energy to the central area and elevate the MLI. However, when taper angle becomes larger than the critical point (35 degree), a considerable portion of the energy is absorbed by the hole wall or repelled out the of the hole which leads to the decline of the MLI. In Fig. 4, it can be observed that the position of the maximum intensity point turns to be closer to the top plane when increasing the taper angle. It is reasonable to explain that the reflection waves from the hole wall at higher taper angle are more likely to focus at a higher position due to a smaller reflection angle (measured from the incoming beam direction). The previous experimental study  reported that the micro-hole wall has a strong influence on the laser induced plasma’s temperature, electron density and position. In the laser material interaction process, it is difficult to avoid the plasma  due to the fact that the laser energy is sufficient high to excite the electron to conduction bond and cause multi-ionization. one researcher  experimentally observed that the plasma plume can be generated by the pre-pulse (originates from the amplified spontaneous emission) several nanoseconds before the starting of the main fs pulse. The extinguishing time of plasma (a few microseconds ) induced by the pre-pulse is much longer that the interval between the pre-pulse and main pulse. Thus, it can be concluded than there is an interaction between the beam and plasma. After the generation of plasma plume near the ablation surface, the laser beam can consider to be the energy source for the heating of the plasma through the inverse bremsstrahlung . Theoretically, it is believed that the larger laser fluence in the spot area can generate more energetic plasma species . Experimentally, Li et al.  used fast photography and optical emission spectroscopy to investigate the dynamic of nanosecond laser induced plasma in the copper. It was found that the position of highest electron density point of plasma is mostly at the axial of laser beam. It thus suggested that on the plane parallel to the target surface, the maximum temperature point of plasma is located at the maximum laser intensity. Wang et al.  experimentally obtained the longitudinally resolved measurement of plasma density along fs laser filament. They reported that the plasma density was lower at the leading and tailing parts of the laser filament and reached a maximum near the geometrical focus of the laser beam. Thus, the maximum temperature point of plasma is mostly located at the maximum laser intensity point. The observation on the varying location of maximum laser intensity point at various taper angles in Fig. 2 and Fig. 4 is of great importance to provide valuable information in predicting the status of plasma plume inside the micro hole. From Fig. 4(i, j, k, l), it can be observed at the higher taper angle the energy absorption by the side wall is much stronger than that at lower taper angle.
In conclusion, theoretical result in analysis of the taper angle’s effect on the laser beam propagation inside the micro-hole of silicon wafer was obtained by numerically solving the time harmonic Maxwell’s wave equation. The MLI at hole bottom firstly increased with increasing taper angle, then reduced with increasing taper angle beyond the critical point. The relationship between the critical taper angle and geometrical dimension of micro-hole was derived based on the reflection behavior of the laser beam. The ellipse entrance hole shape and zonal structure at the hole bottom reported in the literatures have been reasonably explained using the laser intensity distribution obtained in the current model. In addition, the position of the maximum laser intensity point in the vertical direction plane was found to become higher with increasing the taper angle. This result provided important information for the laser induced plasma that is of great significance for both experimental and theoretical studies of the laser material interaction.
This work was supported by A*STAR AOP project (1223600005).
References and links
1. M. Ghoreishi, D. K. Y. Low, and L. Li, “Comparative statistical analysis of hole taper and circularity in laser percussion drilling,” Int. J. Mach. Tools Manuf. 42(9), 985–995 (2002). [CrossRef]
2. B. Tan, “Deep micro hole drilling in a silicon substrate using multi-bursts of nanosecond UV laser pulses,” J. Micromech. Microeng. 16(1), 109–112 (2006). [CrossRef]
3. S. Döring, S. Richter, A. Tünnermann, and S. Nolte, “Evolution of hole depth and shape in ultrashort pulse deep drilling in silicon,” Appl. Phys., A Mater. Sci. Process. 105(1), 69–74 (2011). [CrossRef]
4. M. N. Ruberto, X. Zhang, R. Scarmozzino, A. E. Willner, D. V. Podlesnik, and R. M. Osgood Jr., “Laser-controlled micrometer-scale photoelectrochemical etching of III-V semiconductors,” J. Electrochem. Soc. 138(4), 1174–1185 (1991). [CrossRef]
5. S. Tao, B. Wu, and S. Lei, “Study of laser beam propagation in microholes and the effect on femtosecond laser micromachining,” J. Appl. Phys. 109(12), 123501 (2011). [CrossRef]
6. M. Courtois, M. Carin, P. L. Masson, S. Gaied, and M. Balabane, “A new approach to compute multi-reflections of laser beam in a keyhole for heat transfer and fluid flow modelling in laser welding,” J. Phys. D Appl. Phys. 46(50), 50530501 (2013). [CrossRef]
7. User's Guide of COMSOL Multiphysics software, version 4.3b (2013).
8. G. R. Fowles, Introduction to Modern Optics (Dover Publications, 1989).
9. M. Schwartz, Principles of Electrodynamics (Dover Publications, 1987).
10. J. D. Jackson, Classical Electrodynamics (John Wiley & Sons, 1999).
11. W. Chen and Q. Zhan, “Numerical study of an apertureless near field scanning optical microscope probe under radial polarization illumination,” Opt. Express 15(7), 4106–4111 (2007). [CrossRef] [PubMed]
12. S. Marburg and S. Schneider, “Influence of element types on numeric error for acoustic boundary elements,” J. Comput. Acoust. 11(03), 363–386 (2003). [CrossRef]
13. S. Marburg, “Six boundary elements per wavelength: Is that enough?” J. Comput. Acoust. 10(01), 25–51 (2002). [CrossRef]
14. J. Leuthold, C. Koos, and W. Freude, “Nonlinear silicon photonics,” Nat. Photonics 4(8), 535–544 (2010). [CrossRef]
15. A. D. Bristow, N. Rotenberg, and H. M. van Driel, “Two-photon absorption and Kerr coefficients of silicon for 850–2200nm,” Appl. Phys. Lett. 90(19), 191104 (2007). [CrossRef]
16. A. Matsunawa and V. Semak, “The simulation of front keyhole wall dynamics during laser welding,” J. Phys. D Appl. Phys. 30(5), 798–809 (1997). [CrossRef]
17. V. V. Semak, B. Damkroger, and S. Kempka, “Temporal evolution of the temperature field in the beam interaction zone during laser material processing,” J. Phys. D Appl. Phys. 32(15), 1819–1825 (1999). [CrossRef]
18. H. Ki and J. Mazumder, “Numerical simulation of femtosecond laser interaction with silicon,” J. Laser Appl. 17(2), 110–117 (2005). [CrossRef]
19. E. D. Palik, Handbook of Optical Constants of Solids (Academic Press, 1999).
20. H. Y. Zheng, Y. C. Lam, C. Sundarraman, and D. V. Tran, “Influence of substrate cooling on femtosecond laser machined hole depth and diameter,” Appl. Phys., A Mater. Sci. Process. 89(2), 559–563 (2007). [CrossRef]
21. C. Y. Chien and M. C. Gupta, “Pulse width effect in ultrafast laser processing of materials,” Appl. Phys., A Mater. Sci. Process. 81(6), 1257–1263 (2005). [CrossRef]
22. K. Venkatakrishnan, B. Tan, P. Stanley, and N. R. Sivakumar, “The effect of polarization on ultrashort pulsed laser ablation of thin metal films,” J. Appl. Phys. 92(3), 1604–1607 (2002). [CrossRef]
23. J. Bonse, S. Baudach, J. Krüger, W. Kautek, and M. Lenzner, “Femtosecond laser ablation of silicon-modification thresholds and morphology,” Appl. Phys., A Mater. Sci. Process. 74(1), 19–25 (2002). [CrossRef]
24. T. Y. Choi and C. P. Grigoropoulos, “Plasma and ablation dynamics in ultrafast laser processing of crystalline silicon,” J. Appl. Phys. 92(9), 4918–4925 (2002). [CrossRef]
25. X. Zeng, S. S. Mao, C. Liu, X. Mao, R. Greif, and R. E. Russo, “Plasma diagnostics during laser ablation in a cavity,” Spectrochim. Acta B 58(5), 867–877 (2003). [CrossRef]
26. D. Devaux, R. Fabbro, L. Tollier, and E. Bartnicki, “Generation of shock waves by laser-induced plasma in confined geometry,” J. Appl. Phys. 74(4), 2268–2273 (1993). [CrossRef]
27. D. J. Hwang, C. P. Grigoropoulos, and T. Y. Choi, “Efficiency of silicon micromachining by femtosecond laser pulses in ambient air,” J. Appl. Phys. 99(8), 0831011 (2006). [CrossRef]
28. N. Smijesh and R. Philip, “Emission dynamics of an expanding ultrafast-laser produced Zn plasma under different ambient pressures,” J. Appl. Phys. 114(9), 0933011 (2013). [CrossRef]
29. J. J. Chang and B. E. Warner, “Laser-plasma interaction during visible-laser ablation of methods,” Appl. Phys. Lett. 69(4), 473–475 (1996). [CrossRef]
30. X. Li, W. Wei, J. Wu, S. Jia, and A. Qiu, “The Influence of spot size on the expansion dynamics of nanosecond-laser-produced copper plasmas in atmosphere,” J. Appl. Phys. 113, 243304 (2013).
31. T. J. Wang, J. Ju, Y. Wei, R. Li, Z. Xu, and S. L. Chin, “Longitudinally resolved measurement of plasma density along femtosecond laser filament via terahertz spectroscopy,” Appl. Phys. Lett. 105, 051101 (2014).