Full Explicit Numerical Modeling in Time-Domain for Nonlinear Electromagnetics Simulations in Ultrafast Laser Nanostructuring

: The purpose of this paper is to present a new and accurate, fully explicit ﬁnite-difference time-domain method for modeling nonlinear electromagnetics. The approach relies on a stable algorithm based on a general vector auxiliary differential equation in order to solve the curl Maxwell’s equation in a frequency-dependent and nonlinear medium. The energy conservation and stability of the presented scheme are theoretically proved. The algorithms presented here can accurately describe laser pulse interaction with metals and nonlinear dielectric media interfaces where Kerr and Raman effects, as well as multiphoton ionization and metal dispersion, occur simultaneously. The approach is ﬁnally illustrated by simulating the nonlinear propagation of an ultrafast laser pulse through a dielectric medium transiently turning to inhomogeneous metal-like states by local free-electron plasma formation. This free carrier generation can also be localized in the dielectric region surrounding nanovoids and embedded metallic nanoparticles, and may trigger collective effects depending on the distance between them. The proposed numerical approach can also be applied to deal with full-wave electromagnetic simulations of optical guided systems where nonlinear effects play an important role and cannot be neglected. and Science Eclipse Parallel Tools Platform. Conﬂicts of Interest: The authors declare no conﬂict of interest.


Introduction
The finite difference time-domain (FDTD) method enables the computation and modeling of light propagation and scattering processes in linear and nonlinear dispersive materials possessing wavelength-dependent and intensity-dependent properties. This ranks the method among the most universal and powerful numerical tools for optics, electrodynamics, antennas, and waveguides theory. In this context, the general vector auxiliary differential equation (GVADE) method, which has been reported in [1] and improved in [2], is one of the most popular methods to deal with nonlinear media. However, the method suffers from some drawbacks that limit its speed and stability due to the semi-implicit updating relation for the current polarization densities. In summary, the method is not appropriate or recommended for a broad range of problems, especially those that address the designs with complex 3D topology, where iterative solutions of coupled nonlinear equations can be unstable and resource-demanding [3]. Some direct or explicit methods take into consideration the Kerr effect in the Born-Oppenheimer approximation [4]. The Kerr effect arises from the nonlinear third-order electric susceptibility. In [3], the authors take into consideration the Kerr effect in an explicit form. Other direct methods [5] require excessive computational burdens, especially for the consideration of complex geometries. When metal-dielectric interfaces are treated, especially at nanometric scales, plasmonic effects appear [6,7] and generate high concentrations or high local densities of the electromagnetic field [8]. Since dielectrics experience nonlinear optical processes, the treatment of these regions in the interfaces is particularly challenging.
In this manuscript, we develop an algorithm that permits solving a nonlinear dielectric medium in which Kerr and Raman's effects [9] occur simultaneously. The Raman effect arises from the scattering caused by the conversion of photon energies into vibrational energy of molecules. The scattered photons have less energy than the incident photons corresponding, therefore lowering frequency fields [2]. We also assume that this takes place in an inhomogeneous dielectric environment that can potentially include metallic heterogeneities with specific optical properties. Along with the laser propagation, provided that the intensity is sufficient to ensure nonlinear absorption, the medium can also turn into a metallic state as a free electron population is generated during multiphoton absorption (MPA) or tunneling and avalanche ionization processes. It is well known that the electron transition from the valence band (VB) to the conduction band (CB) is followed by various processes, including heating in the CB, impact ionization, and various kinds of electron collisions. The purpose of this work consists of computing all these events in a full explicit form. The carrier generation rate G( r, t) has been computed owing to the BVkP approach implying Bloch-Volkov states [10]. The algorithms presented in this manuscript are especially suitable for simulations dealing with confinements and guiding high-energy pulsed beams. In summary, in this work, we coupled Maxwell equations and the electric charge carriers conservation equation in which the quantum ionization and recombination are governed by BVkP self-consistently.
In the next sections, by following and enriching some pioneers' works [11][12][13], we detail a fully explicit, stable algorithm based on GVADE FDTD, able to solve the curl Maxwell's equation in a frequency-dependent way for a nonlinear medium. In order to illustrate the combination of all these nonlinear effects occurring simultaneously in a simulation, we have modeled a Bessel beam traveling through dielectric fused silica, where some metal inclusions are modeled in a similar form as in Ref. [14]. In Section 2, we introduce the theoretical model, which describes the nonlinear effects. Section 3 presents the algorithms that discretize in time and space the main magnitudes, electromagnetic fields and electric charge. Appendix A is devoted to the stability condition assessment for these algorithms, which imposed numerical limits to the simulation performance. These limits are studied in Section 4, where we specify the validity framework considering natural limits imposed by physics. In Section 5 we present the simulation outcomes provided by the full numerical modelling and in Section 6 the conclusions are drawn. Finally, the stability conditions are derived in Appendix A and the used computational resources are presented in Appendix B.

Modeling Electromagnetic Fields in Nonlinear Media
The Equations (1)-(3) model the laser beam propagation by means of the electromagnetic fields E( r, t) and H( r, t). Throughout the light propagation path, fields generate and couple with the density of free charges as free electrons n f e − ( r, t), which in turn interact with the local field. In the given form, nonlinearities are determined by the electric displacement field D( r, t) in the Ampère-Maxwell Equation (1). Meanwhile, the evolution of the free-electron density or medium conductivity is described by Equation (3), where G( r, t) is the carrier generation rate and R( r, t) is the carrier recombination rate.
The electric displacement field is composed by the electric field contribution plus the polarization vectors D( r, t) = ε 0 ε ∞ E( r, t) + ∑ P( r, t) [2]. The magnitude ε ∞ is the non-unity high-frequency relative permittivity of the media. On the right hand side of the displacement electric field, the sum of polarization vectors in the model that we are considering ∑ P( r, t) is the contribution of four terms, the metal polarization P Metal ( r, t), the Kerr effect polarization P Kerr ( r, t), the Raman polarization P Raman ( r, t) and the free charge polarization vector P f e − ( r, t) such as ∑ P( r, t) = P Metal ( r, t) + P Kerr ( r, t) + P Raman ( r, t) + P f e − ( r, t) with the detail of each term below: is a coefficient gathering parameters depending on the considered medium. Expressions Equations (5) and (6) providing the Kerr effect polarization and Raman polarization, respectively, are taken from [9]. Here, α is a real-valued constant in the range α ∈ [0, 1] that parameterizes the relative strengths of Kerr and Raman contributions, δ(t) refers to the Dirac delta function that models the Kerr nonresonant virtual electronic transitions, which is medium dependent and in the order of few femtoseconds for fused silica, U(t) is the Heaviside unit step function, e −t τ 2 sin( t τ 1 ) models the impulse response of a single Lorentzian relaxation centered on the optical phonon frequency 1/τ 1 and having a bandwidth of 1/τ 2 . Finally χ (3) 0 is the strength of the third-order nonlinear electric susceptibility.
The upper limit of the integrations in Equations (5)-(7) extends only up to t because the response function must be zero for τ > t to ensure causality [3]. Moreover, freecharge polarization is introduced by using a free-charge density of current, which will be defined later. Finally, to model the metal, we consider a complex electrical permittivity (ω) = a jω+b where a and b are complex numbers that fit experimental data [15]. To model the dispersive media, we assume an electric permittivity, which is frequency-dependent as the calculations are devoted to ultrafast laser propagation, with bandwidth-limited pulses that exhibit wavelength dispersion.
In ref. [2], the general vector auxiliary differential equation technique allows considering nonlinearities or dispersive characteristics of a medium. They are introduced into the Ampère-Maxwell equation by means of current densities. These currents correspond intrinsically to the temporal variation of the polarization vectors. Thereby, the following Equations (8)- (11) can be derived to correlate the currents J Metal ( r, t) = ∂ P Metal ( r,t) with the electric field: We define a free-charge density of current as an Ohmic current, where the conductivity σ f e − ( r, t) = qµ f e − n f e − ( r, t) is given by the following frequency dependent expression: where q is the electron charge, µ f e − the free carriers of charge mobility, n f e − ( r, t) the freecharge density and n ph−e − ( r, t) the photo-ionized electron density. This definition for the density of current Equation (11) is similar to Equation (9) in [10]. In Equation (13), we employ free-charge which is governed by Equation (3), or generation Equation (13) and recombination Equation (14) rates. In [10], the authors explain that electron density n ph−e − ( r, t) produced in the conduction band only describes the interaction of an electron with the laser electric field and accounts for the lattice periodicity, but does not account for possible collisions with phonons, ions, or other electrons. Due to these collisions, the coherence between the excited electron to the CB and its parent ion is destroyed and may introduce deviations from the above-predicted evolution of the electron density in the CB with respect to time. Indeed, without collisions, the produced electron density was shown to oscillate with time as a conduction electron may go back and forth to the VB through the action of the electromagnetic field [10,16]. Obviously, in an oscillation, the energy stays constant [17]. Hence, the authors of [10] defined a free state as a state where the electrons can be heated by inverse bremsstrahlung (IB) and are decorrelated from the parent ion. The density of charge associated with this free state is called here density of freecharge n f e − ( r, t), and its evolution along the time is modeled by Equation (3). Nevertheless, the µ f e − , the free carriers of charge mobility, then we employ experimental data, which allow the determination of this parameter. The method is explained in detail in Section 5.
R( r, t) = f r ( r, t)n f e − ( r, t) In the carrier generation rate, the term associated with the photo-ionized electron density n ph−e − ( r, t) produced in the conduction band is a function of | E( r, t)|.
Finally the functions related to the generation by ionization f d ( r, t), impact f η ( r, t) and avalanche f a ( r, t), as well as the function related to the recombination of the free-charge f r ( r, t), are given by the Equations (16)- (19), respectively.
In Equations (16)- (19), we find that the factor of maximum electric field reduction, which is suggested twenty times by [10], is taken into consideration. Nevertheless, when the photoionization process is considered in a truncated region where the magnitude of | E( r, t)| is from the beginning higher than the threshold E max 20 reference value, these logical conditions can be suppressed in the nested loops, and this increases the algorithm performance. In other words, this simplification is possible when we assume the description with nonlinear effects in a particular region of the computational domain. So, the electromagnetic field that arrives at that region has a magnitude higher than the threshold E max 20 .
In Equation (18), we assume that is the fluence with ζ being the medium impedance. However, this is a particular case only valid when the electromagnetic field is a harmonic function in time. In general, and assuming a propagation along the z-axis, the power per surface is given by the Poynting vector.

Algorithms
The polarization vector due to the Kerr effect at the time t = ∆ t s − 1 2 can be written ; however, afterwards, at the instant ∆ t s + 1 2 , in order to define a full explicit approach, we have to make the approximation ( This approximation was demonstrated in [3] and leads to: This approximation P 1 2 is valid under particular physical conditions explained in Section 4. To compute the nonlinear Raman effect, we make an analogous approximation. Considering that in this case there is no instantaneous cause-effect relationship, it is evident that the approximation is justified [18]. Therefore, by defining the scalar magnitude | Q a ( r, t)| that is defined from the vector Q a ( r, t), which could be considered as an instantaneous anisotropic Raman conductivity, we can discretize the Raman polarization vectors at the instants ∆ t s − 1 2 and ∆ t s + 1 2 as follows: where the Raman response function for a particular medium is modeled by i n i , depends on the refractive index of the media n i as well as the location of the electromagnetic field in the media ( Figure 1). Here u(t) is the Heaviside step function: Scattering body TF/SF CPML Figure 1. Scheme of the calculation box, illustrating the delay nt 0 (k) that depends on the refractive index of the media n i as well as the location of the electromagnetic field in the media.
In this way, we obtain the Raman current density and the Kerr current density: The current density updating associated to the metal is given by the equation: where g a = R 2−b∆ t 2+b∆ t and g b = R a 0 2+b∆ t . The operator R[] calculates the real part of these complex numbers. Finally, the electric field updating is given by the algorithm:

Metal
where the parameter C s B has different time-dependent values at the different locations: where ∂¯ implicates any region in metals or CPMLs, meanwhile D − ∂¯ represents any computational region filled by a dielectric medium. Finally we have two parameters, g c and Υ( r), that are defined as follows: Υ permits to switch the algorithm from dielectric to metal and vice versa. g c is a coefficient gathering parameters for convenience in the metal current. Table 1 summarizes the algo-rithms presented in this section, and Figure 2 sketches the general flow chart. In the next section, we discuss the stability conditions for these algorithms. Table 1. The list of the new explicit algorithms added to the classical FDTD scheme in order to deal with the nonlinear effects. The current list of algorithms is depicted in the form of a flow chart and it is complemented by the global integration view shown in Figure 2.
Step 0: Update 6D-Discrete Incident Electric Field Array ↓ Step 1: Total field / Scattering field algorithm on E ↓ Step 2: Step 9: CPML on Magnetic Field ↓ Step 10: Update 6D-Discrete Incident Magnetic Field Array ↓ Step 11: Total field / Scattering field algorithm on H ↓ Step 12: In Appendix A, we determine the stability conditions for the presented algorithms. Table 2 summarizes the stability conditions of all these effects together or separately.

Performance and Framework
Some approximations made here are not always applicable and, in this section, we will establish the validity framework of the proposed algorithms. In particular, we should con- Both expression are derived in [9] under the Born-Oppenheimer approximation [4] where the Kerr non-resonant virtual electronic transitions are in the order of 1fs, so it can be considered an instantaneous event, which is modeled by the Dirac delta function. However, while the time stepping keeps beneath 1 fs, we could model the medium reaction by using a time narrow Gaussian function with a standard deviation of 1 fs. In our simulations, ∆ t ∼ 0.1 fs, so we are adding in a fraction of 1 10 the standard deviation. In this particular situation, the approximation is assumed to be valid. We employ the full explicit form | Q This expression is valid for all instant ∀t. Hence, we fix the limit for this full explicit approach in the relation ∆ t < 1 fs. The advantage of a fully explicit approach concerns the simulation time. Due to the updating connection in an implicit scheme, most of the algorithms should be re-updating to advance a time step. In classical FDTD, it is usually fixing the limits of the index in all nested loops. In a significant computational domain, for instance, 1.2 G cells, the updating demands 30 s. If the algorithm has to be repeated at least three times, the simulation should need around three times the the implicit scheme. Figure 3 shows the number of simulation status updates on the abscissa axis and on the ordinate axis the time needed in order to update these simulation states. Each step or update corresponds to executing the 14 algorithms shown by Figure 2. In addition to these algorithms that are necessary to perform a system status update, two numerical processes are employed. One exports simulation data or results (every 40 updates) and the other exports the complete simulation status (every 500 updates). The latter is done for two reasons. The first is to recover a state of a simulation in case of computer collapse (for example a failure in the computer electric power source). The second reason concerns the alteration of the computational domain (electromagnetic properties or/and size of the computational domain) without the need to restart the simulation from the first time stepping. In Appendix B the Table A1 lists the computational resources employed in the simulations. Under these conditions we reach a performance which can be summarized by the computational burdens and some performance parameters. Those relies on a computational domain size of 121.  Figure 3 depicts the time performance for a simulation carried out with dynamics nested loops. In our scheme, the size of the nested loops is provided by the information created. We know that our implementation violates the second rule of NASA's Ten-Rules for Developing Safety Critical Code [19]. In general, it is a good practice to arrange the loop indexes in the simulation beginning and keep these ranges. By doing this, we fix, from the start, the nested loop size. Nevertheless, the employment of a dynamic nested loop saves simulation time. Figure 3 depicts that we save more than forty hours. In a simulation with the same characteristics and static nested loop, we should have a constant duration per time step updating of 32 seconds which is equal to the asymptotic value in Figure 3. The total simulation duration is the accumulated time or area under the plot schematized in Figure 3. Hence, the triangle in gray color is the time saved by means of dynamic nested loops.  Finally, we close this section with some consideration on Equation (3). This equation is a derivation from the charge/electrons continuity equation [20], which has the general form where we neglect the term ∇ · J( r, t) qn f e − ( r, t). We know from theory that ρ( r, t) = qn f e − ( r, t) and J( r, t) = qµ f e − n f e − ( r, t) E( r, t). The key is the extremely low mobility µ f e − ∼ 10 −22 , which is deduced from [21]. We can consider a highly intense laser beam with | E( r, t)| 2 ∼ 10 13 , which induces a free electron concentration in fused silica of n f e − ( r, t) ∼ 10 26  Therefore, we can conclude that n f e − ( r, t) ∼ 10 8 ∇· J( r,t) q and Equation (3) is a valid approximation in the present context.

Results
An efficient technique to compute scattered fields in the context of FDTD modeling is the total-field/scattered-field (TF/SF) incident wave source, which is employed by most all current commercial FDTD solvers [2]. Fundamentally, the TF/SF technique is an application of the well-known electromagnetic field equivalence principle [22]. By this principle, the original incident wave of infinite extent and arbitrary propagation direction, polarization, and time-waveform is replaced by electric and magnetic current sources appropriately defined on a finite closed surface, called the Huygens surface [23], containing the object of interest. The reformulated problem confines the incident illumination to a compact total-field region and provides a finite scattered-field region external to the total-field region that is terminated by an absorbing boundary condition (ABC) to simulate the FDTD grid extending to infinity. In particular, we use convolutional perfectly matched layers, a technique already used in other works [24]. In this work, we have used the finite-difference time-domain discrete plane wave technique (FDTD-DPW), which is a numerical approach based on TF/SF technique that allows one to propagate plane waves quasi-perfectly isolated [25]. This means a propagation in the total field domain without reflections to the scattered field domain on the order of machine precision (∼300 dB) [26]. This technique is valid for any angle of propagation, for any grid cell aspect ratio and even for nonuniform grids [27]. Figure 4 illustrates the total field region where there are two important elements, an axicon lens that generates the Bessel beam from the initial plane waves and a block of fused silica. In this domain, the incident plane wave is generated by the FDTD-DPW scheme. The calculation domain contains a block of fused silica glass where the nonlinear effects apply. The interaction between the incident plane wave and the axicon lens results in a Bessel beam, which impinges the SiO 2 block. We assume a pulsed laser source of wavelength (800 ± 11) nm and a time full width at half maximum of σ t = 40 fs. The pulsed incident time-waveform is modeled by with ς ≥ 1, is an appropriate way to avoid numerical instabilities in the simulation beginning. The incident field six-dimensional array (IFA) that feeds the FDTD-DPW is quite simple E( r IFA,t ) = E 0 h(t)x. In addition, there are other parameters that characterized the media properties and are summarize in Table 3.

Name
Symbol Value Unit Ref.
Gap energy E g 9 eV [28] Third-order elec-tric susceptibility χ Axicon refractive index n a 1.45 - [30] By using the developed irradiation strategy, the computational domain elements and their properties as well as the algorithm we carried out, the obtained results can be divided into three groups:  I( r, t)) is calculated as the maximum power that flows through a section of σ 2 w , being σ w the pulsed laser transverse waist, σ the pulsed length in the propagation direction and z p the central pulse location at the calculation time. Hence, To validate the algorithm, which includes nonlinear effects during light propagation, we have investigated the features of the propagation regime depending on the laser intensity. Figure 5a,b redraw the distribution evolution of the pulsed square electric field module for simulation performed for different laser intensity sources. Figure 5a plots the iso-slices that depicts | E| 2 inside the fused silica bulk for high laser power in the system. From top to bottom, the time evolution sequence reveals the defocusing process of the Bessel beam, which ends by forming some filamentation structures in the light distribution. At the time t = 5720, ∆ t = 1.52 ps, the pulsed square electric field module begins to penetrate the interface air-fused silica. All successive time plots correspond to a propagation inside the fused silica block. Figure 5b illustrates, for the same time moments, but a higher initial power source, the square electric field module during the propagation inside SiO 2 material.
We complement the picture with a simulation carried out in the medium power regime. From the instant of time, t 1 = 1.67 ps to the instant t 18 = 2.56 ps, Figure 5c illustrates the antagonistic trade-off between the focusing and the defocusing process experienced by the Bessel beam along with its travel through the bulk fused silica target. Figure 6b shows the history of the power density at the critical threshold separating medium and high power regimes, along with a few steps of time. In this figure we can identify the interface air-fused silica due to the rings which record the history of a stationary wave formed after the interaction. We can observe a strong focusing process in the entrance of the fused silica block followed by an alternation of refocusing and the defocusing events induced and governed by the photoionization and Kerr effect competition, and we know that these are not diffraction artifacts because in the air, at the same laser intensity, we do not observe them. At this stage, we have to give a clear definition of the pulsed power density and the power history. The magnitude of pulsed power density is defined as ∇ · S( r, t), where S( r, t) = E( r, t) ∧ H( r, t) is the Poynting vector [31]. Figure 6a sketches the pulsed power density. From this magnitude ∇ · S( r, t), we can calculate the delivered pulsed electromagnetic density of energy t 0 ∇ · S( r, t)dt. Besides this, it is interesting to consider the range of values in Figure 6a. We can identify a negative power density that requires clarification. As the Poynting vector represents the flux of power through a given surface, the divergence of this magnitude represents the power density that comes into/out a particular point in a given region. In particular, in FDTD calculations, we are describing the power flux through the computational nodes that discretized the computational domain. In the density power accounting, the reference system is placed inside the cell that encloses the nodes. In this way, a negative power density, in a particular instance of time, in a local region, means a local sink of power ∇ · S( r, t) < 0 in that region, and vice versa, we have a local source of power as ∇ · S( r, t) > 0. These oscillations in the density power magnitude indicate a redistribution that is able to support the pulsed density power and the pulsed density energy traveling.   Turning to the matter at hand, we define the history of the power density by applying the absolute value to the divergence of the Poynting vector |∇ · S( r, t)|.
The bi-dimensional longitudinal cut of the history of the power density depicted in Figure 6b permits one to identify the air-fused silica interface. There, in that interface, we can appreciate the reflected beam and the strong focusing in the entrance due to the Kerr effect. We assume the same refractive index for both the axicon lens and the fused silica block in which the wave Bessel beam impinges.
Although we assume nonlinear media, in the simulation results presented up to this point in the manuscript, we have considered that the fused silica is isotropic and homogeneous. However, in order to study the effect of micro-and nanopores in the fused silica, 650 regions or voids have been introduced that range in size from 200 to 650 nanometers in radius. These spherical embedded cavities, that are plotted in Figure 7a, emulate a nanoporous fused silica bulk. The traditional method of linear combined congruence random number generators has been used to assign the size and location of these pores of air/void particles. The linear congruential generator yields a sequence of pseudorandomized numbers calculated with a discontinuous piecewise linear equation. The most important thing about the method is that, under the same seeds, the same random number distributions will be obtained, which allows us to reproduce both size distributions and particle locations. This will be useful to evaluate the effects resulting from another kind of media, such as metal particles, which will allow us to make a point-to-point comparison and infer conclusions from the results. In this setup, Figure 7b shows the history of the power in the region where the particles are located. The simulation is stopped when the nonlinear Bessel beam exits the porous fused silica region. Looking at Figure 7c, it appears that a correlation can be made between the distribution of photoionized electrons and the laser power density, in full consistency with our model. Figure 7c shows a history of the density of the ionized electrons as a consequence of the Bessel beam, which travels through the fused silica block. This result leads us to conjecture that both quantities look similar in terms of spatial structure.
It is relevant to study the distribution or history of power in a non-homogeneous dielectric, as is the present case. If we observe Figure 7b in detail, we see that the history of power reveals the history of field intensity as it passes through the nanoporous medium of SiO 2 . Due to a simple matter of impedance, the wavefront of the pulse is redistributed, avoiding entering inside the voids. As we can see in the longitudinal section shown by Figure 7b, the power vanishes inside the voids and the maximum power is not found at the interface fused silica-void. The maximum is reached in the glass medium between the pores. There is a visual explanation for this behavior. In some lenses, an anti-coating multi-layer that matches the impedance of the lens with the surrounding media (in general, glasses could be done for water or air), is coupled in order to remove the reflections. This gradient in the refractive index is necessary to facilitate the flowing energy from the lens to the environment. In the porous fused silica, there is an abrupt change in the refractive index that explains our result.  As stated above, we have followed identical reasoning and repeat the simulation, this time replacing nanovoids with metal particles. The results are shown in Figure 8. Figure 8a illustrates the power history under the same intensity of light. Along the photoexcitation path, the maximum power occurs around the metal-fused silica interface, enhancing and localizing the light coupling to the nanocomposite material. Even if the non-diffractive behavior is still present, it is less observable, as the beam is strongly attenuated by metal absorption. Moreover, the arrangement of the higher region of light absorption follows the randomly distributed composition of the metallic sphere. This indicates that contrary to dielectric-void interfaces, the dielectric-metal region fosters local effects, preventing any collective response of the medium. This statement may differ if the particles are judiciously defined in terms of concentration and size in order to enhance resonance effects as plasmonic modes. In Figure 9, we observe a pattern in the electron photoionized history. In particular, Figure 9a depicts a weak structure of parallel planes, which is more intense in some regions. Figure 9b shows a stronger pattern. The present results, sketched by Figure 9a, reveal that ultrashort laser pulse interacting with distributed performed scattering nanovoids induces local field enhancement that results in enhanced nonlinear absorption and localized optical breakdown. At the void-SiO 2 interface, the multiphotonic ionization process transforms the material into a thin layer of absorbing plasma with metallic-like properties. The light coupling initiates the growth of nanoplates of enhanced energy deposition, perpendicular to laser propagation direction, but parallel to the polarization. This kind of arranged structure exhibiting a periodicity approximately the wavelength of light in fused silica was commonly experimentally reported in conditions of tightly focused and high repetition rate [32,33]. They directly emerge from the spatial coherence of the waves scattered from the inhomogeneous material in the plane perpendicular to the beam propagation direction. Note the relatively used large nanovoids associated to the grid resolution do not allow to discriminate nanogratings of higher frequency perpendicular to the electric field. For propagation inside laser-induced nanoporous media, the fused silica glass remains underdense, below the critical plasma concentration, and the standing wave effect remains weak. However, the other simulation conditions, portrayed in Figure 9b), clearly indicate that for metallic nanoparticles, local field enhancement occurs near the metal-dielectric interface [34]. Upon electron-matrix energy coupling and void deformation, the light pattern may result in elongated nanostructures [35]. This shows that nanoplasmonic behavior dominates and opens the route for advanced applications as embedded micro-reflectors for a stronger near-field interaction regime. Propagation and scattering through sphere inclusions with different optical properties demonstrate that the current algorithms are stable and can model efficiently nonlinear media with sharp dielectric interfaces among them. Moreover, nonlinear metal-dielectric interface investigations are far reaching for technology-based photonics. For instance, in a wider frequency spectrum, resonators, metasurfaces with designed cavities, metaldielectric waveguides or even coating optical fibers are composed of metal-dielectric interfaces that may involve simultaneously nonlinear mechanisms of different nature. If we imagine for a sake of simplicity an harmonic electromagnetic pulse, which is guided, for instance, inside of a optical fiber, we should solve two coupled inhomogeneous Helmholtz's equations in the frequency domain ∇ 2 E( r) + ω 2 µ ( r) E( r) = −∇ E( r) · ∇(ln( ( r))) and ∇ 2 H( r) + ω 2 µ ( r) H( r)( r) = (∇ ∧ H( r)) ∧ ∇ ln( ( r)). In these equations, the refractive index n( r) = r ( r)µ r , which is a dynamic magnitude along the pulsed beam propagation or simply a nonlinear magnitude, has to be taken into consideration. In this way, the algorithms present here can be used to deal with this kind of problems.

Conclusions
In this work, we have proposed a numerical approach for nonlinear laser beam propagation based on the time-domain discretization of Maxwell equations in a fully explicit scheme, including four nonlinear optical effects. The developed algorithms have been benchmarked in three different ways. First, we have studied the algorithm's stability conditions by combining the von Neumann method with the Routh-Hurwitz criterion. Second, we have established the physical framework for the algorithms, with special attention given to the approximation carried out in the algorithm derivation. Third, we have employed the algorithms in order to determine the transient optical properties of a photoexcited media when a laser beam is focused inside a fused silica bulk. In particular, we have considered a Bessel beam as a free-diffraction beam due to its stable and rebuildablelike energy distribution along with electromagnetic pulse travel path. This non-diffractive behavior has been observed numerically by tracking the power/energy history. This has been determined by evaluating the absolute value of the power |∇ · S( r, t)| and the density energy history as t 0 |∇ · S( r, t)|dt. The main advantage of the present method consists of its remarkable numerical performance (see Figure 3), that together with the tolerant stability conditions found out and the accuracy results, under an also permissive physical framework, made it valid for a wide range of applications (such as ultrafast laser processing which has been extensively applied to fabricate 3D photonic devices based on embedded optical waveguides inscribed in glass materials by inducing permanent refractive index changes in the focal volumes [36], optical couplers and splitters [37], volume Bragg gratings [38], or even diffractive lenses [39] which have been efficiently inscribed by refractive-index changes during light propagation) in which optical/electromagnetic nonlinear effects have to be taken into consideration. Finally, our approach is generalizable to any nonlinear electromagnetic problem with dielectric-metal interfaces. The proposed work elucidates the simultaneous treatment of nonlinear effects and provides new routes toward the simulation of light-matter interaction with sub-wavelength features for the optimization of guided structures or loss-free metasurfaces.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Stability
In this appendix, we study the stability of the algorithm by means of the combination of the von Neumann method with the Routh-Hurwitz criterion [40]. The von Neumann method employs a Fourier series expansion of the error at the mesh nodes at a given time instant t = s∆ t . Due to linearity, only a single term of this expansion needs to be considered, i.e., E s | i,j,k = E 0 y s e j(i∆ x k x +j∆ y k y +k∆ z k z ) (A1) D s | i,j,k = D 0 y s e j(i∆ x k x +j∆ y k y +k∆ z k z ) where E 0 and D 0 are complex amplitudes, indexes i, j, k denote the position of the nodes in the mesh, ∆ x , ∆ y , ∆ z are the sizes of the discretization cell, and k x , k y and k z are the wavenumbers of the discrete modes. The factor y in (A1) and (A2) is a complex variable, often called the amplification factor, which gives the growth of the error in a time updating.
To ensure that a finite-difference algorithm is stable, the error must not grow as the time increases, thus, the condition |y| ≤ 1 must be satisfied. The stability condition for a particular FDTD scheme leads to a characteristic polynomial in y: Assuming that y are the roots of Φ(y), the condition for stability implicates that the roots y must lie inside or on the unit circle in the Y-plane. In this stability purpose, the Routh-Hurwitz criterion provides a method that ensures the roots y on the unit circle in the Y-plane. It establishes that the polynomial Φ(r) = ∑ N p =0 b y N p − with constant real coefficients b has no roots in the right-half of the r-plane if all the entries of the first column of the Routh table are non-negative quantities. In [40], there is a error that makes it impossible to build the Routh table, so we explain here how to build up the Routh table. The first and second rows of the coefficients correspond to even and odd powers of r, respectively. For simplicity in the notation, and without loss of generality, it is usual assuming N p to be an even number, therefore c 1, = b 2 , with = 0, 1, 2, . . . , . The Routh-Hurwitz criterion can be used to determine if any root of the stability polynomial Φ(y) lies outside the unit circle in the y-plane. To this end, the bilinear transformation y = r+1 r−1 is carried out in Φ(y). This transformation maps the exterior of the unit circle in the y-plane onto the right-half of the r-plane. To satisfy the von Neumann stability condition for a different scheme, as a function of the parameters of interest, all the entries of the first column of the Routh table are forced to be non-negative quantities. Nevertheless, before beginning the search for stability conditions, some definitions are necessary. We define the following operators applied to Ξ s :

Appendix A.2. Stability Condition on Kerr Effect
In addition to the wave equation, we have to deal with the differential equation associated to the Kerr effect ∂ ∂t D( r, t) = 0 ∞ ∂ ∂t E( r, t) + ∂ ∂t P Kerr ( r, t). The discretization of this equation by employing the operators results tō It is evident that | E s− 1 2 | 2 = E 0 E * 0 y 2s−1 and this result leads to the characteristic polynomial Φ(y) = (y − 1) 2 ( 0 ( ∞ + αχ (3) 0 | E 0 | 2 y 2s−1 )) + 4y 0 ∞ ν 2 Comparing this characteristic polynomial with this one obtained from the wave difference differential equation Φ (y) = (y − 1) 2 0 ∞ + 4y 0 ∞ ν 2 , we see that they are similar when ∞ αχ (3) 0 | E 0 | 2 y 2s−1 . Obviously, we assume that |y| < 1, but notice here that we are not assuming the same thing that we want to demonstrate. We know that Φ (y), under the classical stability condition, leads to |y| < 1 and in general for ∀s ≥ 1 2 we find that Φ(y) ≡ Φ (y) when ∞ αχ (3) 0 | E 0 | 2 . We can identify the problem at a glance. The stability depends on the laser intensity | E 0 | through its propagation. In that case, we find the classical condition for stability already presented. Therefore, we can impose a condition to the laser intensity in order to preserve the stability; for instance, we can define an arbitrary factor D f such that ∞ D f αχ (3) 0 | E 0 | 2 . Therefore, to ensure the stability . From the last expression, we see that stability is related to the medium, optical properties as well as laser intensity. Hence a heuristic search helps to establish an asymptotic limit for this factor D f ≥ 10 4 .

Appendix A.3. Stability Condition on Raman Effect
Again, we consider the wave equation plus the below difference equation for the Raman effect, By applying the operators, we determine the discretized expression: We calculate the characteristic polynomial replacing equations Equations (A1) and (A2) in the discretized wave equation and in the above expression. Then, we combine both derivation in the same formula, Φ(y) = (y − 1) 2 ∞ + y 2s−1 Ψ y s∆ t 0 g(s∆ t − τ)dτ Under an approach analogous to the one utilizes in the previous case we have to link the power laser to the medium. That means that the scheme should verify the inequality, ∞ y 2s−1 Ψ y s∆ t 0 g(s∆ t − τ)dτ + ∆ t g(s∆ t ) 2 (y 2 + y + 1) | E 0 | 2 .
We know that we observe that stability is related to the medium optical properties as well as the laser intensity when the Raman nonlinear effect is present.

Appendix A.4. Stability Condition on Photoionization Effect
Together the wave equation we take into consideration the differential equation ∂ ∂t D( r, t) = ∞ 0 ∂ ∂t E( r, t) + σ f e − ( r, t) E( r, t), which by means of the defined operators can be discretized as follows,∂ t D s = ∞ 0∂t E s + σ f e −∆ E s Again we combine the equation with the discretized wave differential equation with the aim to determine the characteristic polynomial Φ(y) = 2 0 ∞ (1 + y 4ν 2 − 2 + y 2 ) + ∆ t σ f e − y 2 − 1 in which we apply the bilinear transformation y = r+1 r−1 . Therefore, we calculate the characteristic polynomial Φ(r) = 2 0 ∞ ν 2 r 2 + σ f e − ∆ t r + 2 0 ∞ 1 − ν 2 for which the Routh-Hurwitz criterion can be applied. The first column of the Routh table are formed by the coefficients c 1,0 = 2 0 ∞ ν 2 , c 2,0 = σ f e − ∆ t and c 3,0 = 2 0 ∞ 1 − ν 2 . The later coefficient provides the classical condition for stability. In addition, we know that σ f e − must be a positive magnitude so as to address the stability in the algorithm updating. Table A1 lists the computational resources employed in the simulations. The operative system in which the simulations were carried out and runs on this hardware is Ubuntu Bionic Beaver.