Frequency-banded nonlinear Schrödinger equation with inclusion of Raman nonlinearity

The well-established generalized nonlinear Schrödinger equation (GNLSE) to simulate nonlinear pulse propagation in optical fibers and waveguides becomes inefficient if only narrow spectral bands are occupied that are widely separated in frequency/wavelength, for example in parametric amplifiers. Here we present a solution to this in the form of a coupled frequencybanded nonlinear Schrödinger equation (BNLSE) that only simulates selected narrow frequency bands while still including all dispersive and nonlinear effects, in particular the inter-band Raman and Kerr nonlinearities. This allows for high accuracy spectral resolution in regions of interest while omitting spectral ranges between the selected frequency bands, thus providing an efficient and accurate way for simulating the nonlinear interaction of pulses at widely different carrier frequencies. We derive and test our BNLSE by comparison with the GNLSE. We finally demonstrate the accuracy of the BNLSE and compare the computational execution times for the


Introduction
Numerical simulations of pulse propagation are an integral part of investigating dispersive and nonlinear effects within optical fibers.One of the most popular models used for this is the generalized nonlinear Schrödinger equation (GNLSE), which in combination with the split-step Fourier method (SSFM) provides an accurate and efficient tool for investigating dispersive, Kerr and Raman effects of pulses throughout propagation in optical fibers [1].
The efficiency of simulating the GNLSE using the SSFM is attributed to the execution speed of the fast Fourier transform (FFT) algorithm.However, for it to be effective the time-frequency grid is required to be linearly spaced.This constraint can prove detrimental to the execution times when the frequency bandwidth considered is very large.This is particularly adverse to simulating fiber optical parametric amplifiers (FOPA) where the effects that are of interest are concentrated in small frequency bands that can be spaced over large separations [2].Because of this, the GNLSE is usually split between multiple coupled nonlinear Schrödinger equations (CNLE) that simulate pulse propagation within certain frequency bands [3].However, these coupled equations usually disregard the Raman effect by setting its fractional contribution to the total nonlinearity to zero.This results in the entire fiber nonlinearity being attributed to the Kerr effect which is a misrepresentation of the optical system and leads to inaccurate results depending on whether the frequency bands are within or outside the Raman gain bandwidth as it ignores the frequency dependence of the nonlinearity [4].
There have been a variety of studies that take into account the Raman nonlinearity when considering coupled nonlinear equations.An investigation of stimulated Raman scattering (SRS) [5] proved that the influence of the Raman susceptibility affects the strength of cross-phase modulation in the coupled wave equations.Following from this, [6,7] produced coupled equations that incorporated Raman associated factors when investigating four-wave mixing and SRS in continuous waves (CW).Variations of these equations were used to investigate the combined Kerr and Raman effects for CW and quasi-CW pulses [8][9][10][11].The coupled equations can also be simplified when the pump is considered to be undepleted throughout propagation [4,[12][13][14] or when its power is much larger than the probes considered [15].These models have been used to explain influences of Raman scattering on noise figures [9,16] and demonstrate the need for the Raman factor to be considered when investigating parametric amplification.
Here we present a numerical method capable of simulating the propagation of pulses that nonlinearly interact with each other but occupy widely separated frequency bands.The model provides a unified framework for investigating self-phase modulation (SPM), cross-phase modulation (XPM), four-wave mixing (FWM), self-steepening and dispersive effects over multiple frequency bands.Additionally, the associated coefficients of the nonlinear effects are derived while taking the Raman nonlinearity into consideration.First the resulting equation, dubbed the banded nonlinear Schrödinger equation (BNLSE), is derived from the GNLSE in Section 2. Following this, the model is solved by using the split-step Fourier method (SSFM) and the results are compared to the GNLSE and previously published CNLSEs in Section 3.

Theoretical model
The evolution of the amplitude A(z, t) of a pulse that is propagating through an optical fiber can be found by solving the GNLSE [1]: where β n is the n th Taylor expansion coefficient of the propagation constant β calculated at the central frequency ω 0 , γ the nonlinear coefficient, f r the fractional contribution of the Raman nonlinearity and h (t) the Raman response function.
The BNLSE can be derived for any number of frequency bands.For clarity here we only consider three such bands, l = −1, 0, 1, centered at angular frequencies Ω −1 < 0, Ω 0 = 0, and Ω +1 > 0, respectively, relative to the central frequency band.The amplitude of the GNLSE is then split as: We apply Eq. (2) to Eq. ( 1) and on both sides of the resulting equation identify the terms with time dependencies close to exp(−iΩ l t), l = −1, 0, 1, to obtain the evolution equations of the amplitudes A −1 , A 0 , and A +1 , respectively.Two approximations are made at this stage.First, terms rotating at frequencies 2Ω ±1 , |Ω +1 | + |Ω −1 |, etc. are neglected as they correspond to non-energy conserving FWM terms.Secondly, we simplify terms involving the Raman response function by This assumes that the individual amplitudes A l are varying slowly compared to the Raman response function.In other words, intra-band Raman effects are neglected.The second line in Eq. ( 3) is obtained by identifying the integral in the first line as the Fourier transform h(ω) of the function h(t).For silica optical fibers, neglecting intra-band Raman effects is justified for pulse lengths larger than 1 ps [1].With these assumption, the amplitudes A l can then be found by solving the BNLSE: with where γ l , ω l are the nonlinear coefficient and central angular frequency of band l.
In contrast to the GNLSE, the BNLSE has two types of nonlinear interactions, intra-band nonlinearities which are incorporated by the SPM term and inter-band nonlinearities which are dictated by XPM and FWM terms.The inclusion of the inter-band Raman response constitutes a distinct enhancement to previously presented CNSLEs [3]; it allows for arbitrary positioning of the bands anywhere inside or outside the Raman gain bandwidth and never disregards the fractional contribution of the Raman nonlinearity.As can be seen in Eq. ( 4), when the bands are separated by frequencies much larger than the Raman gain bandwidth (and therefore h (ω) ≈ 0) there is no Raman contribution to the inter-band nonlinearity, while the correct fractional contribution is applied for smaller frequency separations.It is also worth mentioning that, like the GNLSE, the BNLSE conserves the number of photons throughout propagation.Finally, if f r is set to zero, the BNLSE reverts to the CNLSE.
As already discussed, some approximations were made in the derivation of the BNLSE.In particular, intra-band Raman effects were neglected.For shorter pulses this is not always justified and in principle it is possible to extend the BNLSE to also include intra-band Raman effects.For example, the simulation of Raman-induced soliton self-frequency shifting would require this modification.Moreover, spontaneous Raman scattering which generates frequency components of light outside the selected bands is not included in this theory.
Additionally, the BNLSE can only include FWM processes that occur within or between the discrete bands.Equations ( 4)-( 7) were derived allowing for different frequency spacings between the bands.This leads to the explicit time dependence of the inter-band FWM terms in Eq. ( 7) and may lead to the generation of signal or idler waves outside the chosen bands.However, such processes depend on the fiber group velocity dispersion profile as they require phase matching and thus in many situations are not of practical importance.In such cases the form of the equations given above is the most generic form possible.
For the numerical investigation of wideband degenerate FWM, which we consider one of the most important applications of the BNLSE, one would generally choose bands that are equally spaced in frequency, i.e., Ω +1 = −Ω −1 = Ω > 0. Then the explicit time dependence in Eq. ( 7) is removed and the signal and idler waves generated by a strong pump in band 0 will be located in the two bands ±1.Note, however, that cascaded FWM effects or the generation of phase matched signals outside the chosen frequency bands is still disregarded.
From this discussion it is clear that a careful selection of the bands' positions and bandwidths (and potentially increasing the number of bands considered) is required for efficient and accurate computation of nonlinear interactions over wide bandwidths.To some degree, this requires prior understanding of the fundamental physics of the system under investigation.We emphasize, however, that this is also true for the full GNLSE where appropriate discrete grids in time and frequency must be chosen to encapsulate the correct pulse propagation dynamics.
In the following section we validate the BNLSE by calculating the phase modulation and FWM between bands and comparing the results with those obtained by the GNLSE.

Comparison of numerical models
For the purposes of comparing the GNLSE, CNLSE and BNLSE all models are used to simulate pulse propagation of a silica based photonic crystal fiber (LMA5-PM from NKT Photonics).While the fiber is birefringent we only consider propagation on its slow axis.The nonlinear coefficient is γ = 10 W km −1 , the zero dispersion wavelength is at λ z = 1051.85nmand at that wavelength the dispersion coefficients are found to be β 3 = 6.756 [17].Since the fiber is silica based the Raman response function h(t) is interpolated from the data provided in [18].We focus here solely on the case of equally spaced bands with frequency separation F (angular frequency separation Ω = 2πF) as discussed at the end of Sec. 2.
The three models solve Eq. ( 1) with f r = 0.18 for the GNLSE and Eq. ( 4) with f r = 0 and 0.18 for the CNLSE and BNLSE, respectively.While the models could be solved by a general ordinary differential equation solver, it was found that the symmetrized SSFM was more efficient in execution time if the required numerical accuracy is not too high.Therefore in all cases the pulse propagation is undertaken using the symmetrized SSFM with a 4 th order Runge-Kutta method [19] applied to the nonlinear step.In all the simulations considered the pulse is propagated for 14 m at a constant propagation step dz = 10 −3 m which was found to be small enough for each individual model to converge.By selecting a constant propagation step the numerical models considered can be directly compared to each other throughout the fiber propagation.Finally, in order for the methods to have the same frequency grid spacing d f , frequency bands from the frequency spectrum of the GNLSE were selected for the CNLE and BNLSE and the Raman response shifts, h(Ω) etc., were calculated accordingly.
The source code used to compare the models is available under the MIT license at [20].The software is written in Python 3.6 while making use of the NumPy, SciPy and Matplotlib libraries from the Python scientific stack [21,22].Additionally, the software makes selective use of the Numba compiler to enhance its execution speed [23].Finally, the Intel distribution for Python is used for the increased efficiency of the FFT algorithm [24].
In the following subsections, the accuracy of the BNLSE and CNLSE is compared with respect to the GNLSE.In particular the effects the Raman factor has on XPM and FWM are discussed and, finally, a comparison of the execution times between the models is presented.

Nonlinear phase modulation
The addition of the Raman dependent factors in Eq. ( 4) indicates that the nonlinear phase of the pulses will vary depending on the value of the Raman contribution.Here, the nonlinear phase of the pulses is investigated and the accuracy of the BNLSE/CNLSE is compared to the GNLSE.The total phase of a pulse can be calculated from φ (z, t) = arg (A (z, t)) which encompasses both the linear and the nonlinear phase.Since the linear phase of the pulses is φ L = exp (iβz), where β is the propagation constant of the pulse at its carrier frequency, the nonlinear phase can be calculated from φ N L = arg (A (z, t) exp (−iβz)).
In all models two pulses are launched at the fiber input, a strong pump pulse with peak power P 0 = 10 W and a weak probe pulse with peak power P 1 = 100 mW.Both pulses have a half width duration of T 0 = 50 ps.The pump wavelength is set at λ p = 1051.35nm (−0.5 nm from the zero dispersion wavelength).For these parameters, degenerate FWM is phase matched at F =15.47 T Hz (997.24nm and 1111.67 nm wavelength), i.e., within the Raman gain bandwidth.
However, for the length of fiber under consideration and the chosen pump power, no significant energy transfer to these phase matched wavelengths is observed.Instead, we investigate the nonlinear interaction of the pump with the weak seeded probe pulse at a non-phase matched wavelength.Initially the probe is placed outside the Raman gain bandwidth at F = 60 THz (868.59 nm wavelength), so that the Raman gain and FWM effects have only negligible effects on the phase and amplitude of the pulses, Fig. 1.
Because of the orders of magnitude difference between the power of the two pulses, the dominant nonlinear phase modulation on the pump will be from SPM.It was found that regardless of the value of the Raman factor the GNLSE, CNLSE and BNLSE produced the same results.This was expected because the SPM factor in Eq. ( 4) is not dependent on f r .On the other hand, the dominant effect on the phase of the probe will be XPM from the pump.Fig. 1 depicts (from left to right) the (nonlinear) phase of the probe pulse calculated by using the GNLSE and the relative error between this solution and those calculated using the CNLSE and BNLSE.As can be seen, the relative error of the CNLSE exhibits a trend that is similar to the nonlinear phase of the probe.This is a direct result of the XPM factors in Eq. ( 4) which during propagation have a compounding effect upon the nonlinear phase change of the pulse.On the other hand, the relative error of the BNLSE with respect to the GNLSE is much smaller, < 1%.This remaining numerical deviation was found to be related to the frequency grid spacing of both models and was reducing for smaller grid spacing d f .However, the huge number of grid points (and thus computer memory and execution time) required in the GNLSE limited the study of small values of d f .
Equivalent results are depicted in Fig. 2 when the probe is placed at a frequency separation of F = 15.2THz (998.14 nm wavelength) so that the bands fall within each other's Raman bandwidth and thus pulse propagation is affected by SRS, Raman phase modulation and FWM.As in the previous case, the phase of the probe using the CNLSE exhibits a relative error with a trend similar to the nonlinear phase modulation, which again is due to the difference of the factors in Eqs. ( 5) and ( 7) when f r is set to 0 for the CNLSE.In contrast, the relative error of the BNLSE is consistently much smaller and again is computationally limited by the size of the frequency grid of the GNLSE.Thus, Figs. 1 and 2 show the accuracy of the BNLSE in calculating the (nonlinear) phase for bands both outside and inside the Raman gain bandwidth.

Four-wave mixing
The CNLSE has been widely used in FOPA systems to study FWM.However, as is presented in this subsection the disregard of the Raman gain factor substantially reduces the accuracy of the model with respect to the equivalent results produced by the GNLSE.
As in subsection 3.1 the comparison of the models with respect to FWM is conducted with a probe pulse placed inside and outside the Raman gain bandwidth of a strong pump pulse.The powers and pulse widths are the same as in subsection 3.1 but the position of the pulses on the frequency grid is changed.In both cases investigated below the probe position with respect to the pump is determined by solving the phase matching equation [1]: where Since the probe will be phase matched to the pump this will result in the generation of an idler wave at ω i = ω p − Ω = ω p − 2πF where F is the frequency difference between the pump and the probe wave.The pump of the initial simulation is set at λ p = 1051.35nm which results in a phase matched frequency separation of F = 15.5 THz (probe wavelength 997.15 nm and idler wavelength 1111.78 nm).The second simulation has a pump situated at λ p = 1046.85nm which gives F = 42.5 THz (probe wavelength 911.57nm and idler wavelength 1229.28 nm).The probe and resulting idler of the first simulation are situated within the Raman gain bandwidth, while for the second simulation they lie outside.Fig. 3 shows the average power of the idler wave versus propagation length for the probe placed inside (left) and outside (right) the Raman gain bandwidth by using the GNLSE, CNLSE and BNLSE.At the fiber input, the idler power vanishes but it increases during propagation due to FWM and Raman gain.The CNLSE neglects the Raman nonlinearity but overestimates the Kerr nonlinearity, which results in a significant deviation of the idler power with propagation distance compared to the exact GNLSE.On the other hand, the BNLSE is consistently accurate with respect to the GNLSE throughout the pulse propagation independent of the position of the individual frequency bands relative to the Raman gain bandwidth of the pump.An interesting effect that is visible in Fig. 3 is that the results from the CNLSE are closer to those calculated by the GNLSE/BNLSE for bands outside the Raman gain bandwidth rather than inside.This is caused by the negative value of the real part of h(±Ω) for the chosen value of Ω which implies that the Raman nonlinearity reduces the FWM gain in the GNLSE/BNLSE, according to the last line of Eq. (7).F = 42.5 THz, respectively, for different propagation distances.As in the case of the idler energy, the BNLSE consistently produces identical results to the GNLSE throughout propagation.This adds to the conjecture that the BNLSE is able to accurately and efficiently simulate the same dispersive and nonlinear effects as the GNLSE.For example, in Fig. 5, the BNLSE predicts the same time shifting of the idler with respect to that of the pump owing to the presence of group velocity dispersion.

Computational execution times
The main incentive behind the derivation and use of the BNLSE is to avoid the large computational execution times and memory demands of the GNLSE for large frequency grids required for the simulation of narrow bandwidth but widely separated pulses.To support this conjecture the execution times t e for varying frequency grid spacings d f (and thus number of grid points) are presented in Fig. 6 for the three models considered.For each simulation considered the frequency window of the GNLSE is 170.2THz wide and the frequency bandwidths of the CNLSE/BNLSE are 5.3 THz.The simulations are conducted on an Intel Core i7-4790 CPU at 3.60 GHz and each time measurement represents the average over 100 repetitions of the simulation.
Within Fig. 6 the GNLSE is considered with and without a Raman factor ( f r = 0.18 and 0, respectively).This is because in order for the Raman gain to be calculated in Eq. (1) a convolution needs to be evaluated at each step.This significantly increases execution times due to the increased number of FFTs required.Therefore both situations are considered here in order to accurately compare the GNLSE to the BNLSE/CNLSE.For all the frequency grid spacings considered the BNLSE consistently exhibits smaller execution times than the GNLSE.In particular, for the smallest d f considered the BNLSE is found to be over 6 and 10 times faster than the GNLSE for f r = 0 and 0.18, respectively.The increase in efficiency that the BNLSE/CNLSE provide is mainly because the size of the arrays required is much smaller than for the GNLSE.
The CNLSE and BNLSE exhibit similar execution times for the same grid sizes with the CNLSE being slightly faster because in this case f r is set to zero in Eq. ( 4), which allows us to simplify the equation analytically and thus make computation faster.However, as d f becomes smaller (i.e. the grid sizes larger) the difference is decreased because the FFTs and matrix multiplications shown in Eq. ( 4) have an increasingly larger effect than the incorporation of the additional factors.Nevertheless, this small increase in execution times is nominal compared to the increase in accuracy that the BNLSE provides.

Conclusion
The derivation of a frequency-banded nonlinear Schrödinger equation (BNLSE) with the inclusion of inter-band Raman effects has been presented.Its accuracy in predicting nonlinear phase modulation and FWM against the GNLSE and the CNLSE was investigated and the BNLSE was found to be substantially more accurate than the CNLSE in describing nonlinear phenomena in fibers.Furthermore, the numerical execution times of all the models have been presented and compared to each other with the BNLSE portraying similar execution times as the CNLSE but substantially faster than the GNLSE.
By reducing the size of the frequency grid considered and introducing the Raman nonlinearity dependent corrections, the BNLSE is an efficient and accurate alternative to the GNLSE.Immediate applications of the model are found in fiber optical parametric amplifiers and oscillators where energy is converted over large bandwidths [2,25].Furthermore the equation can be used to simulate mid-IR generation [26,27].
It is also worth mentioning that the BNLSE can be extended to incorporate more than three bands by increasing the number of amplitudes considered in Eq. ( 2).This would be beneficial to high power applications where cascade FWM processes are of importance [28].Finally, the equation could also be extended to a multi-mode version of the BNLSE to efficiently simulate parametric conversion between multiple fiber modes over extended frequency bands [29,30].

Fig. 1 :
Fig. 1: Nonlinear phase of the probe at a frequency separation of F = 60 THz (outside the Raman gain bandwidth of the pump) throughout propagation calculated using the GNLSE (left) and the relative error of the CNLSE (middle) and BNLSE (right) compared to the GNLSE.All parameters are given in the text.

Fig. 2 :
Fig. 2: Nonlinear phase of the probe at a frequency separation of F = 15.2THz (inside the Raman gain bandwidth of the pump) throughout propagation calculated using the GNLSE (left) and the relative error of the CNLSE (middle) and BNLSE (right) compared to the GNLSE.

Fig. 3 :
Fig. 3: Average power of the idler wave with respect to the propagation distance with probe inside, at F = 15.5 THz, (left) and outside, at F = 42.5 THz, (right) the Raman gain bandwidth.

Fig. 4 :
Fig. 4: Temporal (left) and spectral (right, on a logarithmic scale) distribution of the idler wave with probe inside the Raman gain bandwidth, at F = 15.5 THz.