Rank-1 lattices and higher-order exponential splitting for the time-dependent Schr\"odinger equation

TL;DR

This paper introduces a high-dimensional numerical method combining Fourier pseudo-spectral discretization on rank-1 lattices with higher-order exponential splitting for solving the time-dependent Schrödinger equation, achieving dimension-independent convergence.

Contribution

The paper presents a novel approach that enables higher-order time convergence for high-dimensional Schrödinger equations using rank-1 lattice spatial discretization, independent of the problem dimension.

Findings

Higher-order convergence demonstrated in 2 to 8 dimensions.

Spatial discretization condition is independent of dimension.

Numerical results confirm practical effectiveness.

Abstract

In this paper, we propose a numerical method to approximate the solution of the time-dependent Schr\"odinger equation with periodic boundary condition in a high-dimensional setting. We discretize space by using the Fourier pseudo-spectral method on rank-$1$ lattice points, and then discretize time by using a higher-order exponential operator splitting method. In this scheme the convergence rate of the time discretization depends on properties of the spatial discretization. We prove that the proposed method, using rank-$1$ lattice points in space, allows to obtain higher-order time convergence, and, additionally, that the necessary condition on the space discretization can be independent of the problem dimension $d$. We illustrate our method by numerical results from 2 to 8 dimensions which show that such higher-order convergence can really be obtained in practice.

Figures

Tables

Table 1. Table 1: Parameters of the rank- 1 1 1 lattice points for our numerical results.

$d$	$n$	${𝒛}^{\top }$
$2$	$2^{16}$	$(1,100135)$
$4$	$2^{20}$	$(1,443165,95693,34519)$
$6$	$2^{24}$	$(1,6422017,7370323,2765761,8055041,2959639)$
$8$	$2^{24}$	$(1,6422017,7370323,2765761,8055041,2959639,7161203,4074015)$

Table 2. Table 2: Coefficients for the sixth-order method, calculated based on [ 10 ] .

	$a_j$		$b_j$
$j=1,9$	0.392161444007314	$j=1,10$	0.196080722003657
$j=2,8$	0.332599136789359	$j=2,9$	0.362380290398337
$j=3,7$	-0.706246172557639	$j=3,8$	-0.186823517884140
$j=4,6$	0.0822135962935508	$j=4,7$	-0.312016288132044
$j=5$	0.798543990934830	$j=5,6$	0.440378793614190
$j=10$	0

Table 3. Table 3: Coefficients for the eighth-order method, calculated based on [ 10 ] .

	$a_j$		$b_j$
$j=1,17$	0.130202483088890	$j=1,18$	0.0651012415444450
$j=2,16$	0.561162981775108	$j=2,17$	0.345682732431999
$j=3,15$	-0.389474962644847	$j=3,16$	0.0858440095651306
$j=4,14$	0.158841906555156	$j=4,15$	-0.115316528044846
$j=5,13$	-0.395903894133238	$j=5,14$	-0.118530993789041
$j=6,12$	0.184539640978316	$j=6,13$	-0.105682126577461
$j=7,11$	0.258374387686322	$j=7,12$	0.221457014332319
$j=8,10$	0.295011723609310	$j=8,11$	0.276693055647816
$j=9$	-0.605508533830035	$j=9,10$	-0.155248405110362
$j=18$	0

Equations

$$\mathrm{i}\gamma \frac{\partial u}{\partial t}(\mathbf{x},t)$$

$$u(\mathbf{x},0)$$

$$\mathrm{i}\hslash \frac{\partial \psi }{\partial t}(\mathbf{x},t)$$

$$\Lambda (\mathbf{z},n):=\left\{\frac{\mathbf{z}k}{n}\mathrm{mod}1:k\in \mathbb{Z}\right\}.$$

$$\mathbf{z}\cdot \mathbf{h}\not\equiv \mathbf{z}\cdot {\mathbf{h}}^{\prime }(\mathrm{mod}n)\text{for all }\mathbf{h},{\mathbf{h}}^{\prime }\in \mathcal{A}(\mathbf{z},n),\mathbf{h}\ne {\mathbf{h}}^{\prime }.$$

$$\displaystyle\begin{split}\mathbb{Z}^{d}&=\biguplus_{\makebox[36.15898pt]{$\scriptstyle{\boldsymbol{h}}\in\Lambda^{\perp}(\boldsymbol{z},n)$}}\left({\boldsymbol{h}}+\mathcal{A}(\boldsymbol{z},n)\right)\\ &=\biguplus_{\makebox[36.15898pt]{$\scriptstyle{\boldsymbol{h}}\in\mathcal{A}(\boldsymbol{z},n)$}}\{{\boldsymbol{h}}^{\prime}\in\mathbb{Z}^{d}:\boldsymbol{z}\cdot{\boldsymbol{h}}^{\prime}\equiv\boldsymbol{z}\cdot{\boldsymbol{h}}\pmod{n}\}\\ &=\biguplus_{\makebox[36.15898pt]{$\scriptstyle j\in\mathbb{Z}_{n}$}}\{{\boldsymbol{h}}\in\mathbb{Z}^{d}:\boldsymbol{z}\cdot{\boldsymbol{h}}\equiv j\pmod{n}\},\end{split}$$

$$E_{\alpha }({\mathbb{T}}^d):=\left\{f\in L_2({\mathbb{T}}^d):\Vert f{\Vert }_{E_{\alpha }({\mathbb{T}}^d)}^2:=\sum _{\mathbf{h}\in {\mathbb{Z}}^d}|\hat{f}(\mathbf{h}){|}^2{r}_{\alpha }^2(\mathbf{h})<\infty \right\},$$

$$r_{\alpha }^2(\mathbf{h}):=\prod _{j=1}^d\max (|h_j{|}^{2\alpha },1).$$

$$A({\mathbb{T}}^d):=\{f\in L_2({\mathbb{T}}^d):\Vert f{\Vert }_{A({\mathbb{T}}^d)}:=\sum _{\mathbf{h}\in {\mathbb{Z}}^d}|\hat{f}(\mathbf{h})|<\infty \}.$$

$$u(\mathbf{x},t)=\sum _{\mathbf{h}\in {\mathbb{Z}}^d}\hat{u}(\mathbf{h},t)\exp (2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{x}),$$

$$\mathrm{i}\gamma {\hat{u}}^{\prime }(\mathbf{h},t)=2{\pi }^2{\gamma }^2\Vert \mathbf{h}{\Vert }_2^2\hat{u}(\mathbf{h},t)+\hat{f}(\mathbf{h},t),$$

$$u_a(\mathbf{x},t):=\sum _{\mathbf{h}\in \mathcal{A}(\mathbf{z},n)}{\hat{u}}_a(\mathbf{h},t)\exp (2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{x}),$$

$${\hat{u}}_a(\mathbf{h},t):=\frac{1}{n}\sum _{\mathbf{p}\in \Lambda (\mathbf{z},n)}u(\mathbf{p},t)\exp (-2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{p}).$$

$$\frac{1}{n}\sum _{\mathbf{p}\in \Lambda (\mathbf{z},n)}\exp (2\pi \mathrm{i}(\mathbf{h}-{\mathbf{h}}^{\prime })\cdot \mathbf{p})={\delta }_{\mathbf{h},{\mathbf{h}}^{\prime }},$$

$$\frac{1}{n}\sum _{\mathbf{h}\in \mathcal{A}(\mathbf{z},n)}\exp (2\pi \mathrm{i}\mathbf{h}\cdot (\mathbf{p}-{\mathbf{p}}^{\prime }))={\delta }_{\mathbf{p},{\mathbf{p}}^{\prime }}.$$

$$u_a(\mathbf{p})=u(\mathbf{p}).$$

$$\mathbf{u}$$

$${\mathbf{u}}_a$$

$${\hat{\mathbf{u}}}_a$$

$${\hat{u}}_a(\mathbf{h})=\sum _{{\mathbf{h}}^{\prime }\in {\Lambda }^{\perp }(\mathbf{z},n)}\hat{u}(\mathbf{h}+{\mathbf{h}}^{\prime })=\hat{u}(\mathbf{h})+\sum _{0\ne {\mathbf{h}}^{\prime }\in {\Lambda }^{\perp }(\mathbf{z},n)}\hat{u}(\mathbf{h}+{\mathbf{h}}^{\prime }).$$

$$\sum _{\mathbf{p}\in \Lambda (\mathbf{z},n)}f(\mathbf{p})\overline{g(\mathbf{p})}=\sum _{\mathbf{h}\in \mathcal{A}(\mathbf{z},n)}\widehat{f_a}(\mathbf{h})\overline{\widehat{g_a}(\mathbf{h})}$$

$$\hat{f}(h)=\frac{1}{n}\sum _{k=0}^{n-1}f(k/n)\exp (-2\pi \mathrm{i}hk/n),$$

$$f(k/n)=\sum _{h=0}^{n-1}\hat{f}(h)\exp (2\pi \mathrm{i}hk/n).$$

$$\hat{f}({\mathbf{h}}_{\xi })=\frac{1}{n}\sum _{k=0}^{n-1}f({\mathbf{p}}_k)\exp (-2\pi \mathrm{i}{\mathbf{h}}_{\xi }\cdot {\mathbf{p}}_k)=\frac{1}{n}\sum _{k=0}^{n-1}f({\mathbf{p}}_k)\exp (-2\pi \mathrm{i}\xi k/n),$$

$$f({\mathbf{p}}_k)=\sum _{\xi =0}^{n-1}\hat{f}({\mathbf{h}}_{\xi })\exp (2\pi \mathrm{i}{\mathbf{h}}_{\xi }\cdot {\mathbf{p}}_k)=\sum _{k=0}^{n-1}\hat{f}({\mathbf{h}}_{\xi })\exp (2\pi \mathrm{i}\xi k/n),$$

$$y^{\prime }(t)=(A+B)y(t),y(0)=y_0,$$

$$y(t+\Delta t)\approx {\mathrm{e}}^{b_{1}B\Delta t}{\mathrm{e}}^{a_{1}A\Delta t}\cdots {\mathrm{e}}^{b_{s}B\Delta t}{\mathrm{e}}^{a_{s}A\Delta t}y(t),$$

$$\Vert {\mathrm{e}}^{b_{1}B\Delta t}{\mathrm{e}}^{a_{1}A\Delta t}\cdots {\mathrm{e}}^{b_{s}B\Delta t}{\mathrm{e}}^{a_{s}A\Delta t}y(t)-{\mathrm{e}}^{(A+B)\Delta t}y(t){\Vert }_X\le C(\Delta t{)}^{p+1},$$

$$y_{k+1}={\mathrm{e}}^{b_{1}B\Delta t}{\mathrm{e}}^{a_{1}A\Delta t}\cdots {\mathrm{e}}^{b_{s}B\Delta t}{\mathrm{e}}^{a_{s}A\Delta t}{y}_k,y_{\{k=0\}}=y_0.$$

$$\Vert y_m-y(t){\Vert }_X\le Cm\Delta t(\Delta t{)}^p=Ct(\Delta t{)}^p.$$

$${\mathrm{ad}}_A^p(B)=[A,{\mathrm{ad}}_A^{p-1}(B)],{\mathrm{ad}}_A^0(B)=B,$$

$$\mathrm{i}\gamma {\hat{\mathbf{u}}}_t^{\prime }$$

Full text

1 Introduction

Rank-$1$ lattice points have been widely used in the context of high-dimensional problems. Their traditional usage is in numerical integration, see, e.g., [5, 19] and references therein. In this work, we use rank-$1$ lattice points for function approximation, to approximate the solution of the time-dependent Schrödinger equation (TDSE). Function approximation using rank-$1$ lattice points has recently received more attention, see, e.g., [3, 11, 12, 13, 21, 14]. In [13], Li and Hickernell introduced the pseudo-spectral Fourier collocation method using rank-$1$ lattice rules. Due to the rank-$1$ lattice structure, Fourier pseudo-spectral methods can be efficiently implemented using one-dimensional Fast Fourier transformations (FFTs). This is well known, and we state the exact form in Theorem 2 together with other useful properties of approximations on rank-$1$ lattice points.

To simulate many particles in the quantum world is a computationally challenging problem. For the TDSE, the dimensionality of the problem increases with the number of particles of the system. In the present paper, the following form is considered:

[TABLE]

where $\mathrm{i}$ represents the imaginary unit, $\boldsymbol{x}$ is the spatial position in the $d$-dimensional torus $\mathbb{T}^{d}=\mathbb{T}([0,1]^{d})\simeq[0,1]^{d}$, the time $t$ is positive valued, and $\gamma$ is a small positive parameter. The function $u(\boldsymbol{x},t)$ is the sought solution, while $v(\boldsymbol{x})$ and $g(\boldsymbol{x})$ are the potential and initial conditions respectively. The Laplacian can be interpreted as $\Delta=\sum_{i=1}^{M}\sum_{j=1}^{D}\partial^{2}/\partial x_{i,j}^{2}=\sum_{i=1}^{d}\partial^{2}/\partial x_{i}^{2}$ where $M$ is the number of particles and $D$ is the physical dimensionality with $MD=d$. We note that the above form of the TDSE becomes equivalent after substitution to the following form which is common in the context of physics:

[TABLE]

where $\hbar$ is the reduced Planck constant and $m$ is the mass.

The form (1) of the TDSE has been studied from various perspectives of numerical analysis [9, 7, 22, 15]. In the present paper, we focus on two perspectives; high-dimensionality and higher-order convergence in time stepping. For the first point, Gradinaru [7] proposed to use sparse grids for the physical space. In [21], the current authors used rank-$1$ lattice points to prove second order convergence for the time discretization using Strang splitting and numerically compared results with the sparse grid approach from [7]. The numerical result using rank-$1$ lattices showed the expected second order convergence even up to $12$ dimensions. Hence rank-$1$ lattice points perform thereby much better than the sparse grid approach. The second point, higher-order convergence in time stepping, is successfully achieved by Thalhammer [22] using higher-order exponential operator splitting. In that paper, the spatial discretization was done by a full grid and therefore was limited to lower dimensional cases ($d\leq 3$).

The rest of this paper is organized as follows: Section 2 describes the proposed method consisting of the higher-order exponential splitting method and Fourier pseudo-spectral method using rank-$1$ lattices. Section 3 shows numerical results with various settings. The main aim here is to show higher-order time stepping convergence in higher-dimensional cases. Finally, Section 4 concludes the present paper with a short summary.

Throughout the present paper, we denote the set of integer numbers by $\mathbb{Z}$ and the ring of integers modulo $n$ by $\mathbb{Z}_{n}:=\{0,1,...,n-1\}$. We distinguish between the normal equivalence in congruence modulo $n$ as $a\equiv b\pmod{n}$ and the binary operation modulo $n$ denoted by $\bmod{n}$ which returns the corresponding value in $\mathbb{Z}_{n}$ for $\bmod{n}$ and in $\mathbb{T}$ for $\bmod{1}$.

2 The numerical method

In this section, we describe necessary ingredients of our method. For the conciseness, we restrict ourselves to the rank-$1$ lattice points instead of general rank-$r$ lattice points. However, our method is indeed possible to generalize to rank-$r$ lattice points, similar as in [21].

We use a rank-$1$ lattice point set and an associated anti-aliasing set for the Fourier pseudo-spectral method. For using the Fourier pseudo-spectral method, one obvious choice is regular grids [22], but the number of points increases too quickly in terms of the number of dimensions. To mitigate this problem, Gradinaru [7, 6] proposed to use sparse grids. For the same reason we introduced lattice points in [21] to get first and second order time convergence, and obtained much better results compared to [7, 6].

2.1 Rank-$1$ lattice point sets and the associated anti-aliasing sets

A rank-$1$ lattice point set $\Lambda(\boldsymbol{z},n)$ is fully determined by the modulus $n$ and a generating vector $\boldsymbol{z}\in\mathbb{Z}_{n}^{d}$:

[TABLE]

Usually, all components of the generating vector are chosen to be relatively prime to $n$ which means all points have different values in each coordinate and the number of points is exactly $n$. The generating vector determines the quality of the rank-$1$ lattice points. Of course, the quality criterion needs to take into account what the lattice points will be used for. A well studied setting is numerical integration, e.g., [19, 17] and [16, Chapter 5]. Function approximation using lattice points is relatively new. In that context, we refer to [11, 12, 3]. We call $\mathcal{A}(\boldsymbol{z},n)\subset\mathbb{Z}^{d}$ an anti-aliasing set for the lattice point set $\Lambda(\boldsymbol{z},n)$ if

[TABLE]

We remark that the anti-aliasing set is not uniquely determined and the cardinality $|\mathcal{A}(\boldsymbol{z},n)|\leq n$. By using the dual lattice $\Lambda^{\bot}(\boldsymbol{z},n):=\{{\boldsymbol{h}}\in\mathbb{Z}^{d}:\boldsymbol{z}\cdot{\boldsymbol{h}}\equiv 0\pmod{n}\}$, we can rewrite the condition as ${\boldsymbol{h}}-{\boldsymbol{h}}^{\prime}\not\in\Lambda^{\bot}(\boldsymbol{z},n)$ for ${\boldsymbol{h}}\neq{\boldsymbol{h}}^{\prime}$. If we have the full cardinality $|\mathcal{A}(\boldsymbol{z},n)|=n$, we can divide $\mathbb{Z}^{d}$ into conjugacy classes:

[TABLE]

where $\uplus$ is the union of conjugacy classes.

2.2 Korobov spaces

Rank-$1$ lattices are closely related to Korobov spaces which are reproducing kernel Hilbert spaces of Fourier series. The Korobov space $E_{\alpha}(\mathbb{T}^{d})$ is given by

[TABLE]

where

[TABLE]

The parameter $\alpha\geq 1/2$ is called the smoothness parameter which determines the rate of decay of the Fourier coefficients. To ensure regularity of the solution of the TDSE (1) and to prove that our method gives higher-order convergence for the temporal discretization, we will assume that the initial condition $g(\boldsymbol{x})$ and the potential function $v(\boldsymbol{x})$ are in the Korobov space with given smoothness, see Lemma 1 and Theorem 4.

2.3 Fourier pseudo-spectral methods using rank-$1$ lattices

We approximate the solution of the TDSE (1) by the truncated Fourier series. To ensure the solution to be regular enough so that the Fourier expansion makes sense (e.g., uniqueness, continuity, point-wise convergence), we require all functions to be in Wiener algebra $A(\mathbb{T}^{d})$:

[TABLE]

For $\alpha>1/2$, we have $E_{\alpha}(\mathbb{T}^{d})\subset A(\mathbb{T}^{d})$. The following lemma shows the regularity of the solution, and the TDSE (1) in terms of Fourier coefficients and was already stated and proven in [21].

Lemma 1** (Regularity of solution and Fourier expansion).**

Given the TDSE (1) with $v,g\in E_{\alpha}(\mathbb{T}^{d})$ and $\alpha\geq 2$, then the solution $u(\boldsymbol{x},t)\in E_{\alpha}(\mathbb{T}^{d})$ for all finite $t\geq 0$ and therefore

[TABLE]

with

[TABLE]

for all ${\boldsymbol{h}}\in\mathbb{Z}^{d}$, with $\widehat{u}^{\prime}({\boldsymbol{h}},t)=(\partial/\partial t)\,\widehat{u}({\boldsymbol{h}},t)$ and $\widehat{f}({\boldsymbol{h}},t)$ the Fourier coefficients of $f(\boldsymbol{x},t):=u(\boldsymbol{x},t)\,v(\boldsymbol{x})$.

We then truncate the Fourier series (4) to a finite sum on an anti-aliasing set $\mathcal{A}(\boldsymbol{z},n)$ associated to a rank-$1$ lattice $\Lambda(\boldsymbol{z},n)$ to get the approximation

[TABLE]

with the approximated coefficients calculated by the rank-$1$ lattice rule

[TABLE]

The subscript $a$ of $u_{a}(\boldsymbol{x},t)$ and $\widehat{u}_{a}({\boldsymbol{h}},t)$ indicates that these are approximations of $u(\boldsymbol{x},t)$ and $\widehat{u}({\boldsymbol{h}},t)$ respectively. For simplicity of notation, we omit the time $t$ in the rest of this section. Due to the rank-$1$ lattice structure and by choosing the anti-aliasing set to be of full size, we have the following properties:

Theorem 2**.**

*Given a rank-$1$ lattice point set $\Lambda(\boldsymbol{z},n)$ and a corresponding anti-aliasing set $\mathcal{A}(\boldsymbol{z},n)$ with $|\mathcal{A}(\boldsymbol{z},n)|=n$, the following properties hold.

(i) (Character property and dual character property) For any two vectors ${\boldsymbol{h}},{\boldsymbol{h}}^{\prime}\in\mathcal{A}(\boldsymbol{z},n)$*

[TABLE]

where $\delta_{\boldsymbol{p},\boldsymbol{p}^{\prime}}$ is the Kronecker delta function that is $1$ if $\boldsymbol{p}=\boldsymbol{p}^{\prime}$ and [math] otherwise. Also, for any two lattice points $\boldsymbol{p},\boldsymbol{p}^{\prime}\in\Lambda(\boldsymbol{z},n)$

[TABLE]

(ii) (Interpolation condition) If $u_{a}$ is the approximation of a function $u\in A(\mathbb{T}^{d})$ by truncating its Fourier series expansion to the anti-aliasing set $\mathcal{A}(\boldsymbol{z},n)$ and by calculating the coefficients by the rank-$r$ lattice rule, cfr. (6) and (7), then for any $\boldsymbol{p}\in\Lambda(\boldsymbol{z},n)$

[TABLE]

(iii) (Mapping through FFT) Define the following vectors:

[TABLE]

*with $\boldsymbol{p}_{k}=\boldsymbol{z}k/n\bmod{1}\in\Lambda(\boldsymbol{z},n)$, and where ${\boldsymbol{h}}_{\xi}\in\mathcal{A}(\boldsymbol{z},n)$ is chosen such that ${\boldsymbol{h}}\cdot\boldsymbol{z}_{\xi}\equiv\xi\pmod{n}$. Then $\boldsymbol{u}=\boldsymbol{u}_{a}$ (by (ii)) is the collection of function values $u(\boldsymbol{p})$ on the lattice points $\boldsymbol{p}\in\Lambda(\boldsymbol{z},n)$ and $\widehat{\boldsymbol{u}}_{a}$ is the collection of Fourier coefficients $\widehat{u}_{a}({\boldsymbol{h}})$ (by using the lattice rule, cfr. (6) and (7)) on the anti-aliasing indices ${\boldsymbol{h}}\in\mathcal{A}(\boldsymbol{z},n)$. The $1$-dimensional discrete Fourier transform and its inverse now maps $\boldsymbol{u}_{a}\in\mathbb{C}^{n}$ to $\widehat{\boldsymbol{u}}_{a}\in\mathbb{C}^{n}$ and back.

(iv) (Aliasing) The approximated Fourier coefficients (7) through the lattice rule $\Lambda(\boldsymbol{z},n)$ alias the true Fourier coefficients in the following way*

[TABLE]

Proof.

We refer to [21, Theorem 2 and Lemma 3] where more general statement for rank-$r$ lattices can be found. ∎

We remark that the above theorem can also be understood in terms of Fourier analysis on a finite Abelian group where the group, normally denoted as $G$, is the rank-$1$ lattice point set $\Lambda(\boldsymbol{z},n)$ and the associated character group (Pontryagin dual) $\widehat{G}:=\{\exp(2\pi\mathrm{i}\,{\boldsymbol{h}}\cdot\circ):{\boldsymbol{h}}\in\mathcal{A}(\boldsymbol{z},n)\}$ with $|\mathcal{A}(\boldsymbol{z},n)|=n$. The (dual) character property is then to be understood as orthonormality of $\widehat{G}$ on $L_{2}(G)$. The interpolation condition can be seen as the representability of functions by using Fourier series. Due to this structure, the Plancherel theorem also holds:

[TABLE]

for $f,g\in L_{2}(G)$.

For readers who are not familiar with Fourier transforms on a rank-$1$ lattice, one intuitive way of seeing why one-dimensional FFTs are available is the following. The usual one-dimensional Fourier transform for equidistant points which is a scalar multiple of a unitary Fourier transform, for a function $f:\mathbb{T}\to\mathbb{C}$, can be written as

[TABLE]

and the inverse

[TABLE]

Now we see that the Fourier transform on a rank-$1$ lattice has the exact same structure for a function $f:\mathbb{T}^{d}\to\mathbb{C}$,

[TABLE]

and

[TABLE]

where we note $\boldsymbol{p}_{k}=\boldsymbol{z}k/n\bmod{1}$ and ${\boldsymbol{h}}_{\xi}\cdot\boldsymbol{z}\equiv\xi\pmod{n}$. Hence we only need one-dimensional FFTs to transform functions on $\mathbb{T}^{d}$.

2.4 Higher-order exponential splitting

For the temporal discretization, we employ a higher-order exponential splitting scheme (also called an exponential propagator), see, e.g., [1, 22, 20]. To describe the higher-order exponential splitting, let us consider the following ordinary differential equation:

[TABLE]

where $A$ and $B$ are differential operators. The solution for the equation (11) is $y(t)=\mathrm{e}^{(A+B)t}y_{0}$. However, often it is not possible to compute this exactly, and one needs to approximate the quantity with cheap computational cost. When both $\mathrm{e}^{At}$ and $\mathrm{e}^{Bt}$ can be computed easily, the higher-order exponential splitting is a powerful tool to approximate the solution $\mathrm{e}^{(A+B)t}y_{0}$. The approximated solution for this case is given by:

[TABLE]

where $a_{i}$ and $b_{i}$, $i=1,\ldots,s$, are coefficients determined by the desired order of convergence $p$. In other words, if the splitting (12) satisfies

[TABLE]

for some normed space $X$, where the constant $C$ is independent of ${\Delta t}$, then the splitting is said to have $p$-th order. The number of steps $s$ and the coefficients $a_{i}$, $b_{i}$ can be determined according to the order $p$, see [8] for details. We evolve the time using this discretization from time [math], i.e.,

[TABLE]

By summing up the local errors (13) of each step $k=1,\ldots,m$, where $t=m\Delta t$, gives the total error:

[TABLE]

We call this quantity the total error in the $L_{2}$ sense, and this is the reason why the splitting is called to be of $p$-th order. The error coming from the exponential splitting can be related to commutators of two operators $A$ and $B$, namely $[A,B]:=AB-BA$, $[A,[A,B]]:=A^{2}B-2ABA+BA^{2}$, etc. We introduce the notation for the $p$-th commutator by following [22]:

[TABLE]

where $p\geq 1$. When the $p$-th commutator is bounded, it is known that the $p$-th order exponential splitting gives the desired order, see [22, Lemma 1 and Theorem 1]. We also refer to [9, Theorem 2.1] for the second-order splitting (namely, Strang splitting) in a more abstract setting.

2.5 Higher-order exponential splitting on rank-$1$ lattices

We apply the higher-order exponential splitting to the space discretized TDSE in this section. For solving the TDSE (1) in the dual space with finite number of Fourier basis functions, we will rewrite the problem in vector form. We let \widehat{\boldsymbol{u}}_{t}:=\bigl{(}\widehat{u}_{a}({\boldsymbol{h}}_{0},t), \dots,\widehat{u}_{a}({\boldsymbol{h}}_{n-1},t)\bigr{)} the approximated solution at time $t$. Throughout time evolution, we use a fixed anti-aliasing set $\mathcal{A}(\boldsymbol{z},n)=\{{\boldsymbol{h}}_{\xi}\mathrel{:}\xi=0,\ldots,n-1\}$ of full size $|\mathcal{A}(\boldsymbol{z},n)|=n$, where we denote ${\boldsymbol{h}}_{\xi}\in\mathcal{A}(\boldsymbol{z},n)$ as such a vector that ${\boldsymbol{h}}_{\xi}\cdot\boldsymbol{z}\equiv\xi\pmod{n}$. We obtain the following relation by imposing that (5) holds for all ${\boldsymbol{h}}\in\mathcal{A}(\boldsymbol{z},n)$,

[TABLE]

with the initial condition $\widehat{\boldsymbol{u}}_{0}=\widehat{\boldsymbol{g}}_{a}:=(\widehat{g}_{a}({\boldsymbol{h}}_{0}),\dots,\widehat{g}_{a}({\boldsymbol{h}}_{n-1}))$,

[TABLE]

and the potential multiplication operator $W_{n}:=F_{n}V_{n}F_{n}^{-1}$ with

[TABLE]

where $F_{n}$ is the unitary Fourier matrix

[TABLE]

The approximation of the multiplication operator, $W_{n}$, is justified by the following lemma which is taken from [21].

Lemma 3** (Multiplication operator on rank-$1$ lattices).**

Given a rank-$1$ lattice point set $\Lambda(\boldsymbol{z},n)$ and corresponding anti-aliasing set $\mathcal{A}(\boldsymbol{z},n)$ of full size, a potential function $v\in E_{\alpha}(\mathbb{T}^{d})$ with $\alpha\geq 2$ and a function $u_{a}\in E_{\beta}(\mathbb{T}^{d})$ with $\beta\geq 2$ with Fourier coefficients only supported on $\mathcal{A}(\boldsymbol{z},n)$. Then the action in the Fourier domain restricted to $\mathcal{A}(\boldsymbol{z},n)$ of multiplying with $v$, that is $f_{a}(\boldsymbol{x})=v(\boldsymbol{x})\,u_{a}(\boldsymbol{x})$, on the nodes of the rank-$1$ lattice, and with $f_{a}$ having Fourier coefficients restricted to the set $\mathcal{A}(\boldsymbol{z},n)$, can be described by a circulant matrix $W_{n}\in\mathbb{C}^{n\times n}$ with $W_{n}=F_{n}\,V_{n}\,F_{n}^{-1}$, with $V_{n}$ given by (16) and $F_{n}$ the unitary Fourier matrix, where the element at position $(\xi,\xi^{\prime})$ of $W_{n}$ is given by

[TABLE]

Proof.

We refer to [21, Lemma 5]. ∎

We approximate the solution of the ordinary differential equation (14)

[TABLE]

by applying the higher-order exponential splitting method (12):

[TABLE]

where again the coefficients $a_{i},b_{i}$ are determined according to the desired order of convergence, and

[TABLE]

The approximated solution at the time $t=k{\Delta t}$ is then obtained by stepping time ${\Delta t}$ iteratively by (18). In the following we show the commutator bounds which correspond to [22, Hypothesis 3] and lead us to the total bound as in [22, Theorem 1].

Theorem 4** ($p$-th commutator bound and total error bound).**

Given a rank-$1$ lattice with generating vector $\boldsymbol{z}\in\mathbb{Z}^{d}$ and modulus $n$ and a TDSE with a potential function $v\in E_{\alpha}(\mathbb{T}^{d})$ with $\alpha>2p+1/2$ and an initial condition $g\in E_{\beta}(\mathbb{T}^{d})$ with $\beta\geq 2$. Let $D=\tfrac{\gamma}{2}D_{n}$ and $W=\frac{1}{\gamma}W_{n}$ with $D_{n}$ and $W_{n}=F_{n}V_{n}F_{n}^{-1}$ as defined in (15) and (17), and with $V_{n}$ as defined in (16) using the potential function $v$.

If the anti-aliasing set $\mathcal{A}(\boldsymbol{z},n)=\{{\boldsymbol{h}}_{\xi}\in\mathbb{Z}^{d}:{\boldsymbol{h}}_{\xi}\cdot\boldsymbol{z}\equiv\xi\pmod{n}\text{ for }\xi=0,\ldots,n-1\}$, with full cardinality, is chosen such that it has minimal $\ell_{2}$ norm, i.e.,

[TABLE]

with

[TABLE]

then for all $\boldsymbol{y}\in\mathbb{R}^{n}$ we have the following bound for the $p$-th commutator:

[TABLE]

*where $c$ is a constant independent of $n$ and $\boldsymbol{y}$.

This commutator condition and [22, Theorem 1] directly give us the total error bound for (14):*

[TABLE]

where $m{\Delta t}=t$ and the constants depend on $t$ but not on $m$ or ${\Delta t}$.

Proof.

Let $M:=\mathrm{ad}_{D}^{p}(W)\;(D+I)^{-p}$. Since $(D+I)^{p}$ is is non-singular, the claim of the theorem is equivalent to the assertion that the induced $\ell_{2}$ norm of the matrix $\|M\|_{2}:=\sup_{\boldsymbol{0}\neq\boldsymbol{y}\in\mathbb{R}^{n}}\|M\,\boldsymbol{y}\|_{2}/\|\boldsymbol{y}\|_{2}$ is bounded independent of $n$. Each element of the matrix $M$ is given by,

[TABLE]

where the constant $c_{1}=1/(2\pi\gamma)^{p}>0$. Now we bound $\|M\|_{2}$ by using $\|M\|_{2}\leq\sqrt{\|M\|_{1}\|M\|_{\infty}}$. First we bound $\|M\|_{1}$:

[TABLE]

We notice that the diagonal components of $M$ ($\xi=\xi^{\prime}$) is always [math], hence we exclude such cases in the following argument. Because we collect the anti-aliasing set by minimizing the $\ell_{2}$ norm (19), we have $\|{\boldsymbol{h}}_{\xi}\|_{2}\leq\|{\boldsymbol{h}}_{\xi}^{\prime}\|_{2}$ for any ${\boldsymbol{h}}_{\xi}^{\prime}\in A(\boldsymbol{z},n,\xi)$. In particular, this holds for ${\boldsymbol{h}}_{\xi}^{\prime}={\boldsymbol{h}}_{\xi-\xi^{\prime}}+{\boldsymbol{h}}_{\xi^{\prime}}$ since $({\boldsymbol{h}}_{\xi-\xi^{\prime}}+{\boldsymbol{h}}_{\xi^{\prime}})\cdot\boldsymbol{z}\equiv\xi\pmod{n}$ for any choice of $\xi^{\prime}=0,\ldots,n-1$. This gives us the connection between $\|{\boldsymbol{h}}_{\xi}\|_{2}$ and $\|{\boldsymbol{h}}_{\xi^{\prime}}\|_{2}$ using $\|{\boldsymbol{h}}_{\xi-\xi^{\prime}}\|_{2}$:

[TABLE]

for $\xi\neq\xi^{\prime}$. We continue from the above bound of $\|M\|_{1}$,

[TABLE]

For the last inequality, we used the conjugacy decomposition (2). By using Cauchy–Schwarz inequality and multiplying and dividing by $r_{\alpha}$, we have

[TABLE]

This means we have bounded $\|M\|_{1}$ independent of $n$. For $\|M\|_{\infty}$ we can proceed in a similar way to obtain

[TABLE]

Therefore, we have $\|M\|_{2}<\infty$ independent of $n$. The total error bound directly follows from this commutator bound and [22, Theorem 1]. ∎

3 Numerical results

We demonstrate our method by showing some numerical results in this section. We construct rank-$1$ lattices by using the component-by-component (CBC) construction [4, 17]. The code for producing the rank-$1$ lattice is available online [18], fastrank1expt.m. With the script, we choose $n$ being a power of $2$ and generate the vector $\boldsymbol{z}$ which is optimized for integration in (unweighted) Korobov space with first order mixed derivatives, i.e., $\alpha=1$. In Table 1 we display the generating vector $\boldsymbol{z}$ and the number of points $n$ for the following numerical results. Using given $n$ and $\boldsymbol{z}$, we construct the anti-aliasing set in accordance with Theorem 4 in the following manner: (i) first we generate all integer vector ${\boldsymbol{h}}\in\mathbb{Z}^{d}$ in a bounded region $\|{\boldsymbol{h}}\|\leq R$ for a well chosen $R$; (ii) then we sort the obtained set according to the $\ell_{2}$ distance in ascending order; (iii) we calculate the value $m_{{\boldsymbol{h}}}:={\boldsymbol{h}}\cdot\boldsymbol{z}\bmod{n}$ in the sorted order and store ${\boldsymbol{h}}$ in $\mathcal{A}(\boldsymbol{z},n)$ if the value $m_{{\boldsymbol{h}}}$ has not appeared before. We repeat this step (iii) until we have the full cardinality $|\mathcal{A}(\boldsymbol{z},n)|=n$.

3.1 Convergence with respect to time step size

We consider a common numerical setting as it is considered in [7, 9, 22] where Fourier pseudo-spectral methods are used. We calculate the error with different value of time steps against a reference solution. For the initial condition $g(\boldsymbol{x})$, we choose the Gaussian wave packet given by:

[TABLE]

where the constant $c$ is a normalizing constant to make $\|g\|_{L_{2}}=1$. For the potential function $v$, we consider a smooth potential function

[TABLE]

and a harmonic potential function

[TABLE]

Our aim is to show the temporal discretization error $\|u_{a}(\boldsymbol{x},t)-u_{a}^{m}(\boldsymbol{x})\|_{L_{2}}$ at fixed time $t=m\,{\Delta t}=1$, for that sake we calculate a reference solution $u^{M}_{a}(\boldsymbol{x})$ with the finest time step size ${\Delta t}=1/M=1/10000$, as an approximation of $u_{a}(\boldsymbol{x},t)$. We then vary the time step size ${\Delta t}=1/m=1/5,\ldots,1/1000$ and calculate $u^{m}_{a}(\boldsymbol{x})$ to see the convergence plot of $\|u^{M}_{a}(\boldsymbol{x})-u_{a}^{m}(\boldsymbol{x})\|_{L_{2}}$.

3.2 Sixth-order splitting

We recall that the higher-order exponential splitting is written as

[TABLE]

For the sixth-order method, we employ the coefficients $a_{j}$ and $b_{j}$ from [10] denoted as “s9odr6a” therein. We exhibit the coefficients in Table 2. We plot the results for dimension $2$ to $8$ in Fig 1. The potential $v_{1}$ is not smooth enough on the boundary of $[0,1]^{d}$ so it does not satisfy the required condition in the strict sense. The initial condition $g$ and the potential $v_{2}$ meet all the required conditions. The expected sixth-order convergence is consistent in every plot. When the error reaches to the machine precision, the plot becomes flat. For the $2$-dimensional case with the potential $v_{2}$, we see the convergence happening when the time step size is very small. This can be explained by a phenomenon, called instability of exponential splitting; this is caused by negative coefficients of the exponential splitting $a_{j}$ and $b_{j}$, and is discussed in e.g., [2]. Especially in [2], commutator-free quasi-Magnus exponential integrators are proposed to avoid the issue, however, this is out of the scope of the present paper. The instability issue does not happen in a higher-dimensional settings.

3.3 Eighth-order splitting

For the eighth-order method, we employ the coefficients again from [10] denoted as “s17odr8a”. The coefficients are shown in Table 3. The results are shown in Fig 2 and we again see that the convergence rate is consistently eighth order in each plot. Most of the plot seems to be similar to Fig 1 but with faster convergence, therefore they reach to the machine precision more quickly.

4 Conclusion

We proposed a numerical method to solve the TDSE. With our method using the time step size ${\Delta t}$, the temporal discretization error converges like $\mathcal{O}(({\Delta t})^{p})$ given that the potential function is in Korobov space of smoothness greater than $2p+1/2$. The numerical results (which are performed from $2$ up to $8$ dimensions) confirmed the theory and the rate of error convergence is consistent. By using rank-$1$ lattices, calculations of the time stepping operator and multiplications are efficiently done by only using one-dimensional FFTs.

Pseudo-spectral methods are widely used technique for solving partial differential equations. It is a common choice to use regular grids, but the number of nodes increases exponential with $d$. We have shown an alternative, rank-$1$ lattice pseudo-spectral methods where the number of points can be chosen freely by the user. In combination with higher-order splitting methods, the proposed method solves the TDSE with higher-order convergence in time.