Analysis of supercontinuum generation under general dispersion characteristics and beyond the slowly varying envelope approximation

The generation of broadband supercontinua (SC) in air-silica microstructured fibers results from a delicate balance of dispersion and nonlinearity. We analyze two models aimed at better understanding SC. In the first one, we characterize linear dispersion in the Fourier domain from the calculated group velocity dispersion (GVD) without using a Taylor approximation for the propagation constant. Results of our numerical simulations are in a good agreement with experiments. A novel relevant length scale, namely the length for shock formation is introduced and its role is discussed. The second part shows similar dynamics for a model that goes beyond the slowly varying approximation for optical pulse propagation.


1.
Introduction. Since the experimental demonstration of optical soliton propagation in single mode fibers some 20 plus years ago, the investigation of pulse dynamics in nonlinear optical fibers has evolved due to the introduction of novel structures with complex properties, such as photonic crystal and holey fibers [19]. In essence these are examples of engineered dielectric structures aimed at tailoring dispersive characteristics and enhancing nonlinear behavior. A direct outcome in terms of the pulse dynamics that has brought much attention from several experimental groups [11,14,13,10,7,29,26,22,17,24] is the ability to generate broadband supercontinuum spectra. Scientifically this is a departure from soliton dynamics that requires careful analytical and numerical modeling in parallel with the experiments. From the applications point of view, it has opened possibilities never seen before in areas such as frequency metrology [28] and medical diagnostics [16,25,23].
Supercontinuum generation (SCG) can arise from various physical processes such as self-and cross-phase modulation, and amplitude modulation [2]. Due to the complex interplay of linear and nonlinear phenomena in SCG dynamics, the theoretical formulation of the SCG mechanism imposes considerable challenges, in particular if this process happens in bulk media [21]. The major recent theory that explains the SCG for relatively low intensities in confined waveguides rests on the evolution and fission of higher-order solitons near the zero-dispersion wavelength in PCFs [17,18,12]. If the input wavelength is close to the zero-dispersion wavelength, then the influence of third-order dispersion is strong, thus a higher-order soliton with number N splits into its constituent solitons with the emission of blueshifted nonsolitonic radiation [30]. Since each soliton and its corresponding radiation has a different central frequency, the width of the generated total spectrum increases with increasing soliton number.
Recent experimental observations of supercontina in soft glass (schott SF6) PCF [22,24], however, suggest an interesting physical mechanism of SCG that cannot be fully explained by the previously known theories. In these experiments, SCG occurs in a dramatic fashion in the very early states of propagation, in particular at a length scale where solitons start forming. Such a phenomenon can only be explained if, initially, nonlinear effects other than soliton fission dominate the physics. Indeed, the underpinning mechanisms that generate supercontinua as reported in most theoretical and experimental studies are shock generation and its dispersive regularization in combination with multisoliton fission. The shock generation is a well known classical phenomenon in fluids and gas dynamics [31]. It also appears in ultra-short pulse propagation in fibers [1,3,32]. For ultrashort pulses, the refractive index depends on the pulse intensity, thus the center of the pulse envelope travels with a different speed than that of the trailing and leading edges of the pulse; this leads to an asymmetric shape of the pulse, which invokes shock formation. However, in optical propagation, dispersion plays an important role, preventing a sharp discontinuity. On the other hand, multisoliton generation resulting from small dispersion effects is a consequence of the integrability of the nonlinear Schrödinger equation (NLSE) [33]. Its eventual fission is the result of perturbations to the NLSE such as third order dispersion. Altogether, a universal feature of nonlinear dispersive wave phenomena is that the long term dynamics results from the delicate balance between linear and nonlinear effects.
To better illustrate this delicate balance we begin by studying a simpler model. Here we do not aim at studying a particular fiber; instead we highlight with a minimal number of linear and nonlinear terms different output scenarios. For this, consider the equation where A is the envelope of the pulse, c 3 , ǫ are respectively coefficients for the second and third linear dispersion terms, and c 2 is the coefficient of a self-steepening term. We can then characterize the dynamics of pulses under different regimes. In all instances we assume an input pulse A inp = √ 2/ cosh(t). (i) 0 < ǫ, c 2 , c 3 ≪ 1 for which N ≈  (ii) 0 < ǫ, c 3 ≪ 1 and both much smaller than c 2 , for which the dominant effect is shock formation at a propagation length proportional to 1 2c2 . The dispersion terms are viewed as perturbations which may or may not prevent shock formation. If they  In contrast, if linear dispersion is dominant and of order one, two cases which illustrate the overall dynamics are, first: (iii) c 3 = 1, c 2 ≪ ǫ ≪ 1 for which a soliton emerges and its dynamics under higher order corrections is well described by soliton perturbation theory [8] (Figure 3). In addition, third order dispersion accounts for linear wave resonances which explain the peak in the spectrum (figure 3b). Finally, (iv) ǫ = 0.1, c 3 ≪ c 2 ≪ 1 for which there is no soliton (Figure 4).
Notice the striking difference in the spectrum at the output. While a broadening of the spectrum is achieved in the first two scenarios, this is not the case for the last two, even if, as the temporal profile displayed in case (iv) suggests, the initial pulse is destroyed.
In the next section, we recognize the aformentioned outcomes for a more realistic model describing the pulse dynamics in photonic crystal fibers. The extended model accounts for all competing effects including self-steepening, which we believe is as important as the effect from the fully detailed linear dispersion. As we hoped the first part showed, in order to understand and exploit these phenomena, it is essential to obtain and analyze better these mathematical models. This in addition could explain for each instance in a real experiment what triggered SC generation. To begin with, the accurate broadband modeling of the dispersion relation is required to make sure one does not obtain spurious results, and to do so here we depart from the commonly used approach where a Taylor series expansion of the propagation constant β models the dispersive properties in a generalized nonlinear Schrödinger equation (gNLSE). Instead, we develop a mathematical model starting from calculated group velocity dispersion (GVD) curves. Then, we construct the function β(ω) over a broad frequency window and integrate the gNLSE preserving the spectral dependence of the propagation constant. As an illustration, we present our numerical results based on the calculated GVD for an LP 01 mode in an air-silica microstructured fiber studied by Dudley et al. [13]. Then, similar to what we did with equation (1), we carry out a careful numerical analysis. We find that if the nonlinear self-steepening term is strong enough, the model as it stands produces a shock that is not arrested by dispersion, whereas for weaker nonlinearity the pulse propagates the full extent of the fiber with the generation of a supercontinuum.
Recent studies, in particular for ultrashort pulse dynamics [9], have recognized that the slowly varying envelope approximation may not hold as a good model. In the last part of this paper we discuss how spectral broadening arises without invoking the slowly varying envelope approximation (SVEA).
2. Formulation of the model. The propagation of an electric field wave packet through an optical fiber can be described by gNLSE [1], Here, the variables z, T, ω represent propagation distance, time and optical frequency, respectively. The envelope of the wave packet is A, and c, λ, ω 0 , β(ω) represent the velocity of light in vacuum, wavelength, central frequency and wave number, respectively. F −1 denotes the inverse Fourier transform, andÂ is the Fourier transform of the pulse envelope. Finally the self-steepening term models the instantaneous nonlinear response function of the medium, which is a good approximation given the temporal lengths of the pulses. The inclusion of a Raman  (non-instantaneous) term we believe will only introduce a shift in the peak location in the spectrum (see Fig. 6 and compare with Fig. 5b in [13]).
The effects of fiber dispersion are accounted for by the propagation constant β(ω) which we calculate based on the dispersion profile presented in [13], without  performing a Taylor series expansion around the carrier frequency. Using two highprecision numerical integrations of an accurate rational interpolant of the GVD curve, we obtain the GVD function D(s). Then, the group velocity ν g (ω) is derived By setting F (λ) = λ λ0 D(s)ds, we obtain Since ∂β ∂ω = 1 νg , it follows that ∂β ∂ω By integrating Eq. (5) with respect to ω and using the relation λ = 2πc ω , we obtain We employ a frame of reference moving with the pulse at the group velocity ν g by making the transformation t = T − z/ν g . In the end, we obtain The resulting equation preserves the complete structure of fiber dispersion which is indeed utilized in experiments. In addition, the equation is valid not only for broad pulses, but also short pulses since the derivation is carried out without the assumption of a pulse centered around a specific carrier frequency (without Taylor series expansion of the propagation constant around a carrier frequency). In the remainder of this section, we present analytical and numerical results obtained from the gNLS (7).
2.1. Optical shock formation. In order to first pay attention to the nonlinear effects governing the mechanism of shock formation [1,3], we consider the dispersionless case by setting F (λ) = 0 in Eq. (7). In the absence of dispersion, we first split Eq. (7) into an intensity-phase system by adding and subtracting where A * is the complex conjugate of A. By defining I = |A| 2 , the addition of Eqs. (8), (9) gives The general solution of Eq. (10) is where f (t) is determined by the initial pulse shape, namely, f (t) = I(0, t). The solution form Eq. (11) implies that asymmetric distortion of the pulse will occur eventually. From Eq. (11), we also find The resulting equation shows that after a distance z s = − ω0 3γ 1 f ′ , a singularity in the pulse intensity will be generated, namely, the formation of an optical shock. This shock does play an important role in the spectral broadening once dispersion regularizes it. In other words, the effects of fiber dispersion cannot be ignored. Moreover, the effect of GVD becomes more important as the pulse steepening becomes significant. This phenomenon prevents further steepening of pulse shape, i.e.,an appropriate strength of linear dispersion results in a mechanism that may prevent (or regularize) the shock.

2.2.
Numerical solutions of the generalized nonlinear Schrödinger equation. Using the GVD profile employed in experiments [13], we perform our numerical simulations of pulse dynamics based on Eq. (7). In particular, we consider the propagation of 100f s pulses at 780nm in a 1-m length air-silica microstructured fiber with γ = 0.1W −1 m −1 . We assume that the input pulse has a form of A 2 0 / cosh 2 t t0 . As a reference, we find from the previous analysis that the shock length for this input pulse is z s = 22.1cm which is much shorter than the actual fiber length, thus in this case shock regularization is a likely scenario for SC generation.
The integration spectrally covers the 400 − 900nm range and, as stated above, it does not require a Taylor expansion of β. Instead we calculate β via two highprecision numerical integrations of the GVD. The numerical evolution is then completed using a standard adaptive ode solver. The results displayed below (figures 5, 6) are in clear qualitative agreement with those in [13]. We should point out two important distinctions between figure 5 in [13] and figure 5 here: we do not capture the peak in the spectrum at wavelengths close to 1200nm in 5a,b of [13] and the corresponding pulse (labeled C) shown in figure 5c. This is because in our approach we computed the dispersion profile based on the calculated GVD curve shown in figure 2 of [13]. This calculation did not extend to wavelengths beyond 900nm and we did not extrapolate such curves. This explains the sharp decay in the spectrum of figure 5 here. Calculations based on a Taylor expansion of β which is commonly used, can in principle be extended to any desired spectral range. On the other hand, our result better reproduces the observed supercontinuum spectrum (figure 5a of [13]) in the short (less than 600nm) wavelength portion and is as good as the Taylor expansion in the intermediate regime.
Finally, figure 5 (right) shows five distinguishable pulses at the output. In figure  6, we spectrally isolate each pulse and find their peaks centered approximately at: 867nm (peak 1), 913nm (peak 2), 848nm (peak 3), 840nm (peak 4) 858nm (peak 5). As we state below, the spectral separation should be accentuated by the presence of the Raman shift which we did not include in the model. What is most important here is we corroborate spectral broadening and that splitting of pulses occurs, with the spectral shift accounting for differences in soliton velocities.
In trying to understand the critical balance between linear and nonlinear effects, we now depart from the concrete example to further illustrate this interplay in a series of simulations of the equation shown below, which is no different than equation (7) except that we placed two adjustable constants c 2 (c 3 ) in front of the self-steepening (linear) term to account for their respective strengths.
By proper re-scaling of the propagation variable z and the pulse peak amplitude, one can eliminate the parameter c 3 but for clarity we analyze our simulations in the 2-parameter space while using the same initial condition. Observe that from (12), an increase of c 2 effectively means the shock length is reduced. In practice, this shock length reduction can be induced by having shorter input pulses. By allowing ourselves to modulate dispersion and nonlinearity through these two parameters, we hope to highlight how delicate this balance is. Figure 7 which presents a separation between two distinct outcomes was obtained by careful simulations in the (c 2 , c 3 ) parameter space. For the region above the curve the shock is not arrested and below the curve dispersion regularizes the shock. It is important to point out the results shown below will strongly depend on the dispersion profile. Nonetheless it is intriguing to see from figure 7 that a universal critical value ( c2 c3 ) C ≈ 2.1 emerges. While one could argue that by modifying c 3 one departs from a particular photonic structure, what matters is that for every value c 3 (that is moving vertically in figure  7), this transition always occurs. At this time, we do not have an explanation for it. Furthermore, this property should be tested for different dispersion profiles. To summarize, by performing a series of careful numerical simulations where we look at the relative strengths of the dispersion (measured by a parameter c 3 that multiplies F (λ) in Eq. (1)) and of the self-steepening term (measured by a parameter c 2 ) we clearly demonstrate two dynamical regions: one where the singularity due to the shock is not suppressed by dispersion (the region above the curve). In this regime the spectral broadening does not saturate and the numerical solution blows up, clearly suggesting that additional physical mechanisms must be considered. In the second region (the region below the curve) propagation leading to supercontinuum generation or other dynamics similar to those highlighted in the first part of the paper is observed. Although we did not show the curve beyond c 3 = 0.8, it should be clear that the point c 2 = c 3 = 1 corresponding to the experimental parameters in [3], is as expected below the curve.
3. SC generation beyond the SVEA. As research into SC generation extends into regimes exhibiting the formation of ultra-short temporal pulses, the validity of the slowly varying envelope approximation (SVEA), even if we incorporate a detailed dispersion profile over a broad spectral range, requires careful consideration. Deriving unidirectional propagation equations for short pulses that would start from Maxwell's equations and account for dynamics such as higher harmonic generation and sub-cycle shock structure leads to complex envelope equations which can only be solved numerically [20,15]. One can instead build computational schemes that solve Maxwell's equations for ultra-short pulses [27]. Careful integration of generalized envelope equations can relove pulses as shot as sub-50 attoseconds [15]. Alternatively, in this section, we illustrate the universality of the features discussed above and observed in many experiments; optical shock formation, dispersion regularization and spectral broadening arise, employing models that do not use the SVEA. Consider where χ (1) , χ (3) are linear and nonlinear susceptibilities (for a detailed derivation of this equation, see [9]). Due to the characteristics of its derivation, this equation again preserves the nonlocal nature of the pulse, and thus is valid for both broad and ultrashort pulses. Moreover, it was reported in [9] that A(t) = α 1+(βt) 2 is a stationary solution provided that χ (3) = 1 and the Fourier transform of linear susceptibility takes the parabolic form, for some real constants α and β. While this may not be a realistic dispersion profile, we believe it illustrates well phenomena observed in realistic scenarios.
Here we assume thatχ (1) has the form of Eq. (15) with α = 0.2, β = 0.75, and an initial pulse of 0.4 1+0.75 2 t 2 . Note that the amplitude of the initial pulse is larger than α, and thus the solution will not be stationary and we expect nonlinear effects will be dominant.
We first illustrate the role of nonlinear effects by turning off the linear susceptibility, i.e., χ (1) = 0, and setting the nonlinear coefficient χ (3) = 1. Numerical simulations of pulse propagation based on Eq. (14) show the output pulse at propagation length z = 5 and its spectrum (Figure 8). In this case, the nonlinear effects lead to steepening of the pulse shape, which results in the broad spectrum of the output pulse. In the other extreme, Fig. 9 illustrates the pure linear effects by turning off the nonlinearity, i.e., χ (3) = 0, and using the linear susceptibility as in (15). At propagation length z = 5, the pulse has slightly changed its form and also its spectrum. Finally, Fig. 10 presents the combined effects of both nonlinearity and linearity, i.e., χ (3) = 1 and χ (1) as in (15). Compared to Fig. (8), the shock is suppressed by linear dispersion. However, with our choice of input conditions, the nonlinear effect still dominates the dynamics, which induces spectral broadening at finite propagation distances.

4.
Conclusions. Supercontinuum generation is a fascinating and important phenomenon observed in certain nonlinear wave systems. In this work, we first discussed several simple models where we tuned dispersion and nonlinearity so that we could showcase different outcomes. In particular we showed shock driven SC generation as well as soliton-fission driven SC. Next we moved to a model closely related to an existing photonic crystal fiber and showed both SC generation as well as critical shock formation. It is important to emphasize that by properly integrating the dispersive terms for a given photonic microstructured fiber, we capture supercontinuum generation as observed in experiments, likely to greater accuracy than the more common expansion to a finite order of the linear dispersion relation. Our numerical simulations illustrate that for some input conditions, shocks rather than soliton fission appear to be dominant and become the major source of spectral broadening. It is true that soliton fission as seen in many works could be the leading mechanism towards SC generation. Which effect is more dominant and what signatures (if any) of the spectral picture can explain the hierarchy of effects coming into the dynamics remains unclear. Interestingly, two recent theoretical papers using a wave-turbulence approach [4,6] propose an explanation of SC generation. In [4], the authors claim that coherent structures (solitons, shocks) no longer play a significant role and instead SC results from a nonequilibrium thermodynamic process. This was based in an NLSE-type model which included higher order dispersion, but excluded self-steepening. On the other hand, if self-steepening is included [6], it plays a critical role in thermalization towards a two peak SC spectral profile.An experimental signature of such optical thermalization was recently reported in [5] Overall, as these works clearly illustrate, an accurate mathematical model is essential to better reflect the experimental outcomes. In particular we have a numerical approach at our disposal to study any photonic fiber structure for which GVD profiles have been or can be computed. Any time one models SC with a slowly varying envelope approximation the question of validity comes to mind. In the last section we presented a simple 1D nonlinear Maxwell's equation model where spectral broadening is achieved, thus re-enforcing the overall principle that SC generation is universal. While more realistic models that go beyond the SVEA are more difficult to analyze, careful numerical simulations will lead to a better understanding of this very rich phenomenon.