1. Introduction

Stimulated Raman scattering is a central process of nonlinear fiber optics, and has been extensively investigated since the early days of research in this field [1–3]. Raman scattering introduces a delayed contribution to the nonlinear response of the medium and, for ultrashort pulse propagation, leads to well-known effects such as the Raman soliton self-frequency shift and Raman induced soliton fission [4–7]. From a theoretical perspective, the effect of Raman scattering on pulse propagation in fibers can sometimes be described using coupled amplitude equations [8], but a more general approach is to incorporate the Raman response directly within a generalized nonlinear Schrödinger equation (GNLSE) model [6,9,10]. This is particularly the case for broadband supercontinuum (SC) generation where Raman-induced effects strongly influence the spectral and temporal evolution and where accurate modeling of dynamics is essential for applications such as telecommunications [11,12].

Many GNLSE-based studies of Raman scattering model the Raman response using experimental gain measurements or accurate lorentzian approximations [13–15]. On the other hand, a convenient first-order linear approximation to the Raman response has also often been used, and applied to both theoretical and numerical investigations [16–28]. The analytic simplification of a linear approximation has led to important insights into many different propagation effects: the soliton self-frequency shift [5]; the formation of quasi-solitons that couple a NLSE soliton to a Raman tail described by an Airy function [16]; soliton-dispersive wave interactions [17]; novel classes of nonlinear localization [18]; and optical shock phenomena [19]. The linear Raman model has also been used in numerical simulations to provide additional important insights into nonlinear spectral broadening in fibers; this work includes early studies of the basic dynamics of spectral broadening in fibers [20,21] as well as more recent studies of the detailed mechanisms of supercontinuum generation [22–28].

The main objective of this paper, however, is to stress that whilst a linear gain approximation may possess significant advantage in analytical studies, its use in numerical simulations of SC generation must be considered with great care. Certainly there are parameter regimes where the linear approximation will capture the essential physics, but this is by no means guaranteed and the validity of the linear model in simulations needs to be examined very carefully on a case by case basis.

Here, we examine these issues in detail by comparing simulations of SC generation that assume a linear Raman response with simulations using the full experimental Raman response curve. We consider differences in the spectral, temporal and stability characteristics for excitation using pulses in the 50 fs–4 ps range, and we identify regimes where the linear Raman approximation predicts the output SC characteristics to good accuracy, and regimes where it does not. An important conclusion is that for conditions typical of many experiments, although the input pulse may satisfy the criteria appropriate for the use of the linear Raman gain approximation, the breakup of the input pulse through soliton fission or modulation instability can rapidly lead to the appearance of ejected solitons whose temporal width is so short that the subsequent propagation dynamics are inaccurately described by the linear Raman approximation. As Raman effects have recently been shown to be important in influencing the statistics of highly incoherent SC generation [29–31], we also compare the spectral characteristics predicted by the full experimental Raman response and the linear Raman response in this regime. Finally we comment on a particular form of numerical artifact associated with pulse collapse that we have identified when using the linear model.

2. Numerical modeling

We begin by reviewing the relevant nonlinear propagation equations. We use dimensional forms of the propagation equations following the normalization that is usual in nonlinear fiber optics [8]. Our numerical modeling is based on the GNLSE written as follows:

We now consider the two approaches used to include the Raman response within the GNLSE given by Eq. (1). The most realistic approach uses the experimentally-measured Raman cross-section to determine the frequency domain form of the finite-bandwidth response function given by ${\tilde{h}}_R(\omega )$ = FT[h_R(t)]. Here FT denotes Fourier transform. For fused silica, this experimental data (adapted from Ref. 13) is shown as the solid line in Fig. 1(a) . This curve exhibits the well known gain peak around ω/2π = 13.2 THz, with the response becoming negligible at frequencies greater than 30 THz. In what follows we refer to this for convenience as “the full model”. Note that in our simulations, we use the experimental form of the response curve taken from Ref. 13, but nearly identical results are obtained using the multiple vibrational mode model of Ref. 15.

 

Fig. 1 (a) Comparison of the two forms of the Raman response function: the experimental response (solid line) is adapted from Ref. 13. The linear response (dashed line) is calculated as explained in the text. (b) As a function of initial pulse bandwidth, we compare the rate of Raman self-frequency shift experienced by an N = 1 soliton for: the experimental (full) model (black) and the linear model (red).

Download Full Size | PDF

The linear approximation to the frequency-domain Raman response function can be written as ${\tilde{h}}_R(\omega )$ = 1 + iωΒ, where B is a slope parameter used to fit the measured Raman gain. Subject to this approximation and also assuming that the effect of shock term on the Raman response is negligible [32], the GNLSE then takes a simpler form [8]:

The characteristic Raman time scale T_R is derived from the Raman gain slope using T_R = f_R B, and can be determined from the slope at ω = 0 such that T_R assumes the interpretation as the first moment of h_R(t) [8]. However, the linear approximation can be used to fit the experimental curve in other ways [33–35], and the dashed line in Fig. 1(a) plots the linear Raman response fitted over the region between the pump and the gain peak. In this case T_R = 3.0 fs and this is the value used in all subsequent simulations. In what follows we refer to this approach described by Eq. (2) as “the linear model”.

3. Soliton propagation and coherent supercontinuum generation

It is first instructive to consider the limitation of the linear model when describing the basic dynamics of the soliton self-frequency shift for solitons of different bandwidth injected into a fiber. To this end we consider the propagation of a fundamental N = 1 soliton in a fiber with only a constant dispersion coefficient β ₂ = −4.10 ´ 10⁻¹ ps² km⁻¹ and nonlinearity γ = 0.01066 W⁻¹m⁻¹, neglecting higher order dispersion and shock terms. These dispersion and nonlinearity parameters are typical of highly nonlinear photonic crystal fibers with zero dispersion wavelength around 1 μm (see below for further details of the fiber considered). For this case, Fig. 1(b) plots the rate of self frequency shift |dν/dz| as a function of the initial pulse spectral bandwidth. Note that the bandwidth of the soliton actually evolves negligibly in the absence of higher order dispersion, remaining constant with propagation distance. From this figure, we see how the rate of soliton self-frequency shift predicted by the linear model with T_R = 3 fs is in good agreement with that obtained from the full model only when the pulse bandwidth is less that ~5 THz (corresponding to a pulse duration FWHM τ > 60 fs). As the bandwidth increases above 5 THz, it is clear that the deviation between the two models start increasing, with the linear model significantly overestimating the rate of self-frequency shift. As we shall see in what follows, it is this fundamental difference in the way in which the two models predict the soliton dynamics that imposes the main restriction on the use of the linear Raman model for SC simulations.

We are now in a position to compare the use of the two models to simulate the more complex process of SC generation. We focus on propagation dynamics near the fiber zero dispersion wavelength where Raman effects play a central role in both the initial stages of pulse break up and the subsequent extension of the spectrum to longer wavelengths through the soliton self-frequency shift [11]. We first consider simulations of SC generation under conditions of coherent spectral broadening [36]. We use parameters typical of experimental conditions, modeling propagation in a commercially-available photonic crystal fiber (SC-5.0-1040) with zero dispersion wavelength at 1030 nm. We consider hyperbolic secant input pulses at 1035 nm, and noise is included through a one-photon-per mode background. Note that accurate comparison with experiments under more general conditions may require more careful consideration of the noise source [37]. At a wavelength of 1035 nm, the fiber nonlinear and shock parameters are γ = 0.01066 W⁻¹m⁻¹ and τ_shock = 0.55 fs. The fiber dispersion curve was included directly in the numerical simulations in the frequency domain, but near-identical results can be obtained using dispersion coefficients (at 1035 nm): β ₂ = −7.10 ´ 10⁻¹ ps² km⁻¹, β ₃ = 6.90 ´ 10⁻² ps³ km⁻¹, β ₄ = −1.02 ´ 10⁻⁴ ps⁴ km⁻¹, β ₅ = 2.7 ´ 10⁻⁷ ps⁵ km⁻¹, β ₆ = −8.3 ´ 10⁻¹⁰ ps⁶ km⁻¹, β ₇ = 3.0 ´ 10⁻¹² ps⁷ km⁻¹, β ₈ = −1.2´10⁻¹⁴ ps⁸ km⁻¹, β ₉ = 6.5′10⁻¹⁷ ps⁹ km⁻¹, β ₁₀ = −3.2´ 10⁻²⁰ ps¹⁰ km⁻¹.

We begin our discussion of the validity of the linear Raman model by showing in Fig. 2 selected results of simulations comparing SC spectra obtained using the linear model (red) with those obtained using the full model (black). We consider differences between the two models for input pulses with durations (FWHM) spanning over the range τ = 50 fs – 1 ps. As the pulse duration is varied, we adjust peak power to keep constant soliton number N = (γP ₀ T ₀ ²/|β ₂|)^1/2 with T ₀ = τ /1.763. For a meaningful comparison between results using different input pulses, we compare spectral characteristics after a propagation distance of z = 10 L _fiss where the characteristic soliton fission distance is L _fiss = L _D/N with L _D = T ₀ ²/|β ₂| [11].

 

Fig. 2 Comparison of output spectra for different input pulse parameters as shown, using the realistic finite bandwidth Raman response (black) and the linear Raman model (red). The propagation distance in each case is z = 10 L _fiss. Dimensional parameters are: (i) P ₀ = 5.2 W, z = 900 m; (ii) P ₀ = 129.4 W, z = 36 m; (iii) P ₀ = 2.07 kW, z = 2.26 m; (iv) P ₀ = 20.7 W, z = 450 m; (v) P ₀ = 517.5 W, z = 18 m; (vi) P ₀ = 8.28 kW, z = 1.13 m. Note the change in wavelength span in figures (iv)-(vi).

Download Full Size | PDF

It is clear from inspection of Fig. 2 that in all cases the linear model predicts a SC bandwidth greater than that predicted by the full model, and the discrepancy between the models increases for larger input soliton number and/or shorter pulse duration. Nonetheless, we note that for the results obtained using N = 5 input pulses, there is good qualitative agreement between the spectra predicted by both models, and for the 1 ps case in particular, there is also very good quantitative agreement. On the other hand, as the soliton number increases to N = 10, the degree of quantitative agreement decreases for the 1 ps case, and for the shorter pulse durations of 200 fs and 50 fs we see significant qualitative differences between the spectra predicted by the two models.

These observations can be explained by the fact that, in contrast to the full Raman model, the linear model does not impose any gain bandwidth limitation, leading to significant differences in the self-frequency shift dynamics of ejected solitons from the SC. Specifically, as we have seen above in Fig. 1(b), the linear model overestimates the rate at which solitons shift toward low frequencies. In fact, we can use the results in Fig. 1(b) to obtain quantitative insight into the differences between the SC simulations seen in Fig. 2 using the linear and full models. In particular, careful examination of the results in Fig. 2 (as well as others carried out over a wider parameter range) show that differences between SC spectra computed using the linear and full Raman models occur for cases when the first soliton that is ejected from the input pulse though soliton fission has a bandwidth approaching 5 THz, corresponding to a temporal FWHM τ ~60 fs. Moreover, since the duration of the first ejected soliton can be predicted from the FWHM of the input pulse and soliton number N using τ ₁ ~τ / (2N −1) fs [6], we obtain a convenient quantitative condition describing the regime of validity of the linear model as: τ / (2N −1) > 60 fs. Thus, for large input pulse durations and/or small input soliton numbers SC generation may be modeled accurately by a linear Raman gain model. In this regard, we note that the soliton number-dependence of this validity condition allows the potential accuracy or inaccuracy of the linear approximation to be conveniently tested for different combinations of peak power and/or duration of the input pulses.

It is instructive now to consider in more detail the case where there is significant disagreement between results from the linear and full Raman models. To this end, Fig. 3 compares the simulations using both models for the case of a 50 fs, N = 10 input pulse [as in Fig. 2 (vi)], and we follow the dynamical evolution of the spectrum over a distance of 20L _fiss. The figure shows also corresponding output temporal profiles and temporal spectrograms. There are clear physical differences between the field characteristics predicted using the different models: besides the significantly larger bandwidth predicted by the linear model, the most notable difference is that whilst the full Raman gain model predicts significant spectral content on both sides of the ZDW for the output spectrum, the corresponding spectrum predicted by the linear Raman model lies almost entirely in the anomalous dispersion regime.

 

Fig. 3 Simulated propagation dynamics and output characteristics for a 50 fs input pulse at 1035 nm with N = 10 using (a) the linear Raman model and (b) the full Raman response. For propagation distance of z = 20 L _fiss subfigures (a) and (b) show the output temporal profile (top) and the spectral evolution dynamics (bottom). Results in (c) show the output spectrograms for both cases as indicated. Spectrogram movies of the results in Fig. 3(a) and 3(b) are at (Media 1) and (Media 2) respectively. The spectrograms are calculated using a 200 fs temporal gate function to sample the evolving field (see e.g. Ref. 11 for more details on the spectrogram representation of supercontinuum generation.)

Download Full Size | PDF

When we inspect the spectral evolution in the framework of the linear model as shown in Fig. 3(a), we see that whatever spectral energy that is initially present in the normal dispersion regime is in fact depleted with propagation. This can also be seen more clearly by inspecting the multimedia movie files linked to in the caption of Fig. 3 that follow the spectrogram evolution with distance for the two cases. We attribute this behavior to two effects: (i) the cross-phase modulation mediated interaction between the shifting soliton and dispersive wave [15,38] that results in temporal overlap of the components over the full propagation distance, and (ii) the fact that the absence of bandwidth limitation in the linear model actually allows the normal dispersion dispersive wave component to act as a Raman pump for temporally-overlapping anomalous dispersion regime components. This results in the continuous loss of energy of the dispersive wave component. In this context note that although the dispersive wave is of low amplitude, the Raman gain increases linearly and very significantly over the ~170 THz frequency interval between the edges of the spectrum. Moreover, the energy transfer from the dispersive waves to the solitons increases the rate of soliton self-frequency shift as well. Finally, we remark here on the very distinctive form of the time domain profile in Fig. 3(a) that is predicted by the linear model. As seen in the corresponding spectrogram this intensity profile consists of two distinct soliton peaks (labeled A and B), each peak associated with an extended lower amplitude pedestal on its leading edge. In fact, this structure is a well-known analytic solution to the NLSE with a linear Raman gain and anomalous dispersion (i.e. Eq. (2) with τ _shock = 0, β ₃ = 0, β ₂ < 0 and T _R ≠ 0) and consists of a quasi-soliton state of a hyperbolic secant pulse coupled to an extended pedestal described by an Airy function [16]. Such a pronounced Airy pedestal, however, is not observed in the time domain profile predicted by the full model. Note that simulations carried out over a much wider parameter range suggest that the condition of an initial ejected soliton bandwidth > 5 THz provides an extremely good indicator of when similar differences between the linear and Raman models will be observed in SC simulations. We further note that the differences in the observed dynamics also impact strongly on the initial soliton fission process itself under conditions when it is primarily induced by the Raman effect (pumping far from the zero dispersion wavelength for example). In particular, the fact that the linear model admits quasi-soliton solutions between a soliton and an extended Airy tail means that simulations using the linear model show a reduced number of ejected solitons at a given distance when compared to the full model.

4. Incoherent supercontinuum generation

The results above clearly show different dynamics using the linear Raman and full Raman gain models. And we have carefully checked that with the parameters used above, the effect of noise was negligible using both models. On the other hand, the shot-to-shot stability characteristics of SC generation have recently been receiving a great deal of attention in the context of analogies between optical field fluctuations and the emergence of giant oceanic rogue waves [27–29]. In this section, we consider how noise manifests itself in the incoherent SC regime using linear and full Raman models, focusing specifically on differences in spectral characteristics between the two cases.

Our simulations here use the same fiber parameters as above, but consider a 4 ps pulse at 1035 nm of 150 W peak power. In this case, we are in a high soliton number regime such that N ~107, and propagation would be expected to be highly incoherent with significant shot-to-shot spectral and temporal fluctuations. Figure 4 compares results from simulations using the linear Raman gain and full Raman gain models. The two models predict very different results, with the absence of gain bandwidth limitation in the linear model leading to significantly greater spectral broadening at the same propagation distance as shown in Fig. 4(a) for a distance of 15 m. In Fig. 4(b) we plot the outputs at two different distances such that there is comparable mean spectral broadening on the long wavelength side of ~200 nm at the −20 dB level relative to the residual pump. This level of broadening is reached after only 11 m of propagation with the linear model compared to the much greater propagation distance of 25 m when using the full model.

 

Fig. 4 Simulated propagation dynamics and output characteristics for a 4 ps pulse at 1035 nm of 150 W peak power such that N ~107 and propagation is highly incoherent with significant shot to shot fluctuations. Results from simulations using the linear Raman gain and full Raman gain models are compared (a) at fixed propagation distance of 15 m and (b) for comparable mean spectral broadening on the long wavelength side of ~200 nm at the −20 dB level relative to the residual pump.

Download Full Size | PDF

The figure shows several important features. Firstly, we note that both models exhibit significant shot to shot fluctuations and, indeed, we have confirmed that both models yield near zero degree of spectral coherence across the generated spectra [11]. This is an important observation as it confirms that SC fluctuations can be observed irrespective of the particular Raman model used. On the other hand, the very significant differences in SC bandwidth at fixed propagation distance suggests that for a quantitative agreement with experiment in this regime, the full Raman model must be used because the linear Raman model greatly overestimates the spectral broadening. This is particularly the case when comparing experiments and stochastic simulations when spectral filtering is used to quantify the statistical fluctuations at particular wavelengths in the incoherent SC [29–31].

From a dynamical perspective, it is important to note that the differences between the two models in the long pulse regime at these higher soliton numbers again occurs because of the dramatic effect that Raman scattering has on the emergence of solitons from the initial phase of modulation instability (MI) [39]. In fact, although deterministic soliton fission does not occur in this long pulse regime, the temporal subpulses that appear on the modulated pulse envelope still possess characteristic bandwidths exceeding 5 THz such that significant differences occur between the dynamics predicted between the two models as these subpulses evolve into self-frequency shifting solitons. In addition, in the multi-soliton regime considered with these parameters, the unbounded linear Raman gain can also facilitate significant energy transfer between temporally overlapping spectral components. This results in a smaller number of emerging solitons of increased energy which undergo an enhanced rate of frequency-shift leading to the much increased spectral bandwidth as seen in Fig. 4. Such unphysical energy transfer is typically observed during soliton-soliton, soliton-residual pump and/or soliton-dispersive wave interactions.

The results above clearly show that for the case of long pulses in the picosecond regime, both the linear and full models yield similar results in predicting incoherent supercontinuum spectra. However, we stress again that the dynamics, spectral structure and bandwidths predicted by the two models are very different as is apparent from Fig. 4. These results also raise the natural question of how the two models compare in terms of predicting coherence for various cases. To investigate this, we performed additional simulations for shorter pulses of 100 fs and 200 fs duration, but for varying input soliton numbers in the range 5-40. For each case, the average degree of supercontinuum coherence <|g⁽¹⁾ ₁₂|> [31] was calculated explicitly from an ensemble of 200 simulations using both linear and full models, and the results are shown in Fig. 5 . Here we plot the mean coherence evaluated across the broadened pulse spectrum. We see that for low soliton number N < 10 when noise-seeded MI plays only a minor role in the initial propagation dynamics [11], both models reproduce the expected coherent spectral broadening (<|g⁽¹⁾ ₁₂|> ~1). On the other hand, for higher soliton numbers when noise plays an increasingly important role in the dynamics, we see that the linear model predicts increased decoherence (lower <|g⁽¹⁾ ₁₂|>) when compared to the full model. This behavior can be readily understood because the coupling between the Raman and MI gain leads to an essentially unlimited bandwidth for noise amplification in the linear model that greatly overestimates the degree of decoherence compared to the full model.

 

Fig. 5 Supercontinuum coherence as a function of soliton number for (a) 100 fs and (b) 200 fs pulses. Results from the full model (open circles) and linear model (closed red circles) are compared. The predicted coherence properties from the two models differ significantly in the partially coherent regime.

Download Full Size | PDF

5. Pulse collapse in the linear Raman gain model

Finally in this section, we describe a particular form of numerical artifact that we have found can lead to a dramatic pulse collapse when using the linear Raman model to simulate SC generation. In fact, although we refer to this collapse process as an artifact because inducing it in simulations involves a particularly poor choice of numerical gridding parameters, it may well be possible that certain physical conditions may also yield the same type of behavior leading to dynamical instability and intermittence phenomena [40].

The collapse process is associated with the dramatic loss of energy of solitons ejected from the input pulse through soliton fission, and can be conveniently illustrated by considering the case of coherent SC generation in Fig. 3(a), but using simulations at lower computational time resolution. Figure 6 shows the results of simulations for the physical parameters of Fig. 3(a) but with different numerical resolution at the same timespan: the resolution in Fig. 3(a) was 1.46 fs/point but for Fig. 6 it is reduced to 2.93 fs/point. It is important to note that although the spectral window is reduced with the decrease in resolution in Fig. 6, it still spans from ~650 nm to ~2500 nm, a range large enough to accommodate the expected SC spectrum (see Fig. 3). This is an important point because it allows us to attribute the observed artifact to inadequate numerical time-resolution and to confirm that frequency domain boundary effects do not play a role.

 

Fig. 6 Simulated propagation dynamics of soliton collapse observed in the linear model. (a) Soliton energy vs. distance (b) temporal profile evolution of the soltion plotted in a frame of reference moving at the soliton velocity. (c) and (d) show the output spectrograms before and after collapse at z = 16 L _fiss and z = 20 L _fiss, respectively. The spectrograms were calculated using a 200 fs gate function. Input pulse parameters are identical to that of Fig. 3. A movie of the spectrogram evolution in (c) which shows also the projected temporal and spectral evolution plotted on log scales is at (Media 3).

Download Full Size | PDF

Figure 6 (a) illustrates the evolution of the energy of the first soliton ejected from the input pulse through soliton fission. The energy remains nearly constant as expected up to a distance of approximately 17L _fiss. However, upon subsequent evolution the soliton undergoes an abrupt loss of energy which is associated with a sudden collapse of the soliton in the temporal domain [Fig. 6(b)]. Interestingly, inspection of the spectrogram representation [Fig. 6 (c)] reveals that during the collapse the energy lost by the soliton is in fact transferred into the normal dispersion regime into the vicinity of the trapped dispersive wave. A link to a movie of the evolution in Fig. 6 is provided in the caption.

We associate the collapse with the way in which the bandwidth-unlimited linear Raman gain is influenced by insufficient temporal resolution in the computational grid. More specifically, extensive numerical simulations for a wide range of input parameters show that the collapse occurs when the intensity profile of the soliton is associated with features whose period is larger than the Nyquist frequency 1/(2 Δτ) of the simulation grid. This situation typically occurs when two spectral components whose frequency separation exceeds the Nyquist frequency temporally overlap. For example, in the particular case shown in Fig. 6 the frequency separation of the soliton and the corresponding trapped dispersive wave is approximately 171 THz which is larger than the Nyquist frequency of the simulation grid (170 THz in this case). In conjunction with the unbounded gain of the linear Raman model this violation of the Nyquist criterion yields significant aliasing of the Raman response part of the nonlinear operator which is the numerical mechanism that eventually leads to the collapse of the soliton. We have been unable to reproduce the collapse shown in Fig. 6 using the full model with finite Raman gain bandwidth although other forms of well-known numerical instability can of course be observed with poorly gridded simulations [11].

6. Conclusion

Understanding the effects of Raman scattering on ultrashort pulse propagation in fiber is essential to obtain a clear picture of the complex interactions that can occur during broadband SC generation. The main conclusion of this paper, however, is that the use of a linear approximation to the Raman gain curve in numerical studies may well result in inaccurate prediction of the expected temporal and spectral characteristics of the evolving SC field. By comparing simulations using the linear model and the realistic Raman response over a range of pulse durations and input soliton numbers, we have shown that depending on the parameter regime considered, the use of a linear approximation to the Raman response will significantly overestimate the degree of nonlinear spectral broadening and introduce qualitative and quantitative errors into both the predicted spectral characteristics and noise properties. Using a linear approximation to the Raman gain can of course yield important physical insights in theoretical studies, but it is important to fully appreciate its limitations and not to apply it in cases where it is likely to yield inaccurate results. As a general rule of thumb, when the bandwidth of the evolving spectral field approaches the bandwidth of the Raman gain and/or contains temporal structures shorter than 60 fs, the linear model is likely to lose validity and the simulation results must be considered with care. Ideally, for numerical studies, our conclusion is that the linear approximation should not be used at all; direct integration of an accurate model for the Raman response is numerically straightforward, and we recommend that this become the norm in numerical modeling of SC generation.