Numerical Method for Coupled Nonlinear Schrödinger Equations in Few-Mode Fiber

: This paper discusses novel approaches to the numerical integration of the coupled nonlinear Schrödinger equations system for few-mode wave propagation. The wave propagation assumes the propagation of up to nine modes of light in an optical ﬁber. In this case, the light propagation is described by the non-linear coupled Schrödinger equation system, where propagation of each mode is described by own Schrödinger equation with other modes’ interactions. In this case, the coupled nonlinear Schrödinger equation system (CNSES) solving becomes increasingly complex, because each mode affects the propagation of other modes. The suggested solution is based on the direct numerical integration approach, which is based on a ﬁnite-difference integration scheme. The well-known explicit ﬁnite-difference integration scheme approach fails due to the non-stability of the computing scheme. Owing to this, here we use the combined explicit/implicit ﬁnite-difference integration scheme, which is based on the implicit Crank–Nicolson ﬁnite-difference scheme. It ensures the stability of the computing scheme. Moreover, this approach allows separating the whole equation system on the independent equation system for each wave mode at each integration step. Additionally, the algorithm of numerical solution reﬁning at each step and the integration method with automatic integration step selection are used. The suggested approach has a higher performance (resolution)—up to three times or more in comparison with the split-step Fourier method—since there is no need to produce direct and inverse Fourier transforms at each integration step. The key advantage of the developed approach is the calculation of any number of modes propagated in the ﬁber.


Introduction
This article is a continuation of our previous work [1], where we suggested a novel solution of CNSES for simulation of ultra-short optical pulse propagation in birefringent fibers. The introduction of our previous work [1] completely reveals the reasons, importance, and necessity of the problem formulation. The main motive of our study was because the designers of fiber optic transmitting and receiving devices lack the tools to determine the shape of distortion of a signal that is transmitted at high power or transmitted over very long distances. Even the commonly-used split-step Fourier method, which gives quite satisfactory results, provides a general solution for the task at hand without particularities. The difference between our previous publication [1] and the current is in a number of propagation modes in the fiber. Birefringent fibers propose propagating two wave modes while few-mode fibers propose propagating up to nine modes. Our goal was to propose a method to solve the system of Schrödinger equations, which allows obtaining more detailed information on the shape of distortions of the signal propagating in the few-mode fiber under various conditions.
It is necessary to emphasize that pulse propagation in fibers with abnormal dispersion is an important investigation theme as femtosecond lasers have an influential position in industrial field [2][3][4][5]. Special optical fibers are developed for transmission of high-power ultrashort pulses [6,7] with accent to polarization maintaining fibers. The negative role of chromatic dispersion appears with increase in the length of the optical fiber communication line. The physical reason for dispersion is the frequency dependence of the phase velocity or phase delay during propagation of electromagnetic waves in the medium. This leads to a violation of the phase vectors in the wave package. When these negative effects are studied, only envelope distortion is usually considered. However, phase distortions are also important for certain classes of problems. For example, the problem of chirping during the propagation of an optical signal in an optical fiber with nonlinear dispersion is one of them.
The first part of the work is devoted to the computing of complex envelop calculations, the second is devoted to the phase velocity or phase delay calculation during the propagation of optical waves in the fibers.

CNSES for Few Modes in Dimensionless Form
The evolution of optical waves in a fiber can be described by the CNSES as: with the initial and boundary conditions for each propagation mode: The definitions are used in Equations (1) and (2): W i -the complex envelope of the optical wave of the i-th mode, α (i) -attenuation of the i-th mode; β (i) 1 , β (i) 2 , β (i) 3 -the first, second and third order dispersion of the i-th mode respectively; γ (i) -parameter of nonlinearity for the i-th mode; C i,m , B i,m -coupling coefficients between the i-th and m-th modes; T R -Raman scattering parameter; ω (i) 0 -angular frequency of the i-th mode; z-coordinate along an optical fiber; t-time, j-imaginary one, T N -final time, and f (i) (t)known functions at the initial time.
To transfer Equation (1) into dimensionless form, the specific values of the process are used, namely Llength, T N -time, and P-power, by relations: Replacement of Equation (3) into Equation (1) transforms it into a dimensionless form: where the dimensionless coefficients: and the definitions for the non-linear component: are used.
Equation (1) is specially separated on linear and nonlinear terms, it allows us to realize the numerical integration algorithm based on computing scheme, where all linear terms are written in implicit and all nonlinear terms are written in explicit finite-difference form.

The Finite-Difference Scheme and Computing Scheme
The Equation system (4) can be simplified by using F (i) (x) terms to denote the right parts of Equation (4): The finite-difference equivalences for Equation (7), according to [1], have the form: where bottom index k is used to denote the mesh points along the length ξ k = ∆ξ·(k − 1) and n is used to denote the dimensionless time mesh points τ n = ∆τ·(n − 1). The dimensionless parameter θ allows to attribute the value of right part of Equation (8) to an arbitrary mesh point with fractional indexes (k + θ, n − 1 /2). The parameter θ ∈ [0, 1], defines the explicit θ = 0 or implicit θ ∈ (0, 1] finite-difference computing scheme. The recommendation for θ parameter choosing was given in [1]. There, it was approved that θ = 1 /2, which leads to stable and reasonable results. All variables with the k-index in Equation (8) are known as well as all variables with the (k + 1) indexes are known at n equal to 0, 1, 2, and N from Equation (2). The values of F (i) are attributed to middleware virtual mesh points between (k+ θ, n − 1 /2).
The nonlinearity in Φ1 (i) and Φ2 (i) , and a third-order partial derivative by time in Equation (4), require modification of Crank-Nicolson's method [1,8]. The aim is to use implicit form for linear terms, and explicit form for nonlinear ones. Thus, we separate the F (i) on the sum of linear and nonlinear terms: with the definitions for linear L (i) and non-linear N (i) terms at the k-layer: For the calculation of the linear terms L (i) on the (k + 1)-th layer, we can just replace k on (k + 1) in Equation (10). For the nonlinear terms on (k + 1)-th layer, we must use finite-difference definition in explicit form: Hence, the explicit form in Equation (8) is used only for Φ1 (i) and Φ2 (i) , while the implicit form is used for the linear terms. The nonlinear terms Φ1 (i) and Φ2 (i) in virtual grid point (k, n − 1 2 ) are written as: The desired function values in (k, n − 1 2 ) point are written as an average of two adjacent values. The partial derivative from this functions in the mentioned point can be written as central finite-differences with second-order accuracy. It gives the linear equation system for each mode at each integration step: It is very important to note that the linear Equation system (13) for the i-th mode depends only on x (i) k+1,n unknowns and does not depend on unknowns from other wave modes. It allows solving the linear equation system for the each i-th mode independently. The index n for each linear equation system for the i-th mode begins from 2 up to the N − 1, because the values of x k+1,1 , x k+1,2 , and x k+1,N are known due to Equation (2).
It should be highlighted that the Equation system (13) decomposes into M independent line equation systems, according to grid node variables x (i) k+1,n (where n begins from 3 and increases to N − 1). Each equation system can be solved separately. These linear equation systems have four-diagonal matrix form, where in addition to the main diagonal, there are one "upper" and two "sub" diagonals. We use the modified Thomas tridiagonal algorithm, which allows us to transform this matrix to triangular form. The described refinement algorithm is used at each integration step. The definition of nonlinear terms (from previous k-th integration layer) explicitly contributes to the error in calculating the values of the new (k + 1)-th layer. The main goal of the proposed refinement algorithm is to introduce iteration of the computation process at each integration step, which corrects the explicit form of nonlinear terms [1].

The Ultra-Short Pulse Evolution in Fiber
The ultra-short pulse propagation in optical fiber including third-order dispersion and Raman scattering is described by complete CNSES. The initial data for the calculations were taken from [9][10][11][12][13][14][15][16]: The single chirped Gauss pulse (chirp C = −0.4579) was put into the fiber−s input; pulse duration is 12 fs, peak power is P = 1.75 × 10 5 W. The pulse form is described as: 20,000 grid points along dimensionless time was chosen; approximately 720,000 integration with initial step ∆ξ = 1·10 −4 d.u. were made. Usage of the automatic integration step correction algorithm was included, allowing to calculate the pulse evolution length up to~2.5 mm. The maximum error was set as 10 −30 d.u. All calculations were made using a processor with double precision and 64-bit architecture.
A comparison of the pulses' forms calculated by our method (red line) and experiment data from [9][10][11][12][13][14][15][16] is shown in Figure 1. 20,000 grid points along dimensionless time was chosen; approximately 720,000 integration with initial step ξ = 1·10 -4 d.u. were made. Usage of the automatic integration step correction algorithm was included, allowing to calculate the pulse evolution length up to ~2.5 mm. The maximum error was set as 10 -30 d.u. All calculations were made using a processor with double precision and 64-bit architecture.
The numerical result (red line in Figure 1) fits well with the experimental result (blue line in Figure 1). It shows that the proposed method is effective and can be used for solving similar nonlinear tasks.

The Phase Velocity or Phase Delay Calculation during the Wave Propagation
Recently, the channel model, based on equivalence principle [17][18][19], was used for description of propagation of wave packages in media. According to this model, the task of wave package propagation was replaced by the task of wave package passing through line system with frequency H(jω) and h(τ) impulse characteristics: where (ω) is phase incursion in media or phase-frequency characteristic of the line system, where the frequency set of the package ω can be presented as the sum of average  and difference Ω package frequencies as ω =  + Ω. The F[] and F -1 [] denote direct and inverse Fourier transform operator.
It is convenient to use the expression of time delays instead of velocities of propagation while studying the dispersion. Therefore, the phase delay τp and group delay τg (fast time) terms are involved into consideration, which are defined from phase-frequency characteristics of the system [20]: The minus sign "-" ahead of the phase (ω, L) of frequency characteristic in Equation (15) is introduced because the derivative of phase defines the phase delay. Group The numerical result (red line in Figure 1) fits well with the experimental result (blue line in Figure 1). It shows that the proposed method is effective and can be used for solving similar nonlinear tasks.

The Phase Velocity or Phase Delay Calculation during the Wave Propagation
Recently, the channel model, based on equivalence principle [17][18][19], was used for description of propagation of wave packages in media. According to this model, the task of wave package propagation was replaced by the task of wave package passing through line system with frequency H(jω) and h(τ) impulse characteristics: where ϕ(ω) is phase incursion in media or phase-frequency characteristic of the line system, where the frequency set of the package ω can be presented as the sum of average and difference Ω package frequencies as ω = + Ω. The F[] and F −1 [] denote direct and inverse Fourier transform operator. It is convenient to use the expression of time delays instead of velocities of propagation while studying the dispersion. Therefore, the phase delay τ p and group delay τ g (fast time) terms are involved into consideration, which are defined from phase-frequency characteristics of the system [20]: The minus sign "-" ahead of the phase ϕ(ω, L) of frequency characteristic in Equation (15) is introduced because the derivative of phase defines the phase delay. Group delay is a differential characteristic of phase incursion in the media (in the neighborhood of the selected frequency), and phase delay is an integrated characteristic.
The most common approach in the investigation of dispersion is related to the differential characteristic definition, where frequency characteristic covers the frequency range Ω ch with average frequency , and amplitude-frequency and phase-frequency characteristics (Equation (15)) are presented by Taylor power series of deferential frequencies (ω − ) = Ω at the average frequency point: For a homogeneous medium, among which there is optical fiber, the phase-frequency characteristic is related to the refractive index of a fiber n(ω) and fiber length L by the relation: The coefficients at frequencies in the Taylor power series of phase-frequency characteristic Equation (18) are the parameters of phase dispersion of appropriate order. Thus, the dispersion parameters evolve by increasing the length L in dispersal systems. When some critical length is achieved, dispersal distortions occur due to non-linear term effects on phase-frequency characteristic.
It is convenient to use the group phase delay as a function of frequency for phase dispersion analysis in the channel model: The ϕ ( , L) = τ g ( , L) parameter is commonly referred to as Group-Delay Dispersion (GDD), and the τ" g ( , L) = ϕ ( , L) = TOD parameter is commonly referred to as Third-Oder Dispersion (TOD).
The complex frequency characteristic is proportional to the non-zero function in the frequency channel in case of phase incursion described by Equation (18) without third-order dispersion (ϕ ( ,L) = 0): The chirplet with rectangular duration window T ef in the time domain corresponds to the function exp(-jϕ ( , L)Ω 2 /2) in the frequency domain at Ω ch T ef >> 1 (where T ef ≈ ϕ ( ,L)·Ω ch ), according to Fourier transform property: The chirplets with a rectangular duration window are shown in Figure 1 when Ω ch T ef >> 1. The low frequencies are shown by red lines, and high frequencies are shown by blue lines.
The picture in Figure 2 shows that second-order frequency dispersion leads to chirping of impulse, characterised by fast time. In the case of normal dispersion, the low frequencies outrun high frequencies and the chirp sign is positive (υ( , L) > 0), in the case of anomalous dispersion, the chirp sign is negative (υ( , L) < 0).
The rough equilibrium in Equation (23) is true at ΩchTef >> 1 and other than zero at delays τ ∈ [τg(ϖ, L) − Tef/2; τg(ϖ, L) + Tef/2). According to Equation (23), the pulse characteristic of the channel is a chirplet by fast time τ, that is second-order dispersion leads to chirping of impulse characteristic. The rate of the frequency changing in the chirplet is positive for normal dispersion and negative for anomalous dispersion.
The chirp effect can lead to the compactification of linear modulated pulse in the medium with second-order dispersion. Optical pulse with Gauss form with line frequency modulation ω = ϖ + ω′t, as well as frequency characteristic of communication line has a phase with quadratic dependence from frequency ϕ″(ϖ, L)Ω 2 /2. Therefore, if the condition Ω 2 /2ω′ + ϕ″(ϖ, La)•Ω 2 /2 = 0 is fulfilled at given communication line length L = La, the linear frequency modulated pulse is matched with dispersal radio channel and it compresses over the fast time. This condition can be rewritten in the form ω′•ϕ″(ϖ, La) = -1. It means that the submitted derivatives must be reciprocal in magnitude and opposite in signs. Thus, for fiber with normal dispersion ϕ″(ϖ, L) > 0, the frequency changing rate of the optical pulse must be negative ω′ < 0, while for the fiber with anomalous dispersion, it must have a positive sign.
With the increase of line length at different signs of derivatives, the dispersion parameter ϕ″(ϖ, L) increases by the module. Therefore, the correlation with line frequency-modulated pulse with dispersion in optical fiber is fulfilled for certain line length L = La.
The computing experiments for the effect were carried out, considering that the phase-frequency characteristic of the channel, given as the line system, is dimensionless. Therefore, the terms in Equation (18) are also dimensionless; it allows to present each term as the product of dimensionless values.
Let us get the coherence range of the channel Ωc2 with second-order dispersion from the equivalence |ϕ″(ϖ, L)(Ωc2/2) 2 | = 1, choosing as the line length L = Le2 the value, for which Ωc2 = Ωch. Therefore, the communication line length can be evaluated by the relation: According to Equation (23), the pulse characteristic of the channel is a chirplet by fast time τ, that is second-order dispersion leads to chirping of impulse characteristic. The rate of the frequency changing in the chirplet is positive for normal dispersion and negative for anomalous dispersion.
The chirp effect can lead to the compactification of linear modulated pulse in the medium with second-order dispersion. Optical pulse with Gauss form with line frequency modulation ω = + ω t, as well as frequency characteristic of communication line has a phase with quadratic dependence from frequency ϕ ( , L)Ω 2 /2. Therefore, if the condition Ω 2 /2ω + ϕ ( , L a )·Ω 2 /2 = 0 is fulfilled at given communication line length L = L a , the linear frequency modulated pulse is matched with dispersal radio channel and it compresses over the fast time. This condition can be rewritten in the form ω ·ϕ ( , L a ) = −1. It means that the submitted derivatives must be reciprocal in magnitude and opposite in signs. Thus, for fiber with normal dispersion ϕ ( , L) > 0, the frequency changing rate of the optical pulse must be negative ω < 0, while for the fiber with anomalous dispersion, it must have a positive sign.
With the increase of line length at different signs of derivatives, the dispersion parameter ϕ"( , L) increases by the module. Therefore, the correlation with line frequencymodulated pulse with dispersion in optical fiber is fulfilled for certain line length L = L a .
The computing experiments for the effect were carried out, considering that the phase-frequency characteristic of the channel, given as the line system, is dimensionless. Therefore, the terms in Equation (18) are also dimensionless; it allows to present each term as the product of dimensionless values.
Let us get the coherence range of the channel Ω c2 with second-order dispersion from the equivalence |ϕ ( , L)(Ω c2 /2) 2 | = 1, choosing as the line length L = L e2 the value, for which Ω c2 = Ω ch . Therefore, the communication line length can be evaluated by the relation: Thus, it gives the equations for dimensionless values (dimensionless length m = L/L e2 , dimensionless second-order dispersion coefficient, and dimensionless frequency η = 2Ω/Ω ch ), where second-order dispersion is described by the following relation: The computing results of the compression of line frequency-modulated pulse with the Gauss form in the fiber with second-order dispersion in dimensionless values are presented in Figure 3. The real part of its envelope is given in the form, where its duration defines the frequency range of the channel: Fibers 2021, 9, x FOR PEER REVIEW 8 of 11 Thus, it gives the equations for dimensionless values (dimensionless length m = L/Le2, dimensionless second-order dispersion coefficient, and dimensionless frequency η = 2Ω/Ωch), where second-order dispersion is described by the following relation: (26) The computing results of the compression of line frequency-modulated pulse with the Gauss form in the fiber with second-order dispersion in dimensionless values are presented in Figure 3. The real part of its envelope is given in the form, where its duration defines the frequency range of the channel: The increase in amplitude and increase in optical loss due to second-order dispersion have similar values (the curve is symmetrical about the maximum). To prevent second-order dispersion, methods of pre-distortion of a light pulse during transmission or correction of dispersion during reception are used. This effect illustrates a technique for overcoming second-order dispersion by transmitting an optimal chirp pulse to the channel.
With an increase in the length of the optical fiber, the third-order dispersion can be of fundamental importance. It can be shown that in this case, the analytical solution for the impulse response has the form: where: The increase in amplitude and increase in optical loss due to second-order dispersion have similar values (the curve is symmetrical about the maximum). To prevent second-order dispersion, methods of pre-distortion of a light pulse during transmission or correction of dispersion during reception are used. This effect illustrates a technique for overcoming second-order dispersion by transmitting an optimal chirp pulse to the channel.
With an increase in the length of the optical fiber, the third-order dispersion can be of fundamental importance. It can be shown that in this case, the analytical solution for the impulse response has the form: where: are incomplete Airy integrals.
In the dimensionless variables, the line length was chosen under conditions that a channel with third-order dispersion coherence bandwidth is equal to its frequency bandwidth Ω C3 = Ω ch : L e3 = 6c/|( n( )) |(Ω ch /2) 3 . In this case, the dimensionless communication line length is m = L/L e3 , the dimensionless dispersion coefficient is p e3 = Ω ch /Ω C3 = 1, and dimensionless frequency is η = 2Ω/Ω ch ∈ [−1, 1]. The third-order nonlinear component at given dimensionless values for phase equals: It must be noted that dimensionless length differs from the length, which was taken for second-order dispersion.
The results of calculating the Gaussian pulse distortion for TOD > 0 and TOD < 0 are shown in Figure 4a,b, respectively. In both cases, starting from a certain length of the optical path, the effect of pulse collapse is observed when it disintegrates into many short pulses. The number of additional pulses depends on the path length (the value of the TOD parameter). At TOD > 0, the collapse appears behind the main pulse and at TOD < 0-in front of it. Obviously, the collapse is described by Airy integrals (28). are incomplete Airy integrals.
In the dimensionless variables, the line length was chosen under conditions that a channel with third-order dispersion coherence bandwidth is equal to its frequency bandwidth ΩC3 = Ωch : Le3 = 6c/|(ϖn(ϖ))″′|(Ωch/2) 3 . In this case, the dimensionless communication line length is m = L/Le3, the dimensionless dispersion coefficient is pe3 = Ωch/ΩC3 = 1, and dimensionless frequency is η = 2Ω/Ωch ∈[-1, 1]. The third-order nonlinear component at given dimensionless values for phase equals: It must be noted that dimensionless length differs from the length, which was taken for second-order dispersion.
The results of calculating the Gaussian pulse distortion for TOD > 0 and TOD < 0 are shown in Figure 4a,b, respectively. In both cases, starting from a certain length of the optical path, the effect of pulse collapse is observed when it disintegrates into many short pulses. The number of additional pulses depends on the path length (the value of the TOD parameter). At TOD > 0, the collapse appears behind the main pulse and at TOD < 0-in front of it. Obviously, the collapse is described by Airy integrals (28).  In the case of second-and third-order dispersion, as the line length increases, in the beginning, pulse broadening is observed, then its asymmetry appears, and then collapse occurs.

Conclusions
This paper shows the effectiveness of the proposed method in solving the CNSES for modeling the propagation of pulses formed by a group of highly coupled modes. Developing the results obtained in [1], it is shown that the decomposition of the CNSES at each integration step into a set of independent systems of linear equations makes it possible to include an arbitrary number of propagation modes in the equations and study their mutual influence. A positive comparison of the results obtained using the proposed algorithm with the works of other authors has shown its viability and prospects in order to further improve its characteristics.
Within the framework of the channel model, an analytical solution is obtained. It describes the propagation of an optical pulse in an optical fiber with a varying length. It is shown that the second-order frequency dispersion can lead to the effect of compression of the optical chirp pulse if the dispersion parameter GDD and the rate of change of the frequency ω′ differ in sign, and if their modules are inversely proportional, then the In the case of second-and third-order dispersion, as the line length increases, in the beginning, pulse broadening is observed, then its asymmetry appears, and then collapse occurs.

Conclusions
This paper shows the effectiveness of the proposed method in solving the CNSES for modeling the propagation of pulses formed by a group of highly coupled modes. Developing the results obtained in [1], it is shown that the decomposition of the CNSES at each integration step into a set of independent systems of linear equations makes it possible to include an arbitrary number of propagation modes in the equations and study their mutual influence. A positive comparison of the results obtained using the proposed algorithm with the works of other authors has shown its viability and prospects in order to further improve its characteristics.
Within the framework of the channel model, an analytical solution is obtained. It describes the propagation of an optical pulse in an optical fiber with a varying length. It is shown that the second-order frequency dispersion can lead to the effect of compression of the optical chirp pulse if the dispersion parameter GDD and the rate of change of the frequency ω differ in sign, and if their modules are inversely proportional, then the maximum compression takes place, i.e., such a pulse is optimal for a fiber of a given length. This effect can be used for optical pulse pre-distortion. The third-order dispersion increasing with length at the beginning leads to an asymmetry of a pulse with a Gaussian envelope, and then to a collapsing effect when it decays into many short adjacent pulses. At TOD > 0, the collapse appears behind the main pulse, and at TOD < 0, in front of it. In the case of second and third-order dispersion, as the cable length increases, in the beginning, pulse broadening is observed, then its asymmetry appears, and then collapse occurs.