A Parallel Simulation Method for Nonlinear Interaction between Laser and Matter

Abstract

The present invention provides a parallel simulation method for the nonlinear interaction between laser and matter, comprising the following steps: introducing the plasma rate equation of the nonlinear optical effect involved in the interaction between laser and matter into the Maxwell equation as an additional differential equation to obtain the density of the excited free electron plasma; constructing a constitutive equation according to the Kerr polarization intensity, the multi-photon ionization polarization intensity, the Drude polarization intensity and the electric flux density; deforming the Maxwell equation by using the relationship between the polarization current and the polarization intensity; solving the magnetic field at time t from the electric field at time t, and solving the electric displacement vector at time t from the magnetic field at time t; numerically solving the constitutive equation by the Newton iteration method, and updating the electric field at time t from the electric displacement vector at time t; and solving the free electron plasma density rate equation by the Runge-Kutta method. In the process of solving the constitutive equation by the Newton iteration method in the present invention, matrix reconstruction is performed on the physical quantities related to the optical field.

Description

Technical Field

[0001] The present invention relates to the field of numerical simulation calculations of the interaction between laser and matter during the laser processing, and particularly relates to a parallelized simulation method for the nonlinear interaction between laser and matter. Background Art

[0002] With the continuous development of laser technology, the interaction between laser and matter has attracted people's attention. The interaction process between laser and matter is a nonlinear optical process, and these nonlinear optical processes affect the quality of laser processing. It is necessary to model and simulate it, and use the simulation results to guide the actual laser processing process.

[0003] In the current numerical simulation of the interaction between laser and matter, two common simulation methods are used. The first is to solve the nonlinear Schrödinger equation to simulate the transmission phenomenon of laser inside matter; the second is to use the finite-difference time-domain method to solve the Maxwell's equations to simulate the interaction between laser and matter with various complex structures and compositions.

[0004] When the interaction between laser and matter involves nonlinear optical effects, such as Kerr effect, multi-photon absorption, multi-photon ionization, avalanche ionization, etc., the nonlinear Schrödinger equation is generally used. Numerically solving the nonlinear Schrödinger equation has been successfully used in the simulation of filamentation and supercontinuum radiation phenomena of laser in free space, the simulation of nonlinear transmission of high-power laser in optical fiber media, and the simulation of the bulk processing process of femtosecond laser on dielectric materials such as quartz. In the field of laser processing, the nonlinear Schrödinger equation is usually used to simulate the transmission phenomenon of ultrafast laser inside homogeneous matter, and is rarely used in the simulation of the interaction between laser and matter with multiple interfaces and complex geometries.

[0005] Different from the nonlinear Schrödinger equation, the finite-difference time-domain method for solving the Maxwell's equations is one of the best numerical methods for solving electromagnetic field problems. In theory, it can simulate various complex material structures, and can accurately simulate the inhomogeneity, anisotropy, dispersion characteristics and nonlinear problems of matter. However, in fact, in the current simulation of the interaction between light and matter by using the finite-difference time-domain method (FDTD) to solve the Maxwell's equations, especially in the simulation related to laser processing, most of them do not consider the nonlinear optical effects of the interaction between laser and matter. In addition, in the first stage of the nonlinear optical process of the interaction between laser and matter described above, that is, the stages of multi-photon ionization and avalanche ionization, etc., in this stage, the laser will excite free electron plasma, and the generation of free electron plasma will then affect the light field by changing the optical properties of the material. They are a process of mutual influence and alternating occurrence over time. In the current process of using FDTD to solve the Maxwell's equations to simulate the interaction between laser and material, the alternating interaction process between light and matter is rarely considered.

[0006] In order to more precisely analyze the laser processing process and improve the laser processing quality, an improved algorithm needs to be proposed, which can comprehensively consider the nonlinear optical effects, the interaction between the optical field and the matter, and the influence of the inhomogeneous structure and the pre-existing surface structure on the interaction between the laser and the matter. Summary of the Invention

[0007] Aiming at the deficiencies in the prior art, the present invention provides a parallel simulation method for the nonlinear interaction between a laser and a matter. Multiple nonlinear optical effects involved in the interaction between the laser and the matter and the interaction effect between the optical field and the matter (specifically, the free electron plasma excited in the matter) are added to the Maxwell equations. The nonlinear optical effects are expressed through the nonlinear polarization intensity; the influence of the free electron plasma on the optical field is expressed through the Drude polarization intensity, and the influence of the optical field on the free electron plasma is expressed through the plasma rate equation. In specific operations, the curl equations of the electric field and the magnetic field in the Maxwell equations are retained and corrected, the nonlinear polarization intensity and the Drude polarization intensity are incorporated into the constitutive equation, and then the free electron plasma rate equation is introduced to construct a coupled system of equations composed of the above four differential equations. Using the anisotropic perfectly matched layer (UPML) as the absorbing boundary, the finite-difference time-domain method is used to solve the curl equations of the electric field and the magnetic field, the Newton iteration method is used to solve the constitutive equation, and the Runge-Kutta method is used to solve the free electron plasma rate equation to obtain the free electron plasma excited by the optical field. As time progresses, in each time step, the update iteration of the electric field strength E → magnetic field strength H → electric displacement vector D' → electric field strength E is realized. In particular, in the numerical simulation, during the process of using the Newton iteration method to solve the constitutive equation, matrix reconstruction of the physical quantities related to the optical field can be carried out to achieve parallel computing of multiple CPUs and multiple cores of the computer.

[0008] The present invention achieves the above technical objectives through the following technical means.

[0009] A parallel simulation method for the nonlinear interaction between a laser and a matter includes the following steps:

[0010] S01: Introduce the plasma rate equation containing the nonlinear optical effects in the interaction between the laser and the matter into the Maxwell equations as an additional differential equation to obtain the density of the excited free electron plasma;

[0011] S02: Express the Kerr effect and the multi-photon absorption nonlinear optical effect generated by the material on light as the Kerr polarization intensity P kerr (t) and the multi-photon absorption polarization intensity P MPI (t); and express the influence of the generation of the free electron plasma reflected by the Drude model on the optical properties of the material as the Drude polarization intensity PDrude (t), according to the Kerr polarization intensity P kerr (t), the multi-photon absorption polarization intensity P MPI (t), the Drude polarization intensity P Drude (t) and the electric flux density D(t) to construct the constitutive equation D′(t), and the constitutive equation can be expressed as:

[0012] D′(t) = D(t) + P kerr (t) + P MPI (t) + P Drude (t)

[0013] D(t) = ε 0 ε ∞ E(t)

[0014] Where:

[0015] P Drude (t) is obtained from the following differential equation:

[0016]

[0017] P Kerr (t) is caused by the instantaneous Kerr effect and is expressed as:

[0018]

[0019] P MPI (t) is caused by multi-photon absorption, and the expression is:

[0020]

[0021] In the formula:

[0022] D′(t) is the electric displacement vector; P Drude (t) is the Drude polarization intensity; P kerr (t) is the Kerr polarization intensity; P MPI (t) is the multi-photon absorption polarization intensity; E(t) is the electric field intensity; ε 0 is the vacuum permittivity; ε ∞ is the relative permittivity of the unexcited material; ω p is the plasma frequency; is the third-order polarizability; Γ is the electron collision frequency;

[0023] S03: Use the relationship between the polarization current and the polarization intensity to transform the Maxwell equation:

[0024] Polarization current J:

[0025] The transformed Maxwell equation is:

[0026]

[0027]

[0028] In the formula:

[0029] E(t) is the electric field strength; μ 0 is the magnetic permeability of vacuum; H(t) is the magnetic field strength; J(t) is the polarization current;

[0030] S04: According to the deformed Maxwell's equations, constitutive equations and differential equations, establish a coupled equation set from the laser parameters and material parameters, mesh the simulation area, and assume that there are i max meshes in the x direction, j max meshes in the y direction, and k max meshes in the z direction;

[0031] S05: Set the anisotropic perfectly matched layer as the absorbing boundary, and use the finite-difference time-domain method to numerically solve the coupled equation set, perform unified numerical simulation on the calculation area and the boundary area, solve the magnetic field H(t) at time t from the electric field E(t) at time t, and solve the electric displacement vector D'(t) at time t from the magnetic field H(t) at time t;

[0032] S06: Numerically solve the constitutive equation in parallel by the Newton iteration method, update the electric field E(t) at time t from the electric displacement vector D′(t) at time t, and realize parallel computing of multiple CPUs and multiple cores of the computer through matrix reconstruction of the physical quantities related to the optical field during the numerical solution process;

[0033] S07: Use the Runge-Kutta method to solve the free electron plasma density rate equation, and update the plasma density n (n+1) excited by the laser from the electric field E e , plasma frequency ω p , polarization intensity P drude (t) and P MPI (t);

[0034] With Δt as the time interval, advance forward according to the time step, and repeat steps S05 to S07 until t > the maximum time t max , and end the simulation.

[0035] Furthermore, the differential equation for obtaining the excited free electron plasma density is:

[0036]

[0037] In the formula:

[0038] n 0 is the refractive index of the unexcited material; Eg is the band gap of the material; n a is the saturated electron density; n e is the free electron plasma density excited by the laser; σ is the inverse bremsstrahlung cross section; σ (k) is the multi-electron ionization coefficient; I is the laser intensity; k is the number of laser photons simultaneously absorbed by the material; τ R is the electron-hole recombination time.

[0039] Furthermore, the electric field E(t) at time t is updated from the electric displacement vector D′(t) at time t, specifically as follows:

[0040] Determine the three components E x 、E y 、E z of the optical field vector E that depends on the three-dimensional spatial coordinates on all grids:

[0041] E x (i, j, k) → i = 1 to i max

[0042] E y (i, j, k) → j = 1 to j max

[0043] E z (i, j, k) → k = 1 to k max

[0044] where: i is the grid number in the x direction; j is the grid number in the y direction; k is the grid number in the z direction; i max is the maximum number of grids in the x direction; j max is the maximum number of grids in the y direction; k max is the maximum number of grids in the z direction;

[0045] Construct a one-dimensional column matrix A, with a total of 3N rows, where N is the number of spatial grids N = i max × j max × k max , and every three rows in the one-dimensional column matrix A correspond to the three components of the optical field of a spatial grid:

[0046]

[0047] Write the function values of the three functions [X Y Z] constructed from the constitutive equation on each spatial grid as a one-dimensional row matrix, repeat this row matrix 3 times, thereby obtaining a 3-row 3-column matrix, and then stack all the 3-row 3-column matrices on all N spatial grids in the column direction to form a 3N-row 3-column matrix B:

[0048]

[0049] Among them

[0050]

[0051] In the formula: Δt is the time grid size, and n is the number of time grids; D′ (n+1) (t) is the electric displacement vector at time t = (n + 1)Δt, and E (n) (t) is the electric field strength vector at time t = nΔt, and E (n+1) (t) is the electric field strength vector of k at time t = (n + 1)Δt, is the Drude polarization intensity vector at time t = nΔt, is the Drude polarization intensity vector of the time grid number at t = (n - 1)Δt, is the multi - photon absorption polarization intensity vector at time t = nΔt, and E x is the electric field component in the x - direction, and E y is the electric field component in the y - direction, and E z is the electric field component in the z - direction;

[0052] Determine the inverse matrix of the 3×3 Jacobian matrix corresponding to each spatial grid:

[0053]

[0054] Among them, are the 9 components of the inverse matrix of the Jacobian matrix;

[0055] Construct a 3N - row and 3 - column matrix C, where every three rows correspond to the inverse matrix of the Jacobian matrix of a spatial grid:

[0056]

[0057] Let g be the number of iterations, and let A g represent the electric field value E at the g - th iteration at time t = (n + 1)Δt g , according to A g+1 = A g - C g .*B g , determine A g+1 , let A g+1 represent the electric field value E at the (g + 1) - th iteration at time t = (n + 1)Δt g+1 , if |E g+1 - E g | > eps, then use E g = E g+1 for iteration, and update the integrated Jacobian inverse matrix C and the integrated matrix B with this electric field value, increase the number of iterations, and the iteration process continues until |E g+1 - Eg | <eps, obtain the electric field E at the new moment (n+1) (t) = E g+1 ; where eqs is the error range.

[0058] The beneficial effects of the present invention are as follows:

[0059] The parallel simulation method for the nonlinear interaction between laser and matter according to the present invention adds nonlinear optical effects such as Kerr effect, multi-photon absorption, multi-photon ionization, avalanche ionization, etc. involved in the interaction between laser and matter, as well as the interaction effect between light and matter (specifically, the free electron plasma excited in the matter) to the Maxwell equations, constructs a coupled equation set composed of the curl equations of the modified Maxwell electric and magnetic fields, the constitutive equation considering the nonlinear effect, and the free electron plasma rate equation, uses the uniaxial perfectly matched layer (UPML) as the absorption boundary, solves the curl equations of the electric and magnetic fields by the finite-difference time-domain method, solves the constitutive equation in parallel by the Newton iteration method, that is, realizes the parallel calculation of multiple CPUs and multi-cores of the computer by matrix reconstruction of the physical quantities related to the light field, uses the Runge-Kutta method to solve the plasma rate equation, and repeats the above algorithm until the simulation is completed as the time step progresses. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. The following drawings are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained obviously without creative efforts based on these drawings.

[0061] Figure 1 It is a flowchart of the parallel simulation method for the nonlinear interaction between laser and matter according to the present invention.

[0062] Figure 2 Modeling diagram of the material structure interacting with the laser (fused silica containing bubbles).

[0063] Figure 3 At t = -27 fs, the laser intensity distribution and the free electron plasma density distribution diagram. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0064] The following further describes the present invention in conjunction with the drawings and specific embodiments, but the protection scope of the present invention is not limited thereto.

[0065] Embodiments of the present invention will be described in detail below. Examples of the embodiments are shown in the accompanying drawings, where like or similar reference numerals denote like or similar elements or elements having like or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to explain the present invention, and should not be construed as limiting the present invention.

[0066] In the description of the present invention, it should be understood that the terms "center", "longitudinal", "transverse", "length", "width", "thickness", "upper", "lower", "axial", "radial", "vertical", "horizontal", "inner", "outer", etc. indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as limiting the present invention. In addition, the terms "first" and "second" are only used for descriptive purposes and cannot be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of the present invention, "a plurality of" means two or more unless otherwise specifically defined.

[0067] In the present invention, unless otherwise clearly specified and defined, the terms "mounted", "connected", "connected to", "fixed", etc. shall be understood in a broad sense. For example, it may be a fixed connection, a detachable connection, or an integral connection; it may be a mechanical connection or an electrical connection; it may be directly connected or indirectly connected through an intermediate medium, and it may be the communication inside two elements. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0068] The present invention is based on Maxwell's curl equation. Optical materials are generally non-magnetic, and μ r = 1,

[0069]

[0070]

[0071] E(t) is the electric field strength, μ 0 is the vacuum magnetic permeability, H(t) is the magnetic field strength, D(t) is the electric flux density, and J is the polarization current.

[0072] As Figure 1 shown, the parallel simulation method for the nonlinear interaction between laser and matter according to the present invention includes the following steps:

[0073] S01: Introduce the plasma rate equation containing the nonlinear optical effect in the interaction between laser and matter into the Maxwell equations as an additional differential equation to obtain the density of the excited free electron plasma. The differential equation is:

[0074]

[0075] In the formula:

[0076] n 0 is the refractive index of the unexcited material; E g is the band gap of the material; n a is the saturated electron density; n e is the density of the free electron plasma excited by the laser; σ is the inverse bremsstrahlung cross section; σ (k) is the multi-electron ionization coefficient; I is the laser intensity; k is the number of laser photons simultaneously absorbed by the material; τ R is the electron-hole recombination time;

[0077] S02: Express the Kerr effect and multi-photon absorption of light by the material, these nonlinear optical effects, as the Kerr polarization intensity P kerr (t) and the multi-photon absorption polarization intensity P MPI (t); and express the influence of the generation of the free electron plasma reflected by the Drude model on the optical properties of the material as the Drude polarization intensity P Drude (t). Construct a constitutive equation based on the Kerr polarization intensity P kerr (t), the multi-photon absorption polarization intensity P MPI (t), the Drude polarization intensity P Drude (t) and the electric flux density D(t). The constitutive equation can be expressed as:

[0078] D′(t) = D(t) + P kerr (t) + P MPI (t) + P Drude (t)

[0079] D(t) = ε 0 ε ∞ E(t)

[0080] Where:

[0081] P Drude (t) is obtained from the following differential equation:

[0082]

[0083] P Kerr (t) is caused by the instantaneous Kerr effect and is expressed as:

[0084]

[0085] P MPi (t) is caused by multi - photon absorption and is expressed as:

[0086]

[0087] Where:

[0088] D′(t) is the electric displacement vector; P Drude (t) is the Drude polarization intensity; P kerr (t) is the Kerr polarization intensity; P MPI (t) is the multi - photon absorption polarization intensity; E(t) is the electric field strength; ε 0 is the vacuum permittivity; ε ∞ is the relative permittivity of the unexcited material; ω p is the plasma frequency; is the third - order susceptibility; Γ is the electron collision frequency;

[0089] S03: Deform the Maxwell's equations using the relationship between polarization current and polarization intensity:

[0090] Polarization current J:

[0091] The deformed Maxwell's equation is:

[0092]

[0093]

[0094] Where:

[0095] E(t) is the electric field strength; μ 0 is the vacuum permeability; H(t) is the magnetic field strength; J(t) is the polarization current;

[0096] S04: Establish a coupled equation set according to the deformed Maxwell's equations, constitutive equations and differential equations from laser parameters and material parameters, and grid the simulation area. Suppose there are i max grids in the x - direction, j max grids in the y - direction, and k max grids in the z - direction;

[0097] The said material parameters include the refractive index n of the unexcited material 0 , the relative permittivity ε of the unexcited material ∞ , the band gap E of the material g , the third - order susceptibility the inverse bremsstrahlung cross - section σ, the multi - electron ionization coefficient σ (k) , the electron - hole recombination time τ R, electron collision frequency Γ, electron effective mass and saturated electron density n a ;

[0098] The laser parameters include the frequency ω, pulse width τ 0 , peak intensity E of single-pulse electric field 0 , beam waist size ω(z) and radius of curvature R(z). The laser can be a Gaussian optical field, Bessel optical field, vortex optical field, etc. In the embodiment, a Gaussian optical field is used, that is, the Gaussian electric field distribution of the laser is a focused beam source in the following form:

[0099]

[0100] S05: Set the anisotropic perfectly matched layer as the absorption boundary, and numerically solve the coupled equations by the finite-difference time-domain method, perform unified numerical simulation on the calculation area and the boundary area, solve the magnetic field H(t) at time t from the electric field E(t) at time t, and solve the electric displacement vector D'(t) at time t from the magnetic field H(t) at time t;

[0101] S06: Numerically solve the constitutive equation by the Newton iteration parallelization method, update the electric field E(t) at time t from the electric displacement vector D′(t) at time t, and realize parallel computing of multiple CPUs and multiple cores of the computer through matrix reconstruction of the physical quantities related to the optical field during the numerical solution process. Specifically:

[0102] S06.1: Determine the three components E x , E y , E z matrix of the optical field vector E depending on the three-dimensional space coordinates on all spatial grids:

[0103] E x (i,j,k) → i = 1 ~ i max

[0104] E y (i,j,k) → j = 1 ~ j max

[0105] E z (i,j,k) → k = 1 ~ k max

[0106] In the formula: i is the grid number in the x direction; j is the grid number in the y direction; k is the grid number in the z direction; i max is the maximum number of grids in the x direction; j max is the maximum number of grids in the y direction; k max is the maximum number of grids in the z direction;

[0107] Construct a one-dimensional column matrix A with 3N rows, where N is the number of spatial grids, and N = i max ×j max ×k max , and in the one-dimensional column matrix A, every three rows correspond to the three components of the optical field of a spatial grid:

[0108]

[0109] S06.2 Write the function values of the three functions [X Y Z] constructed by the constitutive equation on each spatial grid into a one-dimensional row matrix, repeat this row matrix 3 times, thus obtaining a 3-row and 3-column matrix, and then stack all the 3-row and 3-column matrices on all N spatial grids in the column direction to form a 3N-row and 3-column matrix B:

[0110]

[0111] where

[0112]

[0113] In the formula: Δt is the time grid size, and n is the number of time grids; D′ (n+1) (t) is the electric displacement vector at t = (n + 1)Δt, and E (n) (t) is the electric field strength vector at t = nΔt, and E (n+1) (t) is the electric field strength vector at t = (n + 1)Δt, is the Drude polarization intensity vector at t = nΔt, is the Drude polarization intensity vector at t = (n - 1)Δt, is the multi-photon absorption polarization intensity vector at t = nΔt, and E x is the electric field component in the x direction, and E y is the electric field component in the y direction, and E z is the electric field component in the z direction;

[0114] S06.3: Determine the inverse matrix of the 3*3 Jacobian matrix corresponding to each spatial grid:

[0115]

[0116] where are the 9 components of the inverse matrix of the Jacobian matrix;

[0117] Construct a 3N-row and 3-column matrix C, where every three rows correspond to the inverse matrix of the Jacobian matrix of a spatial grid:

[0118]

[0119] Let g be the number of iterations, let A g Represents the electric field value E of the g-th iteration at time t = nΔt g , according to A g+1 =A g -C g .*B g , determine A g+1 , let A g+1 represents the electric field value E at the g+1th iteration at time t=(n+1)Δt g+1 , if |E g+1 -E g |>eps, then use E g =E g+1 Iterate and use the electric field value to update the integrated Jacobian inverse matrix C and the integrated matrix B, increase the number of iterations, and iterate until |E g+1 -E g | <eps,获得新的时间网格下的电场E (n+1) (t) = E g+1 ; where eqs is the error range.

[0120] S07: Runge-Kutta method is used to solve the free electron plasma density rate equation, given by the electric field E (n+1) (t) Update the laser excited plasma density n e , plasma frequency ω p , polarization intensity P drude (t) and P MPI (t);

[0121] With Δt as the time interval, move forward according to the time step, repeating steps S05 to S07 until t> maximum time t max , end the simulation.

[0122] Example 1

[0123] In Example 1, the incident laser has a wavelength of λ=800nm ​​(frequency 3.75×10 14 Hz), pulse width is s p =120fs, single pulse energy of 50nJ femtosecond laser. A microscope with numerical aperture NA = 0.15 is used to irradiate the incident beam on the bubble inside the fused silica (beam radius (waist) ω 0 =3.65μm, radius of curvature R(z)→∞).

[0124] Laser parameters and material parameters are shown in the following table:

[0125]

[0126]

[0127] In the simulations of S01 - S07 in the parallel simulation method for the nonlinear interaction between laser and matter according to the present invention, a modeling diagram of the matter structure interacting with the laser (fused quartz containing bubbles) is simulated, as Figure 2 shown.

[0128] At t = -27 fs, on the cross-section at z = 0 (passing through the center of the bubble), a laser intensity distribution diagram and a free electron plasma density distribution diagram are simulated, as Figure 3 shown.

[0129] It should be understood that although this specification is described according to each embodiment, not every embodiment only contains an independent technical solution. This narrative manner of the specification is only for clarity. Those skilled in the art should regard the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments understandable to those skilled in the art.

[0130] The series of detailed descriptions listed above are only specific descriptions of the feasible embodiments of the present invention, and they are not intended to limit the protection scope of the present invention. Any equivalent embodiments or modifications made without departing from the technical spirit of the present invention should be included within the protection scope of the present invention.

Claims

1. A parallel simulation method for the nonlinear interaction between laser and matter, characterized in that, it includes the following steps: S01: Introduce the plasma rate equation containing the nonlinear optical effect in the interaction between laser and matter into the Maxwell equation as an additional differential equation to obtain the density of the excited free electron plasma; S02: Express the Kerr effect and multi - photon absorption non - linear optical effect generated by the material on light as the Kerr polarization intensity P kerr (t) and the multi - photon absorption polarization intensity P MPI (t); and express the influence of the generation of free - electron plasma reflected by the Drude model on the optical properties of matter as the Drude polarization intensity P Drude (t). According to the Kerr polarization intensity P kerr (t), the multi - photon absorption polarization intensity P MPI (t), the Drude polarization intensity P Drude (t) and the electric flux density D(t), construct the constitutive equation D′(t), and the constitutive equation can be expressed as: D′(t) = D(t) + P kerr (t) + P MPI (t) + P Drude (t) D(t) = ε 0 ε ∞ E(t) where: P Drude (t) is obtained from the following differential equation: P Kerr (t) is caused by the instantaneous Kerr effect and is expressed as: P MPI (t) is caused by multi-photon absorption and is expressed as: In the formula: D′(t) is the electric displacement vector; P Drude (t) is the Drude polarization intensity; P kerr (t) is the Kerr polarization intensity; P MPI (t) is the multi - photon absorption polarization intensity; E(t) is the electric field strength; ε 0 is the vacuum permittivity; ε ∞ is the relative permittivity of the unexcited material; ω p is the plasma frequency; is the third - order polarizability; Γ is the electron collision frequency; S03: Deform the Maxwell equation using the relationship between the polarization current and the polarization intensity: Polarization current J: The deformed Maxwell equation is: In the formula: E(t) is the electric field strength; μ 0 is the magnetic permeability of vacuum; H(t) is the magnetic field strength; J(t) is the polarization current; S04: Establish a coupled equation set based on the laser parameters and material parameters according to the deformed Maxwell's equations, constitutive equations, and differential equations. Mesh the simulation area. Assume there are i max meshes in the x direction, j max meshes in the y direction, and k max meshes in the z direction; S05: Set the anisotropic perfectly matched layer as the absorbing boundary, and use the finite-difference time-domain method to numerically solve the coupled equations, perform unified numerical simulation on the calculation area and the boundary area, solve the magnetic field H(t) at time t from the electric field E(t) at time t, and solve the electric displacement vector D′(t) at time t from the magnetic field H(t) at time t; S06: Numerically solve the constitutive equation in parallel by the Newton iteration method, and update the electric field E(t) at time t from the electric displacement vector D′(t) at time t; S07: Solve the rate equation of free electron plasma density using the Runge-Kutta method, and update the plasma density n (n+1) excited by the laser, the plasma frequency ω e , the polarization intensity P p (t) and P drude (t) and P MPI (t); Advance forward in time steps with a time interval of Δt, and repeat steps S05 - S07 until t > the maximum time t max , and end the simulation.

2. The parallel simulation method for the nonlinear interaction between laser and matter according to claim 1, characterized in that, The differential equation for obtaining the density of the excited free electron plasma is: In the formula: n 0 is the refractive index of the unexcited material; E g is the band gap of the material; n a is the saturated electron density; n e is the density of the free electron plasma excited by the laser; σ is the inverse bremsstrahlung cross section; σ (k) is the multi-electron ionization coefficient; I is the laser intensity; k is the number of laser photons simultaneously absorbed by the material; τ R is the electron-hole recombination time.

3. The parallel simulation method for the nonlinear interaction between laser and matter according to claim 1, characterized in that, Updating the electric field E(t) at time t from the electric displacement vector D′(t) at time t is specifically: S06.1: Determine the three components E x , E y , E z matrix of the light field vector E that depends on the three-dimensional spatial coordinates on all spatial grids: E x (i, j, k) → i = 1 to i max E y (i, j, k) → j = 1 to j max E z (i, j, k) → k = 1 to k max Where: i is the grid number in the x direction; j is the grid number in the y direction; k is the grid number in the z direction; i max is the maximum number of grids in the x direction; j max is the maximum number of grids in the y direction; k max is the maximum number of grids in the z direction; Construct a one-dimensional column matrix A with 3N rows, where N is the number of spatial grids, and N = i max ×j max ×k max , and every three rows in the one-dimensional column matrix A correspond to the three components of the optical field of a spatial grid: S06.2 Write the function values of the three functions [X Y Z] constructed from the constitutive equation on each spatial grid as a one-dimensional row matrix, repeat this row matrix 3 times, thereby obtaining a 3-row 3-column matrix, and then stack all the 3-row 3-column matrices on all N spatial grids in the column direction to form a 3N-row 3-column matrix B: where where: Δt is the time grid size, and n is the number of time grids; D′ (n+1) (t) is the electric displacement vector at the time t = (n + 1)Δt, and E (n) (t) is the electric field strength vector at the time t = nΔt, and E (n+1) (t) is the electric field strength vector at the time t = (n + 1)Δt, is the Drude polarization intensity vector at the time t = nΔt, is the Drude polarization intensity vector at the time t = (n - 1)Δt, is the multi-photon absorption polarization intensity vector at the time t = nΔt, and E x is the electric field component in the x direction, and E y is the electric field component in the y direction, and E z is the electric field component in the z direction; S06.3: Determine the inverse matrix of the 3*3 Jacobian matrix corresponding to each spatial grid: Among them, are the 9 components of the inverse matrix of the Jacobian matrix; Construct a 3N-row 3-column matrix C, and every three rows correspond to the inverse matrix of the Jacobian matrix of a spatial grid: Taking.g as the number of iterations, let A g represent the electric field value E at the g-th iteration at time t=(n + 1)Δt g . According to A g+1 =A g -C g .*B g , determine A g+1 . Let A g+1 represent the electric field value E at the (g + 1)-th iteration at time t=(n + 1)Δt g+1 . If |E g+1 -E g |>eps, then use E g =E g+1 for iteration, and update the integrated Jacobi inverse matrix C and the integrated matrix B with this electric field value, increase the number of iterations. The iteration process continues until |E g+1 -E g |<eps, and obtain the electric field E (n+1) (t)=E g+1 ; where eqs is the error range.

Images