High-Order Analytical Formulation of Soliton Self-Frequency Shift

: We derive an analytical formulation of the Raman-induced frequency shift experienced by a fundamental soliton. By including propagation losses, self-steepening, and dispersion slope, the resulting formulation is a high-order (HO) extension of the well-known Gordon ' s formula for soliton self-frequency shift (SSFS). The HO-SSFS formula closely agrees with numerical results of the generalized nonlinear Schrödinger equation, but without the computational complexity and required computation time. The HO-SSFS formula is a useful tool for the design and validation of wavelength conversion systems and supercontinuum generation systems.


Introduction
Raman-induced frequency shift (RIFS) is an intrapulse scattering process in which the short wavelength content of a pulse gradually transfers to long wavelengths [1]. Solitons are particularly subject to RIFS, in which case we refer to soliton self-frequency shift (SSFS). SSFS is the main driving mechanism of several wavelength conversion and supercontinuum generation processes [1][2][3][4]. The amount of SSFS is accurately calculated by numerical simulation of the generalized nonlinear Schrödinger equation (GNLSE), but an analytical expression would be advantageous to simplify the computation task and provide insight about the impact of each parameter that influence SSFS. In 1986, J. P. Gordon reported an analytical description of SSFS derived from perturbation theory, in which high-order linear and nonlinear terms of the GNLSE were neglected [5]. Gordon showed that for an initially unchirped soliton, SSFS is well approximated from [1,5] where Ω is SSFS in / , 2 is group velocity dispersion coefficient, is Raman time, 0 is soliton duration, and is propagation distance. The beauty of Gordon s formula is that it explicitly links the amount of SSFS as a function of pulse and fiber parameters. For instance, it shows that SSFS increases proportionally with propagation distance and with inverse pulse duration to the fourth order. In most optical fibers, SSFS is significant for solitons with subpicosecond duration.
Although Gordon s formula provides a good estimate in many situations, a non-negligible amount of propagation losses, self-steepening and/or dispersion slope ( 3 ) lead to diverging results in between Gordon s formula and high precision estimates or experimental data [6][7][8][9]. While self-steepening leads to an SSFS that decelerates with propagation distance, dispersion slope can act in two opposite manners. For positive values of 3 , SSFS decelerates with propagation distance whereas the opposite occurs for negative values of 3 [10]. Since, selfsteepening and dispersion slope act simultaneously, their impact on the deviation of Gordon s formula depends on their relative strength. Taking the example of a dispersion flattened fiber, the effect of self-steepening dominates over the one from dispersion slope. This case is thus important for supercontinuum generation using dispersion flattened fibers [10][11][12]. On the other hand, dispersion slope has a leading influence over self-steepening e.g. in a high numerical aperture (HNA) fiber or a non-silica (e.g. ZBLAN and chalcogenide) fiber which often serve for wavelength conversion and supercontinuum generation [13][14][15][16][17][18][19]. Finally, Gordon s formula neglects the impact of propagation losses on SSFS which leads to misleading results e.g. when SSFS occurs in HNA silica fiber, ZBLAN fiber, or chalcogenide fiber [13,20].
In 2003, Santhanam et al. provided a general theory for RIFS based on moment method which includes high-order linear and nonlinear effects [21,22]. A limitation of the theory is that it omits the perturbation on the temporal pulse position resulting from the combined effects of RIFS and chromatic dispersion in the moment equation of chirp. As a result, the theory predicts a misleading RIFS when the temporal position of the soliton changes significantly with respect to the residual pulse position. This inconsistency was later resolved by Chen et al. with an improved moment method that considers a Raman delay that is pulse duration dependent [23]. In both the works of Santhanam et al. and Chen et al., SSFS is solved numerically from pulse evolution equations obtained by the moment method. Although a numerical solution of SSFS from the moment equations saves computation time in comparison to solving SSFS from the GNLSE, an analytical expression of SSFS remains a desirable option to quantify the influence of pulse and fiber parameters on SSFS. Indeed, calculation of SSFS from an analytical expression is much simpler and faster to process than any numerical solution.
In this paper, we propose a high-order (HO) SSFS formula that considers propagation losses, self-steepening, and dispersion slope. The HO-SSFS formula is derived from the moment method and an adiabatic approximation of soliton propagation. Solving for SSFS only from the moment equations requires a Taylor series expansion of the pulse parameters which results in an analytical form of SSFS that is cumbersome and diverges from the exact solution as the distance increases. Instead, we consider the adiabatic approximation along with the moment method thanks to the stable nature of soliton in presence of slowly varying fiber parameters along the fiber length [23,24]. The obtained HO-SSFS formula thus reflects acceleration and deceleration of SSFS with propagation distance resulting from propagation losses, self-steepening, and dispersion slope. The condition under which SSFS is dominated whether by self-steepening or by dispersion slope is also quantified. Finally, practical cases leading to SSFS are evaluated with the HO-SSFS, Gordon s formula, and numerical solutions from the moment method and GNLSE. In all cases, the HO-SSFS formula agrees well with numerical solutions, adding precision to Gordon s formula.

Generalized nonlinear Schrödinger equation
The evolution of an optical pulse in an optical fiber is well described by the generalized pulse-propagation equation [1] where is the slowly varying pulse envelope, is the propagation loss at the central frequency of the pulse, is the n th order dispersion coefficient, 0 is the waveguide nonlinear parameter, 0 is the central angular frequency, and ( ) is the nonlinear response function defined as [1] where ℎ ( ) = ( 1 -2 + 2 -2 ) 1 (− / 2 ) ( / 1 ) is the Raman response function including 1/ 1 as the phonon frequency associated with stimulated Raman scattering and 1/ 2 as the bandwidth of the Raman gain spectrum [25]. The parameter represents the fractional contribution of the delayed Raman response to nonlinear polarization. For pulses with duration much larger than 2 , a simplification can be performed by defining the Raman time as Considering a Taylor expansion of the propagation constant up to third order, Eq. (2) becomes [23] where = ( − / ) is a time variable in a frame of reference moving with the pulse at group velocity . At the right side of Eq. (5), the first term represents self-phase modulation, the second term is responsible for self-steepening, the third term accounts for intrapulse Raman scattering, and the last term represents the Stokes loss. Figure(1) shows a typical case of temporal and spectral evolution of a fundamental soliton in a 50 long HNA silica fiber. Fiber parameters are propagation loss = 10 / , 2 = −77.19 2 / , 3 = 3.1 3 / , and waveguide nonlinear parameter 0 = 2.08 −1 / at the initial pulse wavelength of 0 = 1.94 [15]. The pulse duration and peak power are 185 and 1.08 , respectively. The Raman time = 3 in silica as mesaured by Atieh et al. [26]. This simulation is preformed by numerically solving the GNLSE of Eq. (5) using split-step Fourier method [27]. As seen in Fig. 1, the soliton experiences a red-shift due to RIFS and simultaneously slows-down in the time domain because of this red-shift in an anomalously dispersive medium. It is also observed from Fig. 1(b) that SSFS decelerates with propagation distance and thus deviates from the linear evolution predicted by Eq. (1). This example will be further analyzed in section 4.2.1, but in the meantime justifies the need to further develop a more accurate analytical formulation of SSFS.

The moment method
The GNLSE is generally solved by numerical or semi-analytical techniques, with the knowledge of linear and nonlinear terms well-described in the frequency domain. The variational method is one of the semi-analytic techniques often used to solve the GNLSE of Eq. (5). However, this technique is not suitable to simulate a pulse propagation with RIFS because of the dissipative nature of the system in the presence of intrapulse Raman scattering [21]. Another approach to simulate RIFS in optical fibers is the perturbation theory [5]. However, for optical pulses in the subpicosecond range, the perturbation theory does not hold because the generated intrapulse Raman scattering becomes too large to be treated as a perturbation. The moment method is a powerful analytical technique for studying pulse propagation in any dissipative and non-dissipative systems because it does not require a Lagrangian as in the variational method [22]. With the moment method, an optical pulse is treated as a particle which maintains a specific shape irrespective of change in amplitude, width, and chirp during the propagation [28,29]. This is true for a fundamental soliton that maintains its shape thanks to the cancellation between the linear and nonlinear phase shifts. For the moment method, an optical pulse envelope is defined as which is a function of five moments: pulse energy , root-mean-square pulse width , temporal shift , chirp , and frequency shift Ω . The evolution of the pulse envelope at any distance is determined from an evaluation of each moments, following [23] In the anomalous dispersion regime ( 2 < 0), the solution to the GNLSE can be a fundamental soliton. Expressing the spectral dependence of linear and nonlinear parameters using a Taylor series expansion, this condition is still valid when dispersion terms of order larger than 2 and nonlinear terms of order larger than zero are negligible. Therefore, a stable fundamental soliton can also be maintained in presence of relatively weak high-order terms. Following the moment method, a fundamental soliton is expressed as [1] ( where and are the local duration and chirp, respectively of the soliton. They are related to the moments and by constant factors such that [23] By differentiating Eqs. (7)-(11) with respect to propagation distance and then using the GNLSE of Eq. (5) substituting with the field amplitude of Eq. (12) provide the evolution of the five moments as [23] = − In Eq. (5), the Taylor series expansion of ( ) and ( ) is considered up to third order and first order, respectively. However, this approximation is not appropriate for complex wavelength dependent dispersion and nonlinearity and large SSFS [23]. This limitation can be eliminated by considering a moving frequency reference frame proposed by Chen et al. [23]. With a moving frequency reference frame, the soliton central frequency is considered to follow a frequency shift such that 0, Also, the local values of 2 and are evaluated at the soliton central frequency as 2, 0, 0, By using Eqs. (20)- (22) in the moment equations (15)- (19), a new set of coupled moment equations is obtained in the moving frequency frame as , In general, Eqs. (23)- (27) are solved numerically because of their coupled nature, which requires significantly less computation time compared to solving the GNLSE from a split-step Fourier method. After solving for the five moments, the soliton after propagation is retrieved with Eq. (12). It is also possible to solve for the five moments analytically from Eqs. (23)-(27) by using a Taylor series expansion of the five moments as a function of propagation distance [30]. However, one problem with the Taylor series expansion method is that it gives a good fit with the numerical solution for a short propagation distance. For longer propagation distance, high-order terms are required to include in the Taylor series expansion, which results in substantially heavier analytical forms of the five moments.

Adiabatic approximation of soliton propagation
As discussed in the previous section, the primary assumption for the moment method to remain valid is that the pulse maintains its shape as it propagates down the optical fiber. For a fundamental soliton, this condition is satisfied when fiber parameters vary slowly with propagation length. This approximation is known as the adiabatic approximation of soliton propagation [2,23,31]. According to the adiabatic approximation, a propagating soliton maintains a balance in pulse width and peak power that satisfies the condition where is the soliton order and subscripts loc stand for local values. This adiabatic approximation also assumes that when a fundamental soliton experiences e.g. SSFS, its evolution is weakly affected by the high-order dispersive and nonlinear effects [2]. The perturbation of the highorder effects can be translated into a variation in the chirp of the soliton during propagation. Hence, a weak perturbation suggests that the chirp should remain nearly zero during the soliton propagation [23]. An example of numerical solution for the evolution of chirp of a fundamental soliton from the moment equations (23)-(27) is shown in Figure (2). As seen from Fig. (2), the variation of soliton chirp shows an oscillatory behavior with an average value of = −0.0026.This nearly zero average allows to apply the adiabatic approximation to the soliton propagation.
According to the adiabatic approximation of soliton propagation in Eq. (28), the pulse width of a fundamental soliton is written as By using the adiabatic approximation of in Eq. (26), the evolution of chirp becomes As seen from Eq. (30), adiabatic approximation suggests that an initially unchirped soliton = 0 will remain unchirped during its propagation. By using this condition of chirp in Eq. (27), the SSFS evolution becomes   As seen from Fig. (3), the evolution of pulse width from the moment method shows an oscillatory pattern around a monotonously increasing expansion of the fundamental soliton. This evolution is understood from the oscillatory behavior of the chirp as illustrated in Fig. (2). When

2,
< 0, the pulse compresses and the opposite occur when 2, > 0. On the other hand, the evolution of the pulse width according to the adiabatic approximation shows no such oscillatory pattern because of the continuous zero chirp approximation. The most noticeable feature of Fig. (3) is that the evolution of the soliton width from the moment equations oscillates around the evolution of pulse width obtained from the adiabatic approximation. Therefore, the expression of pulse width in Eq. (29) can be used with confidence in Eq. (31) to evaluate SSFS.

Analytic expression of SSFS
The analytical expression of SSFS is derived from the adiabatic approximation Eq. (29) and Eq. (31). After substituting the expression of pulse width from Eq. (29) in Eq. (31), the differential equation for the evolution of SSFS becomes By using Eqs. (21) and (22), Eq. (32) is converted into the integral form The evaluation of the integral on the left-hand side of Eq. (33) with the approximation that Ω << 0 and thus 1 + Ω / 0 ≈ Ω / 0 , as well as the use of the solution of Eq. (23) with Stokes loss ignored in the right-hand side of Eq. (33), result in a quartic equation of the form where the coefficients are as the effective length of the fiber. Because of the quartic form of Eq. (34), the resulting solution for Ω is cumbersome and of limited use in the context of soliton dynamics in optical fibers. However, the complexity in solution of Eq. (34) is greatly reduced if the effects of self-steepening and dispersion slope are taken separately. This approach is valid because in the presence of large dispersion slope the effect of self-steepening to SSFS can be ignored. Conversely, the effect of self-steepening to SSFS can only be observed when 3 is small such as in the case of photonic crystal fibers with ultra-flattened dispersion profile [32][33][34]. Next, the two cases will be considered separately to obtain the respective expressions of SSFS and then will be combined to get the HO-SSFS formula which covers both the effects of self-steepening and dispersion slope.

Case I: Effect of self-steepening
The effect of self-steepening on the SSFS dominates over the effect of dispersion slope when | 3 | << | 2 | / 0 (this inequality will be shown with (60)). In such a case, Eq. (34) reduces into a cubic polynomial equation as where, Searching for roots of Eq. (46), there is only one real root for which SSFS becomes Here, the subscript ss refers to the analytical expression of SSFS when only the effect of self-steepening is considered and 3 neglected. A Taylor series expansion of Ω provides where, Γ represents the Gamma function. To evaluate the SSFS deviation with respect to Gordon s formula due to self-steepening, it is necessary to develop at first a modified Gordon s formula which includes only the propagation loss. Ignoring both the effects of self-steepening and 3 in Eq. (33), the modified Gordon s formula is The modified Gordon s formula is also apparent as the first term of Ω in Eq. (48), while the second term is specific to self-steepening. The SSFS deviation with respect to the modified Gordon s formula is thus

Case II: Effect of dispersion slope
Multiplying the numerator and the denominator of the left-hand side of Eq. (33) with | 2 | 4 leads to When the effect of dispersion slope dominates over self-steepening, then | 3 | >> | 2 | / 0 and Eq. (52) simplifies into By evaluating the integral of the left-hand side of Eq. (53) and by using the solution of Eq. (23) with Stokes loss ignored in the right-hand side, the solution for the SSFS results into Here, the subscript beta3 refers to the analytical expression of SSFS when the dispersion slope dominates, and self-steepening is neglected. A Taylor series expansion of Ω 3 provides Take note that Ω 3 encloses the modified Gordon s formula and terms specific to 3 . The SSFS deviation with respect to the modified Gordon s formula is thus The analytical expression of SSFS including both contributions from self-steepening and dispersion slope is obtained from ΔΩ , ΔΩ 3 and Ω as Where 3 is obtained from the Taylor series expansion of Eq. (32) and consists of the coupled terms due to self-steepening and dispersion slope. This series of converging high-order terms can be ignored for the sake of conciseness of the expression with satisfactory level of precision. By considering SSFS in presence of propagation losses, self-steepening, and dispersion slope, we obtain the HO-SSFS formula (59) It is of interest to quantify a threshold where the high-order terms in HO-SSFS are dominated by self-steepening or by dispersion slope. A subtraction in between Eq. (57) and Eq. (51) leads to (60) where, In both cases, the Taylor series expansion of up to 6 ℎ order is enough to get an approximate value of SSFS. Table 1 shows the values of the coefficient with the order of Taylor expansion. In Eq. (60), when | 3 | = | 2 | / 0 and making the approximation of ∼ 1 for the first 6 orders in the Taylor expansion, then ΔΩ 3, − ΔΩ , ≈ 0 and both the self-steepening and dispersion slope have equal contributions on the evolution of SSFS. This therefore sets the threshold of equal magnitude for both self-steepening and dispersion slope. When | 3 | << | 2 | / 0 , self-steepening dominates the high-order contribution on HO-SSFS, and the effect of 3 can be neglected. On the other hand, when | 3 | >> | 2 | / 0 , dispersion slope dominates the high-order contribution on HO-SSFS and the effect of self-steepening to SSFS can be neglected.

Numerical results
To validate the HO-SSFS formula, a comparison is made with numerical solutions of the GNLSE of Eq. (5). The comparison is performed for the case of SSFS dominated by self-steepening (Case I), and for the case of SSFS dominated by dispersion slope (Case II). SSFS in presence of propagation losses is also considered in both cases.

Case I: Effect of self-steepening
The validity of the HO-SSFS formula dominated by self-steepening (Ω ) is verified by comparison with the GNLSE of Eq. (5), for the propagatioin scenario of a dispersion flattened fiber [35]. The fiber has a negligible 3 so that self-steepening dominates (i.e. | 2 | / 0 >> | 3 |) and SSFS evolution is expected to decelerate with propagation distance. Other fiber parameters are a length of 50 , 2 = −21.85 2 / , 3 where, is the number of frequency samples in the simulation, and˜ ( , ) is the Fourier transform of the pulse envelope ( , ) in Eq. (5). In all situations, Ω is the most accurate reference value and deviation from this value is considered an estimation error. The estimation error is calculated in between a target result Ω and the reference result Ω .  As seen from Fig. (4), numerical calculations of Ω and Ω agree well in this example as well as in most examples that will be evaluated next, at the exception of the case of SSFS dominated by a positive dispersion slope, which will be specifically discussed in section 4.2.2.
An error < 4% is observed for Ω (inset Fig. 4) due to the omission ( 3 ) terms in Eq. (59) and dispersive wave radiation by soliton. However, this error peaks within the few initial meters of propagation distance and converges towards 0 with further propagation. Take note that in this example, self-steepening is the only high-order contribution to the SSFS formula and thus Ω = Ω because 3 = 0. Self-steepening contributes to a deceleration (i.e. curvature, or saturation) of the SSFS evolution observed with Ω , Ω , Ω , and Ω . As observed from a comparison of Ω and Ω , the modified Gordon s formula allows for a deceleration of the SSFS evolution caused uniquely by the propagation losses, in contrast with Gordon's formula that allows SSFS to follow a linear function of propagation distance. Ω thus deviates from Ω in this example because of the combined effects of self-steepening and propagation losses.
To observe solely the effect of self-steepening on the evolution of SSFS, the same example is taken again but setting instead a propagation loss of = 0. Figure (5) shows the resulting evolution of SSFS.  Fig. (5), cancellation of the propagation losses did not cancel the tendency for Ω , Ω , and Ω to decelerate with propagation distance. In a lossless medium however, Ω (not shown) degenerates into Ω . Again, an error < 4% is observed for Ω (inset Fig. 5) due to the omission of ( 3 ) terms in Eq. (59) and dispersive wave radiation by soliton.

Case II: Effect of dispersion slope
To verify the validity of the analytical expression (Ω 3 ) and the HO-SSFS formula (Ω ), three different types of fibers are considered next. The first two types include an HNA silica fiber and a ZBLAN fiber which are widespread for supercontinuum generation. With their large dispersion slope, these two fibers are good candidates to show the impact of 3 > 0 on the evolution of SSFS with propagation distance. The third fiber is a silica-based dispersion flattened and dispersion decreasing fiber (DF-DDF) to show the impact of 3 < 0 on the evolution of SSFS with propagation distance. A first noticeable feature of Fig. (6) is the pronounced deceleration of the SSFS evolution, leading to an increasing deviation of Ω with respect to the linear evolution predicted by Ω . By comparison of Figs. (6a) and (6b), the effect of self-steepening is negligible over a short propagation length but plays a significant role that reduces asymptotically the error as the propagation length increases. It is also concluded that the dispersion slope is the main element responsible for the strong deviation in between of Ω and Ω . Despite this deviation, the error of Ω is maintained to < 2.5%.
Results of Fig. (7) for the ZBLAN fiber lead to conclusions that are similar to the ones of the HNA fiber. That is, Ω has an error that peaks to around 3% after a short propagation length but this error reduces asymptotically towards zero with propagation distance. Dispersion slope is the main factor responsible for the deviation of Ω with respect to the linear evolution predicted by Ω . . In (a), Ω 3 is shown with inset that provides the error of Ω 3 . In (b), Ω is shown with inset that provides the error of Ω .

Effect of dispersion slope with 3 < 0
In some cases, the slope of 2 with respect to frequency is negative such as in a DF-DDF silica fiber with two zero-dispersion wavelengths [10][11][12].This type of fiber has unique properties as described in [14] and is being used for supercontinuum generation. For this case, fiber parameters are a length of 50 , 2 = −5.29 2 / , 3 = −0.13 3 / , 0 = 5 −1 / at 0 = 1.7 . Propagation losses are ignored to specifically observe the effect of 3 < 0.The pulse is a fundamental soliton with central wavelength of 0 = 1.7 , with a duration of 100 and a peak power of 106 . Figure (8) shows the evolution of SSFS in the DF-DDF fiber. As seen in Fig. (8), a negative dispersion slope accelerates the SSFS evolution as observed from Ω , with respect to a linear evolution such as the one predicted by Ω . The HO-SSFS formula Ω remains in good agreement with Ω , providing an error up to 2.5% but diverging with propagation distance. It is worth noting that the SSFS acceleration cannot be maintained over an unlimited propagation distance because the soliton will eventually reach the zero-dispersion wavelength, and thus the soliton itself would not be preserved.
Another important feature of Fig. (8) is the observable discrepancy in between Ω and Ω . It is assumed that such discrepancy must arise from dispersive wave radiation of the soliton. When a soliton is weakly perturbed by high-order linear and nonlinear effects, it is pushed away from its stable condition and sometimes returns to equilibrium by means dispersive wave radiation [38]. The emission of dispersive radiation results in a damped oscillation in the evolution of pulse duration and chirp of the soliton [23]. Such an emission of dispersive radiation is automatically considered in the numerical solution of Ω because the GNLSE can deal with arbitrary and transitory pulse shapes [23]. On the other hand, within the moment method or adiabatic approximations, the soliton is considered to maintain its shape during the propagation. As a result, the moment method or adiabatic approximation omits the emission of dispersive radiation, which ultimately leads to the discrepancy between the results of SSFS calculated from the analytical expressions and numerical simulation of GNLSE. Since the HO-SSFS formula stems from a combination of the moment method and adiabatic approximation, then it is expected that the error on Ω partly originates from the error on Ω , in addition to the omitted 3 terms in Eq. (59).
The impact of propagation losses on the evolution of SSFS from Figs.(4), (6), and (7) depicts the importance of its inclusion when designing wavelength converter based on Raman frequency shift. Since the effective length of the fiber in this case depends on four times the propagation loss as seen from Eq. (40), a small amount of propagation loss can also contribute significantly on the evolution of SSFS even in the absence of high order linear and nonlinear effects.

Conclusion
We have presented a high-order (HO) extension of the well-known Gordon s formula for soliton self-frequency shift (SSFS), using an adiabatic approximation of soliton propagation combined with the moment method. Taking solely into account the propagation losses of a soliton leads to a modified Gordon s formula, into which propagation distance of Gordon s formula is replaced by an effective propagation length. Furthermore, an analytical expression of HO-SSFS was obtained when including the wavelength dependence of the waveguide nonlinear parameter (i.e. self-steepening) and the wavelength dependence of chromatic dispersion (i.e. dispersion slope). To verify the accuracy of the HO-SSFS formula, numerical simulations were carried out for two cases where self-steepening was the dominant HO term of the HO-SSFS formula, as well as three cases where dispersion slope was the dominant HO term of the HO-SSFS formula. Numerical simulations were considered with and without propagation losses to highlight their importance on the evolution of SSFS. The comparison in between results of the HO-SSFS formula and numerical solution solving the generalized nonlinear Schrödinger equation are showing an excellent agreement within an error < 3%.The error is mainly due to the omission of crossed high-order terms that are not considered in the HO-SSFS formula.
The HO-SSFS formula is a useful tool for the design and validation of optical devices making use of SSFS, such as wavelength converters and supercontinuum sources. It provides physical insight about every fiber and pulse parameters leading to SSFS, as well as it leads to an accurate estimate of SSFS without the time and effort of numerical simulations systems.