<figure><img src="image_1.png" /></figure>

HAL
open science

# Multiresolution of the one dimensional free-particle propagator

Evgueni Dinvay, Yuliya Zabelina, Luca Frediani

► To cite this version:

Evgueni Dinvay, Yuliya Zabelina, Luca Frediani. Multiresolution of the one dimensional free-particle propagator. 2024. hal-04574531

HAL Id: hal-04574531

https://hal.science/hal-04574531

Preprint submitted on 14 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
<figure><img src="image_1.png" /></figure>

HAL
open science

# Multiresolution of the one dimensional free-particle propagator

Evgueni Dinvay, Yuliya Zabelina, Luca Frediani

► To cite this version:

Evgueni Dinvay, Yuliya Zabelina, Luca Frediani. Multiresolution of the one dimensional free-particle propagator. 2024. hal-04574531

HAL Id: hal-04574531

https://hal.science/hal-04574531v1

Preprint submitted on 14 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# MULTIRESOLUTION OF THE ONE DIMENSIONAL FREE-PARTICLE PROPAGATOR

EVGUENI DINVAY, YULIYA ZABELINA, AND LUCA FREDIANI

ABSTRACT. Novel methods to integrate the time-dependent Schrödinger equation within the framework of multiscale approximation is presented. The methods are based on symplectic splitting algorithms to separate the kinetic and potential parts of the corresponding propagator. The semigroup associated with the free-particle Schrödinger operator is represented in a multiwavelet basis. The propagator is effectively discretized with a contour deformation technique, which overcomes the challenges presented by previous discretization methods. The discretized operator is then employed in simple numerical simulations to test the validity of the implementation and to benchmark its precision.

## 1. INTRODUCTION

Let us consider the one-dimensional time-dependent Schrödinger equation

$$i\partial_t u = -\partial_x^2 u + V(x, t)u, \qquad (1.1)$$

complemented by the initial condition given at $t = 0$ by

$$u(0) = u_0 \in L^2(\mathbb{R}). \qquad (1.2)$$

The space variable $x \in \mathbb{R}$ and the time variable $t > 0$. In general a real-valued potential $V(x, t)$ can depend on the solution, as for example, in water waves and Hartree-Fock theories.

The main objective of this work is to provide an effective numerical discretization of the exponential operator $\exp(it\partial_x^2)$, where $t \in \mathbb{R}$ is a small parameter, not necessarily positive. In what follows we assume that $t \neq 0$. Our goal is to construct a numerical representation of this operator to solve Equation (1.1) adaptively. It is known that a use of multiwavelet bases [6] may lead to adaptive numerical solutions for some partial differential equations [7]. However, the time-evolution operator $\exp(it\partial_x^2)$ falls outside this operator class, that are known to give sparse representations in this type of bases [9, 21]. Nevertheless, it turns out that $\exp(it\partial_x^2)$ exhibits a unique peculiarity that can be leveraged in adaptive calculations, as we demonstrate below.

We recall that the $t$-dependent exponential operator under consideration here forms a semigroup representing solutions of the free particle equation

$$i\partial_t u + \partial_x^2 u = 0. \qquad (1.3)$$

In other words, for any square integrable $u_0 \in L^2(\mathbb{R})$, the function $u = \exp(it\partial_x^2)u_0$ solves Equation 1.3. In numerical simulations, the parameter $t$ usually plays the role of a time (sub)step. The importance of an effective discretization of this exponential operator, at least for a time-independent potential $V$, can be exemplified as follows. Let $A, B$ be time independent linear operators. One of the simplest second-order splitting methods can be expressed as follows

$$e^{At+Bt} = e^{Bt/2}e^{At}e^{Bt/2} + \mathcal{O}(t^3). \qquad (1.4)$$

Taking $A = i\partial_x^2$ and $B = -iV$, one can evolve Equation (1.1) by applying the multiplication operator $e^{Bt/2}$ twice and the pseudo-differential operator $e^{At}$ once at each time step starting from the initial function (1.2). Such operator splitting simplifies the problem (1.1), (1.2) significantly: the exponential of the multiplicative operator is straightforward and the problem reduces to finding an efficient representation for the exponential of the kinetic energy. Below we will also make use of

<footnote>Date: May 13, 2024.</footnote>

<footnote>Key words and phrases.</footnote>
# MULTIRESOLUTION OF THE ONE DIMENSIONAL FREE-PARTICLE PROPAGATOR

EVGUENI DINVAY, YULIYA ZABELINA, AND LUCA FREDIANI

ABSTRACT. Novel methods to integrate the time-dependent Schrödinger equation within the framework of multiscale approximation is presented. The methods are based on symplectic splitting algorithms to separate the kinetic and potential parts of the corresponding propagator. The semigroup associated with the free-particle Schrödinger operator is represented in a multiwavelet basis. The propagator is effectively discretized with a contour deformation technique, which overcomes the challenges presented by previous discretization methods. The discretized operator is then employed in simple numerical simulations to test the validity of the implementation and to benchmark its precision.

## 1. INTRODUCTION

Let us consider the one-dimensional time-dependent Schrödinger equation

$$i\partial_t u = -\partial_x^2 u + V(x, t)u, \qquad (1.1)$$

complemented by the initial condition given at $t = 0$ by

$$u(0) = u_0 \in L^2(\mathbb{R}). \qquad (1.2)$$

The space variable $x \in \mathbb{R}$ and the time variable $t > 0$. In general a real-valued potential $V(x, t)$ can depend on the solution, as for example, in water waves and Hartree-Fock theories.

The main objective of this work is to provide an effective numerical discretization of the exponential operator $\exp(it\partial_x^2)$, where $t \in \mathbb{R}$ is a small parameter, not necessarily positive. In what follows we assume that $t \neq 0$. Our goal is to construct a numerical representation of this operator to solve Equation (1.1) adaptively. It is known that a use of multiwavelet bases [6] may lead to adaptive numerical solutions for some partial differential equations [7]. However, the time-evolution operator $\exp(it\partial_x^2)$ falls outside this operator class, that are known to give sparse representations in this type of bases [9, 21]. Nevertheless, it turns out that $\exp(it\partial_x^2)$ exhibits a unique peculiarity that can be leveraged in adaptive calculations, as we demonstrate below.

We recall that the $t$-dependent exponential operator under consideration here forms a semigroup representing solutions of the free particle equation

$$i\partial_t u + \partial_x^2 u = 0. \qquad (1.3)$$

In other words, for any square integrable $u_0 \in L^2(\mathbb{R})$, the function $u = \exp(it\partial_x^2)u_0$ solves Equation 1.3. In numerical simulations, the parameter $t$ usually plays the role of a time (sub)step. The importance of an effective discretization of this exponential operator, at least for a time-independent potential $V$, can be exemplified as follows. Let $A, B$ be time independent linear operators. One of the simplest second-order splitting methods can be expressed as follows

$$e^{At+Bt} = e^{Bt/2}e^{At}e^{Bt/2} + \mathcal{O}(t^3). \qquad (1.4)$$

Taking $A = i\partial_x^2$ and $B = -iV$, one can evolve Equation (1.1) by applying the multiplication operator $e^{Bt/2}$ twice and the pseudo-differential operator $e^{At}$ once at each time step starting from the initial function (1.2). Such operator splitting simplifies the problem (1.1), (1.2) significantly: the exponential of the multiplicative operator is straightforward and the problem reduces to finding an efficient representation for the exponential of the kinetic energy. Below we will also make use of

<footnote>Date: May 13, 2024.</footnote>

<footnote>Key words and phrases.</footnote>
higher order schemes than (1.4). They will permit larger time steps, which turns out to be especially crucial for multiwavelet representation of the semigroup $\exp (it\partial_x^2)$, due to some limitations on how small a time step $t$ might be taken.

We are motivated by the desire of simulating attosecond electron dynamics. In recent years there has been a growing interest in studying the electronic structure of molecules under strong and ultra-fast laser pulses [30]. This is clearly testified by the 2018 and 2023 Nobel prizes in Physics[29, 33]. This could potentially lead to reaction control by attosecond laser pulses in the future. From a computational perspective, electron dynamics in atoms and molecules is challenging: standard atomic orbital approaches routinely used in Quantum Chemistry cannot be applied in strong laser pulses, as simple numerical experiments suggest that ionization plays a significant role and cannot be neglected [17]. In other words, the coupling with continuum or unbound states complicates the use of atomic orbitals for the problem discretization. Inevitably it demands the use of grid based methods, such as finite elements, for example.

In the case of hydrogen atom subjected to a laser field we have the following three-dimensional Schrödinger equation

$$i\partial_t\psi = -\frac{1}{2}\Delta\psi - \frac{1}{|x|}\psi + \mathcal{E}(x, t)\psi, \qquad (1.5)$$

where the second term $-|x|^{-1}$ stands for the Coulomb potential and $\mathcal{E}(x, t)$ stands for the external electric field generated by a laser. It is normally complemented by the ground state $\psi_0(x) = e^{-|x|}/\sqrt{\pi}$ serving as an initial condition $\psi(x, 0) = \psi_0(x)$.

It can be advantageous to make use of adaptive methods to describe the singular Coulomb interaction between the electron and the nucleus, as well as the cusp in the eigenfunctions which follows from it. This, however, sets some limitations on the method development. Additionally, for many-body systems (atoms and molecules), the potential becomes nonlinear due to the presence of electron-electron interactions. The corresponding nonlinearities are normally presented in equations as singular integral convolutions. One such method is constituted by a multiwavelet (MW) representation, within the framework of multiresolution analysis [7]. This approach has demonstrated excellent results in achieving fast algorithms and high precision for static quantum chemistry problems [23, 24, 26]. However, discussing these static results is beyond the scope of the current contribution. Furthermore, the research field is rapidly evolving, and we recommend referring to a comprehensive review paper for a broader perspective [10]. An overview of the most recent advancements will also be available in an upcoming paper [37].

We aim to extend the multiwavelet approach to address time-dependent scenarios. The first step in integrating the generic time-dependent Schrödinger equation, similar to (1.5), involves discretizing the semigroup $e^{it\Delta}$. Due to the limitations highlighted above, the fast Fourier transform (FFT) cannot be relied upon for quantum chemistry applications. This is why we find the multiscale approach promising. To our knowledge, the only work on multiwavelet representation of the dynamical Schrödinger type equations was conducted by Vence et al. [38]. They view the semigroup as a convolution and mitigate the errors arising from its oscillatory behavior by damping high-frequency oscillations in Fourier space. The resulting kernel approximation can be seamlessly integrated into the multiwavelet framework developed for static problems.

In contrast, we adopt a different approach by leveraging the smoothing property of the semigroup instead of damping it. It enables us to achieve a highly accurate multiscale approximation of the free-particle propagator. Moreover, it turns out that this representation can be regarded as sparse, as long as particular restrictions are met during numerical simulations. Since we do not treat the semigroup as a convolution suitable for chemistry applications, unlike the mentioned above method employed by Vence et al., further development is necessary to fully exploit our representation for studying electron dynamics in molecules. Therefore, this article focuses exclusively on the one-dimensional case. In the future, we plan to extend our results and simulate the three-dimensional quantum mechanical problems.

The paper is organised as follows. In Section 2 we introduce some preliminary notions and recall the theory behind multiwavelet bases. In Section 3 we lay out the problem with some insights
higher order schemes than (1.4). They will permit larger time steps, which turns out to be especially crucial for multiwavelet representation of the semigroup $\exp (it\partial_x^2)$, due to some limitations on how small a time step $t$ might be taken.

We are motivated by the desire of simulating attosecond electron dynamics. In recent years there has been a growing interest in studying the electronic structure of molecules under strong and ultra-fast laser pulses [30]. This is clearly testified by the 2018 and 2023 Nobel prizes in Physics [29, 33]. This could potentially lead to reaction control by attosecond laser pulses in the future. From a computational perspective, electron dynamics in atoms and molecules is challenging: standard atomic orbital approaches routinely used in Quantum Chemistry cannot be applied in strong laser pulses, as simple numerical experiments suggest that ionization plays a significant role and cannot be neglected [17]. In other words, the coupling with continuum or unbound states complicates the use of atomic orbitals for the problem discretization. Inevitably it demands the use of grid based methods, such as finite elements, for example.

In the case of hydrogen atom subjected to a laser field we have the following three-dimensional Schrödinger equation

$$i\partial_t\psi = -\frac{1}{2}\Delta\psi - \frac{1}{|x|}\psi + \mathcal{E}(x,t)\psi, \quad (1.5)$$

where the second term $-|x|^{-1}$ stands for the Coulomb potential and $\mathcal{E}(x,t)$ stands for the external electric field generated by a laser. It is normally complemented by the ground state $\psi_0(x) = e^{-|x|}/\sqrt{\pi}$ serving as an initial condition $\psi(x, 0) = \psi_0(x)$.

It can be advantageous to make use of adaptive methods to describe the singular Coulomb interaction between the electron and the nucleus, as well as the cusp in the eigenfunctions which follows from it. This, however, sets some limitations on the method development. Additionally, for many-body systems (atoms and molecules), the potential becomes nonlinear due to the presence of electron-electron interactions. The corresponding nonlinearities are normally presented in equations as singular integral convolutions. One such method is constituted by a multiwavelet (MW) representation, within the framework of multiresolution analysis [7]. This approach has demonstrated excellent results in achieving fast algorithms and high precision for static quantum chemistry problems [23, 24, 26]. However, discussing these static results is beyond the scope of the current contribution. Furthermore, the research field is rapidly evolving, and we recommend referring to a comprehensive review paper for a broader perspective [10]. An overview of the most recent advancements will also be available in an upcoming paper [37].

We aim to extend the multiwavelet approach to address time-dependent scenarios. The first step in integrating the generic time-dependent Schrödinger equation, similar to (1.5), involves discretizing the semigroup $e^{it\Delta}$. Due to the limitations highlighted above, the fast Fourier transform (FFT) cannot be relied upon for quantum chemistry applications. This is why we find the multiscale approach promising. To our knowledge, the only work on multiwavelet representation of the dynamical Schrödinger type equations was conducted by Vence et al. [38]. They view the semigroup as a convolution and mitigate the errors arising from its oscillatory behavior by damping high-frequency oscillations in Fourier space. The resulting kernel approximation can be seamlessly integrated into the multiwavelet framework developed for static problems.

In contrast, we adopt a different approach by leveraging the smoothing property of the semigroup instead of damping it. It enables us to achieve a highly accurate multiscale approximation of the free-particle propagator. Moreover, it turns out that this representation can be regarded as sparse, as long as particular restrictions are met during numerical simulations. Since we do not treat the semigroup as a convolution suitable for chemistry applications, unlike the mentioned above method employed by Vence et al., further development is necessary to fully exploit our representation for studying electron dynamics in molecules. Therefore, this article focuses exclusively on the one-dimensional case. In the future, we plan to extend our results and simulate the three-dimensional quantum mechanical problems.

The paper is organised as follows. In Section 2 we introduce some preliminary notions and recall the theory behind multiwavelet bases. In Section 3 we lay out the problem with some insights
into the difficulties that we aim to overcome. In Section 4, the multiresolution representation of $\exp (it\partial_x^2)$ is reduced to the evaluation of a special kind of integrals. In Section 4 these integrals are evaluated by deforming the integration contour in a way that allows to exploit the smoothing property of this exponential operator. Furthermore, we demonstrate that by carefully selecting the time parameter $t$ and polynomial approximation order, its multiresolution representation will display an unusual sparsity pattern. More specifically, the matrices used in calculations have non-zero elements away from their main diagonals. This sparsity depends significantly on the value of $t$, necessitating the use of high-order numerical schemes for accurate wave propagation simulations. To address this requirement, we introduce high-order symplectic integrators in Section 5. Finally, we test our numerical representation of $\exp (it\partial_x^2)$ by running numerical simulations for toy benchmark models in Section 6. A summary is provided in Section 7 with some plans for future works.

## 2. MULTIWAVELET BASES

This section provides a concise theoretical background on multiwavelet bases. We are interested in approximations of functions in $L^2(\mathbb{R})$. However, for practical reasons, we have to restrict ourselves to a sufficiently large interval. Equation (1.1) can always be normalised in such a way that the unit computational domain $[0, 1]$ will suit our needs, while keeping the constant in front of the second derivative to be $-1$. Note that any function $u \in L^2(0, 1)$ can be trivially extended to a function in $L^2(\mathbb{R})$, by setting it to zero outside $(0, 1)$. On the contrary, any function from $L^2(\mathbb{R})$ can be projected into $L^2(0, 1)$, by multiplying it with the characteristic function of the unit interval $(0, 1)$. This remark allows us to make sense of the inclusion $L^2(0, 1) \subset L^2(\mathbb{R})$. Now for any $\mathfrak{k} \in \mathbb{N} = \{1, 2, \dots\}$ and any $n \in \mathbb{N}_0 = \{0, 1, 2, \dots\}$ we introduce a space of piecewise polynomials $V_n^{\mathfrak{k}}$ as follows. A function $f \in V_n^{\mathfrak{k}}$ provided that on each dyadic interval $(l/2^n, (l+1)/2^n)$ with $l = 0, 1, \dots, 2^n - 1$ it is a polynomial of degree less than $\mathfrak{k}$ and it is zero elsewhere. Note that the space $V_n^{\mathfrak{k}}$ has dimension $2^n\mathfrak{k}$ and

$$V_0^{\mathfrak{k}} \subset V_1^{\mathfrak{k}} \subset \dots \subset V_n^{\mathfrak{k}} \subset \dots \subset L^2(0, 1) \subset L^2(\mathbb{R}),$$

which defines a multiresolution analysis (MRA) [27, 28], since the union $\cup_{n=0}^{\infty} V_n^{\mathfrak{k}}$ is dense in $L^2(0, 1)$. We refer to the given fixed number $\mathfrak{k}$ as the order of MRA and to the integer variable $n$ as a scaling level.

Let $\phi_0, \dots, \phi_{\mathfrak{k}-1}$ be a basis of $V_0^{\mathfrak{k}}$ that will be called scaling functions. The space $V_n^{\mathfrak{k}}$ is spanned by $2^n\mathfrak{k}$ functions which are obtained from the scaling functions by dilation and translation

$$\phi_{j,l}^{n}(x) = 2^{n/2} \phi_j(2^n x - l), \quad j = 0, \dots, \mathfrak{k}-1, \quad l = 0, \dots, 2^n - 1. \qquad (2.1)$$

There is some freedom in choosing a basis in $V_0^{\mathfrak{k}}$, at least for $\mathfrak{k} > 1$. The first example appeared in [6], and it consists of the Legendre scaling functions

$$\phi_j(x) = \sqrt{2j+1} P_j(2x-1), \quad x \in (0,1), \quad j = 0, \dots, \mathfrak{k}-1, \qquad (2.2)$$

where $P_j$ are standard Legendre polynomials. Outside the unit interval $(0, 1)$ they are set to zero. An alternative basis, proved to be more efficient that is especially crucial in multi-dimensional numerical calculations, was presented in [7]. It consists of the interpolating scaling functions

$$\varphi_j(x) = \sqrt{w_j} \sum_{m=0}^{\mathfrak{k}-1} \phi_m(x_j) \phi_m(x), \quad j = 0, \dots, \mathfrak{k}-1, \qquad (2.3)$$

where $\phi_m$ are the Legendre scaling functions (2.2). Here $x_0, \dots, x_{\mathfrak{k}-1}$ denote the roots of $P_{\mathfrak{k}}(2x-1)$ and the quadrature weights $w_j = 1/(\mathfrak{k}P'_{\mathfrak{k}}(2x_j - 1)P_{\mathfrak{k}-1}(2x_j - 1))$. The expansion coefficients of a general function $f$ in a scaling basis are the integrals

$$s_{j,l}^{n}(f) = \int_{\mathbb{R}} f(x) \phi_{j,l}^{n}(x) dx = 2^{-n/2} \int_{0}^{1} f(2^{-n}(x+l)) \phi_{j}(x) dx, \qquad (2.4)$$
into the difficulties that we aim to overcome. In Section 4, the multiresolution representation of $\exp (it\partial_x^2)$ is reduced to the evaluation of a special kind of integrals. In Section 4 these integrals are evaluated by deforming the integration contour in a way that allows to exploit the smoothing property of this exponential operator. Furthermore, we demonstrate that by carefully selecting the time parameter $t$ and polynomial approximation order, its multiresolution representation will display an unusual sparsity pattern. More specifically, the matrices used in calculations have non-zero elements away from their main diagonals. This sparsity depends significantly on the value of $t$, necessitating the use of high-order numerical schemes for accurate wave propagation simulations. To address this requirement, we introduce high-order symplectic integrators in Section 5. Finally, we test our numerical representation of $\exp (it\partial_x^2)$ by running numerical simulations for toy benchmark models in Section 6. A summary is provided in Section 7 with some plans for future works.

## 2. MULTIWAVELET BASES

This section provides a concise theoretical background on multiwavelet bases. We are interested in approximations of functions in $L^2(\mathbb{R})$. However, for practical reasons, we have to restrict ourselves to a sufficiently large interval. Equation (1.1) can always be normalised in such a way that the unit computational domain $[0, 1]$ will suit our needs, while keeping the constant in front of the second derivative to be $-1$. Note that any function $u \in L^2(0, 1)$ can be trivially extended to a function in $L^2(\mathbb{R})$, by setting it to zero outside $(0, 1)$. On the contrary, any function from $L^2(\mathbb{R})$ can be projected into $L^2(0, 1)$, by multiplying it with the characteristic function of the unit interval $(0, 1)$. This remark allows us to make sense of the inclusion $L^2(0, 1) \subset L^2(\mathbb{R})$. Now for any $\mathfrak{k} \in \mathbb{N} = \{1, 2, \dots\}$ and any $n \in \mathbb{N}_0 = \{0, 1, 2, \dots\}$ we introduce a space of piecewise polynomials $V_n^{\mathfrak{k}}$ as follows. A function $f \in V_n^{\mathfrak{k}}$ provided that on each dyadic interval $(l/2^n, (l+1)/2^n)$ with $l = 0, 1, \dots, 2^n - 1$ it is a polynomial of degree less than $\mathfrak{k}$ and it is zero elsewhere. Note that the space $V_n^{\mathfrak{k}}$ has dimension $2^n\mathfrak{k}$ and

$$V_0^{\mathfrak{k}} \subset V_1^{\mathfrak{k}} \subset \dots \subset V_n^{\mathfrak{k}} \subset \dots \subset L^2(0, 1) \subset L^2(\mathbb{R}),$$

which defines a multiresolution analysis (MRA) [27, 28], since the union $\cup_{n=0}^{\infty} V_n^{\mathfrak{k}}$ is dense in $L^2(0, 1)$. We refer to the given fixed number $\mathfrak{k}$ as the order of MRA and to the integer variable $n$ as a scaling level.

Let $\phi_0, \dots, \phi_{\mathfrak{k}-1}$ be a basis of $V_0^{\mathfrak{k}}$ that will be called scaling functions. The space $V_n^{\mathfrak{k}}$ is spanned by $2^n\mathfrak{k}$ functions which are obtained from the scaling functions by dilation and translation

$$\phi_{j,l}^{n}(x) = 2^{n/2}\phi_j(2^n x - l), \quad j = 0, \dots, \mathfrak{k}-1, \quad l = 0, \dots, 2^n - 1. \qquad (2.1)$$

There is some freedom in choosing a basis in $V_0^{\mathfrak{k}}$, at least for $\mathfrak{k} > 1$. The first example appeared in [6], and it consists of the Legendre scaling functions

$$\phi_j(x) = \sqrt{2j+1} P_j(2x-1), \quad x \in (0,1), \quad j = 0, \dots, \mathfrak{k}-1, \qquad (2.2)$$

where $P_j$ are standard Legendre polynomials. Outside the unit interval $(0, 1)$ they are set to zero. An alternative basis, proved to be more efficient that is especially crucial in multi-dimensional numerical calculations, was presented in [7]. It consists of the interpolating scaling functions

$$\varphi_j(x) = \sqrt{w_j} \sum_{m=0}^{\mathfrak{k}-1} \phi_m(x_j)\phi_m(x), \quad j = 0, \dots, \mathfrak{k}-1, \qquad (2.3)$$

where $\phi_m$ are the Legendre scaling functions (2.2). Here $x_0, \dots, x_{\mathfrak{k}-1}$ denote the roots of $P_{\mathfrak{k}}(2x-1)$ and the quadrature weights $w_j = 1/(\mathfrak{k}P'_{\mathfrak{k}}(2x_j - 1)P_{\mathfrak{k}-1}(2x_j - 1))$. The expansion coefficients of a general function $f$ in a scaling basis are the integrals

$$s_{jl}^{n}(f) = \int_{\mathbb{R}} f(x) \phi_{jl}^{n}(x) dx = 2^{-n/2} \int_{0}^{1} f(2^{-n}(x+l)) \phi_{j}(x) dx, \qquad (2.4)$$
that can be evaluated numerically using the Gauss-Legendre quadrature

$$s_{jl}^{n}(f) \approx 2^{-n/2} \sum_{m=0}^{\ell-1} w_m f(2^{-n}(x_m + l))\phi_j(x_m), \quad j = 0, \dots, \ell-1, \quad l = 0, \dots, 2^n - 1. \quad (2.5)$$

Note that the interpolating scaling functions (2.3) satisfy $\phi_j(x_m) = w_j^{-1/2}\delta_{jm}$ and so (2.5) simplifies to

$$s_{jl}^{n}(f) \approx 2^{-n/2} \sqrt{w_j} f(2^{-n}(x_m + l)) \quad (2.6)$$

for this particular choice of a basis.

The multiwavelet space $W_n^\ell$ is defined as the orthogonal complement of $V_n^\ell$ in $V_{n+1}^\ell$, so we have

$$V_{n+1}^{\ell} = V_{n}^{\ell} \oplus W_{n}^{\ell} = V_{0}^{\ell} \oplus W_{0}^{\ell} \oplus W_{1}^{\ell} \oplus \dots \oplus W_{n}^{\ell}.$$

As above it is enough to define a basis $\psi_0, \dots, \psi_{\ell-1}$ in $W_0^\ell$ and then transform it in line with the general rule (2.1), in order to get an orthonormal basis in each $W_n^\ell$. The construction of a specific multiwavelet basis that satisfies certain given restrictions is more involved. We are using the one that was constructed in [6], since it provides us with an additional vanishing moment property, namely,

$$\int_0^1 \psi_j(x) x^m dx = 0, \quad m = 0, 1, \dots, j + \ell - 1. \quad (2.7)$$

The multiwavelet expansion coefficients of a general function $f$ are the integrals

$$d_{jl}^{n}(f) = \int_{\mathbb{R}} f(x) \psi_{jl}^{n}(x) dx, \quad (2.8)$$

that could in principal be also evaluated with the help of the Gauss-Legendre quadrature. However, there are precision issues connected to this approach. The multiwavelet coefficients (2.8) are instead obtained from the higher level scaling coefficients $s_{jl}^{n+1}(f)$, where the latter are calculated by the quadrature rule (2.5). This is possible thanks to the so called forward wavelet transform [18]

$$\begin{pmatrix} \phi_l^n \\ \psi_l^n \end{pmatrix} = U \begin{pmatrix} \phi_{2l}^{n+1} \\ \phi_{2l+1}^{n+1} \end{pmatrix} \quad (2.9)$$

where the two bases of $V_{n+1}^\ell$ are grouped in accordance with the short notation agreement

$$\phi_l^n = \begin{pmatrix} \phi_{0,l}^n(x) \\ \dots \\ \phi_{\ell-1,l}^n(x) \end{pmatrix}, \quad \psi_l^n = \begin{pmatrix} \psi_{0,l}^n(x) \\ \dots \\ \psi_{\ell-1,l}^n(x) \end{pmatrix}$$

and $U$ stands for the unitary matrix

$$U = \begin{pmatrix} H^{(0)} & H^{(1)} \\ G^{(0)} & G^{(1)} \end{pmatrix} \quad (2.10)$$

consisting of $\ell \times \ell$-size filter blocks $H^{(0)} = (h_{pj}^{(0)})$, $\dots$, $G^{(1)} = (g_{pj}^{(1)})$ depending only on the type of scaling basis (Legendre or interpolating) in use [7]. Integrating (2.9) together with $f$, we immediately obtain the following relation

$$\begin{pmatrix} s_l^n \\ d_l^n \end{pmatrix} = U \begin{pmatrix} s_{2l}^{n+1} \\ s_{2l+1}^{n+1} \end{pmatrix} \quad (2.11)$$

with vectors $s_l^n, d_l^n$ consisting of coefficients $s_{jl}^n, d_{jl}^n, j = 0, \dots, \ell - 1$, accordingly.

There are two main advantages with multiwavelets, compared to other discretization methods: (1) the disjoint support of the basis functions combined with the vanishing moments of the wavelets enables sparse and adaptive function representations, with a rigorous precision control; (2) some pseudo-differential operators have sparse multiwavelet representations [6], leading to fast (ideally linearly scaling) algorithms.

Let $P^n, Q^n$ be the orthogonal projectors of $L^2(\mathbb{R})$ into $V_n^\ell, W_n^\ell$, respectively. Clearly, $Q^n = P^{n+1} - P^n$. Moreover, $P^n$ converges strongly to the projection $P: L^2(\mathbb{R}) \to L^2(0, 1)$, whereas $Q^n$
that can be evaluated numerically using the Gauss-Legendre quadrature

$$s_{jl}^{n}(f) \approx 2^{-n/2} \sum_{m=0}^{\ell-1} w_{m} f(2^{-n}(x_{m} + l)) \phi_{j}(x_{m}), \quad j = 0, \dots, \ell-1, \quad l = 0, \dots, 2^{n} - 1. \quad (2.5)$$

Note that the interpolating scaling functions (2.3) satisfy $\phi_j(x_m) = w_j^{-1/2}\delta_{jm}$ and so (2.5) simplifies to

$$s_{jl}^{n}(f) \approx 2^{-n/2} \sqrt{w_j} f(2^{-n}(x_m + l)) \quad (2.6)$$

for this particular choice of a basis.

The multiwavelet space $W_n^\ell$ is defined as the orthogonal complement of $V_n^\ell$ in $V_{n+1}^\ell$, so we have

$$V_{n+1}^{\ell} = V_{n}^{\ell} \oplus W_{n}^{\ell} = V_{0}^{\ell} \oplus W_{0}^{\ell} \oplus W_{1}^{\ell} \oplus \dots \oplus W_{n}^{\ell}.$$

As above it is enough to define a basis $\psi_0, \dots, \psi_{\ell-1}$ in $W_0^\ell$ and then transform it in line with the general rule (2.1), in order to get an orthonormal basis in each $W_n^\ell$. The construction of a specific multiwavelet basis that satisfies certain given restrictions is more involved. We are using the one that was constructed in [6], since it provides us with an additional vanishing moment property, namely,

$$\int_{0}^{1} \psi_{j}(x) x^{m} dx = 0, \quad m = 0, 1, \dots, j + \ell - 1. \quad (2.7)$$

The multiwavelet expansion coefficients of a general function $f$ are the integrals

$$d_{jl}^{n}(f) = \int_{\mathbb{R}} f(x) \psi_{jl}^{n}(x) dx, \quad (2.8)$$

that could in principal be also evaluated with the help of the Gauss-Legendre quadrature. However, there are precision issues connected to this approach. The multiwavelet coefficients (2.8) are instead obtained from the higher level scaling coefficients $s_{jl}^{n+1}(f)$, where the latter are calculated by the quadrature rule (2.5). This is possible thanks to the so called forward wavelet transform [18]

$$\begin{pmatrix} \phi_l^n \\ \psi_l^n \end{pmatrix} = U \begin{pmatrix} \phi_{2l}^{n+1} \\ \phi_{2l+1}^{n+1} \end{pmatrix} \quad (2.9)$$

where the two bases of $V_{n+1}^\ell$ are grouped in accordance with the short notation agreement

$$\phi_l^n = \begin{pmatrix} \phi_{0,l}^n(x) \\ \vdots \\ \phi_{\ell-1,l}^n(x) \end{pmatrix}, \quad \psi_l^n = \begin{pmatrix} \psi_{0,l}^n(x) \\ \vdots \\ \psi_{\ell-1,l}^n(x) \end{pmatrix}$$

and $U$ stands for the unitary matrix

$$U = \begin{pmatrix} H^{(0)} & H^{(1)} \\ G^{(0)} & G^{(1)} \end{pmatrix} \quad (2.10)$$

consisting of $\ell \times \ell$-size filter blocks $H^{(0)} = (h_{pj}^{(0)})$, $\dots$, $G^{(1)} = (g_{pj}^{(1)})$ depending only on the type of scaling basis (Legendre or interpolating) in use [7]. Integrating (2.9) together with $f$, we immediately obtain the following relation

$$\begin{pmatrix} s_l^n \\ d_l^n \end{pmatrix} = U \begin{pmatrix} s_{2l}^{n+1} \\ s_{2l+1}^{n+1} \end{pmatrix} \quad (2.11)$$

with vectors $s_l^n, d_l^n$ consisting of coefficients $s_{jl}^n, d_{jl}^n, j = 0, \dots, \ell-1$, accordingly.

There are two main advantages with multiwavelets, compared to other discretization methods: (1) the disjoint support of the basis functions combined with the vanishing moments of the wavelets enables sparse and adaptive function representations, with a rigorous precision control; (2) some pseudo-differential operators have sparse multiwavelet representations [6], leading to fast (ideally linearly scaling) algorithms.

Let $P^n, Q^n$ be the orthogonal projectors of $L^2(\mathbb{R})$ into $V_n^\ell, W_n^\ell$, respectively. Clearly, $Q^n = P^{n+1} - P^n$. Moreover, $P^n$ converges strongly to the projection $P: L^2(\mathbb{R}) \to L^2(0, 1)$, whereas $Q^n$
tends to zero. Then any bounded operator $\mathcal{T}$ can be represented on the computational domain $[0, 1]$ as

$$ P\mathcal{T}P = \lim_{n \to \infty} \mathcal{T}_n = \mathcal{T}_0 + \sum_{n=0}^{\infty} (\mathcal{A}_n + \mathcal{B}_n + \mathcal{C}_n), \qquad (2.12) $$

where we introduced the multiresolution restrictions

$$ \mathcal{A}_n = Q^n \mathcal{T} Q^n, \quad \mathcal{B}_n = Q^n \mathcal{T} P^n, \quad \mathcal{C}_n = P^n \mathcal{T} Q^n, \quad \mathcal{T}_n = P^n \mathcal{T} P^n. \qquad (2.13) $$

In other words, each approximation $\mathcal{T}_n$ of $\mathcal{T}$ on $[0, 1]$ is defined by its restriction $\mathcal{T}_0$ in $V_0^\mathcal{t}$ and by the collection of triplets $(\mathcal{A}_0, \mathcal{B}_0, \mathcal{C}_0), \dots, (\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n)$. We refer to this representation as the nonstandard form and write shortly

$$ \mathcal{T}_n \cong \{\mathcal{T}_0, (\mathcal{A}_0, \mathcal{B}_0, \mathcal{C}_0), \dots, (\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n)\}. \qquad (2.14) $$

This form is very useful for some classes of operators. In particular, Calderón-Zygmund operators have sparse representations for the matrices associated with $\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n$, leading to fast effective algorithms [8]. Another advantage of the nonstandard form is the absence of coupling between scales when the operator is applied (such a coupling results from a post-processing step which relies only on the fast multiwavelet transform [8, 26]). This is a key-feature as one wants to preserve the adaptivity of functions while the operator is applied.

Before proceeding to concrete examples of operators $\mathcal{T}$, we introduce the following matrices, defined in the $V_n^\mathcal{t} \oplus W_n^\mathcal{t}$ basis, as

$$ \begin{align} [\sigma_{l'l}^n]_{j'j} &= \int_0^1 \phi_{j'l'}^n(x) \mathcal{T} \phi_{jl}^n(x) dx, \qquad [\gamma_{l'l}^n]_{j'j} = \int_0^1 \phi_{j'l'}^n(x) \mathcal{T} \psi_{jl}^n(x) dx, \tag{2.15} \\ [\beta_{l'l}^n]_{j'j} &= \int_0^1 \psi_{j'l'}^n(x) \mathcal{T} \phi_{jl}^n(x) dx, \qquad [\alpha_{l'l}^n]_{j'j} = \int_0^1 \psi_{j'l'}^n(x) \mathcal{T} \psi_{jl}^n(x) dx. \notag \end{align} $$

For every scale $n \in \mathbb{N}_0$, symbol $\sigma^n$ stands for a table of $2^n \times 2^n$-size, where each element is itself a matrix of $\mathcal{t} \times \mathcal{t}$-size. Similarly, the blocks $\alpha^n, \beta^n$ and $\gamma^n$, associated with the triple $(\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n)$, have the same forms. Making use of the wavelet transform (2.9) one can easily deduce the following relation.

$$ \begin{pmatrix} \sigma_{l'l}^n & \gamma_{l'l}^n \\ \beta_{l'l}^n & \alpha_{l'l}^n \end{pmatrix} = U \begin{pmatrix} \sigma_{2l',2l}^{n+1} & \sigma_{2l',2l+1}^{n+1} \\ \sigma_{2l'+1,2l}^{n+1} & \sigma_{2l'+1,2l+1}^{n+1} \end{pmatrix} U^T. \qquad (2.16) $$

Finally, we can describe a practical algorithm for the operator application in the nonstandard form. Let us for a given arbitrary $f \in L^2(0, 1)$ consider the sequences $\tilde{d}_l^n, \tilde{s}_l^n, \tilde{s}_l^n$ with $n \in \mathbb{N}_0, l = 0, 1, \dots, 2^n - 1$ of elements from $\mathbb{C}^\mathcal{t}$ defined by the following iterative procedure

$$ \tilde{d}_l^n = \sum_{l=0}^{2^n-1} (\alpha_{l'l}^n d_l^n(f) + \beta_{l'l}^n s_l^n(f)), \quad \tilde{s}_l^n = \sum_{l=0}^{2^n-1} \gamma_{l'l}^n d_l^n(f) $$

and

$$ \begin{align*} \tilde{s}_{2m}^n &= H^{(0)} (\tilde{s}_m^{n-1} + \tilde{s}_m^{n-1}) + G^{(0)} \tilde{d}_m^{n-1}, & n \ge 1, \quad \text{where } \tilde{s}_0^0 = \sigma_{00}^0 s_0^0(f). \\ \tilde{s}_{2m+1}^n &= H^{(1)} (\tilde{s}_m^{n-1} + \tilde{s}_m^{n-1}) + G^{(1)} \tilde{d}_m^{n-1}, & n \ge 1, \quad \text{where } \tilde{s}_0^1 = \sigma_{01}^0 s_0^1(f). \end{align*} $$

**Lemma 1.** Let $f \in L^2(0, 1)$, $n \in \mathbb{N}_0$ and $g = \mathcal{T}_{n+1}f$. Then its coefficients can be obtained by the above iterative procedure as

$$ d_l^n(g) = \tilde{d}_l^n, \quad s_l^n(g) = \tilde{s}_l^n + \tilde{s}_l^n \qquad (2.17) $$

at the finest scale $n$. The rest coefficients, at scales $n = 0, \dots, n-1$, are found by (2.11).

*Proof.* It is enough to prove (2.17), since the scale relation (2.11) is already known. Moreover, the first equality in (2.17) is obvious. Now we can notice that the second equality in (2.17) clearly
tends to zero. Then any bounded operator $\mathcal{T}$ can be represented on the computational domain $[0, 1]$ as

$$ P\mathcal{T}P = \lim_{n \to \infty} \mathcal{T}_n = \mathcal{T}_0 + \sum_{n=0}^{\infty} (\mathcal{A}_n + \mathcal{B}_n + \mathcal{C}_n), \qquad (2.12) $$

where we introduced the multiresolution restrictions

$$ \mathcal{A}_n = Q^n \mathcal{T} Q^n, \quad \mathcal{B}_n = Q^n \mathcal{T} P^n, \quad \mathcal{C}_n = P^n \mathcal{T} Q^n, \quad \mathcal{T}_n = P^n \mathcal{T} P^n. \qquad (2.13) $$

In other words, each approximation $\mathcal{T}_n$ of $\mathcal{T}$ on $[0, 1]$ is defined by its restriction $\mathcal{T}_0$ in $V_0^\mathcal{t}$ and by the collection of triplets $(\mathcal{A}_0, \mathcal{B}_0, \mathcal{C}_0), \dots, (\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n)$. We refer to this representation as the nonstandard form and write shortly

$$ \mathcal{T}_n \cong \{\mathcal{T}_0, (\mathcal{A}_0, \mathcal{B}_0, \mathcal{C}_0), \dots, (\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n)\}. \qquad (2.14) $$

This form is very useful for some classes of operators. In particular, Calderón-Zygmund operators have sparse representations for the matrices associated with $\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n$, leading to fast effective algorithms [8]. Another advantage of the nonstandard form is the absence of coupling between scales when the operator is applied (such a coupling results from a post-processing step which relies only on the fast multiwavelet transform [8, 26]). This is a key-feature as one wants to preserve the adaptivity of functions while the operator is applied.

Before proceeding to concrete examples of operators $\mathcal{T}$, we introduce the following matrices, defined in the $V_n^\mathcal{t} \oplus W_n^\mathcal{t}$ basis, as

$$ \begin{align} [\sigma_{l'l}^n]_{j'j} &= \int_0^1 \phi_{j'l'}^n(x) \mathcal{T} \phi_{jl}^n(x) dx, \qquad [\gamma_{l'l}^n]_{j'j} = \int_0^1 \phi_{j'l'}^n(x) \mathcal{T} \psi_{jl}^n(x) dx, \tag{2.15} \\ [\beta_{l'l}^n]_{j'j} &= \int_0^1 \psi_{j'l'}^n(x) \mathcal{T} \phi_{jl}^n(x) dx, \qquad [\alpha_{l'l}^n]_{j'j} = \int_0^1 \psi_{j'l'}^n(x) \mathcal{T} \psi_{jl}^n(x) dx. \notag \end{align} $$

For every scale $n \in \mathbb{N}_0$, symbol $\sigma^n$ stands for a table of $2^n \times 2^n$-size, where each element is itself a matrix of $\mathcal{t} \times \mathcal{t}$-size. Similarly, the blocks $\alpha^n, \beta^n$ and $\gamma^n$, associated with the triple $(\mathcal{A}_n, \mathcal{B}_n, \mathcal{C}_n)$, have the same forms. Making use of the wavelet transform (2.9) one can easily deduce the following relation.

$$ \begin{pmatrix} \sigma_{l'l}^n & \gamma_{l'l}^n \\ \beta_{l'l}^n & \alpha_{l'l}^n \end{pmatrix} = U \begin{pmatrix} \sigma_{2l',2l}^{n+1} & \sigma_{2l',2l+1}^{n+1} \\ \sigma_{2l'+1,2l}^{n+1} & \sigma_{2l'+1,2l+1}^{n+1} \end{pmatrix} U^T. \qquad (2.16) $$

Finally, we can describe a practical algorithm for the operator application in the nonstandard form. Let us for a given arbitrary $f \in L^2(0, 1)$ consider the sequences $\tilde{d}_l^n, \tilde{s}_l^n, \tilde{s}_l^n$ with $n \in \mathbb{N}_0, l = 0, 1, \dots, 2^n - 1$ of elements from $\mathbb{C}^\mathcal{t}$ defined by the following iterative procedure

$$ \tilde{d}_l^n = \sum_{l=0}^{2^n-1} (\alpha_{l'l}^n d_l^n(f) + \beta_{l'l}^n s_l^n(f)), \quad \tilde{s}_l^n = \sum_{l=0}^{2^n-1} \gamma_{l'l}^n d_l^n(f) $$

and

$$ \begin{align*} \tilde{s}_{2m}^n &= H^{(0)} (\tilde{s}_m^{n-1} + \tilde{s}_m^{n-1}) + G^{(0)} \tilde{d}_m^{n-1}, & n \ge 1, \quad \text{where } \tilde{s}_0^0 = \sigma_{00}^0 s_0^0(f). \\ \tilde{s}_{2m+1}^n &= H^{(1)} (\tilde{s}_m^{n-1} + \tilde{s}_m^{n-1}) + G^{(1)} \tilde{d}_m^{n-1}, & n \ge 1. \end{align*} $$

**Lemma 1.** Let $f \in L^2(0, 1)$, $n \in \mathbb{N}_0$ and $g = \mathcal{T}_{n+1}f$. Then its coefficients can be obtained by the above iterative procedure as

$$ d_l^n(g) = \tilde{d}_l^n, \quad s_l^n(g) = \tilde{s}_l^n + \tilde{s}_l^n \qquad (2.17) $$

at the finest scale $n$. The rest coefficients, at scales $n = 0, \dots, n-1$, are found by (2.11).

*Proof.* It is enough to prove (2.17), since the scale relation (2.11) is already known. Moreover, the first equality in (2.17) is obvious. Now we can notice that the second equality in (2.17) clearly
holds for the finest scale $n = 0$. Thus the statement follows by induction over $n$ from the equalities

$$ \int_0^1 \phi_{pl}^n(x) \phi_{j,2m}^{n+1}(x) dx = h_{pj}^{(0)} \delta_{lm}, \quad \int_0^1 \psi_{pl}^n(x) \phi_{j,2m}^{n+1}(x) dx = g_{pj}^{(0)} \delta_{lm} \\ \int_0^1 \phi_{pl}^n(x) \phi_{j,2m+1}^{n+1}(x) dx = h_{pj}^{(1)} \delta_{lm}, \quad \int_0^1 \psi_{pl}^n(x) \phi_{j,2m+1}^{n+1}(x) dx = g_{pj}^{(1)} \delta_{lm} $$

that are straightforward to check. Further details are omitted. □

### 3. REVIEW ON TIME EVOLUTION OPERATORS IN MULTIWAVELET BASES

There are two main examples of the operator $\mathcal{T}$ that we regard here, namely, the heat exp $(t\partial_x^2)$ and Schrödinger exp $(it\partial_x^2)$ evolution operators. The former one is probably the most used and studied operator in relation to application of multiwavelets. In this section we remind how exp $(t\partial_x^2)$ is descretised in the multiwavelet basis, and what problems arise with attempts of direct adaptation of the same techniques to the Schrödinger semigroup exp $(it\partial_x^2)$.

#### 3.1. Finite interval representation with Dirichlet boundary conditions.
Laplacian $\partial_x^2$ together with the Dirichlet boundary conditions, or more precisely, defined on the domain

$$ \mathrm{Dom} (\partial_x^2) = \{ u \in H_2^2(0, 1) \mid u(0) = u(1) = 0 \}, $$

constitute an unbounded self-adjoint operator in $L^2(0, 1)$.

Collection of functions $\{x \mapsto \sqrt{2} \sin(\pi\nu x)\}_{\nu \in \mathbb{N}}$ form an orthonormal basis in $L^2(0, 1)$. Let $\mathbb{F}$ be the Fourier transform associated with this basis, that is

$$ \mathbb{F}u(\nu) = \int_0^1 u(x)\sqrt{2}\sin(\pi\nu x)dx. $$

It is obvious that elements of this basis are eigenfunctions of the Laplacian $\partial_x^2$ corresponding to eigenvalues $-(\pi\nu)^2$, $\nu = 1, 2, \dots$. Therefore according to the general spectral theory for any function $\Psi$ we can define the following operator

$$ \Psi(\partial_x^2) u(x) = \sum_{\nu=1}^{\infty} \Psi(-(\pi\nu)^2) \int_0^1 u(x')\sqrt{2}\sin(\pi\nu x')dx'\sqrt{2}\sin(\pi\nu x) $$

on functions $u \in L^2(0, 1)$. The same can be written down in terms of the Fourier transform

$$ \Psi(\partial_x^2) = \mathbb{F}^{-1}\Psi(-(\pi\nu)^2)\mathbb{F}, $$

where we have a composition of Fourier operators and multiplication by symbol $\Psi(-(\pi\nu)^2)$. Note that

$$ \mathbb{F}^{-1} : \ell^2(\mathbb{N}) \ni \alpha \mapsto \sum_{\nu=1}^{\infty} \alpha_{\nu} \sqrt{2} \sin(\pi \nu x) \in L^2(0, 1). $$

The main example we are interested in is $\Psi(\lambda) = \exp(it\lambda)$. The heat exponent with $\Psi(\lambda) = \exp(t\lambda)$ was considered in [7]. The propagator $\exp(it\partial_x^2)$ is defined via symbol $\exp(-it(\pi\nu)^2)$ as described above. One would like to represent this operator in the multiwavelet basis in the same spirit as in [7]. Following the same logic we start with

$$ P^n u(x) = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n \phi_{jl}^n(x), \quad N = 2^n, $$

and look for its evolution in time defined by the operator $\mathcal{T} = \Psi(\partial_x^2)$. Function $\mathcal{T}u$ is approximated by $\mathcal{T}_n u$ which lies in the span of piecewise Legendre polynomials. We have

$$ \mathcal{T}_n u = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} \tilde{s}_{jl}^n \phi_{jl}^n, $$
holds for the finest scale $n = 0$. Thus the statement follows by induction over $n$ from the equalities

$$ 
\begin{aligned}
\int_0^1 \phi_{pl}^n(x) \phi_{j,2m}^{n+1}(x) dx &= h_{pj}^{(0)} \delta_{lm}, & \int_0^1 \psi_{pl}^n(x) \phi_{j,2m}^{n+1}(x) dx &= g_{pj}^{(0)} \delta_{lm} \\
\int_0^1 \phi_{pl}^n(x) \phi_{j,2m+1}^{n+1}(x) dx &= h_{pj}^{(1)} \delta_{lm}, & \int_0^1 \psi_{pl}^n(x) \phi_{j,2m+1}^{n+1}(x) dx &= g_{pj}^{(1)} \delta_{lm}
\end{aligned}
$$

that are straightforward to check. Further details are omitted. □

### 3. REVIEW ON TIME EVOLUTION OPERATORS IN MULTIWAVELET BASES

There are two main examples of the operator $\mathcal{T}$ that we regard here, namely, the heat exp $(t\partial_x^2)$ and Schrödinger exp $(it\partial_x^2)$ evolution operators. The former one is probably the most used and studied operator in relation to application of multiwavelets. In this section we remind how exp $(t\partial_x^2)$ is descretised in the multiwavelet basis, and what problems arise with attempts of direct adaptation of the same techniques to the Schrödinger semigroup exp $(it\partial_x^2)$.

#### 3.1. Finite interval representation with Dirichlet boundary conditions.
Laplacian $\partial_x^2$ together with the Dirichlet boundary conditions, or more precisely, defined on the domain

$$ \mathrm{Dom} (\partial_x^2) = \{ u \in H_2^2(0, 1) \mid u(0) = u(1) = 0 \}, $$

constitute an unbounded self-adjoint operator in $L^2(0, 1)$.

Collection of functions $\{x \mapsto \sqrt{2} \sin(\pi\nu x)\}_{\nu \in \mathbb{N}}$ form an orthonormal basis in $L^2(0, 1)$. Let $\mathbb{F}$ be the Fourier transform associated with this basis, that is

$$ \mathbb{F}u(\nu) = \int_0^1 u(x)\sqrt{2}\sin(\pi\nu x)dx. $$

It is obvious that elements of this basis are eigenfunctions of the Laplacian $\partial_x^2$ corresponding to eigenvalues $-(\pi\nu)^2$, $\nu = 1, 2, \dots$. Therefore according to the general spectral theory for any function $\Psi$ we can define the following operator

$$ \Psi(\partial_x^2) u(x) = \sum_{\nu=1}^{\infty} \Psi(-(\pi\nu)^2) \int_0^1 u(x')\sqrt{2}\sin(\pi\nu x')dx'\sqrt{2}\sin(\pi\nu x) $$

on functions $u \in L^2(0, 1)$. The same can be written down in terms of the Fourier transform

$$ \Psi(\partial_x^2) = \mathbb{F}^{-1}\Psi(-(\pi\nu)^2)\mathbb{F}, $$

where we have a composition of Fourier operators and multiplication by symbol $\Psi(-(\pi\nu)^2)$. Note that

$$ \mathbb{F}^{-1} : \ell^2(\mathbb{N}) \ni \alpha \mapsto \sum_{\nu=1}^{\infty} \alpha_{\nu} \sqrt{2} \sin(\pi \nu x) \in L^2(0, 1). $$

The main example we are interested in is $\Psi(\lambda) = \exp(it\lambda)$. The heat exponent with $\Psi(\lambda) = \exp(t\lambda)$ was considered in [7]. The propagator $\exp(it\partial_x^2)$ is defined via symbol $\exp(-it(\pi\nu)^2)$ as described above. One would like to represent this operator in the multiwavelet basis in the same spirit as in [7]. Following the same logic we start with

$$ P^n u(x) = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n \phi_{jl}^n(x), \quad N = 2^n, $$

and look for its evolution in time defined by the operator $\mathcal{T} = \Psi(\partial_x^2)$. Function $\mathcal{T}u$ is approximated by $\mathcal{T}_n u$ which lies in the span of piecewise Legendre polynomials. We have

$$ \mathcal{T}_n u = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} \tilde{s}_{jl}^n \phi_{jl}^n, $$
where

$$ \tilde{s}_{j'l'}^n = \int_0^1 \Psi(\partial_x^2) P^n u(x) \phi_{j'l'}^n(x) dx = \int_0^1 \mathbb{F}^{-1} \left( \Psi(-\pi\nu)^2 \mathbb{F} P^n u(\nu) \right) (x) \phi_{j'l'}^n(x) dx. $$

The latter is the inner product in $L^2(0, 1)$. It can be rewritten in terms of the inner product in $\ell^2(\mathbb{N})$ by noticing that $\mathbb{F}^* = \mathbb{F}^{-1}$ which follows immediately from the Parseval's identity. Thus

$$ \tilde{s}_{j'l'}^n = \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \mathbb{F} P^n u(\nu) \mathbb{F} \phi_{j'l'}^n(\nu) = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \mathbb{F} \phi_{jl}^n(\nu) \mathbb{F} \phi_{j'l'}^n(\nu). $$

Here we make a first remark on the convergence of this series. Our main example $\Psi(-\pi\nu)^2 = \exp(-it(\pi\nu)^2)$ is a function bounded with respect to $\nu \in \mathbb{N}$. Both $\mathbb{F}\phi_{jl}^n$ and $\mathbb{F}\phi_{j'l'}^n$ are $\ell^2$-sequences. Therefore the series is absolutely convergent by the Cauchy-Schwarz's inequality. So far we do not have much information about the rate of convergence, though one can anticipate a significant cancellations at high frequencies $\nu$ due to oscillations caused by the exponent.

Now let us simplify further the expression for $\mathbb{F}\phi_{jl}^n(\nu)$ in the following way

$$ \begin{aligned} \mathbb{F}\phi_{jl}^n(\nu) &= \int_0^1 \phi_{jl}^n(x) \sqrt{2} \sin(\pi\nu x) dx = \frac{1}{\sqrt{N}} \int_0^1 \phi_j(x) \sqrt{2} \sin\left(\pi\nu \frac{x+l}{N}\right) dx \\ &= -\frac{i}{\sqrt{2N}} \int_0^1 \phi_j(x) \left( \exp\left(i\pi\nu \frac{x+l}{N}\right) - \exp\left(-i\pi\nu \frac{x+l}{N}\right) \right) dx. \end{aligned} $$

If we introduce the usual Fourier transform as

$$ \hat{f}(\xi) = \mathcal{F}f(\xi) = \int_{\mathbb{R}} f(x)e^{-i\xi x} dx \qquad (3.1) $$

then after simple algebraic manipulations we get

$$ \begin{aligned} \mathbb{F}\phi_{jl}^n(\nu)\mathbb{F}\phi_{j'l'}^n(\nu) = \frac{1}{N} \operatorname{Re} \left( \mathbb{F}\phi_j\left(\frac{\pi\nu}{N}\right) \mathbb{F}\phi_{j'}\left(-\frac{\pi\nu}{N}\right) \exp\left(\frac{i\pi\nu(l'-l)}{N}\right) \right. \\ \left. - \mathbb{F}\phi_j\left(-\frac{\pi\nu}{N}\right) \mathbb{F}\phi_{j'}\left(-\frac{\pi\nu}{N}\right) \exp\left(\frac{i\pi\nu(l'+l)}{N}\right) \right). \end{aligned} $$

Summing up we obtain

$$ \tilde{s}_{j'l'}^n = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n \left( [\sigma_{l'-l}^{1n}]_{j'j} - [\sigma_{l'+l}^{2n}]_{j'j} \right), \quad N = 2^n, $$

and so $\mathcal{T}_n$ has matrix elements $\sigma_{l'l}^{1n} = \sigma_{l'-l}^{1n} - \sigma_{l'+l}^{2n}$ in $\{\phi_{jl}^n\}$, where we introduced the following notations

$$ [\sigma_l^{1n}]_{j'j} = \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \frac{1}{N} \operatorname{Re} \mathbb{F}\phi_j\left(\frac{\pi\nu}{N}\right) \overline{\mathbb{F}\phi_{j'}\left(\frac{\pi\nu}{N}\right)} \exp\left(\frac{i\pi\nu l}{N}\right), $$

$$ [\sigma_l^{2n}]_{j'j} = \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \frac{1}{N} \operatorname{Re} \overline{\mathbb{F}\phi_j\left(\frac{\pi\nu}{N}\right) \mathbb{F}\phi_{j'}\left(\frac{\pi\nu}{N}\right)} \exp\left(\frac{i\pi\nu l}{N}\right). $$

A big advantage of these formulas is that the Fourier transforms can be calculated exactly, since we integrate polynomials multiplied by exponents. Moreover, a recursion formula for the Legendre polynomials implies a recursion for the corresponding Fourier transforms. Indeed, integrating the following relation [4] with an exponential weight

$$ (2j+1)P_j(x) = \partial_x(P_{j+1}(x) - P_{j-1}(x)) $$

one obtains

$$ \frac{2j+1}{i\xi} \int_{-1}^{1} P_j(x) e^{-i\xi x} dx = \int_{-1}^{1} P_{j+1}(x) e^{-i\xi x} dx - \int_{-1}^{1} P_{j-1}(x) e^{-i\xi x} dx. $$
where

$$ \tilde{s}_{j'l'}^n = \int_0^1 \Psi(\partial_x^2) P^n u(x) \phi_{j'l'}^n(x) dx = \int_0^1 \mathbb{F}^{-1} \left( \Psi(-\pi\nu)^2 \mathbb{F} P^n u(\nu) \right) (x) \phi_{j'l'}^n(x) dx. $$

The latter is the inner product in $L^2(0, 1)$. It can be rewritten in terms of the inner product in $\ell^2(\mathbb{N})$ by noticing that $\mathbb{F}^* = \mathbb{F}^{-1}$ which follows immediately from the Parseval's identity. Thus

$$ \tilde{s}_{j'l'}^n = \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \mathbb{F} P^n u(\nu) \mathbb{F} \phi_{j'l'}^n(\nu) = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \mathbb{F} \phi_{jl}^n(\nu) \mathbb{F} \phi_{j'l'}^n(\nu). $$

Here we make a first remark on the convergence of this series. Our main example $\Psi(-\pi\nu)^2 = \exp(-it(\pi\nu)^2)$ is a function bounded with respect to $\nu \in \mathbb{N}$. Both $\mathbb{F}\phi_{jl}^n$ and $\mathbb{F}\phi_{j'l'}^n$ are $\ell^2$-sequences. Therefore the series is absolutely convergent by the Cauchy-Schwarz's inequality. So far we do not have much information about the rate of convergence, though one can anticipate a significant cancellations at high frequencies $\nu$ due to oscillations caused by the exponent.

Now let us simplify further the expression for $\mathbb{F}\phi_{jl}^n(\nu)$ in the following way

$$ \begin{aligned} \mathbb{F}\phi_{jl}^n(\nu) &= \int_0^1 \phi_{jl}^n(x) \sqrt{2} \sin(\pi\nu x) dx = \frac{1}{\sqrt{N}} \int_0^1 \phi_j(x) \sqrt{2} \sin\left(\pi\nu \frac{x+l}{N}\right) dx \\ &= -\frac{i}{\sqrt{2N}} \int_0^1 \phi_j(x) \left( \exp\left(i\pi\nu \frac{x+l}{N}\right) - \exp\left(-i\pi\nu \frac{x+l}{N}\right) \right) dx. \end{aligned} $$

If we introduce the usual Fourier transform as

$$ \hat{f}(\xi) = \mathcal{F}f(\xi) = \int_{\mathbb{R}} f(x)e^{-i\xi x} dx \quad (3.1) $$

then after simple algebraic manipulations we get

$$ \begin{aligned} \mathbb{F}\phi_{jl}^n(\nu)\mathbb{F}\phi_{j'l'}^n(\nu) = \frac{1}{N} \operatorname{Re} \left( \mathbb{F}\phi_j\left(\frac{\pi\nu}{N}\right) \mathbb{F}\phi_{j'}\left(-\frac{\pi\nu}{N}\right) \exp\left(\frac{i\pi\nu(l'-l)}{N}\right) \right. \\ \left. - \mathbb{F}\phi_j\left(-\frac{\pi\nu}{N}\right) \mathbb{F}\phi_{j'}\left(-\frac{\pi\nu}{N}\right) \exp\left(\frac{i\pi\nu(l'+l)}{N}\right) \right). \end{aligned} $$

Summing up we obtain

$$ \tilde{s}_{j'l'}^n = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n \left( [\sigma_{l'-l}^{1n}]_{j'j} - [\sigma_{l'+l}^{2n}]_{j'j} \right), \quad N = 2^n, $$

and so $\mathcal{T}_n$ has matrix elements $\sigma_{l'l}^{1n} = \sigma_{l'-l}^{1n} - \sigma_{l'+l}^{2n}$ in $\{\phi_{jl}^n\}$, where we introduced the following notations

$$ [\sigma_l^{1n}]_{j'j} = \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \frac{1}{N} \operatorname{Re} \mathbb{F}\phi_j\left(\frac{\pi\nu}{N}\right) \overline{\mathbb{F}\phi_{j'}\left(\frac{\pi\nu}{N}\right)} \exp\left(\frac{i\pi\nu l}{N}\right), $$

$$ [\sigma_l^{2n}]_{j'j} = \sum_{\nu=1}^{\infty} \Psi(-\pi\nu)^2 \frac{1}{N} \operatorname{Re} \overline{\mathbb{F}\phi_j\left(\frac{\pi\nu}{N}\right) \mathbb{F}\phi_{j'}\left(\frac{\pi\nu}{N}\right)} \exp\left(\frac{i\pi\nu l}{N}\right). $$

A big advantage of these formulas is that the Fourier transforms can be calculated exactly, since we integrate polynomials multiplied by exponents. Moreover, a recursion formula for the Legendre polynomials implies a recursion for the corresponding Fourier transforms. Indeed, integrating the following relation [4] with an exponential weight

$$ (2j + 1)P_j(x) = \partial_x(P_{j+1}(x) - P_{j-1}(x)) $$

one obtains

$$ \frac{2j+1}{i\xi} \int_{-1}^{1} P_j(x) e^{-i\xi x} dx = \int_{-1}^{1} P_{j+1}(x) e^{-i\xi x} dx - \int_{-1}^{1} P_{j-1}(x) e^{-i\xi x} dx. $$
Hence the Legendre scaling functions (2.2) enjoy

$$ \mathcal{F}\phi_{j+1}(\xi) = \sqrt{\frac{2j+3}{2j-1}} \mathcal{F}\phi_{j-1}(\xi) + \frac{2\sqrt{(2j+1)(2j+3)}}{i\xi} \mathcal{F}\phi_j(\xi), \quad (3.2) $$

whereas the first two Fourier transforms

$$ \mathcal{F}\phi_0(\xi) = \frac{1}{i\xi} (1 - e^{-i\xi}), \quad (3.3) $$

$$ \mathcal{F}\phi_1(\xi) = \frac{i\sqrt{3}}{\xi} (1 + e^{-i\xi}) + \frac{2\sqrt{3}}{\xi^2} (e^{-i\xi} - 1) \quad (3.4) $$

can be obtained by direct calculations.

We remark again here that each Fourier transform $\mathcal{F}\phi_j(\frac{\pi\nu}{N}) = \mathcal{O}(\frac{N}{\nu})$, when $\nu \to \infty$, as follows from the calculation of Fourier transforms. This is prohibitively slow for our main example of bounded $\Psi$. Certainly, keeping only the first $N^m$ terms of the sums we get

$$ [\sigma_l^{1n}]_{j'j} = \sum_{\nu=1}^{N^m} \Psi(-(\pi\nu)^2) \frac{1}{N} \operatorname{Re} \mathcal{F}\phi_j\left(\frac{\pi\nu}{N}\right) \overline{\mathcal{F}\phi_{j'}\left(\frac{\pi\nu}{N}\right)} \exp\left(\frac{i\pi\nu l}{N}\right) + \mathcal{O}\left(\frac{1}{N^{m-1}}\right), $$

and similarly for $[\sigma_l^{2n}]_{j'j}$. Indeed, the rest is bounded up to a constant by

$$ \frac{1}{N} \sum_{\nu=N^m+1}^{\infty} \frac{1}{(\nu/N)^2} \le \int_{N^{m-1}}^{\infty} \frac{dx}{x^2} = \frac{1}{N^{m-1}}. $$

This is problematic since in practice an acceptable precision demands setting $m = 3$ and probably even more, which corresponds to the tolerance $\mathcal{O}(\frac{1}{N^2}) = \mathcal{O}(2^{-2n})$. This requires at least $N^3$ multiplications. In fact, the rest can be estimated more thoroughly and we could get a better estimate of it, making use of high oscillation of $\Psi$ mentioned above. Advantage of the smoothing property of the exponential operator will be exploited in Section 4.

For the heat equation the series $[\sigma_l^{1n}]_{j'j}$ and $[\sigma_l^{2n}]_{j'j}$ converge significantly faster thanks to the decay $\Psi(-(\pi\nu)^2) = \exp(-t(\pi\nu)^2)$. Moreover, the corresponding non-standard form blocks $\alpha^n, \beta^n$ and $\gamma^n$ defined in (2.15) turn out to be effectively sparse. We remark here that entries of matrices $[\sigma_l^{1n}]_{j'j}, [\sigma_l^{2n}]_{j'j}$ are calculated only once for fixed given $t$, or several times for high order temporal numerical schemes. Then each of these matrices applied to a vector $\{s_{jl}^n\}$ following a fast numerical procedure developed in [8]. More precisely, the sparse approximations of $\alpha^n, \beta^n$ and $\gamma^n$ are applied to $\{s_{jl}^n\}, \{d_{jl}^n\}$.

### 3.2. Convolution representation on the real line.
On the real line the operator $\mathcal{T} = \Psi(\partial_x^2)$ is a convolution of the form

$$ \mathcal{T}u(x) = \int_{\mathbb{R}} K(x-y)u(y)dy \quad (3.5) $$

with the kernel $K = \mathcal{F}_{\xi}^{-1}\Psi(-\xi^2)$. Its matrix elements (2.15) are the integrals

$$ [\sigma_{l'l}^n]_{j'j} = \int_0^1 \int_0^1 K(x-y)\phi_{j'l}^n(x)\phi_{jl}^n(y)dx dy $$

simplifying to

$$ [\sigma_{l'l}^n]_{j'j} = [\sigma_{l'-l}^n]_{j'j} = \frac{1}{N} \int_{-1}^{1} K\left(\frac{z+l'-l}{N}\right) \Phi_{j'j}(z) dz $$

according to (2.1), where $N = 2^n$ and the so called correlation functions

$$ \Phi_{j'j}(z) = \int_{0}^{1} \phi_{j'}(x) \phi_{j}(x-z) dx $$
Hence the Legendre scaling functions (2.2) enjoy

$$ \mathcal{F}\phi_{j+1}(\xi) = \sqrt{\frac{2j+3}{2j-1}} \mathcal{F}\phi_{j-1}(\xi) + \frac{2\sqrt{(2j+1)(2j+3)}}{i\xi} \mathcal{F}\phi_j(\xi), \qquad (3.2) $$

whereas the first two Fourier transforms

$$ \mathcal{F}\phi_0(\xi) = \frac{1}{i\xi} (1 - e^{-i\xi}), \qquad (3.3) $$

$$ \mathcal{F}\phi_1(\xi) = \frac{i\sqrt{3}}{\xi} (1 + e^{-i\xi}) + \frac{2\sqrt{3}}{\xi^2} (e^{-i\xi} - 1) \qquad (3.4) $$

can be obtained by direct calculations.

We remark again here that each Fourier transform $\mathcal{F}\phi_j(\frac{\pi\nu}{N}) = \mathcal{O}(\frac{N}{\nu})$, when $\nu \to \infty$, as follows from the calculation of Fourier transforms. This is prohibitively slow for our main example of bounded $\Psi$. Certainly, keeping only the first $N^m$ terms of the sums we get

$$ [\sigma_l^{1n}]_{j'j} = \sum_{\nu=1}^{N^m} \Psi(-(\pi\nu)^2) \frac{1}{N} \operatorname{Re} \mathcal{F}\phi_j\left(\frac{\pi\nu}{N}\right) \overline{\mathcal{F}\phi_{j'}\left(\frac{\pi\nu}{N}\right)} \exp\left(\frac{i\pi\nu l}{N}\right) + \mathcal{O}\left(\frac{1}{N^{m-1}}\right), $$

and similarly for $[\sigma_l^{2n}]_{j'j}$. Indeed, the rest is bounded up to a constant by

$$ \frac{1}{N} \sum_{\nu=N^m+1}^{\infty} \frac{1}{(\nu/N)^2} \le \int_{N^{m-1}}^{\infty} \frac{dx}{x^2} = \frac{1}{N^{m-1}}. $$

This is problematic since in practice an acceptable precision demands setting $m = 3$ and probably even more, which corresponds to the tolerance $\mathcal{O}(\frac{1}{N^2}) = \mathcal{O}(2^{-2n})$. This requires at least $N^3$ multiplications. In fact, the rest can be estimated more thoroughly and we could get a better estimate of it, making use of high oscillation of $\Psi$ mentioned above. Advantage of the smoothing property of the exponential operator will be exploited in Section 4.

For the heat equation the series $[\sigma_l^{1n}]_{j'j}$ and $[\sigma_l^{2n}]_{j'j}$ converge significantly faster thanks to the decay $\Psi(-(\pi\nu)^2) = \exp(-t(\pi\nu)^2)$. Moreover, the corresponding non-standard form blocks $\alpha^n, \beta^n$ and $\gamma^n$ defined in (2.15) turn out to be effectively sparse. We remark here that entries of matrices $[\sigma_l^{1n}]_{j'j}, [\sigma_l^{2n}]_{j'j}$ are calculated only once for fixed given $t$, or several times for high order temporal numerical schemes. Then each of these matrices applied to a vector $\{s_{jl}^n\}$ following a fast numerical procedure developed in [8]. More precisely, the sparse approximations of $\alpha^n, \beta^n$ and $\gamma^n$ are applied to $\{s_{jl}^n\}, \{d_{jl}^n\}$.

### 3.2. Convolution representation on the real line.
On the real line the operator $\mathcal{T} = \Psi(\partial_x^2)$ is a convolution of the form

$$ \mathcal{T}u(x) = \int_{\mathbb{R}} K(x-y)u(y)dy \qquad (3.5) $$

with the kernel $K = \mathcal{F}_{\xi}^{-1}\Psi(-\xi^2)$. Its matrix elements (2.15) are the integrals

$$ [\sigma_{l'l}^n]_{j'j} = \int_0^1 \int_0^1 K(x-y)\phi_{j'l}^n(x)\phi_{jl}^n(y)dxdy $$

simplifying to

$$ [\sigma_{l'l}^n]_{j'j} = [\sigma_{l'-l}^n]_{j'j} = \frac{1}{N} \int_{-1}^{1} K\left(\frac{z+l'-l}{N}\right) \Phi_{j'j}(z) dz $$

according to (2.1), where $N = 2^n$ and the so called correlation functions

$$ \Phi_{j'j}(z) = \int_{0}^{1} \phi_{j'}(x)\phi_{j}(x-z)dx $$
have been introduced. Each of them is continuous with the support in [-1, 1]. Their restrictions to either [-1, 0] or [0, 1] are polynomials of order at most 2k - 1. Therefore, one may expand the correlation functions in $L^2(-1, 1)$ as

$$ \Phi_{j'j}(z) = \sum_{p=0}^{2k-1} \left( c_{j'jp}^{(+)} \phi_p(z) + c_{j'jp}^{(-)} \phi_p(z+1) \right), $$

where the cross correlation coefficients

$$ \begin{cases} c_{j'jp}^{(+)} = \int_0^1 \int_0^1 \phi_{j'}(x) \phi_j(x-z) \phi_p(z) dx dz \\ c_{j'jp}^{(-)} = \int_0^1 \int_0^1 \phi_{j'}(x) \phi_j(x-z+1) \phi_p(z) dx dz \end{cases}, \quad j', j = 0, \dots, k-1, \quad p = 0, \dots, 2k-1, $$

are easily tabulated. Thus the convolution operator $\mathcal{T}$ in the orthonormal collection $\{\phi_{jl}^n\}$ has the matrix elements $[\sigma_{l'l}^n]_{j'j} = [\sigma_{l'l-1}^n]_{j'j}$ depending only on the distance $l' - l = 1 - N, \dots, N-1$ to the diagonal. These elements may be evaluated via

$$ [\sigma_l^n]_{j'j} = \frac{1}{\sqrt{N}} \sum_{p=0}^{2k-1} \left( c_{j'jp}^{(+)} s_{p,l}^n(K) + c_{j'jp}^{(-)} s_{p,l-1}^n(K) \right), \quad (3.6) $$

where $N = 2^n$ and the scaling coefficients of the kernel $K$ are defined by (2.4).

Numerical evaluation of convolution type operators is reviewed in [21], for instance. Formula (3.6) gives us only the first integral in (2.15). The rest can be calculated by the decomposition transformation (2.16) that simplifies to

$$ \begin{pmatrix} \sigma_l^n & \gamma_l^n \\ \beta_l^n & \alpha_l^n \end{pmatrix} = U \begin{pmatrix} \sigma_{2l}^{n+1} & \sigma_{2l-1}^{n+1} \\ \sigma_{2l+1}^{n+1} & \sigma_{2l}^{n+1} \end{pmatrix} U^T \quad (3.7) $$

for the convolution operator $\mathcal{T}$ and $l = -N+1, -N+2, \dots, N-1$ standing for the distance to the diagonal. Finally, we are ready to demonstrate how the multiwavelet machinery can be used for the exponential heat operator $\exp(t\partial_x^2)$. It can be regarded as a convolution operator in $L^2(\mathbb{R})$ of the form

$$ \exp(t\partial_x^2) u(x) = \frac{1}{\sqrt{4\pi t}} \int_{\mathbb{R}} \exp\left(-\frac{(x-y)^2}{4t}\right) u(y)dy, \quad t > 0, $$

see [20]. It is associated with the heat equation $\partial_t u = \partial_x^2 u$ on the real line $\mathbb{R}$. The Green's function scaling coefficients $s_{p,l}^n(K)$ can be calculated easily with high precision. Then from (3.6), (3.7) one obtains the nonstandard form matrices $\alpha^n, \beta^n, \gamma^n$ with $n = 0, 1, \dots$. They turn out to be narrow diagonal banded. As a matter of fact, this was used a lot in Chemistry. Indeed, the Green's function for the Helmholtz equation in the three dimension admits an expansion into sum of Gaussian functions [23, 26] similar to the one appearing in this heat semigroup convolution.

An attempt of adaptation of the above approach to the time-dependent Schrödinger equation (1.1) was made in [38]. The exponential operator $\exp(it\partial_x^2)$ can be regarded as a convolution operator in $L^2(\mathbb{R})$ of the form

$$ \exp (it\partial_x^2) u(x) = \frac{\exp(-i\pi/4)}{\sqrt{4\pi t}} \int_{\mathbb{R}} \exp \left( \frac{i(x-y)^2}{4t} \right) u(y)dy, \quad t > 0, $$

see [20]. Oscillatory behaviour of the kernel does not allow us to calculate the convolution with a given high precision. Therefore, in [38] this difficulty is overcome by damping high frequencies, namely, the kernel is approximated by

$$ K(x) \approx \mathcal{F}_{\xi}^{-1} \left( e^{-it\xi^2} \chi(\xi) \right) (x), $$

where $\chi(\xi)$ is a smooth cut off function. The problem with such a brutal approach is that it is difficult to control the error and the cut off $\chi(\xi)$ should be tuned for each particular problem.
have been introduced. Each of them is continuous with the support in [-1, 1]. Their restrictions to either [-1, 0] or [0, 1] are polynomials of order at most 2k - 1. Therefore, one may expand the correlation functions in $L^2(-1, 1)$ as

$$ \Phi_{j'j}(z) = \sum_{p=0}^{2k-1} \left( c_{j'jp}^{(+)} \phi_p(z) + c_{j'jp}^{(-)} \phi_p(z+1) \right), $$

where the cross correlation coefficients

$$ \begin{cases} c_{j'jp}^{(+)} = \int_0^1 \int_0^1 \phi_{j'}(x) \phi_j(x-z) \phi_p(z) dx dz \\ c_{j'jp}^{(-)} = \int_0^1 \int_0^1 \phi_{j'}(x) \phi_j(x-z+1) \phi_p(z) dx dz \end{cases}, \quad j', j = 0, \dots, k-1, \quad p = 0, \dots, 2k-1, $$

are easily tabulated. Thus the convolution operator $\mathcal{T}$ in the orthonormal collection $\{\phi_{jl}^n\}$ has the matrix elements $[\sigma_{l'l}^n]_{j'j} = [\sigma_{l'l-1}^n]_{j'j}$ depending only on the distance $l' - l = 1 - N, \dots, N-1$ to the diagonal. These elements may be evaluated via

$$ [\sigma_l^n]_{j'j} = \frac{1}{\sqrt{N}} \sum_{p=0}^{2k-1} \left( c_{j'jp}^{(+)} s_{p,l}^n(K) + c_{j'jp}^{(-)} s_{p,l-1}^n(K) \right), \quad (3.6) $$

where $N = 2^n$ and the scaling coefficients of the kernel $K$ are defined by (2.4).

Numerical evaluation of convolution type operators is reviewed in [21], for instance. Formula (3.6) gives us only the first integral in (2.15). The rest can be calculated by the decomposition transformation (2.16) that simplifies to

$$ \begin{pmatrix} \sigma_l^n & \gamma_l^n \\ \beta_l^n & \alpha_l^n \end{pmatrix} = U \begin{pmatrix} \sigma_{2l}^{n+1} & \sigma_{2l-1}^{n+1} \\ \sigma_{2l+1}^{n+1} & \sigma_{2l}^{n+1} \end{pmatrix} U^T \quad (3.7) $$

for the convolution operator $\mathcal{T}$ and $l = -N+1, -N+2, \dots, N-1$ standing for the distance to the diagonal. Finally, we are ready to demonstrate how the multiwavelet machinery can be used for the exponential heat operator $\exp(t\partial_x^2)$. It can be regarded as a convolution operator in $L^2(\mathbb{R})$ of the form

$$ \exp(t\partial_x^2) u(x) = \frac{1}{\sqrt{4\pi t}} \int_{\mathbb{R}} \exp\left(-\frac{(x-y)^2}{4t}\right) u(y)dy, \quad t > 0, $$

see [20]. It is associated with the heat equation $\partial_t u = \partial_x^2 u$ on the real line $\mathbb{R}$. The Green's function scaling coefficients $s_{p,l}^n(K)$ can be calculated easily with high precision. Then from (3.6), (3.7) one obtains the nonstandard form matrices $\alpha^n, \beta^n, \gamma^n$ with $n = 0, 1, \dots$. They turn out to be narrow diagonal banded. As a matter of fact, this was used a lot in Chemistry. Indeed, the Green's function for the Helmholtz equation in the three dimension admits an expansion into sum of Gaussian functions [23, 26] similar to the one appearing in this heat semigroup convolution.

An attempt of adaptation of the above approach to the time-dependent Schrödinger equation (1.1) was made in [38]. The exponential operator $\exp(it\partial_x^2)$ can be regarded as a convolution operator in $L^2(\mathbb{R})$ of the form

$$ \exp (it\partial_x^2) u(x) = \frac{\exp(-i\pi/4)}{\sqrt{4\pi t}} \int_{\mathbb{R}} \exp \left( \frac{i(x-y)^2}{4t} \right) u(y)dy, \quad t > 0, $$

see [20]. Oscillatory behaviour of the kernel does not allow us to calculate the convolution with a given high precision. Therefore, in [38] this difficulty is overcome by damping high frequencies, namely, the kernel is approximated by

$$ K(x) \approx \mathcal{F}_{\xi}^{-1} \left( e^{-it\xi^2} \chi(\xi) \right) (x), $$

where $\chi(\xi)$ is a smooth cut off function. The problem with such a brutal approach is that it is difficult to control the error and the cut off $\chi(\xi)$ should be tuned for each particular problem.
Moreover, as in the previous section this approach ignores the smoothing property of the exponential operator.

### 4. CONTOUR DEFORMATION TECHNIQUE

Numerically, a use of the smoothing property of the exponential operator $\mathcal{T} = \exp(it\partial_x^2)$ is made in [25], where they regarded Equation (1.1) in the frequency domain. Returning back to the physical domain, where one applies the potential $V(x, t)$, they calculate the inverse Fourier transform $\mathcal{F}^{-1}$ over a specifically deformed contour $\Gamma$ instead of $\mathbb{R}$, which allows them to increase the precision significantly, while keeping advantage of using fast Fourier transform based schemes. In this section we adopt their approach in order to get a precise multiresolution representation of $\mathcal{T}$. As in Section 3.1, we can write

$$ \mathcal{T}_n u = P^n \exp (it\partial_x^2) P^n u = P^n \mathcal{F}^{-1} e^{-it\xi^2} \mathcal{F} P^n u = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} \tilde{s}_{jl}^n \phi_{jl}^n, $$

where now the operator $\mathcal{T}$ is unitary in $L^2(\mathbb{R})$. Similarly, we obtain

$$ \tilde{s}_{j'l'}^n = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n [\sigma_{l'-l}^n]_{j'j}, \quad N = 2^n, $$

where the matrix $[\sigma_{l'l}^n]_{j'j}$ defined in (2.15) depends only on the distance $l' - l$ to the diagonal, as explained in Section 3.2. And so in terms of the introduced notations we have the following expression

$$ \begin{aligned} [\sigma_l^n]_{j'j} &= \frac{1}{2\pi N} \int_{\mathbb{R}} \exp\left(\frac{i\xi l}{N} - it\xi^2\right) \mathcal{F}\phi_j\left(\frac{\xi}{N}\right) \overline{\mathcal{F}\phi_{j'}\left(\frac{\xi}{N}\right)} d\xi \\ &= \frac{1}{2\pi} \int_{\mathbb{R}} \exp\left(i\xi l - itN^2\xi^2\right) \mathcal{F}\phi_j(\xi) \mathcal{F}\phi_{j'}(-\xi) d\xi. \end{aligned} $$

Note that $[\sigma_l^n]_{j'j} = [\sigma_{-l}^n]_{jj'}$. Fourier transforms of the scaling functions are given by (3.2), (3.3), (3.4). Clearly, the integrand can be extended to an entire function of complex variable $\zeta \in \mathbb{C}$. In two quadrants of $\mathbb{C}$-plane we have $\text{Re}(-itN^2\zeta^2) < 0$, which suggests deforming the integration contour into them for more effective calculation of the integral with respect to $\zeta$. Moreover, we could in principal deform the integration contour exactly as was done in [25], so that the new contour $\Gamma$ would be chosen depending on the sign of $t$ and an accuracy parameter $\varepsilon$. Without loss of generality one can assume $t > 0$ from now on. Their contour $\Gamma_H$ is kept in an $H$-neighbourhood of the real axis $\mathbb{R}$, i.e. in the band $|\text{Im}\,\zeta| \le H$, where the bound $H > 0$ is introduced in order to avoid multiplication of big and small numbers while calculating the integral. A direct repetition of their argument leads to the optimal bound

$$ H = \frac{1}{N} \log \frac{\pi \varepsilon}{2\varepsilon_{\text{mach}}}, \qquad (4.1) $$

where $\varepsilon$ is the desired precision and $\varepsilon_{\text{mach}}$ is the machine epsilon (the precision limit of floating-point arithmetic on a computer, approximately $2^{-52}$ for double type), see the details in Subsection 4.3. However, in our case we can use as an advantage the fact that we are able to calculate Fourier transforms of piecewise polynomials $\phi_j$ exactly. Indeed, all the integrands are sums of exponents (up to powers of $1/\zeta$), and so, intuitively, combining them together one can avoid inaccurate multiplications. Without this restriction we are able to take a contour that will allow us to exploit the smoothing semigroup property at its best. Namely, we choose a contour to be the line $\Gamma = (1 - i)\mathbb{R}$ oriented with angle $-\pi/4$ to the real line. On such contour $\Gamma$ the main exponent part $-itN^2\zeta^2 = -tN^2|\zeta|^2$, which guarantees fast integral convergence.
Moreover, as in the previous section this approach ignores the smoothing property of the exponential operator.

## 4. CONTOUR DEFORMATION TECHNIQUE

Numerically, a use of the smoothing property of the exponential operator $\mathcal{T} = \exp(it\partial_x^2)$ is made in [25], where they regarded Equation (1.1) in the frequency domain. Returning back to the physical domain, where one applies the potential $V(x, t)$, they calculate the inverse Fourier transform $\mathcal{F}^{-1}$ over a specifically deformed contour $\Gamma$ instead of $\mathbb{R}$, which allows them to increase the precision significantly, while keeping advantage of using fast Fourier transform based schemes. In this section we adopt their approach in order to get a precise multiresolution representation of $\mathcal{T}$. As in Section 3.1, we can write

$$ \mathcal{T}_n u = P^n \exp (it\partial_x^2) P^n u = P^n \mathcal{F}^{-1} e^{-it\xi^2} \mathcal{F} P^n u = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} \tilde{s}_{jl}^n \phi_{jl}^n, $$

where now the operator $\mathcal{T}$ is unitary in $L^2(\mathbb{R})$. Similarly, we obtain

$$ \tilde{s}_{j'l'}^n = \sum_{l=0}^{N-1} \sum_{j=0}^{\ell-1} s_{jl}^n [\sigma_{l'-l}^n]_{j'j}, \quad N = 2^n, $$

where the matrix $[\sigma_{l'l}^n]_{j'j}$ defined in (2.15) depends only on the distance $l' - l$ to the diagonal, as explained in Section 3.2. And so in terms of the introduced notations we have the following expression

$$ [\sigma_l^n]_{j'j} = \frac{1}{2\pi N} \int_{\mathbb{R}} \exp\left(\frac{i\xi l}{N} - it\xi^2\right) \mathcal{F}\phi_j\left(\frac{\xi}{N}\right) \overline{\mathcal{F}\phi_{j'}\left(\frac{\xi}{N}\right)} d\xi \\ = \frac{1}{2\pi} \int_{\mathbb{R}} \exp\left(i\xi l - itN^2\xi^2\right) \mathcal{F}\phi_j(\xi) \mathcal{F}\phi_{j'}(-\xi) d\xi. $$

Note that $[\sigma_l^n]_{j'j} = [\sigma_{-l}^n]_{jj'}$. Fourier transforms of the scaling functions are given by (3.2), (3.3), (3.4). Clearly, the integrand can be extended to an entire function of complex variable $\zeta \in \mathbb{C}$. In two quadrants of $\mathbb{C}$-plane we have $\text{Re}(-itN^2\zeta^2) < 0$, which suggests deforming the integration contour into them for more effective calculation of the integral with respect to $\zeta$. Moreover, we could in principal deform the integration contour exactly as was done in [25], so that the new contour $\Gamma$ would be chosen depending on the sign of $t$ and an accuracy parameter $\varepsilon$. Without loss of generality one can assume $t > 0$ from now on. Their contour $\Gamma_H$ is kept in an $H$-neighbourhood of the real axis $\mathbb{R}$, i.e. in the band $|\text{Im}\,\zeta| \le H$, where the bound $H > 0$ is introduced in order to avoid multiplication of big and small numbers while calculating the integral. A direct repetition of their argument leads to the optimal bound

$$ H = \frac{1}{N} \log \frac{\pi \varepsilon}{2\varepsilon_{\text{mach}}}, \quad (4.1) $$

where $\varepsilon$ is the desired precision and $\varepsilon_{\text{mach}}$ is the machine epsilon (the precision limit of floating-point arithmetic on a computer, approximately $2^{-52}$ for double type), see the details in Subsection 4.3. However, in our case we can use as an advantage the fact that we are able to calculate Fourier transforms of piecewise polynomials $\phi_j$ exactly. Indeed, all the integrands are sums of exponents (up to powers of $1/\zeta$), and so, intuitively, combining them together one can avoid inaccurate multiplications. Without this restriction we are able to take a contour that will allow us to exploit the smoothing semigroup property at its best. Namely, we choose a contour to be the line $\Gamma = (1 - i)\mathbb{R}$ oriented with angle $-\pi/4$ to the real line. On such contour $\Gamma$ the main exponent part $-itN^2\zeta^2 = -tN^2|\zeta|^2$, which guarantees fast integral convergence.
Thus the nonstandard form matrices (2.15) take the final forms

$$[\sigma_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\phi_j(\zeta)\mathcal{F}\phi_p(-\zeta)d\zeta, \quad (4.2)$$

$$[\alpha_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\psi_j(\zeta)\mathcal{F}\psi_p(-\zeta)d\zeta, \quad (4.3)$$

$$[\beta_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\phi_j(\zeta)\mathcal{F}\psi_p(-\zeta)d\zeta, \quad (4.4)$$

$$[\gamma_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\psi_j(\zeta)\mathcal{F}\phi_p(-\zeta)d\zeta, \quad (4.5)$$

where the matrices are effectively sparse, as we shall see below, due to $\mathcal{F}\psi_j(0) = 0$ of at least order $\mathfrak{k}$. We point out here that Fourier transform $\mathcal{F}\psi_j(\zeta)$ can be easily found analytically and extended to the complex plane $\zeta \in \mathbb{C}$ as

$$\hat{\psi}_m(\zeta) = \frac{1}{\sqrt{2}} \sum_{j=0}^{\mathfrak{k}-1} \hat{\phi}_j\left(\frac{\zeta}{2}\right) \left[ g_{ij}^{(0)} + g_{ij}^{(1)} \exp\left(-\frac{i\zeta}{2}\right) \right],$$

where $g_{ij}^{(0)}$ and $g_{ij}^{(1)}$ are elements of the filter matrices $G^{(0)}$ and $G^{(1)}$, respectively. Justification of the contour deformation is straightforward, so we omit the proof.

Before we continue with evaluation of these integrals it is worth to make a couple of remarks on behaviour of the integrands in (4.2)-(4.5). As one can see from (3.2), (3.3), (3.4) that the Legendre Fourier transform extensions $\hat{\phi}_j(\zeta)$ have a removable singularity at zero $\zeta = 0$. Clearly, the same is true for the interpolating Fourier transform extensions $\hat{\phi}_j(\zeta)$ and the wavelet Fourier transform extensions $\hat{\psi}_j(\zeta)$. This suggests us to approximate these functions by power series around zero. Moreover, taking into account the Gaussian factor $\exp(-tN^2|\zeta|^2)$, one may expect that these power series will do a good job even away from zero $\zeta = 0$. Below we are mostly focused on further development of this idea. The last remark concerns the price one pays for the contour deformation. As it will be obvious below, apart from the Gaussian fast decreasing factor we also get the increasing factor $\exp((|l|+1)|\zeta|/\sqrt{2})$, that turns out to be crucial, due to round up error, when the time step $t$ is small and the scale $n$ is low. This problem is overcome in practice by the fact that we need the first integral (4.2) to be calculated with high precision at high scales $n$, depending on $t$. The rest are calculated by (3.7).

### 4.1. Legendre scaling functions.

$$[\sigma_l^n]_{pj}(a) = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - ia\zeta^2) \hat{\phi}_j(\zeta) \hat{\phi}_p(-\zeta) d\zeta, \quad 0 \le p, j < \mathfrak{k}, \quad (4.6)$$

where $a = tN^2 = t4^n$. We can represent the multiplication of Fourier transforms as a $\zeta$-power series. This will transform this integral into a series of integrals $J_k(a, l)$, that can be evaluated exactly, with some coefficients depending solely on $p, j, k$. These coefficients need to be computed only once and tabulated, as they are problem-independent. We will call them correlation coefficients. We make an extensive use of (3.2) which is valid for integers $j \ge 1$.

Firstly, we notice that $\hat{\phi}_j(\zeta)$ has a root at zero of order $j$. Indeed,

$$\partial_{\zeta}^{m} \hat{\phi}_{j}(0) = (-i)^{m} \int_{\mathbb{R}} \phi_{j}(x) x^{m} dx$$

that equals zero for any $0 \le m < j$ and a non-zero provided $m = j$. In particular, $\hat{\phi}_0(0) = 1$ and $\hat{\phi}_1(0) = 0$. In other words the Taylor series of the entire function $\hat{\phi}_j(\zeta)$ starts with the power $j$.

Secondly, we can see from (3.2), (3.3), (3.4) that each $\hat{\phi}_j(\zeta)$ is a combination of powers of $1/\zeta$ and of exponents $e^{\pm i\zeta}$. More precisely, $\hat{\phi}_j(\zeta) = \Phi_j(-i\zeta)$, where

$$\Phi_j(x) = A_0^j \frac{1}{x} + \dots + A_j^j \frac{1}{x^{j+1}} + B_0^j \frac{e^x}{x} + \dots + B_j^j \frac{e^x}{x^{j+1}}.$$
Thus the nonstandard form matrices (2.15) take the final forms

$$[\sigma_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\phi_j(\zeta)\mathcal{F}\phi_p(-\zeta)d\zeta, \quad (4.2)$$

$$[\alpha_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\psi_j(\zeta)\mathcal{F}\psi_p(-\zeta)d\zeta, \quad (4.3)$$

$$[\beta_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\phi_j(\zeta)\mathcal{F}\psi_p(-\zeta)d\zeta, \quad (4.4)$$

$$[\gamma_l^n]_{pj} = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - itN^2\zeta^2) \mathcal{F}\psi_j(\zeta)\mathcal{F}\phi_p(-\zeta)d\zeta, \quad (4.5)$$

where the matrices are effectively sparse, as we shall see below, due to $\mathcal{F}\psi_j(0) = 0$ of at least order $\mathfrak{k}$. We point out here that Fourier transform $\mathcal{F}\psi_j(\zeta)$ can be easily found analytically and extended to the complex plane $\zeta \in \mathbb{C}$ as

$$\hat{\psi}_m(\zeta) = \frac{1}{\sqrt{2}} \sum_{j=0}^{\mathfrak{k}-1} \hat{\phi}_j\left(\frac{\zeta}{2}\right) \left[ g_{ij}^{(0)} + g_{ij}^{(1)} \exp\left(-\frac{i\zeta}{2}\right) \right],$$

where $g_{ij}^{(0)}$ and $g_{ij}^{(1)}$ are elements of the filter matrices $G^{(0)}$ and $G^{(1)}$, respectively. Justification of the contour deformation is straightforward, so we omit the proof.

Before we continue with evaluation of these integrals it is worth to make a couple of remarks on behaviour of the integrands in (4.2)-(4.5). As one can see from (3.2), (3.3), (3.4) that the Legendre Fourier transform extensions $\hat{\phi}_j(\zeta)$ have a removable singularity at zero $\zeta = 0$. Clearly, the same is true for the interpolating Fourier transform extensions $\hat{\phi}_j(\zeta)$ and the wavelet Fourier transform extensions $\hat{\psi}_j(\zeta)$. This suggests us to approximate these functions by power series around zero. Moreover, taking into account the Gaussian factor $\exp(-tN^2|\zeta|^2)$, one may expect that these power series will do a good job even away from zero $\zeta = 0$. Below we are mostly focused on further development of this idea. The last remark concerns the price one pays for the contour deformation. As it will be obvious below, apart from the Gaussian fast decreasing factor we also get the increasing factor $\exp((|l|+1)|\zeta|/\sqrt{2})$, that turns out to be crucial, due to round up error, when the time step $t$ is small and the scale $n$ is low. This problem is overcome in practice by the fact that we need the first integral (4.2) to be calculated with high precision at high scales $n$, depending on $t$. The rest are calculated by (3.7).

### 4.1. Legendre scaling functions.

$$[\sigma_l^n]_{pj}(a) = \frac{1}{2\pi} \int_{\Gamma} \exp(i\zeta l - ia\zeta^2) \hat{\phi}_j(\zeta) \hat{\phi}_p(-\zeta) d\zeta, \quad 0 \le p, j < \mathfrak{k}, \quad (4.6)$$

where $a = tN^2 = t4^n$. We can represent the multiplication of Fourier transforms as a $\zeta$-power series. This will transform this integral into a series of integrals $J_k(a, l)$, that can be evaluated exactly, with some coefficients depending solely on $p, j, k$. These coefficients need to be computed only once and tabulated, as they are problem-independent. We will call them correlation coefficients. We make an extensive use of (3.2) which is valid for integers $j \ge 1$.

Firstly, we notice that $\hat{\phi}_j(\zeta)$ has a root at zero of order $j$. Indeed,

$$\partial_{\zeta}^{m} \hat{\phi}_{j}(0) = (-i)^{m} \int_{\mathbb{R}} \phi_{j}(x) x^{m} dx$$

that equals zero for any $0 \le m < j$ and a non-zero provided $m = j$. In particular, $\hat{\phi}_0(0) = 1$ and $\hat{\phi}_1(0) = 0$. In other words the Taylor series of the entire function $\hat{\phi}_j(\zeta)$ starts with the power $j$.

Secondly, we can see from (3.2), (3.3), (3.4) that each $\hat{\phi}_j(\zeta)$ is a combination of powers of $1/\zeta$ and of exponents $e^{\pm i\zeta}$. More precisely, $\hat{\phi}_j(\zeta) = \Phi_j(-i\zeta)$, where

$$\Phi_j(x) = A_0^j \frac{1}{x} + \dots + A_j^j \frac{1}{x^{j+1}} + B_0^j \frac{e^x}{x} + \dots + B_j^j \frac{e^x}{x^{j+1}}.$$
These coefficients are real and can be found exactly. Indeed, from (3.3) we deduce

$$A_0^0 = -1, \quad B_0^0 = 1. \tag{4.7}$$

From (3.4) we deduce

$$A_0^1 = \sqrt{3}, \quad A_1^1 = 2\sqrt{3}, \quad B_0^1 = \sqrt{3}, \quad B_1^1 = -2\sqrt{3}. \tag{4.8}$$

Finally, for $j \ge 1$ from (3.2) we deduce the following relation

$$ \begin{align} A_0^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} A_0^{j-1} \notag \\ A_1^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} A_1^{j-1} - 2\sqrt{(2j+1)(2j+3)} A_0^j \notag \\ \multicolumn{4}{c}{\dots \quad \dots \quad \dots \quad \dots} \tag{4.9} \end{align} $$

$$A_{j-1}^{j+1} = \sqrt{\frac{2j+3}{2j-1}} A_{j-1}^{j-1} - 2\sqrt{(2j+1)(2j+3)} A_{j-2}^{j}$$

$$A_j^{j+1} = -2\sqrt{(2j+1)(2j+3)}A_{j-1}^j$$

$$A_{j+1}^{j+1} = -2\sqrt{(2j+1)(2j+3)}A_j^j$$

and $B_m^j$ obeying the same recurrence for $j \ge 1$.

Now we are ready to provide with the Taylor series about zero $\zeta = 0$ for the entire function $\hat{\phi}_j(\zeta)\hat{\phi}_p(-\zeta) = \Phi_j(-i\zeta)\Phi_p(i\zeta)$. For $j, p \ge 0$ we have

$$\Phi_j(-x)\Phi_p(x) = \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(-1)^{m+1}}{x^{m+q+2}} (A_m^j A_q^p + B_m^j B_q^p + A_m^j B_q^p e^x + B_m^j A_q^p e^{-x}).$$

Expanding the exponents $e^{\pm x}$ in powers of $x$ we obtain

$$\Phi_j(-x)\Phi_p(x) = \sum_{m=0}^{j} \sum_{q=0}^{p} \sum_{k=j+p+2+m+q}^{\infty} \frac{(-1)^{m+1}x^{k-m-q-2}}{k!} (A_m^j B_q^p + (-1)^k B_m^j A_q^p),$$

where we have discarded small powers, since the Taylor series starts from the power $j + p$ in $x$-variable. We can change the summation variable in the third sum as

$$\Phi_j(-x)\Phi_p(x) = \sum_{m=0}^{j} \sum_{q=0}^{p} \sum_{k=0}^{\infty} \frac{(-1)^{m+1}x^{k+j+p}}{(k+2+j+p+m+q)!} (A_m^j B_q^p + (-1)^{k+j+p+m+q} B_m^j A_q^p).$$

This series can be written down as

$$\Phi_j(-x)\Phi_p(x) = \sum_{k=0}^{\infty} \frac{C_{jp}^k x^{k+j+p}}{(k+2+j+p)!},$$

where

$$C_{jp}^k = \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(-1)^{m+1}(k+2+j+p)!}{(k+2+j+p+m+q)!} (A_m^j B_q^p + (-1)^{k+j+p+m+q} B_m^j A_q^p).$$

In particular, $C_{00}^k = 1 + (-1)^k$, that constitute all the cross correlation coefficients for the Haar multiresolution. There is an obvious problem here, namely, $|A_m^j| \sim (4j)^m$ and the same is true for $B_m^j$. So one may try to balance multiplication by factorising it in the following way

$$C_{jp}^k = \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(-1)^{m+1}(k+2+j+p)!(4j)^m(4p)^q}{(k+2+j+p+m+q)!} (\tilde{A}_m^j \tilde{B}_q^p + (-1)^{k+j+p+m+q} \tilde{B}_m^j \tilde{A}_q^p),$$
These coefficients are real and can be found exactly. Indeed, from (3.3) we deduce

$$A_0^0 = -1, \quad B_0^0 = 1. \tag{4.7}$$

From (3.4) we deduce

$$A_0^1 = \sqrt{3}, \quad A_1^1 = 2\sqrt{3}, \quad B_0^1 = \sqrt{3}, \quad B_1^1 = -2\sqrt{3}. \tag{4.8}$$

Finally, for $j \ge 1$ from (3.2) we deduce the following relation

$$ \begin{align} A_0^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} A_0^{j-1} \notag \\ A_1^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} A_1^{j-1} - 2\sqrt{(2j+1)(2j+3)} A_0^j \notag \\ \multicolumn{4}{c}{\dots \quad \dots \quad \dots \quad \dots} \tag{4.9} \end{align} $$

$$A_{j-1}^{j+1} = \sqrt{\frac{2j+3}{2j-1}} A_{j-1}^{j-1} - 2\sqrt{(2j+1)(2j+3)} A_{j-2}^{j}$$

$$A_j^{j+1} = -2\sqrt{(2j+1)(2j+3)}A_{j-1}^j$$

$$A_{j+1}^{j+1} = -2\sqrt{(2j+1)(2j+3)}A_j^j$$

and $B_m^j$ obeying the same recurrence for $j \ge 1$.

Now we are ready to provide with the Taylor series about zero $\zeta = 0$ for the entire function $\hat{\phi}_j(\zeta)\hat{\phi}_p(-\zeta) = \Phi_j(-i\zeta)\Phi_p(i\zeta)$. For $j, p \ge 0$ we have

$$\Phi_j(-x)\Phi_p(x) = \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(-1)^{m+1}}{x^{m+q+2}} (A_m^j A_q^p + B_m^j B_q^p + A_m^j B_q^p e^x + B_m^j A_q^p e^{-x}).$$

Expanding the exponents $e^{\pm x}$ in powers of $x$ we obtain

$$\Phi_j(-x)\Phi_p(x) = \sum_{m=0}^{j} \sum_{q=0}^{p} \sum_{k=j+p+2+m+q}^{\infty} \frac{(-1)^{m+1}x^{k-m-q-2}}{k!} (A_m^j B_q^p + (-1)^k B_m^j A_q^p),$$

where we have discarded small powers, since the Taylor series starts from the power $j + p$ in $x$-variable. We can change the summation variable in the third sum as

$$\Phi_j(-x)\Phi_p(x) = \sum_{m=0}^{j} \sum_{q=0}^{p} \sum_{k=0}^{\infty} \frac{(-1)^{m+1}x^{k+j+p}}{(k+2+j+p+m+q)!} \left(A_m^j B_q^p + (-1)^{k+j+p+m+q} B_m^j A_q^p\right).$$

This series can be written down as

$$\Phi_j(-x)\Phi_p(x) = \sum_{k=0}^{\infty} \frac{C_{jp}^k x^{k+j+p}}{(k+2+j+p)!},$$

where

$$C_{jp}^k = \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(-1)^{m+1}(k+2+j+p)!}{(k+2+j+p+m+q)!} \left(A_m^j B_q^p + (-1)^{k+j+p+m+q} B_m^j A_q^p\right).$$

In particular, $C_{00}^k = 1 + (-1)^k$, that constitute all the cross correlation coefficients for the Haar multiresolution. There is an obvious problem here, namely, $|A_m^j| \sim (4j)^m$ and the same is true for $B_m^j$. So one may try to balance multiplication by factorising it in the following way

$$C_{jp}^k = \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(-1)^{m+1}(k+2+j+p)!(4j)^m(4p)^q}{(k+2+j+p+m+q)!} \left(\tilde{A}_m^j \tilde{B}_q^p + (-1)^{k+j+p+m+q} \tilde{B}_m^j \tilde{A}_q^p\right),$$
with

$$ \tilde{A}_m^j = \frac{A_m^j}{(4j)^m}, \quad \tilde{B}_q^p = \frac{B_q^p}{(4p)^q}, \quad \text{and so on.} $$

These new coefficients satisfy the following relation

$$ \begin{aligned} \tilde{B}_0^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \tilde{B}_0^{j-1} \\ \tilde{B}_1^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \frac{j-1}{j+1} \tilde{B}_1^{j-1} - \sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \tilde{B}_0^j \\ \tilde{B}_2^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \left(\frac{j-1}{j+1}\right)^2 \tilde{B}_2^{j-1} - \sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \frac{j}{j+1} \tilde{B}_1^j \\ \dots & \dots \dots \dots \dots \\ \tilde{B}_{j-1}^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \left(\frac{j-1}{j+1}\right)^{j-1} \tilde{B}_{j-1}^{j-1} - \sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \left(\frac{j}{j+1}\right)^{j-2} \tilde{B}_{j-2}^j \\ \tilde{B}_j^{j+1} &= -\sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \left(\frac{j}{j+1}\right)^{j-1} \tilde{B}_{j-1}^j \\ \tilde{B}_{j+1}^{j+1} &= -\sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \left(\frac{j}{j+1}\right)^j \tilde{B}_j^j \end{aligned} \quad (4.10) $$

and $\tilde{A}_m^j$ obeying the same recurrence for $j \ge 1$. The first coefficients are

$$ \begin{aligned} \tilde{A}_0^0 &= -1, & \tilde{A}_0^1 &= \sqrt{3}, & \tilde{A}_1^1 &= \frac{\sqrt{3}}{2}, \\ \tilde{B}_0^0 &= 1, & \tilde{B}_0^1 &= \sqrt{3}, & \tilde{B}_1^1 &= -\frac{\sqrt{3}}{2}. \end{aligned} \quad (4.11) $$

These formulas are well balanced and can be easily implemented. It turns out that due to multiple self-cancellations, the correlation coefficients stay bounded. For all the ranges of $k, j, p$ of interest our simulations give $|C_{jp}^k| < 10$. Moreover, we still can simplify the expression for $C_{jp}^k$ even more, by exploiting some symmetries of A- and B-coefficients. For example, we can guarantee that $C_{jp}^{2k+1} = 0$ for all $k = 0, 1, \dots$.

**Lemma 2.** For any $j = 0, 1, \dots$ and $m = 0, 1, \dots, j$ it holds true that $A_m^j = (-1)^{1+j-m} B_m^j$.

*Proof.* The proof is split in two steps. Firstly, we claim the statement for $m = 0$, namely,

$$ A_0^j = (-1)^{1+j} B_0^j, \quad j \in \mathbb{N}_0. \quad (4.12) $$

Indeed, for $j = 0$ it follows from (4.7) and for $j = 1$ it follows from (4.8). For the other integers it follows by the induction

$$ A_0^{j+1} = c_1(0, j+1)A_0^{j-1} = c_1(0, j+1)(-1)^j B_0^{j-1} = (-1)^j B_0^{j+1} \quad (4.13) $$

The second step is to prove that

$$ A_z^{z+s} = (-1)^{1+s} B_z^{z+s} \quad (4.14) $$

for all $z, s \in \mathbb{N}_0$. One could proceed by using induction over two variables, but we can easily reduce the problem to the standard one variable induction. Indeed, let us consider a standard bijection between $n \in \mathbb{N}$ and $(z, s) \in \mathbb{N}_0^2$, respecting the increment $n \mapsto n+1$ either as

$$ (z, s) \mapsto (z+1, s-1) \quad (4.15) $$
with

$$ \tilde{A}_m^j = \frac{A_m^j}{(4j)^m}, \quad \tilde{B}_q^p = \frac{B_q^p}{(4p)^q}, \quad \text{and so on.} $$

These new coefficients satisfy the following relation

$$ \begin{aligned} \tilde{B}_0^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \tilde{B}_0^{j-1} \\ \tilde{B}_1^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \frac{j-1}{j+1} \tilde{B}_1^{j-1} - \sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \tilde{B}_0^j \\ \tilde{B}_2^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \left(\frac{j-1}{j+1}\right)^2 \tilde{B}_2^{j-1} - \sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \frac{j}{j+1} \tilde{B}_1^j \\ \dots & \dots \dots \dots \dots \\ \tilde{B}_{j-1}^{j+1} &= \sqrt{\frac{2j+3}{2j-1}} \left(\frac{j-1}{j+1}\right)^{j-1} \tilde{B}_{j-1}^{j-1} - \sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \left(\frac{j}{j+1}\right)^{j-2} \tilde{B}_{j-2}^j \\ \tilde{B}_j^{j+1} &= -\sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \left(\frac{j}{j+1}\right)^{j-1} \tilde{B}_{j-1}^j \\ \tilde{B}_{j+1}^{j+1} &= -\sqrt{\frac{(2j+1)(2j+3)}{(2j+2)(2j+2)}} \left(\frac{j}{j+1}\right)^j \tilde{B}_j^j \end{aligned} \quad (4.10) $$

and $\tilde{A}_m^j$ obeying the same recurrence for $j \ge 1$. The first coefficients are

$$ \begin{aligned} \tilde{A}_0^0 &= -1, & \tilde{A}_0^1 &= \sqrt{3}, & \tilde{A}_1^1 &= \frac{\sqrt{3}}{2}, \\ \tilde{B}_0^0 &= 1, & \tilde{B}_0^1 &= \sqrt{3}, & \tilde{B}_1^1 &= -\frac{\sqrt{3}}{2}. \end{aligned} \quad (4.11) $$

These formulas are well balanced and can be easily implemented. It turns out that due to multiple self-cancellations, the correlation coefficients stay bounded. For all the ranges of $k, j, p$ of interest our simulations give $|C_{jp}^k| < 10$. Moreover, we still can simplify the expression for $C_{jp}^k$ even more, by exploiting some symmetries of A- and B-coefficients. For example, we can guarantee that $C_{jp}^{2k+1} = 0$ for all $k = 0, 1, \dots$.

**Lemma 2.** For any $j = 0, 1, \dots$ and $m = 0, 1, \dots, j$ it holds true that $A_m^j = (-1)^{1+j-m} B_m^j$.

*Proof.* The proof is split in two steps. Firstly, we claim the statement for $m = 0$, namely,

$$ A_0^j = (-1)^{1+j} B_0^j, \quad j \in \mathbb{N}_0. \quad (4.12) $$

Indeed, for $j = 0$ it follows from (4.7) and for $j = 1$ it follows from (4.8). For the other integers it follows by the induction

$$ A_0^{j+1} = c_1(0, j+1)A_0^{j-1} = c_1(0, j+1)(-1)^j B_0^{j-1} = (-1)^j B_0^{j+1} \quad (4.13) $$

The second step is to prove that

$$ A_z^{z+s} = (-1)^{1+s} B_z^{z+s} \quad (4.14) $$

for all $z, s \in \mathbb{N}_0$. One could proceed by using induction over two variables, but we can easily reduce the problem to the standard one variable induction. Indeed, let us consider a standard bijection between $n \in \mathbb{N}$ and $(z, s) \in \mathbb{N}_0^2$, respecting the increment $n \mapsto n+1$ either as

$$ (z, s) \mapsto (z+1, s-1) \quad (4.15) $$
or as

$$z \mapsto 0 \text{ and } z + s \mapsto z + s + 1. \qquad (4.16)$$

Let $\mathcal{P}(n)$ stand for the statement (4.14) with $z(n), s(n)$. The induction base follows from the previous step (4.12). Let $\mathcal{P}(1), \dots, \mathcal{P}(n)$ hold true. We need to check the validity of $\mathcal{P}(n+1)$. If the increment $n \mapsto n+1$ corresponds to (4.16), then $s(n) = 0, s(n+1) = z(n)+1$ and $\mathcal{P}(n+1)$ stands for the statement

$$A_0^{z(n)+s(n)+1} = (-1)^{2+z(n)} B_0^{z(n)+s(n)+1}$$

holding true by (4.12). It is left to check $\mathcal{P}(n+1)$ for the case of the correspondence $n \mapsto n+1$ to (4.15), namely, we need to prove

$$A_{z(n)+1}^{z(n)+s(n)} = (-1)^{s(n)} B_{z(n)+1}^{z(n)+s(n)}.$$

This is obviously true when $s(n) = 1$ or $s(n) = 2$ due to (4.13). For other possible $z = z(n)$ and $s = s(n)$, i.e. $s \ge 3$, it follows from the recurrence relation (4.9) as

$$\begin{aligned} A_{z+1}^{z+s} &= c_1(z+1, z+s)A_{z+1}^{z+s-2} - c_2(z+1, z+s)A_z^{z+s-1} \\ &= c_1(z+1, z+s)(-1)^{s-2}B_{z+1}^{z+s-2} - c_2(z+1, z+s)(-1)^s B_z^{z+s-1} = (-1)^s B_{z+1}^{z+s} \end{aligned}$$

and the induction assumption on validity of the first $n$ statements.

□

From this lemma one can easily derive the final expression

$$C_{jp}^{k} = (-1)^{j} (1 + (-1)^{k}) \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(k + 2 + j + p)!(4j)^{m}(4p)^{q}}{(k + 2 + j + p + m + q)!} \tilde{B}_{m}^{j} \tilde{B}_{q}^{p}. \qquad (4.17)$$

In particular, $C_{jp}^k = 0$ for odd indices $k$. Finally, we show that these coefficients are uniformly bounded with respect to $k$.

**Lemma 3.** For any non-negative integers $k, j, p$ the following bound holds true

$$|C_{jp}^k| \le 2^{j+p+1} \sqrt{(2j+1)(2p+1)} \left\{ \begin{array}{ll} (j+1) \frac{(p/j)^{p+1}-1}{p/j-1}, & \text{provided } p=3(j-1), \\ (p+1) \frac{(j/p)^{j+1}-1}{j/p-1}, & \text{provided } j=3(p-1), \\ \frac{((4j)^{j+1} - (3+j+p)^{j+1})((4p)^{p+1} - (3+j+p)^{p+1})}{(3(j-1)-p)(3(p-1)-j)(3+j+p)^{j+p}}, & \text{other.} \end{array} \right.$$

*Proof.* We start by showing that

$$\max_{m=0,1,\dots,j} |\tilde{B}_m^j| \le 2^j \sqrt{2j+1}$$

for each $j \in \mathbb{N}_0$. Indeed, $\tilde{B}_0^j$ can be calculated directly by the first line in the recurrence relation (4.10) with the first coefficients (4.11). It yields

$$\tilde{B}_0^j = \sqrt{2j+1}, \quad j \in \mathbb{N}_0.$$

The last two lines in (4.10) give

$$|\tilde{B}_{j-1}^j| \le |\tilde{B}_0^1|, \quad |\tilde{B}_j^j| \le |\tilde{B}_1^1|, \quad j \in \mathbb{N}.$$

In particular,

$$|\tilde{B}_{j-1}^j|, |\tilde{B}_j^j| \le \sqrt{2j+1}, \quad j \in \mathbb{N}.$$

Now let us consider an array $M_z^j$ with $z = 0, \dots, j$ and $j \in \mathbb{N}_0$ satisfying

$$M_z^{j+1} = \left(\frac{j-1}{j+1}\right)^z M_z^{j-1} + \frac{2j+1}{2j+2} \left(\frac{j}{j+1}\right)^{z-1} M_{z-1}^j, \quad z = 1, \dots, j-1, \quad j \ge 2,$$
or as

$$z \mapsto 0 \text{ and } z + s \mapsto z + s + 1. \tag{4.16}$$

Let $\mathcal{P}(n)$ stand for the statement (4.14) with $z(n), s(n)$. The induction base follows from the previous step (4.12). Let $\mathcal{P}(1), \dots, \mathcal{P}(n)$ hold true. We need to check the validity of $\mathcal{P}(n+1)$. If the increment $n \mapsto n+1$ corresponds to (4.16), then $s(n) = 0$, $s(n+1) = z(n)+1$ and $\mathcal{P}(n+1)$ stands for the statement

$$A_0^{z(n)+s(n)+1} = (-1)^{2+z(n)} B_0^{z(n)+s(n)+1}$$

holding true by (4.12). It is left to check $\mathcal{P}(n+1)$ for the case of the correspondence $n \mapsto n+1$ to (4.15), namely, we need to prove

$$A_{z(n)+1}^{z(n)+s(n)} = (-1)^{s(n)} B_{z(n)+1}^{z(n)+s(n)}.$$

This is obviously true when $s(n) = 1$ or $s(n) = 2$ due to (4.13). For other possible $z = z(n)$ and $s = s(n)$, i.e. $s \ge 3$, it follows from the recurrence relation (4.9) as

$$\begin{aligned} A_{z+1}^{z+s} &= c_1(z+1, z+s)A_{z+1}^{z+s-2} - c_2(z+1, z+s)A_z^{z+s-1} \\ &= c_1(z+1, z+s)(-1)^{s-2}B_{z+1}^{z+s-2} - c_2(z+1, z+s)(-1)^s B_z^{z+s-1} = (-1)^s B_{z+1}^{z+s} \end{aligned}$$

and the induction assumption on validity of the first $n$ statements.

□

From this lemma one can easily derive the final expression

$$C_{jp}^{k} = (-1)^{j} (1 + (-1)^{k}) \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(k + 2 + j + p)!(4j)^{m}(4p)^{q}}{(k + 2 + j + p + m + q)!} \tilde{B}_{m}^{j} \tilde{B}_{q}^{p}. \tag{4.17}$$

In particular, $C_{jp}^k = 0$ for odd indices $k$. Finally, we show that these coefficients are uniformly bounded with respect to $k$.

**Lemma 3.** For any non-negative integers $k, j, p$ the following bound holds true

$$\left| C_{jp}^{k} \right| \le 2^{j+p+1} \sqrt{(2j+1)(2p+1)} \left\{ \begin{array}{ll} (j+1) \frac{(p/j)^{p+1}-1}{p/j-1}, & \text{provided } p=3(j-1), \\ (p+1) \frac{(j/p)^{j+1}-1}{j/p-1}, & \text{provided } j=3(p-1), \\ \frac{((4j)^{j+1} - (3+j+p)^{j+1}) ((4p)^{p+1} - (3+j+p)^{p+1})}{(3(j-1)-p)(3(p-1)-j)(3+j+p)^{j+p}}, & \text{other.} \end{array} \right.$$

*Proof.* We start by showing that

$$\max_{m=0,1,\dots,j} \left| \tilde{B}_m^j \right| \le 2^j \sqrt{2j+1}$$

for each $j \in \mathbb{N}_0$. Indeed, $\tilde{B}_0^j$ can be calculated directly by the first line in the recurrence relation (4.10) with the first coefficients (4.11). It yields

$$\tilde{B}_0^j = \sqrt{2j+1}, \quad j \in \mathbb{N}_0.$$

The last two lines in (4.10) give

$$\left| \tilde{B}_{j-1}^j \right| \le \left| \tilde{B}_0^1 \right|, \quad \left| \tilde{B}_j^j \right| \le \left| \tilde{B}_1^1 \right|, \quad j \in \mathbb{N}.$$

In particular,

$$\left| \tilde{B}_{j-1}^j \right|, \left| \tilde{B}_j^j \right| \le \sqrt{2j+1}, \quad j \in \mathbb{N}.$$

Now let us consider an array $M_z^j$ with $z = 0, \dots, j$ and $j \in \mathbb{N}_0$ satisfying

$$M_z^{j+1} = \left( \frac{j-1}{j+1} \right)^z M_z^{j-1} + \frac{2j+1}{2j+2} \left( \frac{j}{j+1} \right)^{z-1} M_{z-1}^j, \quad z = 1, \dots, j-1, \quad j \ge 2,$$
where

$$M_0^j = M_j^j = M_j^{j+1} = 1, \quad j \in \mathbb{N}_0.$$

Then it is easy to see that

$$|\tilde{B}_z^j| \le M_z^j \sqrt{2j+1}$$

holds true for all admissible pairs $z, j$. These new parameters can be roughly estimated by

$$M_z^{j+2} \le \max_{z=0, \dots, j+1} M_z^{j+1} + \max_{z=0, \dots, j} M_z^j \le 2^j, \quad z = 1, \dots, j, \quad j \ge 0,$$

which leads to the bound for $\tilde{B}_m^j$ claimed above at the beginning of the proof. We finish now by noticing

$$\sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(k+2+j+p)!(4j)^m (4p)^q}{(k+2+j+p+m+q)!} \le \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(4j)^m (4p)^q}{(3+j+p)^{m+q}}$$

that is a multiplication of two finite geometric series. Summing these geometric series and using the bound for $\tilde{B}_m^j$ we conclude the proof by (4.17). □

Now we can make a use of the obtained power series for $\hat{\phi}_j(\zeta)\hat{\phi}_p(-\zeta)$ standing in (4.6). The contour $\Gamma$ is parameterised by $\zeta = \rho e^{-i\pi/4}$, $\rho \in \mathbb{R}$. One can change the integration variable, and then exchange integration and summation by appealing to the dominated convergence theorem. Thus our main operator (4.6) has the following expansion

$$[\sigma_l^n]_{pj}(a) = \sum_{k=0}^{\infty} C_{jp}^{2k} J_{2k+j+p}(l, a), \qquad (4.18)$$

where $a = tN^2 = t4^n$ and

$$J_m = \frac{e^{i\pi(m-1)/4}}{2\pi(m+2)!} \int_{\mathbb{R}} \exp\left(\rho l \exp\left(i\frac{\pi}{4}\right) - a\rho^2\right) \rho^m d\rho \qquad (4.19)$$

satisfying the following relation

$$J_{m+1} = \frac{i}{2a(m+3)} \left( lJ_m + \frac{m}{(m+2)}J_{m-1} \right), \quad m = 0, 1, 2, \dots, \qquad (4.20)$$

with $J_{-1} = 0$ and

$$J_0 = \frac{e^{-i\frac{\pi}{4}}}{4\sqrt{\pi a}} \exp\left(\frac{il^2}{4a}\right). \qquad (4.21)$$

In other words, the zero-power integral $J_0$ is standard and the rest are obtained from it via integration by parts. Expansion (4.18) turns out to be very efficient as long as the scaling order $n$ is big enough depending on how small is $t$. In other words, the smaller time $t$ the bigger $n$ should be, in order that $|J_m|$ does not reach large values that can destroy the accuracy due to rounding error. For example, $n$ should be bigger than 12 for $t = 10^{-5}$.

Using this recurrence formula we can prove some very useful error estimates. It is worth to point out that $C_{jp}^k \in \ell^\infty(\mathbb{N}_0)$, and so it makes sense to estimate the sequence of power integrals in $\ell^1(\mathbb{N}_0)$. Firstly, it is clear that $J_m \to 0$, since the series converges. Appealing repetitively to the recurrence relation we can get that $J_m$ tends to zero faster than any power of $1/m$. Denoting $b_m = |J_m|$ we can notice that

$$b_{m+1} \le \frac{|l|}{2a(m+3)} b_m + \frac{m}{2a(m+2)(m+3)} b_{m-1}$$

$$b_{m+2} \le \frac{|l|}{2a(m+4)} b_{m+1} + \frac{m+1}{2a(m+3)(m+4)} b_m$$

$$b_{m+3} \le \frac{|l|}{2a(m+5)} b_{m+2} + \frac{m+2}{2a(m+4)(m+5)} b_{m+1}$$

... ... ... ... ...
where

$$M_0^j = M_j^j = M_j^{j+1} = 1, \quad j \in \mathbb{N}_0.$$

Then it is easy to see that

$$|\tilde{B}_z^j| \le M_z^j \sqrt{2j+1}$$

holds true for all admissible pairs $z, j$. These new parameters can be roughly estimated by

$$M_z^{j+2} \le \max_{z=0, \dots, j+1} M_z^{j+1} + \max_{z=0, \dots, j} M_z^j \le 2^j, \quad z = 1, \dots, j, \quad j \ge 0,$$

which leads to the bound for $\tilde{B}_m^j$ claimed above at the beginning of the proof. We finish now by noticing

$$\sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(k+2+j+p)!(4j)^m (4p)^q}{(k+2+j+p+m+q)!} \le \sum_{m=0}^{j} \sum_{q=0}^{p} \frac{(4j)^m (4p)^q}{(3+j+p)^{m+q}}$$

that is a multiplication of two finite geometric series. Summing these geometric series and using the bound for $\tilde{B}_m^j$ we conclude the proof by (4.17). □

Now we can make a use of the obtained power series for $\hat{\phi}_j(\zeta)\hat{\phi}_p(-\zeta)$ standing in (4.6). The contour $\Gamma$ is parameterised by $\zeta = \rho e^{-i\pi/4}$, $\rho \in \mathbb{R}$. One can change the integration variable, and then exchange integration and summation by appealing to the dominated convergence theorem. Thus our main operator (4.6) has the following expansion

$$[\sigma_l^n]_{pj}(a) = \sum_{k=0}^{\infty} C_{jp}^{2k} J_{2k+j+p}(l, a), \qquad (4.18)$$

where $a = tN^2 = t4^n$ and

$$J_m = \frac{e^{i\pi(m-1)/4}}{2\pi(m+2)!} \int_{\mathbb{R}} \exp\left(\rho l \exp\left(i\frac{\pi}{4}\right) - a\rho^2\right) \rho^m d\rho \qquad (4.19)$$

satisfying the following relation

$$J_{m+1} = \frac{i}{2a(m+3)} \left( lJ_m + \frac{m}{(m+2)}J_{m-1} \right), \quad m = 0, 1, 2, \dots, \qquad (4.20)$$

with $J_{-1} = 0$ and

$$J_0 = \frac{e^{-i\frac{\pi}{4}}}{4\sqrt{\pi a}} \exp\left(\frac{il^2}{4a}\right). \qquad (4.21)$$

In other words, the zero-power integral $J_0$ is standard and the rest are obtained from it via integration by parts. Expansion (4.18) turns out to be very efficient as long as the scaling order $n$ is big enough depending on how small is $t$. In other words, the smaller time $t$ the bigger $n$ should be, in order that $|J_m|$ does not reach large values that can destroy the accuracy due to rounding error. For example, $n$ should be bigger than 12 for $t = 10^{-5}$.

Using this recurrence formula we can prove some very useful error estimates. It is worth to point out that $C_{jp}^k \in \ell^\infty(\mathbb{N}_0)$, and so it makes sense to estimate the sequence of power integrals in $\ell^1(\mathbb{N}_0)$. Firstly, it is clear that $J_m \to 0$, since the series converges. Appealing repetitively to the recurrence relation we can get that $J_m$ tends to zero faster than any power of $1/m$. Denoting $b_m = |J_m|$ we can notice that

$$b_{m+1} \le \frac{|l|}{2a(m+3)} b_m + \frac{m}{2a(m+2)(m+3)} b_{m-1}$$

$$b_{m+2} \le \frac{|l|}{2a(m+4)} b_{m+1} + \frac{m+1}{2a(m+3)(m+4)} b_m$$

$$b_{m+3} \le \frac{|l|}{2a(m+5)} b_{m+2} + \frac{m+2}{2a(m+4)(m+5)} b_{m+1}$$

... ... ... ... ...
Summing these inequalities one obtains

$$b_{m+1} + b_{m+2} + \dots \le \frac{m}{2a(m+2)(m+3)} b_{m-1} + \frac{1}{2a(m+3)} \left( |l| + 1 - \frac{3}{m+4} \right) b_m + \frac{1}{2a(m+4)} \left( |l| + 1 - \frac{3}{m+5} \right) b_{m+1} + \dots$$

One can notice that the coefficient standing in front of $b_{m+1}$ forms a decreasing sequence (strictly decreasing for either $m \ge 1$ or $l \ne 0$), and so we can take it out of the brackets. Moreover, for $m \ge \frac{|l|+1}{2a} - \frac{9}{2} + \sqrt{\left(\frac{|l|+1}{2a} + \frac{1}{2}\right)^2 - \frac{3}{a}}$, or any $m$ in case the value under the square root is negative, we have

$$b_{m+1} + b_{m+2} + \dots \le \frac{m}{a(m+2)(m+3)} b_{m-1} + \frac{1}{a(m+3)} \left( |l| + 1 - \frac{3}{m+4} \right) b_m$$

that can be used in practice for series termination. We can precise it even a bit more, since we only need to sum $b_m$ with either even or odd indices $m$, as follows

$$b_{m+1} + b_{m+3} + \dots \le \frac{m}{2a(m+2)(m+3)} b_{m-1} + \frac{|l|}{2a(m+3)} b_m + \frac{m+2}{2a(m+4)(m+5)} b_{m+1} + \frac{|l|}{2a(m+5)} b_{m+2} + \dots \le \frac{m}{2a(m+2)(m+3)} b_{m-1} + \frac{|l|}{2a(m+3)} b_m + \frac{|l|+1}{2a(m+5)} \left[ \frac{m}{a(m+2)(m+3)} b_{m-1} + \frac{1}{a(m+3)} \left( |l| + 1 - \frac{3}{m+4} \right) b_m \right].$$

Given some precision, one may use these inequalities to cut the infinite sum in (4.18) in practical calculations. In principal, we could proceed further and try to obtain an optimal precision-dependent expression for the series cut theoretically as well, though for big $|l|$ we anticipate it to be way too excessive. Moreover, the bound in Lemma 3 turns out to be very rough in practice, due to very fast convergence of the integral series. In normal practical calculations it is enough to know that $|C_{jp}^k| < 10$ for $k \le 50$.

Finally, we turn our attention to the sparsity of matrices (4.3)-(4.5) associated with the non-standard form (2.14) by means of (2.13), (2.15). One can calculate matrices $\sigma_l^{n+1}, l = 1 - 2^{n+1}, \dots, 2^{n+1} - 1$, using the expansion (4.18) and then find matrices $\alpha_l^n, \beta_l^n, \gamma_l^n$ with $l = 1 - 2^n, \dots, 2^n - 1$ exploiting (3.7). The Frobenius norm

$$\|\alpha\| = \left( \sum_{p,j=0}^{\ell-1} |[\alpha]_{pj}|^2 \right)^{1/2} \qquad (4.22)$$

is used for a matrix $\alpha$ of $\ell \times \ell$-size. Figures 1-6 demonstrate the tendency of considerable non-zero values $\|\alpha_l^n\|, \|\beta_l^n\|$ to stay away from the main diagonal $l = 0$. It is in contrast with the fact that the non-standard form matrices of the heat equation are narrow diagonal bounded as we mentioned in Subsection 3.2. Intuitively, this difference can be explained as follows. The heat defuses with time, that is, it transports away from the initial local perturbation point with some smearing. Waves described by the Schrödinger equation (1.3) disperse. Such a wave moves away from the initial local perturbation point forming high frequency ripples that are captured by multiresolution analysis.

Alternatively, the values on the diagonal and away could be roughly estimated as follows. Due to the vanishing moment property (2.7), $\hat{\psi}_j(\zeta)$ has a root at zero of order $\ell + j$. It results in the fact, that the corresponding expansions, similar to (4.18), for the matrix elements $[\alpha_l^n]_{pj}$ and $[\beta_l^n]_{pj}$ start with integrals $J_{2\ell+p+j}$ and $J_{\ell+p+j}$, respectively. These integrals constitute the leading terms, due to their fast convergence to zero. Moreover, the corner elements $[\alpha_l^n]_{00}$ and $[\beta_l^n]_{00}$ are dominant for the same reason. Therefore, $\|\alpha_l^n\| \sim |J_{2\ell}(l)|$ and $\|\beta_l^n\| \sim |J_\ell(l)|$.
Summing these inequalities one obtains

$$b_{m+1} + b_{m+2} + \dots \le \frac{m}{2a(m+2)(m+3)} b_{m-1} + \frac{1}{2a(m+3)} \left( |l| + 1 - \frac{3}{m+4} \right) b_m + \frac{1}{2a(m+4)} \left( |l| + 1 - \frac{3}{m+5} \right) b_{m+1} + \dots$$

One can notice that the coefficient standing in front of $b_{m+1}$ forms a decreasing sequence (strictly decreasing for either $m \ge 1$ or $l \ne 0$), and so we can take it out of the brackets. Moreover, for $m \ge \frac{|l|+1}{2a} - \frac{9}{2} + \sqrt{\left(\frac{|l|+1}{2a} + \frac{1}{2}\right)^2 - \frac{3}{a}}$, or any $m$ in case the value under the square root is negative, we have

$$b_{m+1} + b_{m+2} + \dots \le \frac{m}{a(m+2)(m+3)} b_{m-1} + \frac{1}{a(m+3)} \left( |l| + 1 - \frac{3}{m+4} \right) b_m$$

that can be used in practice for series termination. We can precise it even a bit more, since we only need to sum $b_m$ with either even or odd indices $m$, as follows

$$b_{m+1} + b_{m+3} + \dots \le \frac{m}{2a(m+2)(m+3)} b_{m-1} + \frac{|l|}{2a(m+3)} b_m + \frac{m+2}{2a(m+4)(m+5)} b_{m+1} + \frac{|l|}{2a(m+5)} b_{m+2} + \dots \le \frac{m}{2a(m+2)(m+3)} b_{m-1} + \frac{|l|}{2a(m+3)} b_m + \frac{|l|+1}{2a(m+5)} \left[ \frac{m}{a(m+2)(m+3)} b_{m-1} + \frac{1}{a(m+3)} \left( |l| + 1 - \frac{3}{m+4} \right) b_m \right].$$

Given some precision, one may use these inequalities to cut the infinite sum in (4.18) in practical calculations. In principal, we could proceed further and try to obtain an optimal precision-dependent expression for the series cut theoretically as well, though for big $|l|$ we anticipate it to be way too excessive. Moreover, the bound in Lemma 3 turns out to be very rough in practice, due to very fast convergence of the integral series. In normal practical calculations it is enough to know that $|C_{jp}^k| < 10$ for $k \le 50$.

Finally, we turn our attention to the sparsity of matrices (4.3)-(4.5) associated with the non-standard form (2.14) by means of (2.13), (2.15). One can calculate matrices $\sigma_l^{n+1}, l = 1 - 2^{n+1}, \dots, 2^{n+1} - 1$, using the expansion (4.18) and then find matrices $\alpha_l^n, \beta_l^n, \gamma_l^n$ with $l = 1 - 2^n, \dots, 2^n - 1$ exploiting (3.7). The Frobenius norm

$$\|\alpha\| = \left( \sum_{p,j=0}^{\ell-1} |[\alpha]_{pj}|^2 \right)^{1/2} \qquad (4.22)$$

is used for a matrix $\alpha$ of $\ell \times \ell$-size. Figures 1 - 6 demonstrate the tendency of considerable non-zero values $\|\alpha_l^n\|, \|\beta_l^n\|$ to stay away from the main diagonal $l = 0$. It is in contrast with the fact that the non-standard form matrices of the heat equation are narrow diagonal bounded as we mentioned in Subsection 3.2. Intuitively, this difference can be explained as follows. The heat defuses with time, that is, it transports away from the initial local perturbation point with some smearing. Waves described by the Schrödinger equation (1.3) disperse. Such a wave moves away from the initial local perturbation point forming high frequency ripples that are captured by multiresolution analysis.

Alternatively, the values on the diagonal and away could be roughly estimated as follows. Due to the vanishing moment property (2.7), $\hat{\psi}_j(\zeta)$ has a root at zero of order $\ell + j$. It results in the fact, that the corresponding expansions, similar to (4.18), for the matrix elements $[\alpha_l^n]_{pj}$ and $[\beta_l^n]_{pj}$ start with integrals $J_{2\ell+p+j}$ and $J_{\ell+p+j}$, respectively. These integrals constitute the leading terms, due to their fast convergence to zero. Moreover, the corner elements $[\alpha_l^n]_{00}$ and $[\beta_l^n]_{00}$ are dominant for the same reason. Therefore, $\|\alpha_l^n\| \sim |J_{2\ell}(l)|$ and $\|\beta_l^n\| \sim |J_{\ell}(l)|$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Sparsity of the non-standard form at scale $n = 7$ corresponding to the time step $t = 0.0001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $\mathfrak{k} = 6$, dashed to $\mathfrak{k} = 2$ and dash-dot to $\mathfrak{k} = 11$. For $\mathfrak{k} = 2$ we have the values on the diagonal $||\alpha_0|| = 8.0\text{e-05}$, $||\beta_0|| = 2.4\text{e-03}$, and in the corner $||\alpha_{2n-1}|| = 5.0\text{e-03}$, $||\beta_{2n-1}|| = 2.6\text{e-03}$. For $\mathfrak{k} = 6$ we have the values on the diagonal $||\alpha_0|| = 3.2\text{e-13}$, $||\beta_0|| = 3.9\text{e-08}$, and in the corner $||\alpha_{2n-1}|| = 3.2\text{e-02}$, $||\beta_{2n-1}|| = 1.7\text{e-02}$. For $\mathfrak{k} = 11$ we have the values on the diagonal $||\alpha_0|| = 2.7\text{e-18}$, $||\beta_0|| = 4.0\text{e-15}$, and in the corner $||\alpha_{2n-1}|| = 1.8\text{e-01}$, $||\beta_{2n-1}|| = 8.1\text{e-02}$.</figcaption></figure>

For $l = 0$ from (4.20) one obtains

$$J_{2k}(0) = \frac{e^{i\pi(2k-1)/4}}{2\sqrt{\pi a}(2a)^k(2k+2)!!(2k+1)}, \quad J_{2k+1}(0) = 0. \qquad (4.23)$$

In the case of $|l| \gg 1$, we can admit the following approximation

$$J_{m+1}(l) \approx \frac{il}{2a(m+3)} J_m(l)$$

and so

$$J_m(l) \approx \left(\frac{il}{2a}\right)^m \frac{\exp\left(\frac{i\ell^2}{4a} - \frac{i\pi}{4}\right)}{2\sqrt{\pi a(m+2)!}} \qquad (4.24)$$

which leads to

$$\|\alpha_l^n\| \lesssim \frac{1}{2\sqrt{\pi a(2\mathfrak{k}+2)!}} \left(\frac{l}{2a}\right)^{2\mathfrak{k}}, \quad |l| \gg 1$$

Meanwhile, $\|\alpha_0^n\| \sim |J_{2\mathfrak{k}}(0)|$, where the integral is calculated in (4.23). Thus

$$\|\alpha_l^n\| \lesssim \frac{1}{2\sqrt{\pi a(2a)^{\mathfrak{k}}(2\mathfrak{k}+2)!!(2\mathfrak{k}+1)}} + \frac{1}{2\sqrt{\pi a(2\mathfrak{k}+2)!}} \left(\frac{l}{2a}\right)^{2\mathfrak{k}}$$

Similarly,

$$\|\beta_l^n\| = \|\gamma_l^n\| \lesssim \frac{1}{2\sqrt{\pi a(2a)^{\mathfrak{k}/2}(\mathfrak{k}+2)!!(\mathfrak{k}+1)}} + \frac{1}{2\sqrt{\pi a(\mathfrak{k}+2)!}} \left(\frac{l}{2a}\right)^{\mathfrak{k}}, \quad \text{provided } \mathfrak{k} \text{ is even}$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Sparsity of the non-standard form at scale $n = 7$ corresponding to the time step $t = 0.0001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $\mathfrak{k} = 6$, dashed to $\mathfrak{k} = 2$ and dash-dot to $\mathfrak{k} = 11$. For $\mathfrak{k} = 2$ we have the values on the diagonal $||\alpha_0|| = 8.0e-05$, $||\beta_0|| = 2.4e-03$, and in the corner $||\alpha_{2^{n-1}}|| = 5.0e-03$, $||\beta_{2^{n-1}}|| = 2.6e-03$. For $\mathfrak{k} = 6$ we have the values on the diagonal $||\alpha_0|| = 3.2e-13$, $||\beta_0|| = 3.9e-08$, and in the corner $||\alpha_{2^{n-1}}|| = 3.2e-02$, $||\beta_{2^{n-1}}|| = 1.7e-02$. For $\mathfrak{k} = 11$ we have the values on the diagonal $||\alpha_0|| = 2.7e-18$, $||\beta_0|| = 4.0e-15$, and in the corner $||\alpha_{2^{n-1}}|| = 1.8e-01$, $||\beta_{2^{n-1}}|| = 8.1e-02$.</figcaption></figure>

For $l = 0$ from (4.20) one obtains

$$J_{2k}(0) = \frac{e^{i\pi(2k-1)/4}}{2\sqrt{\pi a}(2a)^k(2k+2)!!(2k+1)}, \quad J_{2k+1}(0) = 0. \qquad (4.23)$$

In the case of $|l| \gg 1$, we can admit the following approximation

$$J_{m+1}(l) \approx \frac{il}{2a(m+3)} J_m(l)$$

and so

$$J_m(l) \approx \left(\frac{il}{2a}\right)^m \frac{\exp\left(\frac{i\ell^2}{4a} - \frac{i\pi}{4}\right)}{2\sqrt{\pi a(m+2)!}} \qquad (4.24)$$

which leads to

$$\|\alpha_l^n\| \lesssim \frac{1}{2\sqrt{\pi a(2\mathfrak{k}+2)!}} \left(\frac{l}{2a}\right)^{2\mathfrak{k}}, \quad |l| \gg 1$$

Meanwhile, $\|\alpha_0^n\| \sim |J_{2\mathfrak{k}}(0)|$, where the integral is calculated in (4.23). Thus

$$\|\alpha_l^n\| \lesssim \frac{1}{2\sqrt{\pi a(2a)^{\mathfrak{k}}(2\mathfrak{k}+2)!!(2\mathfrak{k}+1)}} + \frac{1}{2\sqrt{\pi a(2\mathfrak{k}+2)!}} \left(\frac{l}{2a}\right)^{2\mathfrak{k}}$$

Similarly,

$$\|\beta_l^n\| = \|\gamma_l^n\| \lesssim \frac{1}{2\sqrt{\pi a(2a)^{\mathfrak{k}/2}(\mathfrak{k}+2)!!(\mathfrak{k}+1)}} + \frac{1}{2\sqrt{\pi a(\mathfrak{k}+2)!}} \left(\frac{l}{2a}\right)^{\mathfrak{k}}, \quad \text{provided } \mathfrak{k} \text{ is even}$$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Sparsity of the non-standard form at scale $n = 8$ corresponding to the time step $t = 0.0001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $\mathfrak{k} = 6$, dashed to $\mathfrak{k} = 2$ and dash-dot to $\mathfrak{k} = 11$. For $\mathfrak{k} = 2$ we have the values on the diagonal $||\alpha_0|| = 2.5e-06$, $||\beta_0|| = 3.0e-04$, and in the corner $||\alpha_{2^{n-1}}|| = 3.2e-03$, $||\beta_{2^{n-1}}|| = 3.2e-03$. For $\mathfrak{k} = 6$ we have the values on the diagonal $||\alpha_0|| = 3.9e-17$, $||\beta_0|| = 3.1e-10$, and in the corner $||\alpha_{2^{n-1}}|| = 8.6e-02$, $||\beta_{2^{n-1}}|| = 4.2e-02$. For $\mathfrak{k} = 11$ we have the values on the diagonal $||\alpha_0|| = 5.2e-19$, $||\beta_0|| = 3.5e-18$, and in the corner $||\alpha_{2^{n-1}}|| = 2.9e-03$, $||\beta_{2^{n-1}}|| = 1.8e-02$.</figcaption></figure>

and

$$ ||\beta_l^n|| = ||\gamma_l^n|| \lesssim \frac{1}{2\sqrt{\pi a}(2a)^{(\mathfrak{k}+1)/2}(\mathfrak{k}+3)!!(\mathfrak{k}+2)} + \frac{1}{2\sqrt{\pi a}(\mathfrak{k}+2)!} \left(\frac{l}{2a}\right)^{\mathfrak{k}}, \quad \text{provided } \mathfrak{k} \text{ is odd.} $$

Note that the expansions for matrices $\alpha_l^n, \beta_l^n$ with respect to the power integrals contain cross correlation coefficients depending on the MRA order $\mathfrak{k}$. Therefore, the implicit constants staying in the obtained inequalities for $||\alpha_l^n||, ||\beta_l^n||$ depend on $\mathfrak{k}$ as well. In other words, these inequalities provide only with a qualitative behaviour of the norms. For instance, one can compare

$$ |J_{2\mathfrak{k}}(0)| = \begin{cases} 8.6e-05, & \mathfrak{k} = 2 \\ 2.1e-11, & \mathfrak{k} = 6 \\ 1.0e-20, & \mathfrak{k} = 11 \end{cases} $$

for $n = 7$ and $t = 0.0001$ with the values $||\alpha_0^n||$ given in the caption of Figure 1.

4.2. Interpolating scaling functions. Similarly to the case of Legendre scaling functions, we introduce functions $\Psi_j(-i\zeta) = \hat{\varphi}_j(\zeta)$, so that

$$ \Psi_j(x) = \sqrt{w_j} \sum_{m=0}^{\mathfrak{k}-1} \phi_m(x_j) \Phi_m(x), \quad x \in (0, 1), \quad j = 0, \dots, \mathfrak{k}-1, $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Sparsity of the non-standard form at scale $n = 8$ corresponding to the time step $t = 0.0001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $\mathfrak{k} = 6$, dashed to $\mathfrak{k} = 2$ and dash-dot to $\mathfrak{k} = 11$. For $\mathfrak{k} = 2$ we have the values on the diagonal $||\alpha_0|| = 2.5e-06$, $||\beta_0|| = 3.0e-04$, and in the corner $||\alpha_{2^{n-1}}|| = 3.2e-03$, $||\beta_{2^{n-1}}|| = 3.2e-03$. For $\mathfrak{k} = 6$ we have the values on the diagonal $||\alpha_0|| = 3.9e-17$, $||\beta_0|| = 3.1e-10$, and in the corner $||\alpha_{2^{n-1}}|| = 8.6e-02$, $||\beta_{2^{n-1}}|| = 4.2e-02$. For $\mathfrak{k} = 11$ we have the values on the diagonal $||\alpha_0|| = 5.2e-19$, $||\beta_0|| = 3.5e-18$, and in the corner $||\alpha_{2^{n-1}}|| = 2.9e-03$, $||\beta_{2^{n-1}}|| = 1.8e-02$.</figcaption></figure>

and

$$ ||\beta_l^n|| = ||\gamma_l^n|| \lesssim \frac{1}{2\sqrt{\pi a}(2a)^{(\mathfrak{k}+1)/2}(\mathfrak{k}+3)!!(\mathfrak{k}+2)} + \frac{1}{2\sqrt{\pi a}(\mathfrak{k}+2)!} \left(\frac{l}{2a}\right)^{\mathfrak{k}}, \quad \text{provided } \mathfrak{k} \text{ is odd.} $$

Note that the expansions for matrices $\alpha_l^n, \beta_l^n$ with respect to the power integrals contain cross correlation coefficients depending on the MRA order $\mathfrak{k}$. Therefore, the implicit constants staying in the obtained inequalities for $||\alpha_l^n||, ||\beta_l^n||$ depend on $\mathfrak{k}$ as well. In other words, these inequalities provide only with a qualitative behaviour of the norms. For instance, one can compare

$$ |J_{2\mathfrak{k}}(0)| = \begin{cases} 8.6e-05, & \mathfrak{k} = 2 \\ 2.1e-11, & \mathfrak{k} = 6 \\ 1.0e-20, & \mathfrak{k} = 11 \end{cases} $$

for $n = 7$ and $t = 0.0001$ with the values $||\alpha_0^n||$ given in the caption of Figure 1.

4.2. Interpolating scaling functions. Similarly to the case of Legendre scaling functions, we introduce functions $\Psi_j(-i\zeta) = \hat{\varphi}_j(\zeta)$, so that

$$ \Psi_j(x) = \sqrt{w_j} \sum_{m=0}^{\mathfrak{k}-1} \phi_m(x_j) \Phi_m(x), \quad x \in (0,1), \quad j = 0, \dots, \mathfrak{k}-1, $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Sparsity of the non-standard form at scale $n = 9$ corresponding to the time step $t = 0.0001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $\mathfrak{k} = 6$, dashed to $\mathfrak{k} = 2$ and dash-dot to $\mathfrak{k} = 11$. For $\mathfrak{k} = 2$ we have the values on the diagonal $||\alpha_0|| = 7.8e-08$, $||\beta_0|| = 3.8e-05$, and in the corner $||\alpha_{2^{n-1}}|| = 3.0e-02$, $||\beta_{2^{n-1}}|| = 9.8e-03$. For $\mathfrak{k} = 6$ we have the values on the diagonal $||\alpha_0|| = 5.3e-19$, $||\beta_0|| = 2.4e-12$, and in the corner $||\alpha_{2^{n-1}}|| = 1.6e-03$, $||\beta_{2^{n-1}}|| = 9.1e-03$. For $\mathfrak{k} = 11$ we have the values on the diagonal $||\alpha_0|| = 3.2e-19$, $||\beta_0|| = 1.5e-18$, and in the corner $||\alpha_{2^{n-1}}|| = 6.8e-09$, $||\beta_{2^{n-1}}|| = 1.9e-05$.</figcaption></figure>

where $\phi_m$ are the Legendre scaling functions. Here $x_0, \dots, x_{\mathfrak{k}-1}$ denote the roots of $P_{\mathfrak{k}}(2x-1)$. We are interested in the following combination

$$ \Psi_j(-x) \Psi_p(x) = \sqrt{w_j w_p} \sum_{j', p' = 0}^{\mathfrak{k}-1} \phi_{j'}(x_j) \phi_{p'}(x_p) \Phi_{j'}(-x) \Phi_{p'}(x). $$

Thus it can be presented as the series

$$ \Psi_j(-x) \Psi_p(x) = \sum_{k=0}^{\infty} \frac{D_{jp}^k x^k}{(k+2)!}, $$

where

$$ D_{jp}^k = \sqrt{w_j w_p} \sum_{\substack{0 \le j', p' \le \mathfrak{k}-1 \\ j'+p' \le k}} \phi_{j'}(x_j) \phi_{p'}(x_p) C_{j'p'}^{k-j'-p'}, \quad j, p = 0, \dots, \mathfrak{k}-1. $$

Note that this double sum can be restricted to $k - j' - p'$ being even non-negative. Thus our main operator in the interpolating basis has the following expansion

$$ [\sigma_l^n]_{pj}(a) = \sum_{k=0}^{\infty} D_{jp}^k J_k(l, a), \qquad (4.25) $$

where as above $a = tN^2 = t4^n$.

Calculating matrices $\alpha_l^n, \beta_l^n, \gamma_l^n$ from $\sigma_l^{n+1}$ by (3.7) and then evaluating the norms $||\alpha_l^n||, ||\beta_l^n||$ and $||\gamma_l^n||$, one arrives to qualitatively the same results as above, and so we do not illustrate them on separate figures.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Sparsity of the non-standard form at scale $n = 9$ corresponding to the time step $t = 0.0001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 7.8e-08$, $||\beta_0|| = 3.8e-05$, and in the corner $||\alpha_{2^{n-1}}|| = 3.0e-02$, $||\beta_{2^{n-1}}|| = 9.8e-03$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 5.3e-19$, $||\beta_0|| = 2.4e-12$, and in the corner $||\alpha_{2^{n-1}}|| = 1.6e-03$, $||\beta_{2^{n-1}}|| = 9.1e-03$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 3.2e-19$, $||\beta_0|| = 1.5e-18$, and in the corner $||\alpha_{2^{n-1}}|| = 6.8e-09$, $||\beta_{2^{n-1}}|| = 1.9e-05$.</figcaption></figure>

where $\phi_m$ are the Legendre scaling functions. Here $x_0, \dots, x_{k-1}$ denote the roots of $P_k(2x-1)$. We are interested in the following combination

$$ \Psi_j(-x) \Psi_p(x) = \sqrt{w_j w_p} \sum_{j', p' = 0}^{k-1} \phi_{j'}(x_j) \phi_{p'}(x_p) \Phi_{j'}(-x) \Phi_{p'}(x). $$

Thus it can be presented as the series

$$ \Psi_j(-x) \Psi_p(x) = \sum_{k=0}^{\infty} \frac{D_{jp}^k x^k}{(k+2)!}, $$

where

$$ D_{jp}^k = \sqrt{w_j w_p} \sum_{\substack{0 \le j', p' \le k-1 \\ j' + p' \le k}} \phi_{j'}(x_j) \phi_{p'}(x_p) C_{j'p'}^{k-j'-p'}, \quad j, p = 0, \dots, k-1. $$

Note that this double sum can be restricted to $k-j'-p'$ being even non-negative. Thus our main operator in the interpolating basis has the following expansion

$$ [\sigma_l^n]_{pj}(a) = \sum_{k=0}^{\infty} D_{jp}^k J_k(l, a), \qquad (4.25) $$

where as above $a = tN^2 = t4^n$.

Calculating matrices $\alpha_l^n, \beta_l^n, \gamma_l^n$ from $\sigma_l^{n+1}$ by (3.7) and then evaluating the norms $||\alpha_l^n||, ||\beta_l^n||$ and $||\gamma_l^n||$, one arrives to qualitatively the same results as above, and so we do not illustrate them on separate figures.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Sparsity of the non-standard form at scale $n = 5$ corresponding to the time step $t = 0.001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 2.6\text{e-04}$, $||\beta_0|| = 4.9\text{e-03}$, and in the corner $||\alpha_{2^{n-1}}|| = 2.6\text{e-02}$, $||\beta_{2^{n-1}}|| = 1.1\text{e-02}$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 6.8\text{e-12}$, $||\beta_0|| = 2.0\text{e-07}$, and in the corner $||\alpha_{2^{n-1}}|| = 1.5\text{e-01}$, $||\beta_{2^{n-1}}|| = 1.4\text{e-01}$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 2.4\text{e-18}$, $||\beta_0|| = 8.5\text{e-14}$, and in the corner $||\alpha_{2^{n-1}}|| = 1.5\text{e-04}$, $||\beta_{2^{n-1}}|| = 6.4\text{e-03}$.</figcaption></figure>

4.3. **Haar multiresolution analysis.** The above conclusions simplify significantly in the case $k = 1$. By (4.6), (3.3) we have

$$[\sigma_l^n]_{00} = \int_{\Gamma} F(\zeta)d\zeta, \quad F(\zeta) = \frac{1 - \cos \zeta}{\pi \zeta^2} \exp(i\zeta l - itN^2\zeta^2).$$

Let us first justify (4.1), so we regard $\zeta$ lying in the band $|\operatorname{Im}\zeta| \le H$ and satisfying $\operatorname{Re}\zeta \cdot \operatorname{Im}\zeta \le 0$, then

$$|F(\zeta)| \le \frac{2}{\pi} e^{H(|l|+1)} \le \frac{2}{\pi} e^{HN}.$$

In order to maintain accuracy $\varepsilon$, while discretizing the integral in such domain, the right hand of this bound should not be bigger than $\varepsilon/\varepsilon_{\text{mach}}$. Thus we obtained (4.1).

Now we have only even power integrals in the expansion (4.18),

$$[\sigma_l^n]_{00} = 2 \sum_{k=0}^{\infty} J_{2k}(l, a), \qquad (4.26)$$

since $C_{00}^k = 1 + (-1)^k$ by (4.17), that of course can be checked directly by expanding $\cos \zeta$ about zero. The convergence rate of this series was analysed in Subsection 4.1.

Finally, let us have a closer look at the dependence of $[\sigma_l^n]_{00}$ on the distance $|l|$ to the diagonal. In the case of $l = 0$, one can calculate the sum of the series by means of special functions as

$$[\sigma_0^n]_{00} = -i \sqrt{\frac{2}{\pi}} F_1 \left( \frac{e^{\frac{i\pi}{4}}}{\sqrt{2a}} \right), \qquad (4.27)$$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Sparsity of the non-standard form at scale $n = 5$ corresponding to the time step $t = 0.001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 2.6\text{e-04}$, $||\beta_0|| = 4.9\text{e-03}$, and in the corner $||\alpha_{2n-1}|| = 2.6\text{e-02}$, $||\beta_{2n-1}|| = 1.1\text{e-02}$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 6.8\text{e-12}$, $||\beta_0|| = 2.0\text{e-07}$, and in the corner $||\alpha_{2n-1}|| = 1.5\text{e-01}$, $||\beta_{2n-1}|| = 1.4\text{e-01}$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 2.4\text{e-18}$, $||\beta_0|| = 8.5\text{e-14}$, and in the corner $||\alpha_{2n-1}|| = 1.5\text{e-04}$, $||\beta_{2n-1}|| = 6.4\text{e-03}$.</figcaption></figure>

4.3. **Haar multiresolution analysis.** The above conclusions simplify significantly in the case $k = 1$. By (4.6), (3.3) we have

$$[\sigma_l^n]_{00} = \int_{\Gamma} F(\zeta) d\zeta, \quad F(\zeta) = \frac{1 - \cos \zeta}{\pi \zeta^2} \exp(i\zeta l - itN^2\zeta^2).$$

Let us first justify (4.1), so we regard $\zeta$ lying in the band $|\operatorname{Im}\zeta| \le H$ and satisfying $\operatorname{Re}\zeta \cdot \operatorname{Im}\zeta \le 0$, then

$$|F(\zeta)| \le \frac{2}{\pi} e^{H(|l|+1)} \le \frac{2}{\pi} e^{HN}.$$

In order to maintain accuracy $\varepsilon$, while discretizing the integral in such domain, the right hand of this bound should not be bigger than $\varepsilon/\varepsilon_{\text{mach}}$. Thus we obtained (4.1).

Now we have only even power integrals in the expansion (4.18),

$$[\sigma_l^n]_{00} = 2 \sum_{k=0}^{\infty} J_{2k}(l, a), \qquad (4.26)$$

since $C_{00}^k = 1 + (-1)^k$ by (4.17), that of course can be checked directly by expanding $\cos \zeta$ about zero. The convergence rate of this series was analysed in Subsection 4.1.

Finally, let us have a closer look at the dependence of $[\sigma_l^n]_{00}$ on the distance $|l|$ to the diagonal. In the case of $l = 0$, one can calculate the sum of the series by means of special functions as

$$[\sigma_0^n]_{00} = -i \sqrt{\frac{2}{\pi}} F_1 \left( \frac{e^{\frac{i\pi}{4}}}{\sqrt{2a}} \right), \qquad (4.27)$$
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Sparsity of the non-standard form at scale $n = 6$ corresponding to the time step $t = 0.001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 8.1e-06$, $||\beta_0|| = 6.1e-04$, and in the corner $||\alpha_{2^{n-1}}|| = 9.8e-02$, $||\beta_{2^{n-1}}|| = 3.7e-02$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 8.3e-16$, $||\beta_0|| = 1.6e-09$, and in the corner $||\alpha_{2^{n-1}}|| = 4.1e-04$, $||\beta_{2^{n-1}}|| = 7.5e-03$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 6.3e-19$, $||\beta_0|| = 1.1e-17$, and in the corner $||\alpha_{2^{n-1}}|| = 1.4e-10$, $||\beta_{2^{n-1}}|| = 4.3e-06$.</figcaption></figure>

where

$$F_1(x) = \sum_{k=0}^{\infty} \frac{x^{2k+1}}{(2k+1)(2k+2)!!} = \frac{1-e^{x^2/2}}{x} - i\sqrt{\frac{\pi}{2}} \operatorname{erf}\left(\frac{ix}{\sqrt{2}}\right),$$

using (4.23). The series $F_1\left(\frac{e^{i\pi}}{\sqrt{2a}}\right)$ is very useful in practice for the precision control, since its real and imaginary parts are alternating series with monotonically decreasing coefficients.

In the case of $|l| \gg 1$, we have (4.24) which leads to

$$[\sigma_l^n]_{00} \approx \frac{1}{\sqrt{\pi a}} \left(\frac{2a}{l}\right)^2 \left(1 - \cos\frac{l}{2a}\right) \exp\left(\frac{il^2}{4a} - \frac{i\pi}{4}\right). \quad (4.28)$$

Here the term in front of the exponent changes slowly with $l$, whereas the exponent itself oscillates very fast. This formula describes the operator matrix with a very high precision for any $l \neq 0$. In fact the limit of this expression as $l/2a \to 0$ is identical to $[\sigma_0^n]_{00}$ evaluated by (4.27) with the first order approximation $F_1(x) \approx x/2$. Moreover, (4.27) and (4.28) give us an idea of how many terms are important in Expansion (4.26) for a given precision. The same series cut is reliable in practise for evaluation of the general expansion (4.18).

## 5. INTRODUCTION INTO EXPONENTIAL INTEGRATORS

It is known that a fundamental solution of the general linear equation

$$\partial_t u = Au + Bu \quad (5.1)$$

forms a two-parameter family of operators $\mathcal{T}(t_2, t_1)$ called propagator. It means that for any solution $u = u(t)$ of (5.1) and time moments $t_1, t_2$ we have $u(t_2) = \mathcal{T}(t_2, t_1)u(t_1)$. This section serves as a
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Sparsity of the non-standard form at scale $n = 6$ corresponding to the time step $t = 0.001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 8.1e-06$, $||\beta_0|| = 6.1e-04$, and in the corner $||\alpha_{2^{n-1}}|| = 9.8e-02$, $||\beta_{2^{n-1}}|| = 3.7e-02$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 8.3e-16$, $||\beta_0|| = 1.6e-09$, and in the corner $||\alpha_{2^{n-1}}|| = 4.1e-04$, $||\beta_{2^{n-1}}|| = 7.5e-03$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 6.3e-19$, $||\beta_0|| = 1.1e-17$, and in the corner $||\alpha_{2^{n-1}}|| = 1.4e-10$, $||\beta_{2^{n-1}}|| = 4.3e-06$.</figcaption></figure>

where

$$F_1(x) = \sum_{k=0}^{\infty} \frac{x^{2k+1}}{(2k+1)(2k+2)!!} = \frac{1-e^{x^2/2}}{x} - i\sqrt{\frac{\pi}{2}} \operatorname{erf}\left(\frac{ix}{\sqrt{2}}\right),$$

using (4.23). The series $F_1\left(\frac{e^{i\pi}}{\sqrt{2a}}\right)$ is very useful in practice for the precision control, since its real and imaginary parts are alternating series with monotonically decreasing coefficients.

In the case of $|l| \gg 1$, we have (4.24) which leads to

$$[\sigma_l^n]_{00} \approx \frac{1}{\sqrt{\pi a}} \left(\frac{2a}{l}\right)^2 \left(1 - \cos\frac{l}{2a}\right) \exp\left(\frac{il^2}{4a} - \frac{i\pi}{4}\right). \quad (4.28)$$

Here the term in front of the exponent changes slowly with $l$, whereas the exponent itself oscillates very fast. This formula describes the operator matrix with a very high precision for any $l \neq 0$. In fact the limit of this expression as $l/2a \to 0$ is identical to $[\sigma_0^n]_{00}$ evaluated by (4.27) with the first order approximation $F_1(x) \approx x/2$. Moreover, (4.27) and (4.28) give us an idea of how many terms are important in Expansion (4.26) for a given precision. The same series cut is reliable in practise for evaluation of the general expansion (4.18).

## 5. INTRODUCTION INTO EXPONENTIAL INTEGRATORS

It is known that a fundamental solution of the general linear equation

$$\partial_t u = Au + Bu \quad (5.1)$$

forms a two-parameter family of operators $\mathcal{T}(t_2, t_1)$ called propagator. It means that for any solution $u = u(t)$ of (5.1) and time moments $t_1, t_2$ we have $u(t_2) = \mathcal{T}(t_2, t_1)u(t_1)$. This section serves as a
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Sparsity of the non-standard form at scale $n = 7$ corresponding to the time step $t = 0.001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 2.5e-07$, $||\beta_0|| = 7.7e-05$, and in the corner $||\alpha_{2^{n-1}}|| = 1.3e-02$, $||\beta_{2^{n-1}}|| = 2.7e-02$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 4.2e-19$, $||\beta_0|| = 1.2e-11$, and in the corner $||\alpha_{2^{n-1}}|| = 1.1e-07$, $||\beta_{2^{n-1}}|| = 8.8e-05$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 2.0e-19$, $||\beta_0|| = 1.8e-18$, and in the corner $||\alpha_{2^{n-1}}|| = 3.0e-17$, $||\beta_{2^{n-1}}|| = 1.4e-09$.</figcaption></figure>

brief reminder of how this propagator can be approximated by the exponential operators $e^{At}, e^{Bt}$, providing us with effective numerical schemes for (5.1). Up to the imaginary unit, $A$ and $B$ will stand for kinetic and potential operators below, as in Equation (1.1), for example.

5.1. **Time independent potential.** Let $A, B$ be constant. Then the propagator $\mathcal{T}$ reduces to the one-parameter family $\mathcal{T}(t_2, t_1) = \mathcal{T}(t_2-t_1)$ with $\mathcal{T}(t) = e^{(A+B)t}$ called a semigroup. Since in general operators $A$ and $B$ do not commute, one cannot substitute this exponent with the multiplication of semigroups associated with $A$ and $B$ each. On the contrary, such substitution is valid for small $t$, namely,

$$e^{At+Bt} = e^{At}e^{Bt} + \mathcal{O}(t^2). \qquad (5.2)$$

A subsequent application of the multiplication $e^{At}e^{Bt}$ to a given initial data $u(0)$ would provide us with the solution at a given final time moment with first order $\mathcal{O}(t)$ of accuracy. A more accurate result one can get with the symmetric symplectic scheme (1.4). General higher order factorizations of this type were systematically investigated in [34, 42]. For any given positive integer $n$ one can decompose the semigroup $\mathcal{T}(t)$ as

$$e^{At+Bt} = e^{At_1}e^{Bt_2}e^{At_3}e^{Bt_4}\dots e^{At_M} + \mathcal{O}(t^{n+1}) \qquad (5.3)$$

with some appropriately chosen parameters $t_1, \dots, t_M$. It allows to set bigger time steps and accelerate computations. However, it turns out that $M$ increases with $n$ exponentially, and so in practice the integrator order $n$ is chosen no bigger than 8, since too many exponents in the decomposition formula can on the opposite slow down the calculations. Moreover, the time parameters $t_1, \dots, t_M$ are not unique. Below we try to choose those ones that lead to exponents $\exp(it_j\partial_x^2)$ in the decomposition (5.3) with overall biggest possible absolute values $|t_j|$ to incorporate the multiresolution sparsity at most.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Sparsity of the non-standard form at scale $n = 7$ corresponding to the time step $t = 0.001$. Here $||\alpha_l||$ is blue, $||\beta_l||$ is red, $||\sigma_l||$ is black. Solid line corresponds to $k = 6$, dashed to $k = 2$ and dash-dot to $k = 11$. For $k = 2$ we have the values on the diagonal $||\alpha_0|| = 2.5e-07$, $||\beta_0|| = 7.7e-05$, and in the corner $||\alpha_{2^{n-1}}|| = 1.3e-02$, $||\beta_{2^{n-1}}|| = 2.7e-02$. For $k = 6$ we have the values on the diagonal $||\alpha_0|| = 4.2e-19$, $||\beta_0|| = 1.2e-11$, and in the corner $||\alpha_{2^{n-1}}|| = 1.1e-07$, $||\beta_{2^{n-1}}|| = 8.8e-05$. For $k = 11$ we have the values on the diagonal $||\alpha_0|| = 2.0e-19$, $||\beta_0|| = 1.8e-18$, and in the corner $||\alpha_{2^{n-1}}|| = 3.0e-17$, $||\beta_{2^{n-1}}|| = 1.4e-09$.</figcaption></figure>

brief reminder of how this propagator can be approximated by the exponential operators $e^{At}, e^{Bt}$, providing us with effective numerical schemes for (5.1). Up to the imaginary unit, $A$ and $B$ will stand for kinetic and potential operators below, as in Equation (1.1), for example.

5.1. **Time independent potential.** Let $A, B$ be constant. Then the propagator $\mathcal{T}$ reduces to the one-parameter family $\mathcal{T}(t_2, t_1) = \mathcal{T}(t_2-t_1)$ with $\mathcal{T}(t) = e^{(A+B)t}$ called a semigroup. Since in general operators $A$ and $B$ do not commute, one cannot substitute this exponent with the multiplication of semigroups associated with $A$ and $B$ each. On the contrary, such substitution is valid for small $t$, namely,

$$e^{At+Bt} = e^{At}e^{Bt} + \mathcal{O}(t^2). \qquad (5.2)$$

A subsequent application of the multiplication $e^{At}e^{Bt}$ to a given initial data $u(0)$ would provide us with the solution at a given final time moment with first order $\mathcal{O}(t)$ of accuracy. A more accurate result one can get with the symmetric symplectic scheme (1.4). General higher order factorizations of this type were systematically investigated in [34, 42]. For any given positive integer $n$ one can decompose the semigroup $\mathcal{T}(t)$ as

$$e^{At+Bt} = e^{At_1}e^{Bt_2}e^{At_3}e^{Bt_4}\dots e^{At_M} + \mathcal{O}(t^{n+1}) \qquad (5.3)$$

with some appropriately chosen parameters $t_1, \dots, t_M$. It allows to set bigger time steps and accelerate computations. However, it turns out that $M$ increases with $n$ exponentially, and so in practice the integrator order $n$ is chosen no bigger than 8, since too many exponents in the decomposition formula can on the opposite slow down the calculations. Moreover, the time parameters $t_1, \dots, t_M$ are not unique. Below we try to choose those ones that lead to exponents $\exp(it_j\partial_x^2)$ in the decomposition (5.3) with overall biggest possible absolute values $|t_j|$ to incorporate the multiresolution sparsity at most.
So far we have not employed the fact that A stands for the kinetic and B for the potential energy operators. An advantage of this particular choice was taken in [14, 15, 36], where several fourth order schemes were introduced allowing minimal amount of exponents in the decomposition (5.3). Among them we extensively use the following

$$e^{At+Bt} = \exp\left(\frac{t}{6}B\right) \exp\left(\frac{t}{2}A\right) \exp\left(\frac{2t}{3}\tilde{B}\right) \exp\left(\frac{t}{2}A\right) \exp\left(\frac{t}{6}B\right) + \mathcal{O}(t^5), \quad (5.4)$$

where

$$\tilde{B} = \tilde{B}(t) = B + \frac{t^2}{48}[B, [A, B]].$$

This double commutator is proportional to the squared gradient of the potential $(\partial_x V)^2$, which makes $\tilde{B}$ a potential operator when $B = -iV$. Indeed,

$$\tilde{B} = -i\tilde{V}, \quad \tilde{V} = \begin{cases} V - \frac{t^2}{24}(\partial_x V)^2, & A = i\partial_x^2, \quad B = -iV \\ V - \frac{t^2}{48}(\partial_x V)^2, & A = \frac{i}{2}\partial_x^2, \quad B = -iV \end{cases},$$

where we regarded two of the most used normalisations for the kinetic energy operator. The fourth order scheme (5.4) requires only two applications of the free-particle semigroup operator. On the contrary it demands to know the gradient of the potential. For comparison we mention an alternative scheme from [15] reading

$$e^{At+Bt} = \exp\left(\frac{t}{8}\tilde{B}\right) \exp\left(\frac{t}{3}A\right) \exp\left(\frac{3t}{8}B\right) \exp\left(\frac{t}{3}A\right) \exp\left(\frac{3t}{8}B\right) \exp\left(\frac{t}{3}A\right) \exp\left(\frac{t}{8}\tilde{B}\right) + \mathcal{O}(t^5), \quad (5.5)$$

where $\tilde{B}$ is the same. Apart from the fact that the free-particle operator $\exp(\frac{t}{3}A)$ appears here 3 times, the factor $1/3$ in front of $tA$ makes its multiresolution less sparse comparing to $\exp(\frac{t}{2}A)$ standing in (5.4). Scheme (5.5) could potentially be of use only together with (5.4) for precision control. The operator splitting (5.4) is far superior to any other existing fourth order symplectic algorithm, to our knowledge.

These algorithms can be recursively extended to higher order schemes [15]. Denoting by $\mathcal{T}$ the evolution operator $\mathcal{T}(t) = e^{At+Bt}$ and by $\mathcal{T}^{(n)}$ its $n$-th order approximation, that is $\mathcal{T}(t) = \mathcal{T}^{(n)}(t) + \mathcal{O}(t^{n+1})$, we can express an $n+2$ order scheme as

$$\mathcal{T}^{(n+2)}(t) = \mathcal{T}^{(n)}(t_n)\mathcal{T}^{(n)}(-s_n t_n)\mathcal{T}^{(n)}(t_n), \quad t_n = \frac{t}{2-s_n}, \quad s_n = 2^{1/(n+1)}. \quad (5.6)$$

A proof can be found in [42]. We make a couple of remarks concerning this recursive formula. Firstly, $\mathcal{T}(-t) = \mathcal{T}^{-1}(t)$ that is in connection to the Schrödinger equation coincides with the adjoint $\mathcal{T}^*(t)$. Secondly, $s_n t_n > t_n > t$, which provides us with a little advantage of exploiting (5.6), since the precision and sparsity of the free particle propagator multiresolution is increasing together with the time step. Although there is an obvious drawback behind (5.6), as the amount of work is increasing three times per time step whenever one goes for a higher order scheme. As a working example of a sixth order scheme, we iterate $\mathcal{T}^{(4)}$ given by (5.4) using (5.6) to get

$$\begin{aligned} \mathcal{T}^{(6)}(t) = & \exp\left(\frac{t_4}{6}B\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{2t_4}{3}\tilde{B}(t_4)\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{1-s_4}{6}B\right) \\ & \exp\left(-\frac{s_4t_4}{2}A\right) \exp\left(-\frac{2s_4t_4}{3}\tilde{B}(s_4t_4)\right) \exp\left(-\frac{s_4t_4}{2}A\right) \\ & \exp\left(\frac{(1-s_4)t_4}{6}B\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{2t_4}{3}\tilde{B}(t_4)\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{t_4}{6}B\right). \end{aligned} \quad (5.7)$$

Here the last 6 operators are the first 6 operators following in reverse, so that the whole operator $\mathcal{T}^{(6)}(t)$ is symmetric. Note that this scheme demands 6 applications of kinetic energy exponent per time step. If we continue the recursive procedure (5.6) we will arrive to an eighth order scheme
So far we have not employed the fact that A stands for the kinetic and B for the potential energy operators. An advantage of this particular choice was taken in [14, 15, 36], where several fourth order schemes were introduced allowing minimal amount of exponents in the decomposition (5.3). Among them we extensively use the following

$$e^{At+Bt} = \exp\left(\frac{t}{6}B\right) \exp\left(\frac{t}{2}A\right) \exp\left(\frac{2t}{3}\tilde{B}\right) \exp\left(\frac{t}{2}A\right) \exp\left(\frac{t}{6}B\right) + \mathcal{O}(t^5), \quad (5.4)$$

where

$$\tilde{B} = \tilde{B}(t) = B + \frac{t^2}{48}[B, [A, B]].$$

This double commutator is proportional to the squared gradient of the potential $(\partial_x V)^2$, which makes $\tilde{B}$ a potential operator when $B = -iV$. Indeed,

$$\tilde{B} = -i\tilde{V}, \quad \tilde{V} = \begin{cases} V - \frac{t^2}{24}(\partial_x V)^2, & A = i\partial_x^2, \quad B = -iV \\ V - \frac{t^2}{48}(\partial_x V)^2, & A = \frac{i}{2}\partial_x^2, \quad B = -iV \end{cases},$$

where we regarded two of the most used normalisations for the kinetic energy operator. The fourth order scheme (5.4) requires only two applications of the free-particle semigroup operator. On the contrary it demands to know the gradient of the potential. For comparison we mention an alternative scheme from [15] reading

$$e^{At+Bt} = \exp\left(\frac{t}{8}\tilde{B}\right) \exp\left(\frac{t}{3}A\right) \exp\left(\frac{3t}{8}B\right) \exp\left(\frac{t}{3}A\right) \exp\left(\frac{3t}{8}B\right) \exp\left(\frac{t}{3}A\right) \exp\left(\frac{t}{8}\tilde{B}\right) + \mathcal{O}(t^5), \quad (5.5)$$

where $\tilde{B}$ is the same. Apart from the fact that the free-particle operator $\exp(\frac{t}{3}A)$ appears here 3 times, the factor $1/3$ in front of $tA$ makes its multiresolution less sparse comparing to $\exp(\frac{t}{2}A)$ standing in (5.4). Scheme (5.5) could potentially be of use only together with (5.4) for precision control. The operator splitting (5.4) is far superior to any other existing fourth order symplectic algorithm, to our knowledge.

These algorithms can be recursively extended to higher order schemes [15]. Denoting by $\mathcal{T}$ the evolution operator $\mathcal{T}(t) = e^{At+Bt}$ and by $\mathcal{T}^{(n)}$ its $n$-th order approximation, that is $\mathcal{T}(t) = \mathcal{T}^{(n)}(t) + \mathcal{O}(t^{n+1})$, we can express an $n+2$ order scheme as

$$\mathcal{T}^{(n+2)}(t) = \mathcal{T}^{(n)}(t_n)\mathcal{T}^{(n)}(-s_n t_n)\mathcal{T}^{(n)}(t_n), \quad t_n = \frac{t}{2-s_n}, \quad s_n = 2^{1/(n+1)}. \quad (5.6)$$

A proof can be found in [42]. We make a couple of remarks concerning this recursive formula. Firstly, $\mathcal{T}(-t) = \mathcal{T}^{-1}(t)$ that is in connection to the Schrödinger equation coincides with the adjoint $\mathcal{T}^*(t)$. Secondly, $s_n t_n > t_n > t$, which provides us with a little advantage of exploiting (5.6), since the precision and sparsity of the free particle propagator multiresolution is increasing together with the time step. Although there is an obvious drawback behind (5.6), as the amount of work is increasing three times per time step whenever one goes for a higher order scheme. As a working example of a sixth order scheme, we iterate $\mathcal{T}^{(4)}$ given by (5.4) using (5.6) to get

$$\begin{aligned} \mathcal{T}^{(6)}(t) = & \exp\left(\frac{t_4}{6}B\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{2t_4}{3}\tilde{B}(t_4)\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{(1-s_4)t_4}{6}B\right) \\ & \exp\left(-\frac{s_4t_4}{2}A\right) \exp\left(-\frac{2s_4t_4}{3}\tilde{B}(s_4t_4)\right) \exp\left(-\frac{s_4t_4}{2}A\right) \\ & \exp\left(\frac{(1-s_4)t_4}{6}B\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{2t_4}{3}\tilde{B}(t_4)\right) \exp\left(\frac{t_4}{2}A\right) \exp\left(\frac{t_4}{6}B\right). \end{aligned} \quad (5.7)$$

Here the last 6 operators are the first 6 operators following in reverse, so that the whole operator $\mathcal{T}^{(6)}(t)$ is symmetric. Note that this scheme demands 6 applications of kinetic energy exponent per time step. If we continue the recursive procedure (5.6) we will arrive to an eighth order scheme
demanding 18 applications, correspondingly. It is still a reasonable amount. However, at this level there is a slightly better alternative derived by Yoshida [42] that we describe next. It is clear from the iterative procedure (5.6) that an even order symmetric integrator can have the form

$$ \mathcal{T}^{(n)}(t) = \mathcal{T}^{(2)}(w_m t) \dots \mathcal{T}^{(2)}(w_0 t) \dots \mathcal{T}^{(2)}(w_m t) $$

with $2m + 1$ second order integrators $\mathcal{T}^{(2)}(t)$ defined by the splitting (1.4) and some weights $w_0, w_1, \dots, w_m$. Number $m$ and these weights depend on the order $n$. However, a direct implementation of (5.6) does not necessarily give the most optimal representation. It turns out that $m$ can be diminished down to $m = 3$ for the order $n = 6$ and to $m = 7$ for the order $n = 8$. The weights are not unique, we choose those providing us with the biggest value $\min_{j=0,\dots,m} |w_j|$, namely, set A from Table 1 and set C from Table 2 in [42]. They are repeated below for convenience.

Case $n = 6$:

$$ \begin{aligned} w_0 &= 1.315186320683906, & w_1 &= -1.17767998417887, \\ w_2 &= 0.235573213359357, & w_3 &= 0.784513610477560. \end{aligned} $$

Case $n = 8$:

$$ \begin{aligned} w_0 &= 1.65899088454396, & w_1 &= 0.311790812418427, \\ w_2 &= -1.55946803821447, & w_3 &= -1.67896928259640, \\ w_4 &= 1.66335809963315, & w_5 &= -1.06458714789183, \\ w_6 &= 1.36934946416871, & w_7 &= 0.629030650210433. \end{aligned} $$

Thus both 6th and 8th order Yoshida integrators can be combined together in the product

$$ \mathcal{T}^{(n)}(t) = \exp(c_1 t B) \exp(d_1 t A) \exp(c_2 t B) \dots \exp(d_{2m+1} t A) \exp(c_{2m+2} t B) \quad (5.8) $$

where $t$ is the time step and $c_i, d_i$ are constants given by

$$ c_1 = c_{2m+2} = \frac{w_m}{2}, \quad c_2 = c_{2m+1} = \frac{w_m + w_{m-1}}{2}, \quad \dots, \quad c_{m+1} = c_{m+2} = \frac{w_1 + w_0}{2} $$

and

$$ d_1 = d_{2m+1} = w_m, \quad d_2 = d_{2m} = w_{m-1}, \quad \dots, \quad d_m = d_{m+2} = w_1, \quad d_{m+1} = w_0. $$

One can notice that the integrator is symmetric. Moreover, it is quite general, one does not need to know the potential gradient $\partial_x V(x)$. The 8th order scheme demands 15 kinetic exponent applications, less than 18 corresponding to the twice iterated (5.4) by (5.6). The weights associated with the latter lie in the interval $[0.65, 0.84]$ in absolute value, as one can easily deduce from (5.4), (5.6). Thus in the case of stationary potential $V(x)$ one may anticipate the best performance for

order $n = 4$: Scheme (5.4),

order $n = 6$: Scheme (5.7),

order $n = 8$: Scheme (5.8).

The difference between Scheme (5.7) and Scheme (5.8) is mild, though in [15] it is claimed that (5.7) outperforms slightly (5.8) at least for the scattering experiment regarded there.

5.2. **Time dependent potential.** Now let $B = B(t)$, the propagator $\mathcal{T} = \mathcal{T}(t + \tau, t)$ can still be approximated at a small time step $\tau = \Delta t$ by exponential operators. An $n$-th order integrator propagating solutions of (5.1) through the interval $[t, t + \tau]$ is denoted by $\mathcal{T}^{(n)}(t + \tau, t)$, meaning $\mathcal{T}(t + \tau, t) = \mathcal{T}^{(n)}(t + \tau, t) + \mathcal{O}(\tau^{n+1})$. Using Suzuki's formal approach [35] we extend the schemes from the previous subsection to the current situation. Let $\mathcal{D}$ be the forward time derivative operator, called also the super-operator in [35], that is a time derivative acting on functions standing on the left, namely,

$$ F(t)e^{\tau \mathcal{D}} G(t) = F(t + \tau)G(t) \quad (5.9) $$

for any time-dependent functions $F(t)$ and $G(t)$. In particular, $e^{\tau \mathcal{D}} G(t) = G(t)$. With the help of this super-operator the propagator is reduced to the exponential representation

$$ \mathcal{T}(t + \tau, t) = \exp(\tau(A + B(t) + \mathcal{D})). \quad (5.10) $$
demanding 18 applications, correspondingly. It is still a reasonable amount. However, at this level there is a slightly better alternative derived by Yoshida [42] that we describe next. It is clear from the iterative procedure (5.6) that an even order symmetric integrator can have the form

$$ \mathcal{T}^{(n)}(t) = \mathcal{T}^{(2)}(w_m t) \dots \mathcal{T}^{(2)}(w_0 t) \dots \mathcal{T}^{(2)}(w_m t) $$

with $2m + 1$ second order integrators $\mathcal{T}^{(2)}(t)$ defined by the splitting (1.4) and some weights $w_0, w_1, \dots, w_m$. Number $m$ and these weights depend on the order $n$. However, a direct implementation of (5.6) does not necessarily give the most optimal representation. It turns out that $m$ can be diminished down to $m = 3$ for the order $n = 6$ and to $m = 7$ for the order $n = 8$. The weights are not unique, we choose those providing us with the biggest value $\min_{j=0,\dots,m} |w_j|$, namely, set A from Table 1 and set C from Table 2 in [42]. They are repeated below for convenience.

Case $n = 6$:

$$ \begin{aligned} w_0 &= 1.315186320683906, & w_1 &= -1.17767998417887, \\ w_2 &= 0.235573213359357, & w_3 &= 0.784513610477560. \end{aligned} $$

Case $n = 8$:

$$ \begin{aligned} w_0 &= 1.65899088454396, & w_1 &= 0.311790812418427, \\ w_2 &= -1.55946803821447, & w_3 &= -1.67896928259640, \\ w_4 &= 1.66335809963315, & w_5 &= -1.06458714789183, \\ w_6 &= 1.36934946416871, & w_7 &= 0.629030650210433. \end{aligned} $$

Thus both 6th and 8th order Yoshida integrators can be combined together in the product

$$ \mathcal{T}^{(n)}(t) = \exp(c_1 t B) \exp(d_1 t A) \exp(c_2 t B) \dots \exp(d_{2m+1} t A) \exp(c_{2m+2} t B) \quad (5.8) $$

where $t$ is the time step and $c_i, d_i$ are constants given by

$$ c_1 = c_{2m+2} = \frac{w_m}{2}, \quad c_2 = c_{2m+1} = \frac{w_m + w_{m-1}}{2}, \quad \dots, \quad c_{m+1} = c_{m+2} = \frac{w_1 + w_0}{2} $$

and

$$ d_1 = d_{2m+1} = w_m, \quad d_2 = d_{2m} = w_{m-1}, \quad \dots, \quad d_m = d_{m+2} = w_1, \quad d_{m+1} = w_0. $$

One can notice that the integrator is symmetric. Moreover, it is quite general, one does not need to know the potential gradient $\partial_x V(x)$. The 8th order scheme demands 15 kinetic exponent applications, less than 18 corresponding to the twice iterated (5.4) by (5.6). The weights associated with the latter lie in the interval $[0.65, 0.84]$ in absolute value, as one can easily deduce from (5.4), (5.6). Thus in the case of stationary potential $V(x)$ one may anticipate the best performance for

order $n = 4$: Scheme (5.4),

order $n = 6$: Scheme (5.7),

order $n = 8$: Scheme (5.8).

The difference between Scheme (5.7) and Scheme (5.8) is mild, though in [15] it is claimed that (5.7) outperforms slightly (5.8) at least for the scattering experiment regarded there.

5.2. **Time dependent potential.** Now let $B = B(t)$, the propagator $\mathcal{T} = \mathcal{T}(t + \tau, t)$ can still be approximated at a small time step $\tau = \Delta t$ by exponential operators. An $n$-th order integrator propagating solutions of (5.1) through the interval $[t, t + \tau]$ is denoted by $\mathcal{T}^{(n)}(t + \tau, t)$, meaning $\mathcal{T}(t + \tau, t) = \mathcal{T}^{(n)}(t + \tau, t) + \mathcal{O}(\tau^{n+1})$. Using Suzuki's formal approach [35] we extend the schemes from the previous subsection to the current situation. Let $\mathcal{D}$ be the forward time derivative operator, called also the super-operator in [35], that is a time derivative acting on functions standing on the left, namely,

$$ F(t)e^{\tau \mathcal{D}} G(t) = F(t + \tau)G(t) \quad (5.9) $$

for any time-dependent functions $F(t)$ and $G(t)$. In particular, $e^{\tau \mathcal{D}} G(t) = G(t)$. With the help of this super-operator the propagator is reduced to the exponential representation

$$ \mathcal{T}(t + \tau, t) = \exp(\tau(A + B(t) + \mathcal{D})). \quad (5.10) $$
Combining $A, \mathcal{D}$ together and approximating this exponent exactly as was done in the previous subsection with $A+\mathcal{D}$ standing now in place of $A$ we can extend all the above schemes. Note that $A$ does not depend on time, so $A$ and $\mathcal{D}$ commute. Therefore, $\exp(\tau(A+\mathcal{D})) = \exp(\tau A)) \exp(\tau \mathcal{D})$. Thus approximating the exponent in (5.10) by the second order scheme (1.4) one deduces

$$ \mathcal{T}(t + \tau, t) = \exp\left(\frac{\tau}{2}B(t)\right) \exp(\tau(A + \mathcal{D})) \exp\left(\frac{\tau}{2}B(t)\right) + \mathcal{O}(t^3) $$

which together with the super-operator property (5.9) leads to the following simple second order splitting

$$ \mathcal{T}^{(2)}(t + \tau, t) = \exp\left(\frac{\tau}{2}B(t + \tau)\right) \exp(\tau A) \exp\left(\frac{\tau}{2}B(t)\right). \quad (5.11) $$

In the same manner we extend both 6th and 8th order Yoshida integrators to

$$ \begin{aligned} \mathcal{T}^{(n)}(t + \tau, t) = & \exp(c_1 \tau B(t + b_1 \tau)) \exp(d_1 \tau A) \exp(c_2 \tau B(t + b_2 \tau)) \\ & \dots \exp(d_{2m+1} \tau A) \exp(c_{2m+2} \tau B(t + b_{2m+2} \tau)), \quad (5.12) \end{aligned} $$

where $b_k$ are constants given by

$$ b_k = \sum_{j=k}^{2m+1} d_j = 1 - \sum_{j=1}^{k-1} d_j, $$

in particular, $b_1 = 1$ and $b_{2m+2} = 0$.

In the case of time dependent potential there is an analogue of (5.4) in [16] reading

$$ \mathcal{T}^{(4)}(t+\tau,t) = \exp\left(\frac{\tau}{6}B(t+\tau)\right) \exp\left(\frac{\tau}{2}A\right) \exp\left(\frac{2\tau}{3}\tilde{B}\left(t+\frac{\tau}{2},\tau\right)\right) \exp\left(\frac{\tau}{2}A\right) \exp\left(\frac{\tau}{6}B(t)\right), \quad (5.13) $$

where the modified potential given at the middle step time moment is defined as

$$ \tilde{B}(t, \tau) = B(t) + \frac{\tau^2}{48} [B(t), [A, B(t)]]. \quad (5.14) $$

This follows from (5.10), (5.4), (5.9) and the identity $[B(t), [\mathcal{D}, B(t)]] = 0$. The latter can be checked directly

$$ F(t)[B(t), [\mathcal{D}, B(t)]] = 2F(t)B(t)\mathcal{D}B(t) - F(t)B^2(t)\mathcal{D} - F(t)\mathcal{D}B^2(t) = 0 $$

on any smooth test operator $F(t)$. Scheme (5.13) is very useful, since the commutator in (5.14) turns out to be proportional to $\partial_x V(x, t)$ for the Schrödinger equation, compare also with the same scheme obtained recently for the Dirac equation [41]. Similarly, the sixth order integrator (5.7) turns into

$$ \begin{aligned} \mathcal{T}^{(6)}(t + \tau, t) = & \exp\left(\frac{\tau_4}{6}B(t + (2 - s_4)\tau_4)\right) \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{2\tau_4}{3}\tilde{B}\left(t + \left(\frac{3}{2} - s_4\right)\tau_4, \tau_4\right)\right) \\ & \cdot \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{1 - s_4}{6}\tau_4\right) B(t + (1 - s_4)\tau_4) \exp\left(-\frac{s_4\tau_4}{2}A\right) \\ & \cdot \exp\left(-\frac{2s_4\tau_4}{3}\tilde{B}\left(t + \left(1 - \frac{s_4}{2}\right)\tau_4, s_4\tau_4\right)\right) \exp\left(-\frac{s_4\tau_4}{2}A\right) \exp\left(\frac{1 - s_4}{6}\tau_4\right) B(t + \tau_4) \\ & \cdot \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{2\tau_4}{3}\tilde{B}\left(t + \frac{\tau_4}{2}, \tau_4\right)\right) \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{\tau_4}{6}B(t)\right). \quad (5.15) \end{aligned} $$

Note that we first extended the time-independent potential scheme (5.4) to (5.7) via recursion (5.6), afterwards we applied the Suzuki decomposition method (5.10) to deduce (5.15) from (5.7). Integrator (5.13) can also be directly extended to higher order schemes by the general Suzuki's recursive procedure

$$ \begin{aligned} \mathcal{T}^{(n+2)}(t + \tau, t) = & \mathcal{T}^{(n)}(t + \tau, t + (1 - \kappa_n)\tau) \mathcal{T}^{(n)}(t + (1 - \kappa_n)\tau, t + (1 - 2\kappa_n)\tau) \\ & \cdot \mathcal{T}^{(n)}(t + (1 - 2\kappa_n)\tau, t + 2\kappa_n\tau) \mathcal{T}^{(n)}(t + 2\kappa_n\tau, t + \kappa_n\tau) \mathcal{T}^{(n)}(t + \kappa_n\tau, t) \end{aligned} \quad (5.16) $$
Combining $A, \mathcal{D}$ together and approximating this exponent exactly as was done in the previous subsection with $A+\mathcal{D}$ standing now in place of $A$ we can extend all the above schemes. Note that $A$ does not depend on time, so $A$ and $\mathcal{D}$ commute. Therefore, $\exp(\tau(A+\mathcal{D})) = \exp(\tau A)) \exp(\tau \mathcal{D})$. Thus approximating the exponent in (5.10) by the second order scheme (1.4) one deduces

$$ \mathcal{T}(t + \tau, t) = \exp\left(\frac{\tau}{2}B(t)\right) \exp(\tau(A + \mathcal{D})) \exp\left(\frac{\tau}{2}B(t)\right) + \mathcal{O}(t^3) $$

which together with the super-operator property (5.9) leads to the following simple second order splitting

$$ \mathcal{T}^{(2)}(t + \tau, t) = \exp\left(\frac{\tau}{2}B(t + \tau)\right) \exp(\tau A) \exp\left(\frac{\tau}{2}B(t)\right). \quad (5.11) $$

In the same manner we extend both 6th and 8th order Yoshida integrators to

$$ \begin{aligned} \mathcal{T}^{(n)}(t + \tau, t) = & \exp(c_1 \tau B(t + b_1 \tau)) \exp(d_1 \tau A) \exp(c_2 \tau B(t + b_2 \tau)) \\ & \dots \exp(d_{2m+1} \tau A) \exp(c_{2m+2} \tau B(t + b_{2m+2} \tau)), \quad (5.12) \end{aligned} $$

where $b_k$ are constants given by

$$ b_k = \sum_{j=k}^{2m+1} d_j = 1 - \sum_{j=1}^{k-1} d_j, $$

in particular, $b_1 = 1$ and $b_{2m+2} = 0$.

In the case of time dependent potential there is an analogue of (5.4) in [16] reading

$$ \mathcal{T}^{(4)}(t+\tau,t) = \exp\left(\frac{\tau}{6}B(t+\tau)\right) \exp\left(\frac{\tau}{2}A\right) \exp\left(\frac{2\tau}{3}\tilde{B}\left(t+\frac{\tau}{2},\tau\right)\right) \exp\left(\frac{\tau}{2}A\right) \exp\left(\frac{\tau}{6}B(t)\right), \quad (5.13) $$

where the modified potential given at the middle step time moment is defined as

$$ \tilde{B}(t, \tau) = B(t) + \frac{\tau^2}{48} [B(t), [A, B(t)]]. \quad (5.14) $$

This follows from (5.10), (5.4), (5.9) and the identity $[B(t), [\mathcal{D}, B(t)]] = 0$. The latter can be checked directly

$$ F(t)[B(t), [\mathcal{D}, B(t)]] = 2F(t)B(t)\mathcal{D}B(t) - F(t)B^2(t)\mathcal{D} - F(t)\mathcal{D}B^2(t) = 0 $$

on any smooth test operator $F(t)$. Scheme (5.13) is very useful, since the commutator in (5.14) turns out to be proportional to $\partial_x V(x, t)$ for the Schrödinger equation, compare also with the same scheme obtained recently for the Dirac equation [41]. Similarly, the sixth order integrator (5.7) turns into

$$ \begin{aligned} \mathcal{T}^{(6)}(t + \tau, t) = & \exp\left(\frac{\tau_4}{6}B(t + (2 - s_4)\tau_4)\right) \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{2\tau_4}{3}\tilde{B}\left(t + \left(\frac{3}{2} - s_4\right)\tau_4, \tau_4\right)\right) \\ & \cdot \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{(1 - s_4)\tau_4}{6}B(t + (1 - s_4)\tau_4)\right) \exp\left(-\frac{s_4\tau_4}{2}A\right) \\ & \cdot \exp\left(-\frac{2s_4\tau_4}{3}\tilde{B}\left(t + \left(1 - \frac{s_4}{2}\right)\tau_4, s_4\tau_4\right)\right) \exp\left(-\frac{s_4\tau_4}{2}A\right) \exp\left(\frac{(1 - s_4)\tau_4}{6}B(t + \tau_4)\right) \\ & \cdot \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{\tau_4}{3}\tilde{B}\left(t + \frac{\tau_4}{2}, \tau_4\right)\right) \exp\left(\frac{\tau_4}{2}A\right) \exp\left(\frac{\tau_4}{6}B(t)\right). \quad (5.15) \end{aligned} $$

Note that we first extended the time-independent potential scheme (5.4) to (5.7) via recursion (5.6), afterwards we applied the Suzuki decomposition method (5.10) to deduce (5.15) from (5.7). Integrator (5.13) can also be directly extended to higher order schemes by the general Suzuki's recursive procedure

$$ \begin{aligned} \mathcal{T}^{(n+2)}(t + \tau, t) = & \mathcal{T}^{(n)}(t + \tau, t + (1 - \kappa_n)\tau) \mathcal{T}^{(n)}(t + (1 - \kappa_n)\tau, t + (1 - 2\kappa_n)\tau) \\ & \cdot \mathcal{T}^{(n)}(t + (1 - 2\kappa_n)\tau, t + 2\kappa_n\tau) \mathcal{T}^{(n)}(t + 2\kappa_n\tau, t + \kappa_n\tau) \mathcal{T}^{(n)}(t + \kappa_n\tau, t) \end{aligned} \quad (5.16) $$
with $\kappa_n = (4 - 4^{1/(n+1)})^{-1}$, mentioned in [40]. However, it very soon becomes expensive, for instance, 8-th order scheme demands 50 applications of kinetic energy exponent per time step. So we will only implement it once to derive $\mathcal{T}^{(6)}$ from $\mathcal{T}^{(4)}$ given in (5.13). The corresponding integrator reads

$$ 
\begin{aligned}
\mathcal{T}^{(6)}(t + \tau, t) = & \exp\left(\frac{\kappa_4 \tau}{6} B(t + \tau)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \left(1 - \frac{\kappa_4}{2}\right)\tau, \kappa_4 \tau\right)\right) \\
& \cdot \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{\kappa_4 \tau}{3} B(t + (1 - \kappa_4)\tau)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \left(1 - \frac{3}{2}\kappa_4\right)\tau, \kappa_4 \tau\right)\right) \\
& \cdot \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{1 - 3\kappa_4}{6}\tau B(t + (1 - 2\kappa_4)\tau)\right) \exp\left(\frac{1 - 4\kappa_4}{2}\tau A\right) \\
& \cdot \exp\left(\frac{2(1 - 4\kappa_4)\tau}{3} \tilde{B}\left(t + \frac{\tau}{2}, (1 - 4\kappa_4)\tau\right)\right) \exp\left(\frac{(1 - 4\kappa_4)\tau}{2} A\right) \exp\left(\frac{(1 - 3\kappa_4)\tau}{6} B(t + 2\kappa_4\tau)\right) \\
& \cdot \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \frac{3\kappa_4 \tau}{2}, \kappa_4 \tau\right)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{\kappa_4 \tau}{3} B(t + \kappa_4 \tau)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \\
& \cdot \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \frac{\kappa_4 \tau}{2}, \kappa_4 \tau\right)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{\kappa_4 \tau}{6} B(t)\right)
\end{aligned}
\quad (5.17) $$

with positive $\kappa_4 \approx 0.4$ and negative $1 - 4\kappa_4 \approx -0.5$ weights staying in front of kinetic operator. These weights together with the amount of the kinetic energy exponents makes (5.17) less desirable for use comparing to (5.15), though it may give better results for some time steps.

## 6. NUMERICAL EXPERIMENTS

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Matrix element $[\alpha_l^n]_{pj}$ with $p = j = 0$ at scale $n = 8$ corresponding to the time step $t = 0.0001$ and the MRA order $k = 6$.</figcaption></figure>

The free-particle semigroup $\exp(it\partial_x^2)$ is encoded in MRCPP (Multiresolution Computation Program Package) [2]. MRCPP is a high-performance C++ library, which provides various tools for working with Multiwavelets in connection to computational chemistry. It is also made available in VAMPyR (Very Accurate Multiwavelets Python Routines) [3], which builds on the MRCPP capabilities, providing an easy to use interface. As an introduction to this software, we refer to our recent article
with $\kappa_n = (4 - 4^{1/(n+1)})^{-1}$, mentioned in [40]. However, it very soon becomes expensive, for instance, 8-th order scheme demands 50 applications of kinetic energy exponent per time step. So we will only implement it once to derive $\mathcal{T}^{(6)}$ from $\mathcal{T}^{(4)}$ given in (5.13). The corresponding integrator reads

$$ 
\begin{aligned}
\mathcal{T}^{(6)}(t + \tau, t) = & \exp\left(\frac{\kappa_4 \tau}{6} B(t + \tau)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \left(1 - \frac{\kappa_4}{2}\right)\tau, \kappa_4 \tau\right)\right) \\
& \cdot \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{\kappa_4 \tau}{3} B(t + (1 - \kappa_4)\tau)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \left(1 - \frac{3}{2}\kappa_4\right)\tau, \kappa_4 \tau\right)\right) \\
& \cdot \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{(1 - 3\kappa_4)\tau}{6} B(t + (1 - 2\kappa_4)\tau)\right) \exp\left(\frac{(1 - 4\kappa_4)\tau}{2} A\right) \\
& \cdot \exp\left(\frac{2(1 - 4\kappa_4)\tau}{3} \tilde{B}\left(t + \frac{\tau}{2}, (1 - 4\kappa_4)\tau\right)\right) \exp\left(\frac{(1 - 4\kappa_4)\tau}{2} A\right) \exp\left(\frac{(1 - 3\kappa_4)\tau}{6} B(t + 2\kappa_4\tau)\right) \\
& \cdot \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \frac{3\kappa_4 \tau}{2}, \kappa_4 \tau\right)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{\kappa_4 \tau}{3} B(t + \kappa_4 \tau)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \\
& \cdot \exp\left(\frac{2\kappa_4 \tau}{3} \tilde{B}\left(t + \frac{\kappa_4 \tau}{2}, \kappa_4 \tau\right)\right) \exp\left(\frac{\kappa_4 \tau}{2} A\right) \exp\left(\frac{\kappa_4 \tau}{6} B(t)\right)
\end{aligned}
\quad (5.17) $$

with positive $\kappa_4 \approx 0.4$ and negative $1 - 4\kappa_4 \approx -0.5$ weights staying in front of kinetic operator. These weights together with the amount of the kinetic energy exponents makes (5.17) less desirable for use comparing to (5.15), though it may give better results for some time steps.

## 6. NUMERICAL EXPERIMENTS

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Matrix element $[\alpha_l^n]_{pj}$ with $p = j = 0$ at scale $n = 8$ corresponding to the time step $t = 0.0001$ and the MRA order $k = 6$.</figcaption></figure>

The free-particle semigroup $\exp(it\partial_x^2)$ is encoded in MRCPP (Multiresolution Computation Program Package) [2]. MRCPP is a high-performance C++ library, which provides various tools for working with Multiwavelets in connection to computational chemistry. It is also made available in VAMPyR (Very Accurate Multiwavelets Python Routines) [3], which builds on the MRCPP capabilities, providing an easy to use interface. As an introduction to this software, we refer to our recent article
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Matrix element $[\beta_l^n]_{pj}$ with $p = 0, j = 5$ at scale $n = 8$ corresponding to the time step $t = 0.0001$ and the MRA order $\ell = 6$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Difference between the analytical solution (6.1) and the numerically propagated one at the time moment $t = 0.0001$ for the MRA order $\ell = 20$ and the precision $\varepsilon = 10^{-10}$.</figcaption></figure>

[11], where we have already included a simple numerical experiment with the exponential operator $\exp(it\partial_x^2)$. In [11] we did not exploit the sparsity analysed thoroughly here in Subsection 4.1. In order to construct the multiresolution of $\exp(it\partial_x^2)$ in an adaptive way, the threshold is set to $\varepsilon/10$ for a given precision $\varepsilon$. In other words matrices $\alpha_l^n, \beta_l^n, \gamma_l^n$ having the Frobenius norm (4.22) less or equal than $\varepsilon/10$, are ignored in numerical calculations. This choice of threshold is motivated
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Matrix element $[\beta_l^n]_{pj}$ with $p = 0, j = 5$ at scale $n = 8$ corresponding to the time step $t = 0.0001$ and the MRA order $\ell = 6$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Difference between the analytical solution (6.1) and the numerically propagated one at the time moment $t = 0.0001$ for the MRA order $\ell = 20$ and the precision $\varepsilon = 10^{-10}$.</figcaption></figure>

[11], where we have already included a simple numerical experiment with the exponential operator $\exp(it\partial_x^2)$. In [11] we did not exploit the sparsity analysed thoroughly here in Subsection 4.1. In order to construct the multiresolution of $\exp(it\partial_x^2)$ in an adaptive way, the threshold is set to $\varepsilon/10$ for a given precision $\varepsilon$. In other words matrices $\alpha_l^n, \beta_l^n, \gamma_l^n$ having the Frobenius norm (4.22) less or equal than $\varepsilon/10$, are ignored in numerical calculations. This choice of threshold is motivated
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Convergence comparison between FFT-based schemes (solid lines) and MRA-based schemes (large dots) in the harmonic potential example. The convergence rates are indicated by the dashed lines. The legend corresponds to different integrators: (S2) for (1.4), (A4) for (5.4), (A6) for (5.7), (Y6) for (5.8) with $n = 6$ and (Y8) for (5.8) with $n = 8$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Convergence comparison between FFT-based schemes (solid lines) and MRA-based schemes (large dots) in the Morse potential example. The convergence rates are indicated by the dashed lines. The legend corresponds to different integrators: (S2) for (5.11), (A4) for (5.13), (A6) for (5.15), (W6) for (5.17), (Y6) for (5.12) with $n = 6$ and (Y8) for (5.12) with $n = 8$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Convergence comparison between FFT-based schemes (solid lines) and MRA-based schemes (large dots) in the harmonic potential example. The convergence rates are indicated by the dashed lines. The legend corresponds to different integrators: (S2) for (1.4), (A4) for (5.4), (A6) for (5.7), (Y6) for (5.8) with $n = 6$ and (Y8) for (5.8) with $n = 8$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Convergence comparison between FFT-based schemes (solid lines) and MRA-based schemes (large dots) in the Morse potential example. The convergence rates are indicated by the dashed lines. The legend corresponds to different integrators: (S2) for (5.11), (A4) for (5.13), (A6) for (5.15), (W6) for (5.17), (Y6) for (5.12) with $n = 6$ and (Y8) for (5.12) with $n = 8$.</figcaption></figure>
by oscillatory behaviour of the matrix entries $[\alpha_l^n]_{pj}$, $[\beta_l^n]_{pj}$, $[\gamma_l^n]_{pj}$ with respect to the distance to diagonal $l$, on the one hand, see Figures 7, 8 and compare with Figure 2. On the other hand, the result of operator application was controlled by the free-particle analytical solution to (1.3) of the form

$$u(x,t) = \exp (it\partial_x^2) \psi_0(x) = \left(\frac{1}{2\pi\sigma^2}\right)^{1/4} \sqrt{\frac{4\sigma^2}{4it + 4\sigma^2}} e^{-\frac{(x-x_0)^2}{4it+4\sigma^2}} \quad (6.1)$$

associated with the initial condition

$$\psi_0(x) = \left(\frac{1}{2\pi\sigma^2}\right)^{1/4} \exp\left(-\frac{(x-x_0)^2}{4\sigma^2}\right) \quad (6.2)$$

in the sense $u(x, 0) = \psi_0(x)$. Indeed, let $x_0 = 0.5$ and $\sigma = 0.04$, for example. Then for the MRA order $k = 20$, time parameter $t = 0.0001$ and precision $\epsilon = 10^{-10}$ we have the $L^2$-norm of the difference between analytical and numerically propagated solutions equaled $2.0 \cdot 10^{-11}$, see Figure 9.

This concludes numerical description of the free-particle propagator. Below we illustrate its use on two toy model problems. Its performance is compared to the periodic spectral method.

### 6.1. Harmonic potential.
Our first numerical example takes the Schrödinger equation

$$i\partial_t\psi = -\frac{1}{2}\partial_x^2\psi + V(x)\psi$$

with the harmonic potential

$$V(x) = V_0 \left(x - \frac{1}{2}\right)^2.$$

It is complemented by the initial condition $\psi(0) = \psi_0$ having the Gaussian form (6.2). It is well known that in the harmonic potential the density $|\psi(t)|^2$ oscillates with the period $t_{per} = \pi\sqrt{2/V_0}$, namely, $\psi(t_{per}) = -\psi_0$. This comes from the fact that the eigenvalues for the Hamiltonian are $\sqrt{2V_0}(n + 1/2)$.

Ideally, a choice of parameters like $V_0$ should be physically motivated. Then given some tolerance one can set the computational domain of size $L$ centered, for our example, at $x = 1/2$. Transforming the problem to the unit space interval $[0, 1]$, one will obviously get that the final $V_0$ turns out to be proportional to $L^4$. Therefore, we set the following parameters $V_0 = 98304$, $x_0 = 0.375$ and $\sigma = 0.025$ associated with the domain $[0, 1]$.

We conduct numerical simulations with time step $\tau$ ranging from $t_{per}/10$ to $t_{per}/160$ up to the final time moment $t = t_{per}$ using the machinery developed here. The integrators in use are referred as follows

(S2) the simple second order scheme (1.4),

(A4) the Chin-Chen's scheme A of fourth order (5.4),

(A6) the Chin-Chen's scheme A of sixth order (5.7),

(Y6) the Yoshida's scheme (5.8) with the order $n = 6$,

(Y8) the Yoshida's scheme (5.8) with the order $n = 8$.

The accuracy of the obtained solution $\psi$ for each of these schemes is estimated by the $L^2$-difference

$$\text{error}(t) = \|\psi(t) - \psi_{\text{ref}}(t)\|_{L^2(0,1)}. \quad (6.3)$$

Here we calculate the error at the moment $t = t_{per}$ with respect to the reference $\psi_{\text{ref}}(t_{per}) = -\psi_0$. For the multiwavelet calculations we take the MRA order $k = 18$ and the tolerance needed for adaptive resolution to be $\epsilon = 10^{-10}$. The performance is compared with the Fourier spectral method calculations using 1024 grid points inside of the unit interval $[0, 1]$.

The electron dynamics in the harmonic potential well is smooth and it can be viewed as periodic with very high precision. Therefore, there is no surprise that the spectral method outperforms the multiresolution technique in Figure 10. The accuracy stagnation for the FFT based schemes reveals itself at four times smaller time steps than the corresponding stagnation for the MRA based schemes, with $\tau = t_{per}/640 \approx 2.2 \cdot 10^{-5}$ compared to $\tau = t_{per}/160 \approx 8.9 \cdot 10^{-5}$ for the
by oscillatory behaviour of the matrix entries $[\alpha_l^n]_{pj}$, $[\beta_l^n]_{pj}$, $[\gamma_l^n]_{pj}$ with respect to the distance to diagonal $l$, on the one hand, see Figures 7, 8 and compare with Figure 2. On the other hand, the result of operator application was controlled by the free-particle analytical solution to (1.3) of the form

$$u(x,t) = \exp (it\partial_x^2) \psi_0(x) = \left(\frac{1}{2\pi\sigma^2}\right)^{1/4} \sqrt{\frac{4\sigma^2}{4it + 4\sigma^2}} e^{-\frac{(x-x_0)^2}{4it+4\sigma^2}} \quad (6.1)$$

associated with the initial condition

$$\psi_0(x) = \left(\frac{1}{2\pi\sigma^2}\right)^{1/4} \exp\left(-\frac{(x-x_0)^2}{4\sigma^2}\right) \quad (6.2)$$

in the sense $u(x, 0) = \psi_0(x)$. Indeed, let $x_0 = 0.5$ and $\sigma = 0.04$, for example. Then for the MRA order $k = 20$, time parameter $t = 0.0001$ and precision $\varepsilon = 10^{-10}$ we have the $L^2$-norm of the difference between analytical and numerically propagated solutions equaled $2.0 \cdot 10^{-11}$, see Figure 9.

This concludes numerical description of the free-particle propagator. Below we illustrate its use on two toy model problems. Its performance is compared to the periodic spectral method.

### 6.1. Harmonic potential.
Our first numerical example takes the Schrödinger equation

$$i\partial_t\psi = -\frac{1}{2}\partial_x^2\psi + V(x)\psi$$

with the harmonic potential

$$V(x) = V_0 \left(x - \frac{1}{2}\right)^2.$$

It is complemented by the initial condition $\psi(0) = \psi_0$ having the Gaussian form (6.2). It is well known that in the harmonic potential the density $|\psi(t)|^2$ oscillates with the period $t_{per} = \pi\sqrt{2/V_0}$, namely, $\psi(t_{per}) = -\psi_0$. This comes from the fact that the eigenvalues for the Hamiltonian are $\sqrt{2V_0}(n + 1/2)$.

Ideally, a choice of parameters like $V_0$ should be physically motivated. Then given some tolerance one can set the computational domain of size $L$ centered, for our example, at $x = 1/2$. Transforming the problem to the unit space interval $[0, 1]$, one will obviously get that the final $V_0$ turns out to be proportional to $L^4$. Therefore, we set the following parameters $V_0 = 98304$, $x_0 = 0.375$ and $\sigma = 0.025$ associated with the domain $[0, 1]$.

We conduct numerical simulations with time step $\tau$ ranging from $t_{per}/10$ to $t_{per}/160$ up to the final time moment $t = t_{per}$ using the machinery developed here. The integrators in use are referred as follows

(S2) the simple second order scheme (1.4),

(A4) the Chin-Chen's scheme A of fourth order (5.4),

(A6) the Chin-Chen's scheme A of sixth order (5.7),

(Y6) the Yoshida's scheme (5.8) with the order $n = 6$,

(Y8) the Yoshida's scheme (5.8) with the order $n = 8$.

The accuracy of the obtained solution $\psi$ for each of these schemes is estimated by the $L^2$-difference

$$\text{error}(t) = \|\psi(t) - \psi_{\text{ref}}(t)\|_{L^2(0,1)}. \quad (6.3)$$

Here we calculate the error at the moment $t = t_{per}$ with respect to the reference $\psi_{\text{ref}}(t_{per}) = -\psi_0$. For the multiwavelet calculations we take the MRA order $k = 18$ and the tolerance needed for adaptive resolution to be $\varepsilon = 10^{-10}$. The performance is compared with the Fourier spectral method calculations using 1024 grid points inside of the unit interval $[0, 1]$.

The electron dynamics in the harmonic potential well is smooth and it can be viewed as periodic with very high precision. Therefore, there is no surprise that the spectral method outperforms the multiresolution technique in Figure 10. The accuracy stagnation for the FFT based schemes reveals itself at four times smaller time steps than the corresponding stagnation for the MRA based schemes, with $\tau = t_{per}/640 \approx 2.2 \cdot 10^{-5}$ compared to $\tau = t_{per}/160 \approx 8.9 \cdot 10^{-5}$ for the
integrators (A6), (Y6), (Y8). Otherwise, the precision is consistent and the convergence is optimal for each scheme, as indicated by the slopes of curves in Figure 10 corresponding to their theoretical counterparts. It is also interesting to notice that (Y8) does not reach the accuracy of the cheaper scheme (A6) in the current setup.

**6.2. Walker-Preston model.** In order to test the multiwavelet machinery on the problems with time-dependent Hamiltonian, we consider the Walker-Preston model [39] of a diatomic molecule in a strong laser field. The equation reads

$$i\partial_t\psi = -\frac{1}{2}\partial_x^2\psi + V(x)\psi + \mathcal{E}(x,t)\psi,$$

where $V(x)$ stands for the Morse potential

$$V(x) = V_0(1 - e^{-\alpha x})^2$$

with $V_0 = 0.2251\mu$, $\mu = 1745$ and $\alpha = 1.1741$. These parameters correspond to the HF molecule. Here $\mathcal{E}(x,t) = Ax\cos\omega t$ is the external field with $A = 0.011025\mu$ and $\omega = 0.01787\mu$. The first $0 \to 1$ vibrational frequency difference in the Morse oscillator is $0.01807\mu$, so $\omega$ is slightly below resonance. As the initial condition $\psi(0) = \psi_0$ we consider the ground state wave function

$$\psi_0(x) = \sqrt{\alpha \frac{\xi^{2\lambda-1}(x)}{\Gamma(2\lambda-1)}} e^{-\xi(x)}, \quad \xi(x) = 2\lambda e^{-\alpha x}, \quad \lambda = \frac{\sqrt{2V_0}}{\alpha}$$

having the ground state energy

$$E_0 = \alpha \sqrt{\frac{V_0}{2}} - \frac{\alpha^2}{8}.$$

It is worth to point out that a direct evaluation of $\psi_0(x)$ may lead to huge round up errors. Therefore, we rewrite it in the form suitable for numerical calculations. There is an integer $m \in \mathbb{N}$ such that $2\lambda - 1 = m + \delta\lambda$ with $1/2 \le \delta\lambda < 3/2$. So we can factorise $\psi_0$ as

$$\psi_0(x) = \sqrt{\alpha \frac{\xi^{\delta\lambda}}{\Gamma(\delta\lambda)} \frac{\xi e^{-\xi/m}}{\delta\lambda} \frac{\xi e^{-\xi/m}}{\delta\lambda + 1} \cdots \frac{\xi e^{-\xi/m}}{\delta\lambda + m - 1}},$$

so eliminating possibility of rounding error.

All the introduced parameters are associated with the computational domain $[-0.8, 4.32]$ of the size $L = 5.12$, used by many authors [12, 16, 22, 31, 32] for testing different numerical methods. Thus it is known that a use of the spectral periodic treatment for this problem with 64 grid points gives very accurate results, that can be regarded as a perfect benchmark for our purposes. However, in order to be able to exploit the multiresolution technique one needs to translate the formulation to the computational domain $[0, 1]$ with an affine transformation. Henceforth $x \in [0, 1]$ and the time scales $t = t^{\text{original}}/L^2$. All the other units are changed accordingly, for example, the dissociation energy $V_0 = V_0^{\text{original}}L^2 = 0.2251\mu L^2$.

Spectral FFT-based simulations are conducted with time step $\tau$ ranging from 0.001 down to $0.001 \cdot 2^{-10}$. The final time moment is set to $t_{\text{fin}} = 0.1$. The integrators in use are referred as follows

(S2) the simple second order scheme (5.11),

(A4) the Chin-Chen's scheme A of fourth order (5.13),

(A6) the Chin-Chen's scheme A of sixth order (5.15),

(W6) the Chin-Chen's scheme A of sixth order (5.17) obtained via the procedure (5.16),

(Y6) the Yoshida's scheme (5.12) with the order $n = 6$,

(Y8) the Yoshida's scheme (5.12) with the order $n = 8$.

The accuracy of the obtained solution $\psi$ for each of these schemes is estimated by (6.3) at the final time moment $t = t_{\text{fin}}$. Now we do not have an analytical result that could stand in place of $\psi_{\text{ref}}$ as in the previous experiment. Presuming we do not know a priori which scheme is the most precise, we take different $\psi_{\text{ref}}$ for each scheme. It is defined by the tiniest small step in use and according
integrators (A6), (Y6), (Y8). Otherwise, the precision is consistent and the convergence is optimal for each scheme, as indicated by the slopes of curves in Figure 10 corresponding to their theoretical counterparts. It is also interesting to notice that (Y8) does not reach the accuracy of the cheaper scheme (A6) in the current setup.

6.2. **Walker-Preston model.** In order to test the multiwavelet machinery on the problems with time-dependent Hamiltonian, we consider the Walker-Preston model [39] of a diatomic molecule in a strong laser field. The equation reads

$$i\partial_t\psi = -\frac{1}{2}\partial_x^2\psi + V(x)\psi + \mathcal{E}(x,t)\psi,$$

where $V(x)$ stands for the Morse potential

$$V(x) = V_0(1 - e^{-\alpha x})^2$$

with $V_0 = 0.2251\mu$, $\mu = 1745$ and $\alpha = 1.1741$. These parameters correspond to the HF molecule. Here $\mathcal{E}(x,t) = Ax\cos\omega t$ is the external field with $A = 0.011025\mu$ and $\omega = 0.01787\mu$. The first $0 \to 1$ vibrational frequency difference in the Morse oscillator is $0.01807\mu$, so $\omega$ is slightly below resonance. As the initial condition $\psi(0) = \psi_0$ we consider the ground state wave function

$$\psi_0(x) = \sqrt{\alpha \frac{\xi^{2\lambda-1}(x)}{\Gamma(2\lambda-1)}} e^{-\xi(x)}, \quad \xi(x) = 2\lambda e^{-\alpha x}, \quad \lambda = \frac{\sqrt{2V_0}}{\alpha}$$

having the ground state energy

$$E_0 = \alpha \sqrt{\frac{V_0}{2}} - \frac{\alpha^2}{8}.$$

It is worth to point out that a direct evaluation of $\psi_0(x)$ may lead to huge round up errors. Therefore, we rewrite it in the form suitable for numerical calculations. There is an integer $m \in \mathbb{N}$ such that $2\lambda - 1 = m + \delta\lambda$ with $1/2 \le \delta\lambda < 3/2$. So we can factorise $\psi_0$ as

$$\psi_0(x) = \sqrt{\alpha \frac{\xi^{\delta\lambda}}{\Gamma(\delta\lambda)} \frac{\xi e^{-\xi/m}}{\delta\lambda} \frac{\xi e^{-\xi/m}}{\delta\lambda + 1} \cdots \frac{\xi e^{-\xi/m}}{\delta\lambda + m - 1}},$$

so eliminating possibility of rounding error.

All the introduced parameters are associated with the computational domain $[-0.8, 4.32]$ of the size $L = 5.12$, used by many authors [12, 16, 22, 31, 32] for testing different numerical methods. Thus it is known that a use of the spectral periodic treatment for this problem with 64 grid points gives very accurate results, that can be regarded as a perfect benchmark for our purposes. However, in order to be able to exploit the multiresolution technique one needs to translate the formulation to the computational domain $[0, 1]$ with an affine transformation. Henceforth $x \in [0, 1]$ and the time scales $t = t^{\text{original}}/L^2$. All the other units are changed accordingly, for example, the dissociation energy $V_0 = V_0^{\text{original}}L^2 = 0.2251\mu L^2$.

Spectral FFT-based simulations are conducted with time step $\tau$ ranging from 0.001 down to $0.001 \cdot 2^{-10}$. The final time moment is set to $t_{\text{fin}} = 0.1$. The integrators in use are referred as follows

(S2) the simple second order scheme (5.11),

(A4) the Chin-Chen's scheme A of fourth order (5.13),

(A6) the Chin-Chen's scheme A of sixth order (5.15),

(W6) the Chin-Chen's scheme A of sixth order (5.17) obtained via the procedure (5.16),

(Y6) the Yoshida's scheme (5.12) with the order $n = 6$,

(Y8) the Yoshida's scheme (5.12) with the order $n = 8$.

The accuracy of the obtained solution $\psi$ for each of these schemes is estimated by (6.3) at the final time moment $t = t_{\text{fin}}$. Now we do not have an analytical result that could stand in place of $\psi_{\text{ref}}$ as in the previous experiment. Presuming we do not know a priori which scheme is the most precise, we take different $\psi_{\text{ref}}$ for each scheme. It is defined by the tiniest small step in use and according
to the scheme in question. In other words, $\psi_{\text{ref}}$ for (S2) numerical solutions is the (S2) FFT-based numerical solution associated to the time step $\tau = 0.001 \cdot 2^{-10}$, and so on. The same references are used later on in order to estimate the accuracy of the corresponding multiwavelet simulations.

In fact we can again see an obvious accuracy outperformance of the periodic spectral method in Figure 11. This time, however, we knew a priori that the best FFT calculations are achieved with the rough grid consisting of only 64 points. Indeed, refining the grid to 256 or 1024 one encounters poorer accuracy due to the Gibbs phenomenon leading to high oscillations of solutions close to the boundary points. It does not affect the multiresolution approach, of course. We do not include figures corresponding to refined grid FFT calculations here, as they exhibit behavior similar to Figure 11. However, these FFT results demonstrate poorer convergence, deviating from the theoretically predicted rate, before the MRA calculations approach stagnation. In other words, FFT-based spectral calculations can yield inferior results compared to MRA-based ones in terms of accuracy.

An interesting and noteworthy outcome is that the extensions (5.12) of Yoshida's schemes (5.8) for the time-dependent potential case degenerate to first order. According to [40] it may happen that some exponential integrators lose their theoretically predicted accuracy after the Suzuki's extension described in Subsection 5.2. It seems to us that Yoshida's schemes (5.8) cannot be efficiently extended to non-autonomous equations. The first author encountered a similar problem before [19].

## 7. CONCLUSION

We have developed a multiresolution analysis of the semigroup $\exp(it\partial_x^2)$ and introduced MRA-based numerical methods for solving the time-dependent Schrödinger equation. These methods offer several notable advantages. First of all, they can be seamlessly integrated with existing tools for solving stationary problems. Specifically, the application of singular integral operators arising from potential $V$, as utilized in Kohn-Sham and Hartree-Fock theories, is already well-established in multiwavelet bases. Consequently, these techniques can be readily adapted to address time evolution problems. Secondly, these methods offer spatial adaptivity, which is particularly crucial for quantum chemistry simulations. In these simulations, it is essential to handle Coulomb-type singularities and solutions with cusps. Therefore, an automatic grid refinement based on specified precision can be extremely beneficial.

There are some drawbacks to consider. Achieving sparsity in the multiresolution representation of $\exp(it\partial_x^2)$, which is obviously needed for efficient computations, requires large time steps $t$ and a high MRA order $\ell$. However, a use of a too large MRA order $\ell$ can slow down simulations, while overly large time steps $t$ may compromise the accuracy of the results obtained.

In the future, we plan to focus on developing artificial boundary conditions for the MRA-based methods introduced here. This approach will enable us to use larger time steps while maintaining accuracy, by confining the computational domain to a relatively compact size. Additionally, it will facilitate the use of reasonable values for $\ell$. The techniques developed will be applied to problems where FFT-based methods encounter challenges, specifically in modeling the attosecond dynamics of electrons in molecules exposed to intense laser pulses.

## DECLARATIONS

The CRediT taxonomy of contributor roles [5, 13] is applied. The “Investigation” role also includes the “Methodology”, “Software”, and “Validation” roles. The “Analysis” role also includes the “Formal analysis” and “Visualization” roles. The “Funding acquisition” role also includes the “Resources” role. The contributor roles are visualized in the following authorship attribution matrix, as suggested in Ref. [1].

**Acknowledgments.** The authors are grateful to S. R. Jensen and Ch. Tantardini for numerous helpful discussions. The research was supported by the Research Council of Norway through its Centres of Excellence scheme, Hylleraas-senteret, project number 262695.