Method for simulating and predicting ablation of silicon carbide by femtosecond laser based on dual-temperature model

Abstract

The invention discloses a method for simulating and predicting ablation of silicon carbide by femtosecond laser based on a dual-temperature model, and relates to the technical field of femtosecond laser processing. According to the method, the local carrier density of a silicon carbide material is corrected based on a Fokker-Plank equation; an improved dual-temperature model suitable for interaction of femtosecond laser and silicon carbide is constructed by combining a heat source model of the femtosecond laser; and then, after spatial discretization and time discretization in the depth direction are carried out on the material spatial three-dimensional model, a finite difference method is adopted to calculate the improved dual-temperature model to obtain a physical quantity evolution process in each grid unit, so that simulation prediction is carried out on the interaction process of femtosecond laser and silicon carbide under different laser parameters. The simulation prediction result can reflect the removal mechanism of thermally modified and non-thermally modified materials in the process of processing the silicon carbide by the femtosecond laser, and is favorable for guiding the practical process of processing the silicon carbide by the femtosecond laser, so that the processing quality and efficiency of the silicon carbide by the femtosecond laser can be improved.

Description

Technical Field

[0001] The present application relates to the technical field of femtosecond laser processing, and in particular to a simulation prediction method for femtosecond laser ablation of silicon carbide based on a two-temperature model. Background Art

[0002] Silicon carbide is a wide-bandgap semiconductor material and an ideal material for fabricating high-temperature, high-frequency, high-power, and high-voltage devices, and is widely used in new energy vehicles and solar cells. Due to the extremely high hardness of silicon carbide, the tool wear is severe during the processing, so the processing difficulty is great, which greatly limits the application of silicon carbide.

[0003] Femtosecond laser, due to its ultra-high energy density and extremely short pulse width, can effectively avoid the plasma shielding phenomenon and obtain very good surface quality. Therefore, it is often used as an excellent process for processing silicon carbide. However, the interaction process between femtosecond laser and silicon carbide is extremely complex, and there is no conclusive theory on the removal mechanism of femtosecond laser ablation of silicon carbide. Therefore, at present, in order to achieve better processing effects, it is usually necessary to rely on manual experience combined with a large number of experimental trials to adjust the processing parameters of femtosecond laser to achieve the required ablation effect, with low efficiency and often unsatisfactory results. Summary of the Invention

[0004] In view of the above problems and technical requirements, the present application proposes a simulation prediction method for femtosecond laser ablation of silicon carbide based on a two-temperature model. The technical solution of the present application is as follows:

[0005] A simulation prediction method for femtosecond laser ablation of silicon carbide based on a two-temperature model, the simulation prediction method for femtosecond laser ablation of silicon carbide includes:

[0006] Construct an ablation model, the ablation model includes a three-dimensional material space model of silicon carbide material and a femtosecond laser, and the laser focus of the femtosecond laser is located at the surface of the silicon carbide material;

[0007] Construct a heat source model characterizing the relationship between the laser intensity I of the femtosecond laser, laser parameters, time t, and the depth z below the surface of the silicon carbide material. The laser intensity I conforms to a Gaussian distribution with respect to time t and depth z respectively;

[0008] Introduce a flux multiplication model and correct the local carrier density n of the silicon carbide material based on the photoionization equation derived from the Fokker - Plank equation e and construct a two-temperature model for describing the coupled variation relationship between the electron temperature T e and the lattice temperature T l when the femtosecond laser interacts with the silicon carbide material;

[0009] The three-dimensional model of the material space is spatially discretized along the depth direction into multiple grid cells, and the initial conditions and boundary conditions of the three-dimensional model of the material space are set. The simulation prediction process is discretized into multiple time steps, and the two-temperature model is calculated using the finite difference method at each time step to obtain the electron temperature T e in each grid cell, the lattice temperature T l and the global carrier density N e for ablation simulation prediction.

[0010] A further technical solution thereof is that in the heat source model, the laser parameters of the femtosecond laser include the laser energy density F and the pulse width τ p , and the laser intensity I is also related to the surface reflectivity R of the silicon carbide material, and the surface reflectivity R of the silicon carbide material is dynamically corrected using the Drude model. Then, the surface reflectivity in each grid cell is also calculated at each time step.

[0011] A further technical solution thereof is that the constructed heat source model is:

[0012]

[0013] where α1 is the single-photon absorption coefficient and θ is the free electron absorption coefficient.

[0014] A further technical solution thereof is that the surface reflectivity of the silicon carbide material where the refractive index the extinction coefficient Re represents the real part, Im represents the imaginary part, and:

[0015]

[0016] where ε sic is the initial dielectric function of the silicon carbide material, the plasma frequency the laser frequency e is the electron charge, m e is the electron mass, ε0 is the vacuum permittivity, c v is the vacuum light speed, λ is the laser wavelength, τ D is the reciprocal of the collision frequency, and i represents the imaginary number.

[0017] A further technical solution thereof is that the constructed two-temperature model includes:

[0018]

[0019] where C e is the electron heat capacity, C l is the lattice heat capacity, K e is the electron thermal conductivity, K lis the lattice thermal conductivity; g is the energy coupling rate between carriers and the lattice, S(z,t) is the laser source term and:

[0020] S(z,t) = α1I + α2I 2 + α3I 3 + N e θI

[0021] where α1 is the single - photon absorption coefficient, α2 is the two - photon absorption coefficient, α3 is the three - photon absorption coefficient, θ is the free - electron absorption coefficient, and the global carrier density N e varies with the local carrier density n e changing.

[0022] A further technical solution is that the formula for calculating the local carrier density n of the silicon carbide material in the constructed two - temperature model is: e The calculation formula is:

[0023]

[0024] where h is Planck's constant, v is the laser oscillation frequency, β is the impact ionization coefficient, and γ is the Auger recombination coefficient.

[0025] A further technical solution is that using the finite - difference method to calculate the two - temperature model at each time step includes:

[0026] Calculating the time difference in the two - temperature model using the forward - difference format, calculating the spatial difference in the two - temperature model using the central - difference format, and converting the two - temperature model into a difference - format two - temperature model for solution:

[0027]

[0028] where T e p+1 (m) is the electron temperature of grid cell m at the (p + 1) - th time step, T e p (m) is the electron temperature of grid cell m at the p - th time step, T l p+1 (m) is the lattice temperature of grid cell m at the (p + 1) - th time step, T l p (m) is the lattice temperature of grid cell m at the p - th time step; is the local carrier density of grid cell m at the (p + 1) - th time step, is the local carrier density of grid cell m at the p - th time step, is the global carrier density of grid cell m at the (p + 1) - th time step and varies with the local carrier density changing; S p+1,m(z,t) is the laser source term at the grid cell m at the (p + 1)-th time step, I p+1,m is the laser intensity at the grid cell m at the (p + 1)-th time step; the starting value of the integer parameter p is 1;

[0029] T e p+1 T(m + 1) is the electron temperature at the grid cell m + 1 at the (p + 1)-th time step, e p+1 T(m - 1) is the electron temperature at the grid cell m - 1 at the (p + 1)-th time step; l p+1 T(m + 1) is the lattice temperature at the grid cell m + 1 at the (p + 1)-th time step, l p+1 T(m - 1) is the lattice temperature at the grid cell m - 1 at the (p + 1)-th time step; the grid cells m + 1 and m - 1 are the grid cells adjacent to the grid cell m on both sides along the depth direction of the silicon carbide material; Δt is the time interval between two adjacent time steps, and Δz is the depth interval between two adjacent grid cells along the depth direction.

[0030] A further technical solution thereof is that setting the initial conditions and boundary conditions of the three-dimensional material space model includes:

[0031] Setting the upper surface of the three-dimensional material space model as a free boundary condition, the side surface as a periodic boundary condition, and the lower surface as an absorption boundary condition; setting the electron temperature and lattice temperature of all grid cells in the three-dimensional material space model to the same initial value, and initializing the global carrier density at all grid cells to the same initial value.

[0032] A further technical solution thereof is that the method further includes:

[0033] Performing ablation simulation predictions at different laser energy densities F respectively, obtaining the temperature-time scale curves of the lattice temperature varying with time at each grid cell at different laser energy densities F, and the global carrier density-time scale curves of the global carrier density varying with time at each grid cell at different laser energy densities F; and the temperature-energy scale curves of the maximum value of the lattice temperature varying with the laser energy density F at each grid cell, and the global carrier density-energy scale curves of the maximum value of the global carrier density varying with the laser energy density F at each grid cell;

[0034] Determining the time period when the global carrier density reaches the critical carrier density N th at which non-thermal ablation occurs at the grid cell according to the global carrier density-time scale curve; and determining the thermal ablation threshold T of the lattice temperature at different laser energy densities F according to the temperature-time scale curveth the time period during which the thermal ablation process occurs at the grid cell; determining, according to the global carrier density energy scale curve, the value such that the maximum value of the global carrier density reaches the critical carrier density N th the value range of the laser energy density F during which the non-thermal ablation process occurs at the grid cell, and determining, according to the temperature energy scale curve, the value such that the maximum value of the lattice temperature reaches the thermal ablation threshold T th the value range of the laser energy density F during which the thermal ablation process occurs at the grid cell.

[0035] A further technical solution thereof is that the constructed three-dimensional material space model includes a plurality of laser absorption layers, a heat conduction layer, and a fixed boundary layer stacked in sequence from the surface of the silicon carbide material downward along the depth direction. The method further includes:

[0036] Calculating the equilibrium temperature reached by the electron temperature and the lattice temperature in each grid cell and the time taken to reach the equilibrium temperature, using the reached equilibrium temperature as the temperature load in the molecular dynamics simulation, the time taken to reach the equilibrium temperature as the time for the temperature to rise to the target temperature in the molecular dynamics simulation, and performing the molecular dynamics simulation, and using the OVITO software to visually display the dump file output by the dynamics simulation.

[0037] The beneficial technical effects of this application are:

[0038] This application discloses a femtosecond laser ablation simulation prediction method for silicon carbide based on the two-temperature model. This method considers the carrier excitation mechanism of semiconductors, introduces the flux multiplication model, and derives the photoionization equation based on the Fokker-Plank equation to calculate the change in carrier density, so as to correct the decay rate of electron density in the two-temperature model, enabling it to accurately describe the recombination mechanism of free electrons, thereby constructing an improved two-temperature model applicable to the interaction between femtosecond laser and silicon carbide. Then, the finite difference method is used to calculate and predict the global carrier density, lattice temperature, and surface reflectivity of the silicon carbide surface, so as to simulate and predict the interaction process between femtosecond laser and silicon carbide under different laser parameters. The simulation prediction results can reflect the material removal mechanism of thermal modification during the femtosecond laser processing of silicon carbide, which is beneficial to guiding the practical process of femtosecond laser processing of silicon carbide and is also beneficial to optimizing the laser parameters, thereby improving the processing quality and processing efficiency of femtosecond laser on silicon carbide.

[0039] The two-temperature model constructed in this application also considers the influence of laser radiation on the optical properties of silicon carbide, and corrects the surface reflectivity of silicon carbide by introducing the Drude model to more realistically reflect the optical response during the interaction between femtosecond laser and silicon carbide, thereby making the ablation simulation prediction more accurate.

[0040] The ablation simulation prediction method of the present application utilizes the evolution of the temperatures of electrons and lattices, the global carrier density, etc. in each grid cell at different time points, and can not only infer the distribution characteristics of the thermal ablation process and the non-thermal ablation process on the time and laser energy density scales, but also visually simulate the microscopic removal process of thermal ablation during the interaction between femtosecond laser and silicon carbide by combining molecular dynamics simulation. Description of the Drawings

[0041] Figure 1 It is a schematic flow chart of the femtosecond laser ablation simulation prediction method for silicon carbide according to an embodiment of the present application.

[0042] Figure 2 It is a schematic diagram of the ablation model constructed in the present application.

[0043] Figure 3 It is the temperature-time scale curve of the lattice temperature of the grid cell varying with time at different laser energy densities F and its partial enlarged view obtained by using the ablation simulation prediction method of the present application in an example.

[0044] Figure 4 It is Figure 3 The global carrier density-time scale curve of the global carrier density of the grid cell varying with time at different laser energy densities F and its partial enlarged view obtained by using the ablation simulation prediction method of the present application in an example.

[0045] Figure 5 It is Figure 3 and Figure 4 The temperature-energy scale curve of the maximum value of the lattice temperature of the grid cell varying with the laser energy density F, and the global carrier density-energy scale curve of the maximum value of the global carrier density of the grid cell varying with the laser energy density F obtained by using the ablation simulation prediction method of the present application in an example.

[0046] Figure 6 It is the change of the surface reflectivity of the ablation pit with time observed by using a pump-probe imaging system in the ablation experiment at different laser energy densities in an example.

[0047] Figure 7 It is Figure 6 The comparison result between the simulation curve of the relative change rate ΔR / R of the surface reflectivity R at different laser energy densities F obtained by using the ablation simulation prediction method of the present application in an example and the pump-probe experiment result.

[0048] Figure 8 It is the atomic configuration diagram at different laser energy densities F obtained by performing molecular dynamics simulation by using the ablation simulation prediction method of the present application in an example.

[0049] Figure 9 It is the surface morphology map and partial enlarged view of the ablation pit obtained by SEM observation after the ablation experiment is carried out when the laser energy density F is 0.8 J / cm 2 3.

[0050] Figure 10 It is the surface morphology map and partial enlarged view of the ablation pit obtained by SEM observation after the ablation experiment is carried out when the laser energy density F is 1.8 J / cm 2 3.

[0051] Figure 11 It is the surface morphology map and partial enlarged view of the ablation pit obtained by SEM observation after the ablation experiment is carried out when the laser energy density F is 3 J / cm 2 3. Specific implementation manners

[0052] The following further describes the specific implementation manners of the present application with reference to the accompanying drawings.

[0053] The present application discloses a femtosecond laser ablation simulation and prediction method for silicon carbide based on a two-temperature model. Please refer to Figure 1 the flowchart of

[0054] Step S1: Construct an ablation model when femtosecond laser ablates silicon carbide material. The constructed ablation model includes a three-dimensional material space model of silicon carbide material and a femtosecond laser. The laser focus of the femtosecond laser is located at the surface z = 0 of the silicon carbide material, where z represents the depth of the femtosecond laser beam penetrating into the silicon carbide material.

[0055] In actual implementation, the ablation model can be constructed using lammps software. Please refer to Figure 2 The constructed three-dimensional material space model includes a heat absorption zone, a heat conduction zone, and a fixed zone stacked in sequence from the surface of the silicon carbide material downward along the depth direction. The thickness of the heat conduction zone is much larger than that of the heat absorption zone and the fixed zone. In addition, above the surface of the three-dimensional material space model is a vacuum zone, which allows the surface of the silicon carbide material to expand freely, and clusters, molecules, etc. generated by ablation can also enter the vacuum zone. In one example, the thickness of the heat absorption zone is set to 50 nm and is divided into 10 laser absorption layers along the depth direction. The thickness of the heat conduction zone is set to 30 mm and supports heat dissipation downward. The thickness of the fixed zone is set to 20 nm. The thickness of the vacuum zone above the surface of the three-dimensional material space model is 100 nm.

[0056] Step S2, construct a heat source model that characterizes the relationship between the laser intensity I of femtosecond laser, laser parameters, time t, and the depth z below the surface of the silicon carbide material. The laser intensity I conforms to a Gaussian distribution with respect to time t and depth z respectively. Before the femtosecond laser interacts with the silicon carbide material, its original beam exhibits Gaussian distribution characteristics in both time and space. Therefore, a normal distribution function is used to simulate the time distribution of the laser intensity I. At the same time, the laser focus is located at the surface z = 0 of the silicon carbide material. Based on these conditions, a heat source model that conforms to a Gaussian distribution in both time and space can be constructed.

[0057] In this heat source model, the laser parameters of the femtosecond laser include the laser energy density F and the pulse width τ p , and the laser intensity I is also related to the surface reflectivity R of the silicon carbide material. The constructed heat source model is:

[0058]

[0059] where α1 is the single-photon absorption coefficient, θ is the free-electron absorption coefficient, and N e is the global carrier density of the silicon carbide material.

[0060] Considering that laser irradiation will affect the optical properties of the silicon carbide material, resulting in a change in the surface reflectivity R of the silicon carbide material, the surface reflectivity R of the silicon carbide material in the constructed heat source model is dynamically corrected using the Drude model, so as to more realistically reflect the optical response during the interaction between the femtosecond laser and the silicon carbide material. Then the surface reflectivity R of the silicon carbide material is expressed as:

[0061]

[0062] In the above formula (2), n represents the refractive index, k represents the extinction coefficient, and:

[0063]

[0064]

[0065] In the above formulas (3) and (4), Re represents the real part of the dielectric function, Im represents the imaginary part of the dielectric function, and

[0066]

[0067] In the above formula (5), ε sic is the initial dielectric function of the silicon carbide material, ω p represents the plasma frequency, ω represents the laser frequency, and:

[0068]

[0069]

[0070] Among them, e is the electronic charge, m e is the electron mass, ε0 is the vacuum permittivity, c v is the vacuum luminous intensity, λ is the laser wavelength, τ D is the reciprocal of the collision frequency, and i represents the imaginary number.

[0071] In step S3, introduce the flux multiplication model and correct the local carrier density n of the silicon carbide material based on the photoionization equation derived from the Fokker - Plank equation. Combine the laser intensity I described by the heat source model to construct a two - temperature model for describing the coupled variation relationship between the electron temperature T e and the lattice temperature T e when femtosecond laser interacts with silicon carbide material. l The two - temperature model is constructed as follows:

[0072] Since the sum of the photon - electron coupling time and the electron - phonon coupling time of the femtosecond laser is much less than the phonon - phonon coupling time, this causes that before the silicon carbide material is ablated, the electrons and the lattice cannot fully exchange heat, and there are significant differences in the thermal properties of electrons and the lattice. Therefore, a two - temperature model is needed to describe the temperature changes of electrons and the lattice respectively. The constructed two - temperature model is:

[0073]

[0074] Among them, C e is the electron heat capacity, C l is the lattice heat capacity, K e is the electron thermal conductivity, K l is the lattice thermal conductivity. g is the energy coupling rate between carriers and the lattice, S(z, t) is the laser source term and:

[0075] S(z, t) = α1I + α2I 2 + α3I 3 + N e θI (9)

[0076] Among them, α1 is the single - photon absorption coefficient, α2 is the two - photon absorption coefficient, α3 is the three - photon absorption coefficient, and θ is the free - electron absorption coefficient.

[0077] The N e in formula (9) evaries with the change of. Considering the carrier excitation mechanism of semiconductors, a flux multiplication model is introduced and the photoionization equation derived based on the Fokker - Plank equation is used to calculate the change of local carrier density, so as to correct the decay rate of local carrier density in the two - temperature model and enable it to accurately describe the recombination mechanism of free electrons. In one embodiment, the local carrier density n of the silicon carbide material e The calculation formula is:

[0078]

[0079] where h is the Planck constant, v is the laser oscillation frequency, β is the impact ionization coefficient, and γ is the Auger recombination coefficient. According to equation (10), the local carrier density n e is corrected, thereby correspondingly updating N in equation (9) e .

[0080] Step S4: The three - dimensional material space model is spatially discretized along the depth direction into multiple grid cells, and the initial conditions and boundary conditions of the three - dimensional material space model are set. The simulation prediction process is discretized into multiple time steps, and at each time step, the finite - difference method is used to calculate the two - temperature model to obtain the electron temperature T e in each grid cell, the lattice temperature T l and the global carrier density N e for ablation simulation prediction.

[0081] That is, physical space and time discretization are performed. The three - dimensional material space model of the silicon carbide material is divided into fine grid cells along the depth direction, and the time axis is segmented into tiny time steps. In one example, the spatial division step size Δz = 10 -9 m, and the time division step size Δt = 10 -15 s.

[0082] The initial conditions are set as follows: The electron temperature and lattice temperature of all grid cells in the three - dimensional material space model are set to the same initial value, and the global carrier density at all grid cell locations is initialized to the same initial value. In one example, the electron temperature and lattice temperature of all grid cells are set to 300K, which means that the system is initially in a normal - temperature environment and they are uniformly distributed throughout the space and time. The global carrier density at all grid cell locations is initialized to 10 24 m -3 , providing a basic value for the subsequent simulation of electron - related processes. The laser intensity I is initialized as a matrix of all 1s, and will be updated later according to the laser pulse characteristics and material absorption. The setting of these initial states constructs the starting point of the model. In subsequent iterative calculations, each physical quantity will continuously evolve under the action of the difference format and physical laws to simulate the actual physical process.

[0083] For the boundary conditions: the upper surface of the three-dimensional model of the material space is set as a free boundary condition, the side surfaces are set as periodic boundary conditions, and the lower surface is set as an absorbing boundary condition. Setting the upper surface as a free boundary condition makes the upper surface as if it is actually exposed to the outside world, and particles can freely enter and exit, which can simulate the evaporation and sputtering of atoms during ablation. Using periodic boundary conditions on the side surfaces enables particles to reappear from the opposite side when they reach the side boundaries, as if the region extends infinitely, which can avoid boundary interference and reflect the real process inside the silicon carbide material. Using an absorbing boundary condition on the lower surface can absorb the wave energy when it reaches the lower surface instead of reflecting it back, thereby reducing the boundary reflection error caused by the finite calculation region. The entire three-dimensional model of the material space selects the NVT ensemble to keep the number of particles, volume, and temperature constant, which is conducive to focusing on the relationship between laser energy and internal changes in the material, excluding the interference of temperature fluctuations, and deeply understanding the laser ablation mechanism. In another embodiment based on this, instead of directly taking the first grid cell from top to bottom in the depth direction of the three-dimensional model of the material space as the upper surface, the second grid cell from top to bottom in the depth direction is taken as the upper surface. Similarly, instead of directly taking the bottommost grid cell in the depth direction of the three-dimensional model of the material space as the lower surface, the fourth grid cell from bottom to top in the depth direction of the three-dimensional model of the material space is taken as the lower surface, which can effectively avoid the numerical instability of the edge points. Adopting a form similar to the Dirichlet boundary condition, the global carrier density, electron temperature, and lattice temperature at the boundary are set to be equal to those of adjacent internal points. This means that there is no additional inflow or outflow of electrons at the boundary, and heat conduction makes the temperature at the boundary transition continuously with the interior, and the lattice thermal state also remains continuous.

[0084] Then, according to the time step, the energy input of the laser heat source model is distributed to the corresponding grid cells, and nested loops are run for iterative calculation. The outer loop traverses the space, and the inner loop traverses the time. Inside the loop, the electron temperature T e 、lattice temperature T l 、global carrier density N e and surface reflectivity and other key physical quantities are continuously updated. The global carrier density N e is updated according to light absorption, ionization, and recombination. The electron temperature takes into account heat conduction, energy exchange, and light absorption. The lattice temperature is affected by lattice heat conduction and electron-lattice energy exchange. By continuously iterating, the changes of physical quantities with time and space are simulated.

[0085] The two-temperature models (8), (9), (10) constructed in this application contain time difference terms and spatial difference terms According to the principle of the finite difference method, the time difference in the two-temperature model is calculated using the forward difference scheme, and the spatial difference in the two-temperature model is calculated using the central difference scheme. After converting the partial differential equation of the two-temperature model into a difference scheme in the form of the following algebraic equation, it is solved:

[0086]

[0087] where, T e p+1 (m) is the electron temperature of grid cell m at the (p + 1)-th time step, T e p (m) is the electron temperature of grid cell m at the p-th time step, T l p+1 (m) is the lattice temperature of grid cell m at the (p + 1)-th time step, T l p (m) is the lattice temperature of grid cell m at the p-th time step. is the local carrier density of grid cell m at the (p + 1)-th time step, is the local carrier density of grid cell m at the p-th time step, is the global carrier density of grid cell m at the (p + 1)-th time step and changes with the local carrier density . S p+1,m (z, t) is the laser source term of grid cell m at the (p + 1)-th time step, I p+1,m is the laser intensity of grid cell m at the (p + 1)-th time step and is calculated by substituting the depth z at grid cell m and the time t at the (p + 1)-th time step into Equation (1). The starting value of the integer parameter p is 1.

[0088] T e p+1 (m + 1) is the electron temperature of grid cell m + 1 at the (p + 1)-th time step, T e p+1 (m - 1) is the electron temperature of grid cell m - 1 at the (p + 1)-th time step. T l p+1 (m + 1) is the lattice temperature of grid cell m + 1 at the (p + 1)-th time step, T l p+1 (m - 1) is the lattice temperature of grid cell m - 1 at the (p + 1)-th time step. Grid cell m + 1 and grid cell m - 1 are the grid cells adjacent to both sides of grid cell m along the depth direction of the silicon carbide material. Δt is the time interval between two adjacent time steps, and Δz is the depth interval between two adjacent grid cells along the depth direction.

[0089] Thus, through nested loops, physical quantities can be updated at each space-time point according to the difference format, thereby discretizing the continuous model and iteratively obtaining numerical solutions at different space-time points. Then, ablation simulations and predictions are performed separately under different laser parameters of femtosecond lasers. For femtosecond lasers, the pulse width τ p can be considered fixed and has little impact on the processing process. Therefore, the laser energy density F is mainly considered. That is, ablation simulations and predictions are performed separately under different laser energy densities F of femtosecond lasers, so as to obtain the simulated numerical solutions of each grid cell under different laser energy densities F. Then, after sorting, the temperature-time scale curves of the lattice temperature of each grid cell changing with time under different laser energy densities F can be obtained, as well as the global carrier density-time scale curves of the global carrier density of the grid cell changing with time under different laser energy densities F, and the temperature-energy scale curves of the maximum value of the lattice temperature of the grid cell changing with the laser energy density F, and the global carrier density-energy scale curves of the maximum value of the global carrier density of the grid cell changing with the laser energy density F.

[0090] In one example, for the grid cell with the highest lattice temperature and the highest global carrier density, through the subplot function, the temperature-time scale curves of the lattice temperature of this grid cell at laser energy densities F of 0.6 J / cm 2 , 0.65 J / cm 2 , 0.7 J / cm 2 changing with time are as shown in Figure 3 (a) in, and the local enlarged view within the green box in (a) is as shown in Figure 3 (b) in. The global carrier density-time scale curves of the global carrier density of this grid cell at laser energy densities F of 0.6 J / cm 2 , 0.65 J / cm 2 , 0.7 J / cm 2 changing with time are as shown in Figure 4 (a) in, and the local enlarged view within the green box in (a) is as shown in Figure 4 (b) in. The temperature-energy scale curves of the maximum value of the lattice temperature of the grid cell changing with the laser energy density F, and the global carrier density-energy scale curves of the maximum value of the global carrier density of the grid cell changing with the laser energy density F are summarized as shown in Figure 5 .

[0091] Then, by analyzing the temperature-time scale curve and the global carrier density-time scale curve, the time scale characteristics of the thermal ablation process and the non-thermal ablation process occurring at the grid cell are obtained, including: determining the time when the global carrier density reaches the critical carrier density N at different laser energy densities F according to the global carrier density-time scale curve th and the time period during which the non-thermal ablation process occurs at the grid cell; and determining the time when the lattice temperature reaches the thermal ablation threshold T at different laser energy densities F according to the temperature-time scale curve th and the time period during which the thermal ablation process occurs at the grid cell. Among them, the thermal ablation threshold T th is 3400K, and the critical carrier density N th is 0.71×10 27 / m 3 . It can be seen from Figure 3 that the higher the laser energy density F, the faster the lattice temperature rises and the faster the thermal ablation occurs. When the laser energy density F is 0.6J / cm 2 , the thermal ablation continuously occurs in the time period starting from about 1.9 ps and later. When the laser energy density F is 0.65J / cm 2 , the thermal ablation continuously occurs in the time period starting from about 1.1 ps and later. When the laser energy density F is 0.7J / cm 2 , the thermal ablation continuously occurs in the time period starting from about 0.85 ps and later. It can be seen from Figure 4 that the higher the laser energy density F, the faster the global carrier density increases and the faster the non-thermal ablation occurs. When the laser energy density F is 0.6J / cm 2 , the non-thermal ablation continuously occurs in the time period starting from about 0.197 ps and later. When the laser energy density F is 0.65J / cm 2 , the non-thermal ablation continuously occurs in the time period starting from about 0.186 ps and later. When the laser energy density F is 0.7J / cm 2 , the non-thermal ablation continuously occurs in the time period starting from about 0.178 ps and later. Combining the results of Figure 3 and Figure 4 , it can be seen that the thermal ablation process generally starts from the sub-picosecond to picosecond level and persists. The non-thermal ablation process generally starts from about 100 femtoseconds and persists. Therefore, in terms of the time scale, the non-thermal ablation process always occurs prior to the thermal ablation process, and during the entire thermal ablation process, the non-thermal ablation process always exists.

[0092] Analyzing the temperature energy scale curve and the global carrier density energy scale curve yields the energy scale characteristics of the thermal ablation process and the non-thermal ablation process occurring at this grid cell, including: determining, according to the global carrier density energy scale curve, the energy scale at which the maximum value of the global carrier density reaches the critical carrier density N th and the value range of the laser energy density F at which the non-thermal ablation process occurs at the grid cell, and determining, according to the temperature energy scale curve, the energy scale at which the maximum value of the lattice temperature reaches the thermal ablation threshold T th and the value range of the laser energy density F at which the thermal ablation process occurs at the grid cell. From Figure 5 the two energy scale curves, it can be seen that when the laser energy density F is less than 0.43 J / cm 2 , the global carrier density at this grid cell does not reach the critical carrier density N th , and the lattice temperature does not reach the thermal ablation threshold T th , indicating that neither thermal ablation nor non-thermal ablation occurs when the laser energy density F is less than 0.43 J / cm 2 . When the laser energy density F increases to 0.43 J / cm 2 , the global carrier density at this grid cell reaches the critical carrier density N th and non-thermal ablation occurs, but the lattice temperature still does not reach the thermal ablation threshold T th , that is, non-thermal ablation occurs at this grid cell without thermal ablation. When the laser energy density F continues to increase to 0.6 J / cm 2 and higher, the global carrier density at this grid cell reaches the critical carrier density N th and non-thermal ablation occurs, and the lattice temperature reaches the thermal ablation threshold T th and thermal ablation occurs. From Figure 5 it can be seen that, in terms of energy scale, the laser energy density F for generating the non-thermal ablation process is lower than the laser energy density F for generating the ablation process, resulting in only non-thermal ablation occurring without thermal ablation when the laser energy density F is in the range of 0.43 J / cm 2 to 0.6 J / cm 2 , thus bringing about an ablation window period of 0.17 J / cm 2 , and the existence of this ablation window period provides more possibilities for precision machining.

[0093] In addition, by organizing the predicted surface reflectivity at the grid cell, the surface reflectivity time scale curve of the grid cell varying with time under different laser energy densities F can be obtained. According to the Drude model, a higher surface reflectivity represents a higher global carrier density, while a lower surface reflectivity represents a superheated liquid phase. Therefore, the surface reflectivity time scale curve can also infer the occurrence of the non-thermal ablation process.

[0094] In one example, the typical value of the laser energy density F is selected as 0.43 J / cm 2 , 0.6 J / cm 2 , 1.85 J / cm 2 , 2 J / cm 2 , 2.3 J / cm 2 And ablation experiments are respectively carried out on silicon carbide material samples, and the pump-probe imaging technology is used to observe the change of the surface reflectivity of the ablation pit with time. The pump-probe imaging system is divided into two beams of light, namely the pump beam and the probe beam. During the experiment, the pump beam is focused on the ablation pit of the silicon carbide material sample to cause changes, and the probe beam is irradiated at the same position after being controlled by a delay line at different time delays. The reflectivity of the probe beam is accurately measured by a detector, and multiple measurements are averaged to ensure the accuracy of the data, thereby obtaining the experimental results of the surface reflectivity. The experimental results of the surface reflectivity of the ablation pit observed under different laser energy densities F change with time as Figure 6 shown. From Figure 6 it can be seen that when the laser energy density F is 0.43 J / cm 2 , an oval bright spot that gradually becomes larger and brighter appears on the silicon carbide surface. Since the pump-probe light is incident at a 45° angle, this oval bright spot is caused. The continuous increase in the surface reflectivity represents a relatively high global carrier density. Therefore, it is speculated that Coulomb explosion may be one of the reasons for the increase in the surface reflectivity, and the influence of material densification may also be a reason. However, the removal of materials by Coulomb explosion is usually at the nanoscale, so there is no obvious removal phenomenon on the silicon carbide surface at this time. When the laser energy density F is 0.6 J / cm 2 , only non-thermal ablation occurs before 10 ps, and there are no obvious material removal traces on the material surface. After that, an annular dark area appears around the ablation pit and expands with time. At this time, the decrease in the surface reflectivity is related to the overheated state of the silicon carbide surface, but there are always some alternating bright and dark rings in the ablation pit. This is actually the Newton's rings caused by spot interference. The shedding of the thin and smooth shell on the material surface participates in the generation of Newton's rings. Therefore, the appearance of Newton's rings represents the occurrence of non-thermal ablation phenomenon of photo-induced mechanical fragmentation. When the laser energy density F is large enough (such as 1.85 J / cm 2 , 2 J / cm 2 , 2.3 J / cm 2At this time, non-thermal ablation appears earlier, and the higher the laser energy density F, the higher the surface reflectivity and the corresponding global carrier density at 1 ps. The decrease in surface reflectivity first appears between 1 - 10 ps. On the one hand, it is due to the decrease in global carrier density, and on the other hand, it may be related to the light absorption of the overheated liquid phase. In the range of dozens of picoseconds to hundreds of picoseconds, Newton's rings appear at the edge of the ablation pit, indicating the presence of non-thermal ablation at the edge at this time. It can be seen from this example that the non-thermal ablation process and the thermal ablation process can also be predicted by predicting the surface reflectivity, and the prediction results obtained by using the simulation prediction method of the present application are also consistent with the observation results of the pump-probe imaging system.

[0095] Then, the experimental results of the observed surface reflectivity over time are plotted as a curve, and it is compared with the curve of the surface reflectivity R over time obtained in the present application under the corresponding laser energy density F. The pump-probe imaging system is divided into two beams of light, namely the pump beam and the probe beam. During the experiment, the pump beam is focused on the ablation pit of the silicon carbide material sample to cause changes, and the probe beam is irradiated at the same position after being controlled by a delay line at different time delays. The reflectivity of the probe beam is accurately measured by a detector, and multiple measurements are averaged to ensure the accuracy of the data. The curve of the surface reflectivity R over time is plotted and compared with the curve of the surface reflectivity R over time calculated in the present application. Only taking the cases where the laser energy density F is 1.85 J / cm 2 、2 J / cm 2 、2.3 J / cm 2 as examples, the comparison results of the simulation curve of the relative change rate ΔR / R of the surface reflectivity R calculated in the present application and the experimental curve observed by using the pump-probe imaging technology are as Figure 7 shown. It can be seen from Figure 7 that within the allowable error range, the simulation curve predicted in the present application is basically consistent with the observed experimental curve, indicating the accuracy of the simulation prediction method of the present application.

[0096] In another embodiment, combining the calculation results of the dual-temperature model of the present application with molecular dynamics simulation (MD simulation) can also accurately simulate the thermal removal process of silicon carbide materials and accurately predict the surface optoelectronic properties, including: calculating the equilibrium temperature reached by the electron temperature and the lattice temperature in each grid cell and the time taken to reach the equilibrium temperature, using the reached equilibrium temperature as the temperature load in the molecular dynamics simulation, and the time taken to reach the equilibrium temperature as the time for the temperature to rise to the target temperature in the molecular dynamics simulation and performing the molecular dynamics simulation. When simulating the material removal process, set the interatomic interaction potential of silicon carbide to accurately describe the force situation between silicon carbide atoms. At the same time, use an appropriate time step, monitor the simulation process in real time, and dynamically adjust the parameters to ensure the reliability of the simulation results. After the simulation is completed, collect the structural information and electronic state data of the surface silicon carbide atoms, and use the OVITO software to visually display the dump file output by the kinetic simulation.

[0097] In one example, the laser energy density F is 0.8 J / cm 2 When the method of the present application is used for simulation and prediction, the atomic configuration diagram is as Figure 8 shown in (a) of. In this case, the surface temperature of the silicon carbide material rises sharply to 8000 K. This high temperature causes significant thermal expansion of the silicon carbide material. At about 15 ps, cavities begin to form, and within the next few picoseconds, the cavities continue to expand, directly causing the silicon carbide material to break. When the laser energy density F is 1.8 J / cm 2 When the method of the present application is used for simulation and prediction, the atomic configuration diagram is as Figure 8 shown in (b) of. At this time, the temperature reaches 20000 K. At 15 ps, atoms or clusters have flown out from the surface, indicating that heat accumulates faster, the surface temperature quickly reaches the vaporization temperature and continues to rise. While the surface atoms vaporize, the deep-layer material becomes very loose. At 25 ps, multiple cavities appear. These cavities expand, connect with each other, and finally contribute to the spallation of the liquid phase part. As time goes by, at 50 ps, the deep-layer material has spalled into larger clusters, while the surface material is still removed in the form of droplets. When the laser energy density F is 3 J / cm 2 When the method of the present application is used for simulation and prediction, the atomic configuration diagram is as Figure 8 shown in (c) of. At this time, the laser energy density F is large enough, and the overheated liquid phase region splits into a mixture of a large number of vaporized atoms and droplets in a very short time, a phase explosion forms, and the mixture of vaporized atoms, clusters and droplets is ejected at a very high speed under the combined action of high temperature and high pressure, causing severe material removal. Therefore, the removal mode of silicon carbide at high energy density is phase explosion. It can be intuitively seen that the material removal forms are different when the laser energy density F is different.

[0098] InFigure 8 The ablation experiment was carried out under the same conditions and the surface morphology of the ablation pit after ablation was measured using a scanning tunneling microscope (SEM) and compared with the simulation prediction results of this application. The laser energy density F is 0.8 J / cm 2 The surface morphology of the ablation pit measured at a scale of 10 μm after the ablation is as follows Figure 9 As shown in (a), the surface morphology at the edge of the ablation pit at a scale of 1 μm is shown in Figure 9 As shown in (b), the surface morphology at the edge of the ablation pit at a scale of 500 nm is shown in Figure 9 As shown in (c), the surface morphology at the edge of the ablation pit at a scale of 200 nm is shown in Figure 9 The surface morphology at the center of the ablation pit at a scale of 1 μm is shown in (d). Figure 9 As shown in (e) in the figure, the surface morphology at the center of the ablation pit at a scale of 500 nm is shown in Figure 9 As shown in (f) in the figure, the surface morphology at the center of the ablation pit at a scale of 200 nm is shown in Figure 9 As shown in (g) in . Figure 9 The SEM measurement results show that when the laser energy density F is low, the surface recasts are mostly large droplets with smooth edges, which corresponds to the critical point phase separation at low energy density in molecular dynamics. Figure 8 The simulation prediction results of (a) are basically consistent.

[0099] Laser energy density F is 1.8J / cm 2 The surface morphology of the ablation pit measured at a scale of 10 μm after the ablation is as follows Figure 10 As shown in (a), the surface morphology at the edge of the ablation pit at a scale of 1 μm is shown in Figure 10 As shown in (b), the surface morphology at the edge of the ablation pit at a scale of 500 nm is shown in Figure 10 As shown in (c), the surface morphology at the edge of the ablation pit at a scale of 200 nm is shown in Figure 10 The surface morphology at the center of the ablation pit at a scale of 1 μm is shown in (d). Figure 10 As shown in (e) in the figure, the surface morphology at the center of the ablation pit at a scale of 500 nm is shown in Figure 10 As shown in (f) in the figure, the surface morphology at the center of the ablation pit at a scale of 200 nm is shown in Figure 10 As shown in (g) in . Figure 10 The SEM measurement results show that the volume of the surface recast at the medium energy density F is slightly smaller than that at the low energy density, and the edge of the recast is not clear, which corresponds to the large clusters and small droplets caused by the gasification and liquid phase spalling in the molecular dynamics simulation.Figure 8 The simulation prediction results in (b) are basically consistent.

[0100] When the laser energy density F is 3 J / cm 2 After ablation, the surface topography map of the ablation pit at a scale of 10 μm measured is as shown in Figure 11 (a) in, and the surface topography map at the edge of the ablation pit at a scale of 1 μm is as shown in Figure 11 (b) in, and the surface topography map at the edge of the ablation pit at a scale of 500 nm is as shown in Figure 11 (c) in, and the surface topography map at the edge of the ablation pit at a scale of 200 nm is as shown in Figure 11 (d) in. The surface topography map at the center of the ablation pit at a scale of 1 μm is as shown in Figure 11 (e) in, and the surface topography map at the center of the ablation pit at a scale of 500 nm is as shown in Figure 11 (f) in, and the surface topography map at the center of the ablation pit at a scale of 200 nm is as shown in Figure 11 (g) in. From Figure 11 the SEM measurement results, it can be seen that at a high laser energy density F, the volume of the surface remelted object is significantly reduced and the distribution is significantly denser, corresponding to the extremely fine droplets and clusters generated by phase explosion in the molecular dynamics simulation, which is basically consistent with Figure 8 the simulation prediction results in (c).

[0101] The above are only the preferred embodiments of the present application, and the present application is not limited to the above examples. It can be understood that other improvements and changes directly derived or associated by those skilled in the art without departing from the spirit and concept of the present application should be considered to be included within the protection scope of the present application.

Claims

1. A femtosecond laser ablation simulation prediction method for silicon carbide based on a dual-temperature model, characterized in that: The femtosecond laser ablation simulation prediction method for silicon carbide comprises: Constructing an ablation model, wherein the ablation model includes a material space three-dimensional model of the silicon carbide material and a femtosecond laser, wherein the laser focus of the femtosecond laser is located at the surface of the silicon carbide material; A heat source model is constructed to characterize the relationship between the laser intensity I of the femtosecond laser and the laser parameters, time t and the depth z below the surface of the silicon carbide material. The laser intensity I, time t and depth z respectively conform to Gaussian distribution. The local carrier density n of silicon carbide material is calculated by introducing the flux multiplication model and the photoionization equation derived from the Fokker-Plank equation. e The correction is made and combined with the laser intensity I described by the heat source model to construct a method for describing the electron temperature T when the femtosecond laser interacts with the silicon carbide material. e and the lattice temperature T l Two-temperature model of coupled variation relationship; The material space three-dimensional model is spatially discretized into multiple grid units along the depth direction, and the initial conditions and boundary conditions of the material space three-dimensional model are set. The simulation prediction process is discretized into multiple time steps, and the finite difference method is used to calculate the dual-temperature model at each time step to obtain the electron temperature T in each grid unit. e , lattice temperature T l and the global carrier density N e Perform ablation simulation prediction.

2. The femtosecond laser ablation simulation prediction method for silicon carbide according to claim 1, characterized in that: In the heat source model, the laser parameters of the femtosecond laser include the laser energy density F and the pulse width τ p , and the laser intensity I is also related to the surface reflectivity R of the silicon carbide material, and the surface reflectivity R of the silicon carbide material is dynamically corrected using the Drude model, then the surface reflectivity in each grid unit is also calculated at each time step.

3. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 2, characterized in that: The heat source model constructed is: Among them, α1 is the single-photon absorption coefficient and θ is the free electron absorption coefficient.

4. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 2, characterized in that: Surface reflectivity of silicon carbide material The refractive index Extinction coefficient Re represents the real part, Im represents the imaginary part, and: Among them, ε sic is the initial dielectric function of the SiC material, the plasma frequency Laser frequency e is the electron charge, m e is the mass of the electron, ε0 is the dielectric constant of vacuum, c v is the vacuum luminosity, λ is the laser wavelength, τ D is the reciprocal of the collision frequency, and i represents an imaginary number.

5. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 1, characterized in that: The constructed dual-temperature model includes: Among them, C e is the electron heat capacity, C l is the lattice heat capacity, K e is the electronic thermal conductivity, K l is the lattice thermal conductivity; g is the energy coupling rate between carriers and the lattice, S(z,t) is the laser source term and: S(z,t)=α1I+α2I 2 +α3I 3 +N e θI Among them, α1 is the single-photon absorption coefficient, α2 is the two-photon absorption coefficient, α3 is the three-photon absorption coefficient, θ is the free electron absorption coefficient, and the global carrier density N e As the local carrier density n e changes with the changes of.

6. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 5, characterized in that: The local carrier density n of silicon carbide material in the constructed dual-temperature model e The calculation formula is: Where h is Planck's constant, v is the laser oscillation frequency, β is the impact ionization coefficient, and γ is the Auger recombination coefficient.

7. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 6, characterized in that: Calculating the two-temperature model using the finite difference method at each time step includes: The time difference in the two-temperature model is calculated using the forward difference format, and the spatial difference in the two-temperature model is calculated using the central difference format. The two-temperature model is converted into a difference format two-temperature model and then solved: in, is the electron temperature of grid cell m at the p+1th time step, is the electron temperature of grid cell m at time step p, T l p+1 (m) is the lattice temperature of grid cell m at the p+1th time step, T l p (m) is the lattice temperature of grid cell m at the pth time step; is the local carrier density of grid cell m at the p+1th time step, is the local carrier density of grid cell m at the pth time step, is the global carrier density of grid cell m at the p+1th time step and changes with the local carrier density changes with the changes of S p+1,m (z, t) is the laser source term of grid cell m at the p+1th time step, I p+1,m is the laser intensity of grid cell m at the p+1th time step; the initial value of the integer parameter p is 1; is the electron temperature of grid cell m+1 at the p+1th time step, is the electron temperature of grid cell m-1 at the p+1th time step; T l p+1 (m+1) is the lattice temperature of grid cell m+1 at the p+1th time step, T l p+1 (m-1) is the lattice temperature of grid unit m-1 at the p+1th time step; grid unit m+1 and grid unit m-1 are the grid units adjacent to grid unit m on both sides along the depth direction of silicon carbide material; Δt is the time interval between two adjacent time steps, and Δz is the depth interval between two adjacent grid units along the depth direction.

8. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 1, characterized in that: The initial conditions and boundary conditions for setting the material space 3D model include: The upper surface of the three-dimensional material space model is set as a free boundary condition, the side surface is set as a periodic boundary condition, and the lower surface is set as an absorbing boundary condition; the electron temperature and lattice temperature of all grid cells in the three-dimensional material space model are set to the same initial value, and the global carrier density at all grid cells is initialized to the same initial value.

9. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 2, characterized in that: The method further comprises: Ablation simulation prediction is performed respectively under different laser energy densities F to obtain a temperature-time scale curve of the lattice temperature of each grid unit under different laser energy densities F changing with time, and a global carrier density time scale curve of the global carrier density of the grid unit under different laser energy densities F changing with time; and a temperature energy scale curve of the maximum value of the lattice temperature of the grid unit changing with the laser energy density F, and a global carrier density energy scale curve of the maximum value of the global carrier density of the grid unit changing with the laser energy density F; According to the global carrier density time scale curve, the global carrier density reaches the critical carrier density N at different laser energy densities F. th The time period during which the non-thermal ablation process occurs at the grid unit; and the time period during which the lattice temperature reaches the thermal ablation threshold T at different laser energy densities F according to the temperature time scale curve. th The time period during which the thermal ablation process occurs at the grid unit; the time period during which the maximum value of the global carrier density reaches the critical carrier density N is determined according to the global carrier density energy scale curve. th The range of laser energy density F at which the non-thermal ablation process occurs at the grid unit, and the maximum value of the lattice temperature reaching the thermal ablation threshold T are determined according to the temperature energy scale curve. th The value range of the laser energy density F during the thermal ablation process at the grid unit.

10. The method for simulating and predicting femtosecond laser ablation of silicon carbide according to claim 1, characterized in that: The constructed material space three-dimensional model includes a plurality of laser absorption layers, heat conduction layers and fixed boundary layers stacked in sequence from top to bottom along the depth direction from the surface of the silicon carbide material, and the method further includes: The equilibrium temperature reached by the electron temperature and the lattice temperature in each grid cell and the time taken to reach the equilibrium temperature are calculated. The equilibrium temperature reached is used as the temperature load in the molecular dynamics simulation, and the time taken to reach the equilibrium temperature is used as the time for the temperature to rise to the target temperature in the molecular dynamics simulation. Molecular dynamics simulation is performed, and the dump file output by the dynamics simulation is visualized using OVITO software.

Images