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

HAL
open science

# Fully discrete backward error analysis for the midpoint rule applied to the nonlinear Schroedinger equation

Erwan Faou, Georg Maierhofer, Katharina Schratz

► To cite this version:

Erwan Faou, Georg Maierhofer, Katharina Schratz. Fully discrete backward error analysis for the midpoint rule applied to the nonlinear Schroedinger equation. 2025. hal-05054444

HAL Id: hal-05054444

https://inria.hal.science/hal-05054444v1

Preprint submitted on 2 May 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# FULLY DISCRETE BACKWARD ERROR ANALYSIS FOR THE MIDPOINT RULE APPLIED TO THE NONLINEAR SCHRÖDINGER EQUATION

ERWAN FAOU, GEORG MAIERHOFER, AND KATHARINA SCHRATZ

**ABSTRACT.** The use of symplectic numerical schemes on Hamiltonian systems is widely known to lead to favorable long-time behaviour. While this phenomenon is thoroughly understood in the context of finite-dimensional Hamiltonian systems, much less is known in the context of Hamiltonian PDEs. In this work we provide the first dimension-independent backward error analysis for a Runge-Kutta-type method, the midpoint rule, which shows the existence of a modified energy for this method when applied to nonlinear Schrödinger equations regardless of the level of spatial discretisation. We use this to establish long-time stability of the numerical flow for the midpoint rule.

## 1. INTRODUCTION

In this work we study the symplectic midpoint rule applied to the nonlinear Schrödinger equation (NLSE)

$$ (1.1) \qquad \begin{cases} \partial_t z = -i\Delta z + i\lambda |z|^{2r} z, & (t, x) \in (0, T) \times \mathbb{R}, \\ z(0, x) = z_0, & x \in \mathbb{R}, \end{cases} $$

where $r \in \mathbb{N}_{\ge 1}$ denotes the degree of the nonlinearity, and $\lambda \in \{\pm 1\}$ determines if the equation is *focusing* ($\lambda = -1$) or *defocusing* ($\lambda = 1$). We will consider a spatial discretization of this equation by finite differences and consider the family of fully discrete schemes depending on the time and space discretization parameters.

The equation (1.1) is a Hamiltonian partial differential equation (PDE) associated with the real energy

$$ (1.2) \qquad \mathcal{H}(u, \bar{u}) = \int_{\mathbb{R}} |\nabla u(x)|^2 dx + \frac{\lambda}{r+1} \int_{\mathbb{R}} |u(x)|^{2r+2} dx $$

which is preserved for all times along smooth solutions of (1.1), and we can write this latter equation under the symplectic form $\partial_t z = i \frac{\partial \mathcal{H}}{\partial \bar{u}}(z, \bar{z})$ (cf. [11, Section III.1] and [21,

<footnote n="2020">2020 Mathematics Subject Classification. 65P10, 35Q55, 65M22.</footnote>

<footnote>Key words and phrases. Midpoint rule, symplectic integrators, nonlinear Schrödinger equation, full discretisation, modified energy.</footnote>

<footnote>All authors gratefully acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 850941). GM additionally gratefully acknowledges funding from the European Union's Horizon Europe research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 101064261. EF is supported by the Simons collaboration on wave turbulence, the ANR project KEN ANR-22-CE40-0016, and benefited from the support of the Centre Henri Lebesgue ANR-11-LABX-0020-0.</footnote>
Section 3.2]). Note that this equation also preserves the $L^2$ norm

$$ (1.3) \qquad \mathcal{N}(u) = \int_{\mathbb{R}} |u(x)|^2 dx. $$

In the case of finite-dimensional Hamiltonian systems the existence of a modified energy corresponding to the midpoint rule and, more generally, symplectic Runge-Kutta methods, is well-known since the work by Benettin & Giorgilli [5], Murua [22] and Tang [25] (cf. also [17, Chapter IX.3]). This means, in the finite-dimensional case, that the discrete values given by the midpoint rule correspond to the evaluation of a continuous function which is the solution of a modified Hamiltonian system. This result is one of the central underpinnings of advantageous properties of symplectic integrators and permits a rigorous understanding of their long-time behaviour. Perhaps somewhat surprisingly results of this form (i.e. the existence of a modified energy and control on the long-time behaviour) are much more limited for symplectic integrators applied to partial differential equations (i.e. infinite-dimensional Hamiltonian systems). While in practise symplectic methods often exhibit good long-time behaviour [11, 8, 20] the aforementioned results for finite-dimensional systems do not translate easily to the infinite-dimensional case, essentially because the presence of unbounded operators means that analytic bounds derived for finite-dimensional cases break down when the spatial discretisation is refined. Recent work has provided some initial results resolving this problem by proving the existence of a modified energy for splitting methods for example in the work of Faou & Grébert [13] and Bambusi et al. [4].

In the present work we provide, for the very first time for a Runge-Kutta method, dimension-independent guarantees of the existence of a modified Hamiltonian applied to a discretisation of the NLSE (1.1). This is achieved by formulating the midpoint rule as a modified implicit-explicit splitting method involving pseudo-differential flows, and thus follows a two-step process: (i) firstly the existence of a suitable modified vector field is shown which leads to the splitting formulation of the midpoint rule; (ii) we use the implicit-explicit splitting decomposition and an approach based on [13, 4] to prove the existence of a modified energy for the full midpoint method.

An important point to notice is that numerical resonances can a priori occur, and that the existence of the modified energy requires the use of a CFL (Courant-Friedrichs-Lewy, [9]) restriction between the temporal and spatial discretisation parameters. This requirement is not surprising as it also appears in the context of splitting methods for the nonlinear Schrödinger equation.

The remainder of this manuscript is structured as follows. In Section 2 we introduce the fully discrete NLSE which we consider for the remainder of this work, as well as useful notation for the presentation of later results. In Section 3, we then formulate the midpoint rule as an implicit-explicit (IMEX) splitting in the spirit of [2, 24, 23], see Propositions 3.1 and 3.3. This is followed in Section 4 by the formal construction and statement of our main result, which is given by Theorem 4.3 and which gives the existence of a modified energy under a CFL (4.13) similar to the one used in [13]. Finally, as an application, in Section 5, we prove the almost-global stability of the numerical scheme for small initial data in the energy space, see Theorem 5.2.
## 2. PROBLEM SETTING AND NOTATION

2.1. **The discrete NLSE.** Let us first describe the spatial discretisation which we apply to the NLSE (1.1) for the purpose of our analysis. We begin by approximating $\Delta f(x) \approx (\delta x)^{-2}(f(x + \delta x) - 2f(x) + f(x - \delta x))$ which, together with a Dirichlet cut-off at $x = \pm(K + 1)\delta x$, $K \in \mathbb{N}$ leads to the following system of ODEs called the discrete NLSE (see [19] for the derivation, applications and references about this model)

$$ (2.1) \qquad \begin{cases} \frac{du_{\ell}}{dt} = i \frac{1}{(\delta x)^2} (-u_{\ell+1} + 2u_{\ell} - u_{\ell-1}) + i\lambda |u_{\ell}|^{2r} u_{\ell}, & -K \le \ell \le K, \\ u_{\pm(K+1)} = 0, \\ u_{\ell}(0) = z_0(\ell \delta x), & -K \le \ell \le K, \end{cases} $$

and we expect $u_{\ell}(t)$ to be an approximation of $z(t, \ell \delta x)$ the exact solution of (1.1) at the grid points $\ell \delta x$. In the present work we are interested in studying this family of spatially discrete problems for approximating the NLSE (1.1) for arbitrary values $K, \delta x$. Note that in practice, fixing the length $X := K \delta x$ results in the Dirichlet problem for (1.1), i.e. transforms the problem from the unbounded domain $x \in \mathbb{R}$ to the same operator with Dirichlet boundary conditions on $[-X, X]$. We will not study the effect of this spatial truncation here (see [4, 6] for qualitative estimates in the case of solitons).

For any fixed value of $K, \delta x$ the previous system corresponds to the following ODE system in $2K + 1$ dimensions

$$ (2.2) \qquad \frac{du}{dt} = iAu + if(u), $$

with $f(u)_{\ell} = \lambda |u_{\ell}|^{2r} u_{\ell}$ and

$$ A = \frac{1}{\delta x^2} \begin{pmatrix} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & & \ddots & \ddots & \ddots & \\ & & & -1 & 2 & -1 \\ & & & & -1 & 2 \end{pmatrix}. $$

This equation turns out to be a Hamiltonian system associated with the energy (note here we think of $u$ as a column vector)

$$ (2.3) \qquad \begin{aligned} H(u, \bar{u}) &= \delta x (\bar{u}^T A u) + F(u) \\ &:= 2\delta x \sum_{\ell=-K}^{K-1} \frac{|u_{\ell+1} - u_{\ell}|^2}{\delta x^2} + \frac{\lambda}{r+1} \delta x \sum_{\ell=-K}^{K} |u_{\ell}|^{2r+2}, \end{aligned} $$

which is a discrete approximation of the continuous energy $\mathcal{H}$. Note, moreover, that we can check directly that the system (2.1) preserves the discrete $L^2$ norm

$$ (2.4) \qquad N(u) = \delta x \sum_{\ell=-K}^{K} |u_{\ell}|^2. $$

In particular, this preservation property ensures the global existence of the solution to the discrete NLSE system.
2.2. **Midpoint rule.** To approximate the solution $u(t) = (u_\ell)_{\ell=-K}^K$ of (2.1), we discretize in time using a time step $h > 0$ and we consider the sequence $u^n = (u_\ell^n)_{\ell=-K}^K$ defined by induction $u^{n+1} := \varphi_h(u^n)$ as the solution of the implicit equation

$$ (2.5) \qquad u^{n+1} = u^n + \frac{ih}{2}A(u^{n+1} + u^n) + ihf\left(\frac{u^{n+1} + u^n}{2}\right). $$

The sequence $(u^n)_{n \ge N}$ is then an approximation of the $u(nh)$ of (2.1). Note that this symplectic scheme preserves the $L^2$ norm $N(u^{n+1}) = N(u^n)$ which is a quadratic invariant of the problem (see [17, Section VI.7]).

Our aim is to establish uniform estimates for the existence of a modified Hamiltonian for (2.5), meaning estimates which are valid for all systems in this family, i.e. independent of both $\delta x$ and $K$. We aim at proving the following result: For any given $N$, the midpoint rule coincide with the flow at time $h$ of a modified Hamiltonian system associated with a Hamiltonian function $H_h^{(N)}$ in the sense that

$$ \varphi_h(u) = \Phi_{H_h^{(N)}}^h(u) + \mathcal{O}(h^{N+1}) \text{ as } h \to 0, $$

where $\Phi_P^t$ denote the flow of the Hamiltonian system of energy $P$. Such a result is not a surprise for general Hamiltonian system discretised with symplectic methods, when the spatial discretisation is fixed. The main goal of this work is to make the previous construction and estimates independent of $K$ and $\delta x$. Such a result exist for splitting methods, see [10, 13, 11] and can be used to prove stability results over long times for small solutions [11], solitary waves and plane waves [4, 12]. But this work is the first one concerning more classical symplectic Runge-Kutta methods, exemplified here by the midpoint rule.

Before introducing the notation and the mathematical framework that will be used in this paper, let us remark that the linear case $\lambda = f = 0$ degenerates to the linear equation

$$ u^{n+1} = R(hA)u^n = R(hA)^n u^0. $$

where $R$ is the stability function of the midpoint rule:

$$ R(hA) = \frac{1 + i \frac{hA}{2}}{1 - i \frac{hA}{2}} = \exp\left(2i \arctan\left(\frac{hA}{2}\right)\right). $$

This operator can be defined in several ways for example as

$$ (2.6) \qquad R(hA) = U^{-1} \exp\left(2i \arctan\left(\frac{hD}{2}\right)\right) U, $$

where the action of the functions is understood to be on each element of the diagonal matrix $D$ and $U$ is a unitary matrix such that $A = U^{-1}DU$. In this case, the backward error analysis is straightforwardly done: $u^n$ coincides with the solution at time $t = nh$ of the modified system

$$ \frac{d}{dt}v = i \frac{2}{h} \arctan \left( \frac{hA}{2} \right) v. $$

ensuring the preservation of a modified energy for all times. This fact was used in [10, 13] to obtain long time energy estimates for splitting methods.
2.3. **Functional setting.** The discrete space of functions is

$$V_{\delta x}(=V_{\delta x,K}) = \{u \in \mathbb{C}^{\mathbb{Z}} | u_j = 0, |j| > K\},$$

where the dependence in $K$ will remain implicit in the notation<fnref n="1" />. This space is equipped with the discrete norm

$$\|\psi\|_{\delta x}^2 = 2\delta x \sum_{j \in \mathbb{Z}} \frac{|\psi_{j+1} - \psi_j|^2}{\delta x^2} + \delta x \sum_{j \in \mathbb{Z}} |\psi_j|^2,$$

which is the norm associated with the real scalar product

$$(2.7) \qquad \langle \psi, \varphi \rangle_{\delta x} := \delta x \operatorname{Re} \left[ \overline{\psi}^T (I + A) \varphi \right], \qquad \|\psi\|_{\delta x} = \delta x \left[ \overline{\psi}^T (I + A) \psi \right].$$

We can prove (see for instance [4]), that this norm is an algebra norm on $V_{\delta x}$, uniformly in $\delta x$, see Lemma 2.1 below.

Following [3], we identify $V_{\delta x}$ with a finite element subspace of $H^1(\mathbb{R}; \mathbb{C})$. More precisely, defining the function $s: \mathbb{R} \to \mathbb{R}$ by

$$(2.8) \qquad s(x) = \begin{cases} 0 & \text{if } |x| > 1, \\ x + 1 & \text{if } -1 \le x \le 0, \\ -x + 1 & \text{if } 0 \le x \le 1, \end{cases}$$

the identification is done through the map $i_{\delta x}: V_{\delta x} \to H^1(\mathbb{R}; \mathbb{C})$ defined by

$$(2.9) \qquad \{\psi_j\}_{j \in \mathbb{Z}} \mapsto (i_{\delta x}\psi)(x) := \sum_{j \in \mathbb{Z}} \psi_j s\left(\frac{x}{h} - j\right),$$

which we can easily check to be a continuous isomorphism between the two normed vector spaces i.e. there exists constant $c > 0$ and $C$ independent of $\delta x$ and $K$ such that for all $v \in V_{\delta x}$,

$$(2.10) \qquad c \|v\|_{\delta x} \le \|i_{\delta x}(v)\|_{H^1} \le C \|v\|_{\delta x}.$$

**Lemma 2.1.** $V_{\delta x}$ with the norm $\| \cdot \|_{\delta x}$ is an algebra, with a constant independent of dimension. In particular, for any $v, w \in V_{\delta x}$,

$$\|v \cdot w\|_{\delta x} \le C \|v\|_{\delta x} \|w\|_{\delta x},$$

where $C$ does not depend on $K$ and $\delta x$ and where $\cdot$ denotes the elementwise product of the two vectors. Occasionally we will drop the notation $\cdot$ when it is clear from context.

<footnote n="1">¹This is consistent with practical applications, where we might take $K = X(\delta x)^{-1}$ with $X$ denoting the size of the Dirichlet cut-off, or the *large box* in which the problem on the real line is embedded. Note the case $X = 2\pi$ with periodic boundary conditions could be also tackled with a similar analysis.</footnote>
Proof. We have

$$
\begin{align*}
\|v \cdot w\|_{\delta x}^2 &= 2\delta x \sum_{j \in \mathbb{Z}} \frac{|v_{j+1}w_{j+1} - v_j w_j|^2}{\delta x^2} + \delta x \sum_{j \in \mathbb{Z}} |v_j|^2 |w_j|^2 \\
&\le 2\delta x \sum_{j \in \mathbb{Z}} \frac{(|w_{j+1}| |v_{j+1} - v_j| + |v_j| |w_{j+1} - w_j|)^2}{\delta x^2} + \delta x \sum_{j \in \mathbb{Z}} |v_j|^2 |w_j|^2 \\
&\le 4\delta x \sum_{j \in \mathbb{Z}} \frac{|w_{j+1}|^2 |v_{j+1} - v_j|^2 + |v_j|^2 |w_{j+1} - w_j|^2}{\delta x^2} + \delta x \sum_{j \in \mathbb{Z}} |v_j|^2 |w_j|^2 \\
&\le 2\|v\|_{\ell^\infty}^2 \|w\|_{\delta x}^2 + 2\|w\|_{\ell^\infty}^2 \|v\|_{\delta x}^2.
\end{align*}
$$

Thus we have, using the definition (2.9) of $i_{\delta x}$, that there is a constant $C > 0$ independent of $v, w, K$ such

$$
\|v \cdot w\|_{\delta x} \le 2 \left( \|i_{\delta x}(v)\|_{L^{\infty}(\mathbb{R})} \|w\|_{\delta x} + \|i_{\delta x}(w)\|_{L^{\infty}(\mathbb{R})} \|v\|_{\delta x} \right).
$$

By Morrey's inequality it follows that for some $C > 0$ independent of $v, w, K$ we have

$$
\|v \cdot w\|_{\delta x} \le C \left( \|i_{\delta x}(v)\|_{H^1(\mathbb{R})} \|w\|_{\delta x} + \|i_{\delta x}(w)\|_{H^1(\mathbb{R})} \|v\|_{\delta x} \right).
$$

Thus the result follows by equivalence of the norms $v \mapsto \|v\|_{\delta x}$ and $v \mapsto \|i_{\delta x}(v)\|_{H^1(\mathbb{R})}$ on $V_{\delta x}$. $\square$

2.4. Hamiltonian formulation. We use the following standard Hamiltonian coordinates $u = \frac{1}{\sqrt{2}}(p+iq)$. This translates to a coordinate-wise identification once we consider the values of $u_j$ at the node $j \in \mathbb{Z}$. With this complex representation we associate the derivatives

$$
\frac{\partial}{\partial u_j} = \frac{1}{\sqrt{2}} \left( \frac{\partial}{\partial p_j} - i \frac{\partial}{\partial q_j} \right) \quad \text{and} \quad \frac{\partial}{\partial \bar{u}_j} = \frac{1}{\sqrt{2}} \left( \frac{\partial}{\partial p_j} + i \frac{\partial}{\partial q_j} \right), \quad j \in \mathbb{Z}.
$$

Any function $H(p, q)$ from $\mathbb{R}^{2K+1} \times \mathbb{R}^{2K+1} \to \mathbb{R}$ can be viewed as a function $H(u)$ defined on $\mathbb{C}^{2K+1}$ and taking real values$^2$. For such a Hamiltonian function $H: \mathbb{C}^{2K+1} \to \mathbb{R}$ we can then consider the vector

$$
\nabla_{\bar{u}} H(u) \equiv \frac{1}{\sqrt{2}} \begin{pmatrix} \partial_p H \\ \partial_q H \end{pmatrix}
$$

which allows us to introduce the vector field associated with a Hamiltonian function:

**Definition 2.2 (Hamiltonian formulation).** The vector field, $X_H : \mathbb{C}^{2K+1} \to \mathbb{C}^{2K+1}$, associated with a Hamiltonian function $H$, is given by

$$
X_H(u) := i \delta x^{-1} \nabla_{\bar{u}} H,
$$

which in $(p, q)$-coordinates corresponds to

$$
\delta x^{-1} \frac{1}{\sqrt{2}} J^{-1} \begin{pmatrix} \partial_p H \\ \partial_q H \end{pmatrix}, \quad J = \begin{pmatrix} 0 & I \\ -I & 0 \end{pmatrix} \quad \text{with} \quad J^2 = -I,
$$

<footnote>$^2$Note that with this identification, $H$ is in fact a function of $u$ and $\bar{u}$ and not a holomorphic function of $u$</footnote>
and I denoting the $(2K + 1) \times (2K + 1)$ identity. The Hamiltonian system associated with the function H is

$$ (2.11) \qquad \frac{du}{dt} = X_H(u) = i\delta x^{-1} \nabla_{\bar{u}} H(u), $$

which is equivalent to the real system

$$ \frac{d}{dt} \begin{pmatrix} p \\ q \end{pmatrix} = \delta x^{-1} J^{-1} \begin{pmatrix} \partial_p H \\ \partial_q H \end{pmatrix} = \delta x^{-1} \begin{pmatrix} -\partial_q H \\ \partial_p H \end{pmatrix}. $$

**Remark 2.3.** The scaling factor $\delta x^{-1}$ in the Hamiltonian formulation is included to make $\delta x^{-1} \nabla_{\bar{u}}$ consistent with a variational derivative in the limit $\delta x \to 0$. For example consider the functional $\mathcal{N}(u) := \int_{\mathbb{R}} |u(x)|^2 dx$, whose discrete analogue in our setting (2.1) is $N(u) = \delta x \sum_{j=-K}^{K} |u_j|^2$. The functional derivative of $\mathcal{N}$ is given by

$$ \delta_{\bar{u}} \mathcal{N} = u, $$

while

$$ (\nabla_{\bar{u}} N)_j = \frac{\partial N}{\partial \bar{u}_j} = (\delta x) u_j, $$

ensuring that the correct scaling in the Hamiltonian formulation is indeed (2.11).

**Definition 2.4.** In this notation we then call a map $\Phi : \mathbb{C}^{2K+1} \to \mathbb{C}^{2K+1}$ symplectic if its Jacobian

$$ M = \begin{pmatrix} \frac{\partial \mathrm{Re} \Phi}{\partial p} & \frac{\partial \mathrm{Re} \Phi}{\partial q} \\ \frac{\partial \mathrm{Im} \Phi}{\partial p} & \frac{\partial \mathrm{Im} \Phi}{\partial q} \end{pmatrix} $$

satisfies

$$ (2.12) \qquad M^T J M = J, \quad J = \begin{pmatrix} 0 & I \\ -I & 0 \end{pmatrix}. $$

**Remark 2.5.** We can check (see [17]) that the midpoint rule (2.5) is symplectic in this sense (provided the solution is well defined as solution of an implicit system, see below). Moreover, the composition of symplectic maps is clearly symplectic.

**Definition 2.6** (Commutator of vector fields). For two real vector fields $X, Y$ from $\mathbb{R}^{2K+1} \times \mathbb{R}^{2K+1}$ to itself, we define the usual commutator $[X, Y]$ as follows

$$ [X, Y] = \sum_{j=-K}^{K} \left( X_j^{(p)} \frac{\partial}{\partial p_j} + X_j^{(q)} \frac{\partial}{\partial q_j} \right) \left( Y_j^{(p)} \frac{\partial}{\partial p_j} + Y_j^{(q)} \frac{\partial}{\partial q_j} \right) \\ - \left( Y_j^{(p)} \frac{\partial}{\partial p_j} + Y_j^{(q)} \frac{\partial}{\partial q_j} \right) \left( X_j^{(p)} \frac{\partial}{\partial p_j} + X_j^{(q)} \frac{\partial}{\partial q_j} \right), $$

where $X_j^{(p)}$ and $X_j^{(q)}$ denote the $p_j$ and $q_j$ components of X respectively.

**Definition 2.7.** The natural Poisson bracket in this formulation is given by

$$ \{F, G\} := \delta x^{-1} \sum_{j=-K}^{K} \frac{\partial F}{\partial p_j} \frac{\partial K}{\partial q_j} - \frac{\partial G}{\partial p_j} \frac{\partial F}{\partial q_j}. $$
**Remark 2.8.** The scaling in the Poisson bracket is important for consistency, as it ensures (as can be easily verified in the $(p, q)$-coordinates), that for two Hamiltonian functions $P, Q$,

$$[X_P, X_Q] = X_{\{P,Q\}}.$$

**Remark 2.9.** We can also check that with the aforementioned scalings we have, as usual, that for any smooth function $g: \mathbb{C}^{2K+1} \to \mathbb{C}$

$$\frac{\mathrm{d}g(u)}{\mathrm{d}t} = \{H, g\}.$$

Note we will in the following often switch between the formulation in $u \in \mathbb{C}^{2K+1}$ and in $(q, p) \in \mathbb{R}^{2K+1} \times \mathbb{R}^{2K+1}$. For this it is helpful to keep the identification $iu \equiv J^{-1}(q, p)$ in mind.

**2.5. Estimating polynomial vector fields.** In order to establish the desired truncation bounds on the modified energy for the midpoint rule we have to introduce a suitable framework for estimating commutators of polynomial vector fields. For this we shall use the following notation introduced in [4, Section 7.2]. Suppose $X$ is a vector field on $V_{\delta x}$ which is a homogeneous polynomial of degree $s$. Then we can associate (by polarization) with $X$ a symmetric mulilinear form $\tilde{X}(\psi_1, \dots, \psi_s)$ such that for all $\psi \in V_{\delta x}$, $X(\psi) = \tilde{X}(\psi, \dots, \psi)$. This symmetric multilinear form is given by

$$(2.13) \qquad \tilde{X}(\psi_1, \dots, \psi_s) := \frac{1}{s!} \sum_{k=1}^{s} \sum_{1 \le j_1 < \dots < j_k \le s} (-1)^{s-k} X(\psi_{j_1} + \dots + \psi_{j_k}).$$

We can then define the following operator norm on such homogeneous polynomials

$$(2.14) \qquad \|X\|_{\delta x} = \sup_{\|\psi_j\|_{\delta x} = 1, j=1,\dots,s} \|\tilde{X}(\psi_1, \dots, \psi_s)\|_{\delta x} = \sup_{\|\psi\|_{\delta x} = 1} \|X(\psi)\|_{\delta x},$$

see for instance [7, Proposition 2]. For notational convenience we introduce the following space.

**Definition 2.10.** We denote by $\mathcal{P}_s$ the space of all polynomial vector fields $X$ of degree no larger than $s$ such that $\|X\|_{\delta x}$ is uniformly bounded in both $K > 0$ and $\delta x > 0$.

The norm (2.14) can then be extended to $\mathcal{P}_s$ by simply defining $\| \cdot \|_{\delta x}$ of a general polynomial to be the sum of the norms applied to the homogeneous components. Note when $s = 1$ the norm reduces to the usual operator norm on the normed vector space $V_{\delta x}$. Using (2.14) we can establish the following commutator estimate:

**Lemma 2.11** (See Lemma 7.6 in [4]). Suppose $X \in \mathcal{P}_{s_1}, Y \in \mathcal{P}_{s_2}$ are two polynomial vector fields. Then $[X, Y] \in \mathcal{P}_{s_1+s_2-1}$ and

$$(2.15) \qquad \|[X, Y]\|_{\delta x} \le (s_1 + s_2)\|X\|_{\delta x} \|Y\|_{\delta x}.$$

*Proof.* The proof of this statement is given in Lemma 7.6 [4]. □
**2.6. Vector fields and Hamiltonian functions.** As usual there is a one-to-one correspondence between Hamiltonian functions and associated vector fields. We have already seen the definition of $X_H$ and in the following lemma we shall construct $H$ from $X$.

**Lemma 2.12.** Let $f: \mathbb{C}^{2K+1} \to \mathbb{C}^{2K+1}$ be a homogeneous polynomial vector field such that the matrix

$$ (2.16) \qquad J \begin{pmatrix} \frac{\partial \mathrm{Re}f}{\partial p} & \frac{\partial \mathrm{Re}f}{\partial q} \\ \frac{\partial \mathrm{Im}f}{\partial p} & \frac{\partial \mathrm{Im}f}{\partial q} \end{pmatrix} = \begin{pmatrix} \frac{\partial \mathrm{Im}f}{\partial p} & \frac{\partial \mathrm{Im}f}{\partial q} \\ -\frac{\partial \mathrm{Re}f}{\partial p} & -\frac{\partial \mathrm{Re}f}{\partial q} \end{pmatrix}, $$

represented here in block-matrix notation, where $p = \mathrm{Re}(u) \in \mathbb{R}^{2K+1}$, $q = \mathrm{Im}(u) \in \mathbb{R}^{2K+1}$, is symmetric. Suppose further that, for given $C > 0$ independent of $K$, $d \in \mathbb{N}_{\ge 1}$, $f$ satisfies the bound

$$ \|f(u)\|_{\delta x} \le C \|u\|_{\delta x}^d, \quad \forall u \in V_{\delta x}. $$

Then there is a real-valued homogeneous polynomial $P: \mathbb{C}^{2K+1} \to \mathbb{R}$ such that for all $u \in \mathbb{C}^{2K+1}$

$$ (2.17) \qquad f(u) = X_H = i\delta x^{-1} \nabla_{\vec{\pi}} P(u), $$

and there is $\tilde{C} > 0$, independent of $K$, such that

$$ (2.18) \qquad |P(u)| \le \tilde{C} \|u\|_{\delta x}^{d+1}, \quad \forall u \in V_{\delta x}. $$

**Remark 2.13.** By slight abuse of notation we shall in the following write

$$ (2.19) \qquad \nabla_{\vec{u}}^T f = \frac{1}{\sqrt{2}} \begin{pmatrix} \frac{\partial \mathrm{Re}f}{\partial p} & \frac{\partial \mathrm{Re}f}{\partial q} \\ \frac{\partial \mathrm{Im}f}{\partial p} & \frac{\partial \mathrm{Im}f}{\partial q} \end{pmatrix}. $$

*Proof of Lemma 2.12.* The proof is based on [17, Lemma VI.2.7]. We define

$$ P(u) := 2\delta x \int_0^1 \mathrm{Im}(\vec{u}^T f(tu)) \mathrm{d}t, $$

which clearly is a *real-valued* homogeneous polynomial. Moreover, we have with $u = \frac{1}{\sqrt{2}}(p + iq)$ and using summation convention over repeated indices

$$ \begin{aligned} (i\delta x^{-1} \nabla_{\vec{u}} P)_j &= i \left( \frac{\partial}{\partial p_j} + i \frac{\partial}{\partial q_j} \right) \int_0^1 p_\ell \mathrm{Im} f_\ell(tu) - q_\ell \mathrm{Re} f_\ell(tu) \mathrm{d}t \\ &= \int_0^1 \left( i \frac{\partial}{\partial p_j} - \frac{\partial}{\partial q_j} \right) (p_\ell \mathrm{Im} f_\ell(tu) - q_\ell \mathrm{Re} f_\ell(tu)) \mathrm{d}t \\ &= \int_0^1 i \mathrm{Im} f_j(tu) + it p_\ell \frac{\partial \mathrm{Im} f_\ell}{\partial p_j} \Big|_{tu} - it q_\ell \frac{\partial \mathrm{Re} f_\ell}{\partial p_j} \Big|_{tu} \mathrm{d}t \\ &\quad + \int_0^1 -tp_\ell \frac{\partial \mathrm{Im} f_\ell}{\partial q_j} \Big|_{tu} + \mathrm{Re} f_j(tu) + tq_\ell \frac{\partial \mathrm{Re} f_\ell}{\partial q_j} \Big|_{tu} \mathrm{d}t. \end{aligned} $$
Using the symmetry of (2.16) this simplifies to

$$ (i\delta x^{-1} \nabla \bar{u} P)_j = \int_0^1 i \operatorname{Im} f_j(tu) + itq_\ell \frac{\partial \operatorname{Im} f_j}{\partial q_\ell} \Big|_{tu} + itp_\ell \frac{\partial \operatorname{Im} f_j}{\partial p_\ell} \Big|_{tu} \, dt \\ + \int_0^1 tq_\ell \frac{\partial \operatorname{Re} f_j}{\partial q_\ell} \Big|_{tu} + \operatorname{Re} f_j(tu) + tp_\ell \frac{\partial \operatorname{Re} f_j}{\partial p_\ell} \Big|_{tu} \, dt \\ = \int_0^1 \frac{d}{dt} (tf_j(tu)) \, dt = f_j(tu). $$

This completes the proof of (2.17). For the bound (2.18) we note that for all $u \in B_{\delta x}(R)$

$$ |P(u)| = 2 \left| \int_0^1 \delta x \operatorname{Im}(\bar{u}^T f(tu)) dt \right| \\ \le 2 \sup_{t \in [0,1]} |\delta x \bar{u}^T f(tu)| \le 2 \|\bar{u}\|_{\delta x} \sup_{t \in [0,1]} \|f(tu)\|_{\delta x} \\ \le 2C \|\bar{u}\|_{\delta x} \sup_{t \in [0,1]} t^d \|u^d\|_{\delta x}, $$

where we used the fact that for $u, v \in \mathbb{C}^{2K+1}$,

$$ |\delta x \bar{u}^T v| \le |\langle u, v \rangle_{\delta x}| \le \|u\|_{\delta x} \|v\|_{\delta x}. $$

Thus the result follows from Lemma 2.1. □

The following stability estimate will also prove useful in Section 4.2.

**Lemma 2.14.** Suppose P is a polynomial of degree k such that $|P(u)| \le C \|u\|_{\delta x}^k$ then

$$ \|X_P(z) - X_P(y)\|_{\delta x} \le C_k C \left( \max_{n=1,...,k-2} (\|y\|_{\delta x}^n, \|z\|_{\delta x}^n) \right) \|z - y\|_{\delta x}. $$

*Proof.* It follows directly from the multilinear estimates (2.14) (see also Proposition 2.7 in [11]). □

### 3. THE MIDPOINT RULE AS A SPLITTING METHOD

Our goal is to resort to tools introduced in [13, 11] to establish dimension-independent error estimates in the energy (i.e. estimates which, beyond a CFL constraint, do not depend on $\delta x$ and $K$). In order for this to work, we need to first establish the following result (where we denote $B_{\delta x}(R) := \{u \in V_{\delta x} \mid \|u\|_{\delta x} \le R\}$):

**Proposition 3.1.** Let $R > 0$ be given. There exists $h_0(R) > 0$ and, for all $h \le h_0$, a symplectic map $u \mapsto \Psi^h(u)$ from $B_{\delta x}(R)$ to $B_{\delta x}(2R)$ such that the midpoint rule can be written in the split form

$$ (3.1) \qquad u^{n+1} = R(hA) \circ \Psi^h(u^n) $$

and moreover, we have for $u \in B_{\delta x}(R)$

$$ (3.2) \qquad \Psi^h(u) = u + \sum_{k \ge 1} h^k \Psi_{h,k}(u) $$
where the $\Psi_{h,k}(u)$ are homogeneous polynomial vector fields such that there exists a constant $C$ such that for all $h \le h_0$ and all $u$ in $B_{\delta x}(R)$

$$ \|\Psi_{h,k}(u)\|_{\delta x} \le C^k \|u\|_{\delta x}^{2rk+1}, $$

where $C$ does not depend on $N = 2K + 1$, the dimension of the ODE, nor on $\delta x$ the mesh size of the space discretisation.

To prove the above splitting formulation, we need a few lemmas.

**Lemma 3.2.** $R(hA)$ is an isometry on $\| \cdot \|_{\delta x}$.

*Proof.* The eigenvalues of $R(hA)$ are all on $i\mathbb{R}$, and it can be jointly diagonalised with the norm $\| \cdot \|_{\delta x}$, cf. (2.7). $\square$

*Proof of Proposition 3.1.* We can write the midpoint rule as

$$ \begin{aligned} u^{n+1} &= R(hA)u^n + \frac{ih}{1-ihA/2} f\left(\frac{u^n + u^{n+1}}{2}\right) \\ &= R(hA) \left( u^n + \frac{ih}{1+ihA/2} f\left(\frac{u^n + u^{n+1}}{2}\right) \right). \end{aligned} $$

Let

$$ v^{n+1} = R(hA)^* u^{n+1}. $$

We have

$$ v^{n+1} = u^n + \frac{ih}{1+ihA/2} f\left(\frac{u^n + R(hA)v^{n+1}}{2}\right). $$

For $\varepsilon > 0$ and $u$ fixed, we define the map

$$ (3.3) \qquad v \mapsto F_{\varepsilon,h,u}(v) = u + \frac{i\varepsilon}{1+ihA/2} f\left(\frac{u+R(hA)v}{2}\right). $$

With $R(hA)$ being an isometry of $V_{\delta x}$, we have that $\frac{1}{2}(u + R(hA)v) \in B_{\delta x}(\frac{3\|u\|_{\delta x}}{2})$ for any $v \in B_{\delta x}(2\|u\|_{\delta x})$ and thus $\|\phi(\frac{u+R(hA)v}{2})\|_{\delta x} \le C\|u\|_{\delta x}^{2r+1}$ for any $v \in B_{\delta x}(2\|u\|_{\delta x})$ and for some constant $C > 0$ depending only on $r$. Finally, we have

$$ \left\| \frac{w}{1 + ihA/2} \right\|_{\delta x} \le \|w\|_{\delta x} $$

similarly to the arguments of the proof of lemma 3.2 and, thus, we deduce that

$$ \|F_{\varepsilon,h,u}(v)\|_{\delta x} \le \|u\|_{\delta x} + C\varepsilon\|u\|_{\delta x}^{2r+1} \le 2\|u\|_{\delta x} $$

for $\varepsilon \le \varepsilon_0 = C^{-1}\|u\|_{\delta x}^{2r}$ if $v \in B_{\delta x}(2\|u\|_{\delta x})$. This shows that $F_{\varepsilon,h,u}$ maps $B_{\delta x}(2\|u\|_{\delta x})$ to itself (so long as $0 < \varepsilon \le \varepsilon_0$). Now we estimate for $v$ and $w$ in $B_{\delta x}(2\|u\|_{\delta x})$:

$$ \|F_{\varepsilon,h,u}(v) - F_{\varepsilon,h,u}(w)\|_{\delta x} \le C\varepsilon\|u\|_{\delta x}^{2r} \|v - w\|_{\delta x}, $$

where $C$ depends on $r$ but not on $K$, $\delta x$, and thus for $0 < \varepsilon < \varepsilon_0$, $F_{\varepsilon,h,u}$ is a contraction from $B_{\delta x}(2\|u\|_{\delta x})$ to $B_{\delta x}(2\|u\|_{\delta x})$ ensuring the existence of a fixed point $v \in B_{\delta x}(2\|u\|_{\delta x})$ such that $v = F_{\varepsilon,h,u}(v)$ and we can define

$$ (3.4) \qquad \Psi_{\varepsilon}^{h}(u) := v. $$
Note that this argument also shows the well-posedness of the midpoint rule for h small enough, depending on the size of the numerical solution.

Next we would like to show that $\Psi_{\varepsilon}^h(u)$ has an expansion of the form

$$ (3.5) \qquad \Psi_{\varepsilon}^{h}(u) = \sum_{k \ge 0} \varepsilon^{k} \Psi_{h,k}(u). $$

For this we apply the implicit function theorem. We define

$$ g(\varepsilon, v) := v - F_{\varepsilon, h, u}(v), $$

where $v$ is viewed as a vector in $\mathbb{R}^{2K+1} \times \mathbb{R}^{2K+1}$ through the identification $v = p + iq$. The map is entire in $(\varepsilon, p, q)$ (note that it is an entire function of $(\varepsilon, v, \bar{v})$ and moreover $\Psi_{\varepsilon}^h$ is the solution to $g(\varepsilon, \Psi_{\varepsilon}^h(u)) = 0$). Let us consider the Jacobian of $g$ around $v = u =: \Psi_{h,0}(u)$ applied to a vector $z = r + is \in \mathbb{C}^{2K+1} \equiv \mathbb{R}^{4K+2}$:

$$ dg_{\varepsilon,h,u}(u)(z) = z - \frac{i\varepsilon\lambda(r+1)}{2+ihA} \left[ \left| \frac{u+R(hA)v}{2} \right|^{2r} \bullet (R(hA)z) \right] \\ - \frac{i\lambda\varepsilon r}{2+ihA} \left[ \left| \frac{u+R(hA)v}{2} \right|^{2r-2} \bullet \left( \frac{u+R(hA)v}{2} \right)^2 \bullet (R(hA)^*\bar{z}) \right], $$

where $\bullet$ denotes element-wise multiplication of two vectors. Thus we have by the above lemmas and the fact that we already showed $\|\Psi_{\varepsilon}^h(u)\|_{\delta x} \le 2\|u\|_{\delta x}$,

$$ \|dg_{\varepsilon,h,u}(u)(z) - z\|_{\delta x} \le C\varepsilon\|u\|_{\delta x}^{2r}\|z\|_{\delta x}, $$

for some constant $C$ independent of $K, \delta x$ and $u$. Therefore, so long as $0 < \varepsilon < \varepsilon_0 = \frac{1}{2}C^{-1}\|u\|_{\delta x}^{-2r}$, the Jacobian of the map $g(\varepsilon, v)$ is invertible and thus, by analyticity of $g$ and the implicit function theorem, the map $\varepsilon \mapsto \Psi_{\varepsilon}^h(u)$ is analytic in the region $|\varepsilon| < \varepsilon_0$. Therefore, the expression (3.5) indeed holds for a sequence of vector fields $\Psi_{h,k}(u) \in V_{\delta x}$ which, by Cauchy's estimate satisfy

$$ (3.6) \qquad \|\Psi_{h,k}(u)\|_{\delta x} \le \frac{\sup_{|\varepsilon|=\varepsilon_0} \|\Psi_{\varepsilon}^h(u)\|_{\delta x} \varepsilon_0^{-k}}{2\pi} \le C^k \|u\|_{\delta x}^{2k+1}, $$

where we have increased $C$ appropriately without changing the notation in the interest of simplicity. Finally, to show that each $\Psi_{h,k}$ is a homogeneous polynomial of degree $\le 2k+1$ we proceed by induction on $k$ as follows: The statement is trivially true for $k=0$ since $\Psi_{h,0}(u) = u$. Suppose we have shown the claim for $\Psi_{h,j}$ with $0 \le j \le k-1$. We consider the expansion

$$ (3.7) \qquad \Psi_{\varepsilon}^{h}(u) = \sum_{j=0}^{k} \varepsilon^{j} \Psi_{h,j}(u) + \varepsilon^{k+1} \mathcal{R}_{k} $$

then we have

$$ \Psi_{\varepsilon}^{h}(u) = F_{\varepsilon, h, u}(\Psi_{\varepsilon}^{h}(u)), $$
and expanding the right hand side (F is a composition of polynomials and linear operators) we find

$$ \Psi_{\varepsilon}^{h}(u) = u + \frac{i\lambda\varepsilon}{1 + ihA/2} \sum_{j=0}^{k} \varepsilon^{j} \sum_{\substack{m_1+\cdots+m_r+n_1+\cdots+n_{r+1}=j \\ 0 \le m_1, \dots, m_r, n_1, \dots, n_{r+1} \le j}} \overline{\alpha_{m_1}} \cdots \overline{\alpha_{m_r}} \alpha_{n_1} \cdots \alpha_{n_{r+1}} + \varepsilon^{k+2} \tilde{\mathcal{R}}_k, $$

where

$$ \alpha_j = \begin{cases} \frac{u+R(hA)u}{R(hA)^2}, & j=0, \\ \Psi_{h,j}(u), & j \ge 1, \end{cases} $$

and the map $u \mapsto \tilde{\mathcal{R}}_k(u)$ is bounded. Using (3.7) and comparing the coefficient on both sides immediately shows that each $\Phi_{h,\varepsilon,k}(u)$ is a homogeneous polynomial of degree $\le 2k+1$ in $u$. In particular, this provides a recursive way of computing these expressions and the first two of them are

$$ (3.8) \qquad \Psi_{h,0}(u) = u, \quad \Psi_{h,1}(u) = \frac{i\lambda}{1+ihA/2} \left| \frac{u+R(hA)u}{2} \right|^{2r} \frac{u+R(hA)u}{2}. $$

Finally, we note that the above restrictions on $\varepsilon$ were completely independent of $h > 0$, thus we can take $\varepsilon = h < \varepsilon_0$ and the result follows with $\Psi^h = \Psi_h^h$. $\square$

Note that using (3.8), we can write the expansion in the form

$$ \begin{aligned} \Psi^h(u) &= u + \frac{i\lambda h}{1 + ihA/2} \left| \frac{i}{1 - ihA/2} u \right|^{2r} \frac{i}{1 - ihA/2} u + \mathcal{O}(\|u\|_{\delta x}^{2r+1}) \\ &= u + ih(\delta x)^{-1} \nabla_{\bar{u}} P_{0,h}(u) + \mathcal{O}(\|u\|_{\delta x}^{2r+1}), \end{aligned} $$

where

$$ (3.9) \qquad P_{0,h}(u) = \frac{\lambda}{r+1} \delta x \sum_{\ell=-K}^{K} \left| \left( \frac{1}{1-ihA/2} u \right)_{\ell} \right|^{2r+2}, $$

is the $(2r+2)$-part of the original energy (see (2.3)) composed with the pseudo-differential operator $(1 - ihA/2)^{-1}$. This suggests that $\psi^h$ is the Hamiltonian flow at time $h$ of a modified bounded pseudo-differential operator depending on $h$. For this we recall the general Hamiltonian formulation in the $u$ coordinate (2.11) and exploit Lemma 2.12.

**Proposition 3.3.** The map $u \mapsto \Psi_{\varepsilon}^{h}(u)$ as introduced in (3.4) is symplectic for any permissible choice of $\varepsilon, h$ (in particular for $\varepsilon = h$, i.e. $\Psi^h$ is symplectic). Moreover, there exists a formal real-valued Hamiltonian

$$ P_h = P_{0,h} + \sum_{k \ge 1} h^k P_{k,h}, $$

where the $P_{k,h}$ are real-valued homogeneous polynomials with

$$ |P_{k,h}(u)| \le C^{(k)} \|u\|_{\delta x}^{2r(k+1)+2} $$
with $C^{(k)} > 0$ independent of $K, \delta x$, such that the following holds. For any given $N$ there exists $C_N, \tau_0^{(N)} > 0$ such that if we define

$$P_h^{(N)} = P_{0,h} + \sum_{k=1}^{N} h^k P_{k,h}$$

then we have for any $u \in B_{\delta x}(R), h \in [0, \tau_0^{(N)})$

$$\Psi^h(u) = \Phi_{P_h^{(N)}}^h(u) + R_{h,N}(u),$$

where $\Phi_{P_h^{(N)}}^t$ is the flow associated with the vector field $i\delta x^{-1}\nabla_{\bar{u}}P_h^{(N)}$ and the remainder terms $R_{h,N}$ take the form

$$(3.10) \qquad R_{h,N}(u) = \sum_{\ell=N+2}^{\infty} h^{\ell} \mathcal{R}_{h,\ell}^{(N)}(u),$$

and each $\mathcal{R}_{h,\ell}^{(N)}$ is a polynomial vector field satisfying the bound

$$(3.11) \qquad \|\mathcal{R}_{h,\ell}^{(N)}(u)\|_{\delta x} \le \tilde{C}^{(N)} (C^{(N)})^{\ell} \|u\|_{\delta x}^{2r\ell+1},$$

for some $\tilde{C}^{(N)}, C^{(N)} > 0$ independent of $K, \delta x$.

**Lemma 3.4** (Cauchy-Kovalevskaya). Suppose $\mathcal{P}: \mathbb{C}^{2K+1} \to \mathbb{R}$ is a real-valued homogeneous polynomial such that

$$|\mathcal{P}(u)| \le C \|u\|_{\delta x}^d, \quad \forall u \in \mathbb{C}^{2K+1}$$

for some constant $C > 0$, and some $d \in \mathbb{N}_{\ge 2}$. Let us fix $R > 0$. Then, there is a $\tau_0 > 0$ dependent only on $d, C$ such that for all $0 \le t < \tau_0$ the flow map $\Phi_{\mathcal{P}}^t: B_{\delta x}(R) \subset \mathbb{C}^{2K+1} \to \mathbb{C}^{2K+1}$ associated with the differential equation

$$(3.12) \qquad \frac{du}{dt} = i(\delta x)^{-1} \nabla_{\bar{u}} \mathcal{P}(u),$$

exists, is symplectic in the sense of Definition 2.4, and is locally analytic. In particular there is $\tilde{C}_0 > 0$ (again dependent only on $C, d$) such that for any $u \in B_{\delta x}(R)$ and $\varepsilon \le \tau_0$ we have

$$(3.13) \qquad \Phi_{\mathcal{P}}^{\varepsilon}(u) = u + i\varepsilon(\delta x)^{-1} \nabla_{\bar{u}} \mathcal{P}(u) + \sum_{k \ge 2} \varepsilon^k \mathcal{R}_k(u),$$

where, for each $k \ge 2$, $\mathcal{R}_k(u)$ is a polynomial vector field satisfying the bound

$$(3.14) \qquad \|\mathcal{R}_k(u)\|_{\delta x} \le \tilde{C}_0 d^k C^k \|u\|_{\delta x}^{k(d-2)+1}.$$

**Proof of Lemma 3.4.** The existence of a locally analytic solution is a direct consequence of the standard Cauchy-Kovalevskaya theorem. The form of $\mathcal{R}_k$ and the bound (3.14) can be obtained by plugging (3.13) into (3.12) (we can differentiate the series by absolute convergence of the sum) and comparing terms order by order analogously to the final part of the proof of Proposition 3.1. $\square$
Proof of Proposition 3.3. This proof is based on the proof of Theorem IX.3.1 in [17], but adapted to ensure each estimate is dimension-independent. To begin with, let us denote by $\Phi_{h,\varepsilon}$ the midpoint rule for the Hamiltonian system (with rescaled nonlinearity)

$$ \frac{du}{dt} = iAu + i\frac{\varepsilon}{h}f(u). $$

Clearly, this implies that $\Phi_{h,\varepsilon}$ is a symplectic map for any choice of $h, \varepsilon$ for which $\Phi_{h,\varepsilon}$ is well-defined. By construction we have with the notation of the previous proof, see (3.3)

$$ (3.15) \qquad \Phi_{h,\varepsilon} = R(hA) \circ \Psi_{\varepsilon}^{h} \iff \Psi_{\varepsilon}^{h} = R(hA)^{*} \circ \Phi_{h,\varepsilon}, $$

i.e. $\Psi_{\varepsilon}^{h}$ is the composition of two symplectic maps $(R(hA)^{*})$ is the adjoint of the midpoint rule applied to the linear Schrödinger equation), hence symplectic.

Moreover, the function $\Psi_{h,\varepsilon}$ admits, for a fixed $h > 0$, a convergent expansion in powers of $\varepsilon$, see (3.5) and the bound (3.6) shows that the radius of convergence is independent of $h$. As symplectic $\varepsilon$-perturbation of the identity, we can thus construct a modified vector field at any order in $\varepsilon$ by following the classical method of [17, Chapter IX] and shows that this modified vector is Hamiltonian. We then conclude by taking $\varepsilon = h$ small enough. Let us recall this construction, and show how the estimates obtained are independent of the parameters $K$ and $\delta x$.

We proceed by induction on $N$. To begin with, we note for $P_{0,h}$ given by (3.9), we have $\Psi_{1,h} = i(\delta x)^{-1}\nabla_{\pi}P_{0,h}$ where $\Psi_{h,1}$ is as defined in Proposition 3.1 and is given explicitly by (3.8). Thus, by Lemma 3.4, for any given $R > 0$ there is a $\tau_{0}^{(0)} > 0$ such that the corresponding Hamiltonian flow takes the form

$$ \Phi_{P_{0,h}}^{\varepsilon}(u) = u + \varepsilon \Psi_{h,1}(u) + \sum_{\ell=2}^{\infty} \varepsilon^{\ell} \mathcal{R}_{h,\ell}^{(0)}(u) $$

for all $0 < \varepsilon, h \le \tau_{0}^{(0)}$, $u \in B_{\delta x}(R)$, where each $\mathcal{R}_{h,\ell}^{(0)}$ is a homogeneous polynomial vector field satisfying the bound

$$ \| \mathcal{R}_{h,\ell}^{(0)}(u) \|_{\delta x} \le \tilde{C}^{(0)} (C^{(0)})^{\ell} \|u\|_{\delta x}^{2r\ell+1}, $$

for some $\tilde{C}^{(0)}, C^{(0)} > 0$ independent of $K, \delta x$. Taking $h = \varepsilon < \tau_{0}^{(0)}$ and using Proposition 3.1, this shows Proposition 3.3 for the case $N = 0$. For the induction step we now assume we have shown for a given $N$ that there are real-valued polynomials $P_{0,h}, \dots, P_{N,h}$ such that for any $R > 0$ there is $\tau_{0}^{(N)} > 0$ such that for all $u \in B_{\delta x}(R)$ and $\varepsilon, h < \tau_{0}^{(N)}$

$$ (3.16) \qquad \Phi_{P_{h,\varepsilon}}^{\varepsilon}(u) = \Psi_{\varepsilon}^{h}(u) + \sum_{\ell=N+2}^{\infty} \varepsilon^{\ell} \mathcal{R}_{h,\ell}^{(N)}(u), $$

where $\Psi_{\varepsilon}^{h}$ is as defined in (3.5) and each $\mathcal{R}_{h,\ell}^{(N)}$ is a homogeneous polynomial vector field satisfying the bound

$$ (3.17) \qquad \| \mathcal{R}_{h,\ell}^{(N)}(u) \|_{\delta x} \le \tilde{C}^{(N)} (C^{(N)})^{\ell} \|u\|_{\delta x}^{2r\ell+1}, $$
for some $\tilde{C}^{(N)}$, $C^{(N)} > 0$ independent of $K, \delta x$, and where

$$ (3.18) \qquad P_{h,\varepsilon}^{(N)} = P_{0,h} + \sum_{k=1}^{N} \varepsilon^k P_{k,h}. $$

We now seek to construct $P_{h,N+1}$ such that the analogous statement is true for $P_{h,\varepsilon}^{(N+1)}$. We note that using the estimate (3.11) the series on the right hand side of (3.10) is locally uniformly convergent, thus we can exchange order of differentiation and summation to find for $u \in B_{\delta x}(R)$ and $\varepsilon < \tau_0^{(N)}$ (using the notation introduced in (2.19)):

$$ \nabla_{\bar{u}}^{T} \Phi_{P_{h,\varepsilon}^{(N)}}^{\varepsilon}(u) = \nabla_{\bar{u}}^{T} \Psi_{\varepsilon}^{h}(u) + \sum_{\ell=N+2}^{\infty} \varepsilon^{\ell} \nabla_{\bar{u}}^{T} \mathcal{R}_{h,\ell}^{(N)}(u). $$

According to Lemma 3.4 and (3.15) both $\Phi_{P_{\varepsilon}^{(N)}}^{\varepsilon}$ and $\Psi_{\varepsilon}^{h}$ are symplectic in the sense of Definition 2.4. Thus in slight abuse of notation we have (cf. (2.12)):

$$ iI = \left( \nabla_{\bar{u}}^{T} \Phi_{P_{h,\varepsilon}^{(N)}}^{\varepsilon}(u) \right)^{T} i \nabla_{\bar{u}}^{T} \Phi_{P_{h,\varepsilon}^{(N)}}^{\varepsilon}(u) = \left( \nabla_{\bar{u}}^{T} \Psi_{\varepsilon}^{h}(u) \right)^{T} i \nabla_{\bar{u}}^{T} \Psi_{\varepsilon}^{h}(u), $$

where $I$ denotes the $(2K+1) \times (2K+1)$ identity matrix and $i$ represents multiplication by $-J$. Therefore, using the expression (3.10), we find

$$ (3.19) \qquad \varepsilon^{N+2} \left( \nabla_{\bar{u}}^{T} \Psi_{\varepsilon}^{h}(u) \right)^{T} i \nabla_{\bar{u}}^{T} \mathcal{R}_{h,N+2}^{(N)}(u) + \varepsilon^{N+2} \left( \nabla_{\bar{u}}^{T} \mathcal{R}_{h,N+2}^{(N)}(u) \right)^{T} i \nabla_{\bar{u}}^{T} \Psi_{\varepsilon}^{h}(u) + \varepsilon^{N+3} \mathcal{R}_{h,\varepsilon,N}(u) = 0, $$

for some remainder $\mathcal{R}_{h,\varepsilon,N}$ uniformly bounded on $u \in B_{\delta x}(R)$, $h, \varepsilon \in [0, \tau_0^{(N)})$. We now note that $\Psi_{\varepsilon}^{h}(u) = u + \mathcal{O}(\varepsilon)$

$$ \nabla_{\bar{u}}^{T} \Psi_{\varepsilon}^{h}(u) = I + \mathcal{O}(\varepsilon), $$

where the $\mathcal{O}(\varepsilon)$ represents terms which are uniformly bounded above by $C\varepsilon$ on $u \in B_{\delta x}(R)$ and with $h, \varepsilon \in [0, \tau_0^{(N)})$. Thus we have, by dividing (3.19) by $\varepsilon^{N+2}$ and taking $\varepsilon \to 0$, that

$$ i \nabla_{\bar{u}}^{T} \mathcal{R}_{h,N+2}^{(N)}(u) = \left( i \nabla_{\bar{u}}^{T} \mathcal{R}_{h,N+2}^{(N)}(u) \right)^{T}, $$

i.e. that the $(2K+1) \times (2K+1)$ matrix represented by $i \nabla_{\bar{u}}^{T} \mathcal{R}_{h,N+2}^{(N)}(u)$ is symmetric. Thus $\mathcal{R}_{h,N+2}^{(N)}(u)$ satisfies the assumptions of Lemma 2.12 and there is a real-valued homogeneous polynomial Hamiltonian $P_{N+1,h}$ such that

$$ (3.20) \qquad -\mathcal{R}_{h,N+2}^{(N)}(u) = i(\delta x)^{-1} \nabla_{\bar{u}} P_{N+1,h}, $$

and

$$ |P_{N+1,h}(u)| \leq C \|u\|_{\delta x}^{2r(N+2)+2}, \quad \forall u \in B_R. $$

Let us now consider

$$ P_{\varepsilon}^{(N+1)} := P_{0,h} + \sum_{k=1}^{N+1} \varepsilon^k P_{k,h}. $$
Then, it can be shown analogously to Lemma 3.4, that there is $\tau_0^{(N+1)} \le \tau_0^{(N)}$ such that for all $0 \le t < \tau_0^{(N+1)}$, $u \in B_{\delta x}(R)$ the flows $\Phi_{P_{\epsilon}^{(N+1)}}^t(u)$, $\Phi_{P_{\epsilon}^{(N)}}^t(u)$ exist and in particular that $\Phi_{P_{\epsilon}^{(N+1)}}^t(u)$ takes the following form

$$ (3.21) \qquad \Phi_{P_{\epsilon}^{(N+1)}}^t(u) = u + \sum_{\ell,p=1}^{\infty} t^{\ell} \epsilon^p Q_{h,\ell,p}^{(N+1)}(u), $$

where $Q_{h,\ell,p}^{(N+1)}$ are homogeneous polynomial vector fields satisfying

$$ \|Q_{h,\ell,p}^{(N+1)}(u)\|_{\delta x} \le C_N^{\ell+p} \|u\|_{\delta x}^{2r(\ell+p)+1} $$

for some $C_N > 0$ depending only on $N$. Note that by taking $t = \epsilon$, this implies the existence of the expansion (3.10) at order $N+1$, and we only have to show bounds on the remainder term and the cancellation of the term of order $N+2$.

For $\epsilon \le \tau_0^{(N+1)}$, we introduce two maps $g, f : B_{\delta x}(R) \to \mathbb{C}^{2K+1}$ by

$$ g_{\epsilon}(u) := i(\delta x)^{-1} \nabla_{\bar{u}} P_{\epsilon}^{(N)}, $$

$$ f_{\epsilon}(u) := i(\delta x)^{-1} \nabla_{\bar{u}} P_{\epsilon}^{(N+1)}. $$

For this choice we have

$$ (3.22) \qquad \frac{d\Phi_{P_{\epsilon}^{(N)}}^t}{dt} = g_{\epsilon}(\varphi_{P_{\epsilon}^{(N)}}^t), \quad 0 < t, h, \epsilon < \tau_0^{(N+1)}, $$

$$ (3.23) \qquad \frac{d\Phi_{P_{\epsilon}^{(N+1)}}^t}{dt} = f_{\epsilon}(\varphi_{P_{\epsilon}^{(N+1)}}^t), \quad 0 < t, h, \epsilon < \tau_0^{(N+1)}. $$

We proceed in two steps. Firstly, we shall show that there is $C > 0$ such that for all $u \in B_{\delta x}(R)$, $0 \le t, \epsilon, h < \tau_0^{(N+1)}$:

$$ (3.24) \qquad \| \Phi_{P_{\epsilon}^{(N)}}^t(u) - \Phi_{P_{\epsilon}^{(N+1)}}^t(u) \|_{\delta x} \le t \epsilon^{N+1} C \|u\|_{\delta x}^{2r(N+2)+1}. $$

We note that

$$ \| f_{\epsilon}(\Phi_{P_{\epsilon}^{(N+1)}}^t(u)) - g_{\epsilon}(\Phi_{P_{\epsilon}^{(N+1)}}^t(u)) \|_{\delta x} = \epsilon^{N+1} \| \mathcal{R}_{h,N+2}^{(N)}(\Phi_{P_{\epsilon}^{(N+1)}}^t(u)) \|_{\delta x} \le \epsilon^{N+1} C \|u\|_{\delta x}^{2r(N+2)+1}, $$

for some $C > 0$ so long as $u \in B_{\delta x}(R)$. This holds true because $\mathcal{R}_{h,N+2}^{(N)}$ is a polynomial vector field in its argument with bounded coefficients and the local well-posedness of the Hamiltonian system corresponding to $P_{\epsilon}^{(N+1)}$ established in (3.21). Moreover, $g$ is a polynomial with bounded coefficients, and it follows immediately that it is Lipschitz continuous (see Lemma 2.14) with a constant independent of $K, \delta x$, i.e. there is a $C > 0$ such that for $u, v \in B_{\delta x}(R)$, $\epsilon \in [0, \tau_0^{(N+1)})$ we have

$$ \|g_{\epsilon}(u) - g_{\epsilon}(v)\|_{\delta x} \le C(\|u\|_{\delta x} + \|v\|_{\delta x})^{2r} \|u - v\|_{\delta x} $$

uniformly in $h, \epsilon$ small enough. Thus we have (noting that the flows are equal at $t = 0$)

$$ \| \Phi_{P_{\epsilon}^{(N+1)}}^t(u) - \Phi_{P_{\epsilon}^{(N)}}^t(u) \|_{\delta x} \le C_2 \int_0^t \epsilon^{N+1} \|u\|_{\delta x}^{2r(N+2)+1} \, ds \le C_2 t \epsilon^{N+1} \|u\|_{\delta x}^{2r(N+2)+1}, $$

for some $C_2 > 0$ which depends on $R$ but not on $K, \delta x$, thus completing the estimate (3.24).
In the second step we will show that

$$ (3.25) \qquad \| \Phi_{P_t^{(N+1)}}^t(u) - \Phi_{P_t^{(N)}}^t(u) + t\epsilon^{N+1} \mathcal{R}_{h,N+2}^{(N)}(u) \|_{\delta x} \le Ct^{N+3} \|u\|_{\delta x}^{2(N+3)+1}, $$

for some constant $C > 0$ independent of $K, \delta x$. For this we note that by (3.22) & (3.23) we have

$$ (3.26) \qquad \begin{aligned} \Phi_{P_{\epsilon}^{(N+1)}}^{t}(u) - \Phi_{P_{\epsilon}^{(N)}}^{t}(u) &= \int_{0}^{t} f_{\epsilon}(\Phi_{P_{\epsilon}^{(N+1)}}^{s}(u)) - g_{\epsilon}(\Phi_{P_{\epsilon}^{(N)}}^{s}(u))ds \\ &= \int_{0}^{t} f_{\epsilon}(\Phi_{P_{\epsilon}^{(N+1)}}^{s}(u)) - g_{\epsilon}(\Phi_{P_{\epsilon}^{(N+1)}}^{s}(u))ds \\ &\quad + \int_{0}^{t} g_{\epsilon}(\Phi_{P_{\epsilon}^{(N+1)}}^{s}(u)) - g_{\epsilon}(\Phi_{P_{\epsilon}^{(N)}}^{s}(u))ds. \end{aligned} $$

Now we note that

$$ (3.27) \qquad \begin{aligned} \int_0^t f_\epsilon(\Phi_{P_\epsilon^{(N+1)}}^s(u)) - g_\epsilon(\Phi_{P_\epsilon^{(N+1)}}^s(u))ds &= -\int_0^t \epsilon^{N+1} \mathcal{R}_{h,N+2}^{(N)}(\Phi_{P_\epsilon^{(N+1)}}^s(u))ds \\ &= -\mathcal{R}_{h,N+2}^{(N)}(u) \int_0^t \epsilon^{N+1}ds + t^2 \epsilon^{N+2} Q_{h,\epsilon,t}^{(1)}(u), \end{aligned} $$

where $Q_{h,\epsilon,t}^{(1)}$ is a function of $u$ bounded uniformly on $B_{\delta x}(R)$, $(h, \epsilon, t) \in [0, \tau_0^{(N+1)}]$ and is bounded by $C_N \|u\|_{\delta x}^{2(N+3)+1}$ with a convergent expansion of the form (3.21). In this final line we used the observation from (3.21) that

$$ \varphi_{P_{\epsilon}^{(N+1)}}^{s}(u) = u + tQ_{h,\epsilon,t}^{(2)}(u), $$

where $Q_{h,\epsilon,t}^{(2)}(t, u)$ is some analytic function of $u$ bounded uniformly on $u \in B_{\delta x}(R)$, $h, \epsilon, t \in [0, \tau_0^{(N+1)}]$. Thus, combining (3.26) & (3.27) we obtain that given $R > 0$ there is a constant $C_2$ depending only on $N$ such that

$$ \begin{aligned} \left\| \Phi_{P_{\epsilon}^{(N+1)}}^{t}(u) - \Phi_{P_{\epsilon}^{(N)}}^{t} + t\epsilon^{N+1} \mathcal{R}_{h,N+2}^{(N)}(u) \right\|_{\delta x} &\le t^2 \epsilon^{N+1} \|u\|_{\delta x}^{2(N+3)+1} \\ &\quad + Ct \|u\|_{\delta x}^{2r} \left\| \Phi_{P_{\epsilon}^{(N+1)}}^{t}(u) - \Phi_{P_{\epsilon}^{(N)}}^{t} \right\|_{\delta x} \\ &\le t^2 \epsilon^{N+1} C_2 \|u\|_{\delta x}^{2(N+3)+1} \end{aligned} $$

for any $u \in B_{\delta x}(R)$. This completes the proof of (3.25), hence the induction step and therefore the proof of the desired result (by taking again $h = t = \epsilon < \tau_0^{(N)}$ for any given choice of $N$). $\square$

### 4. MODIFIED ENERGY FOR THE MIDPOINT RULE

Proposition 3.3 essentially shows that the midpoint rule can (to arbitrary desired order) be written as a Hamiltonian splitting method. In this section we will establish our central result, Theorem 4.3 which shows the existence of a modified energy for the midpoint rule, exploiting this “approximate” splitting formulation. We begin by introducing the vector field associated with the linear part of the flow of (2.5).
**Lemma 4.1.** The vector field $X_{A_0}$ associated with the flow $\Phi_{A_0}^1 = R(hA)$ is given by

$$ (4.1) \qquad X_{A_0} = 2iU^{-1} \arctan \left( \frac{hD}{2} \right) U, $$

where $U, D$ are as in (2.6). Furthermore it satisfies the estimate

$$ (4.2) \qquad \|X_{A_0} u\|_{\delta x} \le 2 \arctan \left( \frac{h}{\delta x^2} \right) \|u\|_{\delta x}. $$

*Proof of Lemma 4.1.* We note by linearity we can integrate (4.1) over $t \in [0, 1]$ to directly recover $R(hA)$. By (2.7) the norm $\| \cdot \|_{\delta x}$ and $A$ can be jointly diagonalised meaning that

$$ \|X_{A_0} u\|_{\delta x} \le \max_{1 \le j \le 2K+1} \left| 2 \arctan \left( \frac{h\lambda_j}{2} \right) \right| \|u\|_{\delta x}, $$

where $\lambda_j, j = 1, \dots, 2K + 1$ are the eigenvalues of $A$. Since $A$ is a symmetric tridiagonal Toeplitz matrix, the eigenvalues are well-known to be

$$ (4.3) \qquad \lambda_j = -\frac{1}{\delta x^2} \left( 2 - 2 \cos \left( \frac{j\pi}{2K+2} \right) \right), $$

whence the estimate (4.2) immediately follows. $\square$

**4.1. Formal construction of the modified energy.** Using the above basics we can now formally construct the modified energy for the midpoint rule. We shall rigorously show that the flow corresponding to this modified energy corresponds (up to arbitrary desired order) to the one of the midpoint rule in Section 4.2. We take a similar approach to [13, 4], but note that in our case we have to include a second small parameter $\varepsilon$ which captures the expansion of the modified vector field corresponding to the nonlinear part of the midpoint rule (3.18), to avoid the confusion with the $h$ appearing in the terms depending on $hA$. Thus, we look for a real Hamiltonian function $Z(t, \varepsilon; u)$ such that

$$ (4.4) \qquad \Phi_{Z(t, \varepsilon)}^1 = \Phi_{A_0}^1 \circ \Phi_{P_{h, \varepsilon}^{(N)}}^t, \quad \forall t, \varepsilon < \tilde{\tau}_0^{(N)}, $$

for some threshold $\tilde{\tau}_0^{(N)} > 0$ to be determined, but which should be independent of $K, \delta x$. In the above, $A_0$ is as defined in Lemma 4.1 and $P_{h, \varepsilon}^{(N)}$ is as defined in (3.18). According to [13, Section 3] we have, formally,

$$ (4.5) \qquad \frac{\partial}{\partial t} \Phi_{Z(t, \varepsilon)}^1 = X_{Q(t, \varepsilon)} \circ \Phi_{Z(t, \varepsilon)}^1, $$

where the modified vector field $X_Q$ has the formal series [13, (3.3)]

$$ (4.6) \qquad X_{Q(t, \varepsilon)} = \sum_{k \ge 0} \frac{1}{(k+1)!} \mathrm{ad}_{X_{Z(t, \varepsilon)}}^k \partial_t X_{Z(t, \varepsilon)}. $$

To establish an expression for $Z$ we would also like to differentiate the right hand side of (4.4). We have

$$ (4.7) \qquad \frac{d}{dt} \left( \Phi_{A_0}^1 \circ \Phi_{P_h^{(N)}}^t (u) \right) = d\Phi_{A_0}^1|_{p=\Phi_{P_h^{(N)}}^t(u)} \circ X_{P_h^{(N)}} \circ \Phi_{P_h^{(N)}}^t(u), $$
where, by linearity, we actually have $d\Phi_{A_0}^1|_{p=\Phi_{P_h^{(N)}}^t}(u) = R(hA)$. Moreover, we recall that
$X_P = i\delta x^{-1}\nabla_{\bar{u}}P$, thus we observe that for $\tilde{P} = P \circ R(hA)^*$, i.e.

$$
\tilde{P}(u, \bar{u}) = P(R(hA)^* u, R(hA)^T \bar{u}),
$$

we find in coordinates

$$
\begin{equation} \tag{4.8}
\begin{split}
(i \delta x^{-1} \nabla_{\bar{u}} \tilde{P})_j &= i \delta x^{-1} \frac{\partial \tilde{P}}{\partial \bar{u}_j} = i \delta x^{-1} \sum_{k, \ell} \frac{\partial P}{\partial \bar{u}_k} \frac{\partial R(hA)^T_{kl}}{\partial \bar{u}_j} \bar{u}_l \\
&= \sum_k R(hA)_{jk} i \delta x^{-1} \frac{\partial P}{\partial \bar{u}_k} = (R(hA) \circ X_P)_j.
\end{split}
\end{equation}
$$

And thus, combining this with (4.7) gives

$$
\begin{equation} \tag{4.9}
\begin{split}
\left. \frac{d}{dt} \left( \Phi_{A_0}^1 \circ \Phi_{P_{h,\epsilon}^{(N)}}^t (u) \right) \right|_{t=t} &= X_{P_{h,\epsilon}^{(N)} \circ R(hA)^*} \circ \Phi_{A_0}^1 \circ \Phi_{P_{h,\epsilon}^{(N)}}^t (u) \\
&= X_{P_{h,\epsilon}^{(N)} \circ R(hA)^*} \circ \Phi_{Z(t,\epsilon)}^1.
\end{split}
\end{equation}
$$

Thus combining (4.6) & (4.9) yields the following equation to be satisfied by $X_{Z(t,\epsilon)}$:

$$
(4.10) \qquad \sum_{k \ge 0} \frac{1}{(k+1)!} \mathrm{ad}_{X_{Z(t,\epsilon)}}^k \partial_t X_{Z(t,\epsilon)} = X_{P_{h,\epsilon}^{(N)} \circ R(hA)^*},
$$

where $\mathrm{ad}_{X_H}(X_G) := [X_H, X_G]$. We shall now seek a formal solution to this equation. We start by noting that (4.10) is formally equivalent to

$$
\partial_t X_{Z(t, \varepsilon)} = \sum_{k \ge 0} \frac{B_k}{k!} \mathrm{ad}_{X_{Z(t, \varepsilon)}}^k X_{P_{h, \varepsilon}^{(N)} \circ R(hA)^*},
$$

where $B_k$ are the Bernoulli numbers. Note this inversion is only valid for the full series when $\|\mathrm{ad}_{X_{Z(t,\varepsilon)}}^k \partial_t X_{Z(t,\varepsilon)}\|_{\delta x} < (2\pi)^k$, but the term-by-term conditions are algebraically equivalent regardless of the size of $\|\mathrm{ad}_{X_{Z(t,\varepsilon)}}^k \partial_t X_{Z(t,\varepsilon)}\|_{\delta x}$. In order to find (a suitable truncation) of $X_Z$ we try the following Ansatz

$$
(4.11) \qquad Z_{t,\epsilon} = \sum_{\ell,j=0}^{\infty} t^{\ell} \epsilon^j Z_{\ell,j,h},
$$

with

$$
Z_{0,j,h} := \begin{cases} A_0, & \text{if } j=0, \\ 0, & \text{otherwise.} \end{cases}
$$

This translates to an equivalent Ansatz for $X_{Z(t,\epsilon)}$ by the results in Section 2.6. For notational simplicity we will in the following write $Z_{\ell,j}$ for $Z_{\ell,j,h}$. Matching terms formally order-by-order in both $t$ and $\epsilon$ this leads to the recursion

$$
(4.12) \quad (\ell + 1) X_{Z_{\ell+1,j}} = \sum_{k \ge 0} \frac{B_k}{k!} \sum_{m=0}^{\min\{j,N\}} \sum_{\substack{\ell_1+\ell_2+\cdots+\ell_k=\ell \\ j_1+j_2+\cdots+j_k=j-m}} \mathrm{ad}_{X_{Z_{\ell_1,j_1}}} \cdots \mathrm{ad}_{X_{Z_{\ell_k,j_k}}} X_{P_{m,h} \circ R(hA)^*},
$$

where $P_{m,h}$ are as constructed in Proposition 3.3.
**Lemma 4.2.** Fix $N \in \mathbb{N}$ and suppose we have for some $\tilde{\epsilon} > 0$ that the following CFL condition is satisfied

$$ (4.13) \qquad h \le \delta x^2 \tan \left( \frac{\pi - \tilde{\epsilon}}{2(2r(N+1)+1)} \right). $$

Then, for every $0 \le l, j$ with $l + j \le N$, the Hamiltonian $Z_{l,j}$ and associated vector field $X_{Z_{l,j}}$ is uniquely defined by (4.12), with the infinite series on the right hand side of (4.12) existing and converging. Moreover, $X_{Z_{l,j}} \in \mathcal{P}_{2r(l+j)+1}$, and $X_{Z_{l,j}} = 0$ whenever $j > Nl$ and there is a constant $C_{l,j} > 0$ depending only on $N, l, j$ such that for all $K, \delta x > 0$,

$$ (4.14) \qquad \|X_{Z_{l,j}}\|_{\delta x} \le C_{l,j}. $$

*Proof of Lemma 4.2.* We shall prove the result by induction on $l$. Clearly $X_{Z_{0,j}} \in \mathcal{P}_{2rj+2}$ for all $j \ge 0$ and $X_{Z_{0,j}} = 0$ when $j > 1$. For the induction step we proceed as follows. Firstly, we note that the only non-zero contributions to $Z_{l+1,j}$ in (4.12) arise if $j_i \le Nl_i, i = 1, \dots, k$, i.e. if

$$ j = j_1 + \dots + j_k + m \le l_1 N + \dots + l_k N + m \le l N + N, $$

i.e. if $j \ge N(l+1)$, so $X_{Z_{l+1,j}} = 0$, whenever $j > N(l+1)$. Next we note that any non-zero contribution to $X_{Z_{l,j}}$ is a polynomial of degree $\le 2r(l+j+1)+1$, since

$$ \deg \left( \mathrm{ad}_{X_{Z_{l_1,j_1}}} \cdots \mathrm{ad}_{X_{Z_{l_k,j_k}}} X_{P_{m,h}} \right) = 2r(m+1) + 1 - 2k + \sum_{q=1}^{k} \deg Z_{l_q,j_q} $$

and so, if $l_1 + l_2 + \dots + l_k = l$, $j_1 + j_2 + \dots + j_k = j - m$, then we have

$$ \begin{aligned} \deg \left( \mathrm{ad}_{X_{Z_{l_1,j_1}}} \cdots \mathrm{ad}_{X_{Z_{l_k,j_k}}} X_{P_{m,h}} \right) &\le 2r(m+1) + 1 - 2k + 2r(l+j-m) + 2k \\ &= 2r(l+j+1) + 1. \end{aligned} $$

Thus, so long as the term on the right hand side of (4.12) converges in $\| \cdot \|_{\delta x}$, $X_{Z_{l+1,j}} \in \mathcal{P}_{2r(l+1+j)+1}$. To show this, we can use the estimate (2.15):

$$ \begin{aligned} & (l+1) \|X_{Z_{l+1,j}}\|_{\delta x} \\ &\le \sum_{k \ge 0} \frac{B_k}{k!} \sum_{m=0}^{\min\{j,N\}} \sum_{\substack{l_1+l_2+\dots+l_k=l \\ j_1+j_2+\dots+j_k=j-m}} \left( (2r(l_1+j+1) + 2r(\sum_{i=2}^k l_i + j_i + m + 1) + 3) \right) \cdots \\ &\qquad \cdots \left( (2r(l_{k-1}+j_{k-1}) + 2r(l_k + j_k + m + 1) + 3) \right) (2r(l_k + j_k) + 2r(m+1) + 3) \\ &\qquad \|X_{Z_{l_1,j_1}}\|_{\delta x} \cdots \|X_{Z_{l_k,j_k}}\|_{\delta x} \|X_{P_{m,h} \circ R(hA)^*}\|_{\delta x}. \end{aligned} $$

Let us denote by $d_{\max}$ the largest degree appearing in the above estimate, i.e.

$$ d_{\max} = 2r(l + j - m + m + 1) + 1 = 2r(l + j + 1) + 1. $$
Then we can simplify the estimate to

$$
\begin{align*}
& (\ell + 1) \|X_{Z_{\ell+1,j}}\|_{\delta x} \\
\le & \sum_{k \ge 0} \frac{B_k}{k!} \sum_{m=0}^{\min\{j,N\}} (2d_{\max})^k \sum_{\substack{\ell_1+\ell_2+\dots+\ell_k=l \\ j_1+j_2+\dots+j_k=j-m}} \|X_{Z_{\ell_1,j_1}}\|_{\delta x} \dots \|X_{Z_{\ell_k,j_k}}\|_{\delta x} \|X_{P_{m,h} \circ R(hA)^*}\|_{\delta x} \\
\le & \sum_{k \ge 0} \frac{B_k}{k!} \sum_{m=0}^{\min\{j,N\}} (2d_{\max})^k \sum_{\tilde{r}=1}^{\ell} \binom{k}{k-\tilde{r}} \|X_{A_0}\|_{\delta x}^{k-\tilde{r}} \\
& \times \sum_{\substack{\ell_i>0, \ell_1+\ell_2+\dots+\ell_{\tilde{r}}=\ell \\ j_1+j_2+\dots+j_{\tilde{r}}=j-m}} \|X_{Z_{\ell_1,j_1}}\|_{\delta x} \dots \|X_{Z_{\ell_{\tilde{r}},j_{\tilde{r}}}}\|_{\delta x} \|X_{P_{m,h} \circ R(hA)^*}\|_{\delta x}
\end{align*}
$$

and thus

$$
(\ell + 1) \|X_{Z_{\ell+1,j}}\|_{\delta x} \le \sum_{m=0}^{\min\{j,N\}} \sum_{\tilde{r}=1}^{\ell} \left( \sum_{k \ge 0} \frac{B_k}{k!} (2d_{\max})^k \binom{k}{k-\tilde{r}} \|X_{A_0}\|_{\delta x}^{k-\tilde{r}} \right) \\
\times \sum_{\substack{\ell_i > 0, \ell_1 + \ell_2 + \dots + \ell_{\tilde{r}} = \ell \\ j_1 + j_2 + \dots + j_{\tilde{r}} = j - m}} \|X_{Z_{\ell_1,j_1}}\|_{\delta x} \dots \|X_{Z_{\ell_{\tilde{r}},j_{\tilde{r}}}}\|_{\delta x} \|X_{P_{m,h} \circ R(hA)^*}\|_{\delta x}.
$$

Now we know that

$$
f : x \mapsto \sum_{k \ge 0} \frac{B_k}{k!} (2d_{\max})^k x^k
$$

has radius of convergence $|x| < 2\pi/(2d_{\max})$ and that

$$
\frac{d^{\tilde{r}}}{dx^{\tilde{r}}} f(x) = \sum_{k \ge 0} \frac{B_k}{k!} (2d_{\max})^k \frac{k!}{(k - \tilde{r})!} x^{k - \tilde{r}}.
$$

Thus we have

$$
\begin{equation} \tag{4.15}
\begin{split}
(\ell + 1) \|X_{Z_{\ell+1,j}}\|_{\delta x} &\le \tilde{C} (\ell + 1) \sum_{m=0}^{\min\{j,N\}} \|X_{P_{m,h}}\|_{\delta x} \\
&\quad \times \sum_{\tilde{r}=1}^{\ell} \sum_{\substack{\ell_i>0, \ell_1+\ell_2+\dots+\ell_{\tilde{r}}=\ell \\ j_1+j_2+\dots+j_{\tilde{r}}=j-m}} \underbrace{\|X_{Z_{\ell_1,j_1}}\|_{\delta x} \dots \|X_{Z_{\ell_{\tilde{r}},j_{\tilde{r}}}}\|_{\delta x}}_{\le \max_{\ell_k \le \ell, j_k \le j} C_{\ell_k, j_k}^{\ell+j-m}} \\
&\le (\ell + 1) \tilde{C} \sum_{m=0}^{\min\{j,N\}} \|X_{P_{m,h} \circ R(hA)^*}\|_{\delta x} \max_{\ell_k \le \ell, j_k \le j} C_{\ell_k, j_k}^{\ell+j-m} \sum_{\tilde{r}=1}^{\ell} \sum_{\substack{\ell_i>0, \ell_1+\ell_2+\dots+\ell_{\tilde{r}}=\ell \\ j_1+j_2+\dots+j_{\tilde{r}}=j-m}} 1
\end{split}
\end{equation}
$$

where $\tilde{C} = \sup_{|x|<2\pi-2\tilde{\epsilon}, \tilde{r}=1,...,l} |f^{(\tilde{r})}(x)|$, provided that

$$
(4.16) \qquad \|X_{A_0}\|_{\delta x} \le \frac{\pi - \tilde{\epsilon}}{d_{\max}}.
$$
By Lemma 4.1, and since $d_{\max} \le 2r(\ell + j + 1) + 1$, the CFL condition (4.13) is sufficient to guarantee (4.16). Thus we have by (4.15) that

$$ \|X_{Z_{\ell+1,j}}\|_{\delta x} \le C_{\ell+1,j} \max_{0 \le m \le \min\{j,N\}} \|X_{P_{m,h}}\|_{\delta x}, $$

for some constant $C_{\ell+1,j} > 0$ independent of $K, \delta x$. Moreover, by construction (cf. (3.20)) we have

$$ \|X_{P_{m,h} \circ R(hA)^*}(u)\|_{\delta x} \le C_m \|u\|_{\delta x}^{2r(m+1)+1}, \quad \forall u \in V_{\delta x}, $$

for some constants $C_m > 0$ independent of $\delta x, K$. The result thus immediately follows from Lemma 2.11. □

## 4.2. Final estimates and modified energy.
We are now in a position to show that the modified energy defined in the above is indeed an accurate representation of the midpoint rule, i.e. we are able to prove the following central result of the present work.

**Theorem 4.3.** Let $u^{n+1} = \varphi_h(u^n)$ be the application defined by the midpoint rule (2.5). Fix $N \in \mathbb{N}$. Then there exists $h_0$ such that for $h \le h_0$ and all $K$ and $\delta x$ satisfying the condition CFL condition (4.13) then there exists

$$ (4.17) \qquad H_h = A_0 + B_h, $$

with

$$ A_0(u) = \bar{u}^T \frac{2}{h} \arctan \left( \frac{hA}{2} \right) u, $$

and

$$ (4.18) \qquad B_h(u) = \sum_{k=0}^{N} h^k B_{k,h}, $$

where the $B_{k,h}$ are polynomials of degree $2r(k+1)+2$, and such that there is a constant $C_N > 0$ independent of $K, \delta x > 0$ such that $\forall R \le 1$, for all $u \in B_{\delta x}(R)$, then

$$ (4.19) \qquad \|\varphi_h(u) - \Phi_{H_h}^h(u)\|_{\delta x} \le h^{N+2} C_N R^{2r(N+2)+1}, $$

i.e. $H_h$ is the modified energy up to $O(h^{N+2})$.

*Proof of Theorem 4.3.* The proof of this theorem is largely the same as Theorem 4.2 in [13], but with the double expansion in $t$ and $\varepsilon$. Let us stress the main points:

(i) **Reduction to the *splitted* form.** The combination of Proposition 3.1 and Proposition (3.3) shows that the Theorem is equivalent to proving the same result with the flow

$$ \Phi_{A_0}^1 \circ \Phi_{P_{h,\varepsilon}^{(N)}}^t $$

instead of the midpoint rule $\varphi_h$.

(ii) **Construction of the *modified* Hamiltonian.** We define the truncated expansion $Z^{(N)}(t, \varepsilon)$ of (4.11) by

$$ Z^{(N)}(t, \varepsilon) := \sum_{\ell+j=0}^{N} t^\ell \varepsilon^j Z_{\ell,j,h}. $$
and we set $H_h = Z^{(N)}(h, h)$. Note that we have indeed $A_0 = Z_{0,0,h}$ and that the

$$B_{k,h} := \sum_{\ell+j=k} Z_{k,j,h}$$

are homogeneous polynomials of order $2r(N+1)+2$.

(iii) Flow of the modified equation. We define the Hamiltonian $Q^{(N)}(t, \varepsilon)$ by the formula

$$\sum_{k \ge 0} \frac{1}{(k+1)!} \mathrm{ad}_{X_{Z^{(N)}(t, \varepsilon)}}^k \partial_t X_{Z^{(N)}(t, \varepsilon)} = X_{Q^{(N)}(t, \varepsilon)}.$$

In view of (4.10) and the construction of the $Z$, we can prove that

$$(4.20) \qquad Q^{(N)}(t, \varepsilon)(u) = P_{h,\varepsilon}^{(N)} \circ R(hA)^*(u) + O(h^{N+1}R^{2r(N+2)+2}),$$

for $|t| \le h$, $|\varepsilon| \le h$, and $u \in B_{\delta x}(R)$. Similarly to Lemma 4.3 in [13], we can prove that for $t, \varepsilon \le \tau_0(N)$ and for $R \le 1$, small enough, we have (cf. the formal equation (4.5))

$$\frac{d}{dt} \Phi_{Z^{(N)}(t, \varepsilon)}(u) = X_{Q^{(N)}(t, \varepsilon)} \circ \Phi_{Z^{(N)}(t, \varepsilon)}(u)$$

for $u \in B_{\delta x}(R)$.

(iv) Estimating the difference. We set

$$v(t) := \Phi_{A_0}^1 \circ \Phi_{P_{h,\varepsilon}^{(N)}}^t - \Phi_{Z^{(N)}(t,\varepsilon)}^1,$$

and we want to control this quantity uniformly in $|t| \le h$ and $|\varepsilon| \le h$. To do this, we write

$$\begin{aligned} \|v(t)\|_{\delta x} &\le \int_0^t \left\| X_{P_{h,\varepsilon}^{(N)} \circ R(hA)^*} \circ \Phi_{A_0}^1 \circ \Phi_{P_{h,\varepsilon}^{(N)}}^t (v(s)) - X_{Q(t,\varepsilon)} \circ \Phi_{Z^{(N)}(t,\varepsilon)}^1 (v(s)) \right\|_{\delta x} \mathrm{d}s \\ &\le \underbrace{\int_0^t \left\| X_{P_{h,\varepsilon}^{(N)} \circ R(hA)^*} \circ \Phi_{Z^{(N)}(t,\varepsilon)}^1 (v(s)) - X_{Q(t,\varepsilon)} \circ \Phi_{Z^{(N)}(t,\varepsilon)}^1 (v(s)) \right\|_{\delta x} \mathrm{d}s}_{=:A} \\ &\quad + \underbrace{\int_0^t \left\| X_{P_{h,\varepsilon}^{(N)} \circ R(hA)^*} \circ \Phi_{A_0}^1 \circ \Phi_{P_{h,\varepsilon}^{(N)}}^t (v(s)) - X_{P_{h,\varepsilon}^{(N)} \circ R(hA)^*} \circ \Phi_{Z^{(N)}(t,\varepsilon)}^1 (v(s)) \right\|_{\delta x} \mathrm{d}s}_{=:B} \end{aligned}$$

Then

* A is estimated using (4.20) and we have that

$$A \le C_N h^{N+2} R^{2r(N+2)+2}$$

for $u \in B_{\delta x}(R)$ and $R \le 1$, as in this case $v(s) \in B_{\delta x}(2R)$ for $s \in [0, h]$.

* B is estimated as the Hamiltonian vector field $X_{P_{h,\varepsilon}^{(N)} \circ R(hA)^*}$ is locally Lipschitz (with constant depending on $N$, but uniformly in $|t| \le h$ and $|\varepsilon| \le h$.

The conclusion then follows from classical Gronwall estimates [16, 18].

□
## 5. APPLICATION TO LONG-TIME STABILITY

The existence of the modified energy guarantees the stability of the numerical scheme. The previous result then implies the existence and stability of numerical solitons, as in [4], as well as normal form result as in [14, 15, 1] yielding preservation of the Fourier modes of the solution over long times for small initial data.

As a simple example of an application of Theorem 4.3, we provide here an almost global existence result for small initial date in $V_{\delta x}$. It is based on the idea that for small data, the energy (4.17) controls the norm $|| \cdot ||_{\delta x}$ for small initial data. This is a consequence of the following result:

**Lemma 5.1.** Assume $\delta x$ and $h$ satisfy the condition $h \le C\delta x^2$. Then we have

$$A_0(u) = \bar{u}^T \frac{2}{h} \arctan \left( \frac{hA}{2} \right) u \ge c \bar{u}^T A u,$$

where the constant $c$ only depends on $C$.

*Proof.* We have the existence of $c$ such that for all $x \le C$,

$$\arctan(x) \ge cx.$$

Then we know that there exists a unitary matrix $U$ such that

$$A = U^{-1}DU \quad \text{with} \quad D = \text{diag}(\lambda_j)$$

where the $\lambda_j$ are given by (4.3). Note that we thus have $\frac{1}{2}h\lambda_j \le C$ by assumption. Then we have, with $v = Uu$,

$$\begin{aligned} A_0(u) &= \bar{v}^T \frac{2}{h} \arctan \left( \frac{hD}{2} \right) v \\ &= \sum_j |v_j|^2 \frac{2}{h} \arctan \left( \frac{h\lambda_j}{2} \right) \\ &\ge \sum_j c\lambda_j |v_j|^2 = \bar{u}^T A u, \end{aligned}$$

and this shows the result.

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

**Theorem 5.2.** Let $u^{n+1} = \varphi_h(u^n)$ be the map defined by the midpoint rule (2.5). Fix $N \in \mathbb{N}$ and $\kappa \in (0, \frac{1}{2})$. Then there exists $\epsilon_0$ and $h_0$ such that for $h \le h_0$ and all $K$ and $\delta x$ satisfying the CFL condition (4.13) and all $\epsilon \le \epsilon_0$, we have

$$\|u^0\|_{\delta x} = \epsilon \implies \|u^n\|_{\delta x} \le \epsilon^{1-\kappa}, \quad \text{for } nh \le (h\epsilon^{2r(1-\kappa)})^{-N},$$

i.e. the numerical solution is stable in $|| \cdot ||_{\delta x}$ for long times.

*Proof.* Let $H_h$ be the energy (4.17) given by Theorem 4.3. Assume that $u^n \in B_{\delta x}(\epsilon^{1-\kappa})$. Then we have using (4.19) and the fact that $H_h(\Phi_{H_h}^h(u)) = H_h(u)$ for all $u$,

$$\begin{aligned} |H_h(u^{n+1}) - H_h(u^n)| &\le |H_h(\varphi_h(u^n)) - H_h(\Phi_{H_h}^h(u^n))| + |H_h(\Phi_{H_h}^h(u^n)) - H_h(u^n)| \\ &\le C_N h^{N+2} \epsilon^{(2r(N+2)+1)(1-\kappa)}. \end{aligned}$$

This shows that as long as $u^n \in B_{\delta x}(\epsilon^{1-\kappa})$, we have

$$|H_h(u^{n+1}) - H_h(u^0)| \le (nh)C_N h^{N+1} \epsilon^{(2r(N+2)+1)(1-\kappa)}.$$
But using (4.18), we have

$$ egin{aligned} |A_0(u^n) - A_0(u^0)| &\le |H_h(u^{n+1}) - H_h(u^0)| + C_N \sum_{k=0}^{N} \epsilon^{(2r(k+1)+2)(1-\kappa)} \\ &\le C_N \epsilon^{(2r+2)(1-\kappa)} [1 + (nh)h^{N+1} \epsilon^{2r(1-\kappa)N}] \le 2C_N \epsilon^{(2r+2)(1-\kappa)} \end{aligned} $$

for $(nh)h^{N+1}\epsilon^{2rN(1-\kappa)} \le 1$. As the $L^2$-norm $N(u)$ is preserved (see (2.4)), we deduce that as long as $u^n \in B_{\delta x}(\epsilon^{1-\kappa})$ and for $nh \le (h\epsilon^{2r(1-\kappa)})^{-N}$ we have, using the previous Lemma,

$$ egin{aligned} \|u^\delta\|_{\delta x}^2 &\le \frac{1}{c} (A_0(u^n) + N(u^n)) \le C_N (A_0(u^0) + N(u^0) + \epsilon^{(2r+2)(1-\kappa)}) \\ &\le C_N (\|u^0\|_{\delta x}^2 + \epsilon^{(2r+2)(1-\kappa)}) \le C_N \epsilon^2 \le \epsilon^{2(1-\kappa)} \end{aligned} $$

as $(r+1)(1-\kappa) > 1$, and for $\epsilon$ small enough, and with possible changes of the constant $C_N$ between the lines of the previous estimates. This shows the desired result. $\square$

## REFERENCES

[1] Charbella Abou Khalil and Joackim Bernier. Almost conservation of the harmonic actions for fully discretized nonlinear Klein-Gordon equations at low regularity. https://arxiv.org/abs/2406.12363. 2024.

[2] U. M. Ascher and S. Reich. The midpoint scheme and variants for Hamiltonian systems: advantages and pitfalls. *SIAM J. Sci. Comput.*, 21:1054-1065, 1999.

[3] D. Bambusi and T. Penati. Continuous approximation of breathers in one-and two-dimensional DNLS lattices. *Nonlinearity*, 23(1):143, 2009.

[4] Dario Bambusi, Erwan Faou, and Benoît Grébert. Existence and stability of ground states for fully discrete approximations of the nonlinear Schrödinger equation. *Numerische Mathematik*, 123(3):461-492, 2013.

[5] Giancarlo Benettin and Antonio Giorgilli. On the Hamiltonian interpolation of near-to-the identity symplectic mappings with application to symplectic integration algorithms. *Journal of Statistical Physics*, 74:1117-1143, 1994.

[6] Joackim Bernier and Erwan Faou. Existence and stability of traveling waves for discrete nonlinear Schrödinger equations over long times. *SIAM J. Math. Anal.*, 51:1607-1656, 2019.

[7] Jacek Bochnak and Józef Siciak. Polynomials and multilinear mappings in topological vector-spaces. *Studia Mathematica*, 39:59-76, 1971.

[8] Elena Celledoni, David Cohen, and Brynjulf Owren. Symmetric exponential integrators with an application to the cubic Schrödinger equation. *Found. Comput. Math.*, 8:303-317, 2008.

[9] R. Courant, K. Friedrichs, and H. Lewy. Über die partiellen Differenzengleichungen der mathematischen Physik. *Mathematische Annalen*, 98:32-74, 1928.

[10] Arnaud Debussche and Erwan Faou. Modified energy for split-step methods applied to the linear Schrödinger equation. *SIAM Journal on Numerical Analysis*, 47(5):3705-3719, 2009.

[11] Erwan Faou. *Geometric numerical integration and Schrödinger equations*, volume 15. European Mathematical Society, 2012.

[12] Erwan Faou, Ludwig Gauckler, and Christian Lubich. Plane wave stability of the split-step Fourier method for the nonlinear Schrödinger equation. In *Forum of Mathematics, Sigma*, volume 2, page e5. Cambridge University Press, 2014.

[13] Erwan Faou and Benoît Grébert. Hamiltonian interpolation of splitting approximations for nonlinear PDEs. *Foundations of Computational Mathematics*, 11:381-415, 2011.

[14] Erwan Faou, Benoît Grébert, and Eric Paturel. Birkhoff normal form for splitting methods applied to semi linear Hamiltonian PDEs. Part I: Finite dimensional discretization. *Numer. Math.*, 114:429-458, 2010.
[15] Erwan Faou, Benoît Grébert, and Eric Paturel. Birkhoff normal form for splitting methods applied to semi linear Hamiltonian PDEs. Part II: Abstract splitting. *Numer. Math.*, 114:459–490, 2010.

[16] Thomas Hakon Gronwall. Note on the derivatives with respect to a parameter of the solutions of a system of differential equations. *Annals of Mathematics*, 20(4):292–296, 1919.

[17] E. Hairer, C. Lubich, and G. Wanner. *Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations*. Springer, 2013.

[18] Ralph Howard. The Gronwall Inequality. https://people.math.sc.edu/howard/Notes/gronwall.pdf.
Accessed: 23-08-2024.

[19] Panayotis G Kevrekidis. *The discrete nonlinear Schrödinger equation: mathematical analysis, numerical computations and physical perspectives*, volume 232. Springer Science & Business Media, 2009.

[20] Georg Maierhofer and Katharina Schratz. Bridging the gap: symplecticity and low regularity in Runge-Kutta resonance-based schemes. *Mathematics in Computation*, 2024.

[21] J.E. Marsden and T. Ratiu. *Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems, 2nd Edition*. Texts in Applied Mathematics. Springer New York, 2002.

[22] Ander Murua Uria. *Métodos simplécticos desarrollables en P-series*. PhD thesis, Universidad de Valladolid, 1995.

[23] Mei qing Zhang and Robert D. Skeel. Cheap implicit symplectic integrators. *Applied Numerical Mathematics*, 25:297–302, 1997.

[24] Ari Stern and Eitan Grinspun. Implicit-explicit variational integration of highly oscillatory problems. *Multiscale Modeling & Simulation*, 7(4):1779–1794, 2009.

[25] Y.-F. Tang. Formal energy of a symplectic scheme for Hamiltonian systems and its applications (I). *Computers & Mathematics with Applications*, 27(7):31–39, 1994.