Numerical Simulation of Soliton Tunneling

MAT-VET-F 20028

Examensarbete 15 hp

Juni 2020

Numerical Simulation of Soliton

Tunneling

Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student

Abstract

Numerical Simulation of Soliton Tunneling

Elias Estensen, Amanda Seger, Matilda Tiberg

This project studied two different ways of imposing boundary conditions weakly with the finite difference summation-by-parts (SBP) operators. These operators were combined with the boundary handling methods of simultaneous-approximation-terms (SAT) and the Projection to impose homogeneous Neumann and Dirichlet boundary conditions. The convergence rate of both methods was analyzed for different boundary conditions for the one-dimensional (1D) Schrödinger equation, without potential, which resulted in both methods performing similarly. A multi-block discretization was then implemented and different combinations of SBP-SAT and SBP-Projection were applied to impose inner boundary conditions of continuity between the blocks. A convergence study of the different methods of imposing the inner BC:s was conducted for the 1D Schrödinger equation without potential. The resulting convergence was the same for all methods and it was

concluded that they performed similarly. Methods involving SBP-Projection had the slight advantage of faster computation time.

Finally, the 1D Gross-Pitaevskii equation (GPE) and the 1D Schrödinger equation were analyzed with a step potential. The waves propagating towards the potential barrier were in both cases partially transmitted and partially reflected. The waves simulated with the Schrödinger equation dispersed, while the solitons simulated with the GPE kept their shape due to the equations reinforcing non-linear term. The bright soliton was partly transmitted and partly reflected. The dark soliton was either totally reflected or totally transmitted.

Populärvetenskaplig sammanfattning

Partiella differentialekvationer används för att beskriva många fenomen i naturen. Inom kvantmekaniken används ofta Schrödingerekvationen. En tolkning av denna ekvation är att den beräknar var det är mest sannolikt att en kvantmekanisk partikel befinner sig, då detta ej går att räkna ut exakt om man samtidigt känner till dess rörelsemängd på grund av Heisenbergs osäkerhetsprincip. En annan ekvation av denna typ är Gross-Pitaevkiiekvationen vilket är en ekvation som beskriver solitoner i ett Bose-Einsteinkondensat. En soliton är en typ av våg som inte diffunderar under tiden den propagerar. Till exempel kan en tsunami långt ifrån land modelleras som en soliton. Om en partiell differentialekvation inte kan lösas analytiskt krävs en numerisk lösning. Det finns olika numeriska metoder för att ta fram lösningar. I denna studie används finita differenser där den matematiska konstruktionen av en oändligt liten differens för att approximera en derivata ersätts av små men ändliga differenser. Randvillkor är viktiga inom den numeriska beräkningen då rummet, som beräkningarna utförs i, måste avgränsas. Dessa villkor bestämmer då hur ekvationen ser ut vid dessa ränder.

I denna rapport undersöks den tidsberoende Schrödingerekvationen och Gross-Pitaevskiiekvationen nu-meriskt i MATLAB med finita differenser för en rumsdimension. Särskilt den numeriska hanteringen av randvillkor behandlas där de två randhanteringsmetoderna SBP-SAT och SBP-Projektion undersöks. I denna undersökning används en diskontinuerlig potential, vilket gör att rummet måste delas upp i flera block vilket i sin tur skapar fler randvillkor att uppfylla. På grund av detta tillämpas ett interface med en inre rand-hantering som behandlar detta och återger en fysikalisk lösning. Både yttre och inre randvillkor undersöks med dels rena SBP-SAT- och SBP-Projektion-metoder men även med kombinationer av dessa. Stabilitet av dessa metoder säkerställs matematiskt för att garantera fysikaliska lösningar. Även konvergensen för olika randvillkor och kombinationer av dessa randhanteringsmetoder undersöks.

Contents

1 Introduction 1

1.1 Background . . . 1

1.2 Objective . . . 1

2 Theory 2 2.1 The Schrödinger Equation . . . 2

2.1.1 The GPE . . . 2

2.2 Well-Posedness and the Energy Method . . . 3

2.3 SBP-operators . . . 3

2.4 SAT . . . 4

2.5 The Projection Method . . . 4

2.6 Introducing a Potential . . . 5

2.6.1 Interfaces with SAT and Projection . . . 6

2.7 Discretization in Time . . . 6

2.8 Error and Convergence Calculation . . . 7

3 Results 8 3.1 Stability . . . 8

3.2 Convergence for 1D Schrödinger without Interface . . . 8

3.3 Convergence for 1D Schrödinger with Interface . . . 11

3.4 Bright Soliton Interaction with Barrier . . . 15

3.5 Dark Soliton Interaction with Barrier . . . 20

3.5.1 Simulations with Potential V = 0.07 . . . 22

3.5.2 Simulations with Potential V = 0.063670375 . . . 24

4 Discussion 27

5 Conclusion 29

References 30

1 Introduction

1.1 Background

Partial differential equations (PDE:s) arise in many applications in the world. One phenomenon that can be described by PDE:s is solitons. Solitons are a type of self-reinforcing solitary wave that were first observed on the ocean by J.S Russel in England in 1834 [1]. Solitons also appear in quantum mechanics as solutions to the non-linear Schrödinger equation. One special case of the non-linear Schrödinger equation is the Gross-Pitaevskii equation (GPE). The GPE describes a soliton that appears in a Bose-Einstein condensate (BEC), which is a state matter that can reach when experiencing extremely low temperatures. Except for BEC:s, the GPE is also used in applications related to nonlinear optics, modeling the propagation of light in nonlinear optical fibers.

If a PDE lacks an analytical solution numerical methods need to be applied in order to obtain a solu-tion. A good numerical method needs to be stable and have a high-order of accuracy. Higher order methods are well known for being able to solve PDE:s while maintaining a small error, thus being quite efficient for a given error tolerance compared to first- or second-order methods. Higher-order finite difference methods (HOFDM) are especially great in capturing transient phenomena, but it can be difficult to obtain a stable treatment of the boundaries.[2] To improve the stability of HOFDM, summation-by-parts (SBP) operators were developed and are used together with other methods to impose boundary conditions (BC:s) weakly. One boundary imposing method is the SAT-method, first developed by M. Carpenter et al. in 1994 [3], which is a well-proven finite difference method with high order of accuracy. The method imposes the BC:s while maintaining strict stability, albeit with some loss in the accuracy of the BC:s. This is a great method for coupling boundary conditions in multi-block problems with non-smooth geometries.[4] Another method to impose the BC:s is the Projection method, developed by P. Olsson in 1995 [5, 6]. This method projects the BC:s on a subspace and imposes them more precisely then the SAT-method. The improved Projection method eliminates the risk of instability and is the one to be used in this study. [2]

In this project, the GPE and the linear Schrödinger equation will be used to evaluate numerical meth-ods. In the presented PDE:s, a non-smooth potential will be applied, meaning the potential will take a discontinuous step from V (x) = 0 to some potential V (x) = V at some point x. To be able to handle this discontinuity without losing method accuracy, the space is split into blocks communicating with each other and fulfilling inner BC:s of continuity via an interface. Numerical methods that can handle this are thus required. In this report the numerical finite difference methods of SBP-SAT and SBP-Projection will be used, where the SBP-operators will be of both fourth and sixth order accuracy. The Projection method has yet not been implemented on the GPE and it is thus unexplored which method is the preferred one.

1.2 Objective

The purpose of this project is to implement a stable HOFDM in MATLAB for both the one-dimensional (1D) linear Schrödinger equation and the 1D GPE, utilizing summation by parts (SBP) operators and different methods to impose the BC:s. To achieve this, the well-posedness must be analyzed for both equations. The performance of four numerical methods will be compared: pure SBP-SAT, pure SBP-Projection and two hybrid versions of these.

2 Theory

2.1 The Schrödinger Equation

This study addresses two different Schrödinger equations: the 1D Schrödinger equation (2) and the non-linear Schrödinger equation (3). The latter is also referred to as the Gross Pitaevskii equation (GPE). The key difference is that the GPE has a non-linear term which makes the equation non-dispersive while the Schrödinger equation is dispersive. The non-linear term in the GPE will not affect the treatment of the boundary values numerically.

i_~ut=

⇣ _~2

2m + ˜V (x) ⌘

u (1)

The solutions of the Schrödinger equations are complex-valued. Due to Heisenberg’s uncertainty principle, the exact location of the particle cannot be measured exactly when the momentum is known. The proba-bility density function u(x, t)u⇤_{(x, t) shows the likelihood of finding the particle and is commonly used for}

visualization of solutions to the Schrödinger equation. In quantum physics, particles interacting with finite potentials will always result in a finite probability in finding the particle on either side of the potential. How these probabilities relate to each other depend on how the energy of the particle relates to the potential barrier.

In the 1D Schrödinger equation, ˜V (x) denotes the potential in the room the particle is situated in. For this study, the dimensionless Schrödinger equation (2) was studied.

iut= ↵uxx+ ˜V (x)u (2)

2.1.1 The GPE

When a quantum wave propagates it diffuses due to the inherent dispersiveness of the Schrödinger equation. The GPE, seen in (3), models a single particle wave function in a BEC. ~ is the reduced Planck constant, mis the mass of the particle and ˜V (x)the external potential. The waves described by the GPE keep their velocity and shape due to the non-linear term balancing the dispersion of the wave packet. Solitary waves with these characteristics are referred to as solitons. Solitons will maintain their identities while colliding with other solitons.[7]

The sign of the coupling constant, given in equation (4), will determine whether the soliton is considered a dark or a bright soliton, where the result is a dark soliton when ˜g > 0 and a bright soliton when ˜g < 0. Both bright and dark solitons are local disturbances in a BEC. Bright solitons have a concave shape in a zero background density, and dark solitons have a convex shape in a non-zero background density. Both types are considered in this study. asis the scattering length of two interacting bosons. The dimensionless GPE in 1D

is given in equation (5). It is difficult to solve the GPE (5) analytically due to the non-linear term g|u|2_u,

but analytical solutions have been derived for a particle in free space for both g > 0 and g < 0.[8]

2.2 Well-Posedness and the Energy Method

A well-posed problem has a unique solution that changes continuously with the initial conditions. The energy of a well posed continuous 1D problem is given by the norm ||u||2_{, where u, w 2 L}2_[x

l, xr]with respect to

the inner product of equation (6). The energy ||u||2 _{of a well posed continuous problem is limited in time}

and fulfills equation (7).

The energy method is a way of manipulating PDE:s into this form, here done by multiplying the PDE with its conjugate u⇤_{, integrating spatially by parts, and adding the Hermitian transpose. A well-posed}

problem with a bounded initial condition has non-increasing energy and is thus stable by definition. Apply-ing the energy method on the dimensionless SchrödApply-inger equation with an arbitrary potential ˜V (x) gives an energy estimate of (8), which is well-posed for both homogeneous Dirichlet and Neumann BC:s. Applying the energy method on the dimensionless GPE yields the same energy estimate since the energy terms from the g|u|2_{u-factor cancel each other.}

(u, w) = Z xr xl uw⇤dx (6) d dt||u|| 2  0 (7) d dt||u|| 2_{= i↵(u} xu⇤|xxrl u ⇤ xu|xxrl) (8)

The energy method can also be applied to discrete problems, which is discussed in the next section.

2.3 SBP-operators

SBP-operators of the finite difference type are known to provide stable and high-order accurate solutions to PDE:s. SBP-operators discretely mimic the behavior of integrating by parts by instead summating by parts. Thus, they are used in the semi-discrete energy method together with methods such as SAT or Projection to impose well-posed BC:s. After numerical well-posedness has been achieved for the semi-discrete case, a stable time discretization yields a stable and fully discrete problem. When applying the semi-discrete energy method the discrete norm of equation (9) is used.

||v||H = vTHv (9)

The domain x 2 [xl, xr]is discretized with m equidistant grid points. Each point is then given by xi= h(i 1),

where i = 1, 2, ..., m 1, m, resulting in a mesh-size h = xr xl

m 1 . The second derivative approximated by a

SBP-operator is given by equation (10), where H is the mass matrix, M is the stiffness matrix, ei is the unit

column vector for the point i, and di is a first derivative operator on row vector form for the point i. Note

that both M and H are symmetric positive definite matrices. SBP-operators are described in [10, 11] and the operators used in this project can be found in Appendix A.

D2= H 1( M e1d1+ emdm) (10)

2.4 SAT

SAT is used to impose BC:s weakly by adding penalty terms to the numerical PDE. A Dirichlet BC v(xi) = g(t), applied to a second derivative SPB-operator D2, is implemented with the penalty term

⌧iH 1dTi(eTi v g(t)). The corresponding Neumann BC vx(xi) = g(t) is implemented with the penalty

term ⌧iH 1ei(div g(t)). When all BC:s are applied this way the penalty factors ⌧i are decided such that

the method is well-posed using the semi-discrete energy method.

Applying the energy method to the semi-discrete dimensionless 1D Schrödinger equation (2) with homo-geneous Dirichlet BC:s will yield an energy estimate of (11). Here, ⌧1= i↵and ⌧m= i↵gives a well-posed

problem with constant energy. Applying the energy method to the same problem with homogeneous Neu-mann BC:s yields an energy estimate of (12) where ⌧1= i↵and ⌧m= i↵gives a well-posed problem with

no energy loss. Note that applying the semi-discrete energy method to the GPE yields the same energy estimate since the additional ig|v|2_v⇤_{v-term cancels with its conjugate.}

d dt||v||

2

H = v⇤e1d1v(⌧1⇤ i↵) + v⇤(e1d1)Tv(⌧1+ i↵) + v⇤emdmv(⌧m⇤ + i↵) + v⇤(emdm)Tv(⌧m i↵) (11)

d dt||v||

2

H = v⇤e1d1v(⌧1 i↵) + v⇤(e1d1)Tv(⌧1⇤+ i↵) + v⇤emdmv(⌧m+ i↵) + v⇤(emdm)Tv(⌧m⇤ i↵) (12)

The dimensionless Schrödinger equation in 1D with well-posed homogeneous Dirichlet BC:s is given in equation (13), both in continuous and semi-discrete form on the domain of x 2 [0, 1]. For homogeneous Neumann BC:s ux(0) = 0, ux(1) = 0, the iteration matrix A is replaced by A = ↵(D2 H 1e1d1+H 1emdm).

A fully discrete solution is then obtained by integrating the semi-discrete form numerically using Runge Kutta 4 (RK4). Note that neither the potential terms nor the non-linear term from the GPE contributes to the energy estimate. 8 > > > < > > > : ut= i↵uxx, u(0) = 0, u(1) = 0, x_{2 [0, 1], t 2 (0, T ]} 8 > > > < > > > : vt= i↵(D2v H 1dT1eT1v + H 1dTmeTmv) = iAv, e1v = 0, emv = 0, x = x1, x2, ..., xm 1, xm, t2 (0, T ] (13)

2.5 The Projection Method

The main idea with the Projection method is to impose the boundary conditions more precise than other BC imposing methods. Applying the energy method to a projected problem results in the discrete BC:s corre-sponding well to the continuous BC:s. The Projection method requires less algebraic manipulation than the SAT to make sure that the time derivative of the energy is strictly negative, which makes it less complicated to implement from a mathematical perspective. There are two different versions of the Projection method: the Traditional Projection method [5,6] and the improved Projection method [2]. The difference is that the improved Projection method removes numerical stability issues with the traditional Projection method by not introducing a zero eigenvalue. In this study, the improved Projection method is used (henceforth referred to as the Projection method).

The projection operator is defined in equation (14), where H represents the discrete norm called the Mass matrix, I is the identity matrix and L denotes the boundary operators. The derivative operator is multiplied with the projection operator on both sides; ergo, the D2 operator is replaced with A = P (↵D2+ V (x))P,

which is the iteration matrix.

When using the Projection method, L in the projection operator (14) is defined by equation (15) for ho-mogenous Dirichlet BC:s, and by (16) for homogeneous Neumann BC:s. Note that more conditions can be applied simultaneous by adding more rows to L.

L =  eT 1 eTm (15) L =  d1 dm (16)

The 1D Schrödinger equation approximated by the Projection method then becomes on the semi-discrete form shown in equation (17). RK4 will be used to acquire a fully discrete problem.

vt= iAv (17)

2.6 Introducing a Potential

Figure 1: A schematic figure over the step potential. The potential V (x) = V when x xv.

Solving the Schrödinger equation (2), with a step potential as defined by (18) and shown in Figure 1, with a single block containing both the regions with and without potential will lower the convergence rate significantly as shown in [4]. This can be circumvented by splitting the region without potential into one block and the region with potential into another block.

V (x) = (

V if xv x,

0else (18)

To create one global solution that covers both blocks and does not violate the Schrödinger equation, the blocks need to communicate and fulfill continuity in space and the first spatial derivative. These inner BC:s implying continuity are applied using an interface. The interface is implemented using a global matrix of the form of equation (19), where D(1)

2 represents a matrix containing the first part of the D2operator in block 1,

D₂(2) represents a matrix containing the latter part of the D2 operator in block 2 and C represents a matrix

containing the remaining parts of these D2 operators as well as the interface terms.

A =⇣ D (1) 2 C D(2) 2 ⌘ (19) The new global semi-discrete problem to be solved is given in equation (20) and results in a global solution vector when fully discretized with a stable numerical time integrator (RK4 in this study). The last point in v(1) _{and the first point in v}(2) _{represent the same spatial point. These points are equal when applying both}

the solution has been calculated, these points are replaced by a single point with an average value of the two points. ⇣ v(1) t v(2) t ⌘ = A⇣ v (1) v(2) ⌘ (20) This multi-block interface treatment can be repeated if more potential steps are added. After an interface treatment has been implemented, the potential V (x) can be added as diagonal elements in the global ma-trix. The same can be done for the non-linear g|u|2_{-term in the GPE. The matrix A is then redefined as}

A + diag(V ) + diag(g|u|2₎

2.6.1 Interfaces with SAT and Projection

The SBP-SAT approximation for the 1D Schrödinger equation, in a two-block interface with Dirichlet condi-tions, was acquired from [4] and can be seen in (21) and (22). v(1)_{denotes the Schrödinger equation without}

potential in block 1, and v(2) _{denotes the Schrödinger equation with potential in block 2.}

vt(1)= i↵(D2v(1) H 1dT1eT1v(1)+ 1 2H 1_dT m(eTmv(1) eT1v(2)) 1 2H 1_e m(dmv(1) d1v(2))) (21) v(2)t = i↵(D2v(2) V v(2)+ H 1dTmeTmv(2) 1 2H 1_dT 1(eT1v(2) eTmv(1)) + 1 2H 1_e 1(d1v(2) dmv(1))) (22) The terms H 1_dT

1eT1v(1)and H 1dTmeTmv(2)apply homogeneous Dirichlet conditions on the outer boundaries

and are omitted when using projection on the outer boundaries. All other terms, besides the D2-terms and

the potential term V v(2)_{, are interface terms. The global matrix constructed thus attains the form of equation}

(19).

For Projection, a global matrix A is created consisting of two D2 operators on the diagonal. The L term is

then constructed by concatenating the vectors applying the inner BC:s according to equation (23), where m denotes the number of points for the first block. The global matrix is then processed by multiplying by the P-operator on both sides, resulting in A being replaced by P AP , resulting in interface terms appearing and the global matrix taking the form of equation (19).

L =  eT m eT1 dm d1 (23)

2.7 Discretization in Time

Denote k as the time-step and A as the iteration matrix, then the time-discretizations for equations (2) and (5) are given by RK4 in equation (24). The fully discretized solution is given by v.

How effective a method is depends on what value of the time-step k in RK4 is required: a smaller k will result in more iterations and therefore longer computational time, but if k is too big the system will become unstable. The stability region of RK4, illustrated in Figure 2, can be used to determine the requirements of k. Since the eigenvalues of equation (13) are purely imaginary the stability condition for RK4 is given by |k · ⇢(A)|  2.8. Therefore, a time-step k  2.8

|⇢(A)| gives a stable system where ⇢(A) is the spectral radius

of matrix A and is defined by (25). The stability condition attained by selecting a time step based on the stability region of the method and the spectral radius is referred to as the CFL-condition in this study.

⇢(A) = max(|eig(A)|) (25)

The spectral radius is proportional to the inversed squared mesh-size, ⇢(A) = C · 1

h2, since A approximates a

second derivative. The largest possible time-step can be expressed as k = 2.8 C · h

2_{. Therefore, lower C results}

in higher efficiency.

Figure 2: Stability region for RK4. The eigenvalues need to be situated on or inside the blue line for the method to be stable.

2.8 Error and Convergence Calculation

Letting uN denote the numerical solution and uAthe analytical solution, the error l2can be calculated with

the discrete l2-norm given by equation (26). The convergence rate is given by q in equation (27), where

||uA umN2||h is the discrete l2-norm [4].

3 Results

Firstly, the results of the time discretization are demonstrated. The second subsection presents the results for the two methods applied to the outer boundaries to the 1D Schrödinger equation without any interface implemented. The results for the 1D Schrödinger equation with an implemented interface is presented in the third subsection. Here the Projection method is used at the outer boundaries, and a combination of methods is applied to the inner boundary. The results for both subsections consists of the error calculated with the analytical solution and the convergence rate. In the fourth subsection, the results for the study of soliton tunneling are presented. Tunneling of the Schrödinger equation is also included in this section for qualitative comparison to the soliton tunneling.

3.1 Stability

The constants C = ⇢(A) · h2 _{are demonstrated in Table 1. The Projection has lower values of C, implying}

that larger time steps k can be chosen than with SAT without jeopardizing stability. The constants C were calculated for the 1D Schrödinger equation (2) with ↵ = 1/2 and V (x) = 0.

Table 1: The constant C = ⇢(A) · h2 _{for the Projection and SAT method with Dirichlet and}

Neumann BC:s and SBP4- respective SBP6-operators.

Proj Dirichlet Proj Neumann SAT Dirichlet SAT Neumann

SBP 4 2.665 2.665 4.19 2.665

SBP 6 6.972 4.4196 7.7575 7.0897

3.2 Convergence for 1D Schrödinger without Interface

A convergence study was carried out for homogeneous Dirichlet and Neumann BC:s for the 1D dimensionless Schrödinger equation in (2), with ↵ = 1

2, x 2 [0, 1] and ˜V (x) = 0. The analytical solution is given in

equation (28). The initial condition is u0=p2sin(⇡x)for homogeneous Dirichlet BC:s and u0=p2cos(⇡x)

for homogeneous Neumann BC:s. The error was calculated using the discrete l2-norm given by equation (26)

and the convergence rate given by (27).

uexact= u0e

i⇡2 t

2 (28)

Table 2: SAT and Projection with Dirichlet BC:s and SBP4. logl2is the logarithm of the error

for each method and q is the convergence between the different points m for each method. m logl(SAT )₂ q(SAT ) _logl(P roj)

2 q(P roj) 41 -6.7563 0.00 -6.7602 0.00 81 -7.9477 3.9577 -7.9478 3.9451 161 -9.1258 3.9137 -9.1259 3.9136 321 -10.3289 3.9964 -10.3289 3.9963 641 -11.3795 3.4901 -11.3795 3.4899 1281 -10.5388 -2.7928 -10.5388 -2.7927

Table 3: SAT and Projection with Neumann BC:s and SBP4. logl2is the logarithm of the error

for each method and q is the convergence between the different points m for each method. m logl₂(SAT ) q(SAT ) _logl(P roj)

2 q(P roj) 41 -6.1318 0.00 -6.1662 0.00 81 -7.4551 4.3959 -7.4874 4.3892 161 -8.7861 4.4213 -8.8111 4.3972 321 -10.1029 4.3746 -10.1225 4.3563 641 -11.3329 4.0857 -11.3449 4.0608 1281 -10.5443 -2.6195 -10.54395 -2.6611

Table 4 states the results of the convergence study for SBP6 with Dirichlet BC:s. For SAT, the result is similar to earlier behavior: the convergence rate is around 5.5 before the error starts to increase again. The outcome for Projection is somewhat different and behaves strangely when m > 321. The numerical solution starts to oscillate between a big error and a small error, due to rounding errors, after the minimum error around logl2= 12is reached. In Table 5, the results for SBP6 with Neumann BC:s are shown. The methods

reach their highest convergence rate at fairly low resolved meshes. The error decreases until m = 321 and reaches a minimal error of 10 12_{for Projection and 10} 13_{for SAT. For greater m the error starts to increase}

again due to rounding errors.

Table 4: SAT and Projection with homogeneous Dirichlet BC:s and SBP6. logl2is the logarithm

of the error for each method and q is the convergence between the different points m for each method.

Table 5: SAT and Projection with Neumann BC:s and SBP6. logl2is the logarithm of the error

for each method and q is the convergence between the different points m for each method. m logl₂(SAT ) q(SAT ) _logl(P roj)

2 q(P roj) 41 -7.6804 0 -7.8197 0.00 81 -9.6644 6.5906 -9.7729 6.4885 161 -11.6639 6.6424 -11.7634 6.6123 321 -13.1796 5.0350 -12.4093 2.1457 641 -13.1112 -0.2272 -11.8776 -1.7664 1281 -12.7876 -1.0750 -11.36215 -1.7124

In Figures 3 and 4, the errors calculated by equation (26) are illustrated as a function of time. The plots show how effective a method is: showing the time it takes to reach a desired error tolerance. A finer mesh means the computational time will be longer. m increases from left to right. The error on the y-axis is defined as the l2-norm. For both SBP4 and SBP6 it can be seen that the Projection method is slightly more effective

than the SAT for each BC, especially for Dirichlet BC:s. When m increases the SAT method outperforms the Projection method by reaching a lower error tolerance for both BC:s. The error for m = 641 for the Projection method is not plotted for SBP6 since the minimal error has been reached and the rounding error dominates the total error.

Figure 3: Time required to reach desired l2-error tolerance for Projection and SAT method

Figure 4: Time required to reach desired l2-error tolerance for Projection and SAT method

with homogeneous Neumann and Dirichlet BC:s using SBP6 operators.

3.3 Convergence for 1D Schrödinger with Interface

A convergence study for the interface was carried out by sending a Gaussian wave packet through an interface constructed with two blocks. The wave packet is an analytical solution to the 1D Schrödinger equation in free space. Different combinations of the SAT and Projection methods were used at the inner boundaries, while the Projection method was used to imply homogeneous Dirichlet conditions at the outer boundary for both SBP4 and SBP6. The error was calculated by (26) and the convergence rate by (27). The dimensionless 1D Schrödinger equation (2) was used with ↵ = 1

2, x 2 [ 40, 80] with an interface at x = 20, ˜V (x) = 0 and

initial condition u0= 4

q

2 ⇡e

x2_+ik

0x. The analytical solution to this problem is given by (29). The numerical

simulation was run up to t = 5.

u = 4 r 2 ⇡ e k240 p 1 + 2ite 1 p_1+2it(x ik0₂ )2 (29) Figures 5 and 6 show the numerical solutions for both a low resolved and a high resolved spatial grid, us-ing pure Projection on the inner boundaries as well as the analytical solution. Observe that the numerical solution does not match the analytical solution well for a low resolution and that part of the solution is reflected even though no potential is present. This reflection is present even in the high resolved solutions but to a much lesser extent. Note that both method combinations of Projection and SAT, as well as the pure Projection and pure SAT, result in similar solutions with the same behavior.

In Table 6, the errors and the convergence rate for SBP4 with homogeneous Dirichlet BC:s together with an interface are presented. It can be seen that the convergence rates approach 4 for each combination of methods when m increases. As can be seen, the error only reaches a minimal value of 10 3.6_{. However,}

minimal value of 10 5.2_{. The convergence rate has started to decrease which could indicate that the minimal}

error could have been reached. Note that in both Tables 6 and 7 the method combination Method1M ethod2

refers to Method1 implementing continuity in the solution v and Method2 implementing continuity in the

first spatial derivative vx. The error is plotted against the calculation time in Figures 7 and 8. There is no

prominent behavior in any of the methods for either SBP4 or SPB6.

Figure 6: Solutions for the 1D Schrödinger equation solved with method SBP6 and pure Pro-jection. The interface is situated at x = 20. The solutions are plotted at the same time t = 5.

Table 6: Pure SAT, pure Projection and mixes. Homogeneous Dirichlet BC at the outer bound-ary and SBP4 are used. logl2is the logarithm of the error for each method and q is the

conver-gence rate between the different number of grid points m over both blocks. m logl(SAT SAT )₂ q(SAT SAT ) _logl(P rojP roj)

2 q(P rojP roj) logl

(SAT P roj)

2 q(SAT P roj) logl

Table 7: Pure SAT, pure Projection and mixes. Homogeneous Dirichlet BC at the outer bound-ary and SBP6 are used. logl2is the logarithm of the error for each method and q is the

conver-gence rate between the different number of grid points m over both blocks. m logl(SAT SAT )₂ q(SAT SAT ) _logl(P rojP roj)

2 q(P rojP roj) logl

(SAT P roj)

2 q(SAT P roj) logl

(P rojSAT ) 2 q(P rojSAT ) 402 -0.7477 0.00 -0.7534 0.00 -0.7477 0.00 -0.7534 0.00 802 -0.7032 -0.1477 -0.7147 -0.1285 -0.7070 -0.1351 -0.7108 -0.1415 1602 -1.7660 3.5304 -1.7877 3.5644 -1.7882 3.5916 -1.7655 3.5038 3202 -3.5619 5.9660 -3.5625 5.8957 -3.5638 5.8983 -3.5607 5.9633 4802 -4.6095 5.9491 -4.6091 5.9433 -4.6098 5.9401 -4.6088 5.9523 6402 -5.2084 4.7932 -5.2082 4.7956 -5.2084 4.7914 -5.2082 4.7973

Figure 7: Time required to reach desired l2-error tolerance for pure SAT, pure Projection and

Figure 8: Time required to reach desired l2-error tolerance for pure SAT, pure Projection and

mixes. Homogeneous Dirichlet BC and SBP6 are used.

3.4 Bright Soliton Interaction with Barrier

The dimensionless GPE (5) with g = 1 was considered. A three-block discretization was implemented on the domain x 2 [ 10, 32.5] with V (x) = 1

2(⌫ 2₊1

3)corresponding to half the kinetic energy of the soliton

with mass m = 2 at x 2 [10, 12.5]. A speed of ⌫ = 5 was used. The initial condition was given by the analytical solution to a bright soliton in free space in equation (30) at t = 0 (acquired from [8]). Projection was used to impose homogeneous Dirichlet BC:s at the outer boundaries, whilst pure Projection, pure SAT and combinations of them were used at the inner boundaries. For Figures 13-16, the method notation is the same as in Tables 6 and 7. The blocks consisted of 201, 26 and 201 points respectively. The time step used was half the maximum time step allowed by the CFL-condition.

u = sech(x ⌫t)ei(⌫ x 12(⌫ 2 ₁2_)t)

(30) For qualitative comparison, the 1D Schrödinger equation (2) was numerically simulated in the domain x 2 [ 0.1, 0.325] nm with potential V (x) = 1.1m⌫2

2 , x 2 [0.1, 0.125] nm up to time t = 6 · 10

18 _{s. The}

initial condition is given by the analytical solution (31) where ~ denotes Planck’s reduced constant, r0 is a

localization parameter set to 0.01 nm, p = m⌫ where m is the electron mass. ⌫ = 0.1c is the velocity where c denotes the speed of light. The constant N1 = ((2⇡)

1 2r₀)

1

2 is a normalization constant. The propagating

waves were plotted at times to facilitate the analysis of the qualitative similarities and differences between the waves propagating via the Schrödinger equation and the solitons propagating via the GPE.

u = N1e

x2

4r2_0ei1_~xp (31)

Figures 9-12 and 13-16 illustrate the probability density function |u(x, t)|2 _{over time, for the Schrödinger}

equation consists of several waves, whilst the GPE reflection only results in a single soliton. As for the 1D Schrödinger equation, the results produced are the same as in [4].

Figure 9: The 1D Schrödinger equation at time t = 1.5 · 10 18_{s, with a potential V = 1.1E at}

the interface (situated between the dotted lines), where E is the kinetic energy of the particle.

Figure 10: The 1D Schrödinger equation at time t = 3 · 10 18_{s, with a potential V = 1.1E at}

Figure 11: The 1D Schrödinger equation at time t = 4.5 · 10 18_{s, with a potential V = 1.1E at}

the interface (situated between the dotted lines), where E is the kinetic energy of the particle.

Figure 12: The 1D Schrödinger equation at time t = 6 · 10 18_{s, with a potential V = 1.1E at}

Figure 13: The GPE with g = 1 and a potential V = 0.5E between the dotted lines.

Figure 15: The GPE with g = 1 and a potential V = 0.5E between the dotted lines.

3.5 Dark Soliton Interaction with Barrier

The dimensionless GPE in (5) was considered with g = 1. The domain of x 2 [ 10, 32.5] was discretized with three blocks and a step potential located at x 2 [10, 12.5]. The analytical solution to a dark soliton in free space, given by equation (32), was used as initial conditions at t = 0 with v = k = 3⇡

42.5, a = 1 (acquired

from [9]). The potential was varied and numerically simulated with the simulation shown in Figures 17 - 29. Periodicity in the solution v and the first spatial derivative vxwas applied as outer boundary conditions with

Projection, and the projection term L given by equation (33), where the number of zeros corresponds to the number of grid points in the second block. The blocks consisted of 121, 16 and 121 points respectively. The time step used was half the maximum time step allowed by the CFL-condition.

u = a tanh a(x vt)ei(k(x vt)+(k22 a 2_)t) (32) L =  d10 0...0 0 dm eT 1 0 0...0 0 eTm (33)

Figures 17 - 20 show interaction with the potential V = 0.06, Figures 21 - 24 interaction with the potential V = 0.07 and Figures 25 - 29 interaction with the potential V = 0.063670375. For the Figures 17-29, the method notation is the same as in Tables 6 and 7. Note that the dark soliton reflects for V = 0.07, transmits for V = 0.06 and lingers inside the potential for V = 0.063670375 before transmitting. Some time intervals where the soliton is almost stationary excluded for brevity in the simulation with V = 0.063670375. The background density changes from a uniform background density to a non-uniform background density due to the step potential.

Figure 18: The GPE with g = 1 and a potential V = 0.06 between the dotted lines.

Figure 20: The GPE with g = 1 and a potential V = 0.06 between the dotted lines. 3.5.1 Simulations with Potential V = 0.07

Figure 22: The GPE with g = 1 and a potential V = 0.07 between the dotted lines.

Figure 24: The GPE with g = 1 and a potential V = 0.07 between the dotted lines. 3.5.2 Simulations with Potential V = 0.063670375

Figure 26: The GPE with g = 1 and a potential V = 0.063670375 between the dotted lines.

Figure 28: The GPE with g = 1 and a potential V = 0.063670375 between the dotted lines.

4 Discussion

Regarding the outer boundaries, there seems to be no difference in the convergence performance of SBP-Projection and SBP-SAT for either SBP4 or SBP6 with homogeneous boundary conditions. When looking at the efficiency plots for the outer boundaries in Figures 3 and 4, it becomes clear that there is no significant difference between the methods except for when applying homogeneous Dirichlet BC:s where SBP-Projection performs slightly better than SBP-SAT. This difference can be explained when looking at the spectral radius, in Table 1, where SAT Dirichlet has a notably higher spectral radius. This explains the difference since the time steps used in the efficiency plots of Figures 3 and 4 were half the maximal allowed time step by the CFL-condition. For computations on a larger scale, the time difference to reach a certain error tolerance may become significant and for that reason, SBP-Projection may be preferred in that case. The spectral radius also had an impact on the choice of method for the outer BC:s in the interface convergence study. The Projection method was used instead of SAT on the outer boundaries since the SAT on the outer boundaries dominated the spectral radii when using homogeneous Dirichlet BC:s.

When comparing the different orders of SBP-operators, SBP6-operators are undoubtedly better than SBP4, which was expected. The operators are both more effective and accurate. Both SBP4- and SBP6-operators were used in this study only to conclude that the convergence behaviors for different methods were not exclusive for SBP6. Higher orders of SBP-operators exist and could be considered for future studies. Note that the minimal error obtained with SBP6-SAT is lower than the error obtained with SBP6-Projection in Tables 4 and 5. This might suggest that the minimum possible error to obtain with SBP6-SAT is lower than the minimum possible error with SPB6-Projection. Although this could be the case, further studies with more grid points would need to be carried out to exclude the possibility that there might exist error values of the same magnitude for fewer grid points.

There was a significant difference in how low values the error could reach for the 1D Schrödinger with an interface compared to without it. Without the interface, the error reached a magnitude of 10 11 _or

smaller before the rounding error started to dominate and the convergence rate changed for the worse. With an interface, the error only reached a magnitude of 10 5 _{before the convergence rate started to level out.}

This does not necessarily mean that the methods cannot reach a low error tolerance for an interface, and such a conclusion cannot be drawn without first running the calculations for more amount of grid points. The domain of the simulation with an interface (x 2 [ 40, 80]) differed from the domain of the simulation without an interface (x 2 [0, 1]) which means that the mesh sizes were higher resolved without an interface and lower resolved with an interface. This could have a numerical impact and the calculations should be run again with an adjusted mesh size. Additionally, different initial conditions were used for the two cases which means error comparisons between the two cases could be unjustified. Finally, the Gaussian wave package used is a solution to the 1D dimensionless Schrödinger equation in free space and the domain may have been insufficiently small as the solutions, in Figures 5 and 6, are close to the right boundary. The initial conditions and the BC:s could thus have impacted the minimal achievable error.

As presented in Table 1, the eigenvalues in the iteration matrix associated with the Projection method are lower than the eigenvalues for the SAT method. This results in an opportunity to choose a greater time-step for the Projection method while still maintaining stability. This leads to fewer iterations required, thus making the method faster. This is shown in Figures 3 and 4, where the Projection method reaches the error tolerance faster. On the contrary, SAT reaches a lower error which can be seen in Figures 3 and 4 and Table 2-3. The minimal possible error obtainable with SBP-SAT and SBP-Projection might still be the same and further studies would be necessary to conclude that one method reaches a lower error than another as argued in the first paragraph of the discussion. The mentioned differences are small, hence it is difficult to establish the preferred method. Nevertheless, the Projection method is slightly easier to implement from a mathematical perspective since no penalty parameters ⌧i have to be chosen.

The difference between the behavior regarding reflection/transmission of the wave packet propagating, ei-ther with the Schrödinger equation or the bright soliton propagating with the dimensionless GPE, seen in Figures 9-16, can be explained by the fact that the Schrödinger equation is dispersive while the GPE is non-dispersive. The dispersion of the Schrödinger equation means that solution will spread out over time, while the soliton maintains its shape. This explains why the reflection of the Schrödinger equation consists of several waves and why the GPE reflection and transmission results in single solitons. When running the simulations it can be seen that another soliton appears on the other side of the barrier. The parameters for the GPE were chosen without physical consideration, which means that the validity of the physical prop-erties has not been taken into account. Physically correct parameters may yield another result such as a totally reflected or totally transmitted soliton. Nevertheless, the results show that both SAT, SBP-Projection and combinations can be used to simulate solitons with a step-potential. Note that this is the first time to the author’s knowledge that the Projection method has been successfully implemented with the GPE. The simulations of dark solitons created with g = 1 in the GPE show that the interaction with a po-tential barrier for a dark soliton differs from the bright soliton case. Notably, the dark soliton does not seem to split upon interaction with the potential barrier as the bright soliton did. Instead, the entire dark soliton moves through the potential barrier, getting slowed while doing so, and then either totally reflecting or transmitting. This is shown in Figures 17 - 29.

5 Conclusion

References

[1] Russell, J.S. Report on Waves. Report of the fourteenth meeting of the British Association for the Ad-vancement of Science, York. September 1844; p 311-390. Plates XLVII-LVII

[2] Mattsson K, Olsson P. An improved projection method. Journal of Computational Physics. November 2018; Vol 372: p. 349 – 372.

[3] Carpenter M.H, Gottlieb D, Abarbanel S. The stability of numerical boundary treatments for compact high-order finite-difference schemes. Journal of Computational Physics. October 1993; Vol 108(2): p. 272-295

[4] Almquist M, Mattsson K, Edvinsson T. High-fidelity numerical solution of the time-dependent Dirac equation. J. Comput. Phys. April 2014; Vol 262: p. 86–103.

[5] Olsson P. Summation by parts, projections and stability I. Mathematics of Computation. July 1995; Vol 64 No 211: p. 1035-1065

[6] Olsson P. Summation by parts, projections and stability II. Mathematics of Computation. October 1995; Vol 64 No 212: p. 1473-1493

[7] Wadati M. Introduction to solitons. Pramana - Journal of Physics. November 2001; Vol 57: p. 841–847 [8] Bao W, Tang Q, Xu Z. Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation. Journal of Computational Physics. February 2013; Vol 235: p. 423 – 445.

[9] Sakaguchi H, Tamura H. Scattering of Solitons and Dark Solitons by Potential Walls in the Nonlinear Schrödinger Equation. Journal of the Physical Society of Japan. January 2005; Vol. 74: p. 292–298 [10] Mattsson K, Nordström J. Summation by parts operators for finite difference approximations of second

derivatives. J. Comput. Phys. September 2004; Vol 199(2): p. 503–540.

Appendix A

SBP-operators

All SBP-matrices are of size m ⇥ m. SBP4

The Mass matrix H is given by

H = h · 0 B B B B B B B B B B B B B B B B B B B B B @ 17/48 0 . . . 0 0 59/48 ... ... ... ... 43/48 49/48 1 ... 1 49/48 43/48 ... ... ... ... 59/48 0 0 . . . 0 17/48 1 C C C C C C C C C C C C C C C C C C C C C A

The first derivative difference operators d1 and dm are defined as

d1=_h1 [ ₆11 3 3₂ 1₃ 0 ... 0]

dm= 1_h [0 ... 0 ₃1 3₂ 3 11₆]

SBP6

The Mass matrix H is defined by H1,1 =13649₄₃₂₀₀, H2,2 =12013₈₆₄₀, H3,3=2711₄₃₂₀ H4,4 =5359₄₃₂₀, H5,5=7877₈₆₄₀, H6,6= 43801₄₃₂₀₀ Hm 5,m 5 =43801₄₃₂₀₀, Hm 4,m 4 =7877₈₆₄₀, Hm 3,m 3 =5359₄₃₂₀ Hm 2,m 2 =2711₄₃₂₀, Hm 1,m 1 =12013₈₆₄₀, Hm,m= 13649₄₃₂₀₀ H = h · 0 B B B B B B B B B B B B B B B B B @ H1,1 0 . . . 0 0 ... ... ... ... ... H6,6 1 ... 1 Hm 5,m 5 ... ... ... ... ... 0 0 . . . 0 Hm,m 1 C C C C C C C C C C C C C C C C C A

The first derivative difference operators d1 and dm are defined as

d1=_h1 [ ₁₂25 4 3 4_{3 4}1 0 ... 0]