## Abstract

We describe generalized nonlinear envelope equation modeling of sub-cycle dynamics on the underlying electric field carrier during one-dimensional propagation in fused silica. Generalized envelope equation simulations are shown to be in excellent quantitative agreement with the numerical integration of Maxwell’s equations, even in the presence of shock dynamics and carrier steepening on a sub-50 attosecond timescale. In addition, by separating the effects of self-phase modulation and third harmonic generation, we examine the relative contribution of these effects in supercontinuum generation in fused silica nanowire waveguides.

©2007 Optical Society of America

## 1. Introduction

Ultrashort pulse propagation in highly nonlinear waveguides such as photonic crystal fibers (PCF) remains a subject of intense study, motivated in large part by the many applications of supercontinuum (SC) generation [1]. A recent application of particular interest has exploited the initial nonlinear evolution phase of SC generation for pulse compression, and experiments using tapered PCF nanowires have reported impressive results approaching the single-cycle regime [2]. These experiments were modeled used extended nonlinear Schrödinger equation (NLSE) simulations and, although this approach has successfully described octave-spanning SC generation, the modeling of genuinely single- and sub-cycle propagation dynamics is generally considered to require the solution of Maxwell’s equations [3, 4]. In this regard, however, we note that there has been considerable previous development of different nonlinear propagation models that maintain the fundamental accuracy of Maxwell’s equations whilst still possessing computational simplicity comparable to the NLSE. Many of these models are based on adding correction terms to the NLSE modeling the evolution of a pulse envelope [5–8], while other (more accurate) approaches develop propagation equations for the full electric field, modified with the addition of a unidirectional constraint [9–12].

Our objective here is to present a new nonlinear propagation equation that remains expressed within an envelope formalism, yet which nonetheless accurately models propagation dynamics in highly nonlinear waveguides in the sub-cycle regime. We compare the results of this generalized nonlinear envelope equation (GNEE) with the direct integration of Maxwell’s equations for one-dimensional propagation in fused silica, and obtain excellent agreement modeling spectral broadening over bandwidths many times the optical carrier frequency. Although our simulations solve for the evolution of a generalized envelope, we demonstrate accurate reconstruction of the underlying carrier waveform, even in the presence of carrier-steepening dynamics on a sub-50 attosecond timescale, and for propagation effects sensitive to the carrier-envelope offset phase. To our knowledge, these results represent the first explicit validation of the absence of a bandwidth constraint in nonlinear envelope equation modeling. Moreover, the GNEE facilitates ready comparison with familiar envelope-based formalism and is computationally very convenient, requiring only straightforward modification to existing generalized NLSE simulations. In addition, it allows the separate contributions of different linear and nonlinear effects to be readily identified, and we examine this explicitly for SC generation in tapered fused silica nanowires where we assess the particular importance of third-harmonic generation (THG) on the generated SC bandwidth.

## 2. Generalized Nonlinear Envelope Equation Model

The derivation of the GNEE uses standard techniques [5–17], but we review the main steps for completeness. Because our interest lies in pulse evolution in waveguides, we consider only one-dimensional propagation, but generalization to three dimensions is straightforward [6, 11]. Our starting point is the reduction of Maxwell’s equations to a one-dimensional second-order scalar wave equation. In the frequency domain, this is given by: [5]

Here *E*͂ (*z*,ω) =1/(2π) ∫^{∞}
_{-∞}
*E*(*z*,*t*)exp(*i*ω*t*)*dt*, where *E*(*z*,*t*) is the real electric field and *z* and *t* are laboratory-frame variables. Considering only a third-order nonlinear response, the operator 𝒩͂ has time-domain definition: 𝒩 = χ^{(3)}
*E*
^{2}(*z*,*t*)/*n*
^{2}(ω). Parameters β (ω) and *n*(ω) are the usual wavevector and refractive index.

Obtaining a computationally-efficient first-order envelope equation is now typically carried out using an envelope-carrier decomposition, followed by a factorization procedure subject to approximations such as the neglect of certain derivative terms in *z* or a perturbative treatment of the nonlinear response. It is these approximations that restrict the validity of the resulting equation to optical fields with temporal structure slower than an optical cycle [5–7]. However, it has been known for some time that the reduction of a second-order equation such as Eq. (1) to a first-order equation can actually be carried out without any assumptions on the underlying field timescale, provided that care is taken to correctly separate forward *E*͂_{+}(*z*,ω) and backward *E*͂ -(*z*,ω) propagating fields [14]. Indeed, analytical techniques based on Green’s functions [15], projection operators [16], and directional field variables [17] have all been demonstrated for this purpose. In our case, using the approach in Ref. [15] and defining the full field: *E*͂ (*z*,ω) = *E*͂_{+}(*z*,ω)+ *E*͂_{-}(*z*,ω), we can transform Eq. (1) into a pair of first-order equations:

At this stage we now decompose each component in terms of envelope and carrier e.g. *E*
_{+}(*z*,*t*) = [*A*
_{+}(*z*,*t*)exp(*i*ω_{0}
*t*)+c.c.)]/2, but it is crucial to note that the (generalized) envelope here can contain temporal structure on an arbitrarily fast scale. The field spectrum *E*͂_{+}(*z*,ω) thus retains frequency components at many multiples of ω_{0}, and it is this that allows the reconstructed field waveform *E*
_{+}(*z*,*t*) to accurately model sub-cycle dynamics.

Further simplification requires that nonlinear coupling between the forward and backward propagating components is neglected, but simulations using Maxwell’s equations have previously shown that this is, in fact, an excellent approximation for propagation at intensity levels below dielectric damage thresholds [13, 17]. From Eq. (2), we can then obtain a simple first order equation for each directional envelope component:

We proceed from here following the usual approach in nonlinear fiber optics where the forward propagating envelope [which we write now as *U*(*z*,*t*)] is normalized such that |*U*|^{2} yields the instantaneous power in watts [18]. After some algebra, we obtain the time domain GNEE:

Here, α is linear loss, and the β_{k}’s are the usual dispersion coefficients, although the global dispersion β(ω) can also be directly applied [18]. On the right-hand-side, the terms proportional to |*U*|^{2}
*U* and *U*
^{3} describe self phase modulation and THG respectively, and the nonlinear coefficient γ = ω_{0}
*n*
_{2}/*cA _{eff}*, where the nonlinear refractive index

*n*

_{2}and the effective area

*A*are evaluated at ω

_{eff}_{0}. In the absence of a frequency-dependent effective area, the envelope self-steepening timescale τ

_{ss}= 1/ω

_{0}, but this can be readily generalized for a dispersive

*A*[19]. The fractional Raman contribution to the nonlinearity is

_{eff}*f*, and the generalized Raman function is given by:

_{R}*g*(

*z*,

*t*,

*U*) = 2/3 (

*h*(

_{R}*t*) ∗ |

*U*(

*z*,

*t*)|

^{2})

*U*(

*z*,

*t*)+ 2/3exp(-

*i*ω

_{0}

*t*)(

*h*´

_{R}(

*t*) ∗

*U*

^{2}(

*z*,

*t*)´Re(exp(-iω

_{0}

*t*)

*U*(

*z*,

*t*)) with

*h*´

_{R}(

*t*) =

*h*

_{R}(

*t*)exp(2iω

_{0}

*t*). Note also that this Raman term includes rapidly oscillating components valid beyond slowly-varying envelope approximations [8]. Finally we note that we consider parameters of fused silica with

*n*

^{2}= 2.6×10

^{-20}m

^{2}W

^{-1}and, when Raman is included, we take

*f*= 0.18 and use the experimentally-measured Raman response function of silica [18].

_{R}## 3. Numerical Simulations

The GNEE can be readily solved numerically using the same modified split-step-Fourier scheme commonly used for generalized NLSE modelling, and transformation into a co-moving
5, 18]. The initial conditions on the envelope *U* are determined from the input field ℰ = [*U* exp(iω_{0}
*t*)+ c.c.]/2, and initial conditions on the carrier envelope offset phase can be specified through an additional constant phase term of the type exp(iΔϕ_{CEO}). Here we use normalization as above so that ℰ has units of W^{1/2}. Once the resulting output envelope is determined, the corresponding output field is then reconstructed and, with no restriction on the temporal structure that can be incorporated onto the envelope, this naturally allows for the inclusion of sub-cycle dynamics, if present.

To illustrate and validate the use of Eq. 4 for the modelling of sub-cycle dynamics, Fig. 1 shows results for one-dimensional propagation under conditions of optical shock generation, previously examined only through the direct solution of Maxwell’s equations [20–22]. Here we consider a 5 fs duration (FWHM) 830 nm hyperbolic secant pulse of intensity *I*
_{0} = *P*
_{0}/*A _{eff}* = 4.6×10

^{17}Wm

^{-2}propagating in the absence of dispersion and Raman scattering, and we assume zero initial carrier-envelope offset phase. The initial nonlinear strength is δ

_{n}=

*n*

_{2}

*I*

_{0}= (

*c*/ω

_{0})γ

*P*

_{0}= 1.2×10

^{-2}. GNEE simulations modeled propagation over 2.81 μm where theory predicts carrier shock onset [20–22]. Fig. 1(a) shows the output reconstructed temporal field, with the expected carrier steepening apparent in the exploded view in Fig. 1(b) (solid line). Physically, the shock arises from the generation of multiple odd-order harmonics that are seen in the corresponding spectrum (solid line) in Fig. 1(c). The generation of the oddorder harmonic cascade occurs as a result of the χ

^{(3)}nonlinearity mixing pump and harmonic fields through processes such as ω

_{0}+ω

_{0}+3ω

_{0}→ 5ω

_{0};ω

_{0}+ω

_{0}+5ω

_{0}→ 7ω

_{0}etc.

We have tested the validity of the GNEE simulations through a comparison with results from the direct numerical integration of Maxwell’s equations using the Pseudo Spectral Spatial Domain (PSSD) method [23]. The PSSD simulations are approximation-free and include both forward and backward propagating fields, and the results are shown as the open circles in Figs 1(b) and (c). Remarkable quantitative agreement with the GNEE results is obtained, even for this extreme case where spectral broadening extends over 10 times the carrier frequency, and the carrier steepening observed at the shock onset is on a sub-50 attosecond timescale. These results confirm that the GNEE approach can be used for accurate modeling of sub-cycle dynamics. Moreover, we have also verified that the GNEE model reproduces nonlinear dynamical effects that depend on the carrier envelope offset phase [13, 24, 25]. Specifically, the simulations above studying shock formation after 2.81 μm propagation were repeated with the introduction of an initial carrier envelope offset phase of Δϕ_{CEO} = π/2 on the input field. Although the effect of the initial phase on the dynamics is minor for this set of parameters, a difference is nonetheless manifested in the depth of the spectral minima, and Fig. 1(d) shows the spectral amplitude in the vicinity of the first minimum from the GNEE simulations for Δϕ_{CEO} = 0 (solid line), and for Δϕ_{CEO} = π/2 (dashed line). The PSSD simulation results for the latter case are shown as the circles, and the excellent quantitative agreement again confirms the GNEE validity.

Fig. 2 shows a further test of the GNEE, using parameters as above but also including fused silica dispersion through an appropriate Sellmeier expansion [3]. To focus on the modifications introduced by dispersion, we take τ_{ss} = 1/ω_{0} and neglect Raman scattering. GNEE results showing the output reconstructed temporal field after 100 μm propagation are shown in Fig. 2(a), with an exploded view near the pulse centre in Fig 2(b) (solid lines). Fig. 2(c) (solid line) shows the corresponding spectrum. The principle differences in the dynamics here arise from pump-THG walkoff, and cross-phase modulation on the THG field from the pump which displaces it from the ideal 3ω_{0} value. Despite these additional effects, however, the GNEE simulations are again in excellent agreement with the numerical integration of Maxwell’s equations (circles) shown in Figs. 2(b) and 2(c). We also note that propagation over the 100 μm distance here is convenient to illustrate the broadened pump and THG signatures clearly, but additional simulations confirm agreement between GNEE and Maxwell’s equation simulations over~mm distances when the spectral characteristics are more complex.

We expect that GNEE modeling will be particularly relevant to future SC generation experiments in highly nonlinear waveguides such as PCF nanowires. To this end, we have simulated propagation in a 600 nm diameter fused silica nanowire, including Raman effects and wavelength-dependent dispersion and effective area, but neglecting loss [2]. For pump pulses of 50 fs FWHM and 500 kW peak power at 1060 nm, Fig. 3 (solid line) shows the simulated spectra after 40 μm propagation. Even after this short distance, we see remarkable broadening, extending in fact to the UV transparency edge of fused silica near 180 nm. Significantly, comparison with simulations where THG is not included (dashed line) do not show this degree of broadening. This illustrates the important role that harmonic generation effects play in this regime, even under non-phasematched conditions. We also note that the influence of the novel nanowire dispersion on the propagating field [26] also plays a key role in extending the bandwidth beyond that which would be observed with standard fused silica dispersion (cf. Fig 2).

## 4. Conclusions

The major conclusion of this work is that the GNEE propagation equation presented here accurately models the generation of higher-harmonic frequency components and sub-cycle carrier shock structure arising from nonlinear propagation in a χ^{(3)} medium. This is explicitly confirmed by comparisons of GNEE and PSSD simulations. The GNEE possesses a number of particular features when compared to field-based propagation models, notably its expression in terms of the familiar NLSE-based formalism of nonlinear fiber optics, and its convenient numerical implementation through straightforward modifications to existing generalized NLSE solvers used for SC generation modeling. We anticipate that it will become widely used to study few cycle dynamical processes in nonlinear waveguides.

## References and links

**1. **J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys. **78**, 1135–1184 (2006). [CrossRef]

**2. **M. A. Foster, A. L. Gaeta, Q. Cao, and R. Trebino, “Soliton-effect compression of supercontinuum to few-cycle durations in photonic nanowires,” Opt. Express **13**, 6848–6855 (2005). [CrossRef] [PubMed]

**3. **V. P. Kalosha and J. Herrmann, “Self-phase modulation and compression of few-optical-cycle pulses,” Phys. Rev. A **62**, 011804(R) (2000). [CrossRef]

**4. **N. Karasawa, “Computer simulations of nonlinear propagation of an optical pulse using a finite-difference in the frequency-domain method,” IEEE J. Quantum Electron. **38**, 626–629 (2002). [CrossRef]

**5. **K. J. Blow and D. Wood, “Theoretical description of transient stimulated Raman scattering in optical fibers,” IEEE J. Quantum Electron. **25**, 2665–2673 (1989). [CrossRef]

**6. **T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. **78**, 3282–3285 (1997). [CrossRef]

**7. **N. Karasawa, S. Nakamura, N. Nakagawa, M. Shibata, R. Morita, H. Shigekawa, and M. Yamashita, “Comparison between theory and experiment of nonlinear propagation for a-few-cycle and ultrabroadband optical pulses in a fused-silica fiber,” IEEE J. Quantum Electron. **37**, 398–404 (2001). [CrossRef]

**8. **P. Kinsler and G.H.C. New, “ Wideband pulse propagation: single-field and multi-field approaches to Raman interactions,”Phys.Rev. A **72**, 033804 (2005). [CrossRef]

**9. **A. V. Husakou and J. Herrmann, “Supercontinuum generation of higher-order solitons by fission in photonic crystal fibers,” Phys. Rev. Lett. **87**, 203901 (2001). [CrossRef] [PubMed]

**10. **A. V. Husakou and J. Herrmann, “Supercontinuum generation, four-wave mixing, and fission of higher-order solitons in photonic-crystal fibers,” J. Opt. Soc. Am. B **19**, 2171–2182 (2002). [CrossRef]

**11. **M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E **70**, 036604 (2004). [CrossRef]

**12. **M. Kolesik, E. M. Wright, A. Becker, and J. V. Moloney, “Simulation of third-harmonic and supercontinuum generation for femtosecond pulses in air,” Appl. Phys. B **85**, 531–538 (2006). [CrossRef]

**13. **Y. Mizuta, M. Nagasawa, M. Ohtani, and M. Yamashita, “Nonlinear propagation analysis of few-optical-cycle pulses for subfemtosecond pulse compression and carrier envelope phase effect,” Phys. Rev. A **72**, 063802 (2005). [CrossRef]

**14. **Y. R. Shen, *Principles of Nonlinear Optics* (Wiley, New York, 1984).

**15. **A. Ferrando, M. Zacares, P. F. de Cordoba, D. Binosi, and A. Montero, “Forward-backward equations for nonlinear propagation in axially invariant optical systems,” Phys. Rev. E **71**, 016601 (2005). [CrossRef]

**16. **M. Kolesik, J. V. Moloney, and M. Mlejnek, “Unidirectional optical pulse propagation equation,” Phys. Rev. Lett. **89**, 283902 (2002). [CrossRef]

**17. **P. Kinsler, S. B. P. Radnor, and G. H. C. New, “Theory of directional pulse propagation,” Phys. Rev. A **72**, 063807 (2005). [CrossRef]

**18. **G. P. Agrawal, *Nonlinear Fiber Optics*, 4th ed. (Academic Press, San Diego, 2006).

**19. **B. Kibler, J. M. Dudley, and S. Coen, “Supercontinuum generation and nonlinear pulse propagation in photonic crystal fiber: influence of the frequency-dependent effective mode area,” Appl. Phys. B **81**, 337–342 (2005). [CrossRef]

**20. **G. Rosen, “Electromagnetic shocks and the self-annihilation of intense linearly polarized radiation in an ideal dielectric material,” Phys. Rev. A **139**, A539–A543 (1965).

**21. **R. G. Flesch, A. Pushkarev, and J. V. Moloney, “Carrier wave shocking of femtosecond optical pulses,” Phys. Rev. Lett. **76**, 2488–2491 (1996). [CrossRef] [PubMed]

**22. **L. Gilles, J. V. Moloney, and L. Vazquez, “Electromagnetic shocks on the optical cycle of ultrashort pulses in triple-resonance lorentz dielectric media with subfemtosecond nonlinear electronic debye relaxation,” Phys. Rev. E **60**, 1051–1059 (1999). [CrossRef]

**23. **J. C. A. Tyrrell, P. Kinsler, and G. H. C. New, “Pseudospectral spatial-domain: a new method for nonlinear pulse propagation in the few-cycle regime with arbitrary dispersion,” J. Mod. Opt. **52**, 973–986 (2005). [CrossRef]

**24. **P. Kinsler, G. H. C. New, and J.C.A. Tyrrell, ”Phase sensitivity of nonlinear interactions”, arXiv.org/physics/0611213.

**25. **P. Kinsler, S. B. P. Radnor, J. Tyrrell, and G. H. C. New, “Optical carrier wave shocking and the effect of dispersion,” Phys. Rev. E , submitted (2007). [CrossRef]

**26. **M. A. Foster, J.M. Dudley, B. Kibler, Q. Cao, D. Lee, R. Trebino, and A. L. Gaeta, “Nonlinear pulse propagation and supercontinuum generation in photonic nanowires: experiment and simulation,” Appl. Phys. B **81**, 363–367 (2005). [CrossRef]