1 Introduction
Nonlinear optics is a branch of optics dealing with phenomena that are consequences of lightinduced modifications to the optical properties of matter. Optical fibers are a very attractive medium for studying nonlinear phenomena. Due to the small size of the core and the geometry of the fibers, which allows light to propagate over very long distances, the waveguides can show nonlinearity even for the relatively low power of the initial pulses. Mathematical equations provide an excellent description of the nonlinear effects observed in this type of medium. Unfortunately, usually they are so complicated that they do not have an analytical solution, and the use of numerical integration techniques is unavoidable.
The numerical tools offer an excellent opportunity for deep inspection of light propagation in the computational investigation of nonlinear phenomena occurring in optical fibers. They play a crucial role in supporting various types of laboratory experiments. First of all, numerical tools allow reducing the cost of experimental work. Moreover, they give access to all information about the propagating pulse on the whole propagation distance, not only measured at the output. It is also possible to switch different terms of the equation being solved on and off (including and excluding some specific phenomena) to understand the underlying physics.
The basic mathematical model describing the shaping of the envelope of the pulse amplitude is a NonLinear Schrödinger Equation (NLSE), which depends on the time and distance of propagation. It takes into account effects such as the second order dispersion and phase selfmodulation. It is successfully used in modeling the picosecond propagation of lowenergy pulses propagating in optical networks. However, it is insufficient for femtosecond high energy pulses. To improve the existing NLSE model, the dispersion of higher orders and expressions describing successive nonlinear processes (e.g., Raman scattering) were added, yielding the GNLSE [1]. Numerical modeling of the propagation of light pulses in various optical fibers based on GNLSE is used to describe complex phenomena such as supercontinuum generation. In practice, the equation facilitates the design of optical fibers, allowing for the generation of coherent radiation.
Considering the fact that the numerical solution of this complicated partial differential equation is quite timeconsuming, it is necessary to develop efficient algorithms to solve it. High accuracy and relatively small inference time of the constructed methods is required. In the report, we introduce the open source Python toolbox for solving the GNLSE, based on the splitstep Fourier transform method (SSFM), which takes the least computation time among all compared numerical schemes for GNLSE when the solution varies slowly with time [2].
2 Numerical model
In general, the GNLSE is an example of an equation that generally does not have analytical solutions. Moreover, the use of a numerical approach is necessary to solve the GNLSE and to describe the nonlinear processes occurring in optical fibers. In the time domain, the GNLSE takes the form
(1) 
where is the attenuation constant, and are the higher order dispersion coefficients obtained by a Taylor series expansion of the propagation constant around the center frequency . The term on the second line describes the nonlinear effects – time derivative in this term is responsible for selfsteepening and optical wave breaking, whereas the convolution integral describes the delayed Raman response [3]. This form of the GNLSE is commonly employed for numerical simulations of pulse propagation in a nonlinear medium such as an optical fiber.
Eq. 1
can also be expressed in the frequency domain by overlaying the Fourier transform as
[1](2) 
The first and second terms on the righthand side of Eq. 1 and Eq. 2 indicate the dispersion and nonlinear terms, respectively.
The SSFMs are a class of numerical methods that allow for the study of the behavior of a light pulse propagating in an optical fiber. To better understand the idea behind the algorithm, Eq. 1 can be written in the following form
(3) 
where is the dispersion operator showing the influence of higher order dispersion and attenuation on the pulse propagation, while is the nonlinear operator describing nonlinear phenomena occurring in the optical fiber. They are properly defined in both the time and frequency domains.
In general, both dispersion and nonlinearity act simultaneously along the entire length of the fiber. However, in the SSFM, we assume that they act independently at a short distance (). In other words, the propagation along the segment can be modeled first by using the nonlinearity operator (we assume that ), and then the dispersion operator . This leads to the approximation [1]
(4) 
The dispersion operator is easier to calculate in the frequency domain to which we transfer using the Fourier transform. The return to the time domain is performed using the inverse Fourier transform. In the case of SSFM, the linear (dispersive) term is exactly accounted for, as it is completely integrable, whereas the nonlinear term is usually integrated using RungeKutta (RK) algorithms. A detailed description about the technique employed here can be found in [3, 4].
3 Toolbox deliverables
In our framework, we proposed an efficient solver, thanks to an adaptive step size implementation of the fourthorder RungeKutta in the Interaction Picture method (RK4IP), adapted from [3, 4]
. Numerical integration of a system of ordinary differential equations is done by
solve_ivp from the scipy.integrate package. Our Python code is partially based on MATLAB code published in [4], and available at http://scgbook.info/. The toolbox prepares the integration using SciPy’s ODE solver, while the transitions between time and frequency domains are accomplished using the FFT and iFFT from pyfftw library.The solver module is divided into three parts:

specific model setup (gnlse.GNLSESetup),

the framework for preparing the splitstep Fourier algorithm (gnlse.GNLSE),

and a class for managing solutions (gnlse.Solution).
3.1 Initial conditions: pulse envelopes
To follow the evolution of the amplitude envelope in time and frequency, it is necessary to set appropriate boundary conditions. In this case, it is the known initial pulse shape described by the functions or . They represent the complex envelope amplitude of the pulse at the fiber input.
In physics, the envelope of an oscillating signal is a smooth curve outlining its extremes. The given envelope is a function of time, describing the amplitude of the input pulse calculated basing on its specific properties, i.e., pulse duration fullwidth halfmaximum and peak power (or only average power for a continuous wave). Fig. 1 illustrates the envelopes of four types of optical signals implemented in the envelopes module. We allow the user to choose between Gaussian, Lorentzian, hyperbolic secantshaped envelopes and a continuous wave with optional noise.
3.2 Dispersion operator
Dispersion is a phenomenon that manifests itself as the dependence of the phase velocity of a wave on other quantities, such as frequency or polarization. One of the most important parameters characterizing optical fibers is the chromatic dispersion, which expresses the dependence of the signal propagation velocity on the wavelength. Chromatic dispersion plays a huge role in the propagation process of a pulse in an optical fiber. Different frequency components have different propagation constants , which can be expressed by the equation [1]
(5) 
Since the frequencies present in timelimited pulses are concentrated around the central frequency , we can use the expansion of into the Taylor series
(6) 
where
(7) 
For such representation, the dispersion operator in the frequency domain is defined as
(8) 
where the second term of the righthand side is responsible for loss.
The and parameters are defined as
(9) 
(10) 
where is the group velocity at which the pulse propagates, and is the group refractive index. Thus, is the reciprocal of the group velocity, and quantifies the broadening of the pulse due to the fact that not all of its components propagate at the same velocity. This phenomenon is called group velocity dispersion (GVD). Higher orders of dispersion are also defined; however, for pulses with narrow spectral bands, they can be omitted. In most cases, it is enough to take into account the parameter, which is responsible for the phenomenon of thirdorder dispersion. Taking higher order dispersion into account is usually necessary in the case of extreme spectral broadening, such as in the supercontinuum generation process.
The linear operator in frequency domain can be also expressed by [1]
(11) 
where is both the center frequency of the defined grid and the central frequency corresponding to the wavelength of the central component of the spectrum of the introduced pulse.
The first term of the expression above corresponds to the Taylor series expansion (starting from the second term) of the propagation constant
(12) 
The dispersion operator defined by Eq. 11 is extremely useful when we know the exact values of the propagation constant , not just its successive derivatives. In the case of a known dependence of the effective refractive index on the frequency, it is enough to make appropriate calculations using the formula
(13) 
The next step is to calculate the value of and its first derivative at .
Both forms of the dispersion operator (namely, Eq. 8, and Eq. 11) have been implemented in the described toolbox. In case of taking into account the full dependence of propagation constants on frequency, we use an extrapolation method for the provided effective refractive indices (and calculate propagation constants from them) from an external file. For this purposes the interp1d function from SciPy library was used. A demonstration of the influence of accounted input parameters is shown in Fig. 2 for the supercontinuum generation example.
3.3 Nonlinear phenomena: Raman responses
The last term in Eq. 1 describes nonlinear phenomena such as: phase selfmodulation, selfsteepening, shock wave formation, and Raman scattering. The function, which determines the last of the listed effects, usually is modelled as [5]
(14) 
where
(15) 
is a Heaviside step function, while , and depend on the type of material. For silica glass, they equal , and .
The response is calculated based on the chosen Raman model. There are three available Raman response models for silica optical fibers implemented:

and based on Q. Lin and Govind P. Agrawal model [6] as a more accurate response function,

based on Dawn Hollenbeck and Cyrus D. Cantrell model [7], which presents the multivibrationalmode function,
Their visualisations in the time domain are presented in Fig. 3.
3.4 Nonlinear phenomena: The nonlinear coefficient
The time derivative term inside GNLSE models the dispersion of the nonlinearity. This is usually associated with effects such as selfsteepening and optical shock formation, characterized by a timescale . In the context of fibre propagation, an additional dispersion of the nonlinearity arises due to the frequency dependence of the effective mode area. The last effect can be accounted in coefficient in an approximate manner [5, 8].
A better – still approximate – approach to include the dispersion of the effective mode area is to describe it directly in the frequency domain [9]. In this case, we can derive a GNLSE for the pulse evolution using defined as
(16) 
This approach is more rigorous than the approximation of () and requires the definition of a pseudoenvelope as
(17) 
Both forms of the equation (namely, Eq. 1, and the form introducing scaling Eq. 17) have been implemented in the described toolbox. In the case of taking into account the full dependence of the nonlinear coefficient on frequency, we used an extrapolation method for the provided effective mode areas (and effective refractive indices) from an external file using interp1d function from SciPy library. Investigation of the influence of accounted input parameters is shown in Fig. 4 for soliton fission example.
3.5 Additional tools
We provide some additional tools and functions inside gnlsepython package. As MATLAB is the most popular environment for nonlinear computations in the scientific community, we provided the module import_export that can handle .mat files. Also the gnlse.Solution class has two methods to load and save the results to an external file of this type. For these purposes we used the hdf5storage package.
In turn, the module visualization is responsible for the visualization of the results. We chose the popular matplotlib library for the task. We decided to show results of calculations as twosided plots presenting a solution both in both the time and spectral domains. Default ranges for the xaxes correspond to the case of supercontinum generation as in work [4].
The common module is the last standalone script, which in the future will keep the physical constants used in the calculations. Currently, we have placed there the value of speed of light in vacuum in nm/ps.
4 Numerical illustrations
We investigated the toolbox performance on some common examples presented in [1, 4]. We observe qualitative compliance with a relatively small calculation time.
4.1 Propagation of higherorder soliton
An electromagnetic wave scattered in all directions has different spectral components. The dominant band has the same frequency as the incident wave – it is the Rayleigh band (elastic scattering). On the other hand, the components with frequencies shifted by the characteristic frequencies of molecules are called Raman scattering bands (inelastic scattering). These are the bands of reduced and increased frequencies, the Stokes and antiStokes bands, respectively. AntiStokes scattering occurs less frequently due to the lower occupation of states with higher oscillatory energy [1].
Raman scattering can be both spontaneous or stimulated. The latter is dominant in the case of propagation of short pulses in the optical fiber. In optical fiber transmission, the stimulated Raman scattering can induce optical Raman gain. In addition, due to the forced Raman scattering, the spectrum shifts towards longer wavelengths, which is caused by the amplification of the Stokes wave with frequencies inside the spectrum. It is the socalled selffrequency shift phenomenon.
The selffrequency shift phenomenon can be observed in the case of solitons – optical pulses propagating in a dispersive medium without changing their temporal and spectral intensity profiles. A higherorder soliton is a soliton with a higher energy than that of a fundamental soliton by a factor that is the square of an integer (soliton order). The temporal shape of such soliton varies periodically during propagation. The soliton can be also affected by the nonlinear increase of the pulse steepness (selfsteepening phenomenon). Fig. 5 shows the evolution of the spectral and temporal characteristics of a higherorder soliton in three cases:

propagation without selfsteepening and Raman response,

soliton fission with selfsteepening, but no Raman response accounted,

soliton fission with selfsteepening, and Raman response accounted,
4.2 Supercontinuum generation in anomalous dispersion regime
To a large extent, the interaction of an electromagnetic wave with the system of atoms in a medium depends on the optical properties of matter and the frequency of the incident light. In other words, quantities such as permittivity, refractive index, and absorption are functions of the wave frequency. These dependencies are called the dispersions, which in a broader sense also means the splitting of light into monochromatic components. It influences the phase matching when interacting with several waves [1].
If we take a closer look at the light pulse, which cannot be treated as a monochromatic wave, it can be seen that the dispersion additionally changes its parameters. Due to the different velocities of the group components of the pulse frequency spectrum, the pulse’s time profile is distorted. We can distinguish between two dispersion regimes: normal and anomalous. In the case of anomalous dispersion, the group refractive index increases with increasing wavelength – the components with lower frequencies have a lower group velocity.
Dispersion is rather undesirable in telecommunications systems, and the compensation of the anomalous dispersion can be achieved by optical nonlinearity, delaying the wave components at the beginning of the pulse and accelerating those at the end of [1]. In this way, we can try to compensate for the elongation of the pulse. Nonlinearity also causes many other useful effects that have found their practical application. One of them is the nonlinear spectral broadening process, which we call supercontinuum generation. Fig. 6 shows an example of supercontinuum generation in the anomalous dispersion regime at a central wavelength of 835 nm in a 15 centimeter long fiber using three different shapes of input pulse envelopes.
4.3 Dispersive wave generation in anomalous dispersion regime – different Raman responses
The optical pulse propagating in a fiber with anomalous chromatic dispersion forms a soliton pulse accompanied with a temporally spreading background – the dispersive wave. The background radiation is generated in a high frequency band. The dispersive wave is spread by chromatic dispersion, which is not compensated by the fiber nonlinearity (due to low power). Fig. 7 shows an example of dispersive wave generation in anomalous dispersion regime at a central wavelength of 835 nm in a 15 centimeter long photonic crystal fiber using three different models of the Raman response.
5 Summary
In this report, we present the open source gnlsepython package written in the Python language. It is a set of Phyton scripts for solving the Generalized Nonlinear Schrodringer Equation. It is one of the WUSTFOG students projects developed by Fiber Optics Group, WUST.
The complete and uptodate documentation of the project is available at https://gnlse.readthedocs.io.
Authors’ contributions statement
PR, MZ, AP, DS, SM, and KT contributed to the final code. PR is the solution architect, who proposed the structure of gnlse integration module based on the existing MATLAB code. MZ prepared a module implementing the initial pulses’ temporal shapes, and a module for importing and exporting data. AP prepared an implementation of the Raman responses models and the dispersion operator interpolation components. DS took responsibility for implementing the visualization of the results. SM implemented the nonlinear module for accounting the dispersion of effective mode area. All authors were involved in the preparation of example visualization of physical phenomena. KT and SM supervised the project. SM was responsible for integrating the existing code with all separate modules and prepared the technical document, while KT watched over the correct implementation from the theoretical point of view.
References
 [1] Govind P. Agrawal. Nonlinear Fiber Optics, 6th Edition. Academic Press, 2019.
 [2] M. Aleshams, A. Zarifkar, and M.H. Sheikhi. Splitstep fourier transform method in modeling of pulse propagation in dispersive nonlinear optical fibers. In Proceedings of CAOL 2005. Second International Conference on Advanced Optoelectronics and Lasers, 2005., volume 2, pages 124–126 vol. 2, 2005.
 [3] J. Hult. A fourthorder runge–kutta in the interaction picture method for simulating supercontinuum generation in optical fibers. Journal of Lightwave Technology, 25(12):3770–3775, 2007.
 [4] J. M. Dudley and J. R. Taylor. Supercontinuum Generation in Optical Fibers. Cambridge University Press, 2010.
 [5] K.J. Blow and D. Wood. Theoretical description of transient stimulated raman scattering in optical fibers. IEEE Journal of Quantum Electronics, 25:2665–2673, 1989.
 [6] Q. Lin and G. P. Agrawal. Raman response function for silica fibers. Optics Letters, 31(21):3086, 2006.
 [7] Dawn Hollenbeck and Cyrus D. Cantrell. Multiplevibrationalmode model for fiberoptic raman gain spectrum and response function. J. Opt. Soc. Am. B, 19(12):2886–2892, Dec 2002.
 [8] Bertrand Kibler, J. Dudley, and S. Coen. Supercontinuum generation and nonlinear pulse propagation in photonic crystal fiber: Influence of the frequencydependent effective mode area. Applied Physics B: Lasers and Optics, 81:337–342, 07 2005.
 [9] J. Lægsgaard. Mode profile dispersion in the generalized nonlinear schrödinger equation. Opt. Express, 15(24):16110–16123, Nov 2007.
Comments
There are no comments yet.