# A symmetric low-regularity integrator for the nonlinear Schrödinger equation

Yvonne Alama Bronsard

► To cite this version:

Yvonne Alama Bronsard. A symmetric low-regularity integrator for the nonlinear Schrödinger equation. IMA Journal of Numerical Analysis, 2023, 10.1093/imanum/drad093. hal-04212978v2

HAL Id: hal-04212978

https://hal.science/hal-04212978v2

Submitted on 12 Feb 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# A symmetric low-regularity integrator for the nonlinear Schrödinger equation

Yvonne Alama Bronsard *

## Abstract

We introduce and analyze a symmetric low-regularity scheme for the nonlinear Schrödinger (NLS) equation beyond classical Fourier-based techniques. We show fractional convergence of the scheme in $L^2$-norm, from first up to second order, both on the torus $T^d$ and on a smooth bounded domain $\Omega \subset \mathbb{R}^d$, $d \le 3$, equipped with homogeneous Dirichlet boundary condition. The new scheme allows for a symmetric approximation to the NLS equation in a more general setting than classical splitting, exponential integrators, and low-regularity schemes (i.e. under lower regularity assumptions, on more general domains, and with fractional rates). We motivate and illustrate our findings through numerical experiments, where we witness better structure preserving properties and an improved error-constant in low-regularity regimes.

## 1 Introduction

We consider the nonlinear Schrödinger (NLS) equation,

$$i\partial_t u(t,x) = -\Delta u(t,x) + |u(t,x)|^2 u(t,x), \quad (t,x) \in \mathbb{R} \times \Omega \quad (1)$$

with $\Omega \subset \mathbb{R}^d$ or $\Omega = \mathbb{T}^d$, $d \le 3$, and an initial condition

$$u|_{t=0} = u_0. \quad (2)$$

When $\partial\Omega \neq \emptyset$, we assume that $\Omega$ is a smooth bounded domain and we assign homogeneous boundary conditions which will be encoded in the choice of the domain of the operator $i\Delta$. In the convergence analysis we will consider either periodic or homogeneous Dirichlet boundary conditions. Nevertheless, one could also consider different types of boundary conditions such as homogeneous Neumann boundary conditions by defining the functional spaces accordingly (see Section 2).

Throughout this article we will be interested in providing a reliable approximation of (1) when the initial data $u_0$ are non-smooth, in the sense that they belong to Sobolev spaces of low order. Namely, we will be interested in studying numerical schemes which approximate the time dynamics of (1) at low regularity. The numerical study of low-regularity approximations to nonlinear evolution equations has gained lots of attention in the past years, and numerous contributions have been made in this direction. The first results were established on the torus for the Korteweg-De Vries (KdV) equation and then the NLS equation with the pioneering works of [26] and [35]. These results could be further greatly extended, see for instance [45, 34, 39, 31] and [37, 43, 30, 15] for the KdV and NLS equations respectively. More types of dispersive equations could be dealt with, including for example the Dirac equation [41] or the “good” Boussinesq

<footnote>*Sorbonne Université, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France.
email: yvonne.alama_bronsard@sorbonne-universite.fr</footnote>
equation ([29]), and a general framework for constructing low-regularity approximations up to arbitrary order and for a class of dispersive equations on the torus was obtained in [14].

The construction of these time integrators (called resonance-based schemes, exponential-type low-regularity integrators, or Fourier integrators) strongly depended on Fourier-based expansions, and hence were restricted to periodic boundary conditions. Recently, this restriction was withdrawn to treat more general domains $\Omega \subset \mathbb{R}^d$ and boundary conditions, as well as more general nonlinearities (see [40, 2, 3, 6, 28]).

The general aim of low-regularity integrators is that they converge under lower regularity assumptions, contrarily to classical methods (see [3] for a comparative analysis on general smooth domains). Their major drawback is that they do not preserve the geometric structure of the underlying system. The NLS equation (1) is time reversible, meaning that $\overline{u(-t, x)}$ is again solution of (1), and both the density and the energy are conserved quantities:

$$ \left\{ \begin{array}{l} \|u(t)\|_{L^2} = \|u_0\|_{L^2}, \\ E(t) = E_0, \end{array} \right. \quad t \in I, \qquad (3) $$

where $I$ is the interval of existence of the solution, and where for $H^1$-solutions we have

$$ E(t) = \frac{1}{2} \int |\nabla u|^2(t, x) dx + \frac{1}{4} \int |u|^4(t, x) dx. $$

Hence, when designing a numerical scheme it is natural to take into account both of these conserved quantities, and to retain (as much as possible) these properties also on the discrete level by introducing so-called structure preserving schemes, see [23] for an extensive introduction on the subject. The latter has received great interest thanks to their good long-time near-preservation of the actions of the integrable properties of the equation, and have been successfully studied in the past for the approximation of the NLS equation (1). Examples of such schemes are splitting schemes ([32], [20]), relaxation finite difference type schemes ([11]), symmetric exponential integrators ([17]) or Crank-Nicolson Galerkin methods ([24]), just to name a couple of them. For an overview of symmetric methods for NLS see [17, 20, 5]. While these classical structure preserving schemes provide excellent approximations to smooth solutions in general even up to long times, they often break down and lead to severe loss of convergence for non-smooth solutions. Low-regularity integrators which are suited for non-smooth solutions on the other hand do not preserve the structure of the underlying equation. The natural question which thus arises is: What about low-regularity structure preserving schemes for solving the NLS equation (1)? Only very little is known in this direction, see the work of [33] on the KdV equation, [44] on the cubic Klein-Gordon equation, and [7] on the isotropic Landau-Lifschitz equation. Also worth to be mentioned is the work of [43] which introduces for the first time a first-order Fourier integrator for the NLS equation (1) set on $\mathbb{T}$ which almost conserves the mass.

In this article, we introduce a symmetric low-regularity integrator for solving the NLS equation (1) which allows for low-regularity approximation while maintaining good long-time preservation of the two conserved properties (3) on the discrete level. We carry-out a rigorous convergence analysis in $L^2(\Omega)$ on smooth domains $\Omega \subset \mathbb{R}^d$ and obtain improved error estimates at low-regularity compared to classical symmetric methods. Our numerical findings not only show better structure preservation properties but also show a much better error constant at low-regularity than previously proposed methods (see Figure 1).

In the finite dimensional ODE setting it is well-known that symmetric methods are of even order. In the context of PDEs this is a much more delicate question as convergence is met only when sufficient regularity assumptions are imposed on the solution. Thanks to the gain of symmetry, we show second order $L^2$-convergence of the symmetric scheme under less regularity assumptions than what is required by classical
symmetric schemes ([32, 11, 12, 24]), while asking for slightly more regularity than asymmetric second-order low-regularity schemes ([2]) which however do not preserve the structure of the system (see Figures 2 and 3). Optimal first order low-regularity convergence rates could be obtained. See Section 1.2 for a detailed discussion on the subject.

The scheme we present here is based on the first-order low-regularity scheme first introduced in [35] which is given by,

$$ \Phi_{\tau}(u^n) := e^{i\tau\Delta} (u^n - i\tau(u^n)^2\varphi_1(-2i\tau\Delta)\overline{u^n}), \quad u^0 = u_0, \qquad (4) $$

where $\varphi_1(z) = \frac{e^z - 1}{z}$, and $\tau$ is the time step. In order to symmetrize the above scheme we introduce the adjoint method as the map

$$ \hat{\Phi}_{\tau} = \Phi_{-\tau}^{-1}, $$

and compute (see [23])

$$ \hat{\Phi}_{\tau/2} \circ \Phi_{\tau/2}. $$

This yields the following implicit symmetric low-regularity scheme,

$$ \begin{aligned} u^{n+1} &= \varphi^{\tau}(u^n) = e^{i\tau\Delta}u^n - i\frac{\tau}{2}e^{i\tau\Delta}((u^n)^2\varphi_1(-i\tau\Delta)\overline{u^n}) - i\frac{\tau}{2}((u^{n+1})^2\varphi_1(i\tau\Delta)\overline{u^{n+1}}) \quad (5) \\ &= e^{i\tau\Delta}u^n + \psi_E^{\tau/2}(u^n) + \psi_I^{\tau/2}(u^{n+1}) \\ &= e^{i\tau\Delta}u^n + \Psi^{\tau}(u^n, u^{n+1}), \end{aligned} $$

which satisfies the discrete analogue of the time-reversible property of (1).

We highlight the properties which the scheme (5) inherits through numerical experiments, where we couple the time-integrators with the standard Fourier pseudo-spectral method which encodes periodic boundary conditions. The case of homogeneous Dirichlet boundary conditions remains very similar, yet for completeness we also include a convergence plot in this case where we expand the solution as a sine series expansion. First, in the case of periodic boundary conditions, we observe in Figures 1a and 1b the favorable convergence properties of the scheme (5) for $H^1$ and $H^2$ data respectively. We notice that the error constant of the symmetric scheme is much better than the asymmetric first-order low-regularity integrator (Low-reg 1), and is also better than the asymmetric second-order low-regularity integrator (Low-reg 2). Figure 1c similarly shows the favorable convergence behaviour when considering homogeneous Dirichlet boundary conditions. Secondly, we study in Figures 2 and 3 the structure preserving properties of the new symmetric low-regularity integrator (5) against previous asymmetric low-regularity integrators. We witness that the asymmetric first and second-order low-regularity integrators (Low-reg 1, Low-reg 2) are unable to preserve the density and energy (see (3)), whereas the symmetric integrator (5) appears to nearly-preserve both conserved properties over long-times. We note that in the finite dimensional ODE setting a general theory for symmetric methods applied to integrable reversible systems has been established in [23] allowing for long-time near-conservation of first-integrals. In the infinite dimensional case the understanding of the long-time behaviour of numerical solutions is an ongoing challenge in the field of geometric integration and few results are known, see for example [20, 21, 22, 18]. We expect that it would be possible to prove long-time near-preservation of the density and energy of the scheme (5) by using the results of [23], and by benefitting of an analysis using modulated Fourier expansions (see [18, 22]) or using normal form techniques ([21, 10]) to show near-conservation of the energy. This delicate analysis is out of scope for this paper, where here we focus on the low-regularity error estimates on the solution itself. Finally, we refer to Figure 4 for a broad indication of the relative
computational cost of each of the three low-regularity integrators, and discuss the added cost of implementing the symmetric implicit scheme (5). We observe that the asymmetric second-order low-regularity integrator (Low-reg 2) costs in CPU-time approximately the same as the symmetric integrator (5). Whereas when comparing with the asymmetric first-order scheme (Low-reg 1) we have that the improved convergence properties of the scheme (5) make up for the extra cost of solving the implicit system (5) at every time step.

<figure><img src="image_2.png" /><figcaption>Figure 1: Convergence plot for data in $H^1$ (Figure (a)) and data in $H^2$ (Figure (b) and (c)) of the asymmetric first and second-order low-regularity integrators (pink, dark blue), the symmetric method (5) (red), the classical Lie splitting, Strang splitting, and Euler Exponential method (light blue, yellow, and green). We observe order reduction of the classical Euler Exponential and Lie and Strang splitting methods (Figure (a)), $H^1(\mathbb{T})$-data), and of the Strang splitting method (Figure (b)), $H^2(\mathbb{T})$-data, and Figure (c)), $H^2([0, 1])$-data). Figures (a) and (b) are with periodic boundary conditions, while Figure (c) is with homogeneous Dirichlet boundary conditions. The slopes of the continuous black lines are one and two, respectively. We took the final time $T = 1$, and the number of Fourier modes $K = 2^{11}$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Plot (a) : We graph the $L^2$-norm $||u(t)||_{L^2(T)}$ up until time $T = 50$ of three low-regularity integrators. The asymmetric first and second-order low-regularity integrators (pink, blue), and the new symmetric low-regularity integrator (red). We also graph the exact value $||u_0||_{L^2(T)}$ (black). Plot (b): We only graph the asymmetric second-order low-regularity integrators (blue), and the symmetric low-regularity integrator (red) together with the exact $L^2$-norm of the initial value (black). We fixed the number of Fourier modes $K = 2^9$, the time step $\tau = 0.05$, and took an initial data $u_0 \in H^2$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Plot (a) : We graph the relative energy $E(t)/E(t_0)$ (see (3)) up to $T = 50$ of the same three low-regularity integrators as in Figures 2a and 2b. Plot (b): We only graph the relative energy of the asymmetric second order low-regularity integrators (blue), and the symmetric low-regularity integrator (red). We again fixed the number of Fourier modes $K = 2^9$, the time step $\tau = 0.05$, and took an initial data $u_0 \in H^2$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: We plot the CPU time versus the L²-error. Namely, we compare the computational cost for running the first and second-order low-regularity scheme (pink, dark blue), with the low-regularity symmetric integrator (5) (red). We took the same parameter values as in Figure 1.</figcaption></figure>

This numerical study motivates the use of the scheme (5), which conserves better the underlying geometric structure of the equation, exhibits a better error constant, and can be implemented at relatively low additional cost despite its implicit nature.

**Remark 1.1** (Implicit versus explicit low-regularity schemes). We make the important remark that unlike the previous (asymmetric) low-regularity integrators [35, 14, 36, 37, 40, 2, 3] the above symmetrized scheme (5) is an implicit one. We have witnessed that the implicit nature of the scheme does not adversely affect the computational cost of the method (see Figure 4). Nevertheless, one could query on the necessity of the implicit nature of the symmetric low-regularity scheme. In the case of second-order wave-type equations, instead of considering implicit symmetrized schemes one could study explicit three-time step symmetric schemes using Gautschi-type methods. Indeed, for the cubic Klein-Gordon equation set on $\mathbb{T}$, an explicit symmetric three time-step low-regularity integrator could be obtained by [44]. While this approach is suited to second order equations, by combining the work of [44] together with the uniformly accurate low-regularity integrator [16] adapted to the non-relativistic regime, an interesting open problem would be to obtain in the non-relativistic limit an explicit symmetric three time-step low-regularity approximation to the NLS equation (1).

Having motivated the scheme (5) we now provide the underlying idea behind its construction. To provide a better intuition to the reader we will work in the twisted variable and place ourselves on the torus $\Omega = \mathbb{T}$ in order to make use of Fourier-based expansions (see also [35]).

## 1.1 Derivation of the scheme

We switch to the twisted variable $v = e^{-it\Delta}u$. We observe that $v$ satisfies

$$i\partial_t v = e^{-it\Delta}(|e^{it\Delta}v|^2e^{it\Delta}v), \quad v_0 = u_0.$$
Equivalently, by integrating the above and mapping Duhamel's formula in Fourier space we have

$$v(t_{n+1}) = v(t_n) - i \sum_{k=-k_1+k_2+k_3} e^{ikx} e^{it_n(k^2+k_1^2-k_2^2-k_3^2)} I^\tau, \quad (6)$$

where the oscillatory integral is given by

$$I^\tau = \int_0^\tau e^{i\omega_1 s} h(s) ds, \quad (7)$$

and $h(s) = e^{i\omega_2 s} g(s)$, $\omega_1 = 2k_1^2$, $\omega_2 = k^2 - k_1^2 - k_2^2 - k_3^2 = -2k_1(k_2 + k_3) + k_2k_3$ and

$$g(s) = \bar{v}_{k_1}(t_n + s)v_{k_2}(t_n + s)v_{k_3}(t_n + s).$$

The central question revolves around making a suitable choice of discretization of the oscillatory integral $I^\tau$, with the aim of minimizing the regularity assumption required by this approximation. The underlying idea behind the construction of the previous (asymmetric) low-regularity integrators (or resonance-based schemes) is to choose an approximation of the integral $I^\tau$ which allows for a practical implementation (by not performing exact integration), while optimizing the local error in the sense of regularity. Namely, by recalling that $2k_1^2$ corresponds to second order derivatives in Fourier space while the terms $k_m k_j$ (for $m \neq j$) correspond to product of first order derivatives, the idea is to separate the dominant ($\omega_1$) and lower-order ($\omega_2$) frequencies. The lower-order and non-oscillatory part $h(s)$ is then approximated by a Taylor series expansion centered at $s = 0$,

$$h(s) = h(0) + O(sw_2g),$$

and the dominant part $e^{i\omega_1 s}$ is integrated exactly. This yields the first-order low-regularity scheme (4) with a local error of $O(\tau^2 \partial_x v)$. At low regularity this is more advantageous than classical techniques (such as exponential integrators [25] or splitting methods [32]) which do not embed the dominant frequency interactions into the scheme and obtain a local error of $O(\tau^2 \partial_x^2 v)$.

The key idea behind obtaining the symmetric scheme (1) is to make a different Taylored discretization of the lower-order and non-oscillatory part $h(s)$. Namely, we again integrate exactly and embed the dominant part $e^{i\omega_1 s}$ into the numerical scheme, while this time approximating the non-dominant part in the following symmetric fashion,

$$h(s) \approx h(0)1_{[0,\tau/2]} + h(\tau)1_{(\tau/2,\tau]}, \quad s \in [0, \tau], \quad (8)$$

where $1_A$ is the indicator function on the set A. By plugging this approximation for $h$ into the oscillatory integral (7) yields two terms: an explicit and an implicit one. The explicit term is given by,

$$\int_0^{\tau/2} e^{i\omega_1 s} h(0) ds = \tau \frac{e^{i\omega_1 \tau/2} - 1}{i\omega_1 \tau} \overline{v_{k_1}^n}(t_n) \overline{v_{k_2}^n}(t_n) \overline{v_{k_3}^n}(t_n) = \frac{\tau}{2} \varphi_1(i k_1^2 \tau) \overline{v_{k_1}^n}(t_n) \overline{v_{k_2}^n}(t_n) \overline{v_{k_3}^n}(t_n).$$

Using the definition of the twisted variable, equation (6), and by mapping the above back to physical space yields the explicit nonlinear term $\psi_E^{\tau/2}(u^n)$ in the scheme (5). Similarly, one obtains the nonlinear implicit term in (5) by using the definitions $\omega_1 = 2k_1^2$, $\omega_2 = k^2 - k_1^2 - k_2^2 - k_3^2$ and noticing that

$$h(\tau) = e^{ik^2\tau} \left( \left( e^{-ik_2^2\tau} v_{k_2}^n(t_{n+1}) \right) \left( e^{-ik_3^2\tau} v_{k_3}^n(t_{n+1}) \right) \left( e^{-ik_1^2\tau} \overline{v_{k_1}^n}(t_{n+1}) \right) \right)$$

and

$$\int_{\tau/2}^{\tau} e^{i\omega_1 s} h(\tau) ds = e^{ik^2\tau} \left( \left( e^{-ik_2^2\tau} v_{k_2}^n(t_{n+1}) \right) \left( e^{-ik_3^2\tau} v_{k_3}^n(t_{n+1}) \right) \int_0^{\tau/2} e^{2ik_1^2(\tau-s)} ds \left( e^{-ik_1^2\tau} \overline{v_{k_1}^n}(t_{n+1}) \right) \right) \\ = e^{ik^2\tau} \left( \left( e^{-ik_2^2\tau} v_{k_2}^n(t_{n+1}) \right) \left( e^{-ik_3^2\tau} v_{k_3}^n(t_{n+1}) \right) \frac{\tau}{2} \varphi_1(-i\tau k_1^2) \left( e^{+ik_1^2\tau} \overline{v_{k_1}^n}(t_{n+1}) \right) \right).$$
We note that a general approach to obtain the approximation (8) is to first give a symmetric approximation to the non-oscillatory part $g(s)$ by iterating Duhamel's formula inside $v_k(t_n + s)$ in a symmetric fashion. Namely, $v_k(t_n + s)$ is approximated on $[0, \frac{\tau}{2}]$ by the linear term in the Duhamel formula centered about $s = 0$ (yielding the approximation $g(0)$, see (6)). While on $(\frac{\tau}{2}, \tau]$, $v_k(t_n + s)$ is approximated by the linear term in the Duhamel formula centered about $s = \tau$ (yielding the approximation $g(\tau)$). We then proceed by approximating the lower-order oscillatory part $e^{i\omega_2 s}$ in a symmetric fashion. In order to obtain higher-order symmetric low-regularity approximations, we would iterate inside $g(s)$ both of these Duhamel expansions (centered about $s = 0$ on $[0, \frac{\tau}{2}]$ and about $s = \tau$ on $(\frac{\tau}{2}, \tau]$) up to higher order. The construction of higher-order low-regularity symmetric integrators will be dealt with in future work.

**Remark 1.2** (Third order local error bound). We make an important point related to the third-order local error structure of the scheme (5). Thanks to the symmetry of the scheme (5) we can expect to have second-order convergence under suitable regularity assumptions on the solution. From the above calculations in Fourier one easily observes that we naturally need three additional derivatives in order to obtain a third-order local error bound of the scheme (5). Indeed, by Taylor expanding around the midpoint one observes that the error induced by the discretization (8) of $I^\tau$ requires the boundedness of a term of the form,

$$ \int_{-\tau}^{0} e^{i\omega_1(\frac{s+\tau}{2})} \left( e^{i\omega_2(\frac{s+\tau}{2})} - 1 \right) \mathrm{dsg}(0) \qquad (9) \\ + \int_{0}^{\tau} e^{i\omega_1(\frac{s+\tau}{2})} \left( e^{i\omega_2(\frac{s+\tau}{2})} - e^{i\omega_2\tau} \right) \mathrm{dsg}(0). $$

For the above to yield a third order term one needs to bound a term of order $O(\tau^3\omega_1\omega_2\bar{v}_{k_1}v_{k_2}v_{k_3})$ which corresponds in physical space to asking for three additional derivatives on $\bar{v}$. This is to be compared with classical (symmetric) schemes which usually have a local error of $O(\tau^3\partial_x^4 v)$ (see for example [32] for splitting schemes), and to asymmetric resonance-based schemes which merely asks for $O(\tau^3\partial_x^2 v)$, ([14, 38]).

While we motivated this symmetric low-regularity integrator on a periodic domain, we show that it also allows for a low-regularity approximation on general smooth domains by establishing its convergence at low-regularity (see Section 1.2). Throughout the remainder of this article we will work on general smooth domains $\Omega \subset \mathbb{R}^d$ and make use of semi-group theory to derive our scheme and establish our convergence result on general domains (see Section 4). This differs from the first structure preserving low-regularity integrators [33, 7] which are restricted to periodic boundary conditions.

We now enter the main bulk of this paper, which answers the question of what can be rigorously proven on the $L^2$-convergence of the scheme (5) when set on a general smooth domain.

We state and prove $L^2$-fractional convergence results, from first to second order, both on the torus $\mathbb{T}^d$ and on a smooth bounded domain $\Omega \subset \mathbb{R}^d$, under moderate regularity assumptions on the solution $u$. These are stronger convergence results than the more typical $H^\sigma(\mathbb{T}^d)$ ($\sigma > \frac{d}{2}$) -convergence analysis, which is restricted to an analysis in smooth Sobolev spaces and to periodic boundary conditions. We state our results in the next subsection.

## 1.2 Result

**Theorem 1.3** ($\Omega = \mathbb{T}^d$). Let $T > 0$, $d \le 3$, and $u_0 \in H^\alpha(\mathbb{T}^d)$ with $\alpha \in [1 + \frac{d}{4}, 3]$. Let $u \in C([0, T], H^\alpha(\mathbb{T}^d))$ be the unique solution of (1). Then there exist $\tau_{\min} > 0$ depending on $T$ and on $\|u_0\|_{H^\alpha}$, and $C_T$ a positive
function depending on T and $\sup_{[0,T]} \|u(t)\|_{H^\alpha}$, such that for every time step size $\tau \le \tau_{\min}$ the numerical solution $u^n$ given in equation (5) has the following error bound:

$$ \|u(n\tau) - u^n\|_{L^2} \le C_T(\sup_{[0,T]} \|u(t)\|_{H^\alpha})\tau^{1+\gamma}, \quad 0 \le n\tau \le T, \quad (10) $$

for $\alpha$ and $\gamma \in [0, 1]$ which satisfy

$$ \begin{cases} \alpha > 1 + \frac{d}{2} \text{ and } 0 \le \gamma \le \frac{\alpha-1}{2}, \\ \alpha = 1 + \frac{d}{2} \text{ and } 0 \le \gamma < \frac{d}{4}, \\ \alpha < 1 + \frac{d}{2} \text{ and } 0 \le \gamma \le \alpha - 1 - \frac{d}{4}. \end{cases} \quad (11) $$

We now consider the case where $\Omega$ is a smooth bounded domain. Given that in this case the space $X^s$ (see Section 2) in which the solution belongs depends not only on Sobolev regularity but also on compatibility conditions which the solution must satisfy on the boundary, we divide the statement of our results depending on the compatibility conditions imposed on $u|_{\partial\Omega}$, (and on the order of convergence).

**Theorem 1.4** ($\Omega \subset \mathbb{R}^d$ smooth bounded domain). Let $\Omega$ be a smooth bounded domain of $\mathbb{R}^d$. We consider the NLS equation (1), equipped with homogeneous Dirichlet boundary conditions. Given any $T > 0$, and $d \le 3$, there exists $\tau_{\min} > 0$ depending on T and the norm of the initial data such that we have the following:

1. Given any $u_0 \in (H^{1+d/4} \cap H_0^1)(\Omega)$ we have first-order convergence of the symmetric scheme (5),

$$ \|u^n - u(n\tau)\|_{L^2} \le C_T \tau, $$

for all $\tau \le \tau_{\min}$, and $0 \le n\tau \le T$.

2. More generally, given any $u_0 \in (H^\alpha \cap H_0^1)(\Omega)$ with $\alpha \in [1 + \frac{d}{4}, 2]$ we have the fractional convergence estimates (10) for $\alpha$ and $\gamma$ which satisfies (26). In particular, we have

$$ \|u(n\tau) - u^n\|_{L^2} \le C_T \begin{cases} \tau^{1+\frac{\alpha-1}{2}} & \text{if } 1 + \frac{d}{2} < \alpha \le 2, \\ \tau^{1+\frac{d}{4}-\epsilon} & \text{for } \alpha = 1 + \frac{d}{2}, \\ \tau^{\alpha-d/4}, & \text{for } 1 + \frac{d}{4} \le \alpha < 1 + \frac{d}{2}, \end{cases} $$

for $0 \le n\tau \le T$, $\tau \le \tau_{\min}$, and for any $\epsilon > 0$.

3. By allowing for more compatibility condition on the boundary we have the following second-order convergence result for an initial data $u_0 \in X^3 = \{u \in H^3(\Omega) : u|_{\partial\Omega} = 0, \Delta u|_{\partial\Omega} = 0 \text{ in } L^2(\partial\Omega)\}$,

$$ \|u^n - u(n\tau)\|_{L^2} \le C_T \tau^2, $$

for all $\tau \le \tau_{\min}$, and $0 \le n\tau \le T$.

We start by making a few remarks on Theorem 1.3, set on the torus $\mathbb{T}^d$. Let us first mention that the symmetric low-regularity integrator (5) requires less regularity assumptions than classical symmetric schemes (see [19, 32, 12, 24]). Indeed, for example the authors [19] require $H^2$ solutions to obtain first order convergence of a Lie splitting scheme for NLS, and the author [32] requires $H^4$-solutions for second-order convergence of a Strang splitting method, whereas we require $H^{1+\frac{d}{4}}$ and $H^3$ to obtain first and resp. second order convergence.

We compare this result to previous convergence results of explicit low-regularity integrators for the NLS equation (1), which are not symmetric and hence do not have good structure preservation properties
(see Figures 2, 3). To the best of our knowledge, this is the first fractional convergence results of a low-regularity schemes to be obtained from first to second order. We compare our full first and second-order convergence result with the work of [2], which also obtains first order convergence in $L^2(\mathbb{T}^d)$ for solutions $u(t) \in H^{1+d/4}(\mathbb{T}^d)$. For the second order convergence in $L^2$ of their asymmetric second-order low-regularity integrator the author [2] asks for solutions $u(t) \in H^{2+d/4}$, whereas the symmetric low-regularity integrator (5) requires a bit more regularity, namely $H^3$ solutions. Moreover, convergence of order $\tau^{1+\gamma}$ in $H^r$-norm, $r > d/2$, for $u_0 \in H^{2\gamma+r+1}(\mathbb{T}^d)$ easily follows from the proof of Theorem 1.3. This is to be compared with asymmetric resonance-based schemes which would typically ask for $u_0 \in H^{\gamma+r+1}(\mathbb{T}^d)$, $r > d/2$. See [35, 38] for a first and resp. second order analysis. We refer to Remark 1.2 which discusses the necessity of requiring three additional derivatives on the solution to obtain second-order convergence of the scheme (5).

We finish by comparing our result with the work of [7], which introduced a symmetrized low-regularity integrator for the Schrödinger map (SM), where they relate the SM flow to the NLS equation set on 1-d torus $\mathbb{T}$ via the Hasimoto transform. The analysis of their scheme is however restricted to the 1-d torus, and to first order convergence in smooth Sobolev spaces $H^r(\mathbb{T})$, $r > 1/2$. The results we present here go beyond the more typical $H^r(\mathbb{T})$ error analysis ($r > 1/2$), by pushing down the error analysis to $L^2$ for first and up to second order convergence. Furthermore, we do not restrict ourselves to Fourier-based techniques, and hence to periodic boundary conditions, as is testified by Theorem 1.4. Using the techniques presented in this article, one can also obtain a symmetric low-regularity approximation to the Schrödinger map in a more general setting than [7].

We now comment upon Theorem 1.4. To our knowledge, this is the first convergence result which goes beyond the first-order convergence analysis of a low-regularity integrator when set on a smooth bounded domain $\Omega \subset \mathbb{R}^d$. We refer to [40, Corollary 20] where the authors show first order convergence in $L^2(\Omega)$ of the asymmetric low-regularity scheme (4) while analogously asking for $(H^{1+d/4} \cap H_0^1)(\Omega)$ solutions.

We also compare our result to the work of [24] which introduces a mass and energy conserving variant of the Crank-Nicolson method as its time-discretization. They show first order convergence on a smooth bounded domain $\Omega \subset \mathbb{R}^d$ under -among other assumptions- $u_t \in L^2(0, T; H^2(\Omega))$, and obtain second order convergence under -among other assumptions- $u_{tt} \in L^2(0, T; H^2(\Omega))$, while assuming $u \in C([0, T], H^2(\Omega))$ throughout their analysis. In contrast to the above classical results Theorem 1.4 permits less regularity assumptions on $u(t)$, namely less than $H^2$-solutions for first order, and less than $H^4$-solutions for second order. We note that the analysis presented here works analogously when adding a potential term $uV$ to equation (1), as is considered in [24]. One would need to ask for the same regularity assumption (and boundary conditions) on $V$ as is required on $u$ in the above theorem. This follows exactly as done in [2]. The case of a rougher potential ($V \in L^\infty(\Omega)$) is dealt in the works of [24, 8, 9].

We also mention that for the 1-d NLS equation with Neumann boundary conditions a low-regularity integrator has been introduced by [6], where using harmonic analysis techniques they could prove up to almost first order convergence with $H^1$-data.

## 1.3 Outline of the paper

In Section 2 we set the scene and introduce the spaces and norms, together with crucial nonlinear estimates, which we will work with throughout the error analysis section. In Section 3 we analyze the implicit nature of the scheme; we show that it is well-defined and establish a crucial a priori estimate on the numerical solution. Finally, in Section 4 we prove the fractional global error estimates presented in Theorems 1.3 and
1.4. First, in Section 4.1 the fractional local error bounds are obtained, followed by Section 4.2 where the stability estimate is shown, and from which the convergence results then naturally follow.

## 2 Norms, spaces, and nonlinear estimates

The norm and space used during the error analysis will depend on the domain $\Omega$ and boundary conditions imposed. We will treat the case where $\Omega = \mathbb{T}^d$ with periodic boundary conditions, and the case of homogeneous Dirichlet boundary conditions when placed on a smooth bounded domain of $\mathbb{R}^d$. In the case where $\Omega = \mathbb{T}^d$ the domain of the operator $\mathcal{L} = -\Delta$ is $D(\mathcal{L}) = H^2(\mathbb{T}^d)$, whereas for Dirichlet boundary conditions we have that $D(\mathcal{L}) = (H^2 \cap H_0^1)(\Omega)$. We can define powers of $\mathcal{L}$, $\mathcal{L}^s$, for $s \ge 0$ using the spectral resolution, and define the space $X^s(\Omega) = \mathcal{D}(\mathcal{L}^{s/2})$ as the domain of the operator $\mathcal{L}^{s/2}$, where $X^0(\Omega) = L^2(\Omega)$. We define the norm on $X^s(\Omega)$ by the usual graph norm

$$ \|u\|_s^2 = \|u\|^2 + \|\mathcal{L}^{s/2}u\|^2, \quad s \ge 0, $$

where $\|u\| = \|u\|_{L^2}$ is the $L^2(\Omega)$-norm. We will be interested in characterizing the space $X^s(\Omega)$ depending on the domain $\Omega$ at study.

### 2.1 The case of periodic boundary conditions

In the case of periodic boundary conditions we have that

$$ X^s(\mathbb{T}^d) = H^s(\mathbb{T}^d) := \left\{ u = \sum_{k \in \mathbb{Z}^d} u_k \frac{e^{ikx}}{\sqrt{(2\pi)^d}} \in L^2(\mathbb{T}^d) : |u|_s^2 \triangleq \sum_{k \in \mathbb{Z}^d} |k|^{2s} |u_k|^2 < \infty \right\} $$

with equivalence of norms

$$ \|u\|_s^2 = \|u\|_{L^2(\mathbb{T}^d)}^2 + \|(-\Delta)^{s/2} u\|_{L^2(\mathbb{T}^d)}^2 = \sum_{k \in \mathbb{Z}^d} (1 + |k|^{2s}) |u_k|^2 = \|u\|_{H^s}^2, $$

$$ \text{where } u_k = \frac{1}{\sqrt{(2\pi)^d}} \int_{\mathbb{T}^d} ue^{-ikx} dx. $$

### 2.2 The case of Dirichlet boundary conditions

We will be interested in characterizing the domain $X^s(\Omega)$ for $s \in [0, 2] \cup \mathbb{N}$ (see Theorem 1.4).

In the case where $s = m \in \mathbb{N}$ we have the following characterization (see [42, Lemma 3.1])

$$ X^m = \{u \in H^m(\Omega) : \Delta^j u = 0 \text{ in } L^2(\partial\Omega) \text{ for } j < m/2\}, $$

with equivalence of the norms on $H^m(\Omega)$ and $X^m$ for functions in $X^m$.

To treat the case where $s$ is not an integer we first introduce the following fractional Sobolev-type spaces known as the Sobolev-Slobodetskij, Gagliardo or Aronszajn space. Given any $s > 0$ of the form $s = m + \sigma$, with $m \in \mathbb{N}$ and $\sigma \in (0, 1)$, we define

$$ H^s(\Omega) = \{u \in H^m(\Omega) : D^\alpha u \in H^\sigma(\Omega) \text{ for any } \alpha \text{ s.t. } |\alpha| = m\}, $$

endowed with the norm

$$ \|u\|_{H^s}^2 = \sum_{|\alpha|=0}^{m} \|D^\alpha u\|^2 + \sum_{|\alpha|=m} |D^\alpha u|_{H^\sigma}^2. $$
For $s = m$ an integer the space $H^s(\Omega)$ coincides with the usual Sobolev space $H^m(\Omega)$, and for $\sigma \in (0, 1)$ we have

$$H^\sigma(\Omega) = \left\{ u \in L^2(\Omega) : |u|_{H^\sigma}^2 := \int_{\Omega} \int_{\Omega} \frac{|u(x) - u(y)|^2}{|x-y|^{d+2\sigma}} dxdy < \infty \right\}.$$

We note that all of the fractional Sobolev spaces which we introduce here can also be defined by using interpolation theory. Indeed, the above space is an intermediary Banach space between $L^2(\Omega)$ and $H^1(\Omega)$, and can be defined by interpolation as

$$H^\sigma(\Omega) = [L^2(\Omega), H^1(\Omega)]_\sigma,$$

see [13, Appendix 1] and [46]. Finally, for $s \in (1/2, 2]$ we define

$$H_D^s(\Omega) = \{u \in H^s(\Omega) : u|_{\partial\Omega} = 0 \text{ in } L^2(\partial\Omega)\},$$

it follows from the above that $D(\mathcal{L}) = H_D^2(\Omega)$. We can now express $X^s$ in terms of Sobolev spaces for $s \in [0, 2] \setminus \frac{1}{2}$ (see [46, Theorem 16.12]),

$$X^s(\Omega) = \begin{cases} H^s(\Omega) & \text{if } 0 \le s < \frac{1}{2} \\ H_D^s(\Omega) & \text{if } \frac{1}{2} < s \le 2 \end{cases},$$

with norm equivalence

$$C^{-1} \|u\|_{H^s} \le \|u\|_s \le C \|u\|_{H^s}, \quad u \in X^s, \qquad (12)$$

for some constant $C > 0$. In the special case where $s = 1/2$ we have that $X^{1/2} = H_{00}^{1/2}$ is the intermediate space defined by

$$H_{00}^{1/2}(\Omega) := \{u \in H^{1/2}(\Omega) : |u|_{H_{00}^{1/2}}^2 := \int_{\Omega} \frac{u^2(x)}{\text{dist}(x, \partial\Omega)} dx < \infty\},$$

with equivalence of norms on $X^{1/2}$ as in (12), see [4, Prop 2.2].

### 2.2.1 Bilinear and nonlinear estimates

In this section we introduce bilinear estimates that are fundamental for the global error analysis, which we now motivate. The results we present in this article go beyond the more typical $H^s$ error analysis ($s > d/2$), by pushing down the analysis to $L^2$ and obtaining fractional rates of convergence, from first up to second order. In particular, to obtain these fractional rates when $\gamma < d/4$, we need to work in the low-order Sobolev spaces $H^{2\gamma}$ (see Section 4). In order to obtain sharp low-regularity error estimates in these spaces we call upon three bilinear estimates (see equations (14), (15), and (16) below) which are taylored to require the least regularity assumptions on $u$ when bounding the local error terms (see also Remark 2.1).

Let $\gamma \ge 0$ and $\epsilon > 0$. Throughout the error analysis we will use the following bilinear estimates, depending on the values of $\gamma$. In the regime $\gamma > d/4$ we call upon the classical bilinear estimate

$$\|uv\|_{H^{2\gamma}} \lesssim \|u\|_{H^{2\gamma}} \|v\|_{H^{2\gamma}}, \text{ for } \gamma > \frac{d}{4}, \qquad (13)$$

whereas in the regime $\gamma \in [0, d/4)$ we exploit the following three bilinear estimates,

$$\|uv\| \lesssim \|u\|_{H^{\frac{d}{4}+\gamma}} \|v\|_{H^{\frac{d}{4}-\gamma}} \text{ for } 0 \le \gamma < \frac{d}{4}, \qquad (14)$$

$$\|uv\|_{H^{2\gamma}} \lesssim \|u\|_{H^{\frac{d}{4}+\gamma}} \|v\|_{H^{\frac{d}{4}+\gamma}} \text{ for } 0 \le \gamma < \frac{d}{4}, \qquad (15)$$
and

$$ \|uv\|_{H^{2\gamma}} \lesssim \|u\|_{H^{\frac{d}{2}+\epsilon}} \|v\|_{H^{2\gamma}} \quad \text{for } 0 \le \gamma \le \frac{d}{4}, \qquad (16) $$

for any $\epsilon > 0$. The above estimates are particular cases of [27, Theorem 8.3.1], valid either on $\mathbb{R}^d$ or $\mathbb{T}^d$. Furthermore, for a smooth bounded domain $\Omega \subset \mathbb{R}^d$, Stein's extension theorem ([1, p.154]) guarantees the existence of a total extension operator, bounded both from $L^2(\Omega)$ to $L^2(\mathbb{R}^d)$ and from $H^m(\Omega)$ to $H^m(\mathbb{R}^d)$, for any $m \in \mathbb{N}$. By interpolation, this operator is bounded from $H^s(\Omega)$ to $H^s(\mathbb{R}^d)$ for any $s \le m$ (see [1, p. 208]). The estimates (13), (14), (15) and (16) consequently hold on $\Omega$ by extending $u$ and $v$ to $\mathbb{R}^d$, applying the estimates on their extensions, and restricting their product to $\Omega$.

**Remark 2.1** (Bilinear estimates in low-order Sobolev spaces $H^{2\gamma}$, $\gamma < d/4$). A natural bilinear estimate which is essential for an analysis in the spaces $H^{2\gamma}$, $\gamma < d/4$, is the estimate (16). This estimate is an analogue of the estimate (13) in the smooth case where $\gamma > d/4$. These two estimates allow to start and fall back on the same space $H^{2\gamma}$. However, the estimate (16) requires more regularity on $u$ than on $v$ ($\frac{d}{2} + \epsilon > 2\gamma$, for $\gamma < d/4$), and asks for $H^{\frac{d}{2}+\epsilon}$-regularity on the solution. One can obtain more optimal bounds which require less regularity assumptions by equally distributing the regularity on $u$ and on $v$. Indeed, by assuming that $u$ and $v$ have the same regularity, applying the estimate (15) requires $\frac{d}{4} + \gamma$ additional derivatives, which is better than $\frac{d}{2} + \epsilon$ for $\gamma < d/4$. While the estimate (14) is used when $v$ requires $2\gamma$ derivatives more than $u$, and balances the regularity requirement to again ask for $\frac{d}{4} + \gamma$ on both $u$ and $v$ (see Proposition 4.2).

We now consider the nonlinearity, which we denote by

$$ f(u, \bar{u})(t, x) = -iu^2(t, x)\bar{u}(t, x). \qquad (17) $$

One can easily deduce from the inequalities (13) and (16) together with the equivalence of norms on $X^s$ the following estimates on the nonlinearity (17)

$$ \begin{gathered} \|f(w, \bar{w})\|_s \le c_{s,\sigma} \|w\|_{\sigma}^2 \|w\|_s \le C_{s,\sigma} (\|w\|_{\sigma}) \|w\|_s \\ \|f(v, \bar{v}) - f(w, \bar{w})\|_s \le c_{s,\sigma} \|v - w\|_s \sum_{k=0}^2 \|v\|_{\sigma}^k \|w\|_{\sigma}^{2-k} \le C_{s,\sigma} (\|v\|_{\sigma}, \|w\|_{\sigma}) \|v - w\|_s \end{gathered} \qquad (18) $$

where $\sigma = \frac{d}{2} + \epsilon$, $c_{s,\sigma} > 0$, and $C_{s,\sigma}(\|u\|, \|v\|)$ denotes a generic constant which depends continuously on the bounded arguments $\|u\|$ and $\|v\|$. In the regime $s > \frac{d}{2}$ the above holds with $\sigma = s$.

**Remark 2.2** (An analysis for very rough solutions). The main ingredient throughout the error analysis section of this article rests upon the crucial bilinear estimates given above, and restricts the solution to belong to the Sobolev space $H^s$, $s > d/2$. In order to consider very rough solutions $u \in H^s$, $s \le d/2$ one needs to call upon more refined tools such as discrete Bourgain spaces when working on the torus ([37]), and discrete Strichartz estimates when working on the full space ([36]). This delicate error-analysis is out of scope for this paper.

Lastly, as we are interested in obtaining fractional error estimates we will call upon the following estimate several times throughout the error analysis section. For $\gamma \in [0, 1]$, we have,

$$ \left\| \frac{(e^{it\Delta} - 1)}{(-t\Delta)^{\gamma}} u \right\| \le 2^{1-\gamma} \|u\|. \qquad (19) $$
The above estimate easily follows from the usual bound

$$|e^{ix} - 1| \le 2^{1-\gamma} |x|^{\gamma}, \quad \gamma \in [0, 1],$$

and using the discrete spectral decomposition of the operator $\mathcal{L} = -\Delta$.

We finish this section by stating the definition of a commutator-type term, which is used in order to obtain low-regularity error estimates (see Section 4.1). For $H(v_1, \dots, v_n)$, $n \ge 1$, a function and $L$ a linear operator, we define the commutator-type term $\mathcal{C}[H, L]$ as

$$\mathcal{C}[H, L](v_1, \dots, v_n) = -L(H(v_1, \dots, v_n)) + \sum_{j=1}^{n} \partial_{v_j} H(v_1, \dots, v_n) \cdot Lv_j.$$

In our setting, $H(v_1, v_2) = f(v_1, v_2) = -iv_1^2 v_2$ is the nonlinearity given in (17), $L = i\Delta$, and hence

$$\partial_{v_1} f(v_1, v_2) = -2iv_1 v_2 \quad \text{and} \quad \partial_{v_2} f(v_1, v_2) = -iv_1^2 \qquad (20)$$

and

$$\mathcal{C}[f, i\Delta](v_1, v_2) = -\Delta(v_1^2 v_2) + 2v_1 v_2 \Delta v_1 + v_1^2 \Delta v_2 = -2(|\nabla v_1|^2 v_2 + 2v_1 \nabla v_1 \cdot \nabla v_2). \qquad (21)$$

## 3 The implicit nature of the scheme

In this section we deal with the question of solving the nonlinear equation (5) at a given time step. We also provide an a priori bound on the numerical solution $\varphi^\tau(v)$ in terms of $v$, which is crucial for the convergence analysis. We recall from Figure 4 that the implicit nature of the scheme does not adversely affect the computational cost of the method.

In the following, we fix $v \in X^\sigma$ for some $\sigma > d/2$. We note that $v$ will play the role of the element $u^n$ in the scheme (5). We then introduce the map

$$z \mapsto S(z) = e^{i\tau\Delta}v + \psi_E^{\tau/2}(v) + \psi_I^{\tau/2}(z),$$

and wish to prove that it admits a unique fixed point given by $\varphi^\tau(v)$.

We start by introducing some useful estimates on the map $S$.

**Proposition 3.1.** Given any $\sigma > d/2$ we have

$$\|S(z_1) - S(z_2)\|_{\sigma} \le \tau M(\|z_1\|_{\sigma}, \|z_2\|_{\sigma}) \|z_1 - z_2\|_{\sigma}, \quad \text{and} \quad \|S(e^{i\tau\Delta}z_1) - e^{i\tau\Delta}z_1\|_{\sigma} \le \tau \tilde{M}(\|z_1\|_{\sigma}),$$
where $M(\|z_1\|, \|z_2\|)$ and $\tilde{M}(\|z_1\|)$ denote generic constants which depend continuously on their arguments $\|z_1\|$ and $\|z_2\|$.

**Proof of Proposition 3.1.** The proof follows directly from the definition of the map $S$, the scheme (5) and of the estimate (18). $\square$

The following theorem shows that the implicit scheme (5) is well-defined, and admits an a priori bound.

**Theorem 3.2.** Let $R > 0$ and $\sigma > d/2$. There exists $\tau_R > 0$ such that, for all $\tau \le \tau_R$ and $v \in X^\sigma$ with $\|v\|_{\sigma} \le R$, we have that $\varphi^\tau(v)$ defined in (5) is given by

$$\varphi^\tau(v) \stackrel{H^\sigma}{=} \lim_{j \to +\infty} S^j(e^{i\tau\Delta}v). \qquad (22)$$

Moreover, under the same conditions, we have

$$\|\varphi^\tau(v)\|_{\sigma} \le 2R. \qquad (23)$$
*Proof of Theorem 3.2.* For notational convenience we let $x_j = S^j(x_0)$, with $x_0 = e^{i\tau\Delta}v$. We first show by induction that for sufficiently small $\tau$ we have the bound

$$ \|x_j\|_{\sigma} \le 2R, \quad j \ge 0. \qquad (24) $$

We choose $\tau_R > 0$ such that $\tau_R M(2R, 2R) \le 1/2$ and $\tau_R \tilde{M}(R) \le R/2$, with $M$ and $\tilde{M}$ from Proposition 3.1. We assume that $\|x_j\|_{\sigma} \le 2R, \forall j \le J$. It follows that for $\tau \le \tau_R$ we have,

$$ \begin{aligned} \|x_{J+1} - x_0\|_{\sigma} &\le \sum_{j=1}^{J} \|S(x_j) - S(x_{j-1})\|_{\sigma} + \|S(x_0) - x_0\|_{\sigma} \\ &\le \sum_{j=0}^{J} \left( \prod_{k=1}^{j} \tau M(\|x_k\|_{\sigma}, \|x_{k-1}\|_{\sigma}) \right) \|S(x_0) - x_0\|_{\sigma} \\ &\le \frac{R}{2} \sum_{j=0}^{J} \frac{1}{2^j} \\ &\le R. \end{aligned} $$

By recalling that, by assumption, $\|x_0\| = \|v\| \le R$, we conclude from the above that

$$ \|x_{J+1}\|_{\sigma} \le 2R, $$

and hence by induction bound (24) holds.

It then follows that for all $\tau \le \tau_R$, $(x_j)_{j \in \mathbb{N}}$ is a Cauchy sequence. Indeed, for $m > p$ we have

$$ \|x_m - x_p\|_{\sigma} \le \sum_{j=p}^{m-1} \|S(x_j) - S(x_{j-1})\|_{\sigma} \le \frac{R}{2} \sum_{j=p}^{\infty} \frac{1}{2^j} \xrightarrow{p \to \infty} 0. $$

This implies that the sequence $(x_j)_{j \in \mathbb{N}}$ converges in $X^{\sigma}$ to the unique fixed-point $\varphi^{\tau}(v)$ of $S$, and the characterization (22) follows. Finally, by passing to the limit in (24) we obtain the desired a priori bound (23) on $\varphi^{\tau}(v)$, which concludes the proof. $\square$

# 4 Error Analysis

In this section we will prove the following proposition.

**Proposition 4.1.** Let $T > 0$, and $\gamma \in [0, 1] \setminus \{\frac{d}{4}\}$. Then there exists $\tau_{\min} > 0$ such that for every time step $\tau \le \tau_{\min}$ the numerical solution $u^n$ given in equation (5) has the following error bound:

$$ \|u(n\tau) - u^n\|_{L^2} \le C_T \left( \sup_{[0,T]} \|u(t)\|_{H^\alpha} \right) \tau^{1+\gamma}, \quad 0 \le n\tau \le T, \qquad (25) $$

where $\alpha$ is given by

$$ \left\{ \begin{array}{ll} \alpha = 2\gamma + 1 & \text{if } \frac{d}{4} < \gamma \le 1 \\ \alpha = \gamma + 1 + \frac{d}{4} & \text{if } 0 \le \gamma < \frac{d}{4}, \end{array} \right. \qquad (26) $$

and where $\tau_{\min}$ depends on $T$ and on $\|u_0\|_{H^\alpha}$, and $C_T$ is a positive function of its argument, depending on $T$.
Given a fixed convergence rate this proposition expresses the regularity assumptions needed in order to obtain this rate, while on the other hand given a fixed regularity assumption on the initial data Theorems 1.3 and 1.4 express the convergence rates one can attain with the method (5). We now link these results.

*Proof of Theorem 1.3 and Theorem 1.4.* By writing the convergence rate $\gamma$ in terms of the regularity assumptions needed on the solution, it directly follows from the above proposition that the convergence rate (10) holds for $\alpha$ and $\gamma$ which satisfy:

$$ \begin{cases} 1 + \frac{d}{2} < \alpha \le 3, & \gamma = \frac{\alpha - 1}{2} \\ 1 + \frac{d}{4} \le \alpha < 1 + \frac{d}{2}, & \gamma = \alpha - 1 - \frac{d}{4}. \end{cases} \qquad (27) $$

The proof of Theorem 1.3 then follows by using the fact that if (10) holds for some $\tilde{\alpha} \in [1 + \frac{d}{4}, 3]$, $\tilde{\gamma} \in [0, 1]$, then the error bound also holds for any $\alpha \ge \tilde{\alpha}$, and $\gamma \le \tilde{\gamma}$. In particular, we recover the case $\alpha = 1 + \frac{d}{2}$ in Theorem 1.3 by applying the second line in (27) with $(\tilde{\alpha}, \tilde{\gamma}) = (1 + \frac{d}{2} - \epsilon, \frac{d}{4} - \epsilon)$ to obtain convergence for $0 \le \gamma < \frac{d}{4}$. See Figure 5 which illustrates graphically the convergence result.

<figure><img src="image_6.png" /><figcaption>Figure 5: Illustration of the interplay between the regularity parameter $\alpha$ and the convergence rate parameter $\gamma$ in the convergence result stated in Proposition (4.1) and Theorem 1.3. We plot the regularity assumption $(u(t) \in H^\alpha)$ needed in order to obtain convergence of order $\tau^{1+\gamma}$.</figcaption></figure>

The proof of Theorem 1.4 follows in the same manner, with the added constraint of the boundary conditions. Namely, we require $u \in C([0, T], X^s)$ and $s = 2\gamma + 1$ or $\gamma + 1 + d/4$, where the boundary conditions are imposed in the definition of the space $X^s$. See Section 2 for the definition of the spaces $X^s$. $\square$

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

In order to prove Proposition 4.1 we combine local error bounds together with a stability argument to conclude via a Lady Windermere's fan argument. We start by showing the local error bound of order $\tau^{2+\gamma}$ with the regularity assumptions stated in Proposition 4.1.

## 4.1 Local error analysis

We decompose the local error term as follows,

$$ \begin{aligned} u(t_n + \tau) - \varphi^\tau(u(t_n)) &= \int_0^\tau e^{i(\tau-s)\Delta} f(u(t_n + s), \bar{u}(t_n + s))ds - \psi_E^{\tau/2}(u(t_n)) - \psi_I^{\tau/2}(\varphi^\tau(u(t_n))) \\ &= \mathcal{R}_1(\tau, t_n) + \mathcal{R}_2(\tau, t_n) + \mathcal{R}_3(\tau, t_n), \end{aligned} $$
with

$$ \mathcal{R}_1(\tau, t_n) = \left( \int_0^{\tau/2} e^{i(\tau-s)\Delta} f(e^{is\Delta}u(t_n), e^{-is\Delta}\bar{u}(t_n))ds - \psi_E^{\tau/2}(u(t_n)) \right) + \left( \int_{\tau/2}^{\tau} e^{i(\tau-s)\Delta} f(e^{is\Delta}u(t_n), e^{-is\Delta}\bar{u}(t_n))ds - \psi_I^{\tau/2}(e^{i\tau\Delta}u(t_n)) \right), \quad (28) $$

$$ \mathcal{R}_2(\tau, t_n) = \psi_I^{\tau/2}(e^{i\tau\Delta}u(t_n)) - \psi_I^{\tau/2}(\varphi_{\tau}(u(t_n))), \quad (29) $$

and

$$ \mathcal{R}_3(\tau, t_n) = \int_0^\tau e^{i(\tau-s)\Delta} \left( f(u(t_n + s), \bar{u}(t_n + s)) - f(e^{is\Delta}u(t_n), e^{-is\Delta}\bar{u}(t_n)) \right) ds. $$

We start by estimating the first error term $\mathcal{R}_1(\tau, t_n)$. This term is the one which asks for the most regularity, and hence dictates the regularity assumptions required on the solution, and thereby the initial data. We then proceed by estimating each of the terms $\mathcal{R}_2(\tau, t_n)$ and $\mathcal{R}_3(\tau, t_n)$, to obtain a cancelation in their sum thanks to the symmetry of the scheme, yielding the desired local error estimate.

**Proposition 4.2.** The error term $\mathcal{R}_1(\tau, t_n)$ satisfies the following bound,

$$ \|\mathcal{R}_1(\tau, t_n)\| \le \begin{cases} C_T \left( \sup_{t \in [0, T]} \|u(t)\|_{X^{2\gamma+1}} \right) \tau^{2+\gamma} & \text{if } \gamma > \frac{d}{4} \\ C_T \left( \sup_{t \in [0, T]} \|u(t)\|_{X^{\gamma+1+\frac{d}{4}}} \right) \tau^{2+\gamma} & \text{if } 0 \le \gamma < \frac{d}{4} \end{cases}, \quad (30) $$

for $0 \le t_n \le T$.

*Proof of Proposition 4.2.* We define the filtered function as

$$ \begin{aligned} \mathcal{N}(\tau, s, \zeta, \Delta, v) &= -ie^{i(\tau-s)\Delta}[(e^{is\Delta}v)^2(e^{is\Delta}e^{-2i\zeta\Delta}\bar{v})] \\ &= e^{i(\tau-s)\Delta}f(e^{is\Delta}v, e^{i(s-2\zeta)\Delta}\bar{v}) \end{aligned} \quad (31) $$

which plays a fundamental role in the derivation and analysis of our scheme on general domains. In the above expression we duplicate the time variable into $s$ and $\zeta$, pulling out a factor $e^{is\Delta}$ in front of the conjugate term $e^{i(s-2\zeta)\Delta}\bar{v}$. Taylor expanding in the variable $s$ yields the right cancellation with the factor $e^{i(\tau-s)\Delta}$ to recover, after integrating in the variable $\zeta$, the explicit term $\psi_E^{\tau/2}(v)$ in the scheme (5), as is detailed below. Similar filtering techniques are used in [40, 2, 3].

Let $v = u(t_n)$. By Taylor expanding the filtering function (31) around $s = 0$, we obtain that the first term in (28) satisfies

$$ \begin{aligned} \int_0^{\tau/2} e^{i(\tau-\zeta)\Delta} f(e^{i\zeta\Delta}v, e^{-i\zeta\Delta}\bar{v}) d\zeta &= \int_0^{\tau/2} \mathcal{N}(\tau, \zeta, \zeta, \Delta, v) d\zeta \\ &= \int_0^{\tau/2} \mathcal{N}(\tau, 0, \zeta, \Delta, v) d\zeta + \int_0^{\tau/2} \int_0^{\zeta} \partial_s \mathcal{N}(\tau, s, \zeta, \Delta, v) ds d\zeta \\ &= \psi_E^{\tau/2}(v) + \int_0^{\tau/2} \int_0^{\zeta} e^{i(\tau-s)\Delta} C[f, i\Delta](e^{is\Delta}v, e^{i(s-2\zeta)\Delta}\bar{v}) ds d\zeta, \end{aligned} \quad (32) $$
where to obtain the last line we used the definition of the $\varphi_1$ function (see (4)) to obtain that

$$ \int_{0}^{\tau/2} e^{-2i\zeta\Delta} d\zeta = \frac{\tau}{2} \frac{e^{-i\tau\Delta} - 1}{-i\tau\Delta} = \frac{\tau}{2}\varphi_1(-i\tau\Delta) $$

and the definition of commutator term $C[f, i\Delta](u, v)$ given in (21).

Similarly, using the filtering function (31), we can treat the second line in (28) by Taylor expanding around $s = \tau$. This yields

$$ \begin{aligned} \int_{\tau/2}^{\tau} e^{i(\tau-\zeta)\Delta} f(e^{i\zeta\Delta}v, e^{-i\zeta\Delta}\bar{v})d\zeta &= \int_{\tau/2}^{\tau} \mathcal{N}(\tau, \zeta, \zeta, \Delta, v)d\zeta \\ &= \int_{\tau/2}^{\tau} \mathcal{N}(\tau, \tau, \zeta, \Delta, v)d\zeta - \int_{\tau/2}^{\tau} \int_{\zeta}^{\tau} \partial_s \mathcal{N}(\tau, s, \zeta, \Delta, v)dsd\zeta \quad (33) \\ &= \psi_I^{\tau/2}(e^{i\tau\Delta}v) - \int_{\tau/2}^{\tau} \int_{\zeta}^{\tau} e^{i(\tau-s)\Delta}C[f, i\Delta](e^{is\Delta}v, e^{i(s-2\zeta)\Delta}\bar{v})dsd\zeta, \end{aligned} $$

where to go from the second line in the above to the third we used that

$$ \begin{aligned} \int_{\tau/2}^{\tau} \mathcal{N}(\tau, \tau, \zeta, \Delta, v)d\zeta &= -i(e^{i\tau\Delta}v)^2 \left( \left( \int_{\tau/2}^{\tau} e^{-2i\zeta\Delta}d\zeta \right) e^{i\tau\Delta\bar{v}} \right) \\ &= -i(e^{i\tau\Delta}v)^2 \left( \left( \int_{0}^{\tau/2} e^{-2i(\tau-\zeta)\Delta}d\zeta \right) e^{i\tau\Delta\bar{v}} \right) \\ &= -i\frac{\tau}{2}(e^{i\tau\Delta}v)^2\varphi_1(i\tau\Delta)(e^{-i\tau\Delta\bar{v}}), \\ &= \psi_I^{\tau/2}(e^{i\tau\Delta}v). \end{aligned} $$

By definition of $\mathcal{R}_1(\tau, t_n)$ and by using (32) and (33) we have that

$$ \begin{aligned} \mathcal{R}_1(\tau, t_n) &= \int_0^{\tau/2} \int_0^{\zeta} e^{i(\tau-s)\Delta} C[f, i\Delta](e^{is\Delta}v, e^{is\Delta}e^{-2i\zeta\Delta}\bar{v})dsd\zeta \\ &\quad - \int_{\tau/2}^{\tau} \int_{\zeta}^{\tau} e^{i(\tau-s)\Delta} C[f, i\Delta](e^{is\Delta}v, e^{is\Delta}e^{-2i\zeta\Delta}\bar{v})dsd\zeta \quad (34) \\ &= \int_0^{\tau/2} \int_0^{\zeta} e^{i(\tau-s)\Delta} C[f, i\Delta](e^{is\Delta}v, e^{is\Delta}e^{-2i\zeta\Delta}\bar{v}) \\ &\quad - e^{is\Delta} C[f, i\Delta](e^{i(\tau-s)\Delta}v, e^{i(\tau-s)\Delta}e^{-2i(\tau-\zeta)\Delta}\bar{v})dsd\zeta. \end{aligned} $$

Note that the local error structure lead by the commutator-type terms in the above expression requires less regularity assumptions than what is required by classical methods, such as exponential integrators of splitting methods (see [25, 32]). Indeed, from the explicit form (21) of the commutator term we see that this error term requires only one additional derivative on the initial datum rather than two (see also [2, 3]).

We now show that thanks to the symmetry of the scheme, we obtain a cancelation in the second-order error term (34) yielding (up to) a third-order remainder. For notational convenience we let $w_1(s) = e^{is\Delta}v$, $w_2(s, \zeta) = e^{i(s-2\zeta)\Delta}\bar{v}$, $z_1(s) = e^{i(\tau-s)\Delta}v$, and $z_2(s, \zeta) = e^{i(2\zeta-\tau-s)\Delta}\bar{v}$ and we denote the integrand by

$$ R_r(\tau, s, \zeta, v) = e^{i(\tau-s)\Delta}C[f, i\Delta](w_1(s), w_2(s, \zeta)) - e^{is\Delta}C[f, i\Delta](z_1(s), z_2(s, \zeta)). \quad (35) $$

It follows from the above equations (34) and (35) that

$$ \| \mathcal{R}_1(\tau, t_n) \| \le \frac{\tau^2}{4} \sup_{s, \zeta \in [0, \tau/2]} \| R_r(\tau, s, \zeta, v) \|. $$
It remains to show that

$$ \sup_{s, \zeta \in [0, \tau/2]} \|R_r(\tau, s, \zeta, v)\| \le \begin{cases} C(\|v\|_{2\gamma+1})\tau^\gamma & \text{if } \gamma > d/4 \\ C(\|v\|_{\gamma+1+\frac{d}{4}})\tau^\gamma & \text{if } \gamma < d/4 \end{cases}. \quad (36) $$

We first approximate the exponentials appearing in front of both commutator terms in (35) to obtain the following first approximation result on $R_r$.

**Lemma 4.3.** We have

$$ R_r(\tau, s, \zeta, v) = R_r^1(\tau, s, \zeta, v) + R_r^2(\tau, s, \zeta, v), $$

with

$$ R_r^1(\tau, s, \zeta, v) = (e^{i(\tau-s)\Delta} - 1)C[f, i\Delta](w_1(s), w_2(s, \zeta)) - (e^{is\Delta} - 1)C[f, i\Delta](z_1(s), z_2(s, \zeta)) $$

and

$$ R_r^2(\tau, s, \zeta, v) = C[f, i\Delta](w_1(s), w_2(s, \zeta)) - C[f, i\Delta](z_1(s), z_2(s, \zeta)), $$

which satisfy

$$ \sup_{s, \zeta \in [0, \tau/2]} \|R_r^i(\tau, s, \zeta, \tau, v)\| \le \begin{cases} C(\|v\|_{2\gamma+1})\tau^\gamma & \text{if } \gamma > d/4, \\ C(\|v\|_{\gamma+1+\frac{d}{4}})\tau^\gamma & \text{if } \gamma < d/4, \end{cases} \quad (37) $$

for $i = 1, 2$.

**Proof of Lemma 4.3.** We write the error term $R_r^1$ as

$$ \begin{aligned} R_r^1(s, \zeta, \tau, v, \bar{v}) = & (\tau - s)^\gamma \left( \frac{e^{i(\tau-s)\Delta} - 1}{(\tau - s)^\gamma (-\Delta)^\gamma} \right) (-\Delta)^\gamma C[f, i\Delta](w_1(s), w_2(s, \zeta)) \\ & - s^\gamma \left( \frac{e^{is\Delta} - 1}{s^\gamma (-\Delta)^\gamma} \right) (-\Delta)^\gamma C[f, i\Delta](z_1(s), z_2(s, \zeta)). \end{aligned} \quad (38) $$

Using the bound given in equation (19), and the boundedness of $e^{it\Delta}$ on Sobolev spaces, it follows from equation (38) that we are left to provide a bound on $C[f, i\Delta](v, \bar{v})$ of the form

$$ \|C[f, i\Delta](v, \bar{v})\|_{2\gamma} \le \begin{cases} C(\|v\|_{2\gamma+1}) & \text{if } \gamma > d/4 \\ C(\|v\|_{\gamma+1+\frac{d}{4}}) & \text{if } \gamma < d/4 \end{cases}. \quad (39) $$

From the definition of the commutator (21), and by using the equivalence of norms $\| \cdot \|_s$ and $\| \cdot \|_{H^s}$ on $X^s$ (see Section 2), it follows that for $d/4 < \gamma \le 1$

$$ \begin{aligned} \|\nabla v\|^2 \bar{v} \|_{2\gamma} & \lesssim \|\nabla v\|^2 \bar{v} \|_{H^{2\gamma}} \\ & \lesssim \|v\|_{H^{2\gamma}} \|\nabla v \cdot \nabla v\|_{H^{2\gamma}} \\ & \lesssim \|v\|_{H^{2\gamma}} \|v\|_{H^{2\gamma+1}}^2 \\ & \lesssim C(\|v\|_{2\gamma+1}), \end{aligned} $$

where we used the estimate (13). Similarly, in the case $0 \le \gamma < d/4$ we have

$$ \begin{aligned} \|\nabla v\|^2 \bar{v} \|_{2\gamma} & \lesssim \|v\|_{H^{d/2+\epsilon}} \|\nabla v \cdot \nabla v\|_{H^{2\gamma}} \\ & \lesssim \|v\|_{H^{d/2+\epsilon}} \|v\|_{H^{\gamma+1+\frac{d}{4}}}^2 \\ & \lesssim C(\|v\|_{\gamma+1+\frac{d}{4}}), \end{aligned} $$
where to obtain the first line we used estimate (16), to go from the first to the second line we used the estimate (15) and concluded using the fact that $d/2 + \epsilon < \gamma + 1 + \frac{d}{4}$ together with the equivalence of norms on $X^{\gamma+1+\frac{d}{4}}$.

We can bound the second term in the commutator-type term (21) in the same manner to obtain the desired bound (39). The estimate (37) on $R_r^1$ follows immediately.

We now deal with the approximation of the second error term $R_r^2$ by approximating each of the exponentials appearing in the arguments of the commutator terms, namely on $w_1, w_2, z_1$ and $z_2$. Given the form of the commutator (21) and of $w_1, w_2, z_1, z_2$, each term to approximate will either be of the form

$$ \|w \nabla ((e^{i\xi\Delta} - 1)u) \cdot \nabla z\| \qquad (40) $$

or

$$ \|(e^{i\xi\Delta} - 1)w) \nabla u \cdot \nabla z\| \qquad (41) $$

for $\xi \in [0, \tau]$, and where we use the boundedness of $e^{it\Delta}$ on Sobolev spaces ($t \in \mathbb{R}$). We can approximate (40) as follows, given any $\epsilon > 0$,

$$ \begin{aligned} \|w \nabla ((e^{i\xi\Delta} - 1)u) \cdot \nabla z\| & \lesssim \xi^\gamma \left\{ \begin{array}{ll} \|w\|_{H^{2\gamma}} \|\nabla(-\Delta)^\gamma \left(\frac{e^{i\xi\Delta}-1}{(-\xi\Delta)^\gamma}u\right)\| \|\nabla z\|_{H^{2\gamma}} & \text{if } \gamma > d/4 \\ \|w\|_{H^{\frac{d}{2}+\epsilon}} \|\nabla(-\Delta)^\gamma \left(\frac{e^{i\xi\Delta}-1}{(-\xi\Delta)^\gamma}u\right)\|_{H^{\frac{d}{4}-\gamma}} \|\nabla z\|_{H^{\frac{d}{4}+\gamma}} & \text{if } \gamma < d/4 \end{array} \right. \\ & \le \xi^\gamma \left\{ \begin{array}{ll} C(\|w\|_{2\gamma}, \|u\|_{2\gamma+1}, \|z\|_{2\gamma+1}) & \text{if } \gamma > d/4 \\ C(\|w\|_{\frac{d}{2}+\epsilon}, \|u\|_{\gamma+1+\frac{d}{4}}, \|z\|_{\gamma+1+\frac{d}{4}}) & \text{if } \gamma < d/4 \end{array} \right. , \end{aligned} $$

where we used the Sobolev embedding $H^\sigma \hookrightarrow L^\infty$, for $\sigma > \frac{d}{2}$, and the estimate (14) to obtain the first inequality. To obtain the second line in the above we used the equivalence of norms, thanks to the fact that $u, v, w$ belong to $X^{2\gamma+1}$ or $X^{\gamma+1+d/4}$ respectively, as well as the estimate (19). Hence, given that $\frac{d}{2} + \epsilon < \gamma + 1 + \frac{d}{4}$, the term above satisfies the desired bound of the form (37).

Furthermore, for the expression (41) we have,

$$ \begin{aligned} \|((e^{i\xi\Delta} - 1)w) \nabla u \cdot \nabla z\| & \lesssim \xi^\gamma \left\{ \begin{array}{ll} \|(-\Delta)^\gamma \frac{e^{i\xi\Delta}-1}{(-\xi\Delta)^\gamma} w\| \|\nabla u \cdot \nabla z\|_{H^{2\gamma}} & \text{if } \gamma > d/4 \\ \|(-\Delta)^\gamma \frac{e^{i\xi\Delta}-1}{(-\xi\Delta)^\gamma} w\|_{H^{1+\frac{d}{4}-\gamma}} \|\nabla u \cdot \nabla z\|_{H^{2\gamma}} & \text{if } \gamma < d/4 \end{array} \right. \\ & \lesssim \xi^\gamma \left\{ \begin{array}{ll} \|(-\Delta)^\gamma \frac{e^{i\xi\Delta}-1}{(-\xi\Delta)^\gamma} w\| \|\nabla u\|_{H^{2\gamma}} \|\nabla z\|_{H^{2\gamma}} & \text{if } \gamma > d/4 \\ \|(-\Delta)^\gamma \frac{e^{i\xi\Delta}-1}{(-\xi\Delta)^\gamma} w\|_{H^{1+\frac{d}{4}-\gamma}} \|\nabla u\|_{H^{\frac{d}{4}+\gamma}} \|\nabla z\|_{H^{\frac{d}{4}+\gamma}} & \text{if } \gamma < d/4 \end{array} \right. \\ & \le \xi^\gamma \left\{ \begin{array}{ll} C(\|w\|_{2\gamma}, \|u\|_{2\gamma+1}, \|z\|_{2\gamma+1}) & \text{if } \gamma > d/4 \\ C(\|w\|_{\gamma+1+\frac{d}{4}}, \|u\|_{\gamma+1+\frac{d}{4}}, \|z\|_{\gamma+1+\frac{d}{4}}) & \text{if } \gamma < d/4, \end{array} \right. \end{aligned} $$

where to obtain the first line we used the Sobolev embedding $H^\sigma \hookrightarrow L^\infty$, $\sigma > \frac{d}{2}$, and the estimate $\|uz\| \lesssim \|u\|_{H^{1+\frac{d}{4}-\gamma}} \|z\|_{H^{2\gamma}}$, $\gamma \in [0, 1]$, (see [27, Theorem 8.3.1]). In order to obtain the second line in the above we again used the estimates (13) and (15), and to obtain the third line we used the equivalence of norms on the spaces $X^s$ together with the estimate (19).

By approximating each of the exponentials in the commutator terms defining $R_r^2$ and by collecting the error terms which are either of the form (40) or (41) we recuperate the desired $\tau^\gamma$ bound in (37).

□
We can conclude from Lemma 4.3 that we have the desired bound (36), and hence the bound (30) on $\mathcal{R}_1(\tau, t_n)$. $\square$

We now continue with the bound on the two remaining terms $(\mathcal{R}_2 + \mathcal{R}_3)(\tau, t_n)$, which asks for less regularity assumptions than the boundedness of $\mathcal{R}_1(\tau, t_n)$, as shown below.

**Proposition 4.4.** For $\gamma \in [0, 1]$, we have the following fractional bound,

$$ \|(\mathcal{R}_2 + \mathcal{R}_3)(\tau, t_n)\| \le C \left( \sup_{[0,T]} \|u(t)\|_{2\gamma} \cdot \sup_{[0,T]} \|u(t)\|_{\sigma} \right) \tau^{2+\gamma}, \quad (42) $$

given any $\sigma > \frac{d}{2}$. In particular we have the bound,

$$ \|(\mathcal{R}_2 + \mathcal{R}_3)(\tau, t_n)\| \le \begin{cases} C_{T,\gamma} \left( \sup_{t \in [0,T]} \|u(t)\|_{2\gamma+1} \right) \tau^{2+\gamma}, & \text{if } \gamma > \frac{d}{4} \\ C_{T,\gamma} \left( \sup_{t \in [0,T]} \|u(t)\|_{\gamma+1+\frac{d}{4}} \right) \tau^{2+\gamma}, & \text{if } 0 \le \gamma < \frac{d}{4} \end{cases}, \quad (43) $$

for $0 \le t_n \le T$.

*Proof of Proposition 4.4.* First, we rewrite the error term $\mathcal{R}_3(\tau, t_n)$ by making suitable Taylor expansions on $f$. We start by expanding $u(t_n + \zeta)$ locally up to second order :

$$ u(t_n + \zeta) = e^{i\zeta\Delta}u(t_n) + \zeta f^n + \tilde{R}(\zeta, t_n) \quad (44) $$

where $f^n = f(u(t_n), \bar{u}(t_n))$ and

$$ \tilde{R}(\zeta, t_n) = \int_0^\zeta e^{i(\zeta-s)\Delta} f(u(t_n + s), \bar{u}(t_n + s))ds - \zeta f^n. \quad (45) $$

Using the above expansion for $u$ we rewrite the error term as

$$ \mathcal{R}_3(\tau, t_n) = \int_0^\tau e^{i(\tau-\zeta)\Delta} \left( f \left( e^{i\zeta\Delta}u(t_n) + \zeta f^n + \tilde{R}(\zeta, t_n), e^{-i\zeta\Delta}\bar{u}(t_n) + \zeta \overline{f^n} + \overline{\tilde{R}(\zeta, t_n)} \right) - f \left( e^{i\zeta\Delta}u(t_n), e^{-i\zeta\Delta}\bar{u}(t_n) \right) \right) d\zeta. \quad (46) $$

For notational convenience we let $a_1 := e^{i\zeta\Delta}u(t_n) + \zeta f^n$. By Taylor expanding $f$ around $(a_1, \bar{a}_1)$ and $(e^{i\zeta\Delta}u(t_n), e^{-i\zeta\Delta}\bar{u}(t_n))$ respectively we obtain,

$$ \begin{aligned} f(a_1 + \tilde{R}(\zeta, t_n), \bar{a}_1 + \overline{\tilde{R}(\zeta, t_n)}) &= f(a_1, \bar{a}_1) + E_1(\zeta) \\ f(a_1, \bar{a}_1) &= f(e^{i\zeta\Delta}u(t_n), e^{-i\zeta\Delta}\bar{u}(t_n)) + E_2(\zeta) \end{aligned} \quad (47) $$

where

$$ E_1(\zeta) = \int_0^1 \partial_{v_1} f \left( a_1 + \theta \tilde{R}(\zeta, t_n), \bar{a}_1 + \theta \overline{\tilde{R}(\zeta, t_n)} \right) \cdot \tilde{R}(\zeta, t_n) + \partial_{v_2} f \left( a_1 + \theta \tilde{R}(\zeta, t_n), \bar{a}_1 + \theta \overline{\tilde{R}(\zeta, t_n)} \right) \cdot \overline{\tilde{R}(\zeta, t_n)} d\theta \quad (48) $$

$$ E_2(\zeta) = \zeta \int_0^1 \left[ \partial_{v_1} f \left( e^{i\zeta\Delta}u(t_n) + \theta\zeta f^n, e^{-i\zeta\Delta}\bar{u}(t_n) + \theta\zeta\overline{f^n} \right) \cdot f^n + \partial_{v_2} f \left( e^{i\zeta\Delta}u(t_n) + \theta\zeta f^n, e^{-i\zeta\Delta}\bar{u}(t_n) + \theta\zeta\overline{f^n} \right) \cdot \overline{f^n} \right] d\theta, \quad (49) $$
and $\partial_{v_1}f$ together with $\partial_{v_2}f$ are given in (20). Hence, plugging the above into equation (46) yields,

$$ \begin{aligned} \mathcal{R}_3(\tau, t_n) &= \int_0^\tau e^{i(\tau-\zeta)\Delta} E_1(\zeta) d\zeta + \int_0^\tau e^{i(\tau-\zeta)\Delta} E_2(\zeta) d\zeta \quad (50) \\ &= \mathcal{E}_1(\tau, t_n) + \mathcal{E}_2(\tau, t_n). \end{aligned} $$

We first deal with the term in the decomposition above which is of highest order and hence is the simplest to bound, namely the third order term $\mathcal{E}_1(\tau, t_n)$. In view of obtaining the bound (42) on $(\mathcal{R}_2 + \mathcal{R}_3)(\tau, t_n)$, we first show that $\mathcal{E}_1(\tau, t_n)$ satisfies this bound.

**Lemma 4.5.** We have the following fractional bound on $\mathcal{E}_1(\tau, t_n)$,

$$ \|\mathcal{E}_1(\tau, t_n)\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_{2\gamma})\tau^{2+\gamma}, $$

for any $\gamma \in [0, 1]$, and $\sigma > d/2$.

*Proof.* It follows from (48) that in order to obtain the above bound we need to show that

$$ \|E_1(\zeta)\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_{2\gamma})\zeta^{1+\gamma}. $$

By equation (48) and by using the Sobolev embedding $H^\sigma \hookrightarrow L^\infty$ we have that for all $\sigma > d/2$,

$$ \begin{aligned} \|E_1(\zeta)\| &\le \sup_{\theta \in [0,1]} \left( \|\partial_{v_1} f (a_1 + \theta \tilde{R}(\zeta, t_n), \bar{a}_1 + \theta \overline{\tilde{R}(\zeta, t_n)})\|_{\sigma} \right. \\ &\qquad \left. + \|\partial_{v_2} f (a_1 + \theta \tilde{R}(\zeta, t_n), \bar{a}_1 + \theta \overline{\tilde{R}(\zeta, t_n)})\|_{\sigma} \right) \|\tilde{R}(\zeta, t_n)\| \quad (51) \\ &\le C \left( \sup_{t \in [0,T]} \|u(t)\|_{\sigma}, \sup_{(\zeta,t) \in [0,\tau] \times [0,T]} \|\tilde{R}(\zeta, t)\|_{\sigma} \right) \|\tilde{R}(\zeta, t_n)\|, \end{aligned} $$

where the last inequality follows by using the explicit form of the derivatives (20), the bilinear inequality (13), the first estimate of equation (18), and the fact that $e^{i\zeta\Delta}$ is an isometry on Sobolev spaces. Next, we show that

$$ \sup_{(\zeta,t) \in [0,\tau] \times [0,T]} \|\tilde{R}(\zeta,t)\|_{\sigma} < \infty, \quad \text{and} \quad \|\tilde{R}(\zeta,t_n)\| \le C_T(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_{2\gamma})\zeta^{1+\gamma}. \quad (52) $$

We obtain the first bound by using the first estimate of equation (18) on $f$ with $r = \sigma$,

$$ \sup_{(\zeta,t) \in [0,\tau] \times [0,T]} \|\tilde{R}(\zeta,t)\|_{\sigma} \le \tau C(\sup_{t \in [0,T]} \|u(t)\|_{\sigma}) < +\infty. $$

Next, we obtain the second fractional estimate of equation (52) by making the following decomposition,

$$ \tilde{R}(\zeta, t_n) = \tilde{R}_1(\zeta, t_n) + \tilde{R}_2(\zeta, t_n), \quad (53) $$

with

$$ \tilde{R}_1(\zeta, t_n) = \int_0^\zeta (\zeta - s)^\gamma \frac{(e^{i(\zeta - s)\Delta} - 1)}{(- (\zeta - s)\Delta)^\gamma} (-\Delta)^\gamma f(u(t_n + s), \bar{u}(t_n + s))ds $$

and

$$ \tilde{R}_2(\zeta, t_n) = \int_0^\zeta f(u(t_n + s), \bar{u}(t_n + s))ds - \zeta f^n. $$
Using the fractional bound (19) and the nonlinear estimate (18) we have that $\tilde{R}_1(\zeta, t_n)$ is bounded by

$$ \|\tilde{R}_1(\zeta, t_n)\| \le \zeta^{1+\gamma} C(\sup_{t \in [0,T]} \|u(t)\|_{\sigma}, \sup_{t \in [0,T]} \|u(t)\|_{2\gamma}). $$

Next, by iterating Duhamel's formula in the first term of $\tilde{R}_2(\zeta, t_n)$ we obtain the following expansion for $\tilde{R}_2(\zeta, t_n)$,

$$ \tilde{R}_2(\zeta, t_n) = \tilde{R}_{2,1}(\zeta, t_n) + \tilde{R}_{2,2}(\zeta, t_n), \quad (54) $$

where

$$ \tilde{R}_{2,1}(\zeta, t_n) = \int_0^\zeta f(e^{is\Delta}u(t_n), e^{-is\Delta}\bar{u}(t_n))ds - \zeta f^n, $$

$$ \tilde{R}_{2,2}(\zeta, t_n) = \int_0^\zeta \int_0^1 \partial_{v_1} f \left( e^{is\Delta} u(t_n) + \theta \tilde{R}_{2,2}^r(s), e^{-is\Delta} \bar{u}(t_n) + \theta \overline{\tilde{R}_{2,2}^r(\zeta)} \right) \cdot \tilde{R}_{2,2}^r(s) \\ + \partial_{v_2} f \left( e^{is\Delta} u(t_n) + \theta \tilde{R}_{2,2}^r(s), e^{-is\Delta} \bar{u}(t_n) + \theta \overline{\tilde{R}_{2,2}^r(s)} \right) \cdot \overline{\tilde{R}_{2,2}^r(s)} d\theta ds, $$

and $\tilde{R}_{2,2}^r(s) = \int_0^s e^{i(s-s_1)\Delta} f(u(t_n + s_1), \bar{u}(t_n + s_1))ds_1$.

Using the nonlinear estimate (18), one easily obtains the bound

$$ \|\tilde{R}_{2,2}(\zeta, t_n)\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma})\zeta^2, $$

and hence in particular the $\zeta^{1+\gamma}$ bound for $\gamma \in [0, 1]$.

In order to deal with the first term in the decomposition (54), we Taylor expand the exponentials appearing in $\tilde{R}_{2,1}(\zeta, t_n)$ which yields,

$$ \tilde{R}_{2,1}(\zeta, t_n) = \int_0^\zeta \int_0^1 \partial_{v_1} f \left( u(t_n) + \theta(e^{is\Delta} - 1)u(t_n), \bar{u}(t_n) + \theta(e^{-is\Delta} - 1)\bar{u}(t_n) \right) \cdot \tilde{R}_{2,1}^r(s) \\ + \partial_{v_2} f \left( u(t_n) + \theta(e^{is\Delta} - 1)u(t_n), \bar{u}(t_n) + \theta(e^{-is\Delta} - 1)\bar{u}(t_n) \right) \cdot \overline{\tilde{R}_{2,1}^r(s)} d\theta ds, $$

with $\tilde{R}_{2,1}^r(s) = s^{\gamma} \frac{e^{is\Delta} - 1}{(-s\Delta)^{\gamma}} (-\Delta)^{\gamma} u(t_n)$. Using the fractional estimate (19) we obtain the bound

$$ \|\tilde{R}_{2,1}^r(s)\| \le C_T(\sup_{[0,T]} \|u(t)\|_{2\gamma})s^{\gamma}. $$

Therefore, by using the usual bilinear inequality (13) we achieve the desired bound on $\tilde{R}_{2,1}$;

$$ \|\tilde{R}_{2,1}(\zeta)\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_{2\gamma})\zeta^{1+\gamma}, $$

which concludes the proof of Lemma 4.5.

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

Now that we have dealt with the third order term $\mathcal{E}_1(t_n, \tau)$ in the decomposition (50) of $\mathcal{R}_3(t_n, \tau)$, we are left to consider the term $\mathcal{E}_2(t_n, \tau)$, together with the term $\mathcal{R}_2(t_n, \tau)$ defined at equation (29). First, we rewrite $\mathcal{E}_2(t_n, \tau)$ as a second order term with a third order remainder. The goal being that this second order term cancels with the second order part of the term $\mathcal{R}_2(t_n, \tau)$, thereby only leaving third order remainders.
By using the definition of $\partial_{v_1}f$ and $\partial_{v_2}f$ given in (20) we have that

$$E_2(\zeta) = -i\zeta \int_0^1 [2(e^{i\zeta\Delta}u(t_n) + \theta\zeta f^n)(e^{-i\zeta\Delta}\bar{u}(t_n) + \theta\zeta \overline{f^n})f^n + (e^{i\zeta\Delta}u(t_n) + \theta\zeta f^n)^2 \overline{f^n}]d\theta.$$

We can separate the first order terms in the above with the higher order ones to obtain the following decomposition for $E_2(\zeta)$,

$$E_2(\zeta) = \tilde{E}_2(\zeta) + E_2^r(\zeta),$$

with

$$\tilde{E}_2(\zeta) = -i\zeta (2(e^{i\zeta\Delta}u(t_n))(e^{-i\zeta\Delta}\bar{u}(t_n))f^n + (e^{i\zeta\Delta}u(t_n))^2 \overline{f^n}),$$

and where one can easily show using the nonlinear estimate (18) with $r = \sigma$ that $E_2^r(\zeta)$ has the following bound: $\|E_2^r(\zeta)\| \le C_T(\sup_{[0,T]} \|u(t)\|_\sigma)\zeta^2$.

We let

$$\tilde{\mathcal{E}}_2(\tau, t_n) = \int_0^\tau e^{i(\tau-\zeta)\Delta} \tilde{\mathcal{E}}_2(\zeta) d\zeta \quad \text{and} \quad \tilde{\mathcal{E}}_2^r(\tau, t_n) = \int_0^\tau e^{i(\tau-\zeta)\Delta} \tilde{\mathcal{E}}_2^r(\zeta) d\zeta,$$

where the error term produced has the bound $\|\tilde{\mathcal{E}}_2^r(\tau, t_n)\| \le C(\sup_{[0,T]} \|u(t)\|_\sigma)\tau^3$, which in particular satisfies the $\tau^{2+\gamma}$ bound for $\gamma \in [0, 1]$.

It remains to show that the sum of the remaining terms to bound $(\mathcal{R}_2 + \tilde{\mathcal{E}}_2)(\tau, t_n)$ also satisfy the $\tau^{2+\gamma}$ bound given in equation (42). In view of this, one last approximation step is made on the term $\tilde{\mathcal{E}}_2(\tau, t_n)$ before estimating its sum with the term $\mathcal{R}_2(\tau, t_n)$. By Taylor expanding around $\zeta = \tau$ the function $\zeta \mapsto e^{i(\tau-\zeta)\Delta}(\tilde{\mathcal{E}}_2(\zeta)/\zeta)$ we obtain the following approximation of $\tilde{\mathcal{E}}_2$,

$$\tilde{\mathcal{E}}_2(\tau, t_n) = \tilde{\mathcal{E}}_{2,1}(\tau, t_n) + \tilde{\mathcal{E}}_{2,1}^r(\tau, t_n),$$

with

$$\tilde{\mathcal{E}}_{2,1}(\tau, t_n) = -i \frac{\tau^2}{2} (2(e^{i\tau\Delta}u(t_n))(e^{-i\tau\Delta}\bar{u}(t_n))f^n + (e^{i\tau\Delta}u(t_n))^2 \overline{f^n}), \quad (55)$$

and where $\tilde{\mathcal{E}}_{2,1}^r(\tau)$ satisfies

$$\|\tilde{\mathcal{E}}_{2,1}^r(\tau)\| \le C(\sup_{[0,T]} \|u(t)\|_\sigma, \sup_{[0,T]} \|u(t)\|_{2\gamma})\tau^{2+\gamma}.$$

The above estimate follows from the definition of $\tilde{\mathcal{E}}_2$, the estimate (19), and the expansion

$$e^{i\tau\Delta}v + (e^{i\zeta\Delta} - e^{i\tau\Delta})v = e^{i\tau\Delta}v + \left(\zeta^\gamma \frac{e^{i\zeta\Delta} - 1}{(-\zeta\Delta)^\gamma} + \tau^\gamma \frac{1 - e^{i\tau\Delta}}{(-\tau\Delta)^\gamma}\right) (-\Delta)^\gamma v.$$

The first estimate (42) of Proposition 4.4 then follows directly once the following lemma is established.

**Lemma 4.6.** The remaining error terms have the following bound,

$$\|\mathcal{R}_2(\tau, t_n) + \tilde{\mathcal{E}}_{2,1}(\tau, t_n)\| \le C(\sup_{[0,T]} \|u(t)\|_\sigma, \sup_{[0,T]} \|u(t)\|_{2\gamma})\tau^{2+\gamma}, \quad (56)$$

for any $\sigma > d/2$.

*Proof.* We perform a very similar analysis as was done on the term $\mathcal{R}_3(\tau, t_n)$ to the term $\mathcal{R}_2(\tau, t_n)$ to show that it can be decomposed as a second and third order term. We then conclude by showing that this second order part coincides with the second order term $-\tilde{\mathcal{E}}_{2,1}(\tau, t_n)$.
First, we expand $\varphi^\tau(u(t_n))$ as follows,

$$ \varphi^\tau(u(t_n)) = e^{i\tau\Delta}u(t_n) + \tau f^n + R^1(\tau, t_n), \quad (57) $$

where

$$ R^1(\tau, t_n) = \Psi^\tau(u(t_n), \varphi^\tau(u(t_n))) - \tau f^n, $$

and $\Psi^\tau$ is the nonlinear part of the numerical scheme (5). We show the same bounds in $X^\sigma$ and $L^2$ given in equation (52) on the error term $R^1(\tau, t_n)$. First, using the estimate (18) with $r = \sigma$ we have that

$$ \sup_{t \in [0,T]} \|R^1(\tau, t)\|_{\sigma} \le \tau C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|\varphi^{\tau}(u(t))\|_{\sigma}) \le C(\sup_{[0,T]} \|u(t)\|_{\sigma}) < \infty, $$

where we use the bound

$$ \|\varphi^\tau(u(t))\|_{\sigma} \le 2\|u(t)\|_{\sigma}, \quad (58) $$

which follows from equation (23). We note that in what follows we will always use the above a priori bound (58) when bounding the term $\varphi^\tau(u(t_n))$, and hence will not show its explicit dependence. In order to obtain the $L^2$-bound on $R^1(\tau, t_n)$ we use the following expansion

$$ R^1(\tau, t_n) = -i \frac{\tau}{2} e^{i\tau\Delta} \left( (u(t_n))^2 \varphi_1(-i\tau\Delta) \overline{u(t_n)} \right) - i \frac{\tau}{2} \left( \left( e^{i\tau\Delta} u(t_n) + \Psi^\tau(u(t_n), \varphi^\tau(u(t_n))) \right)^2 \right. \\ \left. \qquad \varphi_1(i\tau\Delta) \left( e^{-i\tau\Delta} \overline{u(t_n)} + \overline{\Psi^\tau(u(t_n), \varphi^\tau(u(t_n))}} \right) \right) - \tau f^n \quad (59) $$

where we used equation (5) and simply inserted into the term $\Psi^\tau(u(t_n), \varphi^\tau(u(t_n)))$ the definition of the scheme

$$ \varphi^\tau(u(t_n)) = e^{i\tau\Delta}u(t_n) + \Psi^\tau(u(t_n), \varphi^\tau(u(t_n))). $$

By expanding about $e^{i\tau\Delta}u(t_n)$ the second term in (59), and then by approximating the remaining exponentials using the usual fractional estimate (19), we obtain that the first-order terms cancel leaving the following bound on $R^1(\tau, t_n)$,

$$ \begin{aligned} \|R^1(\tau, t_n)\| &\le \tau^{1+\gamma} C(\sup_{[0,T]} \|u(t)\|_{2\gamma}, \sup_{[0,T]} \|u(t)\|_{\sigma}) + \tau C(\sup_{n\tau \le T} \|\Psi^\tau(u(t_n), \varphi^\tau(u(t_n)))\|_{\sigma}) \\ &\le \tau^{1+\gamma} C(\sup_{[0,T]} \|u(t)\|_{2\gamma}, \sup_{[0,T]} \|u(t)\|_{\sigma}), \end{aligned} $$

where in order to obtain the second line we used the estimate

$$ \|\Psi^\tau(u(t_n), \varphi^\tau(u(t_n)))\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma})\tau. $$

We make note that in order to approximate the $\varphi_1$ functions appearing in (59) we use the following expansion

$$ \tau \varphi_1(i \tau \Delta) v = \int_0^\tau e^{is\Delta} ds v = \tau v + \int_0^\tau \frac{(e^{is\Delta} - 1)}{(-s\Delta)^\gamma} s^\gamma ds (-\Delta)^\gamma v. \quad (60) $$

We conclude from the above calculations that the error term $R^1(\tau, t_n)$ satisfies the $X^\sigma$ and $L^2$ bounds,

$$ \sup_{t \in [0,T]} \|R^1(\tau, t)\|_{\sigma} < \infty \quad \text{and} \quad \|R^1(\tau, t_n)\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_{2\gamma}) \tau^{1+\gamma}. \quad (61) $$
We now return to the definition (29) of $\mathcal{R}_2(\tau, t_n)$ and to the expansion (57). By letting $b_1 = e^{i\tau\Delta}u(t_n) + \tau f^n$ we have

$$ \begin{aligned} \psi_I^{\tau/2}(b_1 + R^1(\tau, t_n)) &= \psi_I^{\tau/2}(b_1) + \mathcal{E}_3(\tau, t_n) \\ \psi_I^{\tau/2}(b_1) &= \psi_I^{\tau/2}(e^{i\tau\Delta}u(t_n)) + \mathcal{E}_4(\tau, t_n), \end{aligned} \quad (62) $$

where using the estimates in equation (61) and the definition of $\psi_I^{\tau/2}$ it follows that

$$ \|\mathcal{E}_3(\tau, t_n)\| \le \tau^{2+\gamma} C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_{2\gamma}). $$

Furthermore, it follows from equation (62) that by isolating the second order terms with the higher order ones we have the following expansion for $\mathcal{E}_4(\tau, t_n)$,

$$ \begin{aligned} \mathcal{E}_4(\tau, t_n) = &-i\frac{\tau^2}{2} \left( 2f^n(e^{i\tau\Delta}u(t_n))(\varphi_1(i\tau\Delta)(e^{-i\tau\Delta}\overline{u(t_n)})) \right. \\ &\left. \qquad + (e^{i\tau\Delta}u(t_n))^2\varphi_1(i\tau\Delta)\overline{f^n} \right) + \mathcal{E}_4^r(\tau, t_n), \end{aligned} \quad (63) $$

where from a simple calculation one obtains that $\mathcal{E}_4^r(\tau, t_n)$ satisfies

$$ \|\mathcal{E}_4^r(\tau, t_n)\| \le C(\sup_{[0,T]} \|u(t)\|_{\sigma})\tau^3 \le C_T \tau^{2+\gamma}. $$

By approximating the $\varphi_1$ functions in (63) following the expansion given in equation (60), and by using once again the fractional estimate (19) we conclude from the above equations together with definition (55) of $\tilde{\mathcal{E}}_{2,1}$ that the bound (56) is met. This concludes the proof of Lemma 4.6. $\square$

The proof of the above lemma concludes the proof of the first estimate (42) on $(\mathcal{R}_2 + \mathcal{R}_3)(\tau, t_n)$ of Proposition 4.4.

The second estimate (43) of Proposition 4.4 follows directly from the first estimate (42) by noticing that for some small $\epsilon > 0$ and with $\sigma = d/2 + \epsilon$, we have that $2\gamma$ and $\sigma$ are smaller than $2\gamma + 1$ for $\gamma > d/4$, and are also smaller than $\gamma + 1 + d/4$ for $d \le 3$ (and $\gamma \in [0, 1]$). $\square$

**Remark 4.7.** Another way of writing the local error terms is as follows,

$$ u(t_n + \tau) - \varphi^\tau(u(t_n)) = \mathcal{R}(\tau, t_n) + \tilde{\mathcal{R}}(\tau, t_n), \quad (64) $$

with

$$ \begin{aligned} \mathcal{R}(\tau, t_n) = &\int_0^{\tau/2} e^{i(\tau-s)\Delta} \left( f(u(t_n+s), \bar{u}(t_n+s)) - f(e^{is\Delta}u(t_n), e^{-is\Delta}\bar{u}(t_n)) \right) ds \\ &+ \int_{\tau/2}^{\tau} e^{i(\tau-s)\Delta} \left( f(u(t_n+s), \bar{u}(t_n+s)) - f(e^{i(s-\tau)\Delta}u(t_{n+1}), e^{-i(s-\tau)\Delta}\bar{u}(t_{n+1})) \right) ds \end{aligned} $$

and

$$ \begin{aligned} \tilde{\mathcal{R}}(\tau, t_n) = &\left( \int_0^{\tau/2} e^{i(\tau-s)\Delta} f(e^{is\Delta}u(t_n), e^{-is\Delta}\bar{u}(t_n)) ds - \psi_E^{\tau/2}(u(t_n)) \right) \\ &+ \left( \int_{\tau/2}^{\tau} e^{i(\tau-s)\Delta} f(e^{i(s-\tau)\Delta}u(t_{n+1}), e^{-i(s-\tau)\Delta}\bar{u}(t_{n+1})) ds - \psi_I^{\tau/2}(\varphi^\tau(u(t_n))) \right). \end{aligned} $$
The above error decomposition uses the fact that on $[0, \frac{\tau}{2}]$ we center the approximation at the left-end point and on $(\frac{\tau}{2}, \tau]$ at the right-end point. Hence, on each interval respectively we iterate the Duhamel expansions

$$u(t_n + s) = e^{is\Delta}u(t_n) + \int_0^s e^{i(s-s_1)\Delta}f(u(t_n + s_1), \bar{u}(t_n + s_1))ds_1, \quad s \in [0, \frac{\tau}{2}],$$

$$u(t_n + s) = e^{i(s-\tau)\Delta}u(t_{n+1}) - \int_s^\tau e^{i(s-s_1)\Delta}f(u(t_n + s_1), \bar{u}(t_n + s_1))ds_1, \quad s \in (\frac{\tau}{2}, \tau].$$

Using the tools in Section 4.1 one can bound the local error terms (64) in an analogous manner.

## 4.2 Stability

**Theorem 4.8.** Let $R > 0$, $s \ge 0$. There exists $\tau_R > 0$ and $\sigma > d/2$ such that for any $\tau \le \tau_R$ and $w, v \in X^\sigma$, such that $\|w\|_\sigma \le R$ and $\|v\|_\sigma \le R$ we have,

$$\|\varphi^\tau(v) - \varphi^\tau(w)\|_s \le e^{\tau C_R} \|v - w\|_s,$$

where $C_R$ denotes a generic constant depending on $R$ (and on $s$).

**Proof of Theorem 4.8.** Using the second estimate in equation (18) we have

$$\|\varphi^\tau(v) - \varphi^\tau(w)\|_s \le (1 + \tau C_s(\|v\|_\sigma, \|w\|_\sigma)) \|v - w\|_s + \tau C_s(\|\varphi^\tau(v)\|_\sigma, \|\varphi^\tau(w)\|_\sigma) \|\varphi^\tau(v) - \varphi^\tau(w)\|_s,$$

with $\sigma = d/2 + \epsilon$ if $s \le d/2$ and $\sigma = s$ if $s > d/2$. By Theorem 3.2 we have that there exists $\tau_R > 0$ such that for all $\tau \le \tau_R$ we have the bounds: $\|\varphi^\tau(v)\|_\sigma \le 2R$ and $\|\varphi^\tau(w)\|_\sigma \le 2R$. Hence, it follows from the above that,

$$\|\varphi^\tau(v) - \varphi^\tau(w)\|_s \le \frac{(1 + \tau C_s(R, R))}{(1 - \tau C_s(2R, 2R))} \|v - w\|_s \le e^{\tau C_R} \|v - w\|_s,$$

for some $C_R > 0$. $\square$

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

It remains to combine the stability argument presented in Section 4.2 together with the local error bounds of Section 4.1 to prove the global convergence result stated in Proposition 4.1.

**Proof of Proposition 4.1.** We let $e^n = u^n - u(t_n)$, with $e^0 = 0$. First, thanks to Proposition 4.2 and 4.4 we have the local error bound in $L^2$ required for the global convergence analysis, where the regularity requirements on $u$ have been optimized depending on the fractional order of convergence desired. In order to apply the stability bound stated in Theorem 4.8 and to conclude with a Lady Windermere's fan type argument, one needs to show the following uniform bound on the numerical solution,

$$\|u^n\|_\sigma \le M_T, \quad \forall n\tau \le T, \qquad (65)$$

for some $\sigma > d/2$ and $M_T > 0$.

We let $\sigma = d/2 + \epsilon$, for some small $\epsilon > 0$. To obtain the bound (65) we show that there exists $\delta > 0$, a constant $C_{R_n} = C_T(\|u^n\|_\sigma)$ depending on $\|u^n\|_\sigma$ and on $\sup_{[0,T]} \|u(t)\|_\sigma$, and some $\tau_{R_n} > 0$ also depending on $\|u^n\|_\sigma$ such that the following global error bound is met,

$$\|e^{n+1}\|_{\sigma} \le \|\varphi^{\tau}(u(t_n)) - u(t_{n+1})\|_{\sigma} + \|\varphi^{\tau}(u^n) - \varphi^{\tau}(u(t_n))\|_{\sigma} \le C_{T,\gamma}\tau^{1+\delta} + e^{\tau C_{R_n}} \|e^n\|_{\sigma} \quad (66)$$

for all $\tau \le \tau_{R_n}$. One can obtain the second term in the above estimate for $\tau \le \tau_{R_n}$, where $\tau_{R_n}$ depends on $\|u_n\|$ and on $\sup_{[0,T]} \|u(t)\|_\sigma$, by applying Theorem 4.8 with $s = \sigma$. Hence, it remains to obtain the
first term in the above estimate, which corresponds to the local error bound in $X^\sigma$. Namely, by letting $\mathcal{R}(\tau, t_n) = u(t_{n+1}) - \varphi^\tau(u(t_n))$ we show that there exists $\delta > 0$ such that

$$ \|\mathcal{R}(\tau, t_n)\|_{\sigma} \le C_{T,\gamma} \tau^{1+\delta}, \qquad (67) $$

where $C_{T,\gamma}$ is a constant depending on $T$ and on the regularity assumptions on $u$ (which in turn depend on $\gamma$, the fractional order of convergence required). We establish the above local error estimate (67) by using the following interpolation bound,

$$ \|\mathcal{R}(\tau, t_n)\|_{\sigma} \le \begin{cases} \|\mathcal{R}(\tau, t_n)\|^\theta \|\mathcal{R}(\tau, t_n)\|_{2\gamma+1}^{1-\theta} & \text{if } \gamma > \frac{d}{4} \\ \|\mathcal{R}(\tau, t_n)\|^{\tilde{\theta}} \|\mathcal{R}(\tau, t_n)\|_{\gamma+1+\frac{d}{4}}^{1-\tilde{\theta}} & \text{if } 0 \le \gamma < \frac{d}{4} \end{cases}, \qquad (68) $$

where $(\theta, \tilde{\theta}) \in (0, 1)^2$ satisfies $\sigma = (1 - \theta)(2\gamma + 1) = (1 - \tilde{\theta})(\gamma + 1 + d/4)$. We have already established the $L^2$-bound on $\mathcal{R}(\tau, t_n)$ in Section 4.1, which is given by,

$$ \|\mathcal{R}(\tau, t_n)\| \le \tau^{2+\gamma} \begin{cases} C_{T,\gamma} \left( \sup_{t \in [0,T]} \|u(t)\|_{2\gamma+1} \right) & \text{if } \gamma > \frac{d}{4} \\ C_{T,\gamma} \left( \sup_{t \in [0,T]} \|u(t)\|_{\gamma+1+\frac{d}{4}} \right) & \text{if } 0 \le \gamma < \frac{d}{4} \end{cases}. \qquad (69) $$

To obtain the $X^{2\gamma+1}$ and $X^{\gamma+1+\frac{d}{4}}$ bound on $\mathcal{R}(\tau, t_n)$ we simply express the local error using Duhamel's formula and the scheme (5),

$$ \mathcal{R}(\tau, t_n) = \int_0^\tau e^{i(\tau-s)\Delta} f(u(t_n + s), \bar{u}(t_n + s))ds - i\frac{\tau}{2}e^{i\tau\Delta}((u(t_n))^2\varphi_1(-i\tau\Delta)\bar{u}(t_n)) \\ - i\frac{\tau}{2}\left((\varphi^\tau(u(t_n)))^2\varphi_1(i\tau\Delta)\overline{\varphi^\tau(u(t_n))}\right). $$

One can bound each of the above terms separately using the first estimate in equation (18) with $r = 2\gamma + 1$ and $r = \gamma + 1 + \frac{d}{4}$ together with equation (23) (with $R = \sup_{[0,T]} \|u(t)\|_{\sigma}$) to obtain that there exists some $\tilde{\tau}_0 > 0$ depending on $u_0$ and $T$ such that for all $\tau \le \tilde{\tau}_0$

$$ \|\mathcal{R}(\tau, t_n)\|_r \le C(\sup_{[0,T]} \|u(t)\|_{\sigma}, \sup_{[0,T]} \|u(t)\|_r)\tau \le C_{T,r}\tau. \qquad (70) $$

We conclude that the bound (67) follows from equation (68) (with $\delta = (1 + \gamma)\theta$ for $\gamma > d/4$ and $\delta = (1 + \gamma)\tilde{\theta}$ for $\gamma < d/4$), and where the constant $C_{T,\gamma}$ is given in equation (69) by the $L^2$ local error bound. We then proceed by induction on (66) to obtain that there exists a $\tau_0 > 0$ which depends on $T$ and $u_0$ for which the uniform bound (65) is true for all $\tau \le \tau_0$.

Finally, by taking $s = 0$, $\sigma = d/2 + \epsilon$, and $R = \max\{M_T, \sup_{[0,T]} \|u(t)\|_{\sigma}\}$ in Theorem 4.8 yields the existence of a $\tau_R$ which depends on $T$ and $u_0$ such that for all $\tau \le \tau_R$,

$$ \|e^{n+1}\| \le C_{T,\gamma} \tau^{2+\gamma} + e^{\tau C_R} \|e^n\|, \quad n\tau \le T, $$

where $C_{T,\gamma}$ is given in equation (69). The global error bound of Proposition 4.1 follows by iterating the above estimate and taking $\tau_{\min} = \min\{\tau_0, \tau_R\}$, which concludes the proof. $\square$

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 850941). The author would like to express her thanks to Katharina Schratz for her guidance, as well as David Lee, Georg Maierhofer, and Frédéric Rousset for helpful discussions.
References

[1] R. Adams, J.J.F. Fournier Sobolev Spaces, Springer, 2003.

[2] Y. Alama Bronsard, Error analysis of a class of semi-discrete schemes for solving the Gross-Pitaevskii equation at low regularity, J. Comp. App. Math, 114632, July 2022, doi: 10.1016/j.cam.2022.114632.

[3] Y. Alama Bronsard, Y. Bruned, K. Schratz, *Low regularity integrators via decorated trees*, preprint (2022), arxiv.org/abs/2202.01171.

[4] H. Antil, J. Pfefferer, S. Rogovs, *Fractional Operators with Inhomogeneous Boundary Conditions: Analysis, Control, and Discretization*, Communication in Mathematical Sciences, 2017, 16(5): doi:10.4310/CMS.2018.v16.n5.a11.

[5] X. Antoine, W. Bao, C. Besse, *Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations*. Computer Physics Communications, 2013, 184(12), 2621-2633, doi:10.1016/j.cpc.2013.07.012.

[6] G. Bai, B. Li, and Y. Wu, *A constructive low-regularity integrator for the 1d cubic nonlinear Schrödinger equation under the Neumann boundary condition*. IMA J. Numer. Anal. doi: 10.1093/imanum/drac075.

[7] V. Banica, G. Maierhofer, K. Schratz, *Numerical integration of Schrödinger maps via the Hasimoto transform* preprint (2022), arxiv.org/abs/2205.05024.

[8] W. Bao, C. Wang, *Optimal error bounds on the exponential wave integrator for the nonlinear Schrödinger equation with low regularity potential and nonlinearity* preprint (2023), arXiv:2302.09262.

[9] W. Bao, Y. Ma, C. Wang, *Optimal error bounds on time-splitting methods for the nonlinear Schrödinger equation with low regularity potential and nonlinearity* preprint (2023), arxiv.org/abs/2308.15089.

[10] J. Bernier, B. Grébert, *Birkhoff normal forms for Hamiltonian PDEs in their energy space*. Journal de l'École polytechnique Mathématiques, Tome 9 (2022), pp. 681-745. doi:10.5802/jep.193.

[11] C. Besse, *A relaxation scheme for the nonlinear Schrödinger equation*, SIAM J. Numer. Anal. 42(3) (2004), 934-952.

[12] C. Besse, S. Descombes, G. Dujardin, I. Lacroix-Violet, *Energy-preserving methods for nonlinear Schrödinger equations*, IMA J. Numer. Anal., Volume 41, Issue 1, 2021, pages 618-653, doi:org/10.1093/imanum/drz067.

[13] M. Bonforte, Y. Sire, J. L. Vázquez, *Existence, Uniqueness and Asymptotic behaviour for fractional porous medium equations on bounded domains*, Discrete and Continuous Dynamical Systems, 2015, 35(12): 5725-5767. doi: 10.3934/dcds.2015.35.5725.

[14] Y. Bruned, K. Schratz, *Resonance based schemes for dispersive equations via decorated trees*, to appear in Forum of Mathematics, Pi (2022).

[15] J. Cao, B. Li, and Y. Lin. *A new second-order low-regularity integrator for the cubic nonlinear Schrödinger equation*, IMA Journal of Numerical Analysis (2023): drad017.
[16] M. Cabrera Calvo, K. Schratz. *Uniformly accurate low regularity integrators for the Klein-Gordon equation from the classical to non-relativistic limit regime*, SIAM Journal on Numerical Analysis, 60(2), 888-912.

[17] E. Celledoni, D. Cohen, B. Owren, *Symmetric exponential integrators with an application to the cubic Schrödinger equation*, Foundations of Computational Mathematics, 8 (2008), pp. 303–317.

[18] D. Cohen, E. Hairer, C. Lubich, *Conservation of energy, momentum and actions in numerical discretizations of non-linear wave equations*. Numerische Mathematik, 110(2):113–143, 2008.

[19] J. Eilinghoff, R. Schnaubelt, K. Schratz, *Fractional error estimates of splitting schemes for the nonlinear Schrödinger equation*. J. Math. Anal. Appl. 442, 740-760 (2016).

[20] E. Faou, *Geometric Numerical Integration and Schrödinger Equations*. European Math. Soc. Publishing House, Zürich, 2012.

[21] E. Faou, B. Grébert, E. Paturel, *Birkhoff normal form for splitting methods applied to semilinear Hamiltonian PDEs. Part I. Finite-dimensional discretization*. Numer. Math. 114, 429–458 (2010). doi: 10.1007/s00211-009-0258-y.

[22] L. Gauckler, C. Lubich, *Splitting Integrators for Nonlinear Schrödinger Equations Over Long Times*. FoCM, 10: 275–302, 2010.

[23] E. Hairer, C. Lubich, G. Wanner, *Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations*. Second Edition. Springer Berlin (2006).

[24] P. Henning and D. Peterseim, *Crank-Nicolson Galerkin approximations to nonlinear Schrödinger equations with rough potentials*. M3AS Math. Models Methods Appl. Sci. 27(11):2147-2184, (2017).

[25] M. Hochbruck, A. Ostermann, *Exponential integrators*, Acta Numerica, 19, 209-286. doi:10.1017/S0962492910000048 (2010).

[26] M. Hofmanová, K. Schratz, *An oscillatory integrator for the KdV equation*, Numer. Math. 136:1117-1137 (2017).

[27] L. Hörmander, *Lectures on Nonlinear Hyperbolic Differential Equations*, Springer, 2003.

[28] B. Li, S. Ma, K. Schratz. *A semi-implicit low-regularity integrator for Navier-Stokes equations*, to appear in SIAM J. Numer. Anal., arXiv:2107.13427 (2021).

[29] H. Li, C. Su, *Low regularity exponential-type integrators for the “good” Boussinesq equation*, IMA Journal of Numer. Anal., doi.org/10.1093/imanum/drac081 (2023).

[30] B. Li, Y. Wu, *A fully discrete low-regularity integrator for the 1D periodic cubic nonlinear Schrödinger equation*, Numerische Mathematik 149.1 (2021): 151-183.

[31] B. Li, Y. Wu, *An unfiltered low-regularity integrator for the KdV equation with solutions below $H^1$*. arxiv.org/abs/2206.09320.

[32] C. Lubich, *On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations*, Math. Comput., 77, pp. 2141–2153 (2008).
[33] G. Maierhofer, K. Schratz, *Bridging the gap: symplecticity and low regularity on the example of the KdV equation*. preprint (2022), arxiv.org/abs/2205.05024.

[34] C. Ning, Y. Wu, X. Zhao, *An Embedded Exponential-Type Low-Regularity Integrator for mKdV Equation*. SIAM Journal on Numerical Analysis, 60(3), 999-1025 (2022).

[35] A. Ostermann, K. Schratz, *Low regularity exponential-type integrators for semilinear Schrödinger equations*, Found Comput Math 18, 731–755 (2018). https://doi.org/10.1007/s10208-017-9352-1

[36] A. Ostermann, F. Rousset, K. Schratz, *Error estimates of a Fourier integrator for the cubic Schrödinger equation at low regularity*, Found Comput Math 21, 725–765 (2021). https://doi.org/10.1007/s10208-020-09468-7.

[37] A. Ostermann, F. Rousset, K. Schratz, *Fourier integrator for periodic NLS: low regularity estimates via discrete Bourgain spaces*, to appear in J. Eur. Math. Soc. (JEMS). http://arxiv.org/abs/2006.12785

[38] A. Ostermann, Y. Wu, F. Yao, *A second-order low-regularity integrator for the nonlinear Schrödinger equation*, Adv Cont Discr Mod 2022, 23 (2022). https://doi.org/10.1186/s13662-022-03695-8.

[39] F. Rousset, K. Schratz, *Convergence error estimates at low regularity for time discretizations of KdV*, to appear in *Pure and Applied Analysis*. https://arxiv.org/abs/2102.11125.

[40] F. Rousset, K. Schratz, *A general framework of low regularity integrators*, SIAM J. Numer. Anal., 2021, 59:1735-1768. doi: 10.1137/20M1371506.

[41] K. Schratz, Y. Wang, X. Zhao, *Low-regularity integrators for nonlinear Dirac equations*, Math. Comp. 90:189-214 (2021).

[42] V. Thomée, *Galerkin Finite Element Methods for Parabolic Problems*. Springer-Verlag, Heidelberg, 2006.

[43] Y. Wu, F. Yao. *A first-order Fourier integrator for the nonlinear Schrödinger equation on $\mathbb{T}$ without loss of regularity*, Mathematics of Computation 91.335 (2022): 1213-1235.

[44] Y. Wang, X. Zhao, *A symmetric low-regularity integrator for nonlinear Klein-Gordon equation*, Math. Comp. 91, (2022) 2215-2245. https://doi.org/10.1090/mcom/3751.

[45] Y. Wu, X. Zhao, *Embedded exponential-type low-regularity integrators for KdV equation under rough data*, BIT Numer. Math. 62, (2022) 1049–1090. https://doi.org/10.1007/s10543-021-00895-8.

[46] A. Yagi, *Abstract Parabolic Evolution Equations and their Applications*. Springer-Verlag Berlin Heidelberg, 2010.