Spontaneous periodic orbits in the Navier-Stokes flow

TL;DR

This paper introduces a computer-assisted method to prove the existence of spontaneous periodic orbits in the three-dimensional Navier-Stokes equations with specific forcing, utilizing a zero-finding approach on Fourier coefficients.

Contribution

It develops a general constructive framework using a Newton-Kantorovich theorem for proving periodic orbits in Navier-Stokes flows, incorporating symmetry reductions and analytic estimates.

Findings

Proofs of spontaneous periodic orbits in Navier-Stokes with Taylor-Green forcing

Application of a zero-finding problem on Fourier coefficients

Use of symmetry to reduce computational complexity

Abstract

In this paper, a general method to obtain constructive proofs of existence of periodic orbits in the forced autonomous Navier-Stokes equations on the three-torus is proposed. After introducing a zero finding problem posed on a Banach space of geometrically decaying Fourier coefficients, a Newton-Kantorovich theorem is applied to obtain the (computer-assisted) proofs of existence. The required analytic estimates to verify the contractibility of the operator are presented in full generality and symmetries from the model are used to reduce the size of the problem to be solved. As applications, we present proofs of existence of spontaneous periodic orbits in the Navier-Stokes equations with Taylor-Green forcing.

Figures

Tables

Table 1. Table 1: Parameters for the two rigorously computed solutions of the Navier-Stokes equations with Taylor-Green forcing ( 1.2 ). The number N x 1 subscript 𝑁 subscript 𝑥 1 N_{x_{1}} , N x 2 subscript 𝑁 subscript 𝑥 2 N_{x_{2}} , N x 3 subscript 𝑁 subscript 𝑥 3 N_{x_{3}} and N t subscript 𝑁 𝑡 N_{t} define the set 𝒮 sol superscript 𝒮 sol \mathcal{S}^{\text{{sol}}} in ( 5.7 ). The solutions are indicated by labels p 1 subscript p 1 \mathrm{p}_{1} and p 2 subscript p 2 \mathrm{p}_{2} in bifurcation diagram 3 and are computed for ν = 0.286 𝜈 0.286 \nu=0.286 and ν = 0.265 𝜈 0.265 \nu=0.265 respectively. The computations were timed on an Intel Xeon E5-1620v2 with a 3.7GHz clock speed.

	$\eta$	$N_{x_1}$	$N_{x_2}$	$N_{x_3}$	$N_t$	$N^{†}$	$\tilde{N}$	RAM (GB)	CPU days
$p_1$	$1$	17	17	0	11	130	265	10	6
$p_2$	$1$	21	21	0	16	210	425	110	95

Equations

$$\{\begin{array}{rl}{\partial }_{t}u+(u\cdot \nabla )u-\nu \Delta u+\nabla p& =f\\ \nabla \cdot u& =0,\end{array}$$

$$f=f(x)=\left(\begin{array}{c}\text{hfill}2\sin x_1\cos x_2\\ -2\cos x_1\sin x_2\\ 0\end{array}\right),$$

$$\Vert u-{\bar{u}}_{\Omega /\bar{\Omega }}{\Vert }_{{\mathcal{C}}^0}\le r_{\text{sol}}^u\text{and}\Vert p-{\bar{p}}_{\Omega /\bar{\Omega }}{\Vert }_{{\mathcal{C}}^0}\le r_{\text{sol}}^p.$$

$$u(x,t)=\sum _{n\in {\mathbb{Z}}^4}{u}_n{e}^{i(\tilde{n}\cdot x+n_4\Omega t)},$$

$${\left(D_{l}a\right)}_n\stackrel{\text{mbox}[}{=}{n}_l{a}_n\text{for all }n\in {\mathbb{Z}}^4.$$

$$(a\ast b{)}_n\stackrel{\text{mbox}[}{=}\sum _{k\in {\mathbb{Z}}^4}{a}_k{b}_{n-k}.$$

$${\left[\left(a\star \tilde{D}\right)b\right]}^{(l)}\stackrel{\text{mbox}[}{=}\sum _{m=1}^3{a}^{(m)}\ast \left(D_m{b}^{(l)}\right)\text{for all }1\le l\le 3,$$

$${\left[\left(a\star \tilde{D}\right)b\right]}^{(l)}=\sum _{m=1}^3{D}_m\left[a^{(m)}\ast b^{(l)}\right].$$

$${\left[\left(a\star \tilde{D}\right)b\right]}^{(l)}$$

$$=\sum _{m=1}^3\left(D_m\left[a^{(m)}\ast b^{(l)}\right]-\left[D_m{a}^{(m)}\ast b^{(l)}\right]\right)$$

$$=\sum _{m=1}^3{D}_m\left[a^{(m)}\ast b^{(l)}\right].\text{qed}$$

$$(u\cdot \nabla )u=\nabla \left(\frac{|u{|}^2}{2}\right)-u\times \omega ,$$

$$\nabla \times \left((u\cdot \nabla )u\right)$$

$$=\left(u\cdot \nabla \right)\omega -\left(\omega \cdot \nabla \right)u+\omega \left(\nabla \cdot u\right)-u\left(\nabla \cdot \omega \right),$$

$$\nabla \times \left((u\cdot \nabla )u\right)=\left(u\cdot \nabla \right)\omega -\left(\omega \cdot \nabla \right)u.$$

$${\partial }_t\omega +\left(u\cdot \nabla \right)\omega -\left(\omega \cdot \nabla \right)u-\nu \Delta \omega =f^{\omega }\text{on }{\mathbb{T}}^3\times \mathbb{R},$$

$$\{\begin{array}{rl}& \nabla \times u=\omega \\ & \nabla \cdot u=0.\end{array}$$

$$-\Delta u=\nabla \times \omega ,$$

$$u=(-\Delta {)}^{-1}\nabla \times \omega .$$

$${\int }_{{\mathbb{T}}^3}u(x,t)dx=\sum _{n_4\in \mathbb{Z}}{u}_{0,n_4}{e}^{i{n}_4\Omega t}$$

$$\frac{d}{dt}{\int }_{{\mathbb{T}}^3}u(x,t)dx={\int }_{{\mathbb{T}}^3}f(x)dx.$$

$$M_n\stackrel{\text{mbox}[}{=}\{\begin{array}{ll}\frac{i}{{\tilde{n}}^2}\left(\begin{array}{ccc}0& -n_3& n_2\\ n_3& 0& -n_1\\ -n_2& n_1& 0\end{array}\right),& \tilde{n}\ne 0,\\ 0,& \tilde{n}=0,\end{array}$$

$$M\omega \stackrel{\text{mbox}[}{=}{\left(M_n{\omega }_n\right)}_{n\in {\mathbb{Z}}^4},$$

$$\left\{\begin{array}[]{l}\nabla\times u=\omega\\ \nabla\cdot u=0\\ \int_{\mathbb{T}^{3}}u=0\end{array}\right.\qquad\Longleftrightarrow\qquad\left\{\begin{array}[]{l}u=M\omega\\ \nabla\cdot\omega=0.\end{array}\right.$$

$${\partial }_t\omega +\left(M\omega \cdot \nabla \right)\omega -\left(\omega \cdot \nabla \right)M\omega -\nu \Delta \omega =f^{\omega }\text{on }{\mathbb{T}}^3\times \mathbb{R},$$

$$W=\left(\begin{array}{c}\Omega \\ {\left({\omega }_n\right)}_{n\in {\mathbb{Z}}^4\setminus \{0\}}\end{array}\right),$$

$$\omega (x,t)=\sum _{n\in {\mathbb{Z}}_{\ast }^4}{\omega }_n{e}^{i(\tilde{n}\cdot x+n_4\Omega t)}=\sum _{n\in {\mathbb{Z}}^4}{\omega }_n{e}^{i(\tilde{n}\cdot x+n_4\Omega t)}.$$

$$F_n(W)\stackrel{\text{mbox}[}{=}i\Omega n_4{\omega }_n+i{\left[\left(M\omega \star \tilde{D}\right)\omega \right]}_n-i{\left[\left(\omega \star \tilde{D}\right)M\omega \right]}_n+\nu {\tilde{n}}^2{\omega }_n-f_n^{\omega },$$

$${\Psi }_n(\omega )\stackrel{\text{mbox}[}{=}i{\left[\left(M\omega \star \tilde{D}\right)\omega \right]}_n-i{\left[\left(\omega \star \tilde{D}\right)M\omega \right]}_n\text{for all }n\in {\mathbb{Z}}_{\ast }^4.$$

$$D_1{\left(M\omega \right)}^{(1)}+D_2{\left(M\omega \right)}^{(2)}+D_3{\left(M\omega \right)}^{(3)}=0,$$

$${\left[\left(M\omega \star \tilde{D}\right)\omega \right]}_{0,n_4}=0={\left[\left(\omega \star \tilde{D}\right)M\omega \right]}_{0,n_4}\text{for all }{n}_4\in \mathbb{Z}.$$

$$\Phi \stackrel{\text{mbox}[}{=}{\partial }_{t}u+(u\cdot \nabla )u-\nu \Delta u-f.$$

Full text

Keywords

Navier-Stokes equations $\cdot$ periodic orbits $\cdot$ symmetry breaking $\cdot$ computer-assisted proofs

Mathematics Subject Classification (2010)

35Q30 $\cdot$ 35B06 $\cdot$ 35B10 $\cdot$ 35B36 $\cdot$ 65G20 $\cdot$ 76D17

1 Introduction

The Navier-Stokes equations for a fluid of constant density $\rho$ can be expressed as

[TABLE]

where $u=u(x,t)$ is the velocity, $p(x,t)=P(x,t)/\rho$ is the pressure scaled by the density, $\nu$ is the kinematic viscosity and $f=f(x,t)$ is an external forcing term. These equations can be considered on compact or unbounded domains, complemented by boundary and initial conditions. The first equation, which expresses momentum balance, has a quadratically nonlinear advection term. While the presence of the nonlinearity generally obstructs obtaining closed form solutions, there are some notable exceptions. In parallel shear flows, the advection term vanishes identically and analytic solutions are available. Examples include Hagen-Poiseuille flow in pipes [34] and Taylor-Couette flow between co-axial cylinders [36]. In Beltrami flows, the nonlinearity takes the form of a gradient and can be absorbed in the pressure. An example of an explicit solution with this property is the ABC flow [7]. Explicit solutions with non-trivial nonlinearities are rare, but some are known. For instance, exact vortical solutions include Burgers’ vortex on $\mathbb{R}^{3}$ [3] and the periodic vortex array of Taylor and Green [37]. However, it has been known for a long time that, even in fluids with strong viscous damping, more complicated, time-periodic motions can occur. A classical experiment is that of a fluid flowing past a stationary cylinder. Experiments started by von Kármán at the beginning of the twentieth century, and carried on by his students, showed that, at a well-defined flow rate, the motion in the wake of the cylinder becomes time-periodic, as alternating clockwise and counter clockwise vortices travel downstream [19].

While there is little hope of writing down explicit solutions that describe such oscillatory behaviour, a number of authors have attempted to at least establish the existence of time-periodic solutions. The earliest contribution was likely the work of James Serrin. In 1959, he published two papers on the existence and stability of certain solutions to the Navier-Stokes equations in the limit of large viscosity. In the first one, he established the existence of globally stable equilibrium solutions by finding bounds for the nonlinear and forcing terms, and by showing that a certain energy decays [32]. In the second one, he considered large viscosity and gave a criterion for the existence of periodic solutions on a three-dimensional bounded domain subject to time-periodic boundary data and body forces [31].

Many authors followed Serrin in studying the periodically forced non-autonomous Navier-Stokes system dominated by viscosity. Kaniel and Shinbrot [15] considered bounded domains with fixed boundaries and showed the existence of periodic strong solutions for small time-periodic forcing $f$. Without making any assumption about the size of $f$, Takeshita [35] showed the same result as Kaniel and Shinbrot. Some time later, Teramoto [38] proved the existence of time-periodic solutions for domains with slowly moving boundaries. Then, Maremonti [23] and Kozono and Nakao [20] extended the results from bounded domains to $\mathbb{R}^{3}$. The latter made use of the $L^{p}$ theory of the Stokes operator rather than the energy method. A similar result, relying on a milder condition on the forcing function, was derived by Kato [16]. Other extensions were those to inhomogeneous boundary conditions on compact domains by Farwig and Okabe [8] and to the case of a rotating fluid in two dimensions by Hsia et al. [14]. The latter paper also contains a fairly extensive list of references of which only a fraction is discussed here.

Thus, our understanding of periodic flows in response to time-periodic forcing is rather advanced. The same cannot be said about spontaneous periodic motions, which we refer to as being periodic flows driven by a time-independent forcing. In other words, spontaneous periodic motions are periodic orbits in the autonomous Navier-Stokes equations. The regular vortex shedding in the wake of a cylinder, for instance, arises in the absence of a body force and as a consequence of the nonlinearity of the Navier-Stokes equation, not by virtue of the advection being dominated by viscous damping. In an attempt to address the difficulties in studying spontaneous motions, the present paper proposes a general (computer-assisted) approach to prove existence of time-periodic Navier-Stokes flows on the three-torus for given time-independent forcing terms $f=f(x)$.

The novelty of our paper is threefold. Foremost, it provides the first computer-assisted proof of existence of spontaneous periodic Navier-Stokes flows. Second, it introduces general analytic bounds applicable to prove existence of three-dimensional time-periodic solutions for any time-independent forcing term. Third, it uses the symmetries present in Navier-Stokes to significantly reduce the size of the problem to work with.

A few comments on the symmetries are in order. Our approach permits one to take advantage of symmetries of the forcing $f$, and in particular of those (subgroup of) symmetries that are also obeyed by the examined solution. We allow general time-independent forcings with zero spatial average (so that periodic solutions are not ruled out a priori, see Equation (2.5)). Any time-periodic solution thus spontaneously breaks the shift symmetry in time, and other symmetries of $f$ may also be broken by the solution. However, the bigger the symmetry group of the solution is, the more we can reduce the computational cost (in terms of time and, especially, memory).

While all of the analysis is performed in full generality on the 3-torus, the solutions we present in Theorem 1.1 below are two-dimensional (in space) time-periodic solutions. Indeed, they are homogeneous in one spatial variable and can thus be interpreted as solutions on the 2-torus. The only reason for this reduction is that the physical memory requirements for a three-dimensional solution are, for now, prohibitive in our current implementation. To be precise, we consider the so-called Taylor-Green forcing

[TABLE]

which corresponds to counter rotating vortex columns. Clearly, this forcing allows one to restrict to the first two spatial variables. It is expected that some periodic solutions in fact break the 2D symmetry of the forcing (1.2), see also Section 5. While we aim to investigate such solutions in future work, the solutions obtained in the current paper respect the 2D symmetry: they are independent of $x_{3}$ and the third component of the velocity vanishes. We call such a solution an (essentially) 2D solution. In addition, the solutions that we find here are invariant under a symmetry group with 16 elements, see Section 5 for details. This allows us to reduce the number of independent Fourier modes on which we perform the computational analysis by a factor 16, which represents considerable savings in memory requirements. Finally, due to the shift-invariance of the torus, in general it may be appropriate to look for solutions which are shift-periodic, but such complications do not arise when studying 2D solutions for the forcing (1.2).

Before we state a representative sample result, we note that determining the period of the solution is part of the problem. Hence the frequency $\bar{\Omega}$ of a numerical approximate solution $(\bar{u},\bar{p})$ only approximates the true frequency $\Omega$. The solution of (1.1) will therefore be close to a slightly time-dilated version $(\bar{u}_{\theta},\bar{p}_{\theta})$ of the numerical data, where $\theta$ is the dilation factor, see Remark 2.17. As outlined below in more detail, we use a Newton operator to show that, under computable conditions, there is a solution to (1.1) near the numerically obtained approximation $(\bar{u},\bar{p})$, where the error is bounded explicitly. As an example, we prove the following result.

Theorem 1.1**.**

Consider (1.1) defined on the three-torus $\mathbb{T}^{3}$ (with size length $L=2\pi$) and consider the time-independent forcing term (1.2). Let $\nu=0.265$ and $(\bar{u},\bar{p})$ be the numerical solution whose Fourier coefficients and time frequency $\bar{\Omega}$ are given in the file dataorbit2.mat and can be downloaded at [41] (and whose vorticity is represented in Figure 1). Let $r_{\text{{sol}}}^{\Omega}=2.2491\cdot 10^{-6}$, $r_{\text{{sol}}}^{u}=2.2491\cdot 10^{-6}$, and $r_{\text{{sol}}}^{p}=5.6486\cdot 10^{-5}$. There exists a $\frac{2\pi}{\Omega}$-periodic solution $(u,p)$ of (1.1) with $|\Omega-\bar{\Omega}|\leq r_{\text{{sol}}}^{\Omega}$ and such that

[TABLE]

We point out that the $\mathcal{C}^{0}$-norm is only used here to get a simple statement. A more general version of Theorem 1.1, with a stronger norm which is the one actually used in the analysis, is presented in Theorem 5.2 in Section 5.

It is important to recognize that in the last forty years, important open problems were settled with computer-assisted proofs: the universality of the Feigenbaum constant [22], the four-colour theorem [27], the existence of the strange attractor in the Lorenz system [39] (i.e. Smale’s 14th problem) and Kepler’s densest sphere packing problem [12]. We refer the interested to reader to the expository works [42, 11, 18, 24, 25, 26, 29, 40] and the references therein, for a more complete overview of the field of rigorously verified numerics. Let us however mention some results related to the present work. In [46], Watanabe proposes an approach to obtain computer-assisted proofs of existence of stationary solution in the Navier-Stokes equation, which then lead to the rigorous computation of equilibria in a three-dimensional thermal convection problem [17] and in Kolmogorov flow, i.e. flow with periodic boundary conditions and a constant body force with a simple structure [44, 45]. Independently, Heywood et al. [13] established fixed-point theorems for steady, two-dimensional Kolmogorov flows. Their results fall short of a proof of existence only because of the presence of round-off error, which Watanabe avoided by using interval arithmetic. The rigorous computation of time-dependent solutions to the autonomous Navier-Stokes equation has so far been out of reach. We note that computer-assisted proofs for periodic orbits, along lines similar to the current paper, have been obtained for the Kuramoto-Shivashinsky PDE [1, 9, 10, 48] and the ill-posed Boussinesq equation [5].

Our strategy begins by identifying a problem of the form $\mathcal{F}(W)=0$ posed on a Banach algebra of geometrically decaying Fourier coefficients, whose solutions yield the time-periodic orbits. This zero finding problem is derived by applying the curl operator to (1.1) and solving for the periodic orbits in the vorticity equation. Expressing a periodic orbit using a space-time Fourier series and plugging the series in the vorticity equation yields the infinite dimensional nonlinear problem $\mathcal{F}(W)=0$, where $W$ corresponds to the sequence of Fourier coefficients of the vorticity $\omega=\nabla\times u$. The details of the derivation of the map $\mathcal{F}$ are given in Section 2.1. A proof that the solutions of $\mathcal{F}=0$ correspond to time-periodic Navier-Stokes flows is presented in Lemma 2.5. The next step is to consider a finite dimensional projection of $\mathcal{F}$ and to numerically obtain an approximation $\bar{W}$ of a zero of $\mathcal{F}$, that is $\mathcal{F}(\bar{W})\approx 0$. Next, we turn the problem $\mathcal{F}(W)=0$ into an equivalent fixed point problem of the form $T(W)=W-D\mathcal{F}(\bar{W})^{-1}\mathcal{F}(W)$. We then set out to prove that $T$ is a contraction on a neighborhood of $\bar{W}$. The advantage is that instead of trying to prove equalities in the formulation $\mathcal{F}(W)=0$, contractivity involves inequalities only. The proof then proceeds by a Newton-Kantorovich type argument (see Theorem 2.15 and Theorem 4.23) to find a ball centered at $\bar{W}$ on which the map $T$ is a contraction mapping. Having done the hard work in the analysis of reducing the problem to finitely many explicit inequalities, one therefore resorts to interval arithmetic computer calculations for this final step of the proof.

The paper is organized as follows. In Section 2, we introduce the rigorous computational approach and the zero finding problem $\mathcal{F}(W)=0$, as well as the Banach space in which we solve for the zeros of $\mathcal{F}$. The Newton-Kantorovich type theorem is presented in Theorem 2.15. In Section 3, we introduce the general bounds necessary to verify the hypotheses of Theorem 2.15. Then in Section 4 we describe how the symmetries of the model can be used to simplify solving the zero-finding problem, by reducing significantly its size. Using this reduction based on symmetries a modified Newton-Kantorovich theorem is proved (Theorem 4.23) and the associated symmetry-adapted estimates are derived in Section 4.6. Sample results are then presented in Section 5. All the estimates obtained in this paper culminate in Theorems 5.1 and 5.2, which allows us to validate periodic solutions $\omega$ of the vorticity equation, with explicit error bounds. In the Appendix we describe how to recover errors bounds for the associated velocity $u$ and pressure $p$ that solve the Navier-Stokes equations.

2 The rigorous computational approach

This section is devoted to the presentation of the framework that is needed to study periodic solutions of (1.1) by computer-assisted means. We first derive a suitable $\mathcal{F}=0$ problem in Section 2.1 and introduce the proper Banach spaces to study that problem in Section 2.2. Well chosen approximations of $D\mathcal{F}$ and $D\mathcal{F}^{-1}$ are then introduced in Section 2.3, and used in Section 2.4 to state Theorem 2.15, which provides us with sufficient conditions for the existence of non trivial zeros of $\mathcal{F}$.

2.1 The zero finding problem in Fourier space

In this section, we introduce the zero finding problem $\mathcal{F}=0$ that we are going to work with. We start by some (somewhat algebraic) manipulations and then explain in Lemma 2.5 how the zero finding problem is related to the original Navier-Stokes equations.

We consider the 3D incompressible Navier-Stokes equations (1.1) on the three-torus $\mathbb{T}^{3}$ and look for time-periodic solutions. As mentioned in the introduction, both the numerical and the theoretical part of our work are based on Fourier series, for which we will use the following notations. For $n=(n_{1},n_{2},n_{3},n_{4})\in\mathbb{Z}^{4}$, we write $n=(\tilde{n},n_{4})$ where $\tilde{n}=(n_{1},n_{2},n_{3})$ and $\tilde{n}^{2}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}n_{1}^{2}+n_{2}^{2}+n_{3}^{2}$. If $u:\mathbb{T}^{3}\times\mathbb{R}\to\mathbb{R}^{3}$ is a time-periodic function (that is periodic in the fourth variable), we denote by $\left(u_{n}\right)_{n\in\mathbb{Z}^{4}}\in\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}}$ its Fourier coefficients:

[TABLE]

where $\Omega$ is the a-priori unknown angular frequency.

Remark 2.1**.**

In this paper, we are only concerned with smooth (that is analytic) periodic functions. Therefore, we can identify a function with its sequence of Fourier coefficients, and to make the notations lighter we use the same symbol to denote both of them. It should be clear from context whether $u$ (and similarly later for $\omega$, $f$, $f^{\omega}$, etc.) denotes a periodic function or a sequence of Fourier coefficients.

For $1\leq l\leq 3$, we use $u^{(l)}\in\mathbb{C}^{\mathbb{Z}^{4}}$ to denote the Fourier sequence of the $l$-th component of $u$. For any Fourier sequence $a=(a_{n})\in\mathbb{C}^{\mathbb{Z}^{4}}$ and any $1\leq l\leq 3$ we define the sequence $D_{l}a$ corresponding to the partial derivative of $a$ with respect to $x_{l}$ (up to a factor $i$):

[TABLE]

For any Fourier sequence $a=(a_{n})$ and $b=(b_{n})$ in $\mathbb{C}^{\mathbb{Z}^{4}}$ we define their convolution product as

[TABLE]

Finally, given $a,b\in\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}}$ we define

[TABLE]

which is the $l$-th component in Fourier space of $(a\cdot\nabla)b$, again up to a factor $i$. We will frequently use the following lemma.

Lemma 2.2**.**

Let $a,b\in\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}}$ be such that $\sum_{m=1}^{3}D_{m}a^{(m)}=0$. Then, for all $l\in\{1,2,3\}$,

[TABLE]

Proof.

This is a consequence of the product rule:

[TABLE]

We are now almost ready to set up our $\mathcal{F}=0$ problem, but instead of looking directly at periodic solutions of (1.1) we are going to work with the vorticity equation. Namely, we consider the vorticity $\omega=\nabla\times u$ and look for the equation it satisfies. Using

[TABLE]

we get

[TABLE]

and since both $u$ and $\omega$ are divergence free we end up with

[TABLE]

The vorticity equation is then given by

[TABLE]

where $f^{\omega}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\nabla\times f$.

We are going to solve for the Fourier coefficients of $\omega$ satisfying the vorticity equation (2.3). More precisely, our unknowns are the Fourier coefficients $\left(\omega_{n}\right)_{n\in\mathbb{Z}^{4}}$ and the angular frequency $\Omega$. However, in (2.3) the unknowns $\omega$ and $u$ both appear. To obtain an equation depending on the vorticity $\omega$ only, we need to express $u$ in term of $\omega$ by solving

[TABLE]

Applying the curl operator to the first equation, and using that $\nabla\cdot u=0$, we get

[TABLE]

and so, formally,

[TABLE]

Remark 2.3**.**

Expression (2.4) is not completely well defined because the Laplacian has, in general, a non-zero kernel. In particular, the space average velocity

[TABLE]

cannot be recovered from the vorticity $\omega$. However, going back to (1.1) we have

[TABLE]

In this work, we consider a time independent forcing with spatial average equal to zero, therefore the space average velocity is a conserved quantity, which can always be assumed to be zero by Galilean invariance.

To give a well defined version of (2.4), we go through Fourier space and introduce

[TABLE]

and

[TABLE]

As mentioned previously, we can assume the space average velocity to be zero, which is why we define $(M\omega)_{0,n_{4}}$ to be zero for all $n_{4}\in\mathbb{Z}$.

Remark 2.4**.**

The above computations and construction of $M$ can be summarized in the following equivalence

[TABLE]

Going back to (2.3) and replacing $u$ by $M\omega$, we obtain the equation

[TABLE]

which is the one we are going to work with. More precisely, we first define

[TABLE]

which corresponds to all the unknowns we are solving for. Notice that $\omega_{0}$ is not part of the unknowns, as it can always be taken equal to [math]. In the sequel, to simplify the presentation we introduce $\mathbb{Z}^{4}_{*}=\mathbb{Z}^{4}\setminus\{0\}$, and always identify $\left(\omega_{n}\right)_{n\in\mathbb{Z}^{4}_{*}}$ with $\left(\omega_{n}\right)_{n\in\mathbb{Z}^{4}}$ where $\omega_{0}=0$. In particular, the representation in physical space associated to $W$ is given by

[TABLE]

We then define $F=\left(F_{n}\right)_{n\in\mathbb{Z}^{4}_{*}}$ by

[TABLE]

for all $n\in\mathbb{Z}^{4}_{*}$, and aim to show the existence of a nontrivial zero of $F$.

For later use (see Section 3.4 and Section 4.6.3), we also introduce the notation $\Psi=\left(\Psi_{n}\right)_{n\in\mathbb{Z}^{4}_{*}}$ for the nonlinear terms in (2.10):

[TABLE]

Lemma 2.5**.**

Let $W\in\mathbb{R}\times(\mathbb{C}^{3})^{\mathbb{Z}^{4}_{*}}$ be such that the corresponding function $\omega$ is analytic. Assume that $F(W)=0$ and $\nabla\cdot\omega=0$. Assume also that $f$ does not depend on time and has space average zero. Define $u=M\omega$. Then there exists a pressure function $p:\mathbb{T}^{3}\times\mathbb{R}\to\mathbb{R}$ such that $(u,p)$ is a $\frac{2\pi}{\Omega}$-periodic solution of (1.1).

Remark 2.6**.**

Let us mention that the analyticity condition could be considerably weakened: the lemma would hold for any smoothness assumption allowing to justify all the taken derivatives and the switches between functional and Fourier representation. In our case analyticity simply happens to be the most natural assumption, because of the space of Fourier coefficients we end up using, see Section 2.2 and Remark 2.8.

Proof.

First, a straightforward computation gives that

[TABLE]

for any $\omega$, which amounts to saying that $\nabla\cdot M\omega=0$, therefore $\nabla\cdot u=0$.

The next step is to prove that $\nabla\times u=\omega$ (that is $\nabla\times M\omega=\omega$). Using the definition of $M_{n}$ and the fact that $\nabla\cdot\omega=0$, another straightforward computation in Fourier space shows that $\left(\nabla\times M\omega\right)_{n}=\omega_{n}$ for all $\tilde{n}\neq 0$. To prove that $\left(\nabla\times M\omega\right)_{n}=\omega_{n}$ for all $\tilde{n}=0$, first observe that $\left(\nabla\times M\omega\right)_{n}=0$ for $\tilde{n}=0$. Next, since $\omega$ and $M\omega$ are divergence free, we can use Lemma 2.2 on both nonlinear terms of $F$, and get that

[TABLE]

Therefore, $F_{0,n_{4}}(W)=0$ implies that $\omega_{0,n_{4}}=0$ for $n_{4}\neq 0$, which concludes the proof that $\nabla\times u=\omega$.

Finally, we define the periodic function

[TABLE]

Using $\nabla\times u=\omega$ we get that, for all $n\neq 0$, $\left(\nabla\times\Phi\right)_{n}=F_{n}(W)$, and since $(\nabla\times\Phi)_{0}$ vanishes as well, we conclude that $\nabla\times\Phi=0$. Recalling that $u=M\omega$ is divergence free, Lemma 2.2 yields that

[TABLE]

Also using that $\left(M\omega\right)_{0,n_{4}}=0$ and that $f_{0,n_{4}}=0$ (since $f$ does not depend on time and has average zero), we get that $\Phi_{0,n_{4}}=0$ for all $n_{4}\in\mathbb{Z}$. We have shown that

[TABLE]

Therefore (see Lemma 6.1) there exists a $p$ such that $\Phi=-\nabla p$, that is

[TABLE]

Since we had already shown that $u$ is divergence free, this completes the proof. ∎

Finally, in view of the contraction argument that we are going to apply, we add a phase condition in order to isolate the solution. Assume that we have an approximate periodic orbit given by $\hat{\omega}$. A common choice for the phase condition (see e.g. [21]) is to require that the orbit $\omega$ satisfies

[TABLE]

Hence, we define the phase condition

[TABLE]

where the superscript ∗ denotes complex conjugation, and we assumed that the reference orbit given by $\hat{\omega}$ is real-valued, that is $\hat{\omega}_{-n}^{(l)}=\left(\hat{\omega}_{n}^{(l)}\right)^{*}$.

We now consider the enlarged problem

[TABLE]

In view of the similarity between (2.8) and (2.14), we will abuse notation and refer to the set of elements of such variables via

[TABLE]

where $Q_{\leftmoon}\in\mathbb{C}$ and $Q_{n}\in\mathbb{C}^{3}$ for $n\in\mathbb{Z}^{4}_{*}$, rather than introducing notation for the projections onto different components.

2.2 Banach spaces, norms, index sets

In this section we introduce the Banach spaces on which we are going to study $\mathcal{F}$, as well as some additional notations that are going to be used throughout this paper.

For $\eta\geq 1$, we denote by $\ell^{1}_{\eta}(\mathbb{C})$ the subspace of all sequences $a\in\mathbb{C}^{\mathbb{Z}^{4}_{*}}$ such that

[TABLE]

with $\left|n\right|_{1}=\sum_{j=1}^{4}|n_{j}|$. We also introduce the subspaces $\ell^{1}_{\eta,-2,-1}(\mathbb{C})$, $\ell^{1}_{\eta,-1,-1}(\mathbb{C})$ and $\ell^{1}_{\eta,-1,0}(\mathbb{C})$ of $\mathbb{C}^{\mathbb{Z}^{4}_{*}}$, associated to the norms

[TABLE]

and

[TABLE]

with $|n|_{\infty}=\max_{1\leq j\leq 4}|n_{j}|$ and $|\tilde{n}|_{\infty}=\max_{1\leq j\leq 3}|n_{j}|$.

The main space in which we are going to work is the Banach space $\mathcal{X}=\mathbb{C}\times\left(\ell^{1}_{\eta}(\mathbb{C})\right)^{3}$ with the norm

[TABLE]

Remark 2.7**.**

Since $\Omega$ and $\omega$ are incommensurable, $\Omega$ and $\omega$ do not live in the same space it is prudent to introduce an extra weight in the norm for the $|\Omega|$ term. Depending on the application at hand this may indeed be necessary, but for the results in the current paper setting this weight to unity suffices.

Remark 2.8**.**

Notice that, as soon as $\Omega\in\mathbb{R}$ and $\eta>1$, the function $\omega$ associated to an element $W$ of $\mathcal{X}$ via (2.9) is analytic.

Similarly, we introduce the Banach spaces $\mathcal{X}_{-2,-1}=\mathbb{C}\times\left(\ell^{1}_{\eta,-2,-1}(\mathbb{C})\right)^{3}$ and $\mathcal{X}_{-1,-1}=\mathbb{C}\times\left(\ell^{1}_{\eta,-1,-1}(\mathbb{C})\right)^{3}$, respectively endowed with the norms

[TABLE]

and

[TABLE]

Notice that $\mathcal{F}$ defined in (2.14) maps $\mathcal{X}$ into $\mathcal{X}_{-2,-1}$. We are also going to consider subspaces of divergence free sequences of $\mathcal{X}$ and $\mathcal{X}_{-2,-1}$, namely:

[TABLE]

Notice that $\omega\mapsto\sum_{m=1}^{3}D_{m}\omega^{(m)}$ is a bounded linear map from $\left(\ell^{1}_{\eta}\right)^{3}$ to $\ell^{1}_{\eta,-1,0}$, therefore $\mathcal{X}^{\text{{div}}}$ is a closed subspace of $\mathcal{X}$ and thus still a Banach space, with the norm $\left\|\cdot\right\|_{\mathcal{X}}$. Similarly, $\mathcal{X}_{-2,-1}^{\text{{div}}}$ is a Banach space, with the norm $\left\|\cdot\right\|_{\mathcal{X}_{-2,-1}}$.

We recall that, in Lemma 2.5, we need the zero of $F$ to be divergence free in order to prove that it corresponds to a solution of Navier-Stokes equation (1.1). Therefore, when we later prove the existence of a zero of $\mathcal{F}$ in $\mathcal{X}$, it is crucial to be able to show that this zero is actually in $\mathcal{X}^{\text{{div}}}$. To do so in Theorem 2.15, we will make use of the following observation.

Lemma 2.9**.**

$\mathcal{F}$* maps $\mathcal{X}^{\text{{div}}}$ to $\mathcal{X}_{-2,-1}^{\text{{div}}}$.*

Proof.

We need to show that, for all $\omega$ satisfying $\nabla\cdot\omega=0$, we have $\nabla\cdot F(W)=0$. The only term in $F(W)$ that is not obviously divergence free is the nonlinear term, but since $\nabla\cdot\omega=0$ by assumption and $\nabla\cdot M\omega=0$ by construction of $M$, we can proceed as in (2.1)-(2.2) to rewrite the nonlinear term as a curl, therefore concluding that it is indeed divergence free. ∎

Similarly, we want to obtain zero of $\mathcal{F}$ that corresponds to a real-valued function. To ensure this, the following notation and observation are needed.

Definition 2.10**.**

We define the complex conjugation symmetry $\gamma_{*}$, acting on $\mathbb{C}\times\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}_{*}}$, by

[TABLE]

where ∗ denotes complex conjugation, which is to be understood component-wise when applied to $Q_{n}\in\mathbb{C}^{3}$. We still use the symbol $\gamma_{*}$ to denote the complex conjugation symmetry acting on $\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}_{*}}$.

Notice that, thanks to the factor $i$ in the definition (2.13) of $F_{\leftmoon}$, $\mathcal{F}$ is $\gamma_{*}$-equivariant. More precisely, one has the following statement, which will be used for Theorem 2.15.

Lemma 2.11**.**

Assume that $\hat{\omega}$ used in the phase condition (2.13) is such that $\gamma_{*}\hat{\omega}=\hat{\omega}$. Then

[TABLE]

We end this section with a last set of notations that will be useful to describe linear operators.

Notation 2.12**.**

To work with a linear operator $B:\mathcal{X}\to\mathcal{X}$, it is convenient to introduce the block notation

[TABLE]

where

• •

$B^{(l,m)}$* is a linear operator from $\ell^{1}_{\eta}$ to $\ell^{1}_{\eta}$, for all $1\leq l,m\leq 3$,*

• •

$B^{(\leftmoon,m)}$* is a linear operator from $\ell^{1}_{\eta}$ to $\mathbb{C}$, for all $1\leq m\leq 3$,*

• •

$B^{(l,\leftmoon)}$* is a linear operator from $\mathbb{C}$ to $\ell^{1}_{\eta}$, for all $1\leq l\leq 3$,*

• •

$B^{(\leftmoon,\leftmoon)}$* is a linear operator from $\mathbb{C}$ to $\mathbb{C}$.*

For all $1\leq l,m\leq 3$, we write $B^{(l,m)}=\left(B^{(l,m)}_{k,n}\right)_{k,n\in\mathbb{Z}^{4}_{*}}$, so that, for all $a\in\ell^{1}_{\eta}$ and all $k\in\mathbb{Z}^{4}_{*}$,

[TABLE]

Similarly, $B^{(\leftmoon,m)}=\left(B^{(\leftmoon,m)}_{n}\right)_{n\in\mathbb{Z}^{4}_{*}}$ and $B^{(l,\leftmoon)}=\left(B^{(l,\leftmoon)}_{k}\right)_{k\in\mathbb{Z}^{4}_{*}}$, so that, for all $a\in\ell^{1}_{\eta}$, all $\Omega\in\mathbb{C}$ and all $k\in\mathbb{Z}^{4}_{*}$,

[TABLE]

We also use

[TABLE]

to denote all the columns of $B$. Similar notations will be used when $B$ is a linear operator from $\mathcal{X}$ to $\mathcal{X}_{-2,-1}$, or vice versa.

2.3 The operators $\widehat{A}$ and $A$

We now assume that we have an approximate zero $\bar{W}=(\bar{\Omega},\bar{\omega})\in\mathcal{X}$ of $\mathcal{F}$. In practice $\bar{W}$ will only a have finite number of non zero coefficients (see (3.1)), but this is not crucial for the moment.

As mentioned in Section 1, our strategy to obtain the existence of a zero of $\mathcal{F}$ in a neighborhood of $\bar{W}$ is to show that a Newton-like operator of the form

[TABLE]

is a contraction in a neighborhood of $\bar{W}$. To prove this directly, we would need to derive a (computable) estimate of $\left\|D\mathcal{F}(\bar{W})^{-1}\right\|_{B(\mathcal{X}_{-2,-1},\mathcal{X})}$ (see Remark 2.16), which would be a quite formidable task. We circumvent this difficulty by working with an approximate inverse $A$ of $D\mathcal{F}(\bar{W})^{-1}$, which is defined via a finite dimensional numerical approximation, and for which estimates are much easier to obtain. More precisely, we are going to consider a finite dimensional projection of $D\mathcal{F}(\bar{W})$, invert it numerically and then use the numerical inverse to construct $A$. To define $A$ precisely, it will be convenient to use the following notation.

Definition 2.13**.**

Let $\nu$ be the kinematic viscosity used in (1.1) and $\bar{\Omega}$ be the time-frequency of the numerical solution $\bar{W}$. For all $n\in\mathbb{Z}^{4}$ and $n\in\mathbb{N}$, we define

[TABLE]

and the set:

[TABLE]

We then fix $N^{\dagger}\in\mathbb{N}\setminus\{0\}$ (to be chosen later) and consider the set $\mathcal{E}^{\dagger}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\mathcal{E}(N^{\dagger})$ corresponding to the subset of indices that are used in the definitions of $\widehat{A}$ and $A$ below. The reasoning for choosing such a set $\mathcal{E}^{\dagger}$ will become apparent later, when we have to control the dominant linear part of $\mathcal{F}$ (see for instance Section 3.4). This choice is also quite natural from a spectral viewpoint. Indeed $\mathcal{E}^{\dagger}$ is the set of all indices corresponding to eigenvalues of the heat operator $\partial_{t}-\nu\Delta$ with modulus less than $N^{\dagger}$ (see Lemma 3.1).

Definition 2.14**.**

We define the subspace $\mathcal{X}^{\dagger}$ of $\mathcal{X}$ as

[TABLE]

For $a\in\mathbb{C}^{\mathbb{Z}^{4}_{*}}$, we define

[TABLE]

For $W=(\Omega,\omega)\in\mathcal{X}$, this notation naturally extends to

[TABLE]

In the sequel, we identify the finite dimensional vector $\Pi^{\dagger}W$ with its natural injection into $\mathcal{X}^{\dagger}$, and therefore interpret $\Pi^{\dagger}$ as the canonical projection from $\mathcal{X}$ to $\mathcal{X}^{\dagger}$. These notations also naturally extend to $\mathcal{X}_{-1,-2}$.

We are now ready to introduce an approximation $\widehat{A}$ of $D\mathcal{F}(\bar{W})$ and then an approximation $A$ of $D\mathcal{F}(\bar{W})^{-1}$. The bounded linear operator $\widehat{A}:\mathcal{X}\to\mathcal{X}_{-2,-1}$ is defined by

[TABLE]

Notice that $\widehat{A}$ leaves both subspaces $\mathcal{X}^{\dagger}$ and $(I-\Pi^{\dagger})\mathcal{X}$ invariant, and that it acts diagonally on $(I-\Pi^{\dagger})\mathcal{X}$.

Next, we introduce $A^{(N^{\dagger})}$, an approximate inverse of $\Pi^{\dagger}D\mathcal{F}(\bar{W})|_{\mathcal{X}^{\dagger}}$ that is computed numerically (interpreting it as a finite matrix), and the bounded linear operator $A:\mathcal{X}_{-2,-1}\to\mathcal{X}$ defined by

[TABLE]

where

[TABLE]

Notice that $A$ also leaves both subspaces $\mathcal{X}^{\dagger}$ and $(I-\Pi^{\dagger})\mathcal{X}$ invariant, and that it acts diagonally on $(I-\Pi^{\dagger})\mathcal{X}$.

2.4 A posteriori validation framework

We are now ready to give sufficient conditions for the a posteriori validation of the solution $\bar{W}$, that is conditions under which the existence of a zero of $\mathcal{F}$ in a neighborhood of $\bar{W}$ is guaranteed. This is the content of the following theorem.

Theorem 2.15**.**

Let $\eta>1$. With the notations of the previous sections, assume there exist $\bar{W}\in\mathcal{X}$ and non-negative constants $Y_{0}$, $Z_{0}$, $Z_{1}$ and $Z_{2}$ such that

[TABLE]

Assume also that

• •

the forcing term $f$ is time independent and has space average zero;

• •

$\bar{W}$* is in $\mathcal{X}^{\text{{div}}}$;*

• •

$\hat{\omega}$* (used to define the phase condition (2.13)) and $\bar{W}$ are such that $\gamma_{*}\hat{\omega}=\hat{\omega}$ and $\gamma_{*}\bar{W}=\bar{W}$.*

If

[TABLE]

then, for all $r\in[r_{\min},r_{\max})$ there exists a unique $\tilde{W}=(\tilde{\Omega},\tilde{\omega})\in\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ such that $\mathcal{F}(\tilde{W})=0$, where $\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ is the closed ball of ${\mathcal{X}}$, centered at $\bar{W}$ and of radius $r$, and

[TABLE]

Besides, this unique $\tilde{W}$ also lies in $\mathcal{X}^{\text{{div}}}$. Finally, defining $u=M\tilde{\omega}$, there exists a pressure function $p$ such that $(u,p)$ is a $\frac{2\pi}{\tilde{\Omega}}$-periodic, real valued and analytic solution of Navier-Stokes equations (1.1).

Proof.

First notice that we have

[TABLE]

hence $A\widehat{A}$ is a bounded linear and invertible operator from $\mathcal{X}$ to itself by a standard Neumann series argument. Besides, since $A$ and $\widehat{A}$ have diagonal tails that are the exact inverse of one another, the above bound gives that the finite part of $A$ is invertible and thus $A^{-1}$ is a bounded linear operator from $\mathcal{X}_{-2,-1}$ to $\mathcal{X}$. We also have that

[TABLE]

and therefore $Q\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}AD\mathcal{F}(\bar{W})$ is a bounded linear and invertible operator from $\mathcal{X}$ to itself. Furthermore, $\left\|Q^{-1}\right\|_{B(\mathcal{X},\mathcal{X})}\leq(1-(Z_{0}+Z_{1}))^{-1}$. Hence, we finally have that

[TABLE]

is a bounded linear and invertible operator from $\mathcal{X}_{-2,-1}$ to $\mathcal{X}$.

We can thus consider

[TABLE]

which maps $\mathcal{X}$ to itself. We are now going to show that $T$ is a contraction on $\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ for all $r\in[r_{\min},r_{\max})$. It is going to be helpful to introduce the polynomial

[TABLE]

Notice that by (2.21), the quadratic polynomial $P$ has two positive roots, $r_{\min}$ being the smallest, and that $P(r_{\max})<0$ since $r_{\max}$ is the apex of $P$. We estimate, for $r>0$ and $W\in\mathcal{B}_{\mathcal{X}}(\bar{W},r)$,

[TABLE]

For all $r\in[r_{\min},r_{\max}]$, we have $P(r)\leq 0$ and hence $\left\|T(W)-\bar{W}\right\|_{\mathcal{X}}\leq r$, that is $T$ maps $\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ into itself.

Furthermore, for $r>0$ and $W\in\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ we also have

[TABLE]

Hence, from the definition of $r_{\max}$ it follows that $\left\|DT(W)\right\|_{B(\mathcal{X},\mathcal{X})}<1$ for any $r\in[0,r_{\max})$.

Thus we infer that $T$ is a contraction on $\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ for any $r\in[r_{\min},r_{\max})$. Banach’s fixed point Theorem then yields the existence of a unique fixed point $\tilde{W}=(\tilde{\Omega},\tilde{\omega})$ of $T$ in $\mathcal{B}_{\mathcal{X}}(\bar{W},r)$, for all $r\in[r_{\min},r_{\max})$, which corresponds to a unique zero of $\mathcal{F}$ since $D\mathcal{F}(\bar{W})^{-1}$ is injective.

In order to prove that $\tilde{\omega}$ is divergence free, we then make use of the fact that $\mathcal{F}(\mathcal{X}^{\text{{div}}})\subset\mathcal{X}_{-2,-1}^{\text{{div}}}$ (Lemma 2.9). Since $\mathcal{X}_{-2,-1}^{\text{{div}}}$ is a closed subspace of $\mathcal{X}_{-2,-1}$, this implies that for all $W\in\mathcal{X}^{\text{{div}}}$, $D\mathcal{F}(W)\left(\mathcal{X}^{\text{{div}}}\right)\subset\mathcal{X}_{-2,-1}^{\text{{div}}}$. In particular, since $D\mathcal{F}(\bar{W})$ is invertible and $\bar{W}\in\mathcal{X}^{\text{{div}}}$, we have that $D\mathcal{F}(\bar{W})^{-1}\left(\mathcal{X}_{-2,-1}^{\text{{div}}}\right)\subset\mathcal{X}^{\text{{div}}}$, and hence that $T(\mathcal{X}^{\text{{div}}})\subset\mathcal{X}^{\text{{div}}}$. Therefore, $T$ is also a contraction on $\mathcal{B}_{\mathcal{X}^{\text{{div}}}}(\bar{W},r)\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\mathcal{B}_{\mathcal{X}}(\bar{W},r)\cap\mathcal{X}^{\text{{div}}}$ and we now get the existence of a unique fixed point $\tilde{W}$ of $T$ in $\mathcal{B}_{\mathcal{X}^{\text{{div}}}}(\bar{W},r)$ (which must be the same as the one obtained previously by uniqueness).

Finally, concerning real-valuedness, by Lemma 2.11 we have $\mathcal{F}(\gamma_{*}\tilde{W})=\gamma_{*}\mathcal{F}(\tilde{W})=0$, and since $\mathcal{B}_{\mathcal{X}}(\bar{W},r)$ is invariant under $\gamma_{*}$, this “new” zero $\gamma_{*}\tilde{W}$ of $\mathcal{F}$ also belongs to $\mathcal{B}_{\mathcal{X}}$. By uniqueness we must have $\gamma_{*}\tilde{W}=\tilde{W}$, which means that that $\Omega\in\mathbb{R}$ and that the functions associated to $\tilde{\omega}$ and $u$ are real-valued.

The function associated to $\tilde{\omega}$ is therefore divergence free, analytic since $\eta>1$, and $\Omega\in\mathbb{R}$. By Lemma 2.5, $u$ then solves Navier-Stokes equations. ∎

Remark 2.16**.**

Many similar versions of this theorem have been used in the last decades, in a posteriori error analysis and computer-assisted proofs (see for instance [4, 26, 47, 2, 6]). One possible approach (which is maybe the most natural one) to show that the operator $T$ defined in (2.22) is a contraction from a small ball around $\bar{W}$ into itself, is to directly estimate

[TABLE]

instead of (2.17)-(2.20). The main difficulty of this approach in practice is to obtain a bound for the inverse $D\mathcal{F}(\bar{W})^{-1}$, which can sometimes be done using eigenvalue enclosing techniques (see [26] and the references therein). Another possibility is to replace $D\mathcal{F}(\bar{W})^{-1}$ by an approximate inverse $A$ of $D\mathcal{F}(\bar{W})$, and to study the fixed point operator

[TABLE]

instead of $T$. Estimating the quantities in (2.17)-(2.20) is then a good way to prove that $T_{\text{{approx}}}$ is a contraction from a small ball around $\bar{W}$ into itself (see for instance [47, 6]). The price to pay for this approach is that one has to actually compute (partially numerically) a good enough approximate inverse $A$ of $D\mathcal{F}(\bar{W})$, but estimating the norm of $A$ then becomes very straightforward, compared to having to work with the exact inverse $D\mathcal{F}(\bar{W})^{-1}$.

What Theorem 2.15 shows is that, when an approximate inverse $A$ is used and when the estimates (2.17)-(2.20) are good enough to show that $T_{\text{{approx}}}$ is a contraction from $\mathcal{B}(\bar{W},r)$ into itself for some $r>0$, then the same is actually true for $T$. This observation may seem inconsequential, as what we really care about is the existence of a zero of $\mathcal{F}$ near $\bar{W}$, rather than which fixed point operator was used to prove this fact (since both give the same error bound with this approach). However, it turns out to be very advantageous in this work. Indeed, to show that a zero $W$ of $\mathcal{F}$ corresponds to a solution of Navier-Stokes equations, we need to know a priori that $\omega$ is divergence free, see Lemma 2.5. Since $\mathcal{F}$ (and thus $D\mathcal{F}(\bar{W})^{-1}$) preserves the divergence free property, so does $T$, and hence we were able to obtain for free that our fixed point is divergence free (by making sure that the numerical approximate solution $\bar{W}$ itself is divergence free). To obtain the same result with $T_{\text{{approx}}}$, we would have to make sure that the approximate inverse $A$ also preserves exactly the divergence free property. For the approximate inverse $A$ described in Section 2.3 this could likely be achieved, albeit at a considerable computational cost. However, as mentioned in the introduction, in practice we work with symmetry reduced variables (see Section 4). Therefore we use a symmetrically reduced version of $A$ instead of $A$ itself. Making sure that this symmetrically reduced version of $A$ preserves the divergence free property without directly having access to $A$ is a formidable task, which we are able to avoid by working with $T$ instead of $T_{\text{{approx}}}$.

Remark 2.17**.**

We note that the weighted $\ell^{1}$-norm controls the $\mathcal{C}^{0}$-norm. Therefore, when the assumptions of Theorem 2.15 are satisfied, we have the following explicit error control between the exact vorticity $\tilde{\omega}$ and the approximate vorticity $\bar{\omega}$:

[TABLE]

where $\theta=\tilde{\Omega}/\bar{\Omega}$ accounts for the time dilation that occurs when the approximate period is not exactly equal to the true period. We also point out that, for $\eta>1$, we could also obtain explicit error estimates on derivatives of the vorticity. In Lemma 6.2 corresponding error bounds on the velocity field and the pressure are presented.

While Theorem 2.15 is the cornerstone of our approach, the main difficulty still lies ahead of us: we need to derive and implement explicit bounds satisfying (2.17)-(2.20), that are sharp enough for (2.21) to hold. These bounds are obtained in Section 3. However, the computations required to evaluate them are quite prohibitive due to the high dimension of the problem. Therefore, in Section 4 we make use of the symmetries of the solution to reduce the amount of computation needed, and update the bounds obtained in Section 3 accordingly, by showing that they are in a certain precise sense compatible with the symmetries. The implementation in a MATLAB code of the bounds in the symmetric setting can be found at [41], and the results are discussed in Section 5.

3 Estimates without using symmetries

In this section, we derive bounds $Y_{0}$, $Z_{0}$, $Z_{1}$ and $Z_{2}$ satisfying (2.17)-(2.20). We first list a few auxiliary lemmas that are going to be used several times, and then devote one subsection to each bound. We assume throughout that the approximate solution $\bar{W}$ only has a finite number of non-zero modes. More precisely, writing $\bar{W}=(\bar{\Omega},\bar{\omega})\in\mathbb{C}\times\mathbb{C}^{\mathbb{Z}^{4}_{*}}$, we assume there exists a finite set of indices $\mathcal{S}^{\text{{sol}}}\subset\mathbb{Z}^{4}_{*}$ such that

[TABLE]

We also assume that $\hat{\omega}$ used in the phase condition (2.13) is chosen so that $\hat{\omega}_{n}=0$ for all $n\notin\mathcal{S}^{\text{{sol}}}$, that $\mathbb{N}^{\dagger}\in N\setminus\{0\}$ is fixed, and recall that $\mathcal{E}^{\dagger}$ is defined in Section 2.3.

3.1 Uniform estimates involving $\lambda_{n}$

We regroup here some straightforward lemmas that are going to facilitate estimates involving the tail part of $A$, in Sections 3.4 and 3.5. Recall that $\mathcal{E}(N)$ was introduced in Definition 2.13 and $\lambda_{n}$ in (2.16).

Lemma 3.1**.**

For all $N\in\mathbb{N}\setminus\{0\}$

[TABLE]

Proof.

Simply notice that $\mu(n)=\frac{1}{\left|\lambda_{n}\right|}$. ∎

Lemma 3.2**.**

For all $N\in\mathbb{N}\setminus\{0\}$ and all $m\in\{1,2,3\}$,

[TABLE]

Proof.

We note that

[TABLE]

hence

[TABLE]

Lemma 3.3**.**

For all $N\in\mathbb{N}\setminus\{0\}$

[TABLE]

Besides, for all $p\in\{1,2,3\}$, we also have

[TABLE]

Proof.

Just notice that

[TABLE]

where the first inequality follows from the Cauchy-Schwarz inequality. Therefore

[TABLE]

The second estimate is obtained similarly, using that

[TABLE]

Lemma 3.4**.**

For all $N\in\mathbb{N}\setminus\{0\}$

[TABLE]

Proof.

Let $n\notin\mathcal{E}(N)$. If $|n|_{\infty}=|n_{4}|$, then clearly

[TABLE]

Otherwise, $|n|_{\infty}=|n_{m}|$ for $m\in\{1,2,3\}$ and we conclude by Lemma 3.2. ∎

3.2 $Y_{0}$ bound for $A\mathcal{F}(\bar{W})$

We start by giving a computable bound for $\left\|A\mathcal{F}(\bar{W})\right\|_{\mathcal{X}}$.

Proposition 3.5**.**

Assume that $f^{\omega}_{n}=0$ for $n\notin\mathcal{S}^{\text{{sol}}}+\mathcal{S}^{\text{{sol}}}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\{n_{1}+n_{2}\ |\ n_{1},n_{2}\in\mathcal{S}^{\text{{sol}}}\}$, and that $A$ is defined as in Section 2.3. Then (2.17) is satisfied with

[TABLE]

Proof.

The only thing to notice is that the approximate solution $\bar{W}$ only has a finite number of non-zero modes, hence the same is true for $\mathcal{F}(\bar{W})$. More precisely, since $F$ is quadratic we infer that $F_{n}(\bar{W})=0$ for all $n\notin\mathcal{S}^{\text{{sol}}}+\mathcal{S}^{\text{{sol}}}$. Recalling the definition of $A$ and using the splitting

[TABLE]

then directly yields $Y_{0}$. ∎

3.3 ${Z_{0}}$ bound for ${I-A\widehat{A}}$

We now give a computable bound for $\left\|I-A\widehat{A}\right\|_{B(\mathcal{X},\mathcal{X})}$.

Proposition 3.6**.**

Assume that $\widehat{A}$ and $A$ are defined as in Section 2.3. Then (2.18) is satisfied with

[TABLE]

Proof.

We recall that we defined $\widehat{A}$ and $A$ in such a way that their tails are exact inverses of each other, that is

[TABLE]

Therefore, the only non-zero part of $I-A\widehat{A}$ is the finite part $\Pi^{\dagger}\left(I-A\widehat{A}\right)\Pi^{\dagger}$, which yields $Z_{0}$. ∎

Besides, for a linear operator $B:\mathcal{X}\to\mathcal{X}$, the operator norm of $B$ is nothing but the supremum of the norm of each of its column, with a weight, since we use a weighted $\ell^{1}$ norm on $\mathcal{X}$:

[TABLE]

Therefore, if $B$ only has a finite number of non-zero columns $B^{(.,m)}_{.,n}$, and if each of these columns only has a finite number of non-zero terms, (which is the case for the operator involved in $Z_{0}$) we can evaluate such a norm on a computer.

3.4 ${Z_{1}}$ bound for {A\bigl{(}D\mathcal{F}(\bar{W})-\widehat{A}\bigr{)}}

In this section, we give a computable bound for \bigl{\|}A\bigl{(}D\mathcal{F}(\bar{W})-\widehat{A}\bigr{)}\bigr{\|}_{B(\mathcal{X},\mathcal{X})}. Because of the way $\widehat{A}$ and $A$ are defined, and since we only consider here the derivative of $\mathcal{F}$ at a numerical solution $\bar{W}$ (which only has a finite number of nonzero coefficients), each column of A\bigl{(}D\mathcal{F}(\bar{W})-\widehat{A}\bigr{)} also only has a finite number of nonzero coefficients. Therefore, we can numerically evaluate the norm of any finite number of columns of A\bigl{(}D\mathcal{F}(\bar{W})-\widehat{A}\bigr{)}. Our strategy to obtain the $Z_{1}$ bound is thus to compute the norm of a finite (but large enough) number of columns of A\bigl{(}D\mathcal{F}(\bar{W})-\widehat{A}\bigr{)}, and to then get an analytic estimate for the remaining columns. To describe for which columns we compute the norm explicitly, we introduce the following set of indices.

Definition 3.7**.**

Let $\tilde{N}\in\mathbb{N}$. We define the set $\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$ by

[TABLE]

Remark 3.8**.**

In the sequel, we will have to estimate quantities such as

[TABLE]

for a given set $\mathcal{S}$. By definition (3.2) of $\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$, and assuming $-\mathcal{S}^{\text{{sol}}}=\mathcal{S}^{\text{{sol}}}$, we have

[TABLE]

hence for $\mathcal{S}=\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$ we can bound (3.3) by

[TABLE]

and we can thus use the estimates of Section 3.1 to control such terms.

Remark 3.9**.**

Whereas the computational parameter $N^{\dagger}$ determines the size of the “computational/finite” part $A^{(N^{\dagger})}$ of the linear operator $A$, the computational parameter $\tilde{N}$ can be chosen (for fixed choice of $N^{\dagger}$) to balance the quality of the estimates and the computational costs. We will always need to choose $\tilde{N}$ so that $\mathcal{E}(\tilde{N})$ contains $\mathcal{E}^{\dagger}=\mathcal{E}(N^{\dagger})$, i.e. take $\tilde{N}\geq N^{\dagger}$, to ensure that the finite part of $A$ does not influence the tail estimate, see equation (3.4) below.

Proposition 3.10**.**

Assume that $\bar{W}\in\mathcal{X}^{\text{{div}}}$, $\mathcal{E}^{\dagger}\subset\mathcal{E}(\tilde{N})$, and that $\widehat{A}$ and $A$ are defined as in Section 2.3. Define $C=D\mathcal{F}(\bar{W})-\widehat{A}$. Then (2.19) is satisfied with

[TABLE]

where

[TABLE]

and

[TABLE]

where $\sqrt{3}$ can be replaced by $\sqrt{2}$ for an (essentially) 2D solution.

Proof.

We have

[TABLE]

and we want to show that this quantity is bounded by $Z_{1}$. From the definition of $Z_{1}^{\textup{finite}}$ and $Z_{1}^{\textup{tail}}$, and the splitting

[TABLE]

we see that we only have to prove that

[TABLE]

In order to do so, it is helpful to first have a look at the structure of $C$.

By definition of $\widehat{A}$, we have that $\Pi^{\dagger}C\Pi^{\dagger}=0$, and that the remaining coefficients are given by $D\Psi(\bar{\omega})$, since the linear part of $D\mathcal{F}(\bar{W})$ is also present in $\widehat{A}$, that is: if $k\notin\mathcal{E}^{\dagger}$ or $n\notin\mathcal{E}^{\dagger}$ then

[TABLE]

with $\Psi$ as defined in (2.11). Using the block-notation introduced in Section 2.2, we write

[TABLE]

with, for all $1\leq l,m\leq 3$,

[TABLE]

and

[TABLE]

Furthermore, using that for all $k\in\mathbb{Z}^{4}_{*}$

[TABLE]

we obtain

[TABLE]

where $M^{(l,m)}_{n}$ is the coefficient on row $l$ and column $m$ in the matrix $M_{n}$, see (2.6). Each of the four terms in $C^{(l,m)}_{k,n}$ comes from one of the four terms in $D\Psi(\bar{\omega})\omega$, and it should be noted that we used Lemma 2.2 on the first of those four terms to obtain this expression (see Remark 3.11 below).

From the expression of $C$ we see that $\left\|AC^{(.,\leftmoon)}\right\|_{\mathcal{X}}$ and $\max_{1\leq m\leq 3}\left(Z_{1}^{\textup{finite}}\right)^{(m)}$ can indeed be evaluated with a computer, since we only consider a finite number of columns, each having only a finite number of non-zero components. In particular, for any $n\in\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$ the coefficient $C_{k,n}^{(l,m)}$ vanishes for $k$ outside $\mathcal{E}(\tilde{N})+2\mathcal{S}^{\text{{sol}}}$.

We now focus on proving (3.5). First, we recall that $\bar{\omega}_{k-n}=0$ if $k-n\notin\mathcal{S}^{\text{{sol}}}$. By definition of $\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$ this means that for any $1\leq l,m\leq 3$

[TABLE]

In particular, provided $\mathcal{E}^{\dagger}\subset\mathcal{E}(\tilde{N})$, see Remark 3.9, for all $1\leq m\leq 3$ and $n\notin\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$ the non-zero coefficients of the column $C^{(.,m)}_{.,n}$ only get hit by the tail part of $A$. Therefore, we obtain for all $1\leq l,m\leq 3$ and $n\notin\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$

[TABLE]

We are going to bound the supremum over $n\notin\widetilde{\mathcal{S}}^{\text{{sol}}}(\tilde{N})$ of each of the four terms in the sum over $l$ separately. For the first term, using Lemma 3.3 we estimate

[TABLE]

Notice that, if we have an (essentially) $2D$ solution (see Section 1), then one component of $M\bar{\omega}$ is zero, which allows us to use the second estimate of Lemma 3.3, to get

[TABLE]

For the second term, using Lemma 3.1 instead of Lemma 3.3, we get

[TABLE]

For the third term, again using Lemma 3.1, as well as the fact that $M_{n}^{(p,m)}=0$ if $p=m$ and $|M_{n}^{(p,m)}|\leq 1$ if $p\neq m$, we infer that

[TABLE]

Finally, to bound the fourth term, we observe that

[TABLE]

Hence

[TABLE]

Putting everything together, we get that

[TABLE]

where $\sqrt{3}$ can be replaced by $\sqrt{2}$ for a 2D solution, and (3.5) is proven, which concludes the proof of Proposition 3.10. ∎

Remark 3.11**.**

To obtain the $Z_{1}$ bound, we have had to estimate $\left\|AC\right\|_{B(\mathcal{X},\mathcal{X})}$. Since $C$ is unbounded as an operator from $\mathcal{X}$ to itself, we have to carefully handle the unbounded terms appearing in $C$, or equivalently in the derivative $D\Psi(\bar{\omega})$. Indeed, the unboundedness of $C$ is compensated by the decaying tail of $A$. Looking at (3.6), we see that the derivatives in the second and the fourth term are not bothersome, as they are “balanced” by $M$ (which, very roughly speaking, acts as an antiderivative). The derivative in the third term of (3.6) is also not an issue, since it applies to $\bar{\omega}$ which has only a finite number of non zero coefficients (hence its derivative still belongs to $\mathcal{X}$ and its norm can be computed explicitly). However, something extra has to be done to control the derivative appearing in the first term of (3.6). Using Lemma 2.2 on this specific term allows us to factor out the derivative from the convolution product, and this derivative is then canceled by the tail part of $A$, which also acts like a kind of antiderivative (see (3.10) for the actual estimate). This explains why we write $\sum_{p=1}^{3}k_{p}(M\bar{\omega})^{(p)}_{k-n}$ rather than the equivalent (in view of $M\bar{\omega}$ being divergence free) $\sum_{p=1}^{3}n_{p}(M\bar{\omega})^{(p)}_{k-n}$ in the first term of (3.7).

3.5 ${Z_{2}}$ bound for {A\bigl{(}D\mathcal{F}(W)-D\mathcal{F}(\bar{W})\bigr{)}}

In this section, we present the last missing estimate for Theorem 2.15, namely a computable bound for $\left\|A\left(D\mathcal{F}(W)-D\mathcal{F}(\bar{W})\right)\right\|_{B(\mathcal{X},\mathcal{X})}$ in term of $\left\|W-\bar{W}\right\|_{\mathcal{X}}$.

Proposition 3.12**.**

Assume that $A$ is defined as in Section 2.3. Then (2.20) is satisfied with

[TABLE]

where

[TABLE]

Proof.

Let $W^{\prime},W^{\prime\prime}\in\mathcal{X}$ and consider

[TABLE]

Proving the proposition amounts to showing that

[TABLE]

By bi-linearity of $D^{2}\mathcal{F}(\bar{W})$, it is enough to assume that $\left\|W^{\prime}\right\|_{\mathcal{X}},\left\|W^{\prime\prime}\right\|_{\mathcal{X}}\leq 1$, and show that $\left\|Az\right\|_{\mathcal{X}}\leq Z_{2}$. We start by taking a closer look at $z$, which can be expanded as

[TABLE]

where

[TABLE]

Since the first component of $z$ vanishes, we have

[TABLE]

where $\mathcal{X}_{-1,-1}^{0}$ is the subspace of $\mathcal{X}_{-1,-1}$ defined as $\{0\}\times(\ell^{1}_{\eta,-1,-1}(\mathbb{C}))^{3}$ (that is, the first column $A^{(\cdot,\leftmoon)}$ of $A$ does not play any role, since we multiply $A$ with a vector whose first component vanishes).

We now proceed to estimate both terms in the right-hand side of (3.13). For the first one, we infer from Lemma 3.4 that

[TABLE]

To bound the second term, we now use the splitting of $z$ introduced earlier and estimate each term separately. We start by

[TABLE]

Next, we want to estimate

[TABLE]

Thanks to Lemma 2.2, we can rewrite

[TABLE]

where $(M\omega^{\prime\prime})^{(p)}\ast\omega^{\prime}$ must be understood as

[TABLE]

For each $p\in\{1,2,3\}$, we estimate

[TABLE]

and using that

[TABLE]

we get

[TABLE]

Similarly, we get that

[TABLE]

Finally, we have to estimate

[TABLE]

Starting from

[TABLE]

we estimate

[TABLE]

Up to a permutation of the indices $1$, $2$ and $3$, we have for any $l,m,p\in\{1,2,3\}$ that either

[TABLE]

or

[TABLE]

for all $\tilde{n}\in\mathbb{Z}^{3}\setminus\{0\}$. Therefore

[TABLE]

and we end up with

[TABLE]

Similarly, we have

[TABLE]

Adding everything up, we end up with

[TABLE]

which concludes the proof of Proposition 3.12 ∎

Remark 3.13**.**

It should be noted that the slightly unusual step (3.13) is crucial for what is to come in Section 4.6.4. Usually, one would try to postpone this step to directly “cancel” some of the unbounded derivative operators in $z$ with $A$, by first splitting $z$ into several terms and estimating then like

[TABLE]

However, this splitting may not conserve some of the symmetries that $z$ has as a whole, and therefore would not be compatible with the symmetry reduced variables used later on.

4 The symmetries

In this section, we introduce the formalism that allows us to take advantage of the symmetries through the a posteriori validation procedure. In Section 4.1, we introduce some group actions that we use to explicitly describe the symmetries that we are considering, in the physical space (that is symmetries acting on the velocity $u$ or the vorticity $\omega$, seen as functions). In Section 4.2, we proceed to explain why the Navier-Stokes equations are equivariant under these group actions. Equivalent group actions in Fourier space are then defined in Section 4.3, and the equivariance of $\mathcal{F}$ is shown. In Section 4.4, we introduce the subspace $\mathcal{X}^{\text{{sym}}}$ of $\mathcal{X}$ containing the symmetric solutions, as well as a minimal (in term of number of Fourier mode used) subspace $\mathcal{X}^{\text{{red}}}$ of $\mathcal{X}$ that can be used to describe these solutions. A reduced version $\mathcal{F}^{\text{{red}}}$ of $\mathcal{F}$ is then defined on this subspace, and sufficient conditions for the existence of zeros of $\mathcal{F}^{\text{{red}}}$ are then given in Theorem 4.23, which mimics Theorem 2.15. Finally, in Section 4.6, we show that the bounds obtained in Section 3 are compatible with the symmetries, and readily give us the estimates needed to apply Theorem 4.23, which is then done in Section 5.

4.1 Symmetries in physical space

We will restrict our attention to solutions $(u,p)$, $u:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}^{3}$ and $p:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}$ which are $2\pi$-periodic in space and $2\pi/\Omega$-periodic in time. We may thus interpret $u$ and $p$ as a function on $\mathbb{T}^{3}\times\mathbb{S}^{1}$, where $\mathbb{T}^{3}$ is a 3-dimensional torus and $\mathbb{S}^{1}=\mathbb{R}/(2\pi/\Omega)$ is a circle.

We look for solutions with additional symmetries. In particular, we will consider periodic solutions which are invariant under a symmetry group $G$, which acts on $(u,p)$ through a right action $a_{g}$, $g\in G$, of the form

[TABLE]

where $b_{g}$ is a right action on $\mathbb{R}^{3}$, while $c_{g}$ is a left action on $\mathbb{T}^{3}$ and $d_{g}$ is a left action on $\mathbb{S}^{1}$. We abuse the notation slightly by using $a_{g}$ to denote both the action on the pair $(u,p)$ and on each of its components, but it should be clear from context which one is considered.

The symmetries under consideration in this work lead to actions of the form

[TABLE]

where $C_{g}$ is a $3\times 3$ signed permutation matrix (it represents an element of $O_{h}$, the symmetry group of the cube): there is a permutation $\tau_{g}$ of $\{1,2,3\}$ and a map $\rho_{g}:\{1,2,3\}\to\{-1,1\}$ such that $(C_{g})_{mm^{\prime}}=\rho_{g}(m^{\prime})\delta_{m\tau_{g}(m^{\prime})}$, $\tilde{C}_{g}$ is a vector in $\mathbb{R}^{3}$ and $D_{g}$ a real number. An explicit description of the symmetries satisfied by the solutions considered in this work, together with the associated $C_{g}$, $\tilde{C}_{g}$ and $D_{g}$ is given in Section 5.2.

We now describe how these actions behave with respect to differential operators.

Lemma 4.1**.**

*Let $u:\mathbb{R}^{4}\to\mathbb{R}^{3}$ and $p:\mathbb{R}^{4}\to\mathbb{R}$ be smooth functions, and assume $a_{g}$ is defined as in (4.1), with $C_{g}$ an orthogonal matrix. Then

     (a) $\nabla a_{g}p=a_{g}\nabla p$,

     (b) $\nabla\cdot a_{g}u=a_{g}\nabla\cdot u$,

     (c) $\Delta a_{g}u=a_{g}\Delta u$,

     (d) $(a_{g}u\cdot\nabla)a_{g}u=a_{g}(u\cdot\nabla)u$,

     (e) $\nabla\times a_{g}u=\det(C_{g})a_{g}\nabla\times u$.*

Proof.

The first four identities follow directly from the chain rule, using that the change of variable $x\mapsto C_{g}x+2\pi\tilde{C}_{g}$ transforms $\nabla$ into $C^{T}_{g}\nabla$. The last one is slightly less straightforward. Writing $y=C_{g}x+2\pi\tilde{C}_{g}$, we have

[TABLE]

For any $v\in\mathbb{R}^{3}$ we then compute

[TABLE]

Since this holds for all $v\in\mathbb{R}^{3}$, we infer

[TABLE]

which yields identity (e). ∎

From the action $a_{g}$ on $u$, we can define a corresponding action $\tilde{a}_{g}$ on $\omega=\nabla\times u$, given by

[TABLE]

From (4.1) and Lemma 4.1 we infer

[TABLE]

We note that $\det(C_{g})=\pm 1$, since $C_{g}$ is orthonormal, hence the only difference between the action of $a_{g}$ and $\tilde{a}_{g}$ is multiplication by a factor $\pm 1$ (depending on $g$).

4.2 Invariance and equivariance

We say that $u$ is $G$-invariant under the action $a_{g}$ if

[TABLE]

Similarly, we say that $\omega$ is $G$-invariant under the action $\tilde{a}_{g}$ if

[TABLE]

If it is clear which action is meant between $\tilde{a}_{g}$ and $a_{g}$, we only speak about $G$-invariance without mentioning the action. Notice that, by Lemma 4.1(b), the collection of divergence free vector fields is invariant under $G$. Analogously, having zero spatial average, denoted by $\int_{\mathbb{T}^{3}}u=0$, is also a $G$-invariant property. Since the Laplacian is invertible on the set of vector fields with zero spatial average, and it commutes with the group action (Lemma 4.1(c)), it follows from Lemma 4.1(e) that

[TABLE]

Next, we mention the equivariance properties of the Navier-Stokes equations that provide intuition for the sequel. To formalize the discussion, we define the pair $\mathcal{P}=(\mathcal{P}_{1},\mathcal{P}_{2})$ of maps

[TABLE]

on some set $\mathbb{X}$ of smooth vector fields $u$ and scalar functions $p$. Assuming that $f$ is $G$-invariant under the action $a_{g}$, it follows from Lemma 4.1 that $\mathcal{P}$ is $G$-equivariant:

[TABLE]

We note that $\Phi=\Phi(u)$ defined in (2.12) is given by $\mathcal{P}_{1}(u,0)$, hence $\Phi$ is also $G$-equivariant: $\Phi(a_{g}u)=a_{g}\Phi(u)$.

Next, we consider the vorticity formulation and set

[TABLE]

It then follows from (4.2), (4.4) and the equivariance of $\Phi$ that $\mathcal{Q}$ is $G$-equivariant under the action $\tilde{a}_{g}$. As before, we simply refer to $G$-equivariance without mentioning the action, if the latter is clear from the context. The explicit expression of $\mathcal{Q}$ is given by

[TABLE]

Notice that this expression is only equal to $F(W)$ when $\omega$ is divergence free, hence a supplemental argument is required to show that $F$ is $G$-equivariant, which we provide next.

Lemma 4.2**.**

$F$* defined is (2.10) is $G$-equivariant, i.e.*

[TABLE]

where $\tilde{a}_{g}W=\tilde{a}_{g}(\Omega,\omega)$ is defined as $(\Omega,\tilde{a}_{g}\omega)$.

Proof.

Looking at (2.7), the statement follows from (4.4) and Lemma 4.1. ∎

Finally, obviously $\mathcal{P}_{1}(u,p)$ is $G$-invariant if both $u$ and $p$ are $G$-invariant, and similarly for $\mathcal{Q}(\omega)$ and $F(W)$.

Remark 4.3**.**

Concerning the pressure term in the Navier-Stokes equation, in the appendix we describe explicitly a well-defined map $\Gamma$ (see Lemma 6.1) such that

[TABLE]

From this it follows that $[\Gamma a_{g}\Phi](x,t)=[\Gamma\Phi](c_{g}x,d_{g}t)$. Hence, when recovering the pressure from the vorticity, see the proof of Lemma 2.5, we can infer that $G$-invariance of $\omega$ implies $G$-invariance of $\Phi(M\omega)$, from which in turn $G$-invariance of $p=\Gamma\Phi(M\omega)$ follows.

4.3 Representation of the symmetry group in Fourier space

We recall (2.9) that we write the Fourier transform on the $m$-th component of the vorticity as

[TABLE]

where $\omega_{0}$ is assumed to be zero. We note that, as in Section 2, it should be clear from the context whether $\omega$ is to be interpreted in physical or Fourier space. From now on we will denote ${j}=(n,m)\in\mathbb{Z}^{4}_{*}\times\{1,2,3\}$ and write

[TABLE]

We define the index set

[TABLE]

and the set of Fourier coefficients by

[TABLE]

where, for convenience, we introduce the weights

[TABLE]

Remark 4.4**.**

Notice that the space $\mathcal{X}$ defined is Section 2.2 is nothing but $\mathbb{C}\times X$. When looking at the action $\tilde{a}_{g}$ in Fourier space, it is enough to consider $X$ because the action does neither act on $\Omega$ nor depend on it (see Remark 4.6). The reason for grouping $n$ and $m$ as a duo in $j=(n,m)$ is that the group action on the Fourier coefficients (made explicit below) permutes the index duos $j\in J$ rather than $n\in\mathbb{Z}^{4}_{*}$ and $m\in\{1,2,3\}$ individually, see (4.11).

We define the basis vectors $\mathrm{e}_{j}\in X$ by

[TABLE]

For $g\in G$, the group actions $a_{g}$ and $\tilde{a}_{g}$ described in Sections 4.1 and 4.2 correspond to a right group action $\gamma_{g}$ on the Fourier coefficients satisfying

[TABLE]

This action is represented component-wise by

[TABLE]

Here $\beta_{g}$ is itself a left group action on $J$, i.e.,

[TABLE]

whereas $\alpha_{g}({j})\in\{z\in\mathbb{C}:|z|=1\}$ for all ${j}\in J$. The product structure on $\alpha$ is given by

[TABLE]

which follows directly from $\gamma_{g_{1}g_{2}}=\gamma_{g_{2}}\gamma_{g_{1}}$. Note that $\alpha_{g}$ is not a group action.

Obviously, by (4.7) the $G$-equivariance of $F$ in terms of the physical symmetries $\tilde{a}_{g}$ induces a $G$-equivariance in Fourier space, in term of the action $\gamma_{g}$.

Lemma 4.5**.**

Let $g\in G$. Then $F(\gamma_{g}W)=\gamma_{g}F(W)$ for all $W=(\Omega,\omega)\in\mathcal{X}$, where again we extend the action notation: $\gamma_{g}(\Omega,\omega)=(\Omega,\gamma_{g}\omega)$.

Remark 4.6**.**

In terms of the notation introduced in (4.1), the expressions for $\beta_{g}$ and $\alpha_{g}$ are, with ${j}=(n,m)=((\tilde{n},n_{4}),m)$,

[TABLE]

and

[TABLE]

Explicit formulas for the symmetries under consideration in this work are given in Section 5.2.

Two straightforward properties that will be useful later are

[TABLE]

Remark 4.7**.**

The real-valuedness of $u$ provides another symmetry, which in Fourier space corresponds to invariance under the transformation $\gamma_{*}$ introduced in Definition 2.10. One way to deal with this symmetry is to consider the extended symmetry group $\widehat{G}$ generated by $G$ and $\gamma_{*}$. However, one would have to interpret $\alpha_{*}({j})$ as the complex conjugation operator on $S^{1}\subset\mathbb{C}$, which is not a multiplication operator, and it is also not differentiable, which would cause problems later on. One solution is to split $\omega$ into real and imaginary parts, but that is not the way we proceed here. Instead, we stick with the symmetry group $G$ and consider the action of $\gamma_{*}$ separately. We note that $\gamma_{*}$ commutes with $\gamma_{g}$:

[TABLE]

which follows either from straightforward symmetry consideration in physical space, or directly in term of $\alpha_{g}$ and $\beta_{g}$ from (4.12).

4.4 Reduction to symmetry variables

The set of Fourier coefficients that are $G$-invariant is given by

[TABLE]

Again we have a correspondence with the physical space, by (4.7) elements of $X^{\text{{sym}}}$ correspond to functions that are $G$-invariant via the action $\tilde{a}_{g}$.

Lemma 4.8**.**

Let $W=(\Omega,\omega)\in\mathbb{C}\times X^{\text{{sym}}}$. Then $\omega(x,t)$ given by (4.5) satisfies $[\tilde{a}_{g}\omega](x,t)=\omega(x,t)$ for all $g\in G$.

In this section we study the properties of $X^{\text{{sym}}}$. Some of the arguments in this section are essentially the same as the ones in [43, Section 3.2]. Nevertheless we repeat the main arguments and definitions here, since the setting and notation are slightly different. The short proofs of some of the lemmas are also transcribed here to keep the current paper self-contained.

We first list some properties of $\alpha$ that will be useful in what follows.

Lemma 4.9**.**

[43*, Lemma 3.5]** Let ${j}\in J$ and $g,\widetilde{g}\in G$.

     (a) $\alpha_{g}({j})\,\alpha_{g^{-1}}(\beta_{g}({j}))=1$.

     (b) If $\beta_{g}({j})=\beta_{\widetilde{g}}({j})$ and $\alpha_{\widetilde{g}^{-1}g}({j})=1$, then $\alpha_{g}({j})=\alpha_{\widetilde{g}}({j})$.

     (c) If $\beta_{g}({j})={j}$, then $\alpha_{\widetilde{g}g\widetilde{g}^{-1}}(\beta_{\widetilde{g}}({j}))=\alpha_{g}({j})$.*

Proof.

For part (a) we use (4.10) to infer that

[TABLE]

For part (b) we write $g=\widetilde{g}\widetilde{g}^{-1}g$. From the first assumption and (4.9) if follows that $\beta_{\widetilde{g}^{-1}g}({j})=j$. By applying (4.10) to $g_{1}=\widetilde{g}$ and $g_{2}=\widetilde{g}^{-1}g$, we see that the second assumption implies that

[TABLE]

For part (c) we write ${j}^{\prime}=\beta_{\widetilde{g}}({j})$ and apply (4.10) twice to obtain

[TABLE]

where the final equality follows from part (a). ∎

The symmetry group $G$ may be exploited to reduce the number of independent variables of elements in $X^{\text{{sym}}}$. From now on we use the notation

[TABLE]

For any ${j}\in J$ we define the stabilizer

[TABLE]

and the orbit

[TABLE]

Remark 4.10** (Orbit-stabilizer).**

The orbit-stabilizer theorem implies that $|G_{{j}^{\prime}}|=|G_{j}|$ for all ${j}^{\prime}\in G.{j}$, where $|\cdot|$ denotes the cardinality. Indeed, stabilizers of different elements in an orbit are related by conjugacy, and $|G|=|G_{j}|\cdot|G.{j}|$. More generally, when $q$ is a function from $J$, or from subset thereof that is $G$-invariant under the action $\beta_{g}$, to some linear space, then we have

[TABLE]

Lemma 4.11**.**

[43*, Lemma 3.8]** Let ${j}\in J$ be arbitrary. We have the following dichotomy:

     (a) either $\alpha_{g}({j})=1$ for all $g\in G_{{j}}$;

     (b) or $\sum_{g\in G_{{j}}}\alpha_{g}({j})=0$.*

Proof.

Fix ${j}\in J$. We see from (4.10) that $\alpha_{g_{1}g_{2}}({j})=\alpha_{g_{1}}({j})\alpha_{g_{2}}({j})$ for all $g_{1},g_{2}\in G_{j}$. Hence we can interpret $\alpha_{g}({j})$ as a group action of the stabilizer subgroup $G_{j}$, acting by multiplication on the unit circle $S^{1}=\{z\in\mathbb{C}:|z|=1\}$. We consider the stabilizer of $1\in S^{1}$:

[TABLE]

and its orbit

[TABLE]

By the orbit-stabilizer theorem (we use Remark 4.10, but now for the $G_{j}$-action $\alpha_{g}({j})$)

[TABLE]

The set $O_{1}\subset S^{1}$ is invariant under multiplication and division. In particular, if $|O_{1}|=N\in\mathbb{N}$, then $O_{1}=\{e^{2\pi in/N}:n=0,1,\dots,N-1\}$. If $N=1$ we have $H_{1}=G_{j}$ and alternative (a) follows, whereas if $N>1$ then we see that $\sum_{z\in O_{1}}z=\sum_{n=0}^{N-1}e^{2\pi in/N}=0$, hence we conclude from (4.15) that alternative (b) holds. ∎

The indices for which alternative (b) in Lemma 4.11 applies are denoted by

[TABLE]

It follows from the next lemma, which is a slight generalization of Lemma 4.11, that the set $J^{\text{{triv}}}$ is invariant under $G$.

Lemma 4.12**.**

[43*, Lemma 3.10]** Let ${j}\in J$ be arbitrary. We have the following dichotomy:

     (a) either $\alpha_{g}({j}^{\prime})=1$ for all $g\in G_{{j}^{\prime}}$ and all ${j}^{\prime}\in G.{j}$;

     (b) or $\sum_{g\in G_{{j}^{\prime}}}\alpha_{g}({j}^{\prime})=0$ for all ${j}^{\prime}\in G.{j}$.*

Proof.

Let ${j}\in J$ and ${j}^{\prime}\in G.{j}$. Let $\widetilde{g}\in G$ be such that $\widetilde{g}.{j}={j}^{\prime}$. Then a conjugacy between $G_{j}$ and $G_{{j}^{\prime}}$ is given by $g\to\widetilde{g}g\widetilde{g}^{-1}$. It follows from Lemma 4.9(c) that

[TABLE]

The assertion now follows from Lemma 4.11. ∎

When considering $\omega\in X^{\text{{sym}}}$, the indices ${j}$ in $J^{\text{{triv}}}$ are the ones for which $\omega_{j}$ necessarily vanishes.

Lemma 4.13**.**

[43, Lemma 3.11]** Let $\omega\in X^{\text{{sym}}}$. Then $\omega_{j}=0$ for all ${j}\in J^{\text{{triv}}}$.

Proof.

Fix ${j}\in J$. For any $\omega\in X^{\text{{sym}}}$ we have in particular $[\gamma_{g}\omega]_{j}=\omega_{j}$ for all $g\in G_{j}$. By summing over $g\in G_{j}$ and using that $g.{j}={j}$ for $g\in G_{j}$, we obtain

[TABLE]

If ${j}\in J^{\text{{triv}}}$ then the right hand side vanishes, hence $\omega_{j}=0$. ∎

Lemma 4.13 implies that

[TABLE]

In other words, we may restrict attention to the Fourier coefficients corresponding to indices in the complement

[TABLE]

Remark 4.14**.**

On $S^{1}$ we have $z^{-1}=z^{*}$. In particular, Lemma 4.11 implies that

[TABLE]

This identity is used, through Lemma 4.17, to extract both the set $J^{\text{{triv}}}$ and the values $|G_{j}|$ for ${j}\in J^{\text{{sym}}}$ from what is computed in the code, see Remark 4.19.

For $\omega\in X^{\text{{sym}}}$ the coefficients $\omega_{j}$ with ${j}\in J^{\text{{sym}}}$ are not all independent. To take advantage of this, we choose a fundamental domain $J^{\text{{dom}}}$ of $G$ in $J$, i.e., $J^{\text{{dom}}}$ contains precisely one element of each group orbit. The arguments below are independent of which fundamental domain one chooses. We now define the set of symmetry reduced indices as

[TABLE]

and the space of symmetry reduced variables as

[TABLE]

In slight abuse of notation we will also interpret $\mathrm{e}_{j}$ with ${j}\in J^{\text{{red}}}$ as elements of $X^{\text{{red}}}$. We may then interpret $X^{\text{{red}}}$ as a subspace of $X$ by writing $\phi=\sum_{{j}\in J^{\text{{red}}}}\phi_{j}\mathrm{e}_{j}$ for $\phi\in X^{\text{{red}}}$.

The action of $\gamma_{g}$ on basis vectors is

[TABLE]

hence

[TABLE]

Before we specify the dependency of the coefficients $\{\omega_{j}\}_{{j}\in J^{\text{{sym}}}}$ on the symmetry reduced variables $\{\omega_{j}\}_{{j}\in J^{\text{{red}}}}$ for $\omega\in X^{\text{{sym}}}$, we derive some additional properties of $\alpha_{g}({j})$ for ${j}\in J^{\text{{sym}}}$.

Lemma 4.15**.**

[43, Lemma 3.13]** Let $g_{1},g_{2}\in G$ and ${j}\in J^{\text{{sym}}}$. If $g_{1}.{j}=g_{2}.{j}$ then $\alpha_{g_{1}}({j})=\alpha_{g_{2}}({j})$.

Proof.

Since $g_{2}^{-1}g_{1}.{j}={j}$ and ${j}\in J^{\text{{sym}}}$ we have $\alpha_{g_{2}^{-1}g_{1}}({j})=1$ by Lemma 4.11(a). An application of Lemma 4.9(b) concludes the proof. ∎

Definition 4.16**.**

Let ${j}$ by any element of $J^{\text{{sym}}}$ and ${j}^{\prime}$ any element in its orbit $G.{j}$. We can choose an $\widetilde{g}=\widetilde{g}({j},{j}^{\prime})\in G$ such that $\widetilde{g}.{j}={j}^{\prime}$. For such ${j}$ and ${j}^{\prime}$ we define

[TABLE]

This is independent of the choice of $\widetilde{g}$ by Lemma 4.15, and

[TABLE]

This $\widetilde{\alpha}$ will be used to define the natural map $\Sigma$ from $X^{\text{{red}}}$ to $X^{\text{{sym}}}$, defined below in (4.22). It also appears in the description of the action of $\sum_{g\in G}\gamma_{g}$ on unit elements.

Lemma 4.17**.**

We have

[TABLE]

Proof.

It follows from (4.17) and Lemma 4.9(a) that

[TABLE]

Based on the orbit-stabilizer theorem, we use the notation $\widetilde{g}({j},{j}^{\prime})$ from Definition 4.16 to write the right-hand side of (4.20) as

[TABLE]

where we have used (4.10) in the second equality. By combining (4.20) and (4.21) with (4.16) and (4.18) the assertion follows. ∎

We define the projection $\Pi:X\to X^{\text{{red}}}$ by

[TABLE]

Clearly $\Pi^{2}=\Pi$. Next, the group average $\mathcal{A}:X\to X^{\text{{sym}}}$ is

[TABLE]

We note that $\mathcal{A}^{2}=\mathcal{A}$. To relate elements of $X^{\text{{red}}}$ to elements of $X^{\text{{sym}}}$ it turns out that it is useful to introduce the rescaling $\mathcal{R}:X\to X$

[TABLE]

Finally, using $\widetilde{\alpha}$ introduced in Definition 4.16, we define the linear map $\Sigma:X^{\text{{red}}}\to X$ by

[TABLE]

Lemma 4.18**.**

**

*     (a) For $\phi\in X^{\text{{red}}}$ we have $\Sigma\phi=\mathcal{A}\mathcal{R}\phi$, hence $\Sigma\phi\in X^{\text{{sym}}}$.

     (b) $\Pi\Sigma$ is the identity on $X^{\text{{red}}}$.

     (c) $\Sigma\Pi$ is the identity on $X^{\text{{sym}}}$.

     (d) $\Sigma$ is bijective as a map from $X^{\text{{red}}}$ to $X^{\text{{sym}}}$ with inverse $\Pi:X^{\text{{sym}}}\to X^{\text{{red}}}$.*

Proof.

We start with part (a). Let $\phi\in X^{\text{{red}}}$. Then it follows from the definitions of $\mathcal{A}$ and $\mathcal{R}$, as well as Lemma 4.17 that

[TABLE]

Part (b) is an immediate consequence of (4.22) and $\widetilde{\alpha}({j},{j})=1$.

To prove part (c), let $\omega\in X^{\text{{sym}}}$ be arbitrary and set $\widetilde{\omega}=\omega-\Sigma\Pi\omega$. Then from part (b) we get

[TABLE]

hence $\widetilde{\omega}_{{j}}=0$ for all ${j}\in J^{\text{{red}}}$. Consider any fixed ${j}^{\prime}\in J^{\text{{sym}}}$, then there exists a ${j}\in J^{\text{{red}}}\cap G.{j}^{\prime}$. Let $\widetilde{g}\in G$ be such that $\widetilde{g}.{j}={j}^{\prime}$. Since $\widetilde{\omega}\in X^{\text{{sym}}}$ we have

[TABLE]

Since $\alpha_{\widetilde{g}}({j})\in S^{1}$ this implies that $\widetilde{\omega}_{{j}^{\prime}}=0$, and since ${j}^{\prime}\in J^{\text{{sym}}}$ was arbitrary, we conclude that $\widetilde{\omega}_{{j}^{\prime}}=0$ for all ${j}^{\prime}\in J^{\text{{sym}}}$. It then follows from Lemma 4.13 that

[TABLE]

Hence, since $\omega\in X^{\text{{sym}}}$ was arbitrary, $\Sigma\Pi$ is the identity on $X^{\text{{sym}}}$.

Part (d) follows immediately from (b) and (c). ∎

We will need to determine $\widetilde{\alpha}$, particularly for computing the derivative of $\mathcal{F}$, see (4.40). The next remark explains how this is accomplished in the code.

Remark 4.19**.**

In the code we compute on a finite set of indices $J^{\dagger}=\{j=(n,m):n\in\mathcal{E}^{\dagger}\}.$ In particular we determine

[TABLE]

Since each orbit contains precisely one element of $J^{\text{{dom}}}$, it follows from Lemma 4.17 that

[TABLE]

Since $\widetilde{\alpha}(j,j^{\prime})\in S^{1}\subset\mathbb{C}$, it is straightforward to infer the index set $J^{\dagger}\cap J^{\text{{sym}}}$ from S. Furthermore, for $j^{\prime}\in J^{\dagger}\cap J^{\text{{sym}}}$ the value of $\texttt{S}_{j^{\prime}}\in\mathbb{C}$ gives us direct access to both the order of the stabilizer $|G_{j^{\prime}}|$ and the value $\widetilde{\alpha}(j,j^{\prime})$ for the unique $j\in J^{\text{{red}}}\cap G.j^{\prime}$.

We finish this subsection by extending the complex conjugate symmetry $\gamma_{*}$ to the reduced space. We define a corresponding action $\widetilde{\gamma}_{*}$ on $X^{\text{{red}}}$ by

[TABLE]

Clearly $\gamma_{*}$ leaves $X^{\text{{sym}}}$ invariant (see (4.13)). It then follows from Lemma 4.18(d) that

[TABLE]

As one may expect, $\widetilde{\gamma}_{*}$ is an involution. In particular, this is useful in order to symmetrize $\phi$ by $\frac{1}{2}(\phi+\widetilde{\gamma}_{*}\phi)$, which is invariant under $\widetilde{\gamma}_{*}$.

Lemma 4.20**.**

We have $\widetilde{\gamma}_{*}^{2}=\text{{Id}}_{X^{\text{{red}}}}$.

Proof.

Let $\phi\in X^{\text{{red}}}$. Since $\gamma_{*}$ leaves $X^{\text{{sym}}}$ invariant it follows that

[TABLE]

where in the third equality we have used 4.18(c), and in the final equality we have used 4.18(b). ∎

4.5 Functional analytic setup in symmetry reduced variables

We now introduce the framework (with symmetry-reduced variables), where we define the fixed point operator that is actually used in practice to validate the numerical solution.

First, we introduce the space

[TABLE]

Elements in $\mathcal{X}^{\text{{red}}}$ are denoted by

[TABLE]

and the operators $\Pi$ and $\Sigma$ naturally extend to $\mathcal{X}$ and $\mathcal{X}^{\text{{red}}}$ respectively via

[TABLE]

Similarly, we extend the definition of $\widetilde{\gamma}_{*}$ to $\mathcal{X}^{\text{{red}}}$ by

[TABLE]

For later use, we also define the space $\mathcal{X}^{\text{{sym}}}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\mathbb{C}\times X^{\text{{sym}}}$.

Introducing the space $\mathcal{X}^{\text{{red}}}_{-2,-1}$, which is obtained from $\mathcal{X}_{-2,-1}$ the same way $\mathcal{X}^{\text{{red}}}$ is obtained from $\mathcal{X}$, we define the map $\mathcal{F}^{\text{{red}}}:\mathcal{X}^{\text{{red}}}\to\mathcal{X}^{\text{{red}}}_{-2,-1}$ by

[TABLE]

where

[TABLE]

and the phase condition in the symmetrized setting is given by

[TABLE]

where $j=(\tilde{n},n_{4},m)$, and $\hat{\phi}$ is some fixed element of $X^{\text{{red}}}$. Notice that, if we choose $\hat{\omega}=\Sigma\hat{\phi}$ in (2.13), then $F^{\text{{red}}}_{\leftmoon}(\phi)=F_{\leftmoon}(\Sigma\phi)$ and therefore

[TABLE]

Recalling (4.6), we also introduce the weights

[TABLE]

and (see Remark 2.7)

[TABLE]

together with the norm

[TABLE]

with index set

[TABLE]

We recall that $|G.{j}|=|G|/|G_{j}|$ by Remark 4.10. This observation is combined with Remark 4.14 to determine the weights $\xi^{\text{{s}}}_{j}$ in the code. The next lemma shows that the norm $\|\cdot\|_{\mathcal{X}^{\text{{red}}}}$ is compatible with the symmetrization, as well as complex conjugation.

Lemma 4.21**.**

**

*     (a) For all $\varphi\in\mathcal{X}^{\text{{red}}}$ we have $\|\varphi\|_{\mathcal{X}^{\text{{red}}}}=\|\Sigma\varphi\|_{\mathcal{X}}$.

     (b) For all $\varphi\in\mathcal{X}^{\text{{red}}}$ we have $\|\widetilde{\gamma}_{*}\varphi\|_{\mathcal{X}^{\text{{red}}}}=\|\varphi\|_{\mathcal{X}^{\text{{red}}}}$.*

Proof.

We first observe that $\xi_{{j}}$ is invariant under $G$. Namely, when we denote $(n^{\prime},m^{\prime})=g.(n,m)$, then it follows from (4.11) that $|n^{\prime}|_{1}=|n|_{1}$ for any $g\in G$, since $C_{g}$ is a signed permutation matrix.

From the definition (4.22), we then obtain,

[TABLE]

since $\widetilde{\alpha}\in S^{1}$. This proves part (a).

To prove part (b), we first note that, for all $W\in\mathcal{X}$, $\|\gamma_{*}W\|_{\mathcal{X}}=\|W\|_{\mathcal{X}}$ by Definition 2.10 and the observation $\xi_{{j}}$ is invariant under $\gamma_{*}$, in the sense that $|{-n}|_{1}=|n|_{1}$. It then follows from part (a) and (4.23) that

[TABLE]

To solve the zero finding problem $\mathcal{F}^{\text{{red}}}=0$ on $\mathcal{X}^{\text{{red}}}$, we analyse a fixed point operator as in Section 2.4. The role of $\mathcal{X}$ is taken over by $\mathcal{X}^{\text{{red}}}$, and the norm $\|\cdot\|_{\mathcal{X}}$ is replaced by the weighted norm $\|\cdot\|_{\mathcal{X}^{\text{{red}}}}$, which is symmetry compatible in the sense of Lemma 4.21.

Notice that $\mathcal{F}^{\text{{red}}}$ inherits the equivariance property of $\mathcal{F}$ with respect to the complex conjugacy $\gamma_{*}$.

Lemma 4.22**.**

Assume that $\hat{\omega}$ used in (2.13) is such that, $\hat{\omega}=\Sigma\hat{\phi}$ for some $\hat{\phi}\in X^{\text{{red}}}$ such that $\widetilde{\gamma}_{*}\hat{\phi}=\hat{\phi}$. Then, for all $\varphi\in\mathcal{X}^{\text{{red}}}$,

[TABLE]

Proof.

Using successively (4.26), (4.23) and Lemma 2.11, we obtain

[TABLE]

Having chosen $N^{\dagger}$ and $\tilde{N}$, we define the index sets

[TABLE]

The size of the Galerkin projection is the number of elements in $\mathcal{E}^{\dagger}_{\text{red}}(N^{\dagger})$, which is substantially smaller than the number of elements in $\mathcal{E}^{\dagger}(N^{\dagger})$, since we restrict to symmetry reduced variables. Indeed, the number of independent variables is reduced by roughly a factor $|G|$.

The construction of linear operators $A_{\text{{red}}}$ and $\widehat{A}_{\text{{red}}}$ is completely analogous to Section 2.3. To incorporate the incompressibility, we introduce

[TABLE]

Theorem 4.23**.**

Let $\eta>1$. Assume there exist non-negative constants $Y_{0}^{\text{{red}}}$, $Z_{0}^{\text{{red}}}$, $Z_{1}^{\text{{red}}}$ and $Z_{2}^{\text{{red}}}$ such that

[TABLE]

Assume also that

• •

the forcing term $f$ is time independent, $G$-invariant, and has average zero,

• •

$\bar{\varphi}=(\bar{\Omega},\bar{\phi})$* is in $\in\mathcal{X}^{\text{{red}}}_{\text{{div}}}$,*

• •

$\widetilde{\gamma}_{*}\bar{\varphi}=\bar{\varphi}$* and $\hat{\phi}$ (used to define the phase condition (4.25)) is such that $\widetilde{\gamma}_{*}\hat{\phi}=\hat{\phi}$.*

If

[TABLE]

then, for all $r\in[r_{\min},r_{\max})$ there exists a unique $\tilde{\varphi}=(\tilde{\Omega},\tilde{\phi})\in\mathcal{B}_{\mathcal{X}^{\text{{red}}}}(\bar{\varphi},r)$ such that $\mathcal{F}^{\text{{red}}}(\tilde{\varphi})=0$, where $\mathcal{B}_{\mathcal{X}^{\text{{red}}}}(\bar{\varphi},r)$ is the closed ball in ${\mathcal{X}^{\text{{red}}}}$, centered at $\bar{\varphi}$ and of radius $r$, and

[TABLE]

Besides, this unique $\tilde{\varphi}$ also lies in $\in\mathcal{X}^{\text{{red}}}_{\text{{div}}}$. Finally, defining $u=M\Sigma\tilde{\phi}$, there exists a pressure function $p$ such that $(u,p)$ is a $\frac{2\pi}{\tilde{\Omega}}$-periodic, real valued, analytic and $G$-invariant solution of Navier-Stokes equations (1.1).

Proof.

The structure of the proof is the same as the one of Theorem 2.15. Therefore, we only outline the main steps. First, we consider the operator

[TABLE]

and use (4.32) to infer that it is a contraction on $\mathcal{B}_{\mathcal{X}^{\text{{red}}}}(\bar{\varphi},r)$ for all $r\in[r_{\min},r_{\max})$, yielding the existence of a unique zero $\tilde{\varphi}$ of $\mathcal{F}^{\text{{red}}}$ in $\mathcal{B}_{\mathcal{X}^{\text{{red}}}}(\bar{\varphi},r)$. Next, we use that $\bar{\varphi}\in\mathcal{X}^{\text{{red}}}_{\text{{div}}}$ and that $T_{\text{{red}}}\left(\mathcal{X}^{\text{{red}}}_{\text{{div}}}\right)\subset\mathcal{X}^{\text{{red}}}_{\text{{div}}}$, which follows from arguments entirely analogous to the ones in the proof of Theorem 2.15, to infer that $\tilde{\varphi}$ also belongs to $\mathcal{X}^{\text{{red}}}_{\text{{div}}}$. Then, we use Lemmas 4.21 and 4.22, once again as in the proof of Theorem 2.15, to conclude that $\widetilde{\gamma}_{*}\tilde{\varphi}=\tilde{\varphi}$.

Defining $\tilde{W}=(\tilde{\Omega},\tilde{\omega})=\Sigma\tilde{\varphi}$, we then have $\mathcal{F}(\tilde{W})=0$ (with $\hat{\omega}=\Sigma\hat{\phi}$ in (2.13)), $\tilde{\omega}\in\mathcal{X}_{\text{{div}}}$ and $\gamma_{*}\tilde{W}=\tilde{W}$. By Lemma 2.5, there exists $p$ such that $(u,p)$ is a $2\pi/\tilde{\Omega}$-periodic, real valued, analytic solution of Navier-Stokes equations (1.1).

Finally, by construction $\tilde{\omega}\in X^{\text{{sym}}}$, and by Lemma 4.8 $\tilde{a}_{g}\tilde{\omega}=\tilde{\omega}$ for all $g\in G$. From (4.4) we infer that $u=M\tilde{\omega}$ satisfies $a_{g}u=u$, and that $p=\Gamma\Phi(u)$ satisfies $a_{g}p=p$ (see Remark 4.3), i.e., $(u,p)$ is $G$-invariant. ∎

Remark 4.24**.**

To make sure that $\bar{\varphi}\in\mathcal{X}^{\text{{red}}}_{\text{{div}}}$ in practice we replace a candidate $\bar{\phi}$ by $\Pi\mathcal{A}\Pi_{\text{{div}}}\Sigma\bar{\phi}$, where $\Pi_{\text{{div}}}$ is a projection of $X$ onto

[TABLE]

which is $G$-invariant and hence invariant under the averaging operator $\mathcal{A}$.

Similarly, to make sure that $\widetilde{\gamma}_{*}\bar{\varphi}=\bar{\varphi}$ we make sure $\bar{\Omega}$ is real-valued and we replace a candidate $\bar{\phi}$ by $\frac{1}{2}(\bar{\phi}+\Pi\gamma_{*}\Sigma\bar{\phi})$, and analogously $\widetilde{\gamma}_{*}\hat{\phi}=\hat{\phi}$. Since $\mathcal{X}^{\text{{red}}}_{\text{{div}}}$ is invariant under $\widetilde{\gamma}_{*}$, this new $(\bar{\Omega},\frac{1}{2}(\bar{\phi}+\Pi\gamma_{*}\Sigma\bar{\phi}))$ is still in $\mathcal{X}^{\text{{red}}}_{\text{{div}}}$ if $\bar{\varphi}$ is.

The resulting $\bar{\varphi}$ and $\hat{\phi}$ may then be represented by intervals in the code, but that is no impediment whatsoever.

In section 4.6 we explain how explicit expressions for the bounds $Y_{0}^{\text{{red}}}$, $Z_{0}^{\text{{red}}}$, $Z_{1}^{\text{{red}}}$ and $Z_{2}^{\text{{red}}}$ can be obtained from the bounds $Y_{0}$, $Z_{0}$, $Z_{1}$ and $Z_{2}$ of Section 3.

4.6 Bounds in the symmetric setting

As in Section 3, we assume throughout that the approximate solution $\bar{\varphi}=(\bar{\Omega},\bar{\phi})$ only has a finite number of non-zero modes, i.e., $\bar{\phi}_{n}=0$ for all $n\notin\mathcal{S}^{\text{{sol}}}$ for some finite set $\mathcal{S}^{\text{{sol}}}$. Similarly, $\hat{\phi}_{n}=0$ for all $n\notin\mathcal{S}^{\text{{sol}}}$.

4.6.1 bound $Y_{0}^{\text{{red}}}$

There are essentially no changes compared to Section 3.2 in the computation of the bound on the residue, except that we need to take into account the symmetry respecting norm:

[TABLE]

4.6.2 bound $Z_{0}^{\text{{red}}}$

The bound $Z_{0}^{\text{{red}}}$ is completely analogous to the one in Section 3.3. We just need to use the operator norm

[TABLE]

applied to $B_{\text{{red}}}=I-A_{\text{{red}}}\widehat{A}_{\text{{red}}}$, which again has finitely many nonzero components only.

4.6.3 bound $Z_{1}^{\text{{red}}}$

We introduce

[TABLE]

The finite part of the $Z_{1}^{\text{{red}}}$ estimate is analogous to section 3.4:

[TABLE]

For the tail part, we will take advantage of the estimates in section 3.4. We need to estimate

[TABLE]

We introduce, cf. (4.24),

[TABLE]

where $\Psi$ is defined in (2.11). Furthermore, we set $\lambda_{j}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\lambda_{n}$ for $j=(n,m)$, see (2.16). It follows from the definition of $\mathcal{F}^{\text{{red}}}$ and $\widehat{A}_{\text{{red}}}$ that

[TABLE]

In particular, analogous to (3.8),

[TABLE]

Therefore, the tail estimate reduces to

[TABLE]

It follows from (4.38) and the formula (4.22) for the symmetrization $\Sigma$ that

[TABLE]

For any $j\in J^{\text{{red}}}\setminus\widetilde{\mathcal{S}}^{\text{sol}}_{\text{red}}(\tilde{N})$ we then use the triangle inequality to estimate

[TABLE]

Writing $\mathcal{C}_{j,j^{\prime}}\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}D\Psi_{{j}}(\Sigma\bar{\phi})e_{{j}^{\prime}}$ for convenience, we use Remark 4.10, Equation (4.27) and Lemma 4.25 below to obtain

[TABLE]

The right-hand side of (4.41) is exactly the one estimated in Section 3.4 with $\bar{\omega}=\Sigma\bar{\phi}$. Hence (4.39) is bounded by

[TABLE]

where the elements in the righthand side are given in (3.4) with $\bar{\omega}=\Sigma\bar{\phi}$.

Lemma 4.25**.**

Let $\omega\in X^{\text{{sym}}}$. Writing $D\Psi(\omega)=\left(\mathcal{C}_{j^{\prime},j}\right)_{j^{\prime},j\in J}$, we have

[TABLE]

Proof.

Analogous to Lemma 4.2, we observe that $\Psi$ is a $G$-equivariant, which implies that $\Psi(\gamma_{g}\omega)=\gamma_{g}\Psi(\omega)$ for all $\omega\in X$ and $g\in G$. For $\omega\in X^{\text{{sym}}}$ we then obtain

[TABLE]

Using (4.17) twice, we get for all $j^{\prime},j\in J$,

[TABLE]

which yields the assertion, since $|\alpha_{g}(j^{\prime})|=1$ for all $g\in G$ and all $j^{\prime}\in J$. ∎

Finally, we set

[TABLE]

4.6.4 bound $Z_{2}^{\text{{red}}}$

For any $\varphi,\varphi^{\prime}\in\mathcal{X}^{\text{{red}}}$ with $\|\varphi\|_{\mathcal{X}^{\text{{red}}}},\|\varphi^{\prime}\|_{\mathcal{X}^{\text{{red}}}}\leq 1$ we need to estimate

[TABLE]

where $\bar{W}=\Sigma\bar{\varphi}$, $W=\Sigma\varphi$, $W^{\prime}=\Sigma\varphi^{\prime}$, and by Lemma 4.21 $\|W\|_{\mathcal{X}},\|W^{\prime}\|_{\mathcal{X}}\leq 1$.

As in Section 3.5, we start by splitting

[TABLE]

But since $W,W^{\prime}\in\mathcal{X}^{\text{{sym}}}$ and $\mathcal{F}$ is $G$-equivariant, we have that $D^{2}\mathcal{F}(\bar{W})(W,W^{\prime})$ is $G$-invariant, i.e. belongs to $\mathcal{X}^{\text{{sym}}}$, and thus by Lemma 4.21

[TABLE]

Therefore, we can directly use the estimates of Section 3.5 to obtain

[TABLE]

where

[TABLE]

with

[TABLE]

5 Application to Taylor-Green flow

In this section, we present the results obtained with our computer-assisted approach for the Navier-Stokes equations (1.1) with forcing (1.2). In Section 5.1, we set the stage by recalling some know analytic results in this context, complemented by numerical simulations. In Section 5.2, we then show how the estimates of Section 4.6, combined with Theorem 4.23, can be used to obtain rigorous existence results and error bounds about periodic solutions.

5.1 Analytical and numerical background

We consider the forced Navier-Stokes equations for an incompressible, homogeneous fluid

[TABLE]

on a cubical domain of dimension $L$ with periodic boundary conditions. Here $\rho_{0}$ is the density of the fluid and $\nu_{0}$ is its viscosity, while the forcing is chosen to be of the simple form

[TABLE]

where $\gamma_{0}$ parametrizes its intensity. This corresponds to the planar case of a family of flows introduced by Taylor and Green [37] to study the interaction of motions on different spatial scales. Non-dimensionalization (without introducing new notation for the new variables) leads to

[TABLE]

on a cube of dimension $2\pi$, where $\nu=\sqrt{\frac{32\pi^{3}}{\gamma_{0}L^{3}}}\,\nu_{0}$ is a dimensionless parameter and the dimensionless forcing is

[TABLE]

To be able to compare with the literature we use the geometric Reynolds number $Re=\frac{\sqrt{8\pi}}{\nu}$ in the discussion of the bifurcation diagram below, cf. Figure 3.

The Navier-Stokes equations (5.1) under the forcing (5.2) admit an equilibrium solution for which we have the analytic expression

[TABLE]

We will refer to this solution as the viscous equilibrium.

The associated forcing for the vorticity equation is given by

[TABLE]

The viscous equilibrium consists of four counter-rotating vortices. Its vertical vorticity is shown in a plane of constant height in Figure 2. It is straightforward to verify that this solution is invariant under the following symmetry operations:

• •

Translation over any distance $d$ in the vertical direction, $T_{d}$.

• •

Reflection in the $x_{1}$-direction, $S_{x_{1}}$.

• •

Reflection in the $x_{2}$-direction, $S_{x_{2}}$.

• •

Reflection in the $x_{3}$-direction, $S_{x_{3}}$.

• •

Rotation about the axis $x_{1}=x_{2}=0$ over $\pi/2$ followed by a shift over $L/2$ in the $x_{1}$-direction, $R$.

• •

A shift over $L/2$ in both the $x_{1}$ and $x_{2}$ directions, $D$.

These operations generate the group of spatial symmetries of the Navier-Stokes equation with planar Taylor-Green forcing. In addition, the system is equivariant under translations in time. One may reduce this continuous symmetry to a discrete subgroup by considering, for a periodic orbit of period $\tau$, the shift by $\frac{\tau}{k}$ in time, $P_{k}$, for some $k\in\mathbb{N}$ ($k=4$ for the solutions studied below).

The linear instability of the viscous equilibrium at high Reynolds number (that is at low viscosity) has been investigated at length in the literature, for instance by Sipp and Jacquin [33]. In those studies, the emphasis is on the rapid formation of small-scale structures. In the current context, we are interested in the first instabilities that occur when increasing the Reynolds number from zero. Due to a classical result by Serrin [32] the viscous equilibrium is guaranteed to be the unique limit state of the flow for any viscosity greater than

[TABLE]

Below that value, linear instabilities occur, giving rise to branches of solutions for which some of the symmetries are broken. Such bifurcating branches can be numerically approximated by standard methods [30]. Using these methods, we found that the first instability is a Hopf bifurcation that gives rise to a branch of two-dimensional periodic solutions. Subsequently, this branch appears to turn unstable at a point where at least one family of three-dimensional periodic solutions branches off. A partial bifurcation diagram is shown in Figure 3. In order to differentiate between the solutions, we compute the deviation from reflection symmetry by computing the maximum of \|u-S_{x_{1}}u\|\raisebox{-2.15277pt}{\scalebox{0.5}{E}}, where

[TABLE]

is the energy. We normalize this measure by the maximum it can attain, and display

[TABLE]

on the vertical axis of the bifurcation diagram.

5.2 New results for 2D periodic orbits

Using a computer program in MATLAB, we computed a numerical approximation $\bar{W}$ of a periodic orbit and applied Theorem 4.23, together with the bounds of Section 4.6, to validate this solution with explicit error bounds, see Theorems 5.1 and 5.2 below. In the appendix, we describe how to recover errors bounds for the associated velocity $u$ and pressure $p$ that solve the Navier-Stokes equations.

This validation was done for two separate orbits, represented by $\mathrm{p}_{1}$ and $\mathrm{p}_{2}$ in the bifurcation diagram of Figure 3. These two orbits are trivially invariant under $T_{d}$ and $S_{x_{3}}$. In addition, they are invariant under a symmetry group $G$ of order $16$, generated by the following three symmetries:

[TABLE]

The associated actions $a_{g_{i}}$ (see (4.1)) can be represented by the following matrices and vectors:

[TABLE]

These group actions can then be represented at the level of the Fourier coefficients of the vorticity, via $\alpha_{g}$ and $\beta_{g}$ (see (4.7)-(4.8) and Remark 4.6) given by

[TABLE]

The symmetries are used to define the subspace $X^{\text{{red}}}$ (as in Section 4.4), which is then used in Theorem 4.23. In practice, these symmetries are of course first found by inspecting the numerical solution. However, once a true solution is shown to exist via Theorem 4.23, we automatically also get a proof that this solution satisfies these symmetries.

The value of the other parameters ($\eta$, $N^{\dagger}$, $\tilde{N}$, etc) that are needed to defined $A_{\text{{red}}}$ or used to compute the estimates of Section 4.6, are listed in Table 1. We choose a rectangular box of Fourier coefficients to compute the numerical approximate solutions:

[TABLE]

See Remark 5.5 below for a brief discussion about how the parameter values were chosen.

Theorem 5.1**.**

Let $\nu=0.286$ and let $f^{\omega}$ be given by (5.4). Let $G$ be the symmetry group generated by $g_{1}$, $g_{2}$, and $g_{3}$ given in (5.6). Let $\bar{\omega}$ and $\bar{\Omega}$ be the Fourier data time frequency provided in the file dataorbit1.mat (which can be downloaded at [41]). Let $r_{\text{{sol}}}^{\omega}=2.6314\cdot 10^{-5}$. Then there exists a unique pair $(\tilde{\Omega},\tilde{\omega})$ satisfying the error estimate

[TABLE]

such that $\tilde{\omega}$ is an analytic, divergence free, $G$-invariant, $2\pi/\tilde{\Omega}$-periodic solution of the vorticity equation (2.7).

Proof.

Let $\bar{\phi}=\Pi\bar{\omega}$ and $\bar{\varphi}=(\bar{\Omega},\bar{\phi})$. Let $\hat{\phi}=\frac{1}{2}(\bar{\phi}+\Pi\gamma_{*}\Sigma\bar{\phi})$ be used to define $F^{\text{{red}}}_{\leftmoon}$.

We start by using Remark 4.24 to replace $\bar{\phi}$ by a divergence free variant $\bar{\phi}_{0}\in X^{\text{{red}}}_{\text{{div}}}$ such that $\widetilde{\gamma}_{*}\bar{\phi}_{0}=\bar{\phi}_{0}$. This new $\bar{\phi}_{0}$ still has a finite number of nonzero modes: $(\bar{\phi}_{0})_{n}=0$ for $n\notin\mathcal{S}^{\text{{sol}}}$. We take $\eta=1$ and perform an interval arithmetic computation to determine bounds

[TABLE]

This bounds absorbs rounding errors introduced in various symmetrization steps.

Let $F^{\text{{red}}}$ and $\mathcal{X}^{\text{{red}}}$ (with $\eta=1$) be as defined in Section 4.5. We now use $\bar{\varphi}_{0}=(\bar{\Omega},\bar{\phi_{0}})$ instead of $\bar{\varphi}$ as the numerical approximate solution and apply the 2D estimates to obtain the bounds $Y_{0}^{\text{{red}}}$, $Z_{0}^{\text{{red}}}$, $Z_{1}^{\text{{red}}}$ and $Z_{2}^{\text{{red}}}$ defined in Section 4.6. Using the script runprooforbit1, available at [41], with $\eta=1$ and the other parameters chosen as in Table 1, we checked that the inequalities (4.32) are satisfied. In order to make this verification rigorous, all the bounds are evaluated using the interval arithmetic package INTLAB [28].

Recalling from (4.6) that $\xi_{j}=\eta^{|n|_{1}}$, one has that the symmetric weights $\xi^{\text{{s}}}_{j}=\xi^{\text{{s}}}_{j}(\eta)$ given in (4.27) are continuous in the parameter $\eta$. Moreover, $\xi^{\text{{s}}}_{j}(\eta)$ is increasing in $\eta$. Hence, the bounds $Y_{0}^{\text{{red}}}$, $Z_{0}^{\text{{red}}}$, $Z_{1}^{\text{{red}}}$ and $Z_{2}^{\text{{red}}}$ (given respectively in (4.35), (4.36), (4.42) and (4.44)) are also continuous functions of $\eta$. By continuity of the bounds in $\eta$, since the strict inequalities in (4.32) are satisfied for $\eta=1$, they still hold for all $\eta\in[1,\tilde{\eta}]$ for some $\tilde{\eta}>1$. In particular, we can use Theorem 4.23 for any $\eta\in(1,\tilde{\eta}]$, which yields the existence of a unique zero $\tilde{\varphi}=(\tilde{\Omega},\tilde{\phi})$ of $\mathcal{F}^{\text{{red}}}$, so that, by monotonicity of $\xi^{\text{{s}}}_{j}(\eta)$,

[TABLE]

for each $r\in[r_{\min}(\eta),r_{\max}(\eta))$, with $r_{\min}(\eta)$ and $r_{\max}(\eta)$ given in (4.33) and (4.34), respectively. By continuity we find

[TABLE]

where $r_{\min}=r_{\min}(1)$, and have that the solution $\tilde{\varphi}$ is unique in the ball of radius $r$ around $\bar{\varphi_{0}}$ (with $\eta=1$) for any $r<r_{\max}=r_{\max}(1)$.

To translate this existence and uniqueness result to the vorticity, let $\tilde{\omega}=\Sigma\tilde{\phi}\in(\ell^{1}_{\tilde{\eta}})^{3}$. The asserted properties of $\tilde{\omega}$ follow from the arguments given towards the end of the proof of Theorem 4.23, and by applying the triangle inequality after checking that $r_{\min}+r_{0}^{\omega}<r_{\text{{sol}}}^{\omega}<r_{\max}-r_{0}^{\omega}$. ∎

The solution in Theorem 5.1 represents the point $\mathrm{p}_{1}$ in Figure 3. Four snap shots of the third component (the nontrivial one for a 2D solution) of the vorticity field $\bar{\omega}$ from Theorem 5.1 are shown in Figure 4. In exactly the same way, but with the script runprooforbit2 we prove

Theorem 5.2**.**

Let $\nu=0.265$ and $r_{\text{{sol}}}^{\omega}=2.2491\cdot 10^{-6}$. Let $\bar{\omega}$ and $\bar{\Omega}$ be the Fourier data time frequency provided in the file dataorbit2.mat. Then there exists a unique pair $(\tilde{\Omega},\tilde{\omega})$ satisfying the error estimate (5.8), such that $\tilde{\omega}$ is an analytic, divergence free, $G$-invariant, $2\pi/\tilde{\Omega}$-periodic solution of the vorticity equation (2.7).

Remark 5.3**.**

As a direct corollary of Theorem 5.2, we get a validated velocity $u$ and pressure $p$ that solves the Navier-Stokes equations via Lemma 2.5, with explicit error bounds presented in the appendix (Lemma 6.2). In particular, Theorem 1.1 holds.

The solution from Theorem 5.2, corresponding to the point $\mathrm{p}_{2}$ in Figure 3, is depicted in Figure 1. The script mimicproofsfloats may be used to perform faster but non-rigorous versions of the computations using floating point calculations rather than interval arithmetic. The computational cost for Theorem 5.2 is considerably higher than for Theorem 5.1, as more Fourier modes need to be considered to make the proof work.

Remark 5.4**.**

Since the size of the estimates deteriorates for larger values of $\eta$, we choose $\eta=1$. By continuity of the bounds in $\eta$, as explained in the proof of Theorem 5.1, this is enough to prove that we obtain an analytic solution, although we do not get an explicit decay rate of the Fourier coefficients. If obtaining such decay rate was needed, for instance to obtain error estimates on derivatives of the solution, or to obtain an explicit domain of analyticity, this could be achieved by using some explicit $\eta>1$ to compute the bounds $Y_{0}^{\text{{red}}}$, $Z_{0}^{\text{{red}}}$, $Z_{1}^{\text{{red}}}$ and $Z_{2}^{\text{{red}}}$.

Remark 5.5**.**

We briefly discuss the choices made for the other computational parameters. The $Z_{2}^{\text{{red}}}$ bound is relatively insensitive to the computational parameters (as it roughly measures the size of the second derivative at the numerical approximation). We then consider a set $\mathcal{S}^{\text{{sol}}}$ of rectangular shape (as explained in Table 1), with sufficiently many modes $N^{\text{{sol}}}_{x_{1}}=N^{\text{{sol}}}_{x_{2}}$ and $N^{\text{{sol}}}_{t}$ (for the 2D solution $N^{\text{{sol}}}_{x_{3}}=0$) to make the bound $Y_{0}^{\text{{red}}}$ on the residual so small that $2Y_{0}^{\text{{red}}}Z_{2}^{\text{{red}}}$ is roughly of size $10^{-2}$. The choice of this threshold size is based on two criteria. Primarily, we choose it small, so that at the next step (see below) $Z_{1}^{\text{{red}}}$ is allowed to be very close to $1$, as decreasing $Z_{1}^{\text{{red}}}$ is computationally very costly. On the other hand, the threshold should not be too small, as that would require many modes in $\mathcal{S}^{\text{{sol}}}$ to decrease the residual sufficiently.

Subsequently, one may determine from a relatively cheap calculation how large $\tilde{N}$ needs to be for the tail bound to satisfy $1-(Z_{1}^{\text{{red}}})^{\text{{tail}}}>(2Y_{0}^{\text{{red}}}Z_{2}^{\text{{red}}})^{1/2}$. We note that $Z_{0}^{\text{{red}}}$ is negligible in practice. Finally, it requires some experimentation to determine how large one needs to choose $N^{\dagger}$ in order to also satisfy $1-(Z_{1}^{\text{{red}}})^{\text{{finite}}}>(2Y_{0}^{\text{{red}}}Z_{2}^{\text{{red}}})^{1/2}$.

6 Appendix

The estimates obtained in Section 3 and Section 4.6 are used as input for Theorem 4.23, which allows us to validate symmetric periodic solutions $\omega$ of the vorticity equation, with explicit error bounds, as illustrated in Section 5. In this appendix, we describe how to recover errors bounds for the associated velocity $u$ and pressure $p$ that solve the Navier-Stokes equations.

We start with a variation, adapted to our framework, of the classical result stating that a curl-free vector field can be written as a gradient, which was used already in the proof of Lemma 2.5.

Lemma 6.1**.**

Let $\Phi\in\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}}$ satisfy

[TABLE]

Then the map $\Gamma:\left(\mathbb{C}^{3}\right)^{\mathbb{Z}^{4}}\to\mathbb{C}^{\mathbb{Z}^{4}}$ constructed component-wise as

[TABLE]

is well defined, and $p=\Gamma\Phi$ satisfies $\Phi=-\nabla p$.

Proof.

To ensure that $\Gamma$ is well defined, it suffices to show that, for all $n\in\mathbb{Z}^{4}$ and $l,m\in\{1,2,3\}$

[TABLE]

Indeed, since $\nabla\times\Phi=0$, we have that for all $n\in\mathbb{Z}^{4}$ and all $l,m\in\{1,2,3\}$,

[TABLE]

which immediately yields (6.1). Therefore $p=\Gamma\Phi$ is well defined, and we are left to check that $\Phi=-\nabla p$. If $\tilde{n}=0$ then we have

[TABLE]

because we assumed $\Phi_{n}=0$ for all $\tilde{n}=0$. If $\tilde{n}\neq 0$, for any $l\in\{1,2,3\}$ we distinguish between two cases. If $n_{l}\neq 0$, then

[TABLE]

If $n_{l}=0$, then $-\left(\nabla p\right)_{n}^{(l)}=0$, but there exists an $m\neq l$ such that $n_{m}\neq 0$ and thus by (6.2) we find $\Phi_{n}^{(l)}=0$, i.e. $-\left(\nabla p\right)_{n}^{(l)}=\Phi_{n}^{(l)}$ also holds. ∎

The above lemma can be used in the context of Navier-Stokes equations, to recover the pressure from the velocity (we recall that the velocity itself is recovered from the vorticity via $u=M\omega$). We point out that an alternative (arguably more classical) approach is to define $p$ as the solution of the Poisson equation

[TABLE]

satisfying

[TABLE]

Indeed, the latter approach is going to be useful in the sequel, as (6.3) allows to recover sharper error bounds for the pressure (compared to using Lemma 6.1 only).

Our aim is to derive error estimates for the velocity and the pressure, that can be applied as soon as we have validated a divergence-free solution $W$ of the vorticity equation via Theorem 2.15 or Theorem 4.23.

Lemma 6.2**.**

Assume that for some $\bar{W}=(\bar{\Omega},\bar{\omega})\in\mathcal{X}^{\text{{div}}}$, $\eta>1$, we have proven the existence of $r>0$ and of $W=(\Omega,\omega)\in\mathcal{B}_{\mathcal{X}^{\text{{div}}}}(\bar{W},r)$ such that $\mathcal{F}(W)=0$. Define

[TABLE]

where $\Phi$ is defined in (2.12) and $\Gamma$ is defined as in Lemma 6.1. We also consider the sequence $\bar{p}\in\mathbb{C}^{\mathbb{Z}^{4}}$ defined as

[TABLE]

Then, we have the following error estimates for the velocity and the pressure:

[TABLE]

Remark 6.3**.**

As explained in Remark 2.17, these weighted $\ell^{1}$-norms control the $\mathcal{C}^{0}$-norms of the errors (explicitly). Notice also that, even though we used $\eta=1$ to validate the vorticity in Theorem 5.1 and Theorem 5.2, by continuity of the estimates with respect to $\eta$ we get for free a validation for some $\tilde{\eta}>1$ (see the proof of Theorem 5.1), and thus Lemma 6.2 is directly applicable.

Proof.

The error estimate for the velocity simply follows from the definition of $M$:

[TABLE]

To obtain the error estimate for the pressure, we use the fact that $(u,p)$ are smooth solutions of Navier-Stokes equations (see Lemma 2.5), and thus (6.3) holds. In Fourier space, this reduces to

[TABLE]

Since both $u$ and $\bar{u}$ are divergence-free, by Lemma 2.2 we can write, for all $\tilde{n}\neq 0$,

[TABLE]

We then estimate

[TABLE]

where $\left|\cdot\right|$ applied to a sequence must be understood component-wise, i.e. $|u^{(l)}|$ is the sequence whose $n$-th element is equal to $|u^{(l)}_{n}|$. Finally, we obtain

[TABLE]