Pulse propagation in the slow and stopped light regime

Electromagnetic pulse propagation in the slow light regime and near a zero group velocity point is relevant to a plethora of potential applications, and has analogies in numerous other wave systems. Unfortunately, the standard frequency-based formulation for pulse propagation is unsuitable for describing the dynamics in such regimes, due to the divergence of the dispersion coefficients. Moreover, in the presence of absorption, it is not clear how to interpret the propagation dynamics due to the drastic change induced by absorption upon the dispersion curves. As a remedy, we present an alternative momentum-based formulation, which is rapidly converging in these regimes, and naturally suitable for lossy and nonlinear media. It is specialized to a waveguide geometry which provides a significant simplification with respect to existing momentum-based schemes. Doing so, we provide a somewhat alternative, yet intuitive picture of the seeming enhanced absorption and nonlinear response in these regimes, and show that light-matter interactions are not enhanced in the slow/stopped light regimes. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (190.5530) Pulse propagation and temporal solitons; (320.7130) Ultrafast processes in condensed matter, including semiconductors. References and links 1. J. R. Pierce, Traveling-Wave Tubes (D. Van Nostrand Company, Inc., 1950). 2. M. Fleischhauer, A. Imamoglu, and J. P. Marangos, “Electromagnetically induced transparency: optics in coherent media,” Rev. Mod. Phys. 77, 633–673 (2005). 3. P. Siddons, “Light propagation through atomic vapours,” J. Phys. B: At. Mol. Opt. Phys 47, 093001 (2014). 4. J. B. Khurgin, “Slow light in various media: a tutorial,” Adv. Opt. Photonics 2, 287–318 (2010). 5. R. W. Boyd, “Material slow light and structural slow light: similarities and differences for nonlinear optics,” J. Opt. Soc. Am. B 28, A38–A44 (2010). 6. T. Baba, “Slow light in photonic crystals,” Nat. Photonics 2, 465–473 (2008). 7. T. F. Krauss, “Slow light in photonic crystal waveguides,” J. Phys. D 40, 2666–2670 (2007). 8. G. Heinze, C. Hubrich, and T. Halfmann, “Stopped light and image storage by electromagnetically induced transparency up to the regime of one minute,” Phys. Rev. Lett. 111, 033601 (2013). 9. R. Won, “Quantum memory: Extended storage times,” Nat. Photonics 7, 677 (2013). 10. G. Brennen, E. Giacobino, and C. Simon, “Focus on quantum memory,” New J. Phys. 17, 050201 (2007). 11. T. Tanabe, M. Notomi, E. Kuramochi, A. Shinya, and H. Taniyama, “Trapping and delaying photons for one nanosecond in an ultrasmall high-Q photonic-crystal nanocavity,” Nat. Photonics 1, 49–52 (2007). 12. M. Notomi, E. Kuramochi, and T. Tanabe, “Large-scale arrays of ultrahigh-q coupled nanocavities,” Nat. Photonics 2, 741–747 (2008). 13. H. Takesue, N. Matsuda, E. Kuramochi, W. J. Munro, and M. Notomi, “An on-chip coupled resonator optical waveguide single-photon buffer,” Nat. Commun. 4, 2725 (2013). 14. C. Monat, B. Corcoran, D. Pudo, M. E. Heidari, C. Grillet, M. D. Pelusi, D. J. Moss, B. J. Eggleton, T. P. White, L. O’Faolain, and T. F. Krauss, “Slow light enhanced nonlinear optics in silicon photonic crystal waveguides,” IEEE J. Sel. Top. Quantum Electron. 16, 344–356 (2010). 15. S. Lavdas and N. Panoiu, “Theory of pulsed four-wave mixing in one-dimensional silicon photonic crystal slab waveguides,” Phys. Rev. B 93, 115435 (2016). 16. T. Kamalakis, “Derivation of coupled envelope propagation equations in silicon-based photonic crystal waveguides,” J. Opt. Soc. Am. B 33, 2339–2349 (2016). 17. X. Checoury, Z. Han, and P. Boucaud, “Stimulated Raman scattering in silicon photonic crystal waveguides under continuous excitation,” Phys. Rev. B 82, 041308 (2010). 18. J. F. McMillan, M. Yu, D.-L. Kwong, and C. W. Wong, “Observation of spontaneous raman scattering in silicon slow-light photonic crystal waveguides,” Appl. Phys. Lett. 93, 251105 (2008). Vol. 26, No. 15 | 23 Jul 2018 | OPTICS EXPRESS 19294 #319900 https://doi.org/10.1364/OE.26.019294 Journal © 2018 Received 18 Jan 2018; revised 15 Mar 2018; accepted 15 Mar 2018; published 17 Jul 2018 19. B. Corcoran, C. Monat, C. Grillet, D. J. Moss, B. J. Eggleton, T. P. White, L. O’Faolain, and T. F. Krauss, “Green light emission in silicon through slow-light enhanced third-harmonic generation in photonic-crystal waveguides,” Nat. Photonics 3, 206–210 (2009). 20. C. Monat, B. Corcoran, M. Ebnali-Heidari, C. Grillet, B. J. Eggleton, T. P. White, L. O’Faolain, and T. F. Krauss, “Slow light enhancement of nonlinear effects in silicon engineered photonic crystal waveguides,” Opt. Express 17, 2944–2953 (2009). 21. Y. Hamachi, S. Kubo, and T. Baba, “Slow light with low dispersion and nonlinear enhancement in a lattice-shifted photonic crystal waveguide,” Opt. Lett. 34, 1072–1074 (2009). 22. Y. Okawachi, M. S. Bigelow, J. E. Sharping, Z. Zhu, A. Schweinsberg, D. J. Gauthier, R. W. Boyd, and A. L. Gaeta, “Tunable all-optical delays via brillouin slow light in an optical fiber,” Phys. Rev. Lett. 94, 153902 (2005). 23. A. Zadok, A. Eyal, and M. Tur, “Stimulated brillouin scattering slow light in optical fibers,” Appl. Opt. 50, E38–E49 (2011). 24. M. Afzelius, N. Gisin, and H. de Riedmatten, “Quantum memory for photons,” Phys. Today 68, 42–47 (2015). 25. H. de Riedmatten, “A long-term memory for light,” Physics 80, 6 (2013). 26. F. Blatt, L. S. Simeonov, T. Halfmann, and T. Peters, “Stationary light pulses and narrowband light storage in a laser-cooled ensemble loaded into a hollow-core fiber,” Phys. Rev. A 94, 043833 (2016). 27. K. V. Rajitha, T. N. Dey, J. Evers, and M. Kiffner, “Pulse splitting in light propagation through N-type atomic media due to an interplay of Kerr nonlinearity and group-velocity dispersion,” Phys. Rev. A 92, 023840 (2015). 28. G. Grigoryan, V. Chaltykyan, E. Gazazyan, O. Tikhova, and V. Paturyan, “Pulse propagation, population transfer, and light storage in five-level media,” Phys. Rev. A 91, 023802 (2015). 29. Y. Chen, Z. Chen, and G. Huang, “Storage and retrieval of vector optical solitons via double electromagnetically induced transparency,” Phys. Rev. A 91, 023820 (2015). 30. L. V. Hau, S. E. Harris, Z. Dutton, and C. H. Behroozi, “Light speed reduction to 17 metres per second in an ultracold atomic gas,” Nature 397, 594–598 (1999). 31. M. M. Kash, V. A. Sautenkov, A. S. Zibrov, L. Hollberg, G. R. Welch, M. D. Lukin, Y. Rostovtsev, E. S. Fry, and M. O. Scully, “Ultraslow group velocity and enhanced nonlinear optical effects in a coherently driven hot atomic gas,” Phys. Rev. Lett. 82, 5229–5232 (1999). 32. P. C. Kuan, C. Huang, W. S. Chan, S. Kosen, and S. Y. Lan, “Large fizeau’s light-dragging effect in a moving electromagnetically induced transparent medium,” Nat. Commun. 7, 13030 (2016). 33. H. Zhang, M. Sabooni, L. Rippe, C. Kim, S. Kröll, L. V. Wang, and P. R. Hemmer, “Slow light for deep tissue imaging with ultrasound modulation,” Appl. Phys. Lett. 100, 131102 (2012). 34. M. S. Shahriar, G. S. Pati, R. Tripathi, V. Gopal, M. Messall, and K. Salit, “Ultrahigh enhancement in absolute and relative rotation sensing using fast and slow light,” Phys. Rev. A 75, 053807 (2007). 35. S. Wuestner, T. W. Pickering, J. M. Hamm, A. F. Page, A. Pusch, and O. Hess, “Ultrafast dynamics of nanoplasmonic stopped-light lasing,” Faraday Discuss. 178, 307 (2015). 36. A. F. Page, “Phd thesis,” Ph.D. thesis, Imperial College London (2015). 37. K. L. Tsakmakidis, T. W. Pickering, J. M. Hamm, A. F. Page, and O. Hess, “Completely stopped and dispersionless light in plasmonic waveguides,” Phys. Rev. Lett. 112, 167401 (2014). 38. K. L. Tsakmakidis, T. W. Pickering, J. M. Hamm, A. F. Page, and O. Hess, “Cavity-free plasmonic nanolasing enabled by dispersionless stopped light,” Nat. Commun. 5, 4972–4979 (2014). 39. N. Sun, J. Chen, and D. Tang, “Stopping light in an active medium,” The Eur. Phys. J. D 69, 219 (2015). 40. J. I. L. Chen, G. von Freymann, S. Y. Choi, V. Kitaev, , and G. A. Ozin, “Amplified photochemistry with slow photons,” Adv. Mater. 18, 1915–1919 (2006). 41. J. I. L. Chen, G. von Freymann, S. Y. Choi, V. Kitaev, , and G. A. Ozin, “Slow photons in the fast lane in chemistry,” J. Mater. Chem. 18, 369–373 (2008). 42. U. Leonhardt, “A laboratory analogue of the event horizon using slow light in an atomic medium,” Nature. 415, 406–409 (2002). 43. M. Fink, “Acoustic time-reversal mirrors: Imaging of complex media with acoustic and seismic waves,” Top. Appl. Phys. 84, 17–43 (2002). 44. A. Hayrapetyan, K. Grigoryan, R. Petrosyan, and S. Fritzsche, “Propagation of sound waves through a spatially homogeneous but smoothly timedependent medium,” Ann. Phys. 333, 47–65 (2013). 45. H. Cui, W. Lin, H. Zhang, and X. Wang, “Backward waves with double zero-group-velocity points in a liquid-filled pipe,” The J. Acoust. Soc. Am. 139, 1179–1194 (2016). 46. A. Chumak, A. Karenowska, A. Serga, and B. Hillebrands, in Topics in Applied Physics. Magnonics: from fundamentals to applications (Springer, 2013). 47. V. Bacot, M. Labousse, A. Eddi, M. Fink, and E. Fort, “Revisiting time reversal and holography with spacetime transformations,” Nat. Phys. 12, 972–977 (2016). 48. C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases (Cambridge, 2008), 2nd ed. 49. J. Mendonca and A. Guerreiro, “Time refraction and the quantum properties of vacuum,” Phys. Rev. A 72, 063805 (2005). 50. A. V. Turukhin, V. S. Sudarshanam, M. S. Shahriar, J. Musser, B. S. Ham, and P. R. Hemmer, “Observation of ultraslow and stored light pulses in a solid,” Phys. Rev. Lett. 88, 023602 (2002). 51. A. Figotin and I. Vitebskiy, “Slow light in photonic crystals,” Waves Random Complex Media 16, 293–382 (2006). Vol. 26, No. 15 | 23 Jul 2018 | OPTICS EXPRESS 19295

In many studies (especially, in the context of atomic vapours), pulse propagation in the slow light regime was studied via a somewhat phenomenological transport equation for the pulse envelope [27][28][29], i.e., neglecting second-and higher-order dispersion terms.More accurate models (usually, for pulse propagation in waveguides or photonic crystals) were based on envelope equations like the nonlinear Schrödinger equation, see Section 2 and [53].However, the validity of these models is limited to moderate levels of slow light.Otherwise (especially close to a zero group velocity point (ZGVP)), an exceeding number of dispersion terms (in the form of high-order time derivatives) has to be taken into account, because the various dispersion coefficients are proportional to the inverse of the group velocity (and its powers) and/or because of the slowly convergence of the series of higher-order derivatives of the propagation constant with respect to the frequency [54].As a result, standard pulse propagation schemes (see e.g., [53,54]), as well as coupled mode models (see e.g., [55,56], which are all frequency-based expansions) become inefficient.In these cases, one is forced to resort to numerical solutions of the full Maxwell equations (e.g., using the Finite Differences Time Domain (FDTD) approach) which is far more demanding in terms of computational resources and run times [11-13, 37, 38, 57-59]; such an approach also hampers physical insights.This computational burden could have been circumvented by using the efficient envelope formulations developed in the context of Bragg gratings.Initially assuming a shallow grating contrast, a scheme has been derived for pulse propagation near the ZGVP on the band edge, and most consequent work focused on the dynamics of nonlinear phenomena such as gap solitons [60][61][62][63] and their applications [64].This formulation was then extended for deeper grating contrasts by de Sterke et al. [65] (and then, further extended to higher dimensions [66][67][68] and pump-probe configurations [69,70]) by employing a k • p (i.e., momentum-based) expansion of the (diffraction [67] and) structural dispersion (sometimes referred to as waveguide dispersion) based on non-paraxial directional fluxes and Bloch carrier waves; the resulting expansion consists of a series of spatial (longitudinal) derivatives, instead of the standard time derivatives.In particular, in this derivation, the expansion coefficients are proportional to the group velocity (or its powers), rather than to its inverse.As a result, this formulation is especially efficient in the slow light and stopped light regimes, where the derivative series converges rapidly.
Unfortunately, despite the generality and accuracy of these rigorous models [65,67], they are far from being in wide use.This is related, most likely, to their relative complexity.In the same vein, the (momentum-based) envelope formulations of [65,67,69,70] were limited to configurations involving only a small number of not-so-short interacting pulses, in order to avoid further complexity in the formulation arising from the need to account for adjacent bands.Finally, these formulations were limited to material constituents of weak dispersion.However, material dispersion can have an important role for semiconductor and plasmonic structures, and can become comparable to the structural dispersion.Such systems are also frequently lossy, so that there is also an ambiguity regarding the convenient way to include the absorption in the formulation [71].Thus, overall, there are still no efficient (slowly-varying) envelope-based formulations that enable the study of pulse propagation in strongly dispersive and absorptive media, and specifically, in the slow and stopped light regimes.
In this paper, we aim to bridge this gap in the existing pulse propagation schemes.Specifically, in Section 3 we develop a (momentum-based) envelope formulation for pulse propagation in a system with translation invariance along one direction, i.e., a general waveguide system.We, thus, exploit the translation invariance along z to avoid the need for a Bloch carrier wave, hence, we work with the more convenient (and popular) Fourier carrier wave basis.Our approach shares the advantages of previous (momentum-based) formulations, but is more general than previous works [65,67,69,70] as it is suitable also for non-periodic transverse dielectric structure, e.g., the abundant layered structures, photonic crystal fibers [72] etc.; these typically do not involve additional spectrally-adjacent bands.
In Section 4, we demonstrate the convergence of the derived model for the two generic casesstopped and slow light.In order to emphasize the suitability of our approach to the deep slow light regime and for strongly dispersive materials, we model possibly the simplest geometry that supports a ZGVP and a slow light regime, namely, a plasmonic slot waveguide (metal-dielectricmetal (MDM) layer structure), see Fig. 1(a).The dispersion of the metallic constituent is indeed strong, and the metal constituent involves significant absorption.We demonstrate convergence of the derivative series upon the analytical solution.We also show that unlike the standard approach, where the nonlinearity and loss seem to diverge near the ZGVP, in our formulation, they are regular.This is explained by noting that the relevant criterion for the loss and nonlinear effect is the time rather than distance of propagation, such that the propagation time through a finite length waveguide segment diverges near a ZGVP [7,14,15].More generally, we identify effects known from weakly dispersive systems, and discuss the interplay of dispersion, nonlinearity and absorption.We summarize the main findings and discuss future steps in Section 5.

Standard formulation -frequency-based expansion
The starting point of most derivations of pulse propagation schemes is the wave equation for a scalar electric field Here, E is the electric field, In this work, we focus on a waveguide geometry, so that the position vector, ì r, is split into z, the propagation direction, and ì r ⊥ , the transverse direction; thus, the response function of the media R(ì r ⊥ , t) and its Fourier transform, namely, the permittivity (ì r ⊥ , ω), are functions of the transverse coordinate only.Lastly, R N L (ì r, t) is a nonlinear response function (see e.g.[53] for explicit expressions) and the asterisk, ' * ', specifies a convolution.Eq. ( 1) is valid for slab waveguides (including free space) for TE polarization only.In more general cases, the field has more than one component and the first term of the left hand side of Eq. (1) has to be replaced by ì ∇ × ì ∇ × ì E(ì r, t).At this stage, we simplify the derivation by ignoring absorption; this effect will be fully accounted for when we present our alternative expansion.
In order to eliminate the fast oscillations on the scale of the optical period and wavelength, the electric field is written as a convolution of the transverse field profile and a slowly-varying envelope [56], both written as sums of (strictly real) frequency components, i.e., where A(z, t ) is a dimensionless amplitude, and ê(ì r ⊥ , ω) is the transverse mode profile (customarily, associated with units of electric field), i.e., the solution of the Helmholtz equation, Further, ω is the angular frequency, ω 0 is the center frequency of the pulse, β 0 = β(ω 0 ) = ω 0 n e f f /c is the corresponding (potentially complex) propagation constant and F t ω is a Fourier transform from the frequency to the time domain.(Note that in some works, ê(ì r ⊥ , ω) is replaced by ê(ì r ⊥ , ω 0 ).While for single mode propagation, there is just insignificant differences between the two approaches, the latter approach leads to some spurious mode coupling [15].)Eq. ( 1) is now transformed into the frequency domain, yielding where F ω t represents the inverse Fourier transform from the time domain to the frequency domain, and the convolution is transformed to frequency domain as follows Finally, the spectral envelope, â(z, ω − ω 0 ) is related to the time domain envelope via [53,56], Now, one Taylor expands β 2 (ω) around ω 0 as follows Substituting Eqs. ( 3) and (7) in Eq. ( 4), considering instantaneous Kerr nonlinear response, i.e., R N L = N L |E | 2 , where N L is the nonlinear part of the permittivity, and transforming back to time domain gives where the dispersion coefficients Dm are defined as and where is the characteristic length scale for the manifestation of the nonlinear effect, where is the power of the mode at ω 0 [74].Further simplification is achieved by transforming to a frame of reference co-moving with the pulse such that Eq. ( 8) then becomes where β 0 ≡ d 2 β dω 2 ω 0 (not to be confused with the popular notation for an imaginary value; the latter is denoted below explicitly.), and is the group velocity of the mode at ω 0 .
In order to understand the relative magnitude of the various terms in Eq. ( 13) we transform once again the coordinates as follows where T P is the pulse duration and L disp = T 2 P /β 0 is the dispersion length.The new coordinates are now dimensionless, and the derivatives by these coordinates are O(1).Eq. ( 13) now reads Since the wavelength of the center mode is λ 0 = 2πn e f f β 0 and the pulse spatial width is Z ∼ v g T P , the first (non-paraxial) term and third (mixed derivative) term in Eq. ( 16) scale as λ 0 L di s p and λ 0 Z , respectively.When comparing the magnitudes of all terms, the non-paraxial and the mixed derivative terms are negligible in the standard scenario in which λ 0 Z L disp .The dispersion terms are proportional to Dm T m P , such that the high-order terms in Eq. ( 16) can be neglected if In general, this condition is satisfied for sufficiently long pulses.Then, if dispersion is accounted for only up to second-order, Eq. ( 13) reduces to which is the well-known nonlinear Schrödinger equation [53,54].Although this approach is well-accepted and efficient in most cases, it becomes inefficient in the slow light regime and completely breaks down at a zero group velocity point.The reason for the failure of Eq. ( 18) lies in the neglect of the mixed derivative term in Eq. ( 13), as it is inversely proportional to the group velocity.This term grows in magnitude in the slow light regime, and becomes singular at the ZGVP.Additional singularities arise at the high-order derivatives Dm since the m'th term consist of terms which are proportional to v −m g .Therefore, when approaching the slow light regime or a ZGVP, the derivative series cannot be truncated at the second-order time derivative, and higher-order terms have to be considered, making the standard expansion inefficient.
In order to circumvent these problems, we present below an alternative derivation.

Derivation
The starting point of our derivation is the vectorial wave function where ì E = (E x , E y , E z ) is the electric field vector.Eq. ( 19) is the vectorial version of Eq. ( 1); thus, it is valid for any waveguide cross section and any polarization.For lossy media, all coefficients in the following derivation attain complex values.
As in former approaches, the convolution is transformed to the frequency domain.However, unlike the standard approach, we do not proceed by separation of variables, replacement of the permittivity with the propagation constant, and Taylor expanding the latter around the central frequency.Instead, inspired by the momentum-based formulations [65,67], we aim to decompose the fields into a sum of spatial frequencies, i.e., in momentum space.First, however, we expand the permittivity around the central frequency, ω 0 , and then transform back to the time domain, giving Now, we can rewrite the electric field and related physical quantities in terms of their spatial Fourier transforms, namely, where as in the standard derivation (Section 2), ã(β − β 0 , t) is a dimensionless amplitude and ì e(ì r ⊥ , ω) is the transverse mode profile, associated with units of electric field.Here, and in what follows, we denote the spatial Fourier transforms by a tilde symbol (˜), in order to distinguish them from the (more standard) temporal Fourier transforms used in the standard derivation (Section 2), which are denoted by a hat symbol (ˆ).
We emphasize that as we consider a spatial (rather than a temporal) Fourier transform, we inherently work with real β values.Thus, in the presence of absorption (i.e., when the permittivity of the constituents is complex), the frequency ω obtains complex values.In this case, the dispersion curve remains similar to the curve in the loss-free case (see Fig. 1) and therefore, the ZGVP and slow light regime are preserved [37,38,71].This means that the evaluation of the modal frequencies and various other coefficients is done with complex arguments.This is an unavoidable complication of the use of complex frequencies [75].
Substituting Eqs. ( 20) and ( 24) in Eq. ( 19), and transforming into momentum space gives Substituting Eq. ( 25) in Eq. ( 26) allows us to replace the first two terms by (ì r ⊥ , ω)ω 2 a(β − β 0 , t) ì e(ì r ⊥ , β).We now get In order to remove the dependence on the transverse coordinate, we now multiply by ì e * (ì r ⊥ , β) and integrate over ì r ⊥ .Then, Eq. ( 27) becomes where Note that the terms in parentheses are evaluated at ω 0 ≡ ω(β 0 ) (rather than the other way around in the standard approach), but ì e(ì r ⊥ , β) is evaluated at a general β.Thus, one can see that (in the absence of absorption, Otherwise, U [76] becomes a non-trivial complex function.)U(β 0 , β 0 ) is related to the energy density of the carrier mode β 0 (see Appendix A), but U(β, β 0 ) equals this quantity only to leading-order.As noted many times in the past, the energy represents the most convenient (and natural) norm for solutions of the vector Helmholtz equation for the electric field (see e.g., discussion in [77]).
Since ω = ω(β), we can Taylor expand the ω-dependent coefficients in terms of β, in analogy to the expansion of β 2 in terms of ω in the standard approach (7).Thus, Eq. ( 28) becomes where and Note that both Dm and Dm,n account for the dispersion of the material parameters, the mode and the "energy density" U; The index n represents the order of material dispersion, and the index m represents the order of the combined effect of material and structural dispersion.Note that our formulation involves a differentiation of the mode profiles themselves ( ì e(ì r ⊥ , β)), a complication which is avoided in the standard formulation.Such complication arises because the initial temporal Fourier transform (Eq.( 20)) does not allow one to get rid of the mode profiles, as in the standard derivation.
Transforming Eq. ( 30) back to real space gives i D2,0 A(z, t) where and In Eqs. ( 34)-( 35), we used Eqs.( 31)- (32), where it follows from the chain rule that d (ì r ⊥ ,ω) The coefficients ω 0 and ω 0 stand for dω dβ β 0 and d 2 ω dβ 2 β 0 , respectively.Eq. ( 33) is our main result; it is an exact reduction of Maxwell's equations, analogous to the standard derivation, Eq. ( 8), but with the roles of z and t exchanged (as in previous momentumbased expansions [65,67]).Thus, for example, it accounts for non-paraxiality via a second-order derivative in time rather than in space, i.e., in z, and for dispersion via a series of derivatives in space rather than time.(In a 3D configuration, the spatial derivatives also account for diffraction effects [67].)Another manifestation of the exchange of roles of time and space is seen through D2 , which represents the group velocity dispersion (GVD).It consists of two terms.The first is the second derivative of the dispersion curve, whereas the second term is the group velocity at the central wavenumber β 0 in the current approach, similarly to Eq. ( 14), defined as For materials with negligible absorption, the second term in Eq. ( 34) is proportional to v 2 g , in analogy to the GVD coefficient in the standard derivation (where the second term is proportional to v −2 g , see Eq. ( 9) and [53,54]).As in the standard (frequency-based) derivation, the relative magnitude of the two terms in Eq. ( 34) depends on the position of the central mode on the dispersion curve.However, near a ZGVP/in the slow light regime, the second contribution to the GVD vanishes, in contrast to the standard derivation, where the second term diverges, again, due to the exchanged roles of time and space.
The first two orders of dispersion are analogous to those in the standard derivation (8), however, the higher-order dispersion terms in (33) involve additional mixed derivative terms.In that respect, our derivation yields a more complicated model compared with the standard one.
Eq. ( 33) generalizes the final results of the previous momentum-based expansions [65,67] to lossy and strongly dispersive media, including in particular, the slow light regime and near ZGVPs.Importantly, in these regimes, Eq. ( 33) has a major advantage.Indeed, the series of dispersion terms which diverges in the standard derivation (18) as v −m g is replaced by a series of terms that scale with v m g , hence, is rapidly converging.In addition, the high-order derivatives of ω, are also small in the slow light regime.

Dimensional analysis
We now assume that the nonlinear response is spatially-local, i.e., that R N L represents a cubic (Kerr-like) response, R N L = N L | ì E | 2 ; Furthermore, U (29) is evaluated at β 0 in the nonlinear term.In order to further simplify the analysis, we also would like to adopt the slowly-varying envelope approximation.However, this has to be done with care.Indeed, the first term in Eq. ( 33) (the second derivative in time) is not necessarily small, even if the envelope is slowly varying -it consists of slow changes in time (∼ T disp , see below), but also a (v g /Z) 2 scale change due to the propagation; this term is comparable to D2,0 ω 0 2 (which is dominant, as explained in Section 3.1, see Eq. ( 34)), as in the standard derivation [54].Yet, the latter contribution can be removed after the coordinates are transformed to a frame co-moving with the pulse, namely, Accordingly, the derivatives of the coordinates transform as Then, Eq. ( 33) becomes i D2,0 where is the characteristic time for the nonlinear effect to become significant.For lossy media, ω 0 is complex, so that the ∂/∂t term does not vanish, but rather, leaves behind its imaginary part, ω 0 ; note that it does not represent a velocity, but rather the dispersion of the absorption.Similarly, the D2,0 ω 0 2 term does not cancel completely either.In order to reveal the hierarchy of the terms in Eq. ( 39), we normalize the coordinates as follows where Z is the initial pulse width and T disp = Z 2 /ω 0 is the dispersion time, the analogue of the dispersion length L disp .Then, Eq. (39) obtains the dimensionless form where the dimensionless coefficients are The second and fifth terms in Eq. ( 42) are of the order of unity, and typically the dominant.Now, in order to estimate the dispersion terms, we assume that the waveguide includes highly dispersive constituents (e.g., metal layers as in Fig. 1(a)).Thus, the estimates that follow provide an upper bound for the magnitude of these terms, while for dielectric materials, the dispersion terms would be typically smaller.Thus, D2,0 can be estimated as 1/ω 0 (see Appendix B).Then, f ∼ T 0 /T disp and g ∼ T 0 /T P where we defined the period T 0 = 2π/ω 0 , and the pulse duration, T P , which scales as Z/v g .Thus, in the typical case where T 0 T P T disp , we can omit the first (non-paraxiality) and the third (mixed derivative) terms in Eq. ( 42); this is analogous to the standard derivation (Section 2).
The fourth and the sixth terms represent the effect of the complex ω terms on the propagation in the presence of absorption, namely, they account for the first-and second-order dispersion of the absorption.Generically, the absorption dispersion is a weak effect even for strongly absorbing media, like metals; accordingly, it is safe to neglect these terms.
Taking into account the scaling of Dm and Dn,m (see Appendix B), one can see that the seventh and eighth terms scale as (T disp /T 0 ) (β 0 Z) −m and (T disp /T 0 ) (β 0 Z) −m (v g T disp /Z) k (T disp /T 0 ) −n , respectively.Typically, T disp /T 0 is much greater than unity; β 0 Z and (v g T disp /Z) are also larger than unity, but not as much.Thus, since m ≥ 3 in the seventh term, it is typically smaller than unity; further, since n ≥ k and T disp /T 0 is larger than v g T disp /Z, it is possible to show that the eighth term is typically smaller than the seventh term in the slow and stopped light regimes.In the specific case of propagation at the ZGVP (slow light regime), all terms proportional to v g vanish (are small).In all other cases, one has to carefully map the relative magnitudes of the many higher-order terms.
The scaling of the ninth (nonlinear) term turns naturally to be by U(β 0 , β 0 ) (rather than as P 0 , as in the standard derivation, see Eq. ( 40)).In lossy media, T N L becomes complex as it contains the complex values U(β 0 , β 0 ) and ω 0 .Overall, the nonlinear term becomes important in cases where T N L is comparable or shorter than T disp .
If all mixed terms are indeed small, Eq. ( 39) becomes where the * symbol was removed for brevity such that z and t are the (dimensional) coordinates of the moving frame (37).In that sense, the additional complexity of the model does not exceed that of the standard model for pulse propagation (Section 2 and [53]).

Scaling with the group velocity
Many previous studies showed that linear effects such as absorption scale with v −1 g and that cubic nonlinear effects scale with v −2 g , see e.g., [7,14,15,67].The former effect appears only implicitly in the scaling of the linear loss term with the group velocity, e.g., via the dispersion relations [71]; the latter is seen from the natural scaling of the cubic nonlinearity with the power P 0 (see Eq. ( 11)), or explicitly in the momentum-based expansions of de Sterke, Sipe et al. (e.g., [65,67]) where the dynamics is described by directional fluxes.
This inverse proportionality with the group velocity implies that these effects are enhanced in the slow light regime, and diverge near a ZGVP!This behaviour motivated many studies of nonlinear effects in the slow light regime [7,[14][15][16][17][18][19][20][21]67].Yet, rather remarkably, our final equation does not show any particular enhancement of (linear and) nonlinear effects in the slow and stopped light regimes.The solution to this seeming paradox is that the strength of the various effects is related to the time of propagation, rather than to the distance of propagation.Thus, if the dynamics is monitored in time, there is no particular sensitivity of the strength of nonlinear effects on the group velocity.In other words, the nonlinear effects are enhanced for low group velocities since a longer time of propagation is required to traverse a given distance in space, but per unit time, the nonlinear effect has the same strength.This is seen implicitly in Fig. 1(c) that shows that T abs is only weakly dependent on β, such that the attenuation scales linearly with t (and is independent of z).Furthermore, the nonlinear term in our formulation scales with the energy density rather than with the power, hence, is independent of v g .
This result also shows that the enhancement of the nonlinear effects reported before were a result of the experimental configuration -the pulse is injected into a finite length nonlinear media, and the nonlinear effect is measured at the output.Thus, unlike many previous claims, light-matter interactions are not enhanced in the slow/stopped light regimes, and the observed enhancements of the (linear and) nonlinear effects is unrelated to the source of the strong dispersion, be it material or structural dispersion [78,79].An alternative experimental configuration where the absorption or nonlinearity are measured in systems with different group velocities for the same duration is expected to show the insensitivity to the group velocity we predict here.
The disappearance of the dependence on v g is a result of two rather subtle mathematical differences our derivation has in comparison to previous studies.First, our formulation is based on momentum modes (rather than frequency modes).Second, unlike [67], we adopt a frame moving in space (see Eq. ( 37)) rather than in time (12).

Examples
We now demonstrate the validity of our derivation by solving Eq. (44).By its very nature, Eq. ( 44) is solved as an initial value problem (IVP) in time.Specifically, the initial condition is set as where β 0 is the center of the spectrum, ∆β is the spectral width of the pulse and C is a normalization constant such that A(z = 0, t = 0) = 1.This is similar to how the Schrödinger equation appears in quantum mechanics, but is in contrast to the standard approach in optics, where the equation of propagation is solved as an IVP in space, see e.g., [53] and Section 2. Due to the above, it is also not straightforward to compare the two methods directly.
In order to demonstrate the dynamics of short pulses in the slow light regime and near a ZGVP, we consider the simplest geometry that supports such regimes, namely, an MDM waveguide (see Fig. 1(a)).The permittivity of the metal is assumed to be given by the Drude model, where ∞ is a constant which determines the value of the permittivity when ω → ∞, ω P is the plasma frequency and γ is a constant real number which relates to the absorption in the metal.Determination of the coefficients in Eq. ( 44) requires calculating the dispersion relations, ω(β) and the mode profiles.If absorption is neglected (γ = 0), one finds a unique set of solutions which satisfy the dispersion relation, where both β and ω attain real values.However, in the presence of absorption, there are two generic sets of solutions to the dispersion relations [80].
The first possibility is to account for the absorption via complex-β values, as in the standard derivation (Section 2).In this case, the dispersion curve splits at the ZGVP, resulting in two branches of fast light [71]; in addition, the dispersion relation curve backbends near plasmon resonance [80].We note that the standard approach does not totally fail in this case since the terms that include factors such as 1/v m g are large but do not diverge anymore.However, the slow light regime becomes limited, and it is not a-priory clear how to link the resulting complex modal structure and the original ZGVP that disappeared.
The second possibility is to account for the absorption via complex ω values.In this case, the dispersion curve remains similar to the curve in the loss-free case (see Fig. 1(a)-(b)) and therefore, the ZGVP and slow light regime are preserved [37,38].As noted above, our approach is inherently based on complex frequency modes, which thus allows for a convenient description of the pulse dynamics near a ZGVP [37,38].Note that while the mode profiles exhibit outgoing energy flow away from the interface, they are also exponentially decaying due to the negative real part of the metal permittivities in the cladding.In that respect, and unlike the structures studied in [35][36][37][38], the modes we study are below the light line, hence, they do not suffer from radiation losses.As a result, normalization of the modes is trivial, and free from the complications arising frequently in the normalization of exponentially growing complex frequency modes [75].
Once ω(β) was computed, we ensured that the non-paraxiality and various mixed derivative terms are indeed relatively small, so that Eq. ( 44) provides an excellent approximation for the exact model, Eq. ( 42).
When simulating pulse propagation according to Eq. ( 44) and considering absorption, we have to recall that e −iω 0 t in Eq. ( 22) also describes the attenuation of the pulse in time as ω 0 consists of an imaginary part.Therefore, in the following plots, we will present the value | A(z, t)e −iω 0 t | instead of just the amplitude, A(z, t).In that sense, the characteristic time for absorption is T abs = 1/ (ω 0 ).

Simulations of pulse propagation in the slow and stopped light regimes; linear media
We first study the linear dynamics by solving Eq. ( 44) without the nonlinear term.In order to test the convergence of the derivative series in Eq. ( 44), we set Ω (M) (β) ≡ M m=1 Dm (β − β 0 ) m , such that the solution of Eq. ( 44) is Here, M accounts for the order beyond which dispersion terms in Eq. ( 44) are truncated.For a reference, we use the solution of the non-truncated series, found by repeating the procedure with Ω (∞) (β) ≡ ω(β).Since there is no source for backward scattering in our case, it can be shown that this approach is equivalent to solving exactly the vector wave equation ( 19), i.e., taking into account even the non-paraxiality and various mixed-derivative terms; see e.g., [81].Hence, we can refer to the Ω (∞) solution as an exact reference.
Our first example presents pulse propagation in the stopped light regime, where β 0 was chosen to be at the ZGVP, see Fig. 1(b).The spectrum width is chosen to be ∆β ≈ 0.1β 0 , to ensure convergence of the dispersion expansion.In Fig. 2(a), we demonstrate the rapid convergence of Eq. ( 44) upon the exact solution in the absence of absorption.Indeed, the simulation results for the field at t = T disp decently match the exact solution when M = 2, and becomes more accurate when M = 3.Specifically, the relative difference of the outcome was tested according to where A (M ma x ) is the reference (the most accurate simulation).In the linear simulation we set M max = ∞.Here, we measured R 2,∞ = 1.1 and R 3,∞ = 0.035.This demonstrate the rapid convergence of Eq. ( 44) with low-order dispersion.
In the following simulations, we chose a narrower spectral regime such that ∆β = 0.01β 0 .The solution of Eq. ( 44) in the absence of absorption is shown in Fig. 3(a).The pulse is centered at z = 0 at all times, as expected, but broadens due to dispersion.In the presence of realistic absorption (see Fig. 3(b)), T abs turns out to be much shorter than T disp , such that the field decays rapidly within t < T disp .Figs. 4(a)-(b) show the respective dynamics of the spatial spectrum of the pulse in these two cases.The spectrum does not change over time if we only consider dispersion, however, its intensity is attenuated in the presence of absorption, as expected.
In order to show the validity of Eq. ( 44) in the slow light regime, we repeated the process above for β 0 centered at the slow light regime where v g = 5.5 × 10 −4 c (see Fig. 1(b)).Here, as earlier, we choose ∆β ≈ 0.1β 0 .However, although ∆β here is nominally wider than the spectral width chosen in the case of the ZGVP (since the value of β 0 is higher), the GVD and higher-order dispersion coefficients are much smaller.As a result, satisfactory accuracy is obtained even for M = 2, as shown in Fig. 2(b), where R 2,∞ = 0.1 and R 3,∞ = 0.005.This is, again, in contrast to the standard approach (Section 2), where due to the large factors ∼ v −m g in the mth-order terms, many more terms are required for convergence.Fig. 3(c) shows the pulse dynamics in the reference of the moving pulse.The results are remarkably similar to the dynamics in the ZGVP (Fig. 3(a)).In that respect, it is worth emphasizing that one aspect of this similarity is somewhat surprising.Indeed, Fig. 1(c) shows that T abs in the slow light regime is comparable to that at the ZVGP, although the latter is not as close to the SPP resonance.This can be understood by noting that the well-known enhanced absorption near SPP resonance is a result of the longer propagation time through a finite length waveguide segment, i.e., it is a direct result of the fact that the propagation is typically formulated as an IVP in space, an approach that makes the long propagation distance implicit (see also discussion in Section 3.3).Nevertheless, T abs indeed becomes shorter when the pulse spectrum is near resonance.This is a result of the increasing portion of energy residing in the metal side.We will see below a similar behaviour in the context of the nonlinear regime.

Simulations of pulse propagation in the slow and stopped light regimes; nonlinear media
We now repeat the simulations shown in Section 4.1 in the presence of a nonlinear term, now employing the split step method [53].Under these circumstances, one has to consider the interplay between dispersion, absorption and nonlinear dynamics via their characteristic time scales T disp ,   T abs and T N L ; in the presence of absorption, we replace Overall, the nonlinearity results in self-focusing in space and spectral broadening (due to self-phase modulation).As a result, the higher-order dispersion terms become relatively more important than in the linear cases studied above.
We first study pulse propagation at the stopped light regime.We assume a negative nonlinearity which counters the effects of the normal dispersion occurring in this regime.The outcome of this simulation in the absence of absorption, where we set T N L T disp (T abs is infinite) is presented in Fig. 5(a), and the corresponding spectrum, is presented in Fig. 6(a).In this case, we set M = 2, after we found that R 2,3 = 0.008.The corresponding dynamics in the presence of absorption is shown in Fig. 5(b), and Fig. 6(b).A comparison of Fig. 5(a) to Fig. 3(a) shows that since the nonlinear effect is much stronger than the dispersion (for this choice of parameters), the pulse exhibits rapid self-focusing.In momentum space, when comparing Fig. 6(a) to Fig. 4(a), we observe the corresponding behavior where the spectral width broadens, and refocuses again [53].
When considering absorption, in order to observe a significant nonlinear effect, we choose T N L to be much shorter than T abs , thus, also much shorter than T N L in the absence of absorption.If T N L was longer, the absorption would attenuate the pulse before any nonlinear effect can take place.When comparing Fig. 5(b) and Fig. 3(b), it seems that in terms of duration, the nonlinear process does not change the fate of the pulse, as the pulse is attenuated very quickly in both cases.Nevertheless, in momentum space, when comparing Fig.  broadened very quickly due to self-phase modulation.
Next, we study the influence of nonlinearity on pulse propagation in the slow light regime.We choose the same initial pulse as in the slow light regime simulation in Section 4.1.In this regime, the dispersion is anomalous, so we set a positive nonlinear term to counteract the effects of the dispersion.As in the case of the ZGVP, we set T N L such that nonlinear effects are stronger than the dispersion and omit the absorption.Such configuration results in T N L T disp (T abs is infinite).The associated pulse dynamics is shown in Fig. 5(c) and the corresponding spatial spectrum is presented in Fig. 6(c).We note that M = 2 is sufficient in this case as R 2,3 = 0.001.Overall, as in the linear case, the dynamics in the slow light regime is qualitatively (and even quantitatively) similar to that in the ZGVP.Specifically, T N L in both cases is similar, and also similar to cases where the mode has a "normal" velocity, see discussion in Section 3.3.

Discussion
In this paper, we developed an efficient formulation suitable for describing pulse propagation in the slow and stopped light regimes.It enables simple simulations for far more extreme regimes of dispersion than performed before.However, since the linear pulse propagation problem has an analytical solution, the main contribution of this part of our work is in providing a somewhat new point of view on pulse propagation in lossy and nonlinear media, which may in fact be intuitive.Indeed, our formulation shows that the relevant quantity is the propagation time (rather than propagation space) -if measured in time, the rate of phase accumulation is not different than in weakly dispersive systems.Thus, the apparent enhanced phase accumulation (hence, nonlinear effect) and absorption in the slow light regime are a direct result of the longer time spent in the system.This is seen also in the nonlinear simulations (Section 4.2) which do not have an analytic solution.
The rigorous derivation of our main equations is quite subtle, as the lengthy derivation shows, especially, if one has to go beyond second-order dispersion.However, a-posteriory, the final form is quite intuitive as its first-and second-order terms could have even be guessed due to the close analogy to the standard formulation.This familiar form allows us enjoy the wide knowledge accumulated from studies of its standard analogue, including, in particular, the availability of analytical solutions for Gaussian pulses, effects associated with the pulse shape, the understanding of the interplay between dispersion and Kerr nonlinearity (see [53]) and the nonlinear dynamics within the (generalized) nonlinear model, in general.The formulation also shows how to handle material absorption in a convenient manner, via complex frequencies.The advantage of their use was not well understood at the time when the more complex momentum-based formulations were derived.For example, the tedious FDTD calculations of the threshold-less lasing configuration [35][36][37][38] and wave mixing and free carrier generation effects [15,82,83] could now be replaced with simple, even analytic, solutions.
Our formulation can be generalized for a weak longitudinal grating, by adopting the approach in [60][61][62][63].The formulation will also be trivially valid for pulse propagation in (uniform) atomic media [2,3,30,31] and quantum models of stopped light [26][27][28], which as noted, were so far usually modelled with rather simplistic models in which second-order and higher-order dispersion effects were neglected.
In addition, the importance of the current study is in performing the first step towards the extension to more peculiar/demanding problems of (linear and especially nonlinear) pulse propagation that do not have an analytical solutions like, e.g., photonic crystals with defects or degenerate band edges [51], or plasmonic photonic crystals, which are absorptive and highly dispersive [84,85].Further, our formulation is a starting point for the study of more complex configurations involving the coupled dynamics of various modes.In this case, one would be able to derive coupled mode analysis in the slow and stopped light regime, so far treated only with the full Maxwell equations via FDTD algorithms [6, 11-13, 37, 38, 57-59].Finally, our results could also be used in many additional contexts such as acoustic waves [43][44][45], spin waves [46], water waves [47], quantum/matter waves [48,49] We now show that although the constituent materials are dispersive, the relation that stands for the equality of the electric and magnetic energies in dispersionless media, namely, where we used ∫ h r i ẽr j dr k r i =∞ where i, j, k represent any spatial coordinate.We now substitute Faraday-Henri Law in Eq. ( 51) and use Eq.
In the final stage, we substitute Since | ì e(ì r ⊥ , β 0 )| accounts for the total amplitude of the electric field (i.e., max(A) = 1), Eq. ( 57) is the expression for energy of the carrier mode in dispersive, loss-free, non-magnetic media.

Fig. 1 .
Fig. 1.(a) Cross section of a one dimensional plasmonic MDM waveguide of width d.The solid line features a transverse profile of the magnetic field propagating in z direction, corresponding to the antisymmetric branch.(b) Dispersion curve of antisymmetric plasmonic modes in a MDM slab waveguide where the core width is 35 nm, the permittivity of the dielectric is D = 2.5 and the parameters of the permittivity of the silver (see Eq. (45)) are ∞ = 5, ω P = 1.4 × 10 16 s −1 and γ = 3.2 × 10 13 s −1[73].In cases where we ignore absorption, we set γ = 0.The red segment on the left is a stopped light regime where forward and backward modes coexist.Here, β 0 = 1.51 × 10 8 m −1 (λ 0 ≈ 363.6nm) and ∆β = 0.15 × 10 8 m −1 , where β 0 lies exactly at the ZGVP.The second red segment lies in a slow light regime where β 0 = 3 × 10 8 m −1 (λ 0 ≈ 362.9nm) and ∆β = 0.3 × 10 8 m −1 .Here β 0 lies at a point where the group velocity is 5.51 × 10 −4 c.The green line on the left represents the light line, corresponding to c/ √ D .In the presence of absorption, the dispersion curve (and the ZGVP) is shifted slightly (not shown).(c) The values of T abs = 1/ (ω 0 ) corresponding to the mode shown in (b).

Fig. 2 .
Fig. 2. Pulse propagation at (a) ZGVP and (b) slow light regime, where the spectral width of the pulse, ∆β and the center value of the propagation constant, β 0 are the corresponding values shown in Fig. 1(b).The green dashed line is A(z, t = T disp ) where M = 2, the blue dotted line corresponds to M = 3 and the red solid line represents the exact solution M = ∞.

Fig. 5 .
Fig. 5. Spatio-temporal dynamics of the amplitude of the electric field (color represents value | A(z, t)|e − (ω 0 )t ) for a pulse in the non-linear case.Plots (a) and (b) represent pulse propagation in the ZGVP regime in the absence and presence of absorption, respectively.Here, the time scales of the dispersion time and the absorption time are exactly as in Figs 3(a)-(b).The nonlinearity time scale is T N L = 0.2T disp in the absence of absorption (a) and T N L = 0.1T abs when absorption is present (b).Plots (c) and (d) represent pulse propagation in the slow light regime in the absence and presence of absorption, respectively.Here, the dispersion and absorption time scales are exactly as in Figs 3(c)-(d).The nonlinearity time scale is T N L = 0.2T disp in the absence of absorption (c) and T N L = 0.1T abs when absorption is present (d).In both cases, the total time of propagation is T N L and 10T abs , respectively.

Fig. 5 (
b) and (d) (as well as in Fig. 6(b) and (d)) show the dynamics in the presence of absorption.