Adaptive Absorbing Boundary Layer for the Nonlinear Schrödinger Equation

: We present an adaptive absorbing boundary layer technique for the nonlinear Schrödinger equation that is used in combination with the Time-splitting Fourier spectral method (TSSP) as the discretization for the NLS equations. We propose a new complex absorbing potential (CAP) function based on high order polynomials, with the major improvement that an explicit formula for the coefficients in the potential function is employed for adaptive parameter selection. This formula is obtained by an extension of the analysis in [R. Kosloff and D. Kosloff, Absorbing boundaries for wave propagation problems, J. Comput. Phys. 63 (1986), no. 2, 363–376]. We also show that our imaginary potential function is more efficient than what is used in the literature. Numerical examples show that our ansatz is significantly better than existing approaches. We show that our approach can very accurately compute the solutions of the NLS equations in one dimension, including in the case of multi-dominant wave number solutions.


Introduction
In this paper we study a new complex potential absorbing boundary layer approach for the numerical solution of nonlinear Schrödinger (NLS) equations.We focus on the one-dimensional case, where the NLS, in scaled form, reads Here V is a real-valued potential and f is a real-valued function defining a local nonlinear interaction term.A standard case is, e.g., f(z) = qz, the cubic nonlinearity, in which case we have the Gross-Pitaevskii equation (GPE) that models Bose Einstein Condensates.Note that f(z) = qz 2 corresponds to a quintic nonlinearity.The equation is posed as IVP for x on whole space.Applications of NLS type equations can be found in many fields, such as quantum physics [3,13,32] (where TDDFT is a large class of widely used systems of NLS), nonlinear optics [26,45], plasma physics [16,57] or electromagnetic waves [33].
For a numerical computation of the NLS equation (1.1), e.g., with the widely used finite difference schemes, one faces the difficulty of the setting on whole space, whereas the computation can only be performed on a finite domain.Usually the domain of computation is restricted to some finite interval [a, b], and the problem is supplemented with "trivial" boundary conditions, i.e. periodic boundary conditions, with "zero" boundary conditions (i.e.homogeneous Dirichlet boundary conditions) as a special case.These artificially introduced boundaries in general result in errors as soon as the solution is not "small enough" near that boundary.These errors can be dealt with by two basically different approaches.
The first one consists of using a reformulation of (1.1) into a problem which can be actually computed on the bounded interval by applying so-called "Transparent boundary conditions (TBC)" on a, b.In this method, the Schrödinger equation (1.1) posed on the bounded interval is coupled to the same problem posed on the complement of [a, b], and then the two parts are coupled by demanding transparency on the boundary.This results in the usage of (e.g.) "Dirichlet-to-Neumann operators".
Another approach is called "Absorbing Boundary Layer (ABL)" method.Broadly speaking, there problem (1.1) is modified in such a way that errors at the boundary points a, b are either suppressed or totally removed.A "boundary layer" is added outside the original computation domain and the equation is modified on this boundary region in such a way that any parts of the solution inside this region will be damped to zero, while at the same time the solution stays unchanged on the computational interval [a, b].Thus the computation can be cut off outside of the boundary layer by trivial boundary conditions such as periodic boundary conditions, with a "reasonably small" small error with some kind of control.
The TBC method has been widely studied for the linear Schrödinger equation case (f = 0 in equation (1.1)).Exact Transparent Boundary Conditions (TBC) have been derived, implemented and analyzed for the Schrödinger equation in [4,9,19,24,39,46].The TBC equations are challenging to solve numerically, e.g., they contain convolutions in time direction, which are computationally expensive, and the stability of approximations is non-trivial.One way to treat this is to introduce fast algorithms for evaluating the time convolution [11,28,35].Another way is to approximate the transparent boundary conditions via rational approximations to the dispersion relation, resulting in approximate TBC, which are local in time but allow small reflections into the computational domain [1,2,18,21,42].Recently, an approximative method for transparent boundary conditions was described allowing for 3d simulations [52].The approximate TBC method has also been carried over to the density matrix level in [34].
For the nonlinear case, only approximate TBC have been derived so far.Approximate TBCs involving time convolution for the cubic NLS equation were designed by different authors via pseudodifferential, potential and para-linear strategies [5-8, 47, 48, 60].In the particular case of the cubic NLS equation, exact TBC can be constructed by inverse scattering [22,58] which involve time convolution; however, this approach is not applicable to more general nonlinearities.In [53], a split local approximate TBC was proposed for the NLS equation, which involves a time-splitting method and the application of the TBC to the linear sub-problem.An adaptive version of this method was developed in [54].
In the absorbing boundary layer (ABL) category, a number of different methods exist.One method, which we will follow in our work, is the addition of a complex absorbing potential (CAP), sometimes called "Optical potential", which has been studied in [31,36,37].An extensive review of this method class was done by Muga, Palao, Navarro and Egusquiza [36].The imaginary potential ABL has also been used to treat the NLS equation [27], and an extensive study of its application in the context of TDDFT physics calculations was done by DeGiovannini, Larsen and Rubio in [23].The method has been systematically studied in the context of strong field physics in [43,56].
Another method from this category is the perfectly matched layer (PML) method, and the related Exterior Complex Scaling (ECS) method.Both of these methods use a variable re-scaling into the complex plane, see the seminal work of Scrinzi, Stimming and Mauser [41] for a discussion of the relation between these two approaches and an application to the NLS.PML has been originally proposed by Bérenger for hyperbolic problems in electromagnetism [15], and applied to the linear Schrödinger equation, e.g., in [17].The ECS method has been successfully used by Scrinzi [40] and extended to "infinite range" ECS (IRECS), where the domain cutoff in the outer region is no longer needed.See [51] for a study of this method with non-uniform FEM for a laser physics application.
In the work of Soffer and Stucchio [44], the time dependent phase space filter method is introduce for the NLS equation, which conceptually also belongs to the ABL approach.An advantage of this method is its compatibility with the Fourier spectral method for the spatial discretization.
Another recent method to deal with the artificial boundary problem was devised by Kaye and others.They use a transform onto basis functions defined on a particular contour in the complex plane, thus calling the method "contour deformation", [29,30].In [30] the method was successfully applied for TDDFT calculations.
In this paper we study a new complex absorbing potential (CAP) method with a focus on the one-dimensional NLS (1.1) which is coupled to the TSSP (Time Splitting Spectral Propagator, see [12]) scheme for highly accurate computation of the solutions of the NLS.Aside from general computational simplicity, the possibility to use an efficient time splitting method is one advantage of the CAP approach.This is a major advantage in view of the performance difference that exists in comparisons between TSSP and other numerical methods for the NLS equation, for example see [49,50].
Both the PML/ECS approach and the TBC methods can not be combined with TSSP.In PML or ECS, the model equation contains variable coefficients, which would introduce computationally complex convolution terms in a spectral representation, and therefore need to be treated by different space discretization methods, e.g., finite differences or finite elements.Non-uniform FFT algorithms can be applied for this case.So even if the CAP method is less efficient than PML as absorbing boundary layer approach, namely it requires a larger layer width in order to ensure the same accuracy as ECS, its combination with the TSSP as solver/propagator for the original equation (1.1) offers an advantage considering that the PML is solved by finite differences or FEM.From this viewpoint, it is expected that the relatively efficient CAP combined with TSSP can be a satisfying tool for computing NLS equations on an unbounded domain.Note that our algorithm can be straightforwardly applied also in more than one space dimension.
The efficiency of the CAP method (in terms of error suppression) depends very much on the chosen complex potential function.In the works [27,31,36,37,56] different choices for the imaginary potential have been studied.In this paper we propose a new complex potential function based on a high order polynomial.We give an explicit formula for selecting the scale parameter of the imaginary potential function adaptively depending on the time evolution, by the computation of transmission and reflection coefficients for the free Schrödinger imaginary potential model (FSIPM) [31].Numerical examples are presented showing that this new complex potential with adaptively chosen coefficients is much more efficient than existing imaginary potential functions.This work is a step forward designing highly efficient CAP for computing NLS equations.
This paper is organized as follows.In Section 2, we introduce the CAP approach for NLS and a fourth order TSSP for the propagation of NLS equations on a finite computation interval.In Section 3 we discuss the computation of transmission and reflection coefficients for the FSIPM.We give a new transmission coefficient formula which cures the instability problems of the previously used formula and verify that the new formulas for both reflection and transmission coefficients give reliable results for the FSIPM.In Section 4 we present our new complex potential formula together with the existing ones in the literature.We use the FSIPM to guide the adaptive choice of reliable coefficients of the new absorbing potential function.Section 5 contains numerical examples to show that the coefficients determined by the FSIPM are reliable relative to the optimal ones and the new complex potential is much more efficient than existing ones for computing one-dimensional NLS.Section 6 concludes.

The Complex Absorbing Potential Model and the Time-Splitting Spectral Propagator
In this section we introduce the CAP approach for the one-dimensional NLS equation and the TSSP for solving the model equation.Let

The Complex Potential Model
In the complex potential (CAP) approach, [27,31,36,37,43,56] ) In Section 4 we will give concrete examples of absorbing complex potential functions.

A Fourth Order Time-Spitting Spectral Method
We now describe a fourth order TSSP for solving the complex potential equations (2.1)-(2.4).We choose the spatial mesh size Δx = 1 N (d − c) for an even positive integer N, the time step Δt > 0 and let the spatial and time grid points be Let u n j be the approximation of u(x j , t n ), and let u n be the solution vector at time t n with components u n j .By performing first order time-splitting of equation (2.1) one obtains two subproblems ) Equation (2.5) is the original NLS equation (1.1), while equation (2.6) corresponds to the reduction of the amplitude of the solution since the coefficient −iV C on the right-hand side has negative imaginary part.Since equation (2.6)only takes effect in the absorbing layer on which one is not concerned with the accuracy of the solution, it is not necessary to use a high order splitting for the above two equations.
Therefore we simply use first order splitting.Let S 1 (u I ; t), S 2 (u I ; t) be the numerical propagation operators for equations (2.5) resp.(2.6) with initial value u I up to time t.Then the algorithm for computing u n+1 from u n is given by ) Denote S 2,j (u I ; t), j = 0, 1, . . ., N, to be the components of S 2 (u I ; t).The operator S 2 is then given by In order to obtain S 1 , one can use the conventional time splitting for the NLS equation (2.5) ) Let S 11 (u I ; t), S 12 (u I ; t) be numerical propagators for equations (2.10) and (2.11), respectively, and let S 11,j (u I ; t), S 12,j (u I ; t) be their components.Since (2.10) is solved with periodic boundary conditions, it follows that S 11,j (u I ; t) can be computed by the Fourier spectral method, and S 12,j (u I ; t) is given exactly by integration of (2.11).
In summary the TSSP for computing the complex potential model (2.1)-(2.4) is given by (2.7)-(2.9),(2.12), which as a scheme has spectral accuracy in space and fourth order accuracy in time for the solution in the physically interesting interval [a, b].

Choosing the Absorbing Potential, FSIPM
In the work [31], the authors studied the CAP approach for the linear Schrödinger equation.In the case of linear equation and imaginary potential, one has f( ⋅ ) = 0 in (2.1), and puts iV C = V I for some imaginary valued potential V I .Outside of the region of interest [a, b], the potential V(x) can be assumed to be constant, the complex potential model in the boundary layer [c, a] ∪ [b, d] is essentially described by propagating the free Schrödinger equation with added absorption term In [31] the authors propose to choose the parameters in the imaginary potential by examining the transmission and reflection coefficients for time harmonic plane wave solutions of the above equation ("FSIPM").By choosing parameters which result in both small transmission and reflection coefficients, one can obtain an imaginary potential function with good absorption for the linear equation related to the given layer width and wave number of the chosen solution.
In this paper we propose to follow this strategy to choose reliable parameters for complex absorption potentials in the one-dimensional NLS case.Consider the imaginary potential case, namely iV C = V I in (2.1).By adopting time splitting, the equation is split into two sub-problems.The first is the linear problem ("FSIPM") (3.1) on the computational interval, the second is From this viewpoint it is expected that the parameter choosing strategy done in [31] giving small transmission and reflection coefficients in the above linear model also works reliably well for the imaginary potential layer in the NLS.This strategy depends on the computation of transmission and reflection coefficients for time harmonic plane wave solutions of (3.1).For certain particular imaginary potential functions, these can be obtained analytically [31,32].For general imaginary potential functions, however, they need to be computed numerically.This can be achieved by the propagator matrix method [25,31,32].In this section we verify that the propagator matrix method gives reliable results for the free Schrödinger model (FSIPM) by testing an imaginary potential function for which analytical expressions of the transmission and reflection coefficients are available.
Assume the boundary layer is located in the interval [0, D], on which the potential function V I is given, and let V I be zero outside of this interval.For the method, V I is approximated by a piece-wise constant function.Let the interval [0, D] be covered by a uniform mesh 0 = x 0 < x 1 < ⋅ ⋅ ⋅ < x M = D with mesh size h, so h = D M .Then the imaginary potential function V I (x) is approximated by where x j+ 1 2 = 1 2 (x j + x j+1 ).At the left boundary x = 0 we assume the solution has the form e i(Kx−Wt) + Re i(−Kx−Wt) , where K > 0 represents an incoming plane wave with wave number K, and R is the coefficient of the reflection part.R is complex in general.At the right boundary x = D, there is only a transmission wave, thus the solution has the form Te i(Kx−Wt) , with T being the transmission coefficient.On each subinterval [x j , x j+1 ], the solution consists of left and right traveling waves and has the form From the dispersion relation of (3.1), one has (3.4) Now one can establish a recurrence relation for the coefficients A j , B j .Consider a grid point x j+1 for 0 ≤ j ≤ M −2.
By imposing continuity of the solution and its first order derivative at this grid point, one obtains ) Equivalently one has for 0 ≤ j ≤ M − 2.
In the same way, one also has Denote the matrices in (3.7) to be G j , 0 ≤ j ≤ M − 2, and the matrices in (3.8) and (3.9) as G −1 and G M−1 respectively.Denote which is called the propagator matrix.Then the above sums up to The above formulas (3.4), (3.10), (3.12) and (3.13) thus can be used to compute the transmission and reflection coefficients T, R for fixed wave number K and layer width D. We now provide an example to verify that these for-mulas give reliable results for the transmission and reflection coefficients of the FSIPM.Consider the imaginary potential function When D is sufficiently large so that U(cosh 2 ( αD 2 )) −1 is essentially zero, the transmission and reflection coefficients of the FSIPM with the imaginary potential (3.14) can be expressed by (see [31,32]) T, (3.15) where Γ denotes the Gamma function and S is given by We then tested the propagator matrix method with M = 1000.We found that numerical results match very well with the analytical ones.It can also be noticed that the computed transmission coefficient is even reliable when the analytical one is below machine precision.This is due to the structure of the transmission coefficient formula (3.13).When T is very small, then G 22 is very large.In particular, for small values of U, the reflection coefficient is near or below the round-off errors, which is reasonable.
This example illustrates that the propagator matrix method should give reliable results for the transmission and reflection coefficients of the FSIPM.

A New Complex Potential Function with Adaptive Parameter Selection
As stated above, the choice of the complex potential function is essential to the efficiency of the absorbing layer.
In this section we propose a new choice of complex potential function, which will be demonstrated to be more efficient than previously used ones.First we consider complex potential functions in the right absorbing layer [b, d].The corresponding potential functions in the left absorbing layer [c, a] can be treated similarly.In [31,37], the following two imaginary potential functions have been used where U, α, C L are real-valued, positive parameters.Because these potential functions need to vanish in the physically interesting region, they are (approximately) zero at the point b.Within the absorbing layer, the absolute value transitions gradually from zero to the maximum at point d.The gradual transition is necessary to reduce reflection waves in the absorbing layer, as discussed in [31].
In this paper we propose the following complex potential function, a polynomial of fifth order.We first consider the purely imaginary function where C P is a positive number.This function has high regularity at the point b considering that V P I has to vanish on [a, b].
We use the FSIPM (Free equation model) discussed in Section 3 to study the practical performance of the imaginary potential functions.We will employ the propagator matrix method presented in Section 3. We then compare the CAP (4.1)-(4.2) and the new CAP (4.3).We choose K = 10, D = 10 in the Free Schrödinger method and M = 1000 in the propagator matrix method.From Figures 1-3 it can be seen that the minimum of the quantity max(|R|, |T|) for the linear, hyperbolic secans and the new (quintic) imaginary potential functions are approximately 10 −3 , 10 −7 and 10 −9 , respectively.So the new choice (4.3) is more effective in absorbing both transmission and reflection waves than the sech potential (4.1), which is more efficient than the linear function (4.2).In the next section we will present numerical examples showing that for computing one-dimensional NLS, the proposed potential is indeed more efficient than the other choices.the logarithmic curve of these quantities should be similar.Therefore it is not essential which one of these possible minimizing quantities is used in selecting the effective parameters.

Adaptive Parameter Selection by Local Wave Number
Our next task is to find an explicit formula for choosing an optimal coefficient C P in the new imaginary potential function (4.3) depending on wave number K and layer width D for practical computations of NLS equations.This optimal coefficient will then be applied to formulate an ABL in practical computations, where the dominant wave number K of the solution in the ABL is to be detected adaptively.Firstly, for fixed wave number K and layer width D, one can obtain an approximately optimal coefficient C P for the FSIPM by computing max(|R|, |T|) corresponding to the parameter choices made.For example, from Figure 3 one obtains that for K = 10 and D = 10 the optimal coefficient is C P = 10.In the same principle, we can obtain approximately optimal coefficients for different values of D for fixed K = 10.These results are shown in Table 1.
Let C P (K, D) be an explicit formula for choosing C P depending on K and D which we are trying to construct.From Table 1, by using curve fitting one obtains an expression when K = 10.One choice is C(10, D) = min (22, 88.042 To derive the full expression for the approximately optimal coefficient C P (K, D), we will use a property of the optimal coefficient.Consider using the propagator matrix method to compute T and R for the new imaginary potential function (4.3).Now observe that the entries in the propagator matrices (3.7)-(3.9)are unchanged under the following rescaling: let β ∈ ℝ + , and for fixed M, set Also the matrices G −1 , G 0 , . . ., G M−2 , G M−1 defined in (3.10) are unchanged, as well as the products Kh, KD.Therefore the results of the transmission and reflection coefficients (3.12)-(3.13)are unchanged under this rescaling.This fact leads to the following rescaling property of the optimal coefficient C P (K, D): One sees that for a wave number K and a layer width D, the corresponding minimum of the quantity max(|R|, |T|) only depends on the product KD.Therefore the smaller the wave number K is, the larger the layer width D needs to be in order to attain the same accuracy of the absorbing layer.Combining the expression (4.4) and the property (4.5), we get the expression for C P (K, D) as follows: Thus, formula (4.6) gives the approximately optimal coefficient of the new imaginary potential function (4.3) for any wave number K and layer width D for the Free Schrödinger model (FSIPM).In the next section we will show by numerical example that this formula gives practically reliable coefficients for the new imaginary potential function for one-dimensional NLS computations.
We now introduce a complex potential function by multiplying the potential (4.3) by a complex factor.By the discussion above, it is clear that the CAP is more efficient for higher wave numbers.Therefore we consider introducing a monotonically decreasing real potential in the absorbing layer which may increase the phase velocity and wave number of the solution.Thus we propose a complex potential by multiplying (4.3) by 1 − i: The optimal choice for C P is still given by the formula (4.6).In the next section we will show numerical examples which indicate that the complex potential (4.7) is more efficient than the imaginary potential (4.6) when the error tolerance of the absorbing layer is set to be not less than 10 −10 .In general, for practical computations one needs to deal with general solutions which do not always have a global wave number.In order to achieve a good absorption in this situation, the optimal parameter choice of the absorbing potential can be determined by adaptive detection of a dominant local wave number in the solution close to the boundary.We will use the energy-weighted method proposed in [54] to detect a dominant wave number from the solution at each time step and apply it for an adaptive use of the new potentials (4.3) or (4.7).
For the application of the scheme (2.5)-(2.6)on the whole interval [c, d], the absorbing potential (4.3) on both the left and right absorbing layers is formulated as Given some local dominant wave number K, which can be different in each absorbing layer and can change with time, and the ABL width D, the coefficients C P,1 , C P,2 in (4.8) can be chosen adaptively in an optimal way by formula (4.6).We need to determine these coefficients at each time step before applying the exact solver for the second subproblem (2.6) in the splitting (2.5)-(2.6),i.e. the step corresponding to the application of the absorbing potential.The layer width D of the left and right ABLs are a − c and d − b, respectively.The dominant wave number K in the two ABLs will be detected numerically.We propose to adopt the energy-weighted approach used in [54] for adaptively detecting the wave number K in the ABLs.Take the ABL [b, d] for example.Taking the windowed Discrete Fourier transform in this ABL, the discrete spectra of the numerical solution in the ABL are given by where , and [x] represents an integer choosing operation, such as choosing the largest integer not bigger than x.Then the dominant wave number K in this ABL obtained by the energy-weighted approach is given by where we choose p = 4 as suggested in [54].The dominant wave number in the ABL [c, a] can be similarly obtained by considering the weighted average of the negative wave numbers.Now suppose we get the dominant wave numbers K 1 , K 2 in the ABLs [c, a] and [b, d] respectively by the above energy-weighted approach.Then we can get the coefficients C 1 , C 2 in the imaginary potential function (4.8) by where C P ( ⋅ , ⋅ ) is the coefficient formula (4.6).In summary, for the evolution of the numerical solution u n to the next time step u n+1 , the time-splitting spectral process (2.12) on u n is used.We denote the resulting solution vector to be u * with components u * j .We then apply the exact solver for the second subproblem (2.6) on u * to get the numerical solution u n+1 at the next time step.Namely, we have where C 1 , C 2 are given by (4.10).

Numerical Examples
In this section we present numerical examples to verify the accuracy of our method and to demonstrate the efficiency of the absorption potential (4. 3) and the complex potential (4.7) compared with potentials (4.1) and (4.2), and the accuracy of its application with adaptive parameter selection.Our method is applicable to general nonlinearities.In this paper we only test the cubic nonlinearity case, namely f( ⋅ ) in (1.1) is defined by f(z) = qz.
Example 5.1.We first consider single soliton propagation problems tested in [60].We consider the cubic NLS equation with the focusing scale q = 2, and the time evolution of a single soliton solution u ana (x, t) = ρf(x, t)g(x, t), ( where ρ is the amplitude and ν is the group velocity of the soliton wave.We set the interval of interest to be We now introduce our approach to measure the error of the CPLs.Tables 2 and 3 list the L 2 -errors on [−10, 40] of the TSSP described in Section 2 using the new imaginary potential function (4.3) with the coefficient C P = 10 as shown in Table 1, and different choices of the layer width D, time step size Δt and computational time T. In Table 2 one observes that the numerical solutions converge in fourth order with temporal mesh refinement for all the tested layer widths D. This is because the computational time T = 0.24 is relatively short at which the essentially non-zero part of the solution does not reach the absorbing boundary layer.Thus the error of the absorbing layer is not relevant at this time.In Table 3 we present the numerical errors at a large computational time at which the soliton has completely disappeared from the computational interval.The analytical solution in the physically interesting interval is very small at this time.We see that the solutions with boundary layer width D = 30 converges in fourth order with time step refinement.This implies the error of the absorbing boundary layer with layer width D = 30 is less than about 10 −12 , thus the error of the solution computed with the smallest time step size Δt = 0.00125 is not influenced by this absorbing boundary layer.For other smaller layer widths, we see that the error of the solution with relatively fine time step sizes are influenced by the absorbing boundary layer.In particular, the error of the solution with the smallest time step size indicates the error of the corresponding absorbing boundary layer.This gives the approach for measuring the error of the CPL.Namely we take the L 2 -error on the physically interesting interval of the numerical solution by TSSP at the computational time T = 6 with the time step size Δt = 0.00125 as the error of the CPL provided the error is significantly larger than 10 −12 .
We do not plot the numerical solutions versus the analytical solution since the numerical errors are relatively small and the behavior of the analytical solution is easy to understand as given by (5.1)-(5.3).We now test the error of the absorbing boundary layer with different complex potential functions and parameters for the layer width D = 10.As performed for the coefficient of the imaginary potential function (4.3) in Table 1, we can also determine parameters for the hyperbolic secans and linear imaginary potentials (4.1) and (4.2) using FSIPM.For K = 10, D = 10, this results in the parameters U = 480, α = 0.91 for the sech potential; and C L = 33 for the linear one.Denote E L to be the error of the absorbing boundary layer.In Figures 4-7 we plot the quantities log 10 (E L ) for different complex potential functions with parameters chosen around the ones determined by FSIPM.It is seen that the new imaginary potential function and complex potential function are more efficient than the hyperbolic secans function, which in turn is more efficient than the linear function, in term of both the optimal parameters and the ones determined by the FSIPM.This observation is similar to that tested by computing the transmission and reflection coefficients of the FSIPM in Section 4.    In Figure 4, one observes that the complex potential function (4.7) is more efficient than the purely imaginary potential function (4.3) for the optimal parameter choice according to (4.6), and also for a large range of parameters around this value.
One also observes that the parameters determined by the FSIPM are relatively reliable compared with the optimal ones since the errors of the absorbing boundary layer with these parameters are no more than one order of magnitude multiple of those with the optimal parameters.with D being the layer width to be chosen.We choose the spatial grid size to be Δx = 1  15 .The initially left soliton formally catches up with the initially right soliton at the position x = 55, therefore both solitons reach the right ABL approximately at the same time.Figures 8 and 9 show the numerical solutions at T = 1.5 and 2 respectively before and when the two solitons are in interaction.We present in Table 4 the L 2 -errors of the numerical solutions at a large enough time T = 10 when the exact solution in the computational domain is essentially zero.As explained in the previous example, these errors reveal the error caused by using the ABL since at this time the reflection waves generated by the ABL should return to the computational domain.As shown in Table 4, these errors decrease to be sufficiently small with the increment of the layer width.With D = 40, the error decreases to the order of 10 −11 .Considering that in our method one only needs to detect one dominant wave number by the energy-weighted approach described in Section 4 even when the solution actually more than a single wave number, neither of which necessarily is dominant, it is seen that our method is very accurate and convenient in computing the multi-dominant wave number solutions of the cubic NLS equation.

Conclusion
In this paper we studied the coupling of the Time-splitting Fourier spectral method (TSSP) with the imaginary potential absorbing boundary layer(ABL) as a numerical method for computing accurate numerical solutions of the nonlinear Schrödinger (NLS) equation posed on the real axis by an approximation on a finite interval.
The choice of the imaginary potential function is essential in designing an effective imaginary potential ABL.
We proposed a new imaginary potential function based on the high order polynomial, and provided an explicit formula for the coefficient in the imaginary potential function used for adaptive selection of this parameter when using this imaginary potential ABL in practical computations.With this adaptive parameter selecting strategy, our imaginary potential ABL is convenient and reliable to use.The coefficient formula was obtained by using the model analysis adopted in [31].The model analysis also showed that our imaginary potential function is significantly more efficient than the existing ones.Numerical examples were presented showing that our approach works well for the NLS equations in one dimension, including highly accurately computing multidominant wave number solutions of the NLS equations.Our imaginary potential ABL is more efficient for high wave numbers, while it has lower efficiency for small wave number problems in which the ABL width needs to be relatively large.We noticed that an adaptive spatial mesh size selection strategy can be possibly used to compensate the low efficiency of our method for small wave number problems.
Funding: This work was supported by the FWF (Austrian Science Foundation) via project F65 (SFB "Taming complexity in nonlinear PDE systems") and by the WWTF (Viennese Science and Technology Fund), project MA16-066 "SEQUEX".
[a, b] the physically interesting interval on which the solution of the original Schrödinger equation (1.1)-(1.2) is sought.We assume that the initial data u 0 (x) is compactly supported on [a, b].In this situation the equation only contains outgoing waves at the boundaries of the interval [a, b].
, one chooses c, d such that c < a < b < d.The intervals [c, a] and [b, d] are the locations of the absorbing boundary layers.A complex potential term "iV C " is added to the equation, which is chosen to be zero on [a, b] and to have negative imaginary part on [c, a] ∪ [b, d].Then the following equations are solved on the enlarged domain [c, d]: ) with periodic boundary conditions on the ends of the interval [c, d].The model equation (2.1) reduces to the original NLS on [a, b], and on [c, a] ∪ [b, d], its solutions are damped in the time propagation.The errors of using CAP originate from two sources.One is that the solution in the absorbing layer can be reflected back into [a, b], and the other is transmission of parts of the solution not fully suppressed by the damping through the absorbing layer back into [a, b] due to the periodic boundary conditions (2.3)-(2.4).Therefore an effective absorbing boundary layer needs to perform well in absorbing both the transmission and the reflection waves.

. 9 )
which implies that the solution is damped on the absorbing layer [c, a] ∪ [b, d] and unchanged on [a, b].

. 2 )
Equation (3.2) preserves periodic boundary conditions of the solution.Therefore, if the imaginary potential function gives a good absorption for (3.1), it will also perform as absorbing boundary layer for the complete NLS in one time step of the splitting algorithm which involves solving the two sub-problems (3.1)-(3.2) subsequently.

Re iKh ) . ( 3 . 11 )
Using the fact that the determinant of the propagator matrix G is equal to one, from(3.11)one obtains

[− 10 ,
40] on which we seek the solution of the cubic NLS equation.Since this problem describes a soliton traveling from left to right, we do not put an absorbing layer at the left of the interval.Let the right ABL be [40, 40 + D].Then the computational interval is [−10, 40 + D].We choose the parameters for the analytical solution (5.1)-(5.3)x c = 15, ρ = 1.5, ν = 20, thus the dominant wave number of the solution is 10.The spatial mesh size is chosen to be Δx = 0.1.

Example 5 . 2 .
In this example we test the effect of our ABL for absorbing solutions with two different wave numbers.Initial data are taken to be the sum of two different solitons of the form (5.1)-(5.3).The left soliton has the parameters x c = 15, ν = 20, ρ = 1.5, and the right soliton has the parameters x c = 35, ν = 10, ρ = 1.5.The physically interesting interval is chosen to be [−10, 60].The left ABL is [−20, −10], and the right ABL is [60, 60 + D]

Figure 9 :
Figure 9: Example 5.2, computed |u| at T = 2 when the two solitons are in interaction.