Adaptive Absorbing Boundary Layer for the Nonlinear Schrödinger Equation

1 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 ) = q ⁢ z , the cubic nonlinearity, in which case we have the Gross–Pitaevskii equation (GPE) that models Bose Einstein Condensates. Note that f ⁢ ( z ) = q ⁢ z 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, 46, 39]. 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, 21, 18, 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, 6, 7, 8, 48, 47, 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.

The PML approach has been successfully applied to the NLS equation [10, 20, 38, 59].

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.

2 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 [ 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 ] .

3 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 i ⁢ V 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 i ⁢ V 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

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. 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 ⁢ ( K ⁢ x - W ⁢ t ) + R ⁢ e i ⁢ ( - K ⁢ x - W ⁢ t ) , 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 T ⁢ e i ⁢ ( K ⁢ x - W ⁢ t ) , 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

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

Using the fact that the determinant of the propagator matrix G is equal to one, from (3.11) one obtains

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 formulas 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])

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.

4 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 I P 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.

Figures 1–3 plot the computed curve for log 10 ⁡ ( max ⁡ ( | R | , | T | ) ) for different imaginary potential functions with varying parameters. In Figure 3 the quantity M R ⁢ T , U denotes

where α is chosen from I α = { 0.2 , 0.3 , 0.4 , … , 4 } and M R ⁢ T ⁢ ( U , α ) = log 10 ⁡ ( max ⁡ ( | R | , | T | ) ) , for fixed ( U , α ) .

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.

Figure 1

Computed curve log 10 ⁡ ( max ⁡ ( | R | , | T | ) ) versus the coefficient C for the imaginary potential function (4.3).

Figure 2

Computed curve log 10 ⁡ ( max ⁡ ( | R | , | T | ) ) versus the coefficient C for the linear imaginary potential function (4.2).

Figure 3

Computed curve M R ⁢ T , U versus the parameter U for the hyperbolic secans imaginary potential function (4.1).

5 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 ) = q ⁢ z .

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

(5.1) f ⁢ ( x , t ) = sech ⁢ ( ρ ⁢ ( x - x c - ν ⁢ t ) ) ,

(5.2) g ⁢ ( x , t ) = exp ⁡ ( i ⁢ ν 2 ⁢ ( x - ν ⁢ t ) ) ⁢ exp ⁡ ( i ⁢ ( ρ 2 + ν 2 4 ) ⁢ t ) ,

(5.3) u a ⁢ n ⁢ a ⁢ ( 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 [ - 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 .

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 .

Table 2

Example 5.1, ν = 20 , L 2 -errors on [ - 10 , 40 ] of the numerical solutions with different D, Δ ⁢ t , at computationaltime T = 0.24 .

	Δ ⁢ t
D	0.01	0.005	0.0025	0.00125
5	8.76E-7	5.60E-8	3.52E-9	2.22E-10
10	8.76E-7	5.60E-8	3.52E-9	2.22E-10
20	8.76E-7	5.60E-8	3.52E-9	2.22E-10

Table 3

Example 5.1, ν = 20 , L 2 -errors on [ - 10 , 40 ] of the numerical solutions with different D, Δ ⁢ t , at computationaltime T = 6 .

	Δ ⁢ t
D	0.01	0.005	0.0025	0.00125
5	6.20E-6	6.22E-6	6.23E-6	6.23E-6
10	2.45E-8	1.08E-8	1.07E-8	1.07E-8
20	2.19E-8	6.95E-10	4.55E-11	3.77E-12
30	2.17E-8	6.89E-10	4.52E-11	2.96E-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 the 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.

Figure 4

Example 5.1, comparison of errors with imaginary potential (4.3) and complex potential (4.7), □ : imaginary potential, ○ : complex potential.

Figure 5

Example 5.1, errors of absorbing boundary layers with linear potential (4.2).

Figure 6

Example 5.1, errors of absorbing boundary layers with hyperbolic secans potential (4.1), α = 0.91 fixed.

Figure 7

Example 5.1, errors of absorbing boundary layers with hyperbolic secans potential (4.1), U = 180 fixed.

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.

Figure 8

Example 5.2, computed | u | at T = 1.5 when the two solitons begin interacting.

Figure 9

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

6 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 multi-dominant 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.