Pulse propagation in the slow and stopped light regime

Tal A. Weiss; Yonatan Sivan

doi:10.1364/OE.26.019294

1. Introduction

Electromagnetic waves with slow propagation velocities (slow light) were observed and used in the microwave range as early as the 1940s [1]. In recent decades, slow light was demonstrated also at optical frequencies [2–7], and it was even shown that the light can nearly be brought to a standstill (“stopped light”) [8–10]. As a result, various potential applications emerged, including in particular optical memory in optical communication systems [6, 11–13], nonlinear optics [7, 14–21], Brillouin scattering [22,23], and more recently, quantum computing [10,24–29], e.g., based on electromagnetic induced transparency [2, 3, 8, 9, 30–32]; Additional potential applications are tissue imaging [33], sensing [34], low threshold lasing [35–39], photochemistry [40,41] and even tabletop cosmology [42]. Slow and stopped light regimes occur also in other wave systems like acoustic waves [43–45], spin waves [46], water waves [47], quantum/matter waves [48,49] etc..

Slow light has been demonstrated in various platforms, including atomic vapours [2,3,30,31], solids [50], and various photonic structures, such as Bragg gratings, photonic crystal fibers and waveguides etc., see [4–7,51,52].

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–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–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–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 cases - stopped 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.

2. Standard formulation - frequency-based expansion

The starting point of most derivations of pulse propagation schemes is the wave equation for a scalar electric field

\nabla^{2} ℰ (\vec{r}, t) - μ_{0} \frac{\partial^{2}}{\partial t^{2}} (R ({\vec{r}}_{⊥}, t) * ℰ (\vec{r}, t)) - \frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}} (R_{NL} (\vec{r}, t) ℰ (\vec{r}, t)) = 0 .

Here, ℰ is the electric field,

\nabla^{2} = \partial_{z}^{2} + \partial_{⊥}^{2}

, μ₀ is the vacuum permeability, t is time, r⃗ is the position vector.

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_NL(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 ∇⃗ × ∇⃗ × ℰ⃗(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.,

ℰ (\vec{r}, t) = e^{i β_{0} z} ℱ_{ω}^{t} [(\int d t^{'} A (z, t^{'}) e^{i (ω - ω_{0}) t^{'}}) (\int d t^{'} ℱ_{ω}^{t^{″}} [\hat{e} ({\vec{r}}_{⊥}, ω)] e^{i ω t^{″}})],

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,

\nabla_{⊥}^{2} \hat{e} ({\vec{r}}_{⊥}, ω) - β^{2} (ω) \hat{e} ({\vec{r}}_{⊥}, ω) + μ_{0} ∊ ({\vec{r}}_{⊥}, ω) ω^{2} \hat{e} ({\vec{r}}_{⊥}, ω) = 0 .

Further, ω is the angular frequency, ω₀ is the center frequency of the pulse, β₀ = β(ω₀) = ω₀n_eff/c is the corresponding (potentially complex) propagation constant and

ℱ_{ω}^{t}

is a Fourier transform from the frequency to the time domain. (Note that in some works, ê(r⃗_⊥, ω) is replaced by ê(r⃗_⊥, ω₀). 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

\begin{array}{l} \hat{a} (z, ω - ω_{0}) \nabla_{⊥}^{2} \hat{e} ({\vec{r}}_{⊥}, ω) + 2 i β_{0} μ_{0} ∊ ({\vec{r}}_{⊥}, ω) \hat{e} ({\vec{r}}_{⊥}, ω) \frac{\partial \hat{a} (z, ω - ω_{0})}{\partial z} \\ - β_{0}^{2} μ_{0} ∊ ({\vec{r}}_{⊥}, ω) \hat{e} ({\vec{r}}_{⊥}, ω) \hat{a} (z, ω - ω_{0}) + ω^{2} μ_{0} ∊ ({\vec{r}}_{⊥}, ω) \hat{e} ({\vec{r}}_{⊥}, ω) \hat{a} (z, ω - ω_{0}) \\ - \frac{1}{c} {(ω - ω_{0})}^{2} ℱ_{t}^{ω} (R_{NL} (\vec{r}, t) ℰ (\vec{r}, t)) = 0, \end{array}

where

ℱ_{t}^{ω}

represents the inverse Fourier transform from the time domain to the frequency domain, and the convolution is transformed to frequency domain as follows

R ({\vec{r}}_{⊥}, t) * ℰ (\vec{r}, t) = ℱ_{ω}^{t} [∊ ({\vec{r}}_{⊥}, ω) \hat{e} ({\vec{r}}_{⊥}, ω) \hat{a} (z, ω - ω_{0}) e^{i β_{0} z}] .

Finally, the spectral envelope, â(z, ω − ω₀) is related to the time domain envelope via [53,56],

ℱ_{t}^{ω} [A (z, t) e^{- i ω_{0} t}] = \hat{a} (z, ω - ω_{0}) .

Now, one Taylor expands β²(ω) around ω₀ as follows

β^{2} (ω) = {\sum_{m = 0}^{\infty} \frac{1}{m!} \frac{d^{m} (β^{2})}{d ω^{m}} |}_{ω_{0}} {(ω - ω_{0})}^{m} .

Substituting Eqs. (3) and (7) in Eq. (4), considering instantaneous Kerr nonlinear response, i.e., R_NL = ∊_NL |ℰ|², where ∊_NL is the nonlinear part of the permittivity, and transforming back to time domain gives

\frac{i}{2 β_{0}} \frac{\partial^{2} A (z, t)}{\partial z^{2}} + \frac{\partial A (z, t)}{\partial z} + \sum_{m = 1}^{\infty} {\hat{D}}_{m} \frac{\partial^{n} A (z, t)}{\partial t^{m}} + \frac{i}{L_{NL}} {| A (z, t) |}^{2} A (z, t) = 0,

where the dispersion coefficients D̂_m are defined as

{\hat{D}}_{m} = \frac{1}{2 β_{0}} \frac{{(- i)}^{m - 1}}{m!} {\frac{d^{m} (β^{2})}{d ω^{m}} |}_{ω_{0}},

and where

L_{NL} = \frac{2 P_{0}}{∊_{NL} ω_{0} \int {| \hat{e} ({\vec{r}}_{⊥}, ω_{0}) |}^{4} d {\vec{r}}_{⊥},}

is the characteristic length scale for the manifestation of the nonlinear effect, where

P_{0} = \frac{β_{0}}{μ_{0} ω_{0}} \int {| \hat{e} ({\vec{r}}_{⊥}, ω_{0}) |}^{2} d {\vec{r}}_{⊥},

is the power of the mode at ω₀ [74].

Further simplification is achieved by transforming to a frame of reference co-moving with the pulse such that

z^{*} = z, t^{*} = t - z / v_{g} .

Eq. (8) then becomes

\begin{array}{l} \frac{i}{2 β_{0}} \frac{\partial^{2} A (z^{*}, t^{*})}{\partial z^{* 2}} & + \frac{\partial A (z^{*}, t^{*})}{\partial z^{*}} - \frac{i}{2 β_{0} v_{g}} \frac{\partial^{2} A (Z^{*}, t^{*})}{\partial z^{*} \partial t^{*}} - \frac{i}{2} β_{0}^{″} \frac{\partial^{2} A (z^{*}, t^{*})}{\partial t^{* 2}} \\ + \sum_{m = 0}^{\infty} {\hat{D}}_{m} \frac{\partial^{m} A (z^{*}, t^{*})}{\partial t^{* m}} A (z^{*}, t^{*}) + \frac{i}{L_{NL}} {| A (z^{*}, t^{*}) |}^{2} A (z^{*}, t^{*}) = 0, \end{array}

where

{β_{0}^{″} \equiv \frac{d^{2} β}{d ω^{2}} |}_{ω_{0}}

(not to be confused with the popular notation for an imaginary value; the latter is denoted below explicitly.), and

v_{g} \equiv {{(\frac{d β}{d ω})}^{- 1} |}_{ω_{0}} = β_{0}^{' - 1},

is the group velocity of the mode at ω₀.

In order to understand the relative magnitude of the various terms in Eq. (13) we transform once again the coordinates as follows

\hat{z} = \frac{z}{L_{disp}},

\hat{t} = \frac{t^{*}}{T_{P}},

where T_P is the pulse duration and

L_{disp} = T_{P}^{2} / β_{0}^{″}

is the dispersion length. The new coordinates are now dimensionless, and the derivatives by these coordinates are O(1). Eq. (13) now reads

\begin{array}{l} \frac{i}{2 β_{0} L_{disp}} \frac{\partial^{2} A (\hat{z}, \hat{t})}{\partial {\hat{z}}^{2}} + \frac{\partial A (\hat{z}, \hat{t})}{\partial \hat{z}} - \frac{i}{2 β_{0} v_{g} T_{P}} \frac{\partial^{2} A (\hat{z}, \hat{t})}{\partial \hat{z} \partial \hat{t}} - i \frac{\partial^{2} A (\hat{z}, \hat{t})}{\partial {\hat{t}}^{2}} + L_{disp} \sum_{m = 3}^{\infty} {\hat{D}}_{m} \frac{\partial^{m} A (\hat{z}, \hat{t})}{\partial t^{* m}} \\ + i \frac{L_{disp}}{L_{NL}} {| A (\hat{z}, \hat{t}) |}^{2} A (\hat{z}, \hat{t}) = 0 . \end{array}

Since the wavelength of the center mode is

λ_{0} = \frac{2 π n_{eff}}{β_{0}}

and the pulse spatial width is Z ∼ v_gT_P, the first (non-paraxial) term and third (mixed derivative) term in Eq. (16) scale as

\frac{λ_{0}}{L_{disp}}

and

\frac{λ_{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 λ₀ ≪ Z ≪ L_disp. The dispersion terms are proportional to

\frac{{\hat{D}}_{m}}{T_{P}^{m}}

, such that the high-order terms in Eq. (16) can be neglected if

{\hat{D}}_{m} ≪ T_{P} {\hat{D}}_{m - 1} .

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

\frac{\partial A (z^{*}, t^{*})}{\partial z^{*}} - \frac{i}{2} β_{0}^{″} \frac{\partial^{2} A (z^{*}, t^{*})}{\partial t^{* 2}} + \frac{i}{L_{NL}} {| A (z^{*}, t^{*}) |}^{2} A (Z^{*}, t^{*}) = 0,

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 D̂_m since the m’th term consist of terms which are proportional to $v_{g}^{- m}$ . 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.

3. Alternative formulation - momentum-based expansion

3.1. Derivation

The starting point of our derivation is the vectorial wave function

\vec{\nabla} \times \vec{\nabla} \times ℰ (\vec{r}, t) - μ_{0} \frac{\partial^{2}}{\partial t^{2}} (R ({\vec{r}}_{⊥}, t) * \vec{ℰ} (\vec{r}, t)) - \frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}} (R_{NL} (\vec{r}, t) \vec{ℰ} (\vec{r}, t)) = 0 .

where ℰ⃗ = (ℰ_x, ℰ_y, ℰ_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, ω₀, and then transform back to the time domain, giving

R ({\vec{r}}_{⊥}, t) * \vec{ℰ} ({\vec{r}}_{⊥}, t) = \sum_{n = 0}^{\infty} \frac{i^{n}}{n!} {\frac{d^{n} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n}} |}_{ω_{0}} \frac{\partial^{n} \vec{ℰ} ({\vec{r}}_{⊥}, z, t)}{\partial t^{n}} .

Now, we can rewrite the electric field and related physical quantities in terms of their spatial Fourier transforms, namely,

\vec{ℰ} ({\vec{r}}_{⊥}, z, t) = ℱ_{β}^{z} [\tilde{\vec{ℰ}} ({\vec{r}}_{⊥}, β, t)] = ℱ_{β}^{z} [\tilde{a} (β - β_{0}, t) \tilde{e} (r_{⊥}, β) {\tilde{\vec{e}}}^{- i ω (β) t}],

where as in the standard derivation (Section 2), ã(β − β₀, t) is a dimensionless amplitude and

\tilde{\vec{e}} ({\vec{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].

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¹⁶s⁻¹ and γ = 3.2 × 10¹³s⁻¹ [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, β₀ = 1.51 × 10⁸m⁻¹ (λ₀ ≈ 363.6nm) and Δβ = 0.15 × 10⁸m⁻¹, where β₀ lies exactly at the ZGVP. The second red segment lies in a slow light regime where β₀ = 3 × 10⁸m⁻¹(λ₀ ≈ 362.9nm) and Δβ = 0.3 × 10⁸m⁻¹. Here β₀ lies at a point where the group velocity is 5.51 × 10⁻⁴c. The green line on the left represents the light line, corresponding to $c / \sqrt{∊_{D}}$ . In the presence of absorption, the dispersion curve (and the ZGVP) is shifted slightly (not shown). (c) The values of T_abs = 1/ℑ(ω₀) corresponding to the mode shown in (b).

Download Full Size | PDF

Once the amplitude ã(β − β₀, t) is found, the electric field can be written as

\begin{array}{l} \vec{ℰ} (\vec{r}, z, t) & = & e^{- i ω_{0} t} ℱ_{β}^{z} [(\int d z^{'} A (z^{'}, t) e^{- i (β - β_{0}) z^{'}}) (\int d z^{″} ℱ_{β}^{z^{″}} [\tilde{\vec{e}} ({\vec{r}}_{⊥}, β)] e^{- i β z^{″}})] \\ = & e^{- i ω_{0} t} \int d ζ A (z - ζ, t) e^{i β_{0} ζ} ℱ_{β}^{ζ} [\tilde{\vec{e}} (r_{⊥}, β)], \end{array}

where ω₀ = ω(β₀), where

A (z, t) = 2 π e^{- β_{0} z} ℱ_{β}^{z} [a (β - β_{0}, t)],

a (β - β_{0}, t) = \tilde{a} (β - β_{0}, t) e^{- i [ω (β) - ω_{0}] t},

and where

\tilde{\vec{e}} ({\vec{r}}_{⊥}, β) = ({\tilde{e}}_{x}, {\tilde{e}}_{y}, {\tilde{e}}_{z})

satisfies the vector Helmholtz equation

\begin{array}{l} [(\frac{\partial^{2} {\tilde{e}}_{y}}{\partial x \partial y} - \frac{\partial^{2} {\tilde{e}}_{x}}{\partial y^{2}} + β^{2} {\tilde{e}}_{x} + i β \frac{\partial^{2} {\tilde{e}}_{z}}{\partial x}) \hat{x} + (\frac{\partial^{2} {\tilde{e}}_{x}}{\partial x \partial y} - \frac{\partial^{2} {\tilde{e}}_{y}}{\partial x^{2}} + β^{2} {\tilde{e}}_{y} + i β \frac{\partial^{2} {\tilde{e}}_{z}}{\partial y}) \hat{y} \\ + (i β (\frac{\partial {\tilde{e}}_{x}}{\partial x} + \frac{\partial {\tilde{e}}_{y}}{\partial y}) - \frac{\partial^{2} {\tilde{e}}_{x}}{\partial x^{2}} - \frac{\partial^{2} {\tilde{e}}_{y}}{\partial y^{2}}) \hat{z}] + μ_{0} ∊ ({\vec{r}}_{⊥}, β) ω^{2} \tilde{\vec{e}} (\vec{r}, β) = 0 . \end{array}

For simplicity, we omitted in Eq. (25) the arguments in the components of

\tilde{\vec{e}} ({\vec{r}}_{⊥}, β)

.

Substituting Eqs. (20) and (24) in Eq. (19), and transforming into momentum space gives

\begin{array}{l} [(\frac{\partial^{2} {\tilde{e}}_{y}}{\partial x \partial y} - \frac{\partial^{2} {\tilde{e}}_{x}}{\partial y^{2}} + β^{2} {\tilde{e}}_{x} + i β \frac{\partial^{2} {\tilde{e}}_{z}}{\partial x}) \hat{x} + (\frac{\partial^{2} {\tilde{e}}_{x}}{\partial x \partial y} - \frac{\partial^{2} {\tilde{e}}_{y}}{\partial x^{2}} + β^{2} {\tilde{e}}_{y} + i β \frac{\partial^{2} {\tilde{e}}_{z}}{\partial y}) \hat{y} \\ + (i β (\frac{\partial {\tilde{e}}_{x}}{\partial x} + \frac{\partial {\tilde{e}}_{y}}{\partial y}) - \frac{\partial^{2} {\tilde{e}}_{x}}{\partial x^{2}} - \frac{\partial^{2} {\tilde{e}}_{y}}{\partial y^{2}}) \hat{z}] a (β - β_{0}, t) + μ_{0} ω_{0}^{2} ∊ ({\vec{r}}_{⊥}, ω_{0}) \tilde{\vec{e}} ({\tilde{\vec{r}}}_{⊥}, β) a (β - β_{0}, t) \\ + i μ_{0} \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) (2 ω_{0} ∊ ({\vec{r}}_{⊥}) + ω_{0}^{2} \frac{d ∊ ({\vec{r}}_{⊥}, ω)}{d ω}) \frac{\partial a (β - β_{0}, t)}{\partial t} \\ + μ_{0} \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) \sum_{n = 2}^{\infty} i^{n} (\frac{1}{(n - 2)!} {\frac{d^{n - 2} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 2}} |}_{ω_{0}} + \frac{2}{(n - 1)!} ω_{0} {\frac{d^{n - 1} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 1}} |}_{ω_{0}} \\ + \frac{1}{n!} ω_{0}^{2} {\frac{d^{n} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n}} |}_{ω_{0}}) \frac{\partial^{n} a (β - β_{0}, t)}{\partial t^{n}} + e^{i ω_{0} t} \frac{\partial^{2}}{\partial t^{2}} ℱ_{z}^{β} [R_{NL} ({\vec{r}}_{⊥}, z, t) \vec{ℰ} ({\vec{r}}_{⊥}, z, t)] = 0 . \end{array}

Substituting Eq. (25) in Eq. (26) allows us to replace the first two terms by

∊ ({\vec{r}}_{⊥}, ω) ω^{2} a (β - β_{0}, t) \tilde{\vec{e}} ({\vec{r}}_{⊥}, β)

. We now get

\begin{array}{l} (∊ ({\vec{r}}_{⊥}, ω) ω^{2} - ∊ ({\vec{r}}_{⊥}, ω_{0}) ω_{0}^{2}) a (β - β_{0}, t) \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) \\ - i \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) (2 ω_{0} ∊ ({\vec{r}}_{⊥}, ω) + ω_{0}^{2} {\frac{d ∊ ({\vec{r}}_{⊥}, ω)}{d ω} |}_{ω_{0}}) \frac{\partial a (β - β_{0}, t)}{\partial t} \\ - \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) \sum_{n = 2}^{\infty} i^{n} (\frac{1}{(n - 2)!} {\frac{d^{n - 2} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 2}} |}_{ω_{0}} + \frac{2}{(n - 1)!} ω_{0} {\frac{d^{n - 1} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 1}} |}_{ω_{0}} \\ + \frac{1}{n!} ω_{0}^{2} {\frac{d^{n} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n}} |}_{ω_{0}}) \frac{\partial^{n} a (β - β_{0}, t)}{\partial t^{n}} - e^{i ω_{o} t} \frac{\partial^{2}}{\partial t^{2}} ℱ_{z}^{β} [R_{NL} ({\vec{r}}_{⊥}, z, t) \vec{ℰ} ({\vec{r}}_{⊥}, z, t)] = 0 . \end{array}

In order to remove the dependence on the transverse coordinate, we now multiply by ${\tilde{\vec{e}}}^{*} ({\vec{r}}_{⊥}, β)$ and integrate over r⃗_⊥. Then, Eq. (27) becomes

\begin{array}{l} \frac{\partial a (β - β_{0}, t)}{\partial t} + \frac{i}{ω_{0} U (β, β_{0})} [\int (∊ ({\vec{r}}_{⊥}, ω) ω^{2} - ∊ ({\vec{r}}_{⊥}, ω_{0}) ω_{0}^{2}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}] (β - β_{0}, t) \\ + \frac{i}{ω_{0} U (β, β_{0})} \sum_{n = 2}^{\infty} i^{n} [\int (\frac{1}{(n - 2)!} {\frac{d^{n - 2} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 2}} |}_{ω_{0}} + \frac{2}{(n - 1)!} ω_{0} {\frac{d^{n - 1} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 1}} |}_{ω_{0}} \\ + \frac{1}{n!} ω_{0}^{2} {\frac{d^{n} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n}} |}_{ω_{0}}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}] \frac{\partial^{n} a (β - β_{0}, t)}{\partial t^{n}} \\ + \frac{i e^{i ω_{0} t}}{2 ω_{0} U (β, β_{0})} \frac{\partial^{2}}{\partial t^{2}} \int ℱ_{z}^{β} [R_{NL} ({\vec{r}}_{⊥}, z, t) \vec{ℰ} ({\vec{r}}_{⊥}, z, t)] {\tilde{\vec{e}}}^{*} ({\vec{r}}_{⊥}, β) d {\vec{r}}_{⊥} = 0, \end{array}

where

U (β, β_{0}) = \int (∊ {\vec{r}}_{⊥}, ω) + \frac{1}{2} ω {\frac{∊ ({\vec{r}}_{⊥}, ω)}{d ω} |}_{ω (β_{0})} {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥} .

Note that the terms in parentheses are evaluated at ω₀ ≡ ω(β₀) (rather than the other way around in the standard approach), but

\tilde{\vec{e}} ({\vec{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(β₀, β₀) is related to the energy density of the carrier mode β₀ (see Appendix A), but U(β, β₀) 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 β² in terms of ω in the standard approach (7). Thus, Eq. (28) becomes

\begin{array}{l} \frac{\partial a (β, β, t)}{\partial t} & + i \sum_{m = 0}^{\infty} {\tilde{D}}_{m} {(β - β_{0})}^{m} - a (β - β_{0}, t) + \sum_{n = 2}^{\infty} \sum_{m = 0}^{\infty} i^{n + 1} {\tilde{D}}_{n, m} {(β - β_{0})}^{m} \frac{\partial^{n} a (β - β_{0}, t)}{\partial t^{n}} \\ \frac{i e^{i ω_{0} t}}{2 ω_{0} U (β, β_{0})} \frac{\partial^{2}}{\partial t^{2}} \int ℱ_{z}^{β} [R_{NL} ({\vec{r}}_{⊥}, z, t) \vec{ℰ} ({\vec{r}}_{⊥}, z, t)] {\tilde{\vec{e}}}^{*} ({\vec{r}}_{⊥}, β) d {\vec{r}}_{⊥} = 0, \end{array}

where

{\tilde{D}}_{m} = \frac{1}{m!} \frac{1}{2 ω_{0}} \frac{d^{m}}{d β^{m}} {[\frac{\int [ω^{2} (β) ∊ ({\vec{r}}_{⊥}, ω (β)) - ω_{0}^{2} ∊ ({\vec{r}}_{⊥}, ω)] {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}}{U (β, β_{0})}] |}_{β_{0}},

and

\begin{array}{l} {\tilde{D}}_{n, m} & = \frac{1}{m!} \frac{1}{2 ω_{0}} \frac{d^{m}}{d β^{m}} [\frac{1}{U (β, β_{0})} \int (\frac{1}{(n - 2)!} \frac{d^{n - 2} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 2}} + \frac{2}{(n - 1)!} ω_{0} \frac{d^{n - 1} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n - 1}} \\ {+ \frac{1}{n!} ω_{0}^{2} \frac{d^{n} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{n}}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}] |}_{ω_{0}} . \end{array}

Note that both D̃_m and D̃_m,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 (

\tilde{\vec{e}} ({\vec{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

\begin{array}{l} i {\tilde{D}}_{2, 0} \frac{A (z, t)}{\partial t^{2}} + \frac{\partial A (z, t)}{\partial t} + ω_{0}^{'} \frac{\partial A (z, t)}{\partial z} - i {\tilde{D}}_{2} \frac{\partial^{2} A (z, t)}{\partial z^{2}} \\ + \sum_{m = 3}^{\infty} {(- i)}^{m - 1} {\tilde{D}}_{m} \frac{\partial^{m} A (z, t)}{\partial z^{m}} + \sum_{n = 2}^{\infty} \sum_{\begin{array}{l} m = 0 \\ m \neq 0 if n = 2 \end{array}}^{\infty} {(- 1)}^{m} i^{n + m + 1} {\tilde{D}}_{n, m} \frac{\partial^{n + m} A (z, t)}{\partial t^{n} \partial z^{m}} \\ + i \frac{e^{i ω_{0} t}}{2 ω_{0}} ℱ_{β}^{z} [\frac{1}{U (β, β_{0})} \frac{\partial^{2}}{\partial t^{2}} \int ℱ_{z}^{β} [R_{NL} ({\vec{r}}_{⊥}, z, t) \vec{ℰ} ({\vec{r}}_{⊥}, z, t)] {\tilde{\vec{e}}}^{*} ({\vec{r}}_{⊥}, β) d {\vec{r}}_{⊥}] = 0, \end{array}

where

{\tilde{D}}_{2} = \frac{ω_{0}^{″}}{2} + {\tilde{D}}_{2, 0} ω_{0}^{' 2},

and

\begin{array}{l} {\tilde{D}}_{2, 0} = \frac{2}{2 ω_{0} U (β, β_{0})} \int (2 ∊ ({\vec{r}}_{⊥}, ω_{0}) + 4 ω_{0} {\frac{d ∊ ({\vec{r}}_{⊥}, ω)}{d ω} |}_{ω_{0}} \\ {+ ω_{0}^{2} \frac{d^{2} ∊ ({\vec{r}}_{⊥}, ω)}{d ω^{2}} |}_{ω_{0}}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥} . \end{array}

In Eqs. (34)–(35), we used Eqs. (31)–(32), where it follows from the chain rule that

{\frac{d ∊ ({\vec{r}}_{⊥}, ω)}{d β} |}_{β_{0}} = {\frac{d ∊ ({\vec{r}}_{⊥}, ω)}{d ω} |}_{ω_{0}} {\frac{d ω}{d β} |}_{β_{0}}

. The coefficients ω′₀ and ω″₀ stand for

{\frac{d ω}{d β} |}_{β_{0}}

and

{\frac{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 momentum-based 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 D̃₂, 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 β₀ in the current approach, similarly to Eq. (14), defined as

v_{g} \equiv ℜ (ω_{0}^{'}) .

For materials with negligible absorption, the second term in Eq. (34) is proportional to

v_{g}^{2}

, in analogy to the GVD coefficient in the standard derivation (where the second term is proportional to

v_{g}^{- 2}

, 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_{g}^{- m}$ is replaced by a series of terms that scale with $v_{g}^{m}$ , hence, is rapidly converging. In addition, the high-order derivatives of ω, are also small in the slow light regime.

3.2. Dimensional analysis

We now assume that the nonlinear response is spatially-local, i.e., that R_NL represents a cubic (Kerr-like) response, R_NL = ∊_NL|ℰ⃗|²; Furthermore, U (29) is evaluated at β₀ 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)² scale change due to the propagation; this term is comparable to D̃_2,0ω′₀² (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,

z^{*} = z - v_{g} t,

t^{*} = t .

Accordingly, the derivatives of the coordinates transform as

\frac{\partial}{\partial z} = \frac{\partial}{\partial z^{*}},

\frac{\partial}{\partial t} = \frac{\partial}{\partial t^{*}} - v_{g} \frac{\partial}{\partial z^{*}} .

Then, Eq. (33) becomes

\begin{array}{l} i {\tilde{D}}_{2, 0} \frac{\partial^{2} A (z^{*}, t^{*})}{\partial t^{* 2}} + \frac{\partial A (z^{*}, t^{*})}{\partial t^{*}} + i ℑ (ω_{0}^{'}) \frac{\partial A (z^{*}, t^{*})}{\partial z^{*}} \\ - [\frac{i}{2} ω_{0}^{″} - (2 ℑ (ω_{0}^{'}) v_{g} + i ℑ {(ω_{0}^{'})}^{2}) {\tilde{D}}_{2, 0}] \frac{\partial^{2} A (z^{*}, t^{*})}{\partial z^{* 2}} \\ - 2 i v_{g} {\tilde{D}}_{2, 0} \frac{\partial^{2} A (z^{*}, t^{*})}{\partial z^{*} \partial t^{*}} + \sum_{m = 3}^{\infty} {(- i)}^{m - 1} {\tilde{D}}_{m} \frac{\partial^{m} A (z^{*}, t^{*})}{\partial z^{* m}} \\ + \sum_{n = 2}^{\infty} \sum_{\begin{array}{l} m = 0 \\ m \neq 0 if n = 2 \end{array}}^{\infty} \sum_{k = 0}^{n} (\begin{matrix} n \\ k \end{matrix}) i^{n + 1} {(- i)}^{m} {(- v_{g})}^{k} {\tilde{D}}_{n, m} \frac{\partial^{n + m} A (z^{*}, t^{*})}{\partial t^{* n - k} \partial z^{* m + k}} \\ + \frac{i}{T_{NL}} {| A (z^{*}, t^{*}) |}^{2} A (z^{*}, t^{*}) = 0, \end{array}

where

T_{NL} = \frac{2 U (β_{0}, β_{0})}{∊_{NL} ω_{0} \int {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{4} d {\vec{r}}_{⊥}},

is the characteristic time for the nonlinear effect to become significant. For lossy media, ω′₀ is complex, so that the ∂/∂t term does not vanish, but rather, leaves behind its imaginary part, ℑ(ω′₀) note that it does not represent a velocity, but rather the dispersion of the absorption. Similarly, the D̃_2,0ω′₀² term does not cancel completely either.

In order to reveal the hierarchy of the terms in Eq. (39), we normalize the coordinates as follows

\tilde{z} = \frac{z^{*}}{Z},

\tilde{t} = \frac{t^{*}}{T_{disp}},

where Z is the initial pulse width and T_disp = Z²/ω″₀ is the dispersion time, the analogue of the dispersion length L_disp. Then, Eq. (39) obtains the dimensionless form

\begin{array}{l} if \frac{\partial^{2} A (\tilde{z}, \tilde{t})}{\partial {\tilde{t}}^{2}} + \frac{\partial A (\tilde{z}, \tilde{t})}{\partial \tilde{t}} - 2 i g \frac{\partial^{2} A (\tilde{z}, \tilde{t})}{\partial \tilde{z} \partial \tilde{t}} + i q \frac{\partial A (\tilde{z}, \tilde{t})}{\partial \tilde{z}} - \frac{i}{2} \frac{\partial^{2} A (\tilde{z}, \tilde{t})}{\partial {\tilde{z}}^{2}} \\ + s \frac{\partial^{2} A (\tilde{z}, \tilde{t})}{\partial {\tilde{z}}^{2}} + \sum_{m = 3}^{\infty} {(- i)}^{m - 1} T_{disp} \frac{{\tilde{D}}_{m}}{Z^{m}} \frac{\partial^{m} A (\tilde{z}, \tilde{t})}{\partial {\tilde{z}}^{m}} \\ + \sum_{n = 2}^{\infty} \sum_{\begin{array}{l} m = 0 \\ m \neq 0 if n = 2 \end{array}}^{\infty} \sum_{k = 0}^{n} (\begin{matrix} n \\ k \end{matrix}) i^{n + 1} {(- i)}^{m} {(\frac{(- v_{g})}{Z})}^{k} \frac{{\tilde{D}}_{n, m}}{T_{disp}^{n - k - 1}} \frac{1}{Z^{m}} \frac{\partial^{n + m} A (\tilde{z}, \tilde{t})}{\partial {\tilde{t}}^{n - k} \partial {\tilde{z}}^{m + k}} \\ + \frac{T_{disp}}{T_{NL}} {| A (\tilde{z}, \tilde{t}) |}^{2} A (\tilde{z}, \tilde{t}) = 0, \end{array}

where the dimensionless coefficients are

f \equiv \frac{{\tilde{D}}_{2, 0}}{T_{disp}},

g \equiv \frac{v_{g} {\tilde{D}}_{2, 0}}{Z},

q \equiv \frac{ℑ (ω_{0}^{'}) T_{disp}}{Z} = \frac{ℑ (ω_{0}^{'})}{ω_{0}^{″}} Z,

s \equiv \frac{ℑ (ω_{0}^{'})}{ω_{0}^{″}} {\tilde{D}}_{2, 0} (2 v_{g} + i ℑ (ω_{0}^{'})) .

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, D̃_2,0 can be estimated as 1/ω₀ (see Appendix B). Then, f ∼ T₀/T_disp and g ∼ T₀/T_P where we defined the period T₀ = 2π/ω₀, and the pulse duration, T_P, which scales as Z/v_g. Thus, in the typical case where T₀ ≪ 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 D̃_m and D̃_n,m (see Appendix B), one can see that the seventh and eighth terms scale as (T_disp/T₀) (β₀Z)^−m and (T_disp/T₀) (β₀Z)^−m (v_gT_disp/Z)^k(T_disp/T₀)⁻ⁿ, respectively. Typically, T_disp/T₀ is much greater than unity; β₀Z and (v_gT_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₀ is larger than v_gT_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(β₀, β₀) (rather than as P₀, as in the standard derivation, see Eq. (40)). In lossy media, T_NL becomes complex as it contains the complex values U(β₀, β₀) and ω₀. Overall, the nonlinear term becomes important in cases where T_NL is comparable or shorter than T_disp.

If all mixed terms are indeed small, Eq. (39) becomes

\frac{\partial A (z, t)}{\partial t} - i \frac{ω_{0}^{″}}{2} \frac{\partial^{2} A (z, t)}{\partial z^{2}} + \sum_{m = 3}^{\infty} {(- i)}^{m - 1} {\tilde{D}}_{m} \frac{\partial^{m} A (z, t)}{\partial z^{m}} + \frac{i}{T_{NL}} {| A (z, t) |}^{2} A (z, t) = 0,

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]).

3.3. Scaling with the group velocity

Many previous studies showed that linear effects such as absorption scale with $v_{g}^{- 1}$ and that cubic nonlinear effects scale with $v_{g}^{- 2}$ , 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₀ (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–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).

4. 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 $a (β, t = 0) = C e^{- {(\frac{β - β_{0}}{Δ β})}^{2}}$ where β₀ 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,

∊_{M} (ω) = ∊_{\infty} - \frac{ω_{P}^{2}}{ω^{2} + i γ ω},

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_{g}^{m}$ 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–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ω₀t in Eq. (22) also describes the attenuation of the pulse in time as ω₀ consists of an imaginary part. Therefore, in the following plots, we will present the value |A(z, t)e^−iω₀t| instead of just the amplitude, A(z, t). In that sense, the characteristic time for absorption is T_abs = 1/ℑ(ω₀).

4.1. 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)} (β) \equiv \sum_{m = 1}^{M} {\tilde{D}}_{m} {(β - β_{0})}^{m}$ , such that the solution of Eq. (44) is

A (z, t) = 2 π e^{- i β z_{0}} ℱ_{β}^{z} [a (β - β_{0}, t = 0) e^{- i Ω^{M} (β) t}] .

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 β₀ was chosen to be at the ZGVP, see Fig. 1(b). The spectrum width is chosen to be Δβ ≈ 0.1β₀, 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

R_{M, M_{\max}} (t) = \frac{\int {| A^{(M)} (z, t) - A^{(M_{\max})} (z, t) |}^{2} d z}{\int {| A^{(M_{\max})} (z, t) |}^{2} d z},

where A^(M_max) 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.

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, β₀ 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 =∞.

Download Full Size | PDF

In the following simulations, we chose a narrower spectral regime such that Δβ = 0.01β₀. 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.

Fig. 3 Spatio-temporal pulse dynamics of the amplitude of the electric field (color represents the value of |A(z, t)|e^{−ℑ(ω₀)t}) in the linear case. Plots (a) and (b) represent the ZGVP regime case in the absence and presence of absorption, respectively. Here, $β_{0} = 1.51 \times 10^{8} \frac{1}{m}$ , $Δ β = 0.015 \times 10^{8} \frac{1}{m}$ , T_disp = 7.7 × 10⁻¹⁰s; in (b), T_abs = 6.3 × 10⁻¹⁴s, such that T_abs ≪ T_disp. Plots (c) and (d) represent the slow light regime case in the absence and presence of absorption, respectively. Here, $β_{0} = 3 \times 10^{8} \frac{1}{m}$ , $Δ β = 0.03 \times 10^{8} \frac{1}{m}$ , T_disp = 1.39 × 10⁻⁹s , T_abs = 6.26 × 10⁻¹⁴s in (d) such that T_abs ≪ T_disp. In both cases, the total time of propagation is T_disp and 10T_abs, respectively.

Download Full Size | PDF

Fig. 4 The spatial spectrum associated with the results shown in Figs. 3(a)–(d).

Download Full Size | PDF

In order to show the validity of Eq. (44) in the slow light regime, we repeated the process above for β₀ centered at the slow light regime where v_g = 5.5 × 10⁻⁴c (see Fig. 1(b)). Here, as earlier, we choose Δβ ≈ 0.1β₀. However, although Δβ here is nominally wider than the spectral width chosen in the case of the ZGVP (since the value of β₀ 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_{g}^{- m}$ in the m th-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.

4.2. 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_NL; in the presence of absorption, we replace T_NL → ℜ[T_NL].

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_NL ≪ 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].

Fig. 5 Spatio-temporal dynamics of the amplitude of the electric field (color represents the value |A(z, t)|e^{−ℑ(ω₀)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_NL = 0.2T_disp in the absence of absorption (a) and T_NL = 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_NL = 0.2T_disp in the absence of absorption (c) and T_NL = 0.1T_abs when absorption is present (d). In both cases, the total time of propagation is T_NL and 10T_abs, respectively.

Download Full Size | PDF

When considering absorption, in order to observe a significant nonlinear effect, we choose T_NL to be much shorter than T_abs, thus, also much shorter than T_NL in the absence of absorption. If T_NL 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. 6(b) to Fig. 4(b), the spectrum is broadened very quickly due to self-phase modulation.

Fig. 6 The spatial spectrum associated with the results shown in Figs. 5(a)–(d).

Download Full Size | PDF

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_NL such that nonlinear effects are stronger than the dispersion and omit the absorption. Such configuration results in T_NL ≪ 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. Fig. 5(b) and (d) (as well as in Fig. 6(b) and (d)) show the dynamics in the presence of absorption. 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_NL in both cases is similar, and also similar to cases where the mode has a “normal” velocity, see discussion in Section 3.3.

5. 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 to 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–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–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–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–45], spin waves [46], water waves [47], quantum/matter waves [48,49] etc..

Appendix A Finding the energy density of the carrier mode

We start with Eq. (29) for β = β₀, namely,

\begin{array}{l} U (β, β_{0}) & = & s (∊ ({\vec{r}}_{⊥}, ω_{0}) + \frac{1}{2} ω_{0} {\frac{d ∊ ({\vec{r}}_{⊥}, ω)}{d ω} |}_{ω_{0}}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d {\vec{r}}_{⊥}, \\ = & \frac{1}{2} \int (∊ ({\vec{r}}_{⊥}, ω_{0}) + {\frac{d [∊ ({\vec{r}}_{⊥}, ω)]}{d ω} |}_{ω_{0}}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d {\vec{r}}_{⊥} . \end{array}

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, $\int ∊ ({\vec{r}}_{⊥}, ω_{0}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} = \int μ {| \tilde{\vec{h}} ({\vec{r}}_{⊥}, β_{0}) |}^{2}$ , where $\tilde{\vec{h}} ({\vec{r}}_{⊥}, β_{0})$ is the transverse magnetic field of mode β₀, is still valid in the presence of dispersion.

The term $\int ∊ ({\vec{r}}_{⊥}, ω_{0}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d x d y$ in Eq. (48) can be written as

\begin{array}{l} \int ∊ ({\vec{r}}_{⊥}, ω_{0}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d x d y = \int ∊ ({\vec{r}}_{⊥}, ω_{0}) ({| {\tilde{e}}_{x} |}^{2} + {| {\tilde{e}}_{y} |}^{2} + {| {\tilde{e}}_{z} |}^{2}) d x d y \\ = \int ∊ ({\vec{r}}_{⊥}, ω_{0}) (\frac{1}{ω_{0} ∊ ({\vec{r}}_{⊥}, ω_{0})} (- \frac{\partial h_{z}}{\partial y} + i β_{0} h_{y}) {\tilde{e}}_{x}^{*} + \frac{1}{ω_{0} ∊ ({\vec{r}}_{⊥}, ω_{0})} (\frac{\partial h_{z}}{\partial x} - i β_{0} h_{x}) {\tilde{e}}_{y}^{*} \\ + \frac{1}{ω_{0} ∊ ({\vec{r}}_{⊥}, ω_{0})} (\frac{\partial h_{y}}{\partial x} - \frac{\partial h_{x}}{\partial y}) {\tilde{e}}_{z}^{*}) d x d y \\ = \frac{1}{ω_{0}} \int (- \frac{\partial h_{z}}{\partial y} {\tilde{e}}_{x}^{*} + i β_{0} h_{y} {\tilde{e}}_{x}^{*} + \frac{\partial h_{z}}{\partial x} {\tilde{e}}_{y}^{*} - i β_{0} h_{x} {\tilde{e}}_{y}^{*} + \frac{\partial h_{y}}{\partial x} {\tilde{e}}_{z}^{*} - \frac{\partial h_{x}}{\partial y} {\tilde{e}}_{z}^{*}) d x d y, \end{array}

where we have used Maxwell-Ampéré law,

{\tilde{e}}_{x} = \frac{i}{ω_{0} ∊ ({\vec{r}}_{⊥}, ω_{0})} (i β_{0} h_{y} - \frac{\partial h_{z}}{\partial y}),

{\tilde{e}}_{y} = \frac{i}{ω_{0} ∊ ({\vec{r}}_{⊥}, ω_{0})} (\frac{\partial h_{x}}{\partial y} - i β_{0} h_{y}),

{\tilde{e}}_{z} = \frac{i}{ω_{0} ∊ ({\vec{r}}_{⊥}, ω_{0})} (\frac{\partial h_{x}}{\partial y} - \frac{\partial h_{y}}{\partial x}) .

Integrating by parts gives

\begin{array}{l} \int ∊ ({\vec{r}}_{⊥}, ω_{0}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d x d y = \\ \frac{1}{ω_{0}} \int (h_{z} \frac{\partial {\tilde{e}}_{x}^{*}}{\partial y} + i β_{0} h_{y} {\tilde{e}}_{x}^{*} - h_{z} \frac{\partial {\tilde{e}}_{y}^{*}}{\partial x} + i β_{0} h_{y} {\tilde{e}}_{x}^{*} - h_{z} \frac{\partial {\tilde{e}}_{z}^{*}}{\partial x} + h_{x} \frac{\partial {\tilde{e}}_{z}^{*}}{\partial y}) d x d y, \end{array}

where we used

\int h_{r_{i}} {\tilde{e}}_{r_{j}} {d r_{k} |}_{r_{i} = - \infty}^{r_{i} = \infty} = 0,

where i, j, k represent any spatial coordinate. We now substitute Faraday-Henri Law

\frac{\partial {\tilde{e}}_{z}}{\partial x} = i β_{0} {\tilde{e}}_{x} + i ω_{0} μ_{0} h_{y},

\frac{\partial {\tilde{e}}_{z}}{\partial x} = i β_{0} {\tilde{e}}_{x} + i ω_{0} μ_{0} h_{y},

in Eq. (51) and use Eq. (52) again. Eq. (51) becomes

\begin{array}{l} \int ∊ ({\vec{r}}_{⊥}, ω_{0}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d x d y = \\ \frac{1}{ω_{0}} \int (h_{z} (\frac{\partial {\tilde{e}}_{x}^{*}}{\partial y} - \frac{\partial {\tilde{e}}_{y}^{*}}{\partial x}) + i ω_{0} μ_{0} ({| h_{x} |}^{2} + {| h_{y} |}^{2})) d x d y . \end{array}

In the final stage, we substitute

\frac{\partial {\tilde{e}}_{x}^{*}}{\partial y} - \frac{\partial {\tilde{e}}_{y}^{*}}{\partial x} = i ω_{0} μ_{0} h_{z}^{*},

and get

\begin{array}{l} \int ∊ ({\vec{r}}_{⊥}, ω_{0}) {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d x d y & = μ_{0} \int ({| h_{x} |}^{2} + {| h_{y} |}^{2} + {| h_{z} |}^{2}) d x d y \\ = μ_{0} \int {| \vec{h} ({\vec{r}}_{⊥}, β_{0}) |}^{2} d x d y . \end{array}

This proves our claim. Substituting Eq. (56) in Eq. (48) gives

U (β, β_{0}) = \frac{1}{2} \int ({\frac{d (∊ ({\vec{r}}_{⊥}, ω_{0}) ω)}{d ω} |}_{ω_{0}} {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β_{0}) |}^{2} + μ_{0} {| \vec{h} ({\vec{r}}_{⊥}, β_{0}) |}^{2}) d x d y .

Since

| \tilde{\vec{e}} ({\vec{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.

Appendix B Estimate of D̃_n,m and D̃_m

In this Appendix, we estimate the magnitude of the coefficients D̃_n,m (32) and D̃_m (31).

Starting with index n, when n increases, according to Eq. (32), the relative magnitude of consecutive terms in n, namely, D̃_n+1,m/D̃_n,m, can be estimated by

R_{n} ~ \frac{{\int \frac{d ∊ ({\vec{r}}_{⊥}, ω_{0}) ({\vec{r}}_{⊥}, ω)}{d ω} |}_{ω_{0}} d {\vec{r}}_{⊥}}{\int {∊ ({\vec{r}}_{⊥}, ω) |}_{ω_{0}} d {\vec{r}}_{⊥}} .

For a Drude metal, when taking into account Eq. (45) and ignoring the details of the waveguide width, we can estimate,

R_{n} ~ \frac{1}{ω_{0}}

. For a dispersive dielectric, R_n will be smaller, hence, the following analysis can be regarded as an estimate for an upper limit. It follows that

{\tilde{D}}_{n, m} ~ \frac{1}{ω_{0}^{n - 1}} {\tilde{D}}_{1, m} .

Continuing with index m, according to Eq. (32), one can evaluate

R_{m} ~ \frac{\int {\frac{d {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2}}{d β} |}_{β_{0}} d {\vec{r}}_{⊥}}{\int {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}},

where R_m is the ratio of order of magnitude of two adjacent terms of D̃_n,m in m, specifically, D̃_n,1/D̃_n,0. The magnitude of R_m may vary when considering different waveguides, however, in this work it will be reasonable to consider a plasmonic waveguide where

| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) | ~ e^{- | {\vec{k}}_{⊥} \cdot {\vec{r}}_{⊥} |},

| {\vec{k}}_{⊥} | = \sqrt{β^{2} - ∊ ({\vec{r}}_{⊥}, β) {(ω / c)}^{2}} .

Substituting Eqs. (61) in Eq. (60) results in

R_{m} ~ \frac{{\int {\vec{r}}_{⊥} | \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}}{\int {| \tilde{\vec{e}} ({\vec{r}}_{⊥}, β) |}^{2} d {\vec{r}}_{⊥}} .

Notice that Eq. (62) is not a direct result of the substitution but rather a term which maintains the magnitude in Eq. (60). It is deduced from Eq. (62) that R_m ∼ 〈 r⃗_⊥(β)〉, where 〈r⃗_⊥(β)〉 is the characteristic width of the waveguide mode

\tilde{\vec{e}} ({\vec{r}}_{⊥}, β)

.

In cases where β₀ ∼ |k⃗_⊥| at any of the constituents, 〈r_⊥(β)〉 scale as 1/β₀ as β₀ is the dominant value in Eq. (61b). If we consider small values of β₀ with respect to |k⃗_⊥|, due to the nature of the metal cladding in a MDM waveguide, 〈r⃗_⊥(β)〉 cannot exceed a range which is much larger than the core, thus R_m ∼ max(〈r⃗_⊥(β)〉, 1/β₀).

Considering the above, ${\tilde{D}}_{n, 1} ~ \frac{{\tilde{D}}_{n, 0}}{β_{0}}$ , where similar analysis in the high-order terms gives,

{\tilde{D}}_{n, m} ~ \frac{{\tilde{D}}_{n, 0}}{β_{0}^{m}} .

Since D̃_1,0 ≡ 1, according to Eqs. (59) and (63), the magnitude of D̃_n,m is

\frac{1}{ω_{0}^{n - 1} β_{0}^{m}}

. For small values of β₀ the scale of D̃_n,m is even smaller as β₀ < |k⃗_⊥|.

The structure of D̃_m is similar of D̃_n,m. Considering the same arguments as above, it can be shown that the magnitude of D̃_m is $ω_{0} / β_{0}^{m}$ , where it can be smaller with large β₀.

Funding

People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no. 333790.

Acknowledgements

The authors would like to acknowledge useful discussions with I. Bar and P. Y. Chen, N. Gavish and M. Spector. The authors would like to acknowledges financial support from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007–2013) under REA grant agreement no. 333790.

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). [CrossRef]

3. P. Siddons, “Light propagation through atomic vapours,” J. Phys. B: At. Mol. Opt. Phys 47, 093001 (2014). [CrossRef]

4. J. B. Khurgin, “Slow light in various media: a tutorial,” Adv. Opt. Photonics 2, 287–318 (2010). [CrossRef]

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). [CrossRef]

6. T. Baba, “Slow light in photonic crystals,” Nat. Photonics 2, 465–473 (2008). [CrossRef]

7. T. F. Krauss, “Slow light in photonic crystal waveguides,” J. Phys. D 40, 2666–2670 (2007). [CrossRef]

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). [CrossRef] [PubMed]

9. R. Won, “Quantum memory: Extended storage times,” Nat. Photonics 7, 677 (2013). [CrossRef]

10. G. Brennen, E. Giacobino, and C. Simon, “Focus on quantum memory,” New J. Phys. 17, 050201 (2007). [CrossRef]

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). [CrossRef]

12. M. Notomi, E. Kuramochi, and T. Tanabe, “Large-scale arrays of ultrahigh-q coupled nanocavities,” Nat. Photonics 2, 741–747 (2008). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef]

16. T. Kamalakis, “Derivation of coupled envelope propagation equations in silicon-based photonic crystal waveguides,” J. Opt. Soc. Am. B 33, 2339–2349 (2016). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

23. A. Zadok, A. Eyal, and M. Tur, “Stimulated brillouin scattering slow light in optical fibers,” Appl. Opt. 50, E38–E49 (2011). [CrossRef]

24. M. Afzelius, N. Gisin, and H. de Riedmatten, “Quantum memory for photons,” Phys. Today 68, 42–47 (2015). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef] [PubMed]

36. A. F. Page, “Phd thesis,” Ph.D. thesis, Imperial CollegeLondon (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). [CrossRef] [PubMed]

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). [CrossRef]

39. N. Sun, J. Chen, and D. Tang, “Stopping light in an active medium,” The Eur. Phys. J. D 69, 219 (2015). [CrossRef]

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). [CrossRef]

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). [CrossRef]

42. U. Leonhardt, “A laboratory analogue of the event horizon using slow light in an atomic medium,” Nature. 415, 406–409 (2002). [CrossRef] [PubMed]

43. M. Fink, “Acoustic time-reversal mirrors: Imaging of complex media with acoustic and seismic waves,” Top. Appl. Phys. 84, 17–43 (2002). [CrossRef]

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). [CrossRef]

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). [CrossRef] [PubMed]

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). [CrossRef]

48. C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases (Cambridge University Press, 2008), 2nd ed. [CrossRef]

49. J. Mendonca and A. Guerreiro, “Time refraction and the quantum properties of vacuum,” Phys. Rev. A 72, 063805 (2005). [CrossRef]

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). [CrossRef] [PubMed]

51. A. Figotin and I. Vitebskiy, “Slow light in photonic crystals,” Waves Random Complex Media 16, 293–382 (2006). [CrossRef]

52. S. P. Shipman and A. T. Welters, “Pathological scattering by a defect in a slow-light periodic layered medium,” J. Math. Phys. 57, 022902 (2016). [CrossRef]

53. G. P. Agrawal, Nonlinear Fiber Optics (Academics, San Diego, 2011).

54. G. Fibich, The nonlinear Schrödinger equation, singular solutions and optical collapse (Springer, Heidelberg, 2015). [CrossRef]

55. F. Poletti and P. Horak, “Description of ultrashort pulse propagation in multimode optical fibers,” J. Opt. Soc. Am. B 25, 1645–1654 (2008). [CrossRef]

56. Y. Sivan, S. Rosenberg, and A. Halstuch, “Coupled-mode theory for electromagnetic pulse propagation in dispersive media undergoing a spatiotemporal perturbation: Exact derivation, numerical validation and peculiar wave mixing,” Phys. Rev. B 93, 144303 (2016). [CrossRef]

57. P. Tran, “Photonic-band-structure calculation of material possessing kerr nonlinearity,” Phys. Rev. B 52, 10673–10676 (1995). [CrossRef]

58. V. Lousse and J. P. Vigneron, “Self-consistent photonic band structure of dielectric superlattices containing nonlinear optical materials,” Phys. Rev. E 63, 027602 (2001). [CrossRef]

59. B. Maes, P. Bienstman, and R. Baets, “Bloch modes and self-localized waveguides in nonlinear photonic crystals,” J. Opt. Soc. Am. B 22, 613–619 (2005). [CrossRef]

60. C. de Sterke and J. Sipe, “Envelope-function approach for the electrodynamics of nonlinear periodic structures,” Phys. Rev. A 38, 5149–5165 (1988). [CrossRef]

61. C. de Sterke and J. Sipe, “Switching dynamics of finite periodic nonlinear media: A numerical study,” Phys. Rev. A 42, 2858–2869 (1990). [CrossRef] [PubMed]

62. J. Sipe, L. Poladian, and C. M. de Sterke, “Propagation through non-uniform grating structures,” J. Opt. Soc. Am. B 11, 1307–1320 (1994). [CrossRef]

63. C. M. de Sterke, K. R. Jackson, and B. D. Robert, “Nonlinear coupled-mode equations on a finite interval: a numerical procedure,” J. Opt. Soc. Am. B 8, 403–412 (1991). [CrossRef]

64. J. T. Mok, C. M. de Sterke, I. C. M. Littler, and B. J. Eggleton, “Dispersionless slow light using gap solitons,” Nat. Phys. 2, 775–780 (2006). [CrossRef]

65. C. de Sterke, D. Salinas, and J. Sipe, “Coupled-mode theory for light propagation through deep nonlinear gratings,” Phys. Rev. E 54, 1969–1989 (1996). [CrossRef]

66. N. Aközbek and S. John, “Optical solitary waves in two- and three-dimensional nonlinear photonic band-gap structures,” Phys. Rev. E 57, 2287–2319 (1998). [CrossRef]

67. N. Bhat and J. Sipe, “Optical pulse propagation in nonlinear photonic crystals,” Phys. Rev. E 64, 056604 (2001). [CrossRef]

68. J. E. Sipe, “Vector k · p approach for photonic band structures,” Phys. Rev. E 62, 5672–5677 (2000). [CrossRef]

69. Y. Sivan and J. B. Pendry, “Time-reversal in dynamically-tuned zero-gap periodic systems,” Phys. Rev. Lett. 106, 193902 (2011). [CrossRef]

70. Y. Sivan and J. B. Pendry, “Theory of wave-front reversal of short pulses in dynamically-tuned zero-gap periodic systems,” Phys. Rev. A 84, 033822 (2011). [CrossRef]

71. M. Reza, M. M. Dignam, and S. Hughes, “Can light be stopped in realistic metamaterials?” Nature 455, E10–E11 (2008). [CrossRef]

72. J. Joannopoulos, S. Johnson, J. Winn, and R. Meade, Photonic Crystals - Molding the Flow of Light (Princeton University Press, 2008), 2nd ed.

73. W. Cai and V. Shalaev, Optical metamaterials, fundamentals and applications (Springer, Heidelberg, 2010).

74. P. Yeh, Optical Waves in Layered Media (Wiley-Interscience, 2005), 2nd ed.

75. P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin, “Light interaction with photonic and plasmonic resonances,” ArXiv: 170502433 Phys. (2017).

76. L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii, Electrodynamic of Continuous Media (Pergamon Press, Kidlington, 2008).

77. S. G. Johnson, Notes on the algebraic structure of wave equations (http://math.mit.edu/stevenj/18.369/wave-equations.pdf, 2010).

78. L. Thévenaz, S. Chin, I. Dicaire, J.-C. Beugnot, S. F. Mafang, and M. G. Herraez, “Experimental verification of the effect of slow light on molecular absorption,” Proceeding SPIE 7503, 75034W (2009). [CrossRef]

79. I. Dicaire, A. D. Rossi, S. Combrié, and L. Thévenaz, “Probing molecular absorption under slow-light propagation using a photonic crystal waveguide,” Opt. Lett. 37, 4934–4936 (2012). [CrossRef] [PubMed]

80. A. Archambault, T. V. Teperik, F. Marquier, and J. J. Greffet, “Surface plasmon fourier optics,” Phys. Rev. B 79, 195414 (2009). [CrossRef]

81. Z. I. Samson, P. Horak, K. F. MacDonald, and N. I. Zheludev, “Femtosecond surface plasmon pulse propagation,” Opt. Lett. 36, 250–252 (2011). [CrossRef] [PubMed]

82. Q. Lin, O. J. Painter, and G. P. Agrawal, “Nonlinear optical phenomena in silicon waveguides: modeling and applications,” Opt. Express 15, 16604–16644 (2007). [CrossRef] [PubMed]

83. C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express 18, 21427–21448 (2010). [CrossRef] [PubMed]

84. R. Atkinson, W. R. Hendren, G. A. Wurtz, W. Dickson, A. V. Zayats, P. Evans, and R. J. Pollard, “Anisotropic optical properties of arrays of gold nanorods embedded in alumina,” Phys. Rev. B 73, 235402 (2006). [CrossRef]

85. G. A. Wurtz, R. J. Pollard, W. R. Hendren, G. P. Wiederrecht, D. J. Gosztola, V. A. Podolskiy, and A. V. Zayats, “Designed ultrafast optical nonlinearity in a plasmonic nanorod metamaterial enhanced by nonlocality,” Nat. Nanotechnol. 6, 107–111 (2011). [CrossRef] [PubMed]

Pulse propagation in the slow and stopped light regime

Abstract

1. Introduction

2. Standard formulation - frequency-based expansion

3. Alternative formulation - momentum-based expansion

3.1. Derivation

3.2. Dimensional analysis

3.3. Scaling with the group velocity

4. Examples

4.1. Simulations of pulse propagation in the slow and stopped light regimes; linear media

4.2. Simulations of pulse propagation in the slow and stopped light regimes; nonlinear media

5. Discussion

Appendix A Finding the energy density of the carrier mode

Appendix B Estimate of D̃_n,m and D̃_m

Funding

Acknowledgements

References and links

Cited By

Figures (6)

Equations (74)

Optics Express

Abstract

1. Introduction

2. Standard formulation - frequency-based expansion

3. Alternative formulation - momentum-based expansion

3.1. Derivation

3.2. Dimensional analysis

3.3. Scaling with the group velocity

4. Examples

4.1. Simulations of pulse propagation in the slow and stopped light regimes; linear media

4.2. Simulations of pulse propagation in the slow and stopped light regimes; nonlinear media

5. Discussion

Appendix A Finding the energy density of the carrier mode

Appendix B Estimate of D̃n,m and D̃m

Funding

Acknowledgements

References and links

Cited By

Figures (6)

Equations (74)

Optics Express

Appendix B Estimate of D̃_n,m and D̃_m