Metastable energy strata in numerical discretizations of weakly nonlinear wave equations

The quadratic nonlinear wave equation on a one-dimensional torus with small initial values located in a single Fourier mode is considered. In this situation, the formation of metastable energy strata has recently been described and their long-time stability has been shown. The topic of the present paper is the correct reproduction of these metastable energy strata by a numerical method. For symplectic trigonometric integrators applied to the equation, it is shown that these energy strata are reproduced even on long time intervals in a qualitatively correct way.


Introduction
We consider the nonlinear wave equation on a one-dimensional torus, x ∈ T = R/(2πZ), with a positive Klein-Gordon parameter ρ. We assume that the initial value consists of a single Fourier mode and is small. The smallness of the initial value corresponds to a weakly nonlinear setting. In [15], it has been investigated how the energy, which is initially located in this single Fourier mode, is distributed among the other modes in the course of time. More precisely, the formation of energy strata has been shown that are persist on long time intervals, see Section 2 for a precise description of this result.
In the present paper, we discretize the nonlinear wave equation and answer the question whether this long-time property of the exact solution is inherited by the numerical method. As a numerical method, we consider a spectral collocation in space combined with a symplectic trigonometric integrator in time. We show that this numerical method in fact reproduces the energy strata of the exact solution even on a long time interval, provided that a certain numerical non-resonance condition is fulfilled, see Section 2 for a formulation of this main result.
The considered numerical method is already known to behave well on long time intervals with respect to preservation of regularity and near-conservation of actions [6,12], as well as near-conservation of energy and momentum [2,3,6,11]. Our result adds to this list an even more sophisticated long-time property of the exact solution that is reproduced in a qualitatively correct way. In comparison to previous results, the considered situation requires less control of interactions in the nonlinearity, which allows us to exclude numerical resonances under weaker restrictions on the time step-size than needed previously.
The proof of our main result is given in Sections 3-5. As for the exact solution, it is based on a modulated Fourier expansion in time, with a multitude of additional difficulties due to the discrete setting that will be described in detail.
Related questions have been studied for symplectic and non-symplectic methods applied to the nonlinear Schrödinger equation. In this case, an initial value consisting of a single Fourier mode yields a solution that continues to consist of a single Fourier mode for all times (plane wave solution). In particular and in contrast to the nonlinear wave equation considered in the present paper, there is no formation of energy strata in the other modes. The stability in numerical discretizations of these plane wave solutions under small perturbations of the initial value is an old question [22] and has been studied on short time intervals [4,7,19,20,22] and on long time intervals [10].

Statement of the main result 2.1 Metastable energy strata revisited
Writing the solution u = u(x, t) of the nonlinear wave equation (1) as a Fourier series j∈Z u j (t)e ijx with Fourier coefficients u j , the mode energies of the nonlinear wave equation (1) are given by E j (t) = 1 2 |ω j u j (t)| 2 + 1 2 |u j (t)| 2 , j ∈ Z, where ω j are the frequencies Note that E j = E −j for real-valued initial values (and hence real-valued solutions). We are interested in the evolution of these mode energies. Assuming that the initial values are small and concentrated in the first mode, an inspection of interaction of Fourier modes in the nonlinearity g(u) = u 2 suggests that the energy in the first mode is distributed among the other modes in a geometrically decaying way: at least for small times t.
The main result of [15] (Theorem 1 for θ = 1) states that these energy strata are in fact stable on long time intervals in the following sense: for fixed but arbitrary  K ≥ 2 and s > 1 2 , we have on the long time interval This can be written in compact and equivalent form as see also Figure 1 for an illustration. In addition to these bounds on all mode energies, the first mode energy is nearly conserved on this long time interval, |E 1 (t) − E 1 (0)| ≤ C 2 . All constants depend on K and s but not on , and they deteriorate as K or s grow. This result is valid for all but finitely many values of the parameter ρ ∈ [0, ρ 0 ] in the wave equation (1), which excludes some resonant situations. The question that we want to answer in the present paper is whether this stable behaviour on long time intervals is inherited by a typical structure-preserving numerical discretization of the nonlinear wave equation.
In these methods, the solution u = u(x, t) of (1) is approximated, at discrete times t n = nτ with the time step-size τ , by a trigonometric polynomial of degree M : The Fourier coefficients u n = (u n −M , . . . , u n M −1 ) T of this trigonometric polynomials are computed with the trigonometric integrator where Ω denotes the diagonal matrix containing the frequencies ω j , j = −M, . . . , M −1, of (3). In addition, * denotes the discrete convolution defined by which can be computed efficiently using the fast Fourier transform. The method is characterized by the diagonal filter matrices Ψ = ψ(τ Ω) and Φ = φ(τ Ω), which are computed from filter functions ψ and φ. These filter functions are assumed to be real-valued, bounded and even with ψ(0) = φ(0) = 1. The starting approximation is computed by and a velocity approximation by 2τ sinc(τ Ω)u n = u n+1 − u n−1 .
An error analysis of these methods has been given in [14]. We assume here that the considered trigonometric integrator is symplectic, see [18, Chap. XIII, Eq. (2.10)]. Examples for such filter functions are φ = 1 and ψ = sinc, which is the impulse method [16,21] or method of Deuflhard [8], as well as φ = sinc and ψ = sinc 2 , which is the mollified impulse method of García-Archilla, Sanz-Serna & Skeel [13]. Certain choices of filter functions leading to non-symplectic methods could also be handled using the transformation indicated in [5, Remark 3.2], but we do not pursue this here. We are interested in the mode energies (2) along the numerical solution (6). We denote them in the following by We note that E n j = E n −j for real-valued initial values (and hence numerical solutions that take real values in the collocation points x k = πk/M , k = −M, . . . , M − 1), and we set E n M = E n −M . As for the exact solution in Section 2.1, we consider initial values with

Metastable energy strata in trigonometric integrators
We now state our main result which says, roughly speaking, that trigonometric integrators (6) integrate the metastable energy strata in nonlinear wave equations (1) qualitatively correctly. For this result, we need a non-resonance condition on the time step-size τ and on the frequencies ω j , j = −M, . . . , M − 1, of (3). In the statement of this condition, we consider indices j ∈ {−M, . . . , M − 1} and vectors k = (k 0 , . . . , k M ) T of integers k l , and we write with the vector ω = (ω 0 , . . . , ω M ) T of frequencies. We fix 3 ≤ K ≤ M , and we denote, corresponding to [15], by K the set K = (j, k) : max(|j|, µ(k)) < 2K and k l = 0 for all l ≥ K (10) where j = (0, . . . , 0, 1, 0, . . . , 0) T is the |j|th unit vector and In comparison to the set K for the exact solution (see [15,Equation (11)]), we have to consider indices modulo 2M due to the discretization in space. In addition, we correct here a typo in [15,Equation (11)], where |j| ≥ K should be replaced by |j − r| ≥ K.
(12b) We emphasize that only k = − j is excluded in the estimate for sin( 1 2 τ (ω j + k · ω)) in the first condition (12a), whereas k = j and k = − j are both excluded in the estimate for this sine in the second condition (12b).
Under this non-resonance condition, we prove in this paper the following discrete analogon of the analytical result [15, Theorem 1] described in Section 2.1. Theorem 1. Fix an integer 3 ≤ K ≤ M and real numbers 0 ≤ ν < 1 2 and s > 1 2 , and assume that the non-resonance condition (12) holds for this ν. Then there exist 0 > 0 and positive constants c and C such that the mode energies (8) along trigonometric integrators (6) for nonlinear wave equations (1) with initial data (9) with 0 < ≤ 0 satisfy, over long times and −e(1) |E n 1 − E 0 1 | ≤ C 1−2ν . The constants c and C depend on K, ν and s, but not on and the discretization parameters τ and M . Remark 2. The result of Theorem 1 gets stronger for larger K (with worse constants, however), but also the assumption gets stronger, because the set K (see (10)) in the non-resonance condition (12) grows. We can thus also cover the (not so interesting) border case K = 2 as considered in [15] provided that the non-resonance assumption holds for K = 3. Without requiring the non-resonance condition with K = 3, our techniques of proof can be used to show that Theorem 1 holds for K = 2, if we replace −e(1) |E n The same remark applies to [15].
The proof of Theorem 1 is given in Sections 3-5 below. The structure of the proof is similar to the structure of the corresponding proof for the exact solution given in [15], with a multitude of additional difficulties due to the discrete setting.
The main difference of this result in comparison with the corresponding result [15,Theorem 1] for the exact solution is the required non-resonance condition. The non-resonance condition (12) excludes two types of resonances. First, it requires that ω j ± k · ω is bounded away from zero for (j, k) ∈ K with k = ∓ j . This is the non-resonance condition for the exact solution, which can be shown to hold for all except finitely many values of ρ in a fixed interval, see [15,Section 3.2]. In addition, the non-resonance condition (12) requires also that products τ (ω j ± k · ω) are bounded away from nonzero integer multiples of 2π. If this latter condition is violated, we observe numerical resonances as illustrated in the following section.
Numerical resonances can be typically excluded under some CFL-type step-size restriction on the time step-size. A feature of the considered situation is that they can be excluded under less restrictive assumptions than in previous studies on the long-time behaviour of trigonometric integrators for nonlinear wave equations [6,12]. In particular, there are no numerical resonances on the time interval 0 ≤ t ≤ −K/2 under the CFL-type step-size restriction for some constant c (we may take ν = 0 here). This follows from the structure of the set K of (10) which enters the non-resonance condition: this set allows only a single large frequency ω l ≤ M √ 1 + ρ in the linear combination k · ω, and it allows to bound the remaining part of k·ω by 2K √ 1 + ρ (using the definition of µ and the bounds on µ in K). Under the above CFL condition (13) and using ω j ≤ M √ 1 + ρ, the arguments 1 2 τ (ω j ± k · ω) of the sines in the non-resonance condition (12) are thus bounded by c < π, which excludes numerical resonances (but not analytical resonances).
In the situation of [6,12], the stronger CFL-type step-size restriction of the form τ M K ≤ c < π has to be used to exclude numerical resonances.
For nonzero ν, the numerical non-resonance condition (12) resembles the one used in [6], see Equations (23) and (24) therein. We note that the reduction (24) there is of no relevance in our context because of the structure of the set K. We also note that the number of indices (j, k) that have to satisfy the stronger condition (12b), which does not appear in [6], is bounded independently of the spatial discretization parameter M .
There are several possible extensions of Theorem 1. First, we could extend the result, as in [15], to an energy profile e with e(l) = K + (|l| − K)(1 − θ) instead of e(l) = K for |l| ≥ K, where 0 < θ ≤ 1. This shows that also the mode-energies E l for l > K decay geometrically, but only with a smaller power of close to 1. As in [15], the time-scale is then restricted to 0 ≤ t n ≤ −θK(1−2ν)/2 . Second, the result holds for general nonlinearities g(u) instead of u 2 if g is realanalytic near 0 and g(0) = g (0) = 0. In addition, stronger estimates hold if further derivatives of g vanish at 0.
Finally, the result can be extended to more general initial energy profiles, for example to the situation where a whole band E 0 l , |l| ≤ B, of initial mode energies is of order .

Numerical experiment
We consider the nonlinear wave equation (1) with ρ = √ 3 and with initial value satisfying (9) for = 10 −3 . We apply the trigonometric integrator (6) with filter functions φ = 1 and ψ = sinc (the impulse method or method of Deuflhard) to the equation. We take M = 2 5 for the spectral collocation in space, and we use three different time step-sizes for the discretization in time 1 .
The first time step-size is τ = 0.05. For this time step-size, the trigonometric integrator integrates the geometrically decaying energy strata qualitatively correctly, even on long time intervals, see Figure 2 (every tenth time step is plotted there). The time step-size fulfills the CFL-type condition (13) with K = 6, under which numerical resonances in the non-resonance condition (12) of Theorem 1 can be excluded.
The second time step-size is τ = 2π/(ω 1 + ω 6 + ω 7 ) ≈ 0.4212. It is chosen in such a way that it does not satisfy the non-resonance condition (12). The resulting numerical resonance becomes apparent in Figure 3, where we observe in particular an exchange among the first, sixth and seventh mode. For better visibility, we have plotted max m=0,...,99 E n+m l for n = 0, 100, 200, . . . instead of E n l for n = 0, 1, 2, . . . in Figure 3. In the same figure, a similar behaviour can be observed for the resonant time step-size τ = 2π/(−ω 1 + ω 6 + ω 7 ).

Approximation ansatz
We will approximate the numerical solution by a modulated Fourier expansion, where K j = {k : (j, k) ∈ K} with the set K of (10). We require that the modulation functions z k j are polynomials. In contrast to the modulated Fourier expansion of [15, Equation (7)] for the exact solution, the modulation functions considered here vary on a slow time-scale ν/2 t, where ν is the parameter from the non-resonance condition (12). With this slow time scale, derivatives of z k j ( ν/2 t) with respect to t carry additional factors ν/2 , which will be used to compensate for this factor in the weak non-resonance condition (12a).
Having in mind that u n j , j = −M, . . . , M − 1, are the Fourier coefficients of a trigonometric polynomial of degree M , we write in the following u n M = u n −M and z k M = z k −M . Inserting the ansatz (14) into the trigonometric integrator (6a) and comparing the coefficients of e i(k·ω)t , yields a set of equations for the modulation functions z = (z k j ) (j,k)∈K , for which solutions are to be constructed up to a small defect d k j : with the extended potential with Φz = (φ(τ ω j )z k j ) (j,k)∈K and with ∇ −k −j denoting the partial derivative with respect to z −k −j . In (16) and in the following we denote by ≡ the congruence modulo 2M .
With the modulated Fourier expansion (14) at hand, the formula (6c) for the velocity approximation leads to an approximation ofu n j by a modulated Fourier expansion: Finally, we get conditions from the fact that the ansatz (14) should satisfy the initial condition: where (17) was used for the derivation of the second equation.

Norms
where e(j) is the energy profile of (5) and σ j are the weights of (4). For ν = 0, this is the Sobolev H s -norm of the corresponding trigonometric polynomial The rescaling −2e(j)ν originates from the non-resonance condition (12a) that introduces ν/2 .
We will make frequent use of the fact that this norm behaves well with respect to the discrete convolution, This follows from the corresponding property of the Sobolev H s -norm, see, for example, [17,Lemma 4.2], and the fact that e(j) ≤ e(j 1 ) + e(j 2 ) for j ≡ j 1 + j 2 . Similarly, we also have In particular, we have for ν = 0 in (19)

The modulated Fourier expansion on a short time interval
Instead of initial mode energies (9), we consider the more general choice of initial mode energies, which is the expected situation at later times (see Theorem 1). We then have the following discrete counterpart of [15,Theorem 3], whose proof will be given in Section 4.
Theorem 3. Under the assumptions of Theorem 1, but with (22) instead of (9), the numerical solution u n = (u n −M , . . . , u n M −1 ) T admits an expansion (10), where the coefficient functions z k j are polynomials satisfying z −k −j = z k j , and where the remainder term r n = (r n −M , . . . , r n M −1 ) T and the corresponding remainder termṙ n = (ṙ n −M , . . . ,ṙ n M −1 ) T in the velocity approximation (17) are bounded by The constant C is independent of , τ and M , but depends on K, on ν and γ of (12), on ρ of (1), on s of (4) and on C 0 of (22).

Almost-invariant energies
We now derive almost-invariant energies of the modulation system (15) which enable us to consider long time intervals. The derivation of these almost-invariant energies is similar as in the case of the exact solution, see [15, Section 3.5], but it now leads to discrete almost-invariant energies that involve additionally the time step-size τ and are related to those of [6]. Let Our aim is to prove that which shows that, for small defect d k j , the quantity E l is almost invariant. To see this almost-invariance, we use, as in the case of the exact solution, that the extended potential U of (16) is invariant under the transformation S λ (θ)z := (e i(k·λ)θ z k j ) (j,k)∈K for θ ∈ R and real vectors λ = (λ 0 , . . . , λ M ). This invariance shows that We then multiply this equation with τ /2, we replace −∇ −k −j U with the help of (15), and we use the symplecticity (7) of the method. This yields In this relation, we use that which follows from the symmetry in j (recall that z k M = z k −M ) and the asymmetry in k. Choosing λ = l with l = 0, . . . , M and using the definition (23) of E l , we then end up with (24).
Note that, due to the appearance of the factor k l in E l of (23), only modulation functions z k j with k l = 0 are actually relevant in E l (recall that k l is the lth entry of the vector k). This will be used all the time in Section 5 below. In particular, we will use there that only very specific vectors k with k l = 0 for l ≥ K appear in the set K of (10).
Using a Taylor expansion of z k j ( ν/2 (t + τ )) and the property (25), we can derive the following alternative form of E l : Using (25), we see that this quantity coincides, for ν = 0 and in the limits τ → 0 and M → ∞, with the almost-invariant energy of [15,Equation (23)] for the exact solution. We therefore also call E l (t) an almost-invariant energy. In Sections 4 and 5 below, we prove the following discrete counterparts of [15, for this new almost-invariant energy.
Theorem 4 (Almost-invariant energies controlled by mode energies). Under the conditions of Theorem 3 we have, for 0 ≤ t ≤ 1, where C 0 is independent of , τ and M and depends on the initial data only through the constant C 0 of (22).

Theorem 5 (Variation of almost-invariant energies).
Under the conditions of Theorem 3 we have, for 0 ≤ t n ≤ 1, where C is independent of , τ and M and depends on the initial data only through the constant C 0 of (22).
At time t = 1, for which we assume without loss of generality that t N = N τ = 1 for some N , we consider a new modulated Fourier expansion leading to new almostinvariant energies.
Theorem 6 (Transitions in the almost-invariant energies). Let the conditions of Theorem 3 be fulfilled. Let z( ν/2 t) = (z k j ( ν/2 t)) (j,k)∈K for 0 ≤ t ≤ 1 be the coefficient  functions as in Theorem 3 with corresponding almost-invariant energies E l (t) for initial data (u 0 ,u 0 ). Let further z( ν/2 t) = ( z k j ( ν/2 t)) (j,k)∈K be the coefficient functions and E(t) the corresponding almost-invariants of the modulated Fourier expansion for 0 ≤ t ≤ 1 corresponding to initial data (u N ,u N ) with t N = N τ = 1, constructed as in Theorem 3. If (u N ,u N ) also satisfies the bound (22) where C is independent of , τ and M and depends on the initial data only through the constant C 0 of (22). and −e(1) |E n 1 − E 1 (t n )| ≤ C 1−2ν , where C depends on C 0 in (27), but is independent of , τ , M and C 0 of (22) if 1−2ν is sufficiently small.

From short to long time intervals: Proof of Theorem 1
Based on Theorems 3-7, the proof of Theorem 1 is the same as in the case of the exact solution, see [15, Section 3.6]: Theorems 3-5 and 7 yield the statement of Theorem 1 on a short time interval 0 ≤ t n ≤ 1. Theorem 6 can be used to patch many of these short time intervals together, on which the almost-invariant energies E l are still well preserved (Theorems 5 and 6) and allow to control the mode energies E l (Theorem 7). A schematic overview of the proof is given in Figure 4.

Construction of a modulated Fourier expansion: Proofs of Theorems 3 and 4
Throughout this section, we work under the assumptions of Theorem 3. In addition, we let δ = 1/2 , such that u 0 j ,u 0 j = O(δ e(l) ). Moreover, we write for an inequality up to a factor that is independent of δ = 1/2 , τ and M .

Expansion of the modulation functions
The construction of modulation functions z k j of the modulated Fourier expansion (14) is based on an expansion with polynomials z k j,m = z k j,m (δ ν t) and with m(j, k) = max e(j), max l : k l =0 e(l) .
This can be motivated by the fact that on the one hand we expect z k j = O(δ m(j,k) ) from the analysis of the exact solution in [15], but on the other hand we do not expect an expansion in true powers of δ as in [15] because of the non-resonance condition (12) involving δ ν .
Our goal now is to derive equations for the modulation coefficient functions z k j,m in the ansatz (28). Inserting this ansatz into (15), expanding z k j (δ ν (t ± τ )) in a Taylor series and requiring that powers of δ 1−2ν agree on both sides yields, neglecting the defect d k j , 4s j −k s j +k z k j,m + 2iτ δ ν s 2kż k j,m + τ 2 δ 2ν c 2kz k j,m + . . . = τ 2 ψ(τ ω j ) Here, we use the notation s k = sin( τ 2 k · ω) and c k = cos( τ 2 k · ω) and the fact that cos x − cos y = 2 sin( y−x 2 ) sin( y+x 2 ). All functions in (29) are evaluated at δ ν t and dots on the functions z k j,m denote derivatives with respect to the slow time scale δ ν t. Note that m(j, k) ≤ m(j 1 , k 1 ) + m(j 2 , k 2 ) if k = k 1 + k 2 and j ≡ j 1 + j 2 mod 2M , and hence the power of δ on the right-hand side of (29) is small. We recall that ≡ denotes the congruence modulo 2M .

Remark 8.
For ν = 0 we recover in (29), after division by τ 2 and in the limit τ → 0, the system of equations which was used in the case of the exact solution, see [15,Equation (28)]. For ν > 0, the above construction becomes significantly more involved than there. The reason is that the non-resonance condition (12a) only allows us to bound the factor s j −k s j +k on the left-hand side of (29) from below by γ 2 τ 2 δ 2ν . As we will see in the proof of Lemma 10 below, we can compensate this possibly small factor with an additional δ 2ν on the right-hand side (together with τ 2 ), which we gain from the special choice (28) as ansatz for z k j .
In addition to (29), we get from condition (18) that where we use again a Taylor expansion of the modulation functions.

Construction of modulation functions
We construct polynomial modulation functions of the form (28) by solving (29) and (30) consecutively for m = 1, 2, . . . , 2K − 1. For convenience, we set z k j,m = 0 for m < m(j, k). Assuming that we have computed all polynomials z k j,m for m < m, the construction relies on the observation that only already computed functions z k j ,m appear on the right-hand side of (29). This equation is thus, for z(s) = z k j,m (s), of the form with coefficients α 0 , . . . , α L and a polynomial p.
For k ∈ K j with k = ± j , the coefficient α 0 in this equation is nonzero (by the non-resonance condition (12)), and the unique polynomial solution of this equation is given by The modulation coefficient functions constructed in this way are called off-diagonal modulation coefficient functions.
For k = ± j , the coefficient α 0 is zero, and the polynomial solutions of (31) are given by In (33a), the initial value z(0) is still a free parameter. We use (30) to fix it. Adding and subtracting the two equations of (30), after multiplying the first equation with 2iω j and the second one with ω j /s 2 j = 1/(τ sinc(τ ω j )), gives us Note again that this equation becomes for ν = 0 and in the limit τ → 0 the corresponding equation for the exact solution, see [15,Equation (31)]. The modulation coefficient functions constructed with (33) and (34) are called diagonal modulation coefficient functions.
Example 9. We illustrate the above construction for the first two nontrivial values of m. For simplicity, we restrict to the case ν = 0 and φ = 1 (and hence ψ = sinc by (7)). For m = 1, we observe that the polynomial p in (31), which is the right-hand side of (29), vanishes for all j and all k. The above construction (32) thus yields For m = 2, we observe that the polynomial p in (31), which is the right-hand side of (29), is non-zero only for j ∈ {±2, 0} and k = {±2 1 , 0}. In particular, we get from (32) For j = 2 and k = 2 1 , the polynomial p is non-zero but constant, and we get from (32) the non-zero off-diagonal modulation coefficient function  to the numerical solution u n j . For τ = 0.05 and ρ = √ 3, the quality of this approximation is illustrated in Figure 5 for initial values satisfying (9) with different values of . We observe that the approximation is exact for n = 0 (by construction) and that the approximation is of order δ 3 for j = 1 on a time interval of order one. We also observe an improved approximation order of δ 4 for j = 0 and j = 2. This can be explained with the construction for m = 3, which yields z k 2,3 = z k 0,3 ≡ 0 for all k.
Going beyond this example, we verify as in [15,Section 4.1] the following properties of the constructed functions z k j,m . Using e(j) ≤ e(j 1 ) + e(j 2 ) for j ≡ j 1 + j 2 , µ(k) ≤ µ(k 1 ) + µ(k 2 ) for k = k 1 + k 2 and e(j) = µ(k) < µ(k 1 ) + µ(k 2 ) for k = ± j = k 1 + k 2 with µ defined in (11), we see that they are polynomials of degree deg z k j,m ≤ m − max(e(j), µ(k)), where a negative degree corresponds to the zero polynomial. Moreover, the polynomial z k j,m can be different from the zero polynomial only in the two cases case 1: |j| ≤ m, µ(k) ≤ m, k l = 0 for l ≥ min(m + 1, K), This explains how the set K of (10) is built up. The two cases follow from the decomposition of the last line in (29), where the first term is only present for m > K: the functions z k i ji,mi in the second sum and the function z k 2 j2,l in the first sum belong inductively to the first case above, whereas the function z k 1 j1,m−l in the first sum belongs either to the first or to the second case, leading to a function z k j,m belonging to the first or second case, respectively. In addition, we have z k j,m = z −k −j,m and z ± l j,e(l) = 0 if |l| = |j|.

Bounds of the modulation functions
For polynomials z = z(s), we introduce the norm It has the properties that |z · w| t ≤ |z| t · |w| t and |v| t ≤ deg(v) |v| t . With this norm, we have the following discrete counterpart of [15,Lemma 1]. In contrast to the analytical result in [15], we gain a factor |sinc(τ ω j )| −1 in the estimates, which comes from the filter ψ in the trigonometric integrator (6) and will be needed later on to deal with this factor in the almost-invariant energies E l of (23).
This shows that off-diagonal modulation coefficient functions z k j,m are bounded by a constant, since only modulation coefficient functions z k j ,m with m < m can appear on right-hand side of the above estimate for z k j,m . To get the summed estimate of the lemma, we use γ k j γ k 1 j1 γ k 2 j2 for k = k 1 + k 2 and j ≡ j 1 + j 2 , which yields The first estimate then follows inductively from the algebra property (21) and 1 ≤ |sinc(τ ω j )| −1 .
(b) In the case k = ± j , the absolute value of the coefficient α 1 is bounded from below by γτ 2 δ 2ν by the (weaker) non-resonance condition (12a). Also in this case, we have an additional factor δ 2ν in the polynomial on the right-hand side of (31) since m(j, k) = e(j) and m(j 1 , k 1 ) ≥ e(j) or m(j 2 , k 2 ) ≥ e(j) for k = ± j = k 1 + k 2 . This shows that on the modulation functions z k j themselves. In particular, we have by the definition of m(j, k).

Proof of Theorem 4
The proof of Theorem 4 relies on the alternative form (26) of the almost-invariant energies E l (t) of (23). Note again that only modulation functions z k j with k l = 0 can appear in E l (t) due to the factor k l . By (41), this argument shows that E l (t) = O(δ 2e(l) ) has the claimed order in δ.
To get the precise estimate of Theorem 4, we multiply (26) with σ l δ −2e(l) , we use e(l) ≤ m(j, k) for k l = 0, and we sum over l, which yields by the definition of the norm | · | t . For the first term on the right-hand side of this estimate, we use and |k · ω| ≤ γ k j , we use that the 1 -norm dominates the 2 -norm (more precisely, k∈Kj |v k j | 2 ≤ ( k∈Kj |v k j |) 2 ), and we use (39) for k = ± j (the two sincs cancel in this case) and (40) for k = ± j . This shows that this first term is bounded by a constant. To show that also the second term is bounded by a constant, we apply in addition the non-resonance condition (12a) to bound δ ν /|sinc(τ ω j )| and the Cauchy-Schwarz inequality.

Bounds of the defect
When constructing z k j with the expansion (28), the defect d k j , k ∈ K j , in (15) is given by Although we consider system (15) defining the defect d k j only for k ∈ K j , we use this formula to define d k j also for k / ∈ K j . This will be helpful below.
Proof. The result follows with Lemma 10 and the arguments used in its proof.
On a higher level, the approximation u n = ( u n −M , . . . , u n M −1 ) T of (14) to the numerical solution u n then has a defect e n = (e n −M , . . . , e n M −1 ) T when inserted into the numerical method (6a): By construction, this defect is given by where the defect d k j and the modulation functions are evaluated at δ ν t n .
Proof. We use that the set K j was constructed in such a way that k belongs to K j if k = k 1 + k 2 , j ≡ j 1 + j 2 and m 1 + m 2 ≤ 2K − 1 for k 1 , k 2 , j 1 , j 2 , m 1 , m 2 with z k 1 j1,m1 = 0 and z k 2 j2,m2 = 0. Using (43) to define d k j also for k / ∈ K j , this yields the compact form e n j = k∈Z M +1 d k j (δ ν t n )e i(k·ω)tn for the defect in (44). From the definition of the norm · (see Section 3.2), we thus get Using e(j) ≤ m(j, k) and ω j ≤ γ k j , the statement of the lemma thus follows from Lemma 11.

Bounds of the remainder
With the constructed modulation functions z k j , we get the approximation u n j of (14) to the numerical solution u n j . From (17), we also get an approximation˙ u n j tou n j . In the following lemma, we establish a bound for the approximation error.
Lemma 13. The remainders u n j − u n j andu n j −˙ u n j satisfy, for 0 ≤ t n ≤ 1, Proof. (a) In a first step, we write the numerical scheme (6) and its approximation by a modulated Fourier expansion in one-step form. The method in one-step form reads τ Ψg n cos(τ Ω)Φg n + Φg n+1 (45) with g n = (Φu n ) * (Φu n ), see [18,Section XIII.2.2]. The first line is obtained by adding (6a) and (6c), and the second line is obtained by subtracting these equations with n + 1 instead of n and using the first line to replace u n+1 as well as the symplecticity (7). In the same way, we derive for the approximations u n and˙ u n , which satisfy (6c) exactly (by definition (17)) and (6a) up to a small defect given by (44), τ Ψ g n cos(τ Ω)Φ g n + Φ g n+1 (46) with g n = g n + e n .
(b) In the following, we write for a norm that is equivalent to the norm considered in the statement of the lemma. By induction on n, we prove that the numerical solutions stays small in this norm, (22), provided that δ 1−2ν is sufficiently small. This is true for n = 0 by the choice of the initial value (22). For n > 0, we first note that, by (45), Ωu n ≤ Ωu n−1 + u n−1 + 1 2 τ 2 ΩΨg n−1 , and hence Ωu n ≤ Cδ 1−2ν by induction, by (7), by the boundedness of φ and by the algebra property (20). We then observe that the matrix appearing in the one-step formulation (45) of the method preserves the norm (47). Moreover, by the algebra properties (19) and (20), the second term in the one-step formulation (45) can be estimated in the norm (47) by τ δ 2(1−2ν) up to a constant. This implies the stated bound of (u n ,u n ) . Similarly, we get for the modulated Fourier expansion the bound To get this bound, we just repeat the presented argument for the numerical solution given by (45) for the modulated Fourier expansion which satisfies the very similar relation (46). The only difference is the additional term e n in (46), which can be estimated with Lemma 12.
(c) For the difference (u n − u n ,u n −˙ u n ), we subtract the above one-step formulations, and then proceed as in the proof of the smallness of (u n ,u n ) and ( u n ,˙ u n ) in (b). We use the smallness of (u n ,u n ) and ( u n ,˙ u n ) together with the algebra property (19) and Lemma 12 to control the difference g n − g n = −e n . The mentioned lemma on the defect introduces the small parameter δ 2K(1−2ν) .
With the results proven in this section, the proofs of Theorems 3 and 4 are complete. τ k l ω l φ(τ ω j )z −k −j (δ ν tñ)d k j (δ ν tñ) for the lth almost-invariant energy E l . Note again that only indices k with k l = 0 are relevant in the above sum, and hence we have there z −k −j = O(δ e(l) ) by (41) and d k j = O(δ e(l)+(2K−e(l))(1−2ν) ) by Lemma 11 (since m(j, k) ≥ e(l)). This shows that the variation in the almost-invariant energy E l is indeed of the claimed order, and it shows the additional statement for l = 1.