![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00671v1-figure-ed7d9866-b418-46ba-9157-9cbcdb6f65a9.jpg)

## Absorbing Boundary Conditions for Variable Potential Schrödinger Equations via Titchmarsh-Weyl Theory

Matthias Ehrhardt , Chunxiong Zheng ∗ †

IMACM, School of Mathematics and Natural Sciences, University of Wuppertal, Germany

Department of Mathematical Sciences, Tsinghua University, Beijing 100084, P.R. China

August 2, 2024

## Abstract

We propose a novel approach to simulate the solution of the time-dependent Schrödinger equation with a general variable potential. The key idea is to approximate the Titchmarsh-Weyl m-function (exact Dirichlet-to-Neumann operator) by a rational function with respect to an appropriate spectral parameter. By using this method, we overcome the usual high-frequency restriction associated with absorbing boundary conditions in general variable potential problems. The resulting fast computational algorithm for absorbing boundary conditions ensures accuracy over the entire frequency range.

AMS classification: 65M99, 81-08

Keywords: absorbing boundary conditions, variable potential, Schrödinger equation, Titchmarsh-Weyl m-function, unbounded domain

## 1 Introduction

In this paper we consider the linear Schrödinger problem of the form

$$\mathrm i u _ { t } + \partial _ { x } ^ { 2 } u & = V ( x ) u, \quad ( x, t ) \in \mathbb { R } \times ( 0, T ], \\ u ( x, 0 ) & = u _ { 0 } ( x ), \quad x \in \mathbb { R },$$

where T denotes the finite evolution time, and u 0 is an initial wave packet supported in a finite interval Ω int = [ x -, x + ] with x -&lt; x + . It is well known that under mild conditions the Cauchy problem (1) has a unique solution u ∈ C ( R + , L 2 ( R )) , cf. [33], e.g.:

Theorem 1. Let u 0 ∈ L 2 ( R ) and real-valued potential V ∈ L ∞ ( R ) . Then the problem (1) has a unique solution u ∈ C ( R + , L 2 ( R )) . Moreover, the 'energy' is preserved, i.e.

$$\| u (., t ) \| _ { L ^ { 2 } ( \mathbb { R } ) } = \| u _ { 0 } \| _ { L ^ { 2 } ( \mathbb { R } ) }, \quad \forall \, t \geq 0.$$

∗ Corresponding author, ehrhardt@uni-wuppertal.de

† Supported by the Alexander-von-Humboldt Foundation, czheng@tsinghua.edu.cn

The Schrödinger problem (1) is defined on an unbounded domain x ∈ R . To numerically simulate its solution, it is common practice to truncate the domain to a bounded one, for example, Ω int ,T = Ω int × (0 , T ] . Absorbing boundary conditions (ABCs) are thus necessary for well-posedness at the two artificially introduced boundaries, Σ ± ,T = { x ± } × (0 , T ] .

Numerical simulation of the linear Schrödinger equation on unbounded domains with an external potential has been a hot research area for nearly thirty years, cf. the concise review article [7]. An ABC is called exact if the solution of the truncated domain problem remains the same as that of the original unbounded domain problem. The exact ABC is guaranteed to exist due to the wellposedness of the linear Schrödinger problem (1), but it can only be formulated analytically for some special potentials, such as constant potential [18], linear potential [19], symmetric periodic potential [21], isotropic free particle potential, Morse potential, harmonic potential, and Bargeman potential, cf. e.g. [34]. In the more general case, i.e., for general variable potential problems one is led to design approximate analytical ABCs for a given frequency regime with respect to some a priori criterion. Methods in this category include the pseudo-differential calculus method [5, 6, 8], the perfectly matched layer (PML) method [42], and the operator splitting method [40]. To the authors' knowledge, all of them are essentially based on the high frequency approximations . For low-frequency problems, the ABCs would be less accurate by these methods.

This paper proposes a new approach to the design of ABCs for the Schrödinger problem. Inspired by the work of Alpert, Greengard, and Hagstrom [2] on the fast evaluation of nonreflecting boundary kernels for time-domain wave propagation, we approximate the Titchmarsh-Weyl m-function (equivalently, the exact DtN operator) in the frequency domain by a rational function with respect to an appropriate spectral parameter. In the time domain, the nonreflecting boundary kernels are thus approximated by a sum of exponentials, which makes the approximate ABCs easy to implement.

The rationality of the above treatment is due to the analyticity property and the asymptotic behavior of the m-function. Since our approximation is performed in the whole frequency regime, the proposed ABCs are expected to be more versatile and accurate, especially in the low-frequency regime, thus overcoming the typical high-frequency restriction. Note that the Titchmarsh-Weyl mfunction is nothing else but the so-called total symbol in microdifferential calculus, which is treated by an asymptotic expansion to obtain a hierarchy of ABCs, cf. [5, 6, 8]. Note also that Titchmarsh-Weyl theory is used in the analysis of initial value problems for Schrödinger equations operator-valued potentials [25] or strongly singular potentials [30]. It is also used in practical applications in quantum mechanics [15, 26] and in option pricing in mathematical finance [31].

This paper is organized as follows. In Section 2 we review the basic facts of the Titchmarsh-Weyl theory for Schrödinger operators in one dimension for ease of later reference. Then, in Section 3 we discuss the Titchmarsh-Weyl m-function (i.e. the exact Dirichlet-to-Neumann operator) and explain the algorithm used to compute the m-function numerically. Thus, at least from a numerical point of view, the exact ABC is explicitly known, see Section 3. However, when simulating the Schrödinger equation (1), the difficulty does not lie in the computation of the m-function in a frequency domain method, which is presented in Section 4, but in its inverse Laplace transformation, which is too expensive. For this reason, in Section 5 we introduce a rational approximation of the m-function in the frequency domain to obtain an approximate ABC that can be computed efficiently using a fast evaluation technique [41] in the time domain. We discuss practical implementation issues and finally in Section 6 we conclude with numerical results illustrating that our new approach leads to an efficient and reliable algorithm for the time-dependent Schrödinger equation with a general variable potential.

## 2 The Titchmarsh-Weyl theory

We will review here, for convenience of later work, the essentials of Titchmarsh-Weyl (TW) theory for Schrödinger operators in one dimension. The interested reader may consult [24, Section 2] or [12, Section 4.3] for a more detailed presentation.

First, consider the Schrödinger operator L on the real line given by

$$L = - \partial _ { x } ^ { 2 } + V ( x ), \quad x \in \mathbb { R },$$

with a real-valued, locally integrable potential V ∈ L 1 loc ( R ) . Let x 0 ∈ R be an arbitrarily chosen point, called reference point . In the following we will study how the solutions depend on this parameter x 0 .

To do this, we consider θ x x , λ ( ; 0 ) and ϕ x x , λ ( ; 0 ) to be the fundamental solutions of the Schrödinger eigenvalue problem.

$$- u _ { x x } + V ( x ) u = \lambda u, \quad x \in \mathbb { R }, \quad \lambda \in \mathbb { C },$$

with the following initial conditions at the reference point x 0 :

$$\theta ( x _ { 0 } ; x _ { 0 }, \lambda ) = 1, \quad \theta _ { x } ( x _ { 0 } ; x _ { 0 }, \lambda ) = 0,$$

$$\phi ( x _ { 0 } ; x _ { 0 }, \lambda ) = 0, \quad \phi _ { x } ( x _ { 0 } ; x _ { 0 }, \lambda ) = 1.$$

It can be shown that under these assumptions θ x x , λ ( ; 0 ) and ϕ x x , λ ( ; 0 ) exist on the whole real axis, and they are entire functions of λ and real for λ ∈ R . Now, as a basic fact of TW theory, the equation (4) has at least one solution ψ ± , called Weyl's solution with

$$\psi _ { \pm } ( x _ { 0 } ; x _ { 0 }, \lambda ) = 1,$$

and

$$\psi _ { \pm } ( x ; x _ { 0 }, \lambda ) \in L ^ { 2 } ( \mathbb { R } _ { \pm } ^ { x _ { 0 } } )$$

for any λ ∈ C + . Here R x 0 ± stands for the interval [ x , 0 ±∞ ) and C + denotes the upper half of the complex plane, i.e. C + = { z ∈ C | Im z &gt; 0 } . A potential V ( x ) is said to be in the limit-point case at ±∞ if and only if there exists only one Weyl's solution in the corresponding L 2 space. The reader will immediately realize that the assumption of V ( x ) in the limit-point case is necessary for the well-posedness of the Schrödinger problem (1) in a more general setting. At positive infinity, a standard sufficient condition for the limit-point case is given by [35]:

Theorem 2 ([35, Theorem X.8]) . Let V ( x ) be a continuous real-valued function on ( x , 0 ∞ ) and suppose that there exists a positive differentiable function M x ( ) such that

$$\left ( i \right ) \, V ( x ) \geq - M ( x ) \ i f \, x > x _ { 0 } ;$$

$$\left ( i i \right ) \underset { x _ { 1 } } { \overset { \infty } { \int } } \left ( M ( x ) \right ) ^ { - 1 / 2 } d x = \infty \text{ for any } x _ { 1 } > x _ { 0 } ;$$

(iii) M x / M x ′ ( ) ( ( ) ) 3 2 / is bounded near ∞ .

Then V ( x ) is in the limit-point case at ∞ .

An analogous result can be given at negative infinity point.

According to this theorem, a potential V ( x ) is in the limit-point case provided that V ( x ) ≥ -kx 2 for some constant k and for all sufficiently large x . This implies that the restriction on the potential for the limit-point case is very weak: It only excludes some especially strange potentials, which may not be physically relevant at all. Roughly speaking, the limit case does not admit potentials that tend too fast (faster than quadratically) to -∞ for x →±∞ .

Due to the boundary conditions (6a) we can write

$$\psi _ { \pm } ( x ; x _ { 0 }, \lambda ) = \theta ( x ; x _ { 0 }, \lambda ) + m _ { \pm } ( x _ { 0 }, \lambda ) \, \phi ( x ; x _ { 0 }, \lambda ),$$

with some uniquely determined coefficient, the Titchmarsh-Weyl m-function m ± ( x , λ 0 ) . This function plays a fundamental role in the spectral theory of the Schrödinger operator (3) on the half-line R x 0 ± .

We will now summarize some of the most important properties of the Titchmarsh-Weyl mfunction. First,

$$m _ { \pm } ( x _ { 0 }, \lambda ) \text{ is analytic with respect to } \lambda \text{ to } \mathbb { C } \mathbb { \Lambda } \text{ and } m _ { \pm } \colon \mathbb { C } _ { + } \rightarrow \mathbb { C } _ { + }$$

and is therefore called a Herglotz function (or Nevanlinna or Pick function), cf. [24, Lemma 2.3]. It is easy to show that this Herglotz-property is directly related to the positive type of the DtN-map in the sense of memory equations, cf. [18] for the corresponding case of constant external potential. Thus it is an essential ingredient of the stability w.r.t. the L 2 norm.

We also have the symmetry property

$$\overline { m _ { \pm } ( x _ { 0 }, \lambda ) } = m _ { \pm } ( x _ { 0 }, \bar { \lambda } )$$

and the local singularities of m are real and and most of them are first order, i.e.

$$\lim _ { \epsilon \to 0 + } ( - i \epsilon ) \, m _ { \pm } ( x _ { 0 }, \lambda + i \epsilon ) \geq 0, \quad \lambda \in \mathbb { R },$$

cf. [24, Theorem A.2].

Another important property is given by the Borg-Marchenko theorem [13, 32], which states that the Titchmarsh-Weyl m-function m ± ( x , λ 0 ) uniquely determines the potential V ( x ) at x &gt; x 0 (or x &lt; x 0 ). Moreover, since ψ ± ( x x , λ ; 0 ) changes with a simple multiplication when the reference point x 0 changes, one has

$$m _ { \pm } ( x, \lambda ) = \frac { \partial _ { x } \psi _ { \pm } ( x ; x _ { 0 }, \lambda ) } { \psi _ { \pm } ( x ; x _ { 0 }, \lambda ) }.$$

Thus, it is easy to verify that the m function satisfies the following Riccati equation:

$$\partial _ { x } m _ { \pm } ( x, \lambda ) = - m _ { \pm } ^ { 2 } ( x, \lambda ) + V ( x ) - \lambda.$$

## 3 The exact ABC by Titchmarsh-Weyl theory

We apply the Laplace transform

$$\hat { u } ( x, s ) = \mathcal { L } ( u ( x, t ) ) ( s ) = \int _ { 0 } ^ { + \infty } u ( x, t ) \, e ^ { - s t } \, d t, \quad \text{Re } s > 0,$$

to the Schrödinger equation (1) on the right exterior domain Ω + = { x ∈ R | x &gt; x + } and on the left exterior domain Ω -= { x ∈ R | x &lt; x -} . In the frequency domain, the Schrödinger equation is a second order homogeneous ODE

$$- \hat { u } _ { x x } + V ( x ) \hat { u } = \lambda \hat { u }, \quad x \in \Omega _ { \pm },$$

with λ = i s ∈ C + . The exact absorbing boundary condition of the DtN form in the frequency domain is thus

$$\hat { u } _ { x } ( x _ { \pm }, \lambda ) = m _ { \pm } ( x _ { \pm }, \lambda ) \hat { u } ( x _ { \pm }, \lambda ).$$

Only in some special cases, the m-function has a closed analytic form [14, 23]. For example, in the case of a constant potential V ≡ V 0 one gets

$$m _ { + } ( x _ { + }, \lambda ) = - \stackrel { + } { \sqrt { - \lambda + V _ { 0 } } }.$$

If the potential represents a harmonic oscillator , i.e. V ( x ) = x 2 on the interval [0 , ∞ ) , one obtains a meromorphic m-function given by the ratio of two gamma functions:

$$m _ { + } ( 0, \lambda ) = - \frac { 2 \Gamma ( \frac { 3 } { 4 } - \frac { 1 } { 4 } \lambda ) } { \Gamma ( \frac { 1 } { 4 } - \frac { 1 } { 4 } \lambda ) }.$$

Finally, for the Bargmann potential

$$V ( x ) = - 8 \beta ^ { 2 } \, \frac { \beta - \gamma } { \beta + \gamma } \frac { e ^ { - 2 \beta x } } { ( 1 + \frac { \beta - \gamma } { \beta + \gamma } e ^ { - 2 \beta x } ) ^ { 2 } }, \quad \beta > 0, \quad \gamma \geq 0,$$

one obtains the m-function

$$\i \i o n \quad & m _ { + } ( 0, \lambda ) = - \, \sqrt { - \lambda } - \frac { \gamma ^ { 2 } - \beta ^ { 2 } } { \sqrt { - \lambda } + \gamma }.$$

In the general case, however, numerical methods must be considered. This problem has been studied in many papers, e.g. [14, 39, 27, 29]. In this paper we simply compute the m -function by evolving the Riccati equation (12) with the classical fourth-order Runge-Kutta scheme and setting an initial data m ± ( x ± ,λ , λ ) = ∓ + √ -λ at a sufficiently distant point x ± ,λ = ± 200 . This treatment is reasonable since the potentials in our numerical tests actually decay to zero for x →∞ .

## 4 The frequency-domain method

The solution of the time-dependent Schrödinger equation could then be computed using the following frequency-domain method:

Step 1. Fix σ &gt; 0 . For each s = σ + i µ with µ ∈ R , solve the Laplace-transformed Schrödinger equation in the bounded interval [ x -, x + ] :

$$- \hat { u } _ { x x } + V ( x ) \hat { u } & = \text{is} \, \hat { u } - \text{i} u _ { 0 } ( x ), \quad x \in [ x _ { - }, x _ { + } ], \\ \hat { u } _ { x } ( x _ { - } ) & = m _ { - } ( x _ { - }, \text{is} ) \, \hat { u } ( x _ { - } ), \\ \hat { u } _ { x } ( x _ { + } ) & = m _ { - } ( x _ { + }, \text{is} ) \, \hat { u } ( x _ { + } ).$$

Step 2. Perform the inverse Laplace transformation

$$u ( x, t ) & = \mathcal { L } ^ { - 1 } ( \hat { u } ( s, t ) ) ( x ) \\ & = \frac { 1 } { 2 \pi i } \int _ { \sigma - i \infty } ^ { \sigma + i \infty } e ^ { s t } \hat { u } ( s, t ) \, d s = \frac { e ^ { \sigma t } } { 2 \pi } \int _ { - \infty } ^ { \infty } e ^ { i f t } u ( x, \sigma + i f ) \, d f,$$

to derive the wave function u x, t ( ) for any t ∈ (0 , T ] .

In the numerical implementation, some parameters have to be tuned. The function ˆ u is smoother for larger damping factor σ , but the evolution time span is then limited because an exponential term is involved in the inverse Laplace transformation. As a common practice, we set σ = 1 /T , where T is a prescribed evolution time. The integration domain is unbounded in (19) and must be truncated. We introduce a cut-off frequency f c and confine the integration to the interval [ -f , f c c ] . In addition, to get rid of high frequency oscillations in the inverse-transformed function, we should introduce another filtering function χ , which remains 1 over a sufficiently large frequency band with zero frequency as its center, and vanishes smoothly near the endpoints of [ -f , f c c ] . A good candidate (empirically) is

$$\chi = \exp \left ( - ( 1. 2 f / f _ { c } ) ^ { 2 0 } \right ).$$

After these treatments, we then derive an approximate inverse transformation as

$$\frac { e ^ { \sigma t } } { 2 \pi } \int _ { - \infty } ^ { \infty } e ^ { i f t } u ( x, \sigma + i f ) \, d f \approx \frac { e ^ { \sigma t } } { 2 \pi } \int _ { - f _ { c } } ^ { f _ { c } } \chi ( f ) \, e ^ { i f t } u ( x, \sigma + i f ) \, d f.$$

The right side is computed with an appropriate quadrature scheme.

## 5 The time-domain method

It follows from Section 2 that the exact ABC we are looking for is now explicitly known, at least from a numerical point of view. But this is not the whole story for simulating the solution of the time-dependent Schrödinger equation. The difficulty lies not in computing the m-function itself, but in computing its inverse Laplace transformation. Of course, a numerical inverse transformation is possible, but it would be too expensive.

Therefore, in this section we design an approximate ABC based on the rational approximation of the m-function. The kernel functions are of exponential type with respect to the half-order time derivative operator, so the fast evaluation technique proposed in [41] (cf. Appendix) is applicable. For some alternative fast evaluation methods, we refer the reader to [7] and the references therein. The rational approximation is realized by solving a least squares problem, an analogous technique to that used in [2] for fast evaluation of the boundary kernel functions of the hyperbolic wave equation.

In the time domain, the truncated Schrödinger problem reads

$$\mathrm i u _ { t } + \partial _ { x } ^ { 2 } u & = V ( x ) u, \quad ( x, t ) \in [ x _ { - }, x _ { + } ] \times ( 0, T ], \\ u ( x, 0 ) & = u _ { 0 } ( x ), \quad x \in [ x _ { - }, x _ { + } ], \\ u _ { x } ( x _ { \pm }, t ) & = \mathcal { L } ^ { - 1 } ( m _ { \pm } ( x _ { \pm }, \mathrm i s ) \, \hat { u } ( x _ { \pm }, s ) ) ( t ), \quad t \in ( 0, T ].$$

To simplify the notation, we will focus on the right boundary at x = x + . Let us recall that the DtN map in the frequency domain is

$$\hat { u } _ { x } ( x _ { + }, s ) = m _ { + } ( x _ { + }, \mathrm i s ) \, \hat { u } ( x _ { + }, s ).$$

Returning to the time domain we have to consider the convolution

$$u _ { x } ( x _ { + }, t ) = K ( t ) * u ( x _ { + }, t ), \ \text{ with } \ K ( t ) = \mathcal { L } ^ { - 1 } ( m _ { + } ( x _ { + }, i s ) ) ( t ).$$

There are two major difficulties here. First, it is generally hard to compute K( ) t and second, the convolution involved naturally leads to a nonlocal-in-time DtN map .

To get an idea, let us first consider two specific simple examples. In the case of the free Schrödinger ( V ≡ 0 ), cf. (15), we have

$$m _ { + } ( i s ) = - \sqrt { - i s } \quad \text{and thus} \quad \mathcal { L } ^ { - 1 } ( m _ { + } ( i s ) ) = - e ^ { - i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } },$$

with the half-order time derivative defined as

$$\partial _ { t } ^ { \frac { 1 } { 2 } } v ( t ) = \frac { 1 } { \sqrt { \pi } } \, \frac { \mathrm d } { \mathrm d t } \int _ { 0 } ^ { t } \frac { v ( \tau ) } { \sqrt { t - \tau } } \, d \tau,$$

which can be efficiently evaluated with some existing methods, e.g. [9, 28, 41]. Second, we consider the Bargmann potential (17), and in this case the m-function is

$$m _ { + } ( 0, i s ) = - \, \dagger \sqrt { - \text{i} s } - \frac { \gamma ^ { 2 } - \beta ^ { 2 } } { \sqrt { - \text{i} s } + \gamma },$$

thus we have

$$\i a v e { \mathcal { L } ^ { - 1 } ( m _ { + } ( 0, i s ) ) } = - e ^ { - i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } - ( \gamma ^ { 2 } - \beta ^ { 2 } ) ( e ^ { - i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } + \gamma ) ^ { - 1 }.$$

This operator can then be evaluated efficiently by introducing an unknown function and using the fast methods for ∂ 1 2 t .

Inspired by these two examples it naturally leads us to think about the possibility of approximating the m -function with a rational function with respect to a new spectral parameter k = + √ -i s (NOT s ), i.e.,

$$, \, \dots, \\ m _ { + } ( x _ { + }, \mathbf i s ) & \approx \tilde { m } _ { + } ( x _ { + }, \mathbf i s ) = - \, \sqrt { - \mathbf i s } + \sum _ { n = 1 } ^ { d } \frac { \alpha _ { n } } { \sqrt { - \mathbf i s } + \beta _ { n } }.$$

Once this is done, we can then replace the exact m -function with the approximate alternative ˜ m , which leads to the approximate kernel function

$$\mathcal { L } ^ { - 1 } ( \tilde { m } _ { + } ( x _ { + }, \mathbf i s ) ) = - e ^ { - \mathbf i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } + \sum _ { n = 1 } ^ { d } \alpha _ { n } ( e ^ { - \mathbf i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } + \beta _ { n } ) ^ { - 1 }.$$

The analogous idea for the Bargmann potential (17) can then be used to handle this kernel function. The answer for the possibility is affirmative considering that the asymptotic expansion has been given in [17] as √ √

$$m _ { + } ( x _ { + }, \lambda ) = - \stackrel { + } { \sqrt { - \lambda } } + \circ ( 1 / \sqrt { r } ), \ \ r \to \infty,$$

where λ = µr , r ∈ R and the convergence is uniform for µ in any compact subset of C + , cf. [36, Theorem C.4].

Now putting

$$g _ { + } ( \lambda ) = m _ { + } ( x _ { + }, \lambda ) + \, \dagger \sqrt { - \lambda },$$

in view of (25) we know that g + ( λ ) is analytic in C + with respect to λ and it tends to zero for λ → ∞ . We then use the method of Alpert, Greengard and Hagstrom [2] to approximate g + ( λ )

with a rational function with respect to + √ -λ (NOT λ ). In terms of (25) we consider the following nonlinear least squares problem

$$\epsilon = \min _ { P, Q } \int _ { - \infty + i \sigma } ^ { \infty + i \sigma } \left | \frac { P ( \vec { \vee } - \lambda ) } { Q ( \vec { \vee } - \lambda ) } - g _ { + } ( \lambda ) \right | ^ { 2 } | d \vec { \vee } - \lambda |,$$

where P , Q are polynomials with deg( P ) + 1 = deg( Q ) = d , and d is determined by making ϵ ≤ ϵ 0 , where ϵ 0 is a prescribed tolerance number. This nonlinear problem (27) is then solved using the technique of linearization and orthogonalization [2]. Finally, by expressing P/Q as a sum of poles, we arrive at

$$m _ { + } ( x _ { + }, \lambda ) \approx \tilde { m } _ { + } ( x _ { + }, \lambda ) = - \sqrt { - \lambda } + \sum _ { n = 1 } ^ { d } \frac { \alpha _ { n } } { \sqrt { - \lambda } + \beta _ { n } }.$$

Note that the coefficients α n and β n should appear as conjugate pairs due to the symmetry property (9). Unfortunately, it is not clear whether the rational approximation ˜ m + ( x + , λ ) in (28) still has the important Herglotz-property of the m-function. Moreover, the Herglotz-property cannot be checked by some conditions on the poles due to the leading square root in (28).

Figure 1: s -plane and k -plane.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00671v1-figure-8f637158-93cd-4dea-9333-1485b365bd14.jpg)

Applying the same idea to the m -function m -( x -, λ ) and we get the approximate boundary condition in the frequency domain

$$\hat { u } _ { x } ( x _ { \pm }, s ) = \left ( \mp \sqrt { - i s } + \sum _ { n = 1 } ^ { d _ { \pm } } \frac { \alpha _ { n, \pm } } { \sqrt { - i s } + \beta _ { n, \pm } } \right ) \hat { u } ( x _ { \pm }, s ).$$

If we introduce new unknowns ˆ w n, ± as

$$\hat { w } _ { n, \pm } = \frac { \hat { u } ( x _ { \pm }, s ) } { \stackrel { + } { \sqrt { - i } s } + \beta _ { n, \pm } },$$

then we can rewrite (29) as

$$\hat { u } _ { x } ( x _ { \pm }, s ) \pm \, \sqrt { - 1 } s \, \hat { u } = \sum _ { n = 1 } ^ { d _ { \pm } } \alpha _ { n, \pm } \hat { w } _ { n, \pm }, & & ( 3 1 a ) \\ \cdot \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \   & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, &\, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \, & \,$$

$$\stackrel { n = 1 } { \sqrt { - i s } } \, \hat { w } _ { n, \pm } + \beta _ { n, \pm } \hat { w } _ { n, \pm } = \stackrel { n = 1 } { \hat { u } } ( x _ { \pm }, s ), \quad n = 1, \dots, d _ { \pm }.$$

In the time domain, the approximate boundary condition is

$$u _ { x } ( x _ { \pm }, t ) \pm e ^ { - i \pi / 4 } \partial _ { t } ^ { \frac { 1 } { 2 } } u ( x _ { \pm }, t ) = \sum _ { n = 1 } ^ { d _ { \pm } } \alpha _ { n, \pm } w _ { n, \pm } ( t ), \$$

$$e ^ { - i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } w _ { n, \pm } ( t ) + \beta _ { n, \pm } w _ { n, \pm } ( t ) = u ( x _ { \pm }, t ), \quad n = 1, \dots, d _ { \pm }.$$

The final approximate truncated time-domain problem is formulated as

$$\mathrm i u _ { t } + \partial _ { x } ^ { 2 } u & = V ( x ) u, \quad ( x, t ) \in \Omega _ { \text{int} } \times ( 0, T ], \\ u ( x, 0 ) & = u _ { 0 } ( x ), \quad x \in \Omega _ { \text{int} },$$

$$u _ { x } ( x _ { \pm }, t ) \pm e ^ { - i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } u ( x _ { \pm }, t ) = \sum _ { n = 1 } ^ { d _ { \pm } } \alpha _ { n, \pm } w _ { n, \pm } ( t ),$$

$$e ^ { - i \pi / 4 } \, \partial _ { t } ^ { \frac { 1 } { 2 } } w _ { n, \pm } ( t ) + \beta _ { n, \pm } w _ { n, \pm } ( t ) = u ( x _ { \pm }, t ), \quad n = 1, \dots, d _ { \pm }.$$

## 6 Numerical results

In this section, we present some numerical results to test the accuracy of the proposed methods. In each example, the standard Crank-Nicolson scheme for time discretization is used. The fast evaluation of the half-order time derivative operator (23) is performed with the method of Zheng [41]. The computational domain is chosen to be Ω int = [ x -, x + ] = [ -5 5] . and the initial data is a Gaussian beam: u 0 ( x ) = e -x 2 +4i x . We use an 8 th -order FEM method with 1024 elements for the spatial discretization and a uniform time step of size ∆ = 10 t -4 .

## 6.1 The Free Schrödinger Equation

The exact solution of the free Schrödinger equation ( V ( x ) ≡ 0) is

$$u _ { \text{exa} } ( x, t ) = \sqrt { \frac { \text{i} } { - 4 t + \text{i} } } \exp \left ( \frac { - \text{i} x ^ { 2 } - 4 x + 1 6 t } { - 4 t + \text{i} } \right ).$$

We set the cut-off frequency f c in (20) to be f c = 256 and used a filtering function χ f ( ) = exp ( -(1 2 . f/f c ) 20 ) . The following Table 1 shows the relative L 2 -errors to the exact solution at certain time points when using 8097 quadrature points with Simpson's rule. Using this simple ex-

Table 1: Relative L 2 -error at certain time points for the free Schrödinger equation.

| Time points          |      0.5 |      0.6 |     0.7 |     0.8 |      0.9 |
|----------------------|----------|----------|---------|---------|----------|
| Relative L 2 -errors | 2.26e-07 | 3.46e-08 | 7.6e-09 | 5.6e-09 | 6.25e-09 |

ample, we can see that the frequency method with truncation and filtering works quite well: the magnitude of the relative L 2 -errors is at most on the order of 10 -7 .

## 6.2 The Coulomb-like Potential

In the second example we test the time-domain method with a Coulomb-like potential

$$V ( x ) = \frac { 1 } { \sqrt { 1 + x ^ { 2 } } }.$$

We set σ for Step 1 to σ = 1 and the tolerance number for the nonlinear least squares problem (27) to ϵ 0 = 10 -8 . Here we get 4 poles. Figure 2 shows the time evolution in a colored contour plot. One can see how the initial beam spreads out as time increases. Figure 3 shows how the relative

## Modulus

Figure 2: Time evolution of solution with Coulomb potential (34).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00671v1-figure-b2a3745d-ed02-44a1-bfaa-1b2824ff214b.jpg)

L 2 -error evolves in time. In this example with a varying external potential the magnitude of the relative L 2 -errors remains below 10 -5 .

## 6.3 The Gaussian Barrier

Next, we change the potential to a Gaussian barrier

$$V ( x ) = 3 0 e ^ { - 3 6 ( x - 8 ) ^ { 2 } }$$

with a height of 30, located in the exterior domain x &gt; 5 and centered at x = 8 . We set σ = 1 and ϵ 0 = 10 -4 and get 21 poles with the nonlinear least squares algorithm. Figure 4 shows the time evolution of the solution. One can clearly see how the initial beam propagates, spreads, and is (partially) reflected by the Gaussian barrier (35). The time evolution of the corresponding relative L 2 -error is shown in Figure 5. The relative L 2 errors remain below 5 × 10 -4 .

Unfortunately, the nonlinear least squares algorithm used in this paper failed to produce a rational approximation within an error tolerance much smaller than ϵ 0 . A more efficient algorithm is still desired, and this problem is currently under investigation.

Figure 3: Time evolution of the relative L 2 -error.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00671v1-figure-00386a56-a036-4244-9a4f-831bd39de9a2.jpg)

## Conclusion and Outlook

In this work we presented a new approach to simulate the solution to the Schrödinger equation with a general space-dependent potential in unbounded domains. Both frequency-domain and time-domain methods have been developed.

Future work will consist of implementing a more sophisticated algorithm for computing the m -function. Instead of solving the Riccati equation (12), we will consider computing the Weyl circles [14, 29] or the boundary control approach of Avdonin, Mikhaylov and Rybkin [10]. We will also look for a more stable algorithm for its rational approximation. It will also be clarified how this rational approximation can be made to preserve the essential Herglotz-property of the analytic m-function. This study will allow us to perform a rigorous stability analysis of this new approach. Finally, as a future goal, we want to extend our approach to the multi-dimensional Schrödinger problem, following the idea of Amrein and Pearson [3].

## Modulus

Figure 4: Time evolution of solution with Gaussian barrier (35).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00671v1-figure-f896eaee-2edf-409b-a939-9a3120937e1d.jpg)

## Appendix: Fast Evaluation Method

Here we present a short description of the method in [41] for evaluating the half-order time derivative ∂ 1 2 t . For any smooth function v = v t ( ) with v (0) = v ′ (0) = 0 , it is known that the semi-discrete half-order derivative

$$D _ { t } ^ { \frac { 1 } { 2 } } v ( t _ { n } ) \coloneqq \sqrt { \frac { 2 } { \Delta t } } \sum _ { m = 0 } ^ { n } \alpha _ { m } \, v ( t _ { n - m } )$$

with

$$\alpha _ { m } = \begin{cases} \beta _ { k } = \frac { ( 2 k )! } { 2 ^ { 2 k } ( k! ) ^ { 2 } } &, \ m = 2 k, \\ - \beta _ { k } &, \ m = 2 k + 1, \end{cases}$$

1

gives a second-order approximation of ∂ 2 t v t ( n ) (see [4, 41]). Suppose there exists a sum of decaying exponentials satisfying

$$\tilde { \beta } _ { k } = \sum _ { j = 1 } ^ { L } w _ { j } \, e ^ { - s _ { j } k }, \ \ s _ { j } > 0, \ | \beta _ { k } - \tilde { \beta } _ { k } | \leq \epsilon, \ k = 0, 1, \dots, \left [ \frac { N } { 2 } \right ].$$

Here N is the total number of time steps. If ϵ is small enough, it is reasonable to approximate (36) with

$$\tilde { D } _ { t } ^ { \frac { 1 } { 2 } } v ( t _ { n } ) \coloneqq \sqrt { \frac { 2 } { \Delta t } } ( v ( t _ { n } ) - v ( t _ { n - 1 } ) ) + \sqrt { \frac { 2 } { \Delta t } } \sum _ { m = 2 } ^ { n } \tilde { \alpha } _ { m } v ( t _ { n - m } ),$$

Figure 5: Time evolution of the relative L 2 -error.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00671v1-figure-ce1dc0d0-1690-4ff6-95a0-3e85b033b9fb.jpg)

where

$$\tilde { \alpha } _ { m } = \begin{cases} \tilde { \beta } _ { k } &, \ m = 2 k, \\ - \tilde { \beta } _ { k } &, \ m = 2 k + 1. \end{cases}$$

Set v k = v t ( k ) , v = ( v , v 0 1 , . . . ) , and define

$$\mathcal { F } _ { \text{odd} } ( w, s ; \mathbf v, k ) \coloneqq \sum _ { m = 1 } ^ { k } w e ^ { - s m } v _ { 2 k + 1 - 2 m }$$

and

$$\mathcal { F } _ { \text{even} } ( w, s ; \mathbf v, k ) \coloneqq \sum _ { m = 1 } ^ { k } w e ^ { - s m } v _ { 2 k - 2 m }.$$

Thus, F odd ( w,s ; v , 0) = F even ( w,s ; v , 0) = 0 . In addition, we have the following recursions

$$\mathcal { F } _ { \text{odd} } ( w, s ; \mathbf v, k ) \ & = \ e ^ { - s } \left [ w v _ { 2 k - 1 } + \mathcal { F } _ { \text{odd} } ( w, s ; \mathbf v, k - 1 ) \right ], \\ \mathcal { F } _ { \text{even} } ( w, s ; \mathbf v, k ) \ & = \ e ^ { - s } \left [ w v _ { 2 k - 2 } + \mathcal { F } _ { \text{even} } ( w, s ; \mathbf v, k - 1 ) \right ].$$

The summation (39) is then computed within O L ( ) operations as

$$\sum _ { m = 2 } ^ { n } \tilde { \alpha } _ { m } v _ { n - m } = \begin{cases} \sum _ { j = 1 } ^ { L } \mathcal { F } _ { \text{even} } ( w _ { j }, s _ { j } ; \nu, k ) - \sum _ { j = 1 } ^ { L } \mathcal { F } _ { \text{odd} } ( w _ { j }, s _ { j } ; \nu, k - 1 ) \quad, n = 2 k, \\ \sum _ { j = 1 } ^ { L } \mathcal { F } _ { \text{odd} } ( w _ { j }, s _ { j } ; \nu, k ) - \sum _ { j = 1 } ^ { L } \mathcal { F } _ { \text{even} } ( w _ { j }, s _ { j } ; \nu, k ) \quad, n = 2 k + 1. \end{cases}$$

In [41] for N = 1 000 000 , , , the authors found a sum of 81 decaying exponentials that approximates β k with an error of less than 5 0 . × 10 -11 .

## References

- [1] T. Aktosun, Mathematical Studies in Nonlinear Wave Propagation , Contemporary Mathematics 379 , American Mathematical Society, Providence, RI, 2005.
- [2] B. Alpert, L. Greengard and T. Hagstrom, Rapid evaluation of nonreflecting boundary kernels for time-domain wave propagation , SIAM J. Numer. Anal. 37 (2000), 1138-1164.
- [3] W.O. Amrein and D.B. Pearson, M operators: a generalisation of Weyl-Titchmarsh theory , J. Comput. Appl. Math. 171 (2004), 1-26.
- [4] X. Antoine and C. Besse, Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation, J. Comput. Phys. 181 (2003), 157-175.
- [5] X. Antoine, C. Besse and V. Mouysset, Numerical Schemes for the simulation of the twodimensional Schrödinger equation using non-reflecting boundary conditions , Math. Comp. 73 (2004), 1779-1799.
- [6] X. Antoine, C. Besse and S. Descombes, Artificial boundary conditions for one-dimensional cubic nonlinear Schrödinger equations , SIAM J. Numer. Anal. 43 (2006), 2272-2293.
- [7] X. Antoine, A. Arnold, C. Besse, M. Ehrhardt and A. Schädle, A Review of Transparent and Artificial Boundary Conditions Techniques for Linear and Nonlinear Schrödinger Equations , Commun. Comput. Phys. 4 (2008), 729-796.
- [8] X. Antoine, C. Besse and J. Szeftel, Towards accurate artificial boundary conditions for nonlinear PDEs through examples , CUBO 11 (2009), 29-48.
- [9] A. Arnold, M. Ehrhardt and I. Sofronov, Discrete transparent boundary conditions for the Schrödinger equation: fast calculation, approximation, and stability , Comm. Math. Sci. 1 (2003), 501-556.
- [10] S. Avdonin, V. Mikhaylov and A. Rybkin, The boundary control approach to the TitchmarshWeyl m-function , Comm. Math. Phys. 275 (2007), 791-803.
- [11] C. Bennewitz and W.N. Everitt, Some remarks on the Titchmarsh-Weyl m-coefficient. Proc. Pleijel Conference, Uppsala, Sweden, 1979, pp. 49-108.
- [12] C. Bennewitz, R. Weikard and M. Brown, Spectral and scattering theory for ordinary differential equations , Springer, 2020.
- [13] G. Borg, Uniqueness theorems in the spectral theory of y ′′ +( λ -q x ( )) y = 0 , Proc. 11th Scandinavian Congress of Mathematicians, Johan Grundt Tanums Forlag, Oslo, 1952, pp. 276-287.
- [14] B.M. Brown, V.G. Kirby and J.D. Pryce, A numerical method for the determination of the Titchmarsh-Weyl m-coefficient , Proc. R. Soc. Lond. A 435 (1991), 535-549.

| [15]   | Y. Chiba and S. Ohnishi, Resonance-state calculation applying the Weyl-Titchmarsh theory: Application for the quantum-confined Stark effects on excitons in a GaAs-Al x Ga 1 - x As quantum well , Phys. Rev. B 41 (1990), 6065-6068   |
|--------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [16]   | E.A. Coddington and N. Levinson, Theory of Ordinary Differential Equations , McGraw-Hill, New York, 1955.                                                                                                                              |
| [17]   | A.A. Danielyan and B.M. Levitan, On the asymptotic behavior of the Weyl-Titchmarsh m- function , Math. USSR Izvestiya 36 (1991), 487-496.                                                                                              |
| [18]   | M. Ehrhardt and A. Arnold, Discrete transparent boundary conditions for the Schrödinger equa- tion , Riv. Mat. Univ. Parma 6 (2001), 57-108.                                                                                           |
| [19]   | M. Ehrhardt and R.E. Mickens, Solutions to the discrete Airy equation: Application to parabolic equation calculations , J. Comput. Appl. Math. 172 (2004), 183-206.                                                                    |
| [20]   | M. Ehrhardt, Discrete transparent boundary conditions for Schrödinger-type equations for non- compactly supported initial data , Appl. Numer. Math. 58 (2008), 660-673.                                                                |
| [21]   | M. Ehrhardt and C. Zheng, Exact artificial boundary conditions for problems with periodic structures , J. Comput. Phys. 227 (2008), 6877-6894.                                                                                         |
| [22]   | M. Ehrhardt, H. Han and C. Zheng, Numerical simulation of waves in periodic structures , Commun. Comput. Phys. 5 (2009), 849-870.                                                                                                      |
| [23]   | F. Gesztesy and B. Simon, A new approach to inverse spectral theory, II. General real potentials and the connection to the spectral measure , Ann. Math. 152 (2000), 593-643.                                                          |
| [24]   | F. Gesztesy and M. Zinchenko, On spectral theory for Schrödinger operators with strongly sin- gular potentials , Math. Nachr. 279 (2006), 1041-1082.                                                                                   |
| [25]   | F. Gesztesy, R. Weikard and M. Zinchenko, Initial value problems and Weyl-Titchmarsh theory for Schrödinger operators with operator-valued potentials , Oper. Matrices 7(2) (2013), 241-283.                                           |
| [26]   | G.N. Gibson, G. Dunne and K.J. Bergquist, Tunneling Ionization Rates from Arbitrary Potential Wells , Phys. Rev. Lett. 81 (1998), 2663-2666.                                                                                           |
| [27]   | N. Gordon, J.P. Killingbeck and M.R.M. Witwit, Numerical Determination of the Titchmarsh- Weyl m-coefficient , Comput. Math. & Math. Phys. 40 , (2000), 108-117.                                                                       |
| [28]   | S. Jiang and L. Greengard, Fast evaluation of nonreflecting boundary conditions for the Schrödinger equation in one dimension , Comp. Math. Appl. 47 (2004), 955-966.                                                                  |
| [29]   | V.G. Kirby, A numerical method for determining the Titchmarsh-Weyl m-coefficient and its application to certain integro-differential inequalities , Ph.D. Thesis, University of Wales, College of Cardiff, UK, 1990.                   |
| [30]   | A. Kostenko, A. Sakhnovich and G. Teschl, Weyl-Titchmarsh theory for Schrödinger operators with strongly singular potentials , Int. Math. Res. Not. 8 (2012), 1699-1747.                                                               |
| [31]   | Y. Li and J.E Zhang, Option pricing with Weyl-Titchmarsh theory , Quant. Finance 4 (2004), 457-464.                                                                                                                                    |

- [32] V.A. Marchenko, Certain problems in the theory of second-order differential operators , Doklady Akad. Nauk SSSR 72 (1950), 457-460. (Russian).
- [33] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations , Springer, New York 1983.
- [34] A.D. Polyanin, Handbook of Linear Partial Differential Equations for Engineers and Scientists , Chapman &amp; Hall/CRC, 2002.
- [35] M. Reed and B. Simon, Methods of modern mathematical physics II: Fourier analysis, selfadjointness , Academic Press, San Diego 1975.
- [36] R. Sims, Reflectionless Sturm-Liouville equations , J. Comput. Appl. Math. 208 (2007), 207-225.
- [37] E.C. Titchmarsh, Eigenfunction Expansions Associated with Second-Order Differential Equations , Oxford University Press, London, 1946, Volumes 1 and 2.
- [38] H. Weyl, Über gewöhnliche Differentialgleichungen mit Singularitäten und die zugehörigen Entwicklungen willkürlicher Funktionen , Math. Ann. 68 (1910), 220-269.
- [39] M.R.M. Witwit, N. Gordon and J.P. Killingbeck, Numerical computation and analysis of the Titchmarsh-Weyl m(l) function for some simple potentials , J. Comput. Appl. Math. 106 (1999), 131-143.
- [40] J. Zhang, Z. Xu and X. Wu, Unified approach to split absorbing boundary conditions for nonlinear Schrödinger equations , Phys. Rev. E 78 (2008), 026709.
- [41] C. Zheng, Approximation, stability and fast evaluation of exact artificial boundary condition for the one-dimensional heat equation , J. Comput. Math. 25 (2007), 730-745.
- [42] C. Zheng, A perfectly matched layer approach to the nonlinear Schrödinger wave equations , J. Comput. Phys. 227 (2007), 537-556.
Matthias Ehrhardt Chunxiong Zheng
Absorbing Boundary Conditions for Variable Potential Schrödinger Equations via Titchmarsh-Weyl Theory
We propose a novel approach to simulate the solution of the time-dependent
Schr\"odinger equation with a general variable potential. The key idea is to
approximate the Titchmarsh-Weyl m-function (exact Dirichlet-to-Neumann
operator) by a rational function with respect to an appropriate spectral
parameter. By using this method, we overcome the usual high-frequency
restriction associated with absorbing boundary conditions in general variable
potential problems. The resulting fast computational algorithm for absorbing
boundary conditions ensures accuracy over the entire frequency range.