Saturable Lorentz model for fully explicit three-dimensional modeling of nonlinear optics

Inclusion of the instantaneous Kerr nonlinearity in the FDTD framework leads to implicit equations that have to be solved iteratively. In principle, explicit integration can be achieved with the use of anharmonic oscillator equations, but it tends to be unstable and inappropriate for studying strong-field phenomena like laser filamentation. In this paper, we show that nonlinear susceptibility can be provided instead by a harmonic oscillator driven by a nonlinear force, chosen in a way to reproduce the polarization obtained from the solution of the quantum mechanical two level equations. The resulting saturable, nonlinearly-driven, harmonic oscillator model reproduces quantitatively the quantum mechanical solutions of harmonic generation in the under-resonant limit, up to the 9th harmonic. Finally, we demonstrate that fully explicit leapfrog integration of the saturable harmonic oscillator is stable, even for the intense laser fields that characterize laser filamentation and high harmonic generation.


Introduction
The finite-difference time-domain (FDTD) method is one of the most powerful methods available for solving rigorously Maxwell equations for light propagation in complex settings, such as for self-focusing and four-wave mixing in photonic devices [1]. However, there are still open fundamental questions and technical challenges associated with the inclusion of atomic-level contributions to that framework [2]. These questions and challenges are particularly relevant to modeling sub-wavelength nanophotonic devices that can have only a small number of atoms in one (e.g., a nanofilm) or more dimensions (e.g., a nanorod or a nanoparticle) and where continuum optics theory fails to describe light-matter interaction accurately. In this paper, we develop a three-dimensional (3D) nonlinear atomic polarization model that integrates well in the FDTD method, while retaining the most important atomic response features that usually manifest only through the solution of the quantum mechanical problem [3]. Ultimately, this model could be used in the microscopic particle-in-cell (MicPIC) framework [2,[4][5][6] to allow ab initio modeling of nonlinear optics, where optical materials are represented by ensembles of individual atoms.
Phenomenological and semi-classical models have been proposed, both in frequency and time domain, for the simulation of nonlinear light propagation. But in the FDTD framework, the inclusion of optical nonlinearities leads to implicit equations that have to be solved iteratively [7][8][9]. The major drawback with implicit schemes is that they rely on complex implementations and require matrix operations that tend to use a lot of computational resources. The nonlinear polarization model presented in this paper has the advantage of being fully explicit while effectively including more physics than models with, say, only an instantaneous Kerr contribution. Moreover, because it is based on an harmonic oscillator model, it can be combined seamlessly with the ordinary differential equations derived from the Sellmeier equation for a particular medium. Ultimately, the saturable oscillator model we propose can be generalized to deal with 3D vectorial propagation of intense laser pulses in anisotropic nonlinear media.
We recall that there are two possible contributions to the third-order nonlinearity: the electronic Kerr effect, which is considered here, and the vibrational Raman effect [3]. Besides allowing explicit integration of the Kerr nonlinearity in FDTD, our model also allows to describe the Kerr effect in the high-intensity non-perturbative limit of light-matter interaction, with potential application to modeling laser filamentation [10,11]. We emphasize that at these intensities, other effects such as ionization and plasma dynamics become important and have to be considered. Our intention in this paper is not to fully model laser filamentation but to show that the saturable oscillator model allows to treat Kerr nonlinearity at high intensities numerically in an explicit and stable manner, and in agreement with quantum mechanical results. The applicability of our model to the Raman nonlinearity will require further investigation.
The paper is divided as follows. In section 2, we show how explicit integration of nonlinear optical response in FDTD is possible. In section 3, we present the anharmonic oscillator model, an intuitive nonlinear extension to the Lorentz medium model. In section 4, we review the concepts of nonlinear optics in quantum mechanical two-level systems to explain the origin of the nonlinear response in terms of atomic transition saturation. In section 5, we present the saturable oscillator model and show that it behaves much like a quantum two-level system, up to the 9th harmonic, even in the strong electric fields that characterize laser filamentation. Finally, we draw conclusions in section 6.

Explicit nonlinear optics in the FDTD framework
Based on Yee's seminal work [12], modern FDTD methods [1,9] model light propagation in continuous polarizable media by solving the Maxwell's macroscopic curl equations [13,14]. For example, if we assume a nonmagnetic medium (B = µ 0 H) with a polarization density P, the FDTD technique solves using centered finite differences on numerical meshes, staggered in space and time, to achieve second-order accuracy. Finite difference operators and staggering are applied to all electric, magnetic, and electric polarization density components. In linear optics, the polarization vector is proportional to the electric field, i.e., P = ε 0 χ (1) E, where χ (1) is the linear susceptibility, a real constant (dispersion is neglected for simplicity). Then, applying the temporal finite difference operator to Eq. (1b), one can easily isolate E n+1 and express it in terms of known values of E and H at previous times only. Integration is then fully explicit. However, if the polarization is nonlinear in the electric field, for example, if P = ε 0 (χ (1) E + χ (3) |E| 2 E) like in the presence of instantaneous Kerr nonlinearity [3], it is no longer possible to isolate the individual field components. To find E n+1 , it is necessary to rely on iterative methods [7,8]. An alternative to the iterative implicit method is to use auxiliary oscillator-type differential equations (AODE). For example, a linear Lorentz medium is modeled by: where Ω is a coupling frequency, related to the static linear susceptibility χ (1) (ω = 0) = Ω 2 /ω 2 0 . Using leapfrogged finite differences: Then, Eq. (1b) can be written: which depends only on values in the past. From here, it is readily seen that fully explicit nonlinear optical material modeling can be achieved in the FDTD framework by replacing the linear driving force ε 0 Ω 2 E n in Eq. (4a) by a nonlinear force f (E n ).
The advantages of the AODE approach are the following. First, it naturally introduces a specific time dependence to the polarization and current densities to model dispersion. Second, with a second-order equation for the polarization density P, a current density J can be defined at the right temporal staggering points for the calculation of the electric field in the future, in a fully explicit manner. Third, because it is based on the harmonic oscillator, it integrates seamlessly to linear models based on the ordinary differential equations derived from the Sellmeier equation for a given medium. Finally, it can be generalized, with the right form of f (E), to take into account tensorial anisotropic response and nonlinear susceptibility of arbitrary order.
Later, in section 4, we define f (E) based on the solution for light propagation in a quantum mechanical two-level system. In section 5, we present the saturable oscillator model and show that quantitative agreement can be obtained with the quantum mechanical solution, up to the 9th harmonic. The observation of harmonics with order higher than 10 at high laser intensities (> 10 12 W/cm 2 ) suggests that f (E) could ultimately be tailored to model high-harmonic generation (HHG). We emphasize the fact that at such high laser intensities field ionization would have to be included to account for the relevant physical processes.

The anharmonic oscillator model
The optical response of dielectric media is usually well described in the context of the Lorentz model, where atoms are represented by classical harmonic oscillators [14]: Note that we neglected damping for simplicity, but it can readily be included by adding γdP/dt to the left-hand side of the left equation. A natural extension to this model is to add nonlinearity to the restoring force felt by the bound electrons, e.g., for a centrosymmetric material [3], where b is a constant proportional to χ (3) that modulates the strength of the third-order nonlinearity. The anharmonic oscillator model presented above is a perturbative approach, where the atomic potential is expanded in powers of the electronic displacement from non-driven equilibrium. It is valid when field excitation is not too strong. Its fully explicit integration in the FDTD/particle-in-cell (PIC) frameworks was demonstrated recently by Gordon, Helle, and Peñano [15]. However, as the authors reported, it tends to be unstable.
In the next two sections, we show how the Lorentz model of the atom can be extended to deal with nonlinear optics, not by adding inherently unstable nonlinear restoring force contributions, but by using the proper nonlinear driving force. The model we propose also supports fully explicit integration in the FDTD/PIC frameworks and can be seen as a drop-in replacement for the method proposed by Gordon, Helle, and Peñano [15].

Nonlinear saturation in a two-level system
The quantum mechanical two-level model provides a very good approximation to a more complete theory of nonlinear light-matter interaction [3,16]. Solution of the two-level optical Bloch equations in the under-resonant limit yields the polarization density of an ensemble of two-level atoms subjected to an electric field as [3] where χ 0 is a complex-valued ω-dependent weak-field susceptibility, ∆ = ω − ω 0 is the angular frequency relative to the resonance (transition) frequency ω 0 , τ is the induced dipole moment decay time, and E s is the line-center saturation field strength. If the denominator is expanded in a power series in terms of |E| 2 /|E s | 2 , it gives We emphasize that this expansion is valid in the weak field limit only (|E| 2 |E s | 2 ). The most important feature of Eq. (8)-the saturation of the susceptibility in the presence of strong fields-is equally lost in Eq. (9b) and the anharmonic oscillator model.
Finally, quantum mechanics suggests that a driving function of the form introduces a nonlinear polarization that expands, in the weak field limit, to the power series that is usually used to model centrosymmetric nonlinear media [3]. It is straightforward to see that the expansion of has even and odd orders in E and can potentially be used to model noncentrosymmetric nonlinear media. We emphasize that for an accurate description of optical saturation in atoms higher energy levels need to be taken into account. This will change the value of the saturation field strength E s , but not the structure of the expression. We will see in the next section that replacing E by f (E) in Eq. (3b) is actually a very good approximation to solving the quantum mechanical two-level equations, with the advantage that it exposes naturally the vectorial nature of the problem.

The saturable oscillator model
The saturable harmonic oscillator model is characterized by the following equation: Actually, it is the Lorentz model [see Eq. (6)] where the driving force is modulated by a saturation function 1/(1+|E| 2 /E 2 s ) whose expression is inspired by the quantum mechanical solution presented in section 4. The angular frequency Ω s and saturation field strength E s are fitting parameters that we will use below to match a particular medium response. Equation (12) can be solved with leapfrogged finite differences, as we have shown earlier at Eqs. (4a)-(4b), i.e., We have seen in section 2 that this particular form allows fully explicit integration of the electromagnetic field components with the FDTD method.
In the strong field limit (|E n | E s ), the driving force starts saturating, which introduces nonlinearities. In the weak field limit (|E n | E s ), after the power series expansion of the saturation term, one gets an anharmonic oscillator equation similar to Eq. (7).
Next, we test the saturable oscillator model numerically with parameters that characterize ultrashort pulse filamentation in air. In particular, in section 5.1, we take advantage of the fact that the leading term of the weak field expansion is the harmonic oscillator to obtain a value for Ω s from the Sellmeier equation for air. In turn, the critical power for self-focusing allows the definition of χ (3) and E s . Ultimately, in section 5.2, we compare the saturable oscillator model against those presented in the previous sections.

Parameters for nonlinear propagation in air
Laser filamentation is a spectacular manifestation of the richness of non-perturbative nonlinear optics. Besides technical applications for long-distance propagation of intense laser beams in the atmosphere, harmonic generation, and supercontinuum emission, there remain fundamental questions [10,11]. There is, in particular, the debated question about the contribution of nonlinear optical saturation to the stopping of self-focusing, which is usually attributed to plasma generation alone [10]. Recent quantum computations performed in the context of laser filamentation in gases display evidences that only relying on instantaneous higher-order Kerr nonlinearities while ignoring plasma generation is not possible at intensities of the order of 10 13 W/cm 2 and beyond [17][18][19]. However, experimental observations show that during filamentation the medium is only weakly ionized, i.e., most electrons remain bound. These bound electrons contribute to the Kerr nonlinearity, even for the high intensities that characterize the filament core, accepted to be in the 10 13 −10 14 W/cm 2 range [10,11]. The FDTD method, combined with the saturable oscillator model and proper ionization/plasma models [11], is a unique tool to study those questions in a setting that takes into account the temporal, full vectorial, and 3D nature of the phenomenon.
For the following numerical tests, we applied the saturable oscillator model to air, because of the importance of that medium in laser filamentation [10,11]. Our intention is not to fully model laser filamentation but to show that our simple model allows to treat Kerr nonlinearity at high intensities in a numerically stable manner and in agreement with quantum mechanical results. Based on the Sellmeier equation presented in [20], the linear response of air in the infrared range of the electromagnetic spectrum is accurately modeled by the following susceptibility: where Ω 1 = 6.411 × 10 14 rad/s ω 1 = 2.906 × 10 16 rad/s (15a) Ω 2 = 1.092 × 10 14 rad/s ω 2 = 1.427 × 10 16 rad/s (15b) From the critical power for self-focusing of a Gaussian beam in air at the wavelength of 800 nm (∼ 3.2 GW) [11], we inferred n 2 3 × 10 −23 m 2 /W and χ (3) 1.062 × 10 −25 m 2 /WΩ. Then we define the weak field nonlinear polarization density of air as: If, for simplicity, we compare only the first band (ω 1 ) near ω = 0 with the static solution of Eq. (12) we get immediately the following correspondence: With the parameters defined so far, it gives E s 6.77 × 10 10 V/m, which corresponds to the saturation intensity I s = E 2 s /2η 0 6 × 10 14 W/cm 2 , where η 0 is the characteristic impedance of vacuum.
As seen at Eq. (18), the saturable oscillator model gives the freedom to chose the linear and nonlinear refractive index independently, through the parameters Ω s , ω 0 , and E s . Although we applied the model to a single resonance of air, generalization to multiple frequencies is straightforward and done by using a specific saturable oscillator for each term of the Sellmeier equation associated with a given material.

Numerical results
We compared the saturable oscillator model against the anharmonic oscillator model [Eq. (7)], the quantum mechanical two-level model [3,16], and the instantaneous Kerr model [Eq. (16)]. In particular for the quantum mechanical two-level model, we integrated numerically Eqs. (6.5.6) and (6.5.8) found in [3], along with the dipole moment defined at Eq. (6.5.31) of the same reference. In all cases, the driving electric field was that of a Gaussian pulse with angular frequency ω l = 2πc/λ l (λ l = 800 nm) and pulse duration T = 10 fs. The field amplitude is defined in terms of the laser intensity I = E 2 0 /2η 0 . For all models we used the parameters defined in section 5.1.
In particular, we compared the harmonic power spectrum generated by the time-dependent nonlinear polarization density P(t) by inserting the Fourier transform P(ω) into the Larmor radiation formula [13] All curves were normalized by the value at the peak of the linear response (ω/ω l = 1) of the instantaneous Kerr model. Individual model parameters were optimized to reproduce the instantaneous Kerr curve in the perturbation limit where only the first and third harmonics are present, shown in figure 1. In this regime all models agree very well, as expected. At the beginning of the laser filamentation process, a Gaussian beam propagating in a Kerr medium with a power above the critical power will start self-focusing under the action of the  Kerr two-level saturable anharmonic Fig. 2. Comparison between all four models for I = 2.5 × 10 11 W/cm 2 (see also caption of figure 1). In insets, a zoomed view of the curves show that the saturable oscillator reproduces fairly well the phase relationship between frequency components of the quantum mechanical solution, while the anharmonic oscillator does not. Also, the anharmonic oscillator tends to overestimate the 7th and 9th harmonics by a few orders of magnitude.
intensity-dependent refractive index [3]. During that phase, the peak intensity increases rapidly to reach 10 11 − 10 12 W/cm 2 [10,11]. Results obtained in this regime are shown in figure 2. From the numerical results it appears that the 5th harmonic is already stronger than the 3rd harmonic in the lower intensity example of figure 1, suggesting that it could have a similar relative influence. Overall, the two-level, saturable oscillator, and anharmonic oscillator models agree well, although the anharmonic oscillator overestimates the 7th and 9th harmonics by a few orders of magnitude and fails at reproducing the dips of destructive interference predicted by both the two-level and saturable oscillator solutions (some are also predicted by the Kerr model, but are not shown).
In the widely accepted picture of laser filamentation in air, self-focusing is slowed down and stopped by a defocusing plasma formed due to multiphoton ionization of O 2 molecules, when the intensity reaches 10 12 − 10 13 W/cm 2 [10,11]. The so-called filament emerges from the Kerr two-level saturable anharmonic Fig. 3. Comparison between all four models for I = 5 × 10 13 W/cm 2 , a laser intensity comparable to the peak intensity inside a light filament in air [10,11] (see also caption of figure 1). Whereas the Kerr, two-level, and saturable oscillator models agree, the aharmonic oscillator fails at reproducing the dynamics of the two-level system. With an intensity just a bit higher, at I = 6.7 × 10 13 W/cm 2 , the numerical integration of the anharmonic oscillator does not converge anymore. balance between nonlinear self-focusing and plasma defocusing. Peak intensity is then in the 10 13 − 10 14 W/cm 2 range. Atomic polarization in that regime is shown in figure 3, for the four models. Let aside the fact that the Kerr model gives only first and third order contributions, the Kerr, two-level, and saturable oscillator models agree very well, whereas the anharmonic oscillator model is completely off by several orders of magnitude at the 5th, 7th, and 9th harmonics.
One can see a manifestation of the instability reported by Gordon, Helle, and Peñano [15] in figure 3(b), where the 9th harmonic peak is overshooting the 7th and the 5th. Ultimately, at ∼ 6.7 × 10 13 W/cm 2 the 3rd harmonic overshoots the linear peak and the numerical integration of the anharmonic oscillator does not converge anymore, no matter how small the time step is. We emphasize that the anharmonic oscillator instability will show at much lower intensities when modeling solid-density materials because the nonlinear susceptibility is much higher. An expanded view of figure 3(a) is given in figure 4. It is readily observed that the agreement between the quantum mechanical two-level model and the classical saturable oscillator model is almost perfect up to the 9th harmonic, but not beyond. Finally, we emphasize that HHG in gases is commonly performed in the 10 14 − 10 15 W/cm 2 intensity range, where the atomic polarization model presented here needs to be complemented by the ionization-recollision process [21,22].

Conclusions
In this paper, we have demonstrated that the saturable oscillator model behaves effectively much like a quantum mechanical two-level system. The advantage of the saturable oscillator over the two-level model is the possibility to extend-in a simple and intuitive way-the Sellmeierharmonic oscillator formalism to include damping, restoring force, and driving force coupling in a tensorial form to allow effective modeling of the anisotropic nonlinear response and excitation decay of specific crystalline structures and molecules. We have shown that numerical integration can be fully explicit, allowing straightforward integration to the FDTD, PIC, and MicPIC frameworks with, ultimately, the possibility to model nonlinear light propagation with atomic scale resolution. In particular, we have shown that the saturable oscillator model gives accurate results up to the 9th harmonic, for laser intensities relevant to the process of laser filamentation in air.