gnlse-python: Open Source Software to Simulate Nonlinear Light Propagation In Optical Fibers

by   Paweł Redman, et al.

The propagation of pulses in optical fibers is described by the generalized nonlinear Schrodinger equation (GNLSE), which takes into account the fiber losses, nonlinear effects, and higher-order chromatic dispersion. The GNLSE is a partial differential equation, whose order depends on the accounted nonlinear and dispersion effects. We present gnlse-python, a nonlinear optics modeling toolbox that contains a rich set of components and modules to solve the GNLSE using the split-step Fourier transform method (SSFM). The numerical solver is freely available, implemented in Python language, and includes a number of optical fiber analysis tools. Code and data are available at



There are no comments yet.


page 5

page 7

page 8

page 9


Physics-informed Neural Network for Nonlinear Dynamics in Fiber Optics

A physics-informed neural network (PINN) that combines deep learning wit...

pyLLE: a Fast and User Friendly Lugiato-Lefever Equation Solver

We present the development of pyLLE, a freely accessible and cross-platf...

Voxelwise nonlinear regression toolbox for neuroimage analysis: Application to aging and neurodegenerative disease modeling

This paper describes a new neuroimaging analysis toolbox that allows for...

Soliton propagation in lossy optical fibers

In this work, we study the propagation of solitons in lossy optical fibe...

Random problems with R

R (Version 3.5.1 patched) has an issue with its random sampling function...

Fast and accurate waveform modeling of long-haul multi-channel optical fiber transmission using a hybrid model-data driven scheme

The modeling of optical wave propagation in optical fiber is a task of f...

Spiral capacitor calculation using FEniCS

The paper shows how to optimize a water level sensor consisting of a cyl...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Nonlinear optics is a branch of optics dealing with phenomena that are consequences of light-induced 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 self-modulation. It is successfully used in modeling the picosecond propagation of low-energy 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 time-consuming, 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 split-step 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


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 self-steepening 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 



The first and second terms on the right-hand 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


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]


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 Runge-Kutta (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 fourth-order Runge-Kutta 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 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 split-step 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 full-width half-maximum 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 secant-shaped envelopes and a continuous wave with optional noise.

Figure 1: Example input pulse envelops in time domain.

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]


Since the frequencies present in time-limited pulses are concentrated around the central frequency , we can use the expansion of into the Taylor series




For such representation, the dispersion operator in the frequency domain is defined as


where the second term of the right-hand side is responsible for loss.

The and parameters are defined as


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 third-order 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]


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


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


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.

Figure 2: Example of supercontinuum generation using different dispersion operators. In our example we used provided coefficients up to the 10th derivative.

3.3 Nonlinear phenomena: Raman responses

The last term in Eq. 1 describes nonlinear phenomena such as: phase self-modulation, self-steepening, shock wave formation, and Raman scattering. The function, which determines the last of the listed effects, usually is modelled as [5]




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:

  • based on K. J. Blow and D. Wood model [5], presented in Eq. 14,

  • 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 multi-vibrational-mode function,

Their visualisations in the time domain are presented in Fig. 3.

Figure 3: Raman scattering functions for silica optical fibers for three models in time domain.

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 self-steepening 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


This approach is more rigorous than the approximation of () and requires the definition of a pseudo-envelope as


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.

Figure 4: Example of soliton fission using different nonlinear coefficient values. In both cases the dispersion operators given by Taylor expansion of were used.

3.5 Additional tools

We provide some additional tools and functions inside gnlse-python 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 two-sided plots presenting a solution both in both the time and spectral domains. Default ranges for the x-axes correspond to the case of supercontinum generation as in work [4].

The common module is the last stand-alone 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 higher-order 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 anti-Stokes bands, respectively. Anti-Stokes 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 so-called self-frequency shift phenomenon.

The self-frequency 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 higher-order 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 (self-steepening phenomenon). Fig. 5 shows the evolution of the spectral and temporal characteristics of a higher-order soliton in three cases:

  • propagation without self-steepening and Raman response,

  • soliton fission with self-steepening, but no Raman response accounted,

  • soliton fission with self-steepening, and Raman response accounted,

Figure 5: Evolution of the spectral and temporal characteristics of the higher-order N = 3 soliton in three cases.

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.

Figure 6: Example of supercontinuum generation in anomalous dispersion regime.

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.

Figure 7: Example of dispersive wave generation in anomalous dispersion regime.

5 Summary

In this report, we present the open source gnlse-python 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 WUST-FOG students projects developed by Fiber Optics Group, WUST.

The complete and up-to-date documentation of the project is available at

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.


  • [1] Govind P. Agrawal. Nonlinear Fiber Optics, 6th Edition. Academic Press, 2019.
  • [2] M. Aleshams, A. Zarifkar, and M.H. Sheikhi. Split-step 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 fourth-order 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. Multiple-vibrational-mode model for fiber-optic 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 frequency-dependent 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.