## Abstract

The Lorentz model provides a simple and intuitive approach to linear optical dispersion. In this paper, we present a nonlinear formulation of the Lorentz model based on an analytical solution of a quantum-mechanical two-level model atom in the under-resonant limit and show how multiple nonlinear Lorentz equations can be combined to describe nonlinear optical dispersion in both time and frequency domain nanophotonics simulations.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The development of current and next-generation nanophotonics technologies heavily relies on computer modelling (see, e.g., [1]). As rigorous quantum *N*-body calculations of the polarization in atoms, molecules, and solids are numerically too expensive to be used in full-scale electromagnetic simulations, there is a pressing need for advanced material models that can describe the underlying quantum physics with sufficient accuracy and computational efficiency. In particular, there has been significant progress on fully spatio-temporal *ab initio* nanophotonics modeling with multi-level atomic systems, including loss-free plasmonics, SPASERS, optical limiters, plasmon-enhanced nonlinear optics, etc. (see, e.g., Refs. [2–5]). In this paper, we present a nonlinear generalization of the Lorentz dispersion model.

The Lorentz model has been used extensively to describe optical dispersion, but it does not deal with nonlinearity. The standard approach is to use the anharmonic (Duffing) oscillator equation instead [6, 7]. In this case, nonlinearity is introduced *via* nonlinear restoring forces, which results in numerical instabilities in time-domain simulations [7–9] and makes establishing the correspondence between time-domain and frequency-domain expressions rather cumbersome (see, e.g., [6]). In previous works, we proposed nonlinear alternatives to the Lorentz equation that are robust to numerical integration in the finite-difference time-domain (FDTD) framework [8,9]. However, the wavelength-dependence of the nonlinear susceptibility was not addressed. This work improves upon Refs. [8,9] by showing how to use the perturbative form of these nonlinear Lorentz equations to model *nonlinear* optical dispersion in both the time and frequency domains for applications in nanophotonics simulations.

The paper is structured as follows. In Sec. 2, we show how the nonlinear Lorentz equation can be derived from the quantum mechanical two-level atom model. In Sec. 3, we present a detailed spectral analysis. In Sec. 4, we propose an intuitive decomposition that generalizes the Sellmeier approach [10,11] to nonlinear dispersion. Finally, we conclude in Sec. 5.

## 2. The nonlinear Lorentz equation from quantum mechanics

The quantum mechanical two-level atom model provides a good approximation to a more complete theory of optical polarizability (see, e.g., [12]). In this context, optical transitions are conveniently described by the optical Bloch equations [13–15]. To deal with both centrosymmetric and non-centrosymmetric media, we generalized the usual formalism by including permanent dipole moments for the ground state and excited state. Following the procedure found in Sec. 6.4.1 of Ref. [13], it is possible to show that the temporal evolution of the expectation value of the dipole moment (*p*) of the two-level transition is then characterized by the following equations:

*E*is the driving electric field,

*w*is the population inversion parameter with equilibrium value

*w*, $\kappa =2{\omega}_{0}{\mu}_{0}^{2}/\hslash $ is an atom-field coupling parameter,

^{eq}*ħω*

_{0}is the transition energy between the two levels,

*T*

_{1}is the population relaxation time, and

*T*

_{2}is the atomic dephasing time. The atomic polarizability is parameterized by the atomic dipole transition constant

*μ*

_{0}, the difference between the permanent dipole moments of the ground (

*μ*) and excited (

_{g}*μ*) states Δ

_{e}*μ*=

*μ*−

_{g}*μ*and the average value

_{e}*μ̄*= (

*μ*+

_{g}*μ*)/2.

_{e}We will now introduce an approximation for the population inversion parameter *w* to combine Eqs. (1) and (2) into a single equation. We assume for now that the laser field is weak and the laser frequency *ω _{L}* is sufficiently detuned from the transition frequency

*ω*

_{0}so that the electron population remains in the ground state at all time, i.e.,

*w*= −1. Then, Eqs. (1) and (2) simplify to the single equation that follows

*μ*=

_{g}*μ*= 0. For strong laser fields,

_{e}*w*= −1 is no longer true and we have to seek for a time-varying solution for the population inversion parameter. Next, we will show that it is effectively possible to write Eqs. (1) and (2) in a Lorentz-like form where

*f*(

*E*) ≡

*f*[

*E*(

*t*)] is any function of

*E*(

*t*),

*a priori*. The Lorentz equation is a special case where

*f*(

*E*) is a linear function of the electric field [see Eq. (4)]. We will look for an approximate solution for

*f*(

*E*) and its Taylor series around

*E*= 0 to get a form like

*f*(

*E*) ≃

*a*

_{0}+

*a*

_{1}

*E*+

*a*

_{2}

*E*

^{2}+

*a*

_{3}

*E*

^{3}+ . . . that is characteristic of

*perturbative*nonlinear optics.

For dielectrics, the band gap energy is typically larger than 4 eV, whereas*T*_{1} and*T*_{2} are associated with slow relaxation processes, with typical values larger than 1 ns and 1 ps, respectively [13,16]. Thus for femtosecond pulses in dielectrics, we can assume that $2{\omega}_{L}/{T}_{2}\ll {\omega}_{L}^{2}<{\omega}_{0}^{2}$, which leads to this approximate solution for the atomic polarization

*a*= Δ

*μ*/

*ħω*

_{0}. This approximation is a necessary step to find an analytical solution for

*w*. As done above to derive Eq. (4), the resulting solution for

*w*will be reinserted in Eq. (1) to find an expression for a driving function

*f*(

*E*) and get a final equation similar to Eq. (5).

Inserting Eq. (6) into Eq. (2) gives

*w*suggests a fast contribution oscillating on the laser-cycle time scale and a slow envelope modulation evolving with

*T*

_{1}≳ 1 ns. In most situations involving laser pulses, pulse duration

*t*is orders of magnitude less than

_{p}*T*

_{1}and the second term on the left-hand side of Eq. (7) can be safely neglected. After taking the time derivative on the right-hand side and doing basic algebra, it is straightforward to show that

*dw*/

*dt*to the left-hand side of the equation and all terms with

*dE*/

*dt*to the right-hand side we get Equation (9) is solved by integrating both sides of the equation, which leads to the following solution:

*C*is an integration constant. Imposing that the electron population is in the ground state (

*w*= −1) when

*E*= 0 gives

*C*= −1.

Equation (10) is an approximate, time-varying solution for *w* valid in a regime that is known as the under-resonant (*ω _{L}* ≪

*ω*

_{0}) adiabatic-following (

*t*≪

_{p}*T*

_{1}

*, T*

_{2}) limit (see, e.g., Ref. [13] in Sec. 6.4.2 for further details). As done above to obtain Eqs. (3) and (4) with

*w*= −1, the approximate solution for

*w*given by Eq. (10) is reinserted into Eq. (1), which gives

*f*= Δ

*μ*/2, whereas

*a*and

*b*are defined below Eqs. (6) and (8), respectively.

For non-centrosymmetric media where *a* ≠ 0, the term ${\omega}_{0}^{2}(1+aE)p$in Eq. (11) needs a special treatment. For the laser intensities typically involved in perturbative nonlinear optics, we could assume that ${\omega}_{0}^{2}(1+aE)p\simeq {w}_{0}^{2}p$ (e.g., for *I* ≲ 10^{9} W/cm^{2}, *aE* ≲ 10^{−3}). An alternative approach is to transfer ${\omega}_{0}^{2}aEp$ to the right-hand side of the equation and use Eq. (6) to approximate this source-term contribution. Both options give similar results, when the right-hand side of Eq. (11) is expanded in a Taylor series around *E* = 0. For the second case, it gives

A macroscopic formulation of Eq. (12) is obtained by multiplying both sides of the equation by the atomic density *N*. Noting that the first term in the curly brackets on the right-hand side is the permanent dipole moment of the ground state, it makes sense to define the *induced* macroscopic polarization density *P* = *N*(*p* − *μ _{g}*) such that

*γ*= 2/

*T*

_{2},

*χ̄*

^{(1)}=

*Na*

_{1}/

*∊*

_{0},

*χ̄*

^{(2)}=

*Na*

_{2}/

*∊*

_{0}, and

*χ̄*

^{(3)}=

*Na*

_{3}/

*∊*

_{0}. It is observed that in the weak-field limit (

*E*→ 0), this nonlinear Lorentz equation [Eq. (13)] simplifies to the usual, linear Lorentz dispersion model.

Summing up, in this section we derived a nonlinear form of the Lorentz dispersion model from the quantum mechanical equations for a two-level atom model. The microscopic [Eq. (12)] and macroscopic [Eq. (13)] formulations of the nonlinear Lorentz model provide equations that will give results comparable to that from to the quantum mechanical equations [e.g., Eqs. (1) and (2)] when the photon energy is much smaller than the transition energy (i.e., *ħω _{L}* ≪

*ħω*

_{0}) and when no relaxation or Rabi oscillations occur during the laser pulse (i.e.,

*t*≪

_{p}*T*

_{1},

*T*

_{2}and 2

*μ*

_{0}

*E*

_{0}/

*ħ*≪ 2

*π*/

*t*, where

_{p}*E*

_{0}is the electric field peak value). These conditions are typically met for moderately intense femtosecond pulses in dielectrics. In the next section, we test the validity of the nonlinear Lorentz model by comparing numerical solutions obtained with both Eq. (13) and the Schrödinger equation.

## 3. Spectral analysis of the nonlinear Lorentz equation

To obtain the spectral signature of the nonlinear Lorentz equation, we integrated numerically Eq. (13) with the leapfrog method (see, e.g., Sec. 2.1 in Ref. [9] for details) for two reference media with optical properties comparable to air (*n*_{0} = 1.0003 and *n*_{2} = 5 × 10^{−19} cm^{2}/W) and BK-7 glass (*n*_{0} = 1.5 and *n*_{2} = 3.4 × 10^{−16} cm^{2}/W) [13]. The driving electric field is defined as *E*(*t*) = *E*_{0} cos(*ω _{L}t*) exp(−

*t*

^{2}/

*T*

^{2}), where

*T*= 10 fs and ${E}_{0}=\sqrt{2{\eta}_{0}{I}_{0}}$, with

*I*

_{0}= 10

^{6}W/cm

^{2}. We set

*γ*= 0 and the oscillator frequency arbitrarily high (

*ω*

_{0}= 3.14 × 10

^{16}rad/s) to fulfill the under-resonant criteria (

*ω*≪

_{L}*ω*

_{0}) for optical frequencies in the infrared range. We scanned the laser wavelength

*λ*= 2

_{L}*πc*/

*ω*(keeping all other parameters constant) and computed the Fourier transform of the polarization density

_{L}*P*(

*t*) to extract the amplitude of the linear (|

*P*(

*ω*)|) and nonlinear (|

_{L}*P*(3

*ω*)|) scattering signals for each

_{L}*λ*. For comparison, we applied the same methodology to the solutions obtained by direct integration of the two-level atom model in the Schrödinger representation (for details, see Sec. 6.5.1 in [13] and Appendix D.2 in [9]). Results are presented in Fig. 1.

_{L}To support our analysis, we recall that the two-level atom model for a centrosymmetric media (damping neglected) has a perturbative solution of the following form [6,13]:

*χ̄*

^{(1)}and

*χ̄*

^{(3)}are the first-order and third-order static susceptibility parameters, respectively. The dispersion curve for the third harmonic is thus expected to show two distinct resonance peaks: one at

*ω*=

_{L}*ω*

_{0}and another at

*ω*=

_{L}*ω*

_{0}/3 (for a second-order material, the second peak would appear at

*ω*=

_{L}*ω*

_{0}/2 instead). Looking at Fig. 1(a), it is observed that the nonlinear Lorentz model reproduces correctly linear scattering [

*P*

^{(1)}(

*ω*)], as expected from the usual Lorentz model. However, looking at Fig. 1(b), we see that the second nonlinear scattering peak is absent from the solutions obtained with Eq. (13). Nevertheless, this is an advantage for fitting nonlinear dispersion data, as we will see in Sec. 4.

_{L}To get a better insight into the spectral signature of Eq. (13), we split the total polarization density *P* into first and third order contributions and took their Fourier transform. With *E* = *E*_{0} cos(*ω _{L}t*) exp(−

*t*

^{2}/

*T*

^{2}),

*f̂*

_{(ω)}= FT{

*f*(

*x*)}, and the two following identities: FT{

*f*(

*x*) cos(

*ax*)} = [

*f̂*(

*ω*−

*a*) +

*f̂*(

*ω*+

*a*)]/2 and cos

^{3}(

*θ*) = [3 cos(

*θ*) + cos(3

*θ*)]/4, it is possible to show that

*ω*and ±3

_{L}*ω*, associated with self-phase modulation and third harmonic generation, respectively.

_{L}*Ĝ*(

_{n}*ω*) is a narrow function centered at

*ω*= 0 such that

*Ĝ*(

_{n}*ω*≠ 0) ≃ 0 (it is a Dirac delta in the limit

*T*→ ∞):

*ω*=

*ω*) and third-harmonic (

_{L}*ω*= 3

*ω*) scattering for

_{L}*T*≫ 1/

*ω*and

_{L}*γ*= 0:

*P̃*

^{(1)}(

*ω*) is the perturbative result of Eq. (14) times a constant, confirming the agreement observed in Fig. 1(a). On the other hand

_{L}*P̃*

^{(3)}(3

*ω*), as seen in Fig. 1(b), has only a single peak at

_{L}*ω*=

_{L}*ω*

_{0}/3 [because

*Ĝ*

_{3}(3

*ω*±

_{L}*ω*) ≃ 0], thus missing to correctly describe the third-harmonic spectra of the quantum solution beyond

_{L}*ω*

_{0}/3. Nevertheless, having a single pole in

*P̃*

^{(3)}(3

*ω*) is actually an advantage to fit nonlinear dispersion data. We develop this idea below.

_{L}## 4. Modelling nonlinear optical dispersion with the nonlinear Lorentz model

The Lorentz dispersion model is used extensively to model linear scattering. Equally, the nonlinear formulation [Eq. (13)] can be used to fit nonlinear scattering data. Next, we present a formal framework that we apply to examples in Secs. 4.1 and 4.2.

The Fourier transform FT { } of Eq. (13) gives

*ξ*is associated with an independent oscillator equation and susceptibility parameter

*χ̄*

^{(ξ)}, i.e.,

*ξ*= 1) is the usual Lorentz equation, whereas higher-order equations (

*ξ*> 1) correspond to the nonlinear contribution of order

*ξ*.

By summing the contributions from an ensemble of oscillators like Eq. (20)—labeled by the subscript *k*—it is possible to get an effective *ξ*th-order dispersion equation as

*ξ*= 1 and

*γ*= 0, Eq. (21) leads to the Sellmeier dispersion formula for the linear refractive index $n(\omega )=\sqrt{1+{{\sum}_{k}\chi}_{k}^{(1)}(\omega )}$ [10, 11]. The nonlinear Lorentz model thus provides a natural and intuitive generalization to nonlinear optics of the traditional modelling approaches for linear optics in transparent media, while dealing phenomenologically with the underlying optical processes.

_{k}#### 4.1. Example 1: fitting the quantum two-level solution

In Sec. 3, we have seen that third-order scattering predicted by the quantum two-level model is characterized by two resonance peaks, whereas the nonlinear Lorentz model has only one [see Fig. 1(b)]. The analysis provided in Sec. 4 suggests that we can mimic the quantum behaviour by adding an extra third-order (*ξ* = 3) nonlinear Lorentz equations [Eq. (20)] with a resonance frequency at 3*ω*_{0}. With optimized third-order susceptibility (${\overline{\chi}}_{k}^{(3)}$) and damping (*γ _{k}*) parameters, this basic two-oscillator fit gives an impressive agreement with the quantum solution, with an average error around 3% (see Fig. 2). The large error around

*ω*=

_{L}*ω*

_{0}could be reduced by using more nonlinear Lorentz equations for the fit (see Example 2 below).

#### 4.2. Example 2: fitting nonlinear dispersion data

Using four third-order (*ξ* = 3) frequency-domain nonlinear-Lorentz equations [see Eq. (21)], we fitted the data and best fit proposed by Torruellas *et al.* for third-order dispersion *χ*^{(3)}(3*ω*) in poly(3-decylthiophene) reported in [17]. The individual oscillator parameters were optimized until the maximum error went below 5% and the average error below 1% (see Table 1). The fit results shown in Fig. 3 clearly demonstrate that the nonlinear Lorentz model can reproduce with high accuracy the model originally proposed by Torruellas *et al.* [17]. We emphasize that Torruellas *et al*.’s fit model is limited to the spectral domain.

## 5. Conclusion

We have introduced a theoretical framework that extends the scope of the Lorentz linear dispersion model for transparent media to the perturbative nonlinear regime. We derived the underlying equations starting from the quantum mechanical equations of the two-level atom model, in the adiabatic-following, under-resonant approximation. We provided temporal and spectral formulations of the nonlinear Lorentz model, along with compelling evidence that the proposed framework can be used for modelling nonlinear dispersion of arbitrary complexity in both time-domain and frequency-domain nanophotonics simulations.

## Funding

Natural Sciences and Engineering Research Council of Canada (NSERC) through the College and Community Innovation Program - Innovation Enhancement Grants (CCIPE 517932-17); Fonds de recherche du Québec - Nature et technologies (FRQNT) through the Programme de recherche pour les chercheurs et les chercheuses de collège (2019-CO-254385); Deutsche Forschungsgemeinschaft via SPP1840 and a Heisenberg fellowship (ID FE 1120/4-1).

## References

**1. **B. Gallinet, J. Butet, and O. J. F. Martin, “Numerical methods for nanophotonics: standard problems and future challenges,” Laser Photonics Rev. **9**, 577–603 (2015). [CrossRef]

**2. **X. Jiang and C. M. Soukoulis, “Time dependent theory for random lasers,” Phys. Rev. Lett. **85**, 70–73 (2000). [CrossRef]

**3. **S. Wuestner, A. Pusch, K. L. Tsakmakidis, J. M. Hamm, and O. Hess, “Overcoming losses with gain in a negative refractive index metamaterial,” Phys. Rev. Lett. **105**, 127401 (2010). [CrossRef]

**4. **J. Lee, M. Tymchenko, C. Argyropoulos, P.-Y. Chen, F. Lu, F. Demmerle, G. Boehm, M.-C. Amann, A. Alù, and M. A. Belkin, “Giant nonlinear response from plasmonic metasurfaces coupled to intersubband transitions,” Nature **511**, 65–69 (2014). [CrossRef] [PubMed]

**5. **S. I. Azzam and A. V. Kildishev, “Time-domain dynamics of saturation of absorption using multilevel atomic systems,” Opt. Mater. Express **8**, 3829–3834 (2018). [CrossRef]

**6. **M. Scalora, M. A. Vincenti, D. de Ceglia, C. M. Cojocaru, M. Grande, and J. W. Haus, “Nonlinear Duffing oscillator model for third harmonic generation,” J. Opt. Soc. Am. B **32**, 2129 (2015). [CrossRef]

**7. **D. Gordon, M. Helle, and J. Peñano, “Fully explicit nonlinear optics model in a particle-in-cell framework,” J. Comput. Phys. **250**, 388–402 (2013). [CrossRef]

**8. **C. Varin, G. Bart, R. Emms, and T. Brabec, “Saturable Lorentz model for fully explicit three-dimensional modeling of nonlinear optics,” Opt. Express **23**, 2686–2695 (2015). [CrossRef]

**9. **C. Varin, R. Emms, G. Bart, T. Fennel, and T. Brabec, “Explicit formulation of second and third order optical nonlinearity in the FDTD framework,” Comput. Phys. Commun. **222**, 70–83 (2018). [CrossRef]

**10. **W. Sellmeier, “Zur Erklärung der abnormen Farbenfolge im Spectrum einiger Substanzen,” Ann. Phys. Chem. **219**, 272–282 (1871). [CrossRef]

**11. **B. Tatian, “Fitting refractive-index data with the Sellmeier dispersion formula,” Appl. Opt. **23**, 4477 (1984). [CrossRef]

**12. **H. Bethe and E. Salpeter, *Quantum mechanics of one- and two-electron atoms* (Springer, 1957). [CrossRef]

**13. **R. Boyd, *Nonlinear Optics*, 3rd ed. (Academic Press, 2008).

**14. **M. Scully and M. Zubairy, *Quantum Optics* (Cambridge University Press, 1997). [CrossRef]

**15. **L. Allen and J. Eberly, *Optical Resonance and Two-level Atoms* (Dover, 1975).

**16. **A. Siegman, *Lasers* (University Science Books, 1986).

**17. **W. Torruellas, D. Neher, R. Zanoni, G. Stegeman, F. Kajzar, and M. Leclerc, “Dispersion measurements of the third-order nonlinear susceptibility of polythiophene thin films,” Chem. Phys. Lett. **175**, 11–16 (1990). [CrossRef]