Efficient algorithms for solving the spectral scattering problems\\ for the Manakov system of nonlinear Schroedinger equations
EEfficient algorithms for solving the spectral scattering problemsfor the Manakov system of nonlinear Schroedinger equations
L. L. Frumin , ∗ Institute of Automation and Electrometry SB RAS, Novosibirsk 630090, Russian Federation, and Novosibirsk State University, Novosibirsk 630090, Russian Federation (Dated: June 9, 2020)“Vectorial” numerical algorithms are proposed for solving the inverse and direct spectral scatteringproblems for the nonlinear vector Schroedinger equation, taking into account wave polarization,known as the Manakov system. It is shown that a new algebraic group of 4-block matrices withoff-diagonal blocks consisting of special vector-like matrices makes possible the generalization ofnumerical algorithms of the scalar problem to the vector case, both for the focusing and defocusingManakov systems. As in the scalar case, the solution of the inverse scattering problem consistsof inversion of matrices of the discretized system of Gelfand-Levitan-Marchenko integral equationsusing the Toeplitz Inner Bordering algorithm of Levinson’s type. Also similar to the scalar case, thealgorithm for solving the direct scattering problem obtained by inversion of steps of the algorithmfor the inverse scattering problem. Testing of the vector algorithms performed by comparing theresults of the calculations with the known exact analytical solution (the Manakov vector soliton)confirmed the numerical efficiency of the vector algorithms.
I. INTRODUCTION
The nonlinear Schroedinger equation (NLSE) widelyused in modern science and technology as one of themost fundamental mathematical models. NLSE belongsto the nontrivial class of integrable nonlinear partial dif-ferential equations whose solutions can be found by theInverse Scattering Transform method (IST) [1–3]. Thescalar NLSE appears abundantly in theoretical physicsand nonlinear physical optics. It can also use to de-scribe the propagation of information signals throughfiber-optical communication lines [4, 5]. Careful con-sideration of polarization phenomena in a medium withdispersion and Kerr nonlinearity is of paramount impor-tance for the development of modern nonlinear physicsand optics.Manakov [6], when exploring self-focusing of lightbeams and self-induced transparency phenomena witha non-negligent contribution of polarization in nonlin-ear dispersive optical media, was probably the firstwho introduced the vector variant of the NLSE, nowknown as the Manakov System. This system consistsof two interaction-coupled nonlinear Schroedinger equa-tions for two optical polarizations. We will not presenthere the equations of the vector NLSE since the fur-ther consideration is based exclusively on the Gelfand-Levitan-Marchenko system of coupled integral equations(GLME), that applied to the solution of the spectral scat-tering problems in the vector case.Manakov showed that his vector variant of the NLSEbelongs to the class of integrable systems. He constructedthe corresponding L-A Lax operator pair, and using themethod of the IST described the general N-soliton solu-tions, and also found its particular solution, known as ∗ Electronic address: [email protected] the Manakov vector soliton. Subsequently Basharov andMaimistov [7] (see also [8]) discovered that the ManakovSystem could be used to describe many other nonlinearpolarization optical effects, including the propagation ofultrashort polarized optical pulses in a resonant two-levelenvironment.To solve the direct and inverse scattering problemsfor the scalar Schroedinger equation in the frame of ISTnumerical algorithms of Toeplitz Inner Bordering (TIB)[9, 10] have been developed. They based on the direct nu-merical solution of GLM integral equations. The TIB al-gorithm of the inverse scattering problem is efficient, i.e.,it is fast, accurate, and stable because it is a modificationof the well known Levinson algorithm [11]. The TIB algo-rithm for solving the direct scattering problem obtainedby inversion of steps of the algorithm for the inverse prob-lem. The numerical efficiency of the algorithms caused bythe Toeplitz symmetry of the discretized GLME system.TIB algorithms find applications in various optics prob-lems, including Bragg gratings synthesis [9, 12, 13], anddevelopment of new nonlinear approaches to the trans-mission of information in fiber-optic lines [14–17].The aim of this paper is a generalization of TIB al-gorithms for the solution of scattering problems for theManakov system of vector NLSE.In the next Section II, we consider the GLM integralequations and give a short description of its applicationto the solution of the vector NLSE. Section III describesa replacement of variables and discretization of GLME.In Section IV, we introduce vector-like matrices and alsoa new algebraic group of 4-block matrices. On the baseof these constructions, in Section V, we derive a vectorTIB algorithm for solving the inverse scattering problem.In SectionVI, the schematics of the vector algorithms forsolving inverse and direct scattering problems are pre-sented. Section VII contains some results of numericalsimulation and testing of the vector TIB algorithms, andSection VIII is a Conclusion. a r X i v : . [ n li n . S I] J un II. VECTOR GLME
The spectral scattering problems for the vector NLSEreduced in [7] and monograph [18]) to a system of nineintegral GLME, in the same way as Zakharov and Shabatdid it for the scalar NLSE in the famous work [3]. How-ever, this cumbersome system can split into three inde-pendent groups, each of three integral equations. It turnsout that it is enough to consider the only one group ofthree integral equations for spectral scattering problemsfor the Manakov system. In the dimensionless notationclose to the notation of Lam’s monograph [19] the systemof GLME for the left-hand scattering problem consists ofthe next equations:A ∗ ( x, y ) + x (cid:90) −∞ (cid:88) α =1 , A α ( x, z )Ω α ( y + z ) dz = 0 , (1) ± A ∗ α ( x, y )+ x (cid:90) −∞ A ( x, z )Ω α ( y + z ) dz = − Ω α ( x + y ) . (2)Here − x < y, z < x , and Ω , Ω are GLME kernels. Hereand hereinafter α = 1 ,
2. An asterisk means complexconjugate. The top sign of the symbol ± corresponds tothe Manakov defocusing system and the lower sign to thefocusing one. The difference from Lam’s notation is thatthe indices of the functions A , A , A do not begin fromunity, but from zero.The solution of the inverse scattering problem is twocomponents of “potential” vector function for two orthog-onal polarizations: q ( x ) , q ( x ), connected with solutionsof GLME equations by the synthesizing relations:q α ( x ) = ± ∗ α ( x, x − . (3)Two critical remarks should make here. First, equa-tions 1 and 2 present actually to not two-dimensional,but three-dimensional problem, since all the consideredfunctions A , A , A and the components of the vectorkernel Ω α parametrically depend on an additional evo-lutionary variable. Since this dependence arises explic-itly only when considering the evolution of the solutionin time or along the optical line, it usually omitted forbrevity’s sake. In this paper, such an evolutionary vari-able corresponds to time-variable that we outline here-under as t .Secondly, for the focusing case, the components ofthe GLME vector kernel are the sum of the componentsV α ( x, t ) of vector kernel of the continuous spectrum (thevector of pulse response function), and the sum of com-ponents Λ kα ( x, t ) of vector kernel for discrete eigenvaluesof the Lax operator spectrum corresponding to the set ofdiscrete eigenvalues { λ k } of the operator, where the in-dex k numbers the eigenvalues of the discrete spectrum:Ω α ( x, t ) = V α ( x, t ) + (cid:80) k Λ kα ( x, t ) . As in the scalar case, the solution of the Cauchy prob-lem for vector NLSE by IST method consists of a se-quence of three main steps: 1. The direct spectral scattering problem. The com-ponents of known solution vector function q α ( x, t = 0, are used to solve the direct scatteringproblem for the Manakov system, and the scatter-ing data at t = 0 are restored as components of thekernel vector Ω α ( x, α ( x, t = 0,are transformed using the spectral evolution trans-form (see, for example, [1, 3, 19]) into the compo-nents of the vector kernel GLME at time t = T :Ω α ( x, T ). Specifically, for kernel components ofdiscrete spectrum Λ kα , corresponding to the k -theigenvalue λ k , this transformation has the followingform: Λ kα ( x, T ) = Λ kα ( x,
0) exp (cid:8) ( − λ k T ) (cid:9) , where iis imaginary unit.3. The inverse scattering problem. The scatteringdata in the form of vector kernel componentsΩ α ( x, T ), are used to solve the GLM equations andto determine the unknown components of the po-tential vector function q α ( x, T ).To construct a numerical algorithm, the scatteringproblem is considered on a finite interval 0 ≤ x ≤ L :it is assumed that the kernels Ω ( x ) , Ω ( x ) vanish out-side this interval. In this case, the GLM equations takethe following form:A ∗ ( x, y ) + x (cid:90) − y (cid:88) α =1 , A α ( x, z )Ω α ( y + z ) dz = 0 , (4) ± A ∗ α ( x, y ) + x (cid:90) − y A ( x, z )Ω α ( y + z ) dz = − Ω α ( x + y ) , (5)where − x < y, z < x ≤ L , and α = 1 , ( x ) = 0, then the problem becomes scalarfor the potential q = q ( x ), and the reduced system ofGLME should coincide with that for the scalar case. Ifwe put Ω ( x ) = 0, the scalar system of GLME holdssimilarly for the potential q ( x ), and that is a check forthe correctness of our system of GLM equations. III. GLME DISCRETE APPROXIMATION
The first step in the GLME discretization is a replace-ment of variables in GLME. It makes it possible to obtainintegral equations with different arguments of the kernelsthat give matrix blocks with Toeplitz symmetry after thediscretization of these equations.Following [9, 10] we carry out the complex conjugationof equation (4) and replace the variables: z → τ − x, y → x − σ , with 0 ≤ σ, τ < x ≤ L . In equations (5), wesimilarly replace y → τ − x, z → x − σ . We make alsothe replacement of unknown functions:u( x, σ ) = A ( x, x − σ ) , v α ( x, τ ) = ± A ∗ α ( x, τ − x ) . (6)Given these notations, we rewrite (4) and (5) in the form:u( x, σ ) ± x (cid:90) σ (cid:88) α =1 , v α ( x, τ )Ω ∗ α ( τ − σ ) dτ = 0 . (7)v α ( x, τ ) + τ (cid:90) u( x, σ )Ω α ( τ − σ ) dσ = − Ω α ( τ ) . (8)Also we respectively rewrite synthesizing relations (3) asq α ( x ) = v α ( x, x − . (9)We will discretize equations (7) and (8) with the 1storder of approximation accuracy. Let us introduce a dis-crete computational grid: x m = mh/ m = 0 , ..., N ; h = 2 L/N ; (10) σ k = kh ; τ n = nh ; k, n = 0 , ..., m. We replace the integral x (cid:82) σ v α ( x, τ )Ω ∗ α ( τ − σ ) dτ in theequation (7) by the left Riemann sum: m − (cid:88) n = k v α ( x m , τ n ) h Ω ∗ α ( τ n − σ k ) = m − (cid:88) n = k v ( m ) α ; n h Ω ∗ α ; n − k . Here italic letters denotes the grid vectors: v ( m ) α ; n =v α ( x m , τ n ) and Ω ∗ α ; n − k = Ω ∗ α ( τ n − σ k ). Since n ≥ k , thissum is the multiplication of the upper triangular Toeplitzmatrix Q α (hereinafter, we will denote ordinary matricesby capital italics) with elements Q α ; k,n = h Ω ∗ α ; n − k , withsize m × m , by the vector v ( m ) α ; n , with size m . The super-script ( m ) ascribed to vector components here indicatesthat it corresponds to m th step of the algorithm. Forbrevity sake we do not use this superscript for matrices.The discrete analog of equation (7) now can be repre-sented as: u ( m ) k ± m − (cid:88) n = k (cid:88) α =1 , Q α ; k,n v ( m ) α ; n = 0 , (11)where k = 0 , ..., m −
1. Case m = 0 corresponds to theinitial condition: u (0)0 = 0The integrals τ (cid:82) u( x, σ )Ω α ( τ − σ ) dσ in equations (8) canbe represented using the left Riemann sum as products oflow triangular Toeplitz matrices R α with size m × m andwith elements R α ; n,k = h Ω α ; n − k , where n ≥ k , on thevector u ( m ) k = u( x m , σ k ). Note here that matrices Q α are Hermitian conjugations of matrices R α : Q α = R † α . We write the discrete analog of equations (8) in the nextform: v ( m ) α ; n + n − (cid:88) k =0 R α ; n,k u ( m ) k = r α,n , (12)where right part is r α,n = − Ω α,n , and n = 0 , ..., m − m = 0 in the form v (0) α ;0 = − Ω α, again correspondsto the initial condition for the potential vector: q α ;0 =2 v (0) α ;0 = − Ω α, With m changing from 1 to N we get N systems oflinear equations (11)–(12) with size of 3 m × m . Thesesystems are “nested” one into another that resembles abordering numerical algorithm. For the numerical solu-tion of the inverse scattering problem, it is necessary tosolve all the obtained nested systems and determine thecomponents of the potential vector: q α ; m = 2 v ( m ) α ; m , m = 1 , ..., N. (13)The direct numerical solution of the nested systemsof equations (11)–(12) by the Gauss elimination methodrequires O (cid:0) (3 N ) (cid:1) floating-point operations, and for ac-tual problem sizes when N can reach several thousand, itis possible only with supercomputers or computer clus-ters. The best variant of the algorithm for solving sucha series of nested linear systems seems to be a Levinson-type bordering algorithm [11], which in the process ofthis bordering addresses all the systems and required only O ( N ) floating-point operations.Recall that in the case of a scalar NLSE, the system ofGLME consists of two coupled integral equations. Thediscrete form of GLME for the scalar case with first-orderapproximation accuracy is derived from equations (11)and (12) if we put R = 0 and v = 0, and also omit thelower indices of the matrix R and grid vector v ( m )1 . Thediscrete form of GLME for the scalar case has the matrixform of nested systems of linear equations: (cid:18) E ± R † R E (cid:19) (cid:18) u ( m ) v ( m ) (cid:19) = (cid:18) r ( m ) (cid:19) , (14)where, m = 1 , .., N , and the unknown column vector ofsize 2 m is composed of two concatenated column vectors u ( m ) , v ( m ) , each of size m . Block E is the identity matrix, R the lower triangular matrix and Hermitian conjugate R † is upper triangular Toeplitz matrix, all these of a sizeof m × m . The in the right-hand side of Eq. (14) denotesa zero column vector with size m , and column vector r ( m ) is given by the vector of discrete samples of the GLMEkernel: r ( m ) = − ( Ω , ..., Ω m − ) T . (15)The matrix of the system (14) in the scalar case hasthe form of a four-block matrix with Toeplitz symmetry,enabling us to apply the Levinson-type TIB algorithmfor its solution.In the case of the vector NLSE, the discretized GLMequations (11) and (12) consists of three equations. Itleads to a 3 m × m block matrix consisting of nineToeplitz blocks: E ± R † ± R † R E R E , where is zero matrix block with size m × m . Despitethat all these blocks are Toeplitz, the complete matrix ofthe system, unlike the scalar case, does not have Toeplitzsymmetry.Note that if we rewrite vector GLME in the vectorform, we get only two coupled integral equations, onescalar, and the other one vector-like. The main idea ofthis paper is to use the vector notation for reducing thevector case to the scalar one by presenting the discretizedGLME system in the form of not nine, but four blocks, asin the scalar case, but some of them are vector-like blocks.Hereunder we introduce new mathematical constructionsof the vector-like matrices and also 4-block matrices withoff-diagonal vector-like matrices’ that has the algebraicgroup properties and can be used to construct an effi-cient ”vector” algorithm for solving the inverse scatter-ing problem for the vector NLSE, similar to the TIB al-gorithm for the scalar case. IV. 2-MATRICES AND 4-BLOCK MATRICES
Consider a 2-dimensional linear vector space over thefield of complex numbers C, with a scalar product andtwo orthogonal unit vectors (cid:126) e and (cid:126) e . Elements of thisvector space will be called hereinafter c-vectors. We willdenote c-vectors using arrow for it and Roman letters forits components: (cid:126) c = c (cid:126) e + c (cid:126) e .Let us introduce the space of special vector-like matri-ces consisting of pairs of square m × m matrices B , B ,that for brevity sake will be called 2-matrices: B = B (cid:126) e + B (cid:126) e . We will denote the 2-matrices in bold ital-ics. We call the matrices B and B the projections ofthe 2-matrix B onto the unit vectors (cid:126) e and (cid:126) e . Notethat the i, j -th element of the 2-matrix B has a form B i,j (cid:126) e + B i,j (cid:126) e . The 2-matrix B is scalar-wise mul-tiplied by the two-dimensional c-vector (cid:126)b = b (cid:126) e + b (cid:126) e (here b , b are complex numbers), the result is the or-dinary matrix b B + b B .Let us define the left multiplication operation of ma-trix A by 2-matrix B : A B = AB (cid:126) e + AB (cid:126) e , as wellas the right multiplication of 2-matrix B on the matrix A : B A = B A(cid:126) e + B A(cid:126) e . The result in both caseswill be a 2-matrix. Since matrix multiplication is not acommutative operation, the results may be different.Now we consider the scalar product of two 2-matrices A = A (cid:126) e + A (cid:126) e and B = B (cid:126) e + B (cid:126) e . The result is theusual (ordinary) matrix: C = A · B = A B + A B (thedot · hereinafter denotes a scalar product). This scalar product can also be non-commutative. However, scalarmultiplication is associative: A C · B = A · C B , wherethe C matrix is ordinary.We also introduce 2-vectors, the two projections ofwhich are normal m sized vectors, rows, or columns,which we denote in bold italics. The projections, gridvectors, we have already indicated above with simple ital-ics. The fact that the indices i, k, j , and the size m alsoindicated in plain italics should not lead to confusion. Forexample, we consider a column 2-vector b = b (cid:126) e + b (cid:126) e with projections b , b .Note that the k th component of 2-vector is a c-vectorand we will denote it as: (cid:126) b k = b k (cid:126) e + b k (cid:126) e The ordinary matrix A can be multiplied on the left bythe column vector b , and we get the matrix form for twosystems of m linear equations: A b = c : Ab = c , Ab = c , where 2-vectpr c = c (cid:126) e + c (cid:126) e .We can also multiply scalar-wise the 2-matrix B onthe column 2-vector b : B · b = ( B (cid:126) e + B (cid:126) e ) · ( b (cid:126) e + b (cid:126) e ) = B b + B b . The corresponding linear equation B · b = d , where d is an (ordinary) vector of size m , can beinterpreted, for example, as a matrix notation of the sumof two systems of m linear equations. Left multiplicationof the 2-matrix B by the (ordinary) column vector c givesus, as the result, 2-vector f : B c = f , ( f = f (cid:126) e + f (cid:126) e ),that should be interpreted (without discussion of theircompatibility) as a compact vector representation of twolinear systems: B c = f , B c = f . Similarly, one candefine the operations of the right multiplication of 2-rowvectors on 2-matrices.Some known properties of ordinary matrices can be ex-tended to the 2-matrices. In particular, the 2-matrix B can be complexly conjugated by conjugating its projec-tions: B ∗ = B ∗ (cid:126) e + B ∗ (cid:126) e , transposed: B T = B T (cid:126) e + B T (cid:126) e , Hermitian conjugated: B † = B † (cid:126) e + B † (cid:126) e . Wedefine an (anti-) Hermitian 2-matrix if both its projec-tions are (anti-) Hermitian: B † = ∓ B , B † = ∓ B . Thesame is applicable and for 2-vectors. We also define aToeplitz 2-matrix if both of its projections are Toeplitz.Finally, a 2-matrix can be persymmetric if the equality J B J = B T , holds, where J is an m × m exchange ma-trix. We will also consider the zero 2-matrix O , bothprojections of which are zero matrices. Note that thedeterminant and inversion of 2-matrix are not defined.It is well known that ordinary non-singular matricesform a group with respect to the operation of matrix mul-tiplication. This group has a unit element, it is the unitmatrix, and an inverse element, it is the inverse matrix.The mentioned group properties of ordinary non-singularmatrices allowed successfully use them for solving linearsystems of equations. Vector matrices and their gener-alization, multidimensional matrices, have not found thesame acceptance in applied mathematics and mathemati-cal physics as ordinary matrices do, possibly because theydo not have the necessary group properties. In particu-lar, 2-matrices do not form a group with respect to theoperation of generalized (scalar matrix) multiplication,since such multiplication results not in a 2-matrix, butin an ordinary matrix. Since the set of 2-matrices do notform a group, it has not a unit 2-matrix and, accordingly,has not inverse 2-matrices. It turns out, however, thatthe group with respect to the generalized operation ofmultiplication, including ordinary matrix multiplication,multiplication of ordinary matrices by 2-matrices, andscalar multiplication of 2-matrices, forms a more complexconstruction of 4-block matrix whose diagonal blocks areformed by ordinary matrices, for example, A and B , andoff-diagonal blocks are 2-matrices, for example, B and C with projections sizes of m × m : M = (cid:18) A BC D (cid:19) , (16)Hereinafter, we denote such 4-block matrices by capitalRoman bold letters and will also call simply block ma-trices. For such block matrices, there is a unit element: E = (cid:18) E E (cid:19) , where E is the unit matrix of size m × m ,and is the 2-matrix with zero projections of m × m size.For non-singular (it will be clear later what it means) 4-block matrices there exists an inverse matrix having thesame 4-block form. We write the generalized Frobeniusformula [20, 21] for the inversion of our M block matrix(16): M − = (cid:18) A − + A − B · H C A − − A − B H − H C A − H (cid:19) , (17)where H = ( D − C A − · B ) − . It can also see from(17) that the inverse matrix has ordinary matrices onits diagonal, and the off-diagonal blocks are 2-matrices.We emphasize that the Frobenius formula does not re-quire the inversion of the diagonal 2-matrices B and C ,for which the inversion not defined. It follows from theFrobenius formula that an inverse matrix exists if thereexist matrices A − and H . Besides, the scalar productof 2-matrices is associative to ensure the equality of theleft and right inverse block matrices.The fact that non-singular block matrices form a groupwith respect to the generalized multiplication operationallows us to use them to solve systems of linear equa-tions, and also expand the range of applicability of somenumerical algorithms and approaches developed for ordi-nary non-singular matrices.To compose a linear system of equations with a blockmatrix we consider a block (column) vector p , denotedby a Roman bold letter, that is an analog of columnsfrom the left half of the block matrix. Its top part is an(ordinary) column vector c , with m size, and the bottompart is a 2-vector column b : p = (cid:18) c b (cid:19) . It is easy toverify, following the rules described above for multiply-ing ordinary and 2-matrices by an ordinary and 2-vector,that multiplying a unit block matrix E by a block vector p leaves the latter unchanged: Ep = p . There is anotherversion of the block vector, it is a “flip” block vector,that corresponds to the columns of the right half of the block matrix, for example, d = (cid:18) b c (cid:19) . It is also easy toverify that Ed = d . Similarly one can define block rowvectors. V. VECTOR TIB ALGORITHM FOR THEINVERSE SCATTERING PROBLEM
We turn to the system of equations (11) (12) using theblock vector notations described above. First we consider2-matrices R = R (cid:126) e + R (cid:126) e , where R , R are Toeplitz m × m matrices, m = 1 , ..., N . Also we define a 2-vectorcolumn of the solution of the system of equations (11)and (12) v ( m ) = v ( m )1 (cid:126) e + v ( m )2 (cid:126) e , and a 2-vector columnof the right-hand side r ( m ) = r ( m )1 (cid:126) e + r ( m )2 (cid:126) e , where thecolumn vectors r ( m ) are given by the discretized kernelsof the GLM equations, as in the scalar case in Eq. (15).It is required to find the 2-vector of the potential q = q (cid:126) e + q (cid:126) e . We denote its n th component as c-vector (cid:126)q n .In these notations, equations (11) and (12) like in scalarcase (14) can be represented in the form of a 4-blockmatrix: (cid:18) E ± R † R E (cid:19) (cid:18) u ( m ) v ( m ) (cid:19) = (cid:18) r ( m ) (cid:19) (18)From a comparison of (18) and (14) it follows that thevector case differs from the scalar one only by the cor-responding vector notation. The matrix of system (18)is precisely a 4-block matrix of size 2 m × m , the diago-nal blocks of which are formed by ordinary unit matrices E , and the off-diagonal blocks formed by 2-matrices R and ± R † . Since the latter are Toeplitz 2-matrices, theentire block matrix of the system G = (cid:18) E ± R † R E (cid:19) , isalso Toeplitz. Off-diagonal blocks of this matrix have ei-ther Hermitian or anti-Hermitian symmetry, dependingon the sign ± , which makes it possible to develop thevector versions of the TIB algorithms. The Levinson al-gorithm for Toeplitz matrices is not directly applicableto this problem. As in the scalar case, if we increase in-dex m by one an “inner bordering” occurs for 2-matrices R , ± R † and also for identity matrices E , each of themincreases by one column and one row. The block matrix G , in this case, increased by two rows and two columns,while in the Levinson algorithm at each step, the size ofthe matrix incremented by 1. In this case, the analogof the TIB algorithm for the inverse scattering problem,described in [9, 10] becomes applicable.Suppose that at the m th step of the algorithm weknow the solution in the form of a column block vector (cid:18) u ( m ) v ( m ) (cid:19) . It is required at the next step of the algo-rithm to find the solution v ( m +1) corresponding to theembedded system of equations of size 2 m + 2.We will obtain the solution of the system of linearequations if we find the inverse block matrix G − . Usingthe generalized Frobenius formula (17), we write formallythis inverse matrix in the form: G − = (cid:18) E ± R † H · R ∓ R † H − H R H (cid:19) , (19)where H = ( E ∓ R · R † ) − . Note that matrix E ∓ R · R † is Hermitian; therefore, matrix H is also Hermitian. Inaddition, the top diagonal block E ± R † H · R of G − isalso Hermitian, i.e. both diagonal blocks of the inversematrix are Hermitian. The off-diagonal blocks of the in-verse matrix are either Hermitian (for the upper sign,i.e., for the defocusing NLSE), or anti-Hermitian (for thefocusing NLSE). The symmetry properties of the blocksof the inverse matrix are important for constructing avector algorithm similar to the scalar TIB. For this al-gorithm, relations between the left f ( m )1 and right f ( m ) m block column of the inverse matrix, and also its top g ( m )1 and bottom g ( m ) m block row are of particular importance.The symmetry of blocks of the inverse matrix G − al-lows us to establish relationships between columns androws framing the matrix. Note that to solve the inversescattering problem it is required to find only one bottomblock row g ( m ) m of the inverse 4-block matrix. The knowl-edge of this row is sufficient to determine the m th compo-nent of the 2-vector of potential (cid:126)q m since it is determinedby only one last element (cid:126)v ( m ) m (it is c-vector) of the solu-tion 2-vector (cid:126)v ( m ) : (cid:126)q m = 2 (cid:126)v ( m ) m = 2 g ( m ) m (cid:18) r ( m ) (cid:19) .For what follows, it is convenient to introduce for m thstep of the algorithm the column vector y ( m ) and thecolumn 2-vector z ( m ) = z ( m )1 (cid:126) e + z ( m )2 (cid:126) e . Let the first(left) block column f of the inverse matrix consist ofthese vectors: f ( m )1 = (cid:18) y ( m ) z ( m ) (cid:19) . The Toeplitz symmetryof the original G matrix ensures the persymmetry of theinverse matrix. Therefore, the bottom block row g ( m ) m ofthe inverse matrix G − is a symmetric reflection of theleft block column f ( m )1 with respect to the northeast-to-southwest diagonal: g ( m ) m = (cid:32) (cid:103) z ( m ) (cid:103) y ( m ) (cid:33) T , where the tildemeans the inverse of the numbering of the elements of arow or column. Due to the Hermitian symmetry of thediagonal blocks of the inverse matrix, the left part of thetop half of the block row g ( m )1 of the inverse matrix, thatis part of its top diagonal block, is Hermitian conjugateto the first half of the left column y ( m ) and has the form y ( m ) † . For off-diagonal blocks of the inverse matrix inthe Hermitian/anti-Hermitian cases, the 2-vector row ofthe top block row g ( m )1 is Hermitian or anti-Hermitianconjugate to the left column vector z ( m ) . Thus, the topblock row of the inverse matrix can be represented as g ( m )1 = (cid:18) y ( m ) ± z ( m ) (cid:19) † . Taking into account the persym-metry of the inverse matrix, we see that the right (last) block column f m ) m of the inverse matrix is persymmetricto the top (first) row g ( m )1 , that is f ( m ) m = (cid:32) ± (cid:103) z ( m ) (cid:103) y ( m ) (cid:33) † .Comparing the obtained rows and columns of inversematrices for the scalar and vector cases of the GLME, wecome to a vector generalization of the TIB inverse scat-tering algorithm, which differs from the scalar one onlyin that some of the vector arrays (for example, y ( m ) ) inthe algorithm remain the same, but another one( z ( m ) )becomes 2-vector, i.e. as if they became doubled. Theauxiliary vectors y ( m ) and z ( m ) at m th step of the algo-rithm are calculated on the basis of the equations: y ( m +1) = c ( m ) (cid:18) y ( m ) (cid:19) + (cid:126)d ( m ) · (cid:32) (cid:126) ± (cid:94) z ( m ) ∗ (cid:33) (20) z ( m +1) = c ( m ) (cid:18) z ( m ) (cid:126) (cid:19) + (cid:126)d ( m ) (cid:32) ± (cid:93) y ( m ) ∗ (cid:33) (21)Here (cid:126) ( m ) is acomplex scalar, and (cid:126)d ( m ) = d ( m )1 (cid:126) e + d ( m )2 (cid:126) e is a c-vector.Equation (20) is scalar, and (21) is a vector equation,i.e. these are two equations for two components of the2-vector (cid:126)z ( m +1) .For the f ( m )1 and right f ( m ) m block columns of the inverseblock matrix, we can write: Gf ( m )1 = (10 ... ) T , Gf ( m ) m = ( ... T . (22)Here is a zero 2-vector. Comparing equations (22) forstep m and for step m + 1, we arrive at the followingequations for the coefficients c ( m ) , (cid:126)d ( m ) :c ( m ) ± (cid:126)β ( m ) ∗ · (cid:126)d ( m ) = 1 , c ( m ) (cid:126)β ( m ) ± (cid:126)d ( m ) = (cid:126) . (23)Here one equation is also scalar and the other is vector,i.e. compact record of 2 equations. The solution to thesystem (23) has the form:c ( m ) = (1 ± | (cid:126)β ( m ) | ) − , (cid:126)d ( m ) = (cid:126)β ( m ) c ( m ) (24)The main parameter of the TIB algorithm c-vector (cid:126)β ( m ) has the components: (cid:126)β ( m ) = β ( m )1 (cid:126) e + β ( m )2 (cid:126) e , andis given by: (cid:126)β ( m ) = m − (cid:88) k =0 h(cid:126) Ω m − k y ( m ) k (25)This equation is also a compact vector representation ofa pair of equations. As a result, for the m th componentof the 2-vector of the potential, we obtain: (cid:126)q m = 2 (cid:126)v ( m ) m = − (cid:126)β ( m ) /h. (26) VI. SCHEMATIC OF THE VECTOR TIBALGORITHMSA. Vector TIB algorithm for the inverse scatteringproblem
First-order algorithm for the inverse scattering prob-lem includes the following steps:1. Put m = 1 and calculate initial value for 0thcomponent of the solution vector (it is a c-vector) (cid:126)q = − (cid:126) Ω and initial values for auxiliary vectors: y (1)0 = (1 ± h | (cid:126) Ω | ) − , (cid:126)z (1)0 = − y (1)0 h(cid:126) Ω . (27)2. Determine the main parameter of the algorithm c-vector (cid:126)β ( m ) using (25).3. Find m th component of the potential vector (cid:126)q m from Eq. (26); this is the output at every step.4. Calculate coefficients c ( m ) and (cid:126)d ( m ) from Eq. (24).5. Determine auxiliary vector y ( m ) and 2-vector z ( m ) using (20) and (21).6. Increment m and go to the step 2 until m < N . B. Vector TIB algorithm for the direct scatteringproblem
Described above algorithm for the inverse scatteringproblem can be inverted to solve the direct scatteringproblem. Resulting algorithm consists of the followingsteps:1. Calculate initial values the kernel vector (cid:126) Ω = − (cid:126)q /
2, and initial values for the auxiliary vectors(27) and put m = 1.2. Determine the main parameter of the algorithm: (cid:126)β ( m ) = − h(cid:126)q m / m th component of the kernel vector (cid:126) Ω m (thisc-vector is the output at every step): (cid:126) Ω m = (cid:32) (cid:126)β ( m ) − h m − (cid:88) k =1 (cid:126) Ω m − k y ( m ) k (cid:33) /y ( m )0 .
4. Calculate coefficients c ( m ) and (cid:126)d ( m ) from Eq. (24).5. Determine auxiliary vector y ( m ) and 2-vector z ( m ) using (20) and (21).6. Increment m and go to the step 2 until m < N . VII. NUMERICAL SIMULATION:ALGORITHMS VERIFICATION
Numerical simulation was performed to test the pre-sented vector TIB algorithms for the Manakov vectorsoliton as an example of an accurate solution. Recall thatthe Manakov vector soliton corresponds to one eigenvalue λ = ω + i a of the discrete spectrum of the Manakov sys-tem:q α ( x ) = − α a sech(2 ax + δ ) exp( − ωx + i θ ) , (28)where l α are components of c-vector of polarization ofthe soliton, a is amplitude, δ is its center displacement, ω is the frequency, θ is its phase, and as elsewhere in thetext α = 1 , . This soliton corresponds to the next components ofvector kernel of the GLMEΩ α ( x ) = c α exp( − i λx ) = c α exp( − i ωx + ax ) . (29)Here c α = 2l α a exp { ( δ + i θ ) } are components of the com-plex vector constant that determines amplitude, shift,phase and polarization of the soliton.Numerical modeling confirmed the efficiency of the vec-tor TIB algorithms. Some of the calculation results, us-ing the example of the Manakov vector soliton, are shownin Fig. 1–4. The calculations performed by the variantof the program for solving the inverse scattering problemon the interval − L/ ≤ x ≤ L/
2. The solution obtainedon the range [ − ,
10] for N = 2 = 8192 calculationintervals, and soliton shifted from the center of the inter-val to test the asymmetric solution. The polarization ofthe soliton chosen so that the real parts of componentsof the soliton potential have different signs. The exactand restored from the GLM kernel real parts of the Man-akov 2-vector potential (28) are presented in the figure1. One polarization component of the Manakov solitondisplayed above the abscissa axis, and below there is an-other one. The calculations carried out with first-orderapproximation accuracy.Figure 2 presents the distribution of the absolute valueof the solution error for the inverse scattering problemfor N = 2 calculation intervals. The maxima of theabsolute error in the figure correspond to the peaks ofthe derivative solution. It follows that the main is theapproximation error. The error fall at the end of the in-terval confirms that it not accumulated in the algorithm.The integral error of inverse problem solution for N = 2 calculation intervals was 0.0028, and for 2 intervals, itwas 0.0012. The ratio of these values is 2, which indicatesthe 1st order of approximation accuracy. The calcula-tion time increased four times that correspond to O ( N )floating-point operations.Figures 3 and 4 presents calculation result for the di-rect scattering algorithm for N = 2 calculation inter-vals. Figures 3 shows a comparison of the exact and re-stored logarithms of the absolute value of integral kernel2-vector Ω α for both polarization components. Figure - - - - - x R e @ q H x L D , R e @ q H x L D FIG. 1: Inverse scattering problem: comparison of the exact(gray curve) and restored (black strokes) real parts of twopolarization components of the potential 2-vector q for theManakov vector soliton. Real part of the component q ofthe potential vector is placed above the abscissa axis and realpart of the component q lies below the abscissa axis. - - x R e s t o r e E rr o r s FIG. 2: Inverse scattering problem: distribution of absolutecalculation errors for both polarization components of the po-tential 2-vector q of the Manakov vector soliton. Gray curverefers to the component q , and black strokes correspond tothe component q . Ω α . It can see from the figure that,starting from about the middle of the Manakov soliton,the relative calculation error increases almost linearly. Itindicates a moderate error accumulation by the TIB algo-rithm for solving the direct scattering problem. However,when the number of calculation intervals doubled, the in-tegral error also halves, which confirms the first order ofapproximation accuracy. - - - - - x L og È W Α H x L È , Α = , FIG. 3: Direct scattering problem: Comparison of the exact(gray curve) and restored (black strokes) logarithms of theabsolute value of integral kernels Ω α for two its polarizationcomponents. - - x R e l a ti v e R e s t o r e E rr o r s FIG. 4: Direct scattering problem: distribution of the rela-tive calculation errors of integral kernels Ω α , α = 1 , , the black strokes correspond to component Ω ofkernel vector. VIII. CONCLUSION
Based on the discovered group properties of 4-blockmatrices with vector-like off-diagonal matrix blocks, ageneralization of the efficient scalar TIB algorithms forsolving inverse and direct spectral scattering problemsfor the vector nonlinear Schroedinger equation (Manakovsystem) is presented. Similar to the scalar algorithms,the new algorithms are based on solving the discretizedsystem of coupled GLM integral equations for both fo-cusing and defocusing cases. Also, as in the scalar case,the acceleration of calculations is achieved due to theToeplitz symmetry of the matrix of the discretized sys-tem of GLM equations. The algorithms were tested onthe exact NLSE analytical solution, Manakov vector soli-ton, and demonstrated high speed, stability, and accu-racy, which is sufficient for many applications.
IX. ACKNOWLEDGEMENTS
The author is grateful to Professor D. A. Shapiro andProfessor S. K. Turitsyn for helpful discussions, recom- mendations, and interest in this work.Funding: This work was supported by the Ministryof Science and Higher Education of Russian Federa-tion, project AAAA-A17-117062110026-3, and SectionVII (Numerical Simulation) was supported by RussianScience Foundation (RSF) (17-72-30006). [1] S. Novikov, S. Manakov, L. Pitaevskii, and V. Za-kharov,
Theory of solitons: the inverse scattering method (Springer Science & Business Media, 1984).[2] M. J. Ablowitz and H. Segur,
Solitons and the inversescattering transform , Vol. 4 (SIAM Studies in AppliedMathematics, Philadelphia 1981).[3] V. Zakharov and A. Shabat, Soviet Physics JETP , 62(1972).[4] L. F. Mollenauer and J. P. Gordon, Solitons in opticalfibers: fundamentals and applications (Academic Press,2006).[5] Y. S. Kivshar and G. Agrawal,
Optical solitons: fromfibers to photonic crystals (Academic press, 2003).[6] S. V. Manakov, Soviet Physics JETP , 248 (1974).[7] A. M. Basharov and A. I. Maimistov, Soviet PhysicsJETP , 913 1984.[8] A. I. Maimistov, Quantum Electronics , 756 (2010).[9] O. V. Belai, L. L. Frumin, E. V. Podivilov, and D. A.Shapiro, Journal of the Optical Society of America B ,1451 (2007).[10] L. L. Frumin, O. V. Belai, E. V. Podivilov, and D. A.Shapiro, Journal of the Optical Society of America B ,290 (2015).[11] R. E. Blahut, Fast Algorithms for Digital Signal Process- ing (Addison-Wesley, 1985).[12] A. Buryak, J. Bland-Hawthorn, and V. Steblina, Opt.Express , 1995 (2009).[13] O. V. Belai, L. L. Frumin, E. V. Podivilov, and D. A.Shapiro, Laser Physics , 318 (2010).[14] S. K. Turitsyn, J. E. Prilepsky, S. T. Le, S. Wahls, L. L.Frumin, M. Kamalian, and S. A. Derevyanko, Optica ,307 (2017).[15] V. Aref, S. T. Le, and H. Buelow IEEE Journal of light-wave technology , 1289 (2018 ).[16] L. L. Frumin, A. Gelash, and S. K. Turitsyn, PhysicalReview Letters , 223901 (2017).[17] S. Le, Y. Prylepskiy, and S. Turitsyn, Optics Express ,26720 (2014)[18] V. I. Nayanov, Multi-Field Solitons (Fizmztlit, Moscow2006, in Russian).[19] G. L. Lamb Jr.,
Elements of soliton theory (New York,Wiley-Interscience, 1980).[20] F. R. Gantmacher,
Theory of Matrices . (AMS ChelseaPublishing: Reprinted by American Mathematical Soci-ety, 2000).[21] D. Bernstein,