1. Introduction

Femtosecond lasers have high intensities which induce material phase and structure changes suitable for optics and photonics fabrication including laser polishing, optical array fabrication, surface nano-structuring, waveguide writing, and photonics-materials welding [1–5]. These applications are enabled either by non-thermal-ablation-based material removal via solid-vapor phase change with bulk temperatures remaining below the melting point and/or thermal and heat accumulation processes via solid-liquid-vapor phase changes. Each laser pulse induces complex physical phenomena which occur within the timescale of several picoseconds: energy absorption, free-carrier electron generation, and electron-lattice energy/temperature coupling and thermalization [6]. After laser pulse interaction, the induced heat diffuses throughout the material bulk until either a subsequent laser pulse strikes the surface or the initial bulk temperature is achieved, both of which normally occur on the microsecond timescale [7].

It is important to assess, differentiate, and control the non-thermal ablation and/or thermal impact of different combinations of laser parameters to ensure that laser processing conditions are aligned with the required interaction mechanism for a particular fabrication application. For thermal-melting-based applications such as laser micro-welding, heat generated during processing must be controlled to ensure precise, localized melting of the substrate to achieve crack-free, long-standing bonding with minimal thermal-gradient-induced stresses [5,8]. For non-thermal ablation-based surface processing applications, optimal laser parameters must be determined to enable direct material removal via solid-vapor phase change and to avoid compounded bulk heating by repetitive pulse incidence. This heat accumulation arises if the combined effect of laser parameters (fluence, repetition rate, and scanning speed) is not controlled to enable sufficient diffusion of heat between pulses and if it is not optimal for the thermal properties of the material. Uncontrolled heat accumulation can lead to the onset of melting, oxidation, and pileup of material [9,10], which are detrimental to structuring applications and interfere with determination of the true ablation rate and removal volume.

It is not a trivial task to determine optimal laser parameters to distinguish between non-thermal ablation-based material breakdown (solid-vapor transition) and heat-accumulation-based phase change (solid-liquid-vapor transition). Many material-dependent, time-consuming designs of experiments are required to assess both the isolated and combined impact of parameters on both ablation and bulk material heating. Optimal parameters must then be identified to enable precision laser-based fabrication processes.

An effective modeling method can instead be implemented to determine initial processing parameter combinations for achieving femtosecond-laser processing relying on different physical mechanisms. This modeling method must consider the complex phenomena on the pulse timescale which cause non-thermal ablation, as well as the heat accumulation process controlled by bulk heat diffusion between pulses. A variety of models – two-temperature, beam propagation, molecular dynamics, etc. – have been implemented to simulate laser interaction on the pulse timescale to describe lattice heating, ablation, and phase change during ultrafast laser-material interaction [11–13]. However, these models require many calibrated input parameters to describe the material’s response to laser interaction mechanisms. It is impractical to experimentally determine these parameters due to laser-parameter and temperature dependencies. Prediction of ablation and heat accumulation processes using these models requires extensive simulation times and advanced computational resources [14], which are unsuitable for simulating area processing with hundreds to thousands of laser pulses. Heat accumulation models can alternatively be used to predict heat accumulation during laser processing by using classical descriptions of heat generation and diffusion [9,15,16]. Although these models are time-efficient, they require processing-condition-specific input parameters to describe the amount of energy absorbed during the laser interaction process. They also cannot predict the onset of non-thermal ablation because they do not simulate the laser interaction phenomena occurring on the pulse timescale. An effective modeling method for predicting laser parameters to assess, differentiate, and control non-thermal ablation and thermal processing conditions during multi-pulse femtosecond laser processing remains to be developed.

Our overall objective is to establish a modeling method to determine optimal femtosecond laser processing parameters which enable non-thermal ablation or thermal melting and establish a path for process-specific parameter optimization. For non-thermal ablation applications, our modeling method aims to fulfill two goals: (1) predict initial laser parameters which enable non-thermal ablation, (2) predict the impact of laser parameters on heat accumulation and determine parameter combinations to avoid the onset of detrimental melting, oxidation, and material pileup. This method needs to combine accuracy and computational efficiency, be suitable for multi-pulse processing scenarios, and alleviate the need for time-consuming, iterative, experimental investigations of optimized laser-parameter combinations.

To develop this modeling method, we investigate and compare the prediction of heat accumulation in silicon due to multi-pulse femtosecond laser processing using both a three-dimensional two-temperature model (TTM) and a classical heat accumulation model [9]. We use the TTM to predict the laser fluence for achieving non-thermal ablation and investigate optimal laser parameters, e.g., repetition rate and scanning speed, which minimize heat accumulation. We investigate the link between the two models and establish an integrated modeling method to simultaneously achieve accuracy and efficiency. Section 2 presents the TTM formulation, simulation of laser-material interaction on the laser-pulse timescale, and the determination and experimental validation of the fluence threshold for non-thermal ablation. In Section 3, we extend the TTM to three spatial dimensions to simulate heat accumulation during non-thermal, ablation-based scanning laser surface processing and evaluate the time-efficiency of the simulation. We also describe the formulation of a classical heat accumulation model and compare its prediction of heat accumulation and its time efficiency to the TTM. In Section 4, for the first time, we devise an integrated modeling method using the TTM and heat accumulation model to achieve simultaneous accuracy and time-efficiency in predicting heat accumulation and to extend the simulation to the general-parameter case. We show that the integrated modeling method effectively predicts heat accumulation with simulation times on the order of a minute per pulse, establishing a method to determine, control, and optimize laser parameters for application-specific, multi-pulse femtosecond laser processing.

2. Predicting and validating non-thermal femtosecond laser ablation

2.1 Two-temperature model prediction of the silicon ablation threshold

To achieve our first goal of predicting initial laser parameters enabling non-thermal ablation, a two-temperature model (TTM) simulating femtosecond laser-silicon interaction was constructed based upon the methods of Chen et al. and Van Driel [6,12], using MATLAB (MATLAB R2017a). The TTM simulates three major phenomena to describe ultrafast laser-material interaction on the pulse timescale, respectively represented by Eqs. (1)-(3): (1) free-carrier electron generation, (2) carrier temperature change, and (3) lattice heating. The temperature, intensity, and wavelength dependencies of each phenomenon are considered.

Equation (1) calculates the number density of free-carrier electrons, $N_c$, which is influenced by carrier generation ($G_c$), carrier depletion ($D_c$), and carrier flow. Non-thermal ablation occurs when $N_c$ reaches the critical density for silicon breakdown (6.9 × 10²⁰ cm⁻³ for 1030-nm wavelength laser pulses [17]) while the bulk lattice temperature remains below the melting point (1687 K [12]). Equation (2) is used to calculate the carrier system temperature, $T_c$, which is raised by absorption of laser energy, $S$, and respectively decreased by thermal energy coupling to the material lattice (with temperature, $T_l$), ambipolar diffusion of electrons, change in the carrier-system kinetic energy, and change in the material bandgap energy ($E_g=$1.12 eV for silicon at room temperature) [12]. Equation (3) calculates heating of the material lattice due to carrier energy coupling and to bulk heat diffusion. Further descriptions of Eqs. (1)-(3) are provided in the Appendix and detailed studies of TTM influence parameters are provided in [6, 12, 17].

We use the TTM to simulate the impact of laser fluence on the achieved maximum free-carrier number density to determine the threshold for achieving non-thermal ablation. Figure 1(a) shows that the maximum free-carrier number density increases with increasing laser fluence. The determined fluence threshold for achieving ablation is ~0.45 J/cm² at which the carrier density is equal to the critical value of 6.9 × 10²⁰ cm⁻³. For each simulation, the initial carrier number density was set to $N_c$ = 10¹² cm⁻³ and both initial lattice and carrier temperatures were set to 300 K [6,12]. Spatial and temporal sampling are simulation-dependent parameters, generally on the respective sub-micron and femtosecond orders, and were selected to avoid steep temperature and carrier density gradients and to ensure complete interactions.

 

Fig. 1 (a) TTM-simulated maximum electron number density versus laser fluence for laser-silicon interaction (300-fs pulse width, 1030-nm wavelength). The dashed line marks the critical density above which ablation occurs (6.9 × 10²⁰ cm⁻³ [17]). (b) Predicted carrier number density, carrier temperature, and lattice temperature with respective maxima of 9.3 × 10²⁰ cm⁻³, 2.1 × 10⁴ K, and 847 K for a TTM simulation at the laser fluence of 0.53 J/cm². The vertical line marks the time of peak pulse intensity at 0 ps.

Download Full Size | PDF

Figure 1(b) shows an example of the TTM prediction of free-carrier electron density, carrier temperature, and lattice temperature for laser-silicon interaction using a fluence above the determined ablation threshold. The carrier energy and temperature rise as free carriers are generated. The peak number density and the peak temperature of the free-carrier system are achieved 0.3 ps after the arrival of the peak pulse intensity. Over time, carrier energy is transferred to the material lattice via collisional interaction [12], causing the lattice temperature to rise and eventually thermalize (equilibrate) with the decaying carrier temperature. The lattice temperature rises to 95% of the equilibrium temperature (847 K) just two picoseconds after the arrival of the peak pulse intensity. For this set of laser parameters, non-thermal ablation is predicted to occur because the peak electron number density exceeds the critical electron density for ablation while the lattice temperature remains below the silicon melting point (1687 K [12]). Therefore, the TTM is capable of predicting both the onset of ablation and the temperature rise in the material lattice to assess and control laser fluence for achieving non-thermal ablation.

2.2 Experimental validation of TTM predictions

The TTM-predicted single-pulse ablation threshold of silicon was experimentally validated via study of ablation sensitivity to laser fluence. An ultrafast laser beam (Satsuma HP3, Amplitude Systèmes) producing ~300-fs laser pulses at 1030 nm was focused to a 75-μm ($1/e^2$ diameter) spot. The laser fluence was varied by changing the pulse energy. The area of each ablation-generated crater was measured using an optical profiler (Zygo NewView 600). Figure 2(a) shows that the ablated crater area increases logarithmically with increasing laser fluence. Figure 2(b) displays an example of the surface profile of an ablated crater.

 

Fig. 2 (a) Ablated crater area on silicon versus laser pulse fluence. The ablation threshold was determined to be 0.43 J/cm², corresponding to the x-intercept of the regression. (b) Surface height map of single-shot laser ablation of silicon at 0.65 J/cm² fluence with depth and width on the orders of 10 nm and 35 µm, respectively.

Download Full Size | PDF

The ablation threshold fluence for silicon processing at 1030 nm was determined to be 0.43 J/cm² by fitting Eq. (4) to the experimental crater area data [17,18].

In Eq. (4), $w_o$ is the $1/e^2$ radius of the focused Gaussian beam, $F_{laser}$ is the peak laser fluence, and $F_{threshold}$ is the ablation threshold fluence. The determined ablation threshold is consistent with the TTM prediction of 0.45 J/cm² and agrees with the numerical and experimental predictions reported by Thorstensen and Foss [17], validating the capability of the TTM to predict ablation. This result indirectly validates the TTM’s ability to accurately predict lattice temperature because the number density of carriers, which determines the onset of ablation, also influences the carrier system temperature which controls lattice heating.

3. Predicting heat accumulation during multi-pulse femtosecond laser processing

3.1 Extension of the TTM to predict heat accumulation

The second part of our two-fold modeling goal is to predict the impact of laser parameters on heat accumulation and determine parameter combinations to avoid onset of detrimental melting, oxidation, and material pileup during non-thermal ablation processes. To predict heat accumulation, we have extended the TTM to three spatial dimensions to simulate multi-pulse laser-material interaction for area processing – advancing the applicability of existing one-dimensional two-temperature models [6,12]. The TTM offers versatile simulation of heat accumulation for different laser-parameter combinations because it considers the interaction effects of laser wavelength, pulse width, and fluence and can be adjusted to simulate different materials.

The three-dimensional TTM (3D-TTM) calculates laser interaction from irradiation through thermalization of carrier and lattice system temperatures. The Gaussian-distributed laser pulse intensity is divided into spatially discrete values which serve as inputs to computationally parallel one-dimensional TTM simulations. Following interaction, in the time between laser pulses, heat diffusion is calculated throughout the entire material bulk. Immediately prior to the next pulse-incidence, the laser beam is spatially displaced relative to its previous location according to the laser repetition rate and scanning speed. This procedure is repeated to simulate each subsequent laser pulse.

The 3D-TTM was used to predict the maximum surface temperature of silicon during multi-pulse processing to determine laser parameters to control heat accumulation. Simulation parameters were selected to ensure the onset of non-thermal ablation, enable sufficient time for heat diffusion between laser pulses, and take advantage of the maximum speed of a commercial optical laser scanner. Figure 3 shows that the surface temperature rise due to pulse incidence remains below the silicon melting point (1687 K [12]) and that the post-diffusion surface temperatures (settling temperatures) remain below the temperature at which silicon oxidation becomes significant (973 K [19]). These conditions exemplify a thermally-controlled scenario for precision, non-thermal ablation-based surface processing, as the laser parameters enable the onset of ablation and prevent the onset of undesirable thermal effects, i.e. melting and significant oxidation.

 

Fig. 3 3D-TTM prediction of maximum surface temperature versus time for multi-pulse laser processing of silicon using a 350-fs pulse width, a 515-nm wavelength, a 0.52 J/cm² peak laser fluence (~0.2 J/cm² ablation threshold for 515-nm silicon processing [17]), a 500-kHz repetition rate, and a 4 m/s scanning speed. The two horizontal lines mark the melting point (T_m = 1687 K) and the temperature at which oxidation becomes significant (T_o = 973 K) for silicon.

Download Full Size | PDF

The time required to complete each 3D-TTM simulation depends on its spatial and temporal discretization, the number of voxels comprising the simulated material lattice, and the number of timesteps required to achieve thermalization and complete heat diffusion. For example, the five-pulse simulation shown in Fig. 3 takes approximately 3.5 hours (~0.7 hours per pulse) to complete when using 2-µm sampling of the input laser intensity, sub-micron/fs-order spatial/temporal sampling for the parallel one-dimensional TTM simulations, and using cubic-micron-order voxels/nanosecond timesteps for heat diffusion calculations. A more time-efficient method is required to simulate heat accumulation during multi-pulse laser-material interaction for area processing because the 3D-TTM simulation time is unsuitable for simulating hundreds to thousands of laser pulses.

3.2 Time-efficient prediction of heat accumulation via classical modeling and comparison to the 3D-TTM

Figure 1(b) in Section 2.1 demonstrated that the free-carrier electron system achieves its peak temperature and imparts heat to/thermalizes with the lattice system on the picosecond timescale. This timescale is orders of magnitude shorter than the microsecond timescale for bulk heat diffusion between incidence of laser pulses [7]. For example, when processing with a 500-kHz repetition rate, the time allowed for heat diffusion between pulses is 2 µs, six orders of magnitude larger than the picosecond timescale for electron-lattice thermalization. Therefore, the temperature rise in the bulk material caused by a femtosecond laser pulse can be considered instantaneous with respect to the calculation of heat accumulation.

Equation (3) of the TTM, which describes lattice heating, demonstrates that the exchange of thermal energy between carrier and lattice systems becomes zero once thermalization occurs. Equation (3) then becomes Eq. (5), the classical heat conduction equation describing heat diffusion [20]:

In Eq. (5), $\rho$ is the silicon lattice density, $C_l$ is the specific heat capacity, and ${\kappa }_l$ is the thermal conductivity. Thus, the calculation of pulse-induced lattice temperature rise and bulk heat diffusion can be performed in a separable fashion using a time-efficient, classical heat accumulation model [9].

Equations (5) – (7) comprise the classical heat accumulation model. For each laser pulse, Eq. (6) is used to calculate the lattice surface energy induced by the irradiation process once thermalization of the carrier and lattice systems has been achieved.

Equation (6) describes the lattice surface energy, $E_l$, resulting from incidence of a spatial-Gaussian laser pulse with waist, $w_o$, energy, $E_p$, and peak located at $\left(x_c,y_c\right)$. Here, $\text{Δ}x$ and $\text{Δ}y$ represent the lateral sampling of each surface voxel, each with a physical dimension of 2 µm, and $E_i$ is the initial surface energy prior to irradiation. Instantaneous calculation of energy absorption is enabled by, $A$ – a coefficient describing the amount of pulse energy imparted to the material once thermalization of the carrier and lattice systems has been achieved. The classical heat accumulation model is not capable of simulating energy absorption on the interaction timescale, so this coefficient is usually experimentally determined and, therefore, case-specific to interaction parameters, e.g., wavelength, pulse width, and fluence [9]. The coefficient does not separately consider carrier generation, so the classical heat accumulation model is unable to predict non-thermal ablation.

The corresponding temperature of the material lattice following laser interaction is calculated using Eq. (7) [9].

Here, $V$ is the physical volume of a simulated voxel, and $C_l\left(T_l\right)$ is the lattice-temperature-dependent specific heat capacity. Following the determination of temperature rise, heat diffusion is calculated throughout the material bulk via Eq. (5). The classical heat accumulation model iteratively calculates energy rise, induced temperature, and heat diffusion for each laser pulse, adjusting the spatial coordinates of the beam prior to each pulse-incidence to accommodate laser scanning.

Figure 4 shows that the classical heat accumulation model prediction of surface temperature evolution is consistent with the 3D-TTM prediction for the same laser processing parameters as used in Fig. 3. The settling temperatures predicted by the classical heat accumulation model differ by only 20 K from the 3D-TTM-predicted values, but the laser-induced surface temperature rises are approximately 300 K higher due to the use of the interaction-parameter-specific absorption coefficient, $A$ = 0.8 [21]. This experimentally-determined value was the best-available value from a survey of relevant literature. It is the average value over a range of fluences up to 1 J/cm² and was determined using multi-pulse irradiation with an 800-nm wavelength laser so it inherently considers in situ changes in optical properties [21]. The coefficient can be implemented to simulate 515-nm laser interaction because single-photon absorption (see Appendix) dominates the absorption mechanism for both wavelengths – i.e., the corresponding photon energies at these two wavelengths (2.41 eV at 515-nm, 1.55 eV at 800-nm) exceed the silicon bandgap energy (1.12 eV). However, use of this coefficient limits the heat accumulation modeling to specific cases for particular combinations of laser wavelength, pulse width, and fluence.

 

Fig. 4 Surface temperature versus time/arrival of different pulses predicted by the classical heat accumulation model (dotted line) and 3D-TTM (solid line; from Fig. 3). Models consistently predict settling temperatures for the same set of simulation parameters (515-nm wavelength, 350-fs pulse width, 0.52 J/cm² peak fluence, 500 kHz repetition rate, and 4 m/s laser scanning speed).

Download Full Size | PDF

As standalone models, neither the 3D-TTM nor the heat accumulation model is ideal for simulating multi-pulse interaction: the TTM is unsuitable for extension to area laser processing simulations due to lengthy processing times (e.g. 3.5 hours for the 3D-TTM vs. less than 4 minutes for the classical heat accumulation model) and the heat accumulation model is case-specific to laser processing conditions, requiring a priori determination of the input absorption coefficient. To extend the time-efficient heat accumulation model to the general case, an offline method is required to determine the absorption coefficient across all potential combinations of laser processing parameters and for various materials since it is time-consuming and impractical to determine experimentally.

4. Integrating models to achieve accurate, versatile, and time-efficient prediction of heat accumulation

To overcome the tradeoff between model versatility and time-efficiency, we devise an integrated modeling method using the 3D-TTM and classical heat accumulation model to accurately predict heat accumulation with reduced simulation times: (1) use the TTM to numerically determine an operational absorption coefficient describing the lattice energy after thermalization and (2) feed this coefficient to the heat accumulation model in place of the experimentally-determined absorption coefficient, $A$, to improve its accuracy. This methodology eliminates the need for time-consuming, iterative experimental calibration of the absorption coefficient. By using the TTM to determine $A$, the integrated modeling method maintains the meticulous nature of a TTM simulation since all TTM-simulated physical phenomena and in situ changes in optical/material properties are inherently incorporated into the absorption coefficient determination. The integrated modeling method enables simulation times comparable to those of the classical heat accumulation model since the multi-pulse TTM simulations for coefficient determination can be performed offline and fed back, as needed, as inputs to the heat accumulation model.

The TTM can be used to predict an operational absorption coefficient via Eq. (8), which determines the ratio of the energy, $E_a$, absorbed by the material lattice once thermalization is achieved and the energy of a single laser pulse, $E_p$.

Absorbed energy,$E_a$, is calculated as the difference between the lattice energy in its thermalized and pre-irradiation states. The energy in each of these states can be calculated by rearranging Eq. (7) as $E_l=\rho \cdot V\cdot T_l\cdot C\left(T_l\right)$.

The average absorption coefficient for the five-pulse 3D-TTM simulation in Section 3.1 was found to be $\overline{\eta }$ = 0.59. This coefficient was fed to the classical heat accumulation model in place of $A$ in Eq. (6) to re-calculate the heat accumulation using the same laser parameters as in Fig. 4.

Figure 5 shows that the simulation of heat accumulation using the improved, integrated modeling method consistently predicts surface temperatures with the 3D-TTM (differing by less than 65 K) whereas the original heat accumulation model over-predicted temperature rises by up to 300 K. The simulation durations of the integrated modeling method are equivalent to the original heat accumulation model (< 4 min. for five pulses). The surface temperature agreement and the maintained time-efficiency show that the integrated modeling method enables accurate predictions of heat accumulation and is suitable for simulating multi-pulse laser-silicon interaction. The integrated modeling method is capable of predicting optimal laser parameters to control heat accumulation, while only the TTM is capable of predicting laser fluences for achieving non-thermal ablation.

 

Fig. 5 Comparison of the integrated modeling method (dashed line) and 3D-TTM (solid line) prediction of heat accumulation using the same laser parameters as in Fig. 4 (0.52-J/cm² laser fluence, 500-kHz repetition rate, 4-m/s scanning speed, 350-fs pulse width, and 515-nm wavelength). The models consistently predict the rise and settling of surface temperatures.

Download Full Size | PDF

5. Conclusion

Efficient and accurate models are required to determine optimal laser parameters for femtosecond laser processing applications relying on fundamentally different physical mechanisms: non-thermal ablation via solid-vapor phase transitions and thermal processing via solid-liquid-vapor transitions. We accurately predict non-thermal ablation and heat accumulation during multi-pulse femtosecond laser-material interaction using a three-dimensional TTM. We, for the first time to our knowledge, accurately and efficiently predict heat accumulation via an integrated modeling method which links the TTM and a classical heat accumulation model. Only the TTM predicts the laser fluences required to achieve non-thermal ablation, which was experimentally validated. The TTM was extended to simulate heat accumulation during multi-pulse irradiation, but at the cost of lengthy simulation times. The classical heat accumulation model enabled time-efficient prediction of heat accumulation induced by multi-pulse laser interaction, but it is limited to specific cases of laser incidence due to the required input of an experimentally-determined, interaction-specific absorption coefficient. An integrated modeling method for predicting heat accumulation was devised to overcome the tradeoff between model versatility and time-efficiency by using the TTM to predict and feed an absorption coefficient to the heat accumulation model and extend it to the general case. High accuracy and efficiency were simultaneously achieved via the integrated modeling method, demonstrating a path for model reduction and time-efficient simulation of multi-pulse ultrafast laser processing.