A Hopf bifurcation in the planar Navier-Stokes equations

We consider the Navier-Stokes equation for an incompressible viscous fluid on a square, satisfying Navier boundary conditions and being subjected to a time-independent force. As the kinematic viscosity is varied, a branch of stationary solutions is shown to undergo a Hopf bifurcation, where a periodic cycle branches from the stationary solution. Our proof is constructive and uses computer-assisted estimates.


Introduction and main result
We consider the Navier-Stokes equations ∂ t u − ν∆u + (u · ∇)u + ∇p = f , ∇ · u = 0 on Ω , (1.1) for the velocity u = u(t, x, y) of an incompressible fluid on a planar domain Ω, satisfying suitable boundary conditions for (x, y) ∈ ∂Ω and initial conditions at t = 0. Here, p denotes the pressure, and f = f (x, y) is a fixed time-independent external force. Our focus is on solution curves and bifurcations as the kinematic velocity ν is being varied. In order to reduce the complexity of the problem, the domain Ω is chosen to be as simple as possible, namely the square Ω = (0, π) 2 . Following [21], we impose Navier boundary conditions on ∂Ω, which are given by u 1 = ∂ x u 2 = 0 on {0, π} × (0, π) , u 2 = ∂ y u 1 = 0 on (0, π) × {0, π} . (

1.2)
A fair amount is known about the (non)uniqueness of stationary solutions in this case [21]. This includes the existence of a bifurcation between curves of stationary solutions with different symmetries.
Here we prove the existence of a Hopf bifurcation for the equation (1.1) with boundary conditions (1.2), and with a forcing function f that satisfies (∂ x f 2 − ∂ y f 1 )(x, y) = 5 sin(x) sin(2y) − 13 sin(3x) sin(2y) . (1.3) In a Hopf bifurcation, a stationary solution loses stability and a small-amplitude limit cycle branches from the stationary solution [1,3,4]. Among other things, this introduces a time scale in the system and increases its complexity. In this capacity, Hopf bifurcations in the Navier-Stokes equation constitute an important first step in the transition to turbulence in fluids, as was described in the seminal work [5]. Numerically, there is plenty of evidence that Hopf bifurcations occur in the Navier-Stokes equation, but proofs are still very scarce. An explicit example of a Hopf bifurcation was given in [6] for the rotating Bénard problem. A proof exists also for the Couette-Taylor problem [7,9]. Sufficient conditions for the existence of a Hopf bifurcation in a Navier-Stokes setting are presented in [13].
Before giving a precise statement of our result, let us replace the vector field u in the equation (1.1) by ν −1 u. The equation for the rescaled function u is α∂ t u − ∆u + γ(u · ∇)u + ∇p = f , ∇ · u = 0 on Ω , (1.4) where γ = ν −2 . The value of α that corresponds to (1.1) is ν −1 , but this can be changed to any positive value by rescaling time.
Numerically, it is possible to find stationary solutions of (1.4) for a wide range of values of the parameter γ. At a value γ 0 ≈ 83.1733117 . . . we observe a Hopf bifurcation that leads to a branch of periodic solutions for γ > γ 0 .
For a fixed value of α, the time period τ of the solution varies with γ. Instead of looking for τ -periodic solution of (1.4) for fixed α, we look for 2π-periodic solutions, where α = 2π/τ has to be determined. To simplify notation, a 2π-periodic function will be identified with a function on the circle T = R/(2πZ). Our main theorem is the following. Theorem 1.1. There exists a real number γ 0 = 83.1733117 . . ., an open interval I including γ 0 , and a real analytic function (γ, x, y) → u γ (x, y) from I × Ω to R 2 , such that u γ is a stationary solution of (1.4) and (1.2) for each γ ∈ I. In addition, there exists a real number α 0 = 4.66592275 . . ., an open interval J centered at the origin, two real analytic functions γ and α on J that satisfy γ(0) = γ 0 and α(0) = α 0 , respectively, as well as two real analytic functions (s, t, x, y) → u s,e (t, x, y) and (s, t, x, y) → u s,o (t, x, y) from J × T × Ω to R 2 , such that the following holds. For any given β ∈ C satisfying β 2 ∈ J, the vector field u = u s,e + βu s,o with s = β 2 is a solution of (1.4) and (1.2) with γ = γ(s) and α = α(s). Furthermore, u 0,e (t, . , .) = u γ 0 and ∂ t u 0,o (t, . , .) = 0.
To our knowledge, this is the first result establishing the existence of a Hopf bifurcation for the Navier-Stokes equation in a stationary environment.
Our proof of this theorem is computer-assisted. The solutions are obtained by rewriting (1.4) and (1.2) as a suitable fixed point equation for scalar vorticity of u. Here we take advantage of the fact that the domain is two-dimensional. We isolate the periodic branch from the stationary branch by using a scaling that admits two distinct limits at the bifurcation point. This approach is also known as the blow-up method, which is a common tool in the study of singularities and bifurcations [8].
Computer-assisted methods have been applied successfully to many different problems in analysis, mostly in the areas of dynamical systems and partial differential equations. Here we will just mention work that concerns the Navier-Stokes equation or Hopf bifurcations. For the Navier-Stokes equation, the existence of symmetry-breaking bifurcations among stationary solutions has been established in [10,21]. Periodic solutions for the Navier-Stokes flow in a stationary environment have been obtained in [20]. In the case of periodic forcing, the problem of existence and stability of periodic orbits has been investigated in [14]. Concerning the existence of Hopf bifurcations, a computer-assisted proof was given recently in [22] for a finite-dimensional dynamical system; and an extension of their method to the Kuramoto-Sivashinsky PDE is presented in [23]. For other recent computer-assisted proofs we refer to [16,17,18,19] and references therein. As mentioned earlier, a system similar to the one considered here is known to exhibit a symmetry-breaking bifurcation within the class of stationary solutions [21]. The broken symmetry is y → π/2 − y. Based on a numerical computation of eigenvalues, we expect an analogous bifurcation to occur here at γ ≈ 1450. Interestingly, the Hopf bifurcation described here occurs at a significantly smaller value of γ. We have not tried to prove the existence of a symmetry-breaking bifurcation for the forcing (1.3), since such an analysis would duplicate the work in [21] and go beyond the scope of the present paper.
The remaining part of this paper is organized as follows. In Section 2, we first rewrite (1.4) as an equation for the function Φ = ∂ y u 1 − ∂ x u 2 , which is the scalar vorticity of −u. After a suitable scaling Φ = U β φ, the problem of constructing the solution branches described in Theorem 1.1 is reduced to three fixed point problems for the function φ. These fixed point equations are solved in Section 3, based on estimates described in Lemmas 3.3, 3.4, and 3.6. Section 4 is devoted to the proof of these estimates, which involves reducing them to a large number of trivial bounds that can be (and have been) verified with the aid of a computer [24].

Fixed point equations
The goal here is to rewrite the equation (1.4) with boundary conditions (1.2) as a fixed point problem. Applying the operator ∂ : (u 1 , u 2 ) → ∂ 2 u 1 − ∂ 1 u 2 on both sides of the equation (1.4), we obtain Here, we have used that ∂ (u · ∇)u = u · ∇Φ. Using the divergence-free condition ∇ · u = 0, one also finds that ∆u = J∇Φ , If Φ vanishes on the boundary of ∂Ω, then the equation (2.2) can be inverted to yield where ∆ denotes the Dirichlet Laplacean on Ω. In Section 3 we will define a space of real analytic functions Φ that admit a representation It is straightforward to check that the corresponding vector field u = (u 1 , u 1 ) satisfies the Navier boundary conditions (1.2). So a solution u of (1.4) and (1.2) can be obtained via (2.5) from a solution Φ of the equation (2.1). For convenience, we write (2.1) as where L is the symmetric bilinear form defined by The coefficients Φ j,k in the series (2.4) are 2π-periodic functions and thus admit an expansion Φ j,k = n∈Z Φ n,j,k cosi n , cosi n (t) = cos(nt) if n ≥ 0, sin(−nt) if n < 0. (2.8) Denote by N 0 the set of all nonnegative integers. For any subset N ⊂ N 0 we define where sin m (z) = sin(mz). In particular, the even frequency part Φ e (odd frequency part Φ o ) of Φ is defined to be the function E N Φ, where N is the set of all even (odd) nonnegative integers. This leads to the decomposition Φ = Φ e + Φ o that will be used below.
To simplify the discussion, consider first non-stationary periodic solutions. For γ near the bifurcation point γ 0 , we expect Φ to be nearly time-independent. So in particular, Φ o is close to zero. Consider the function φ = φ e + φ o obtained by setting φ e = Φ e and The scaling factor β = 0 will be chosen below, in such a way that φ e and φ o are of comparable size. Substituting into (2.6) yields the equation Finally, we convert (2.11) to a fixed point equation by applying the inverse of α∂ t − ∆ to both sides.
One of the features of the equation (2.11) is that the time-translate of a solution is again a solution. We eliminate this symmetry by imposing the condition φ 1,1,1 = 0. In addition, we choose β = θ −1 Φ −1,1,1 , where θ is some fixed constant that will be specified later. This leads to the normalization conditions (2.14) Notice that β enters our main equationφ = φ only via its square s = β 2 . It is convenient to regard s to be the independent parameter and express γ as a function of s. The functions γ = γ(s) and α = α(s) are determined by the condition thatφ satisfies the normalization conditions (2.14). Applying the functionals A and B to both sides of (2.11), using the identities A∆ = −2A, A∂ t = B, B∆ = −2B, B∂ t = −A, and imposing the conditions Aφ = 0 and Bφ = θ, we find that For a fixed value of s, define F s (φ) =φ, whereφ is given by (2.13), with γ = γ(s, φ) and α = α(s, φ) determined by (2.15). The fixed point equation for F s is used to find non-stationary time-periodic solutions of (2.11).
Remark 1. The choice (2.15) guarantees that Aφ = 0 and Bφ = θ, even if φ does not satisfy the normalization conditions (2.14). Thus, the domain of the map F s can include non-normalized function φ. (The same is true for the map F γ described below.) But a fixed point of this map will be normalized by construction.
In order to determine the bifurcation point γ 0 and the corresponding frequency α 0 , we consider the map F : φ →φ given by (2.13) with s = 0. The values of γ and α are again given by (2.15), so that Aφ = 0 and Bφ = θ. We will show that this map F has a fixed point φ with the property that φ n,j,k = 0 whenever |n| > 1. The values of γ and α for this fixed point define γ 0 and α 0 .
A similar map F γ : φ →φ, given by (2.13) with s = 0, is used to find stationary solutions of the equation (2.6). In this case, the value of γ is being fixed, and φ o is taken to be zero. The goal is to show that this map F γ has a fixed point φ γ that is independent of time t. Then Φ = φ γ is a stationary solution of (2.6).
We finish this section by computing the derivative of the map F s described after (2.15). The resulting expressions will be needed later. Like some of the above, the following is purely formal. A proper formulation will be given in the next section. For simplicity, assume that φ depends on a parameter. The derivative of a quantity q with respect to this parameter will be denoted byq. Define The above expressions forγ andα are obtained by differentiating (2.15).

The associated contractions
In this section, we formulate the fixed point problems for the maps F , F γ , and F s in a suitable functional setting. The goal is to reduce the problems to a point where we can invoke the contraction mapping theorem. After describing the necessary estimates, we give a proof of Theorem 1.1 based on these estimates. We start by defining suitable function spaces. Given a real number ρ > 1, denote by A the space of all functions h ∈ L 2 (T) that have a finite norm h , where Here cosi n are the trigonometric function defined in (2.8). It is straightforward to check that A is a Banach algebra under the pointwise product of functions. That is, gh ≤ g h for any two functions g, h ∈ A. We also identify functions on T with 2π-periodic functions on R. In this sense, a function in A extends analytically to the strip T (ρ) = {z ∈ C : | Im z| < log ρ}.
Given in addition ̺ > 1, we denote by B the space of all function Φ : T 2 → A that admit a representation (2.4) and have a finite norm x, y)) in this space will also be identified with a function (t, x, y) → Φ(t, x, y) on T 3 , or with a function on R 3 that is 2π-periodic in each argument.
In this sense, every function in B extends analytically to T (ρ) × T (̺) 2 .
We consider A and B to be Banach spaces over F ∈ {R, C}. In the case F = R, the functions in these spaces are assumed to take real values for real arguments.
Clearly, a function Φ ∈ B admits an expansion (2.9) with N = N 0 . The sequence of Fourier coefficients Φ n,k,j converges to zero exponentially as |n| + j + k tends to infinity. If all but finitely many of these coefficients vanish, then Φ is called a Fourier polynomial. The equation (2.9) with N ⊂ N 0 non-empty defines a continuous projection E N on B whose operator norm is 1. Using Fourier series, it is straightforward to see that the equation (2.16) defines two bounded linear operators L α and L ′ α on B, for every α ∈ C. The operator L α is in fact compact. Specific estimates will be given in Section 4. The following will be proved in Section 4 as well.
Proposition 3.1. If Φ and φ belong to B, then so does |∆| −1/2 L(Φ)φ, and This estimate implies e.g. that the transformation φ →φ, given by (2.13) for fixed values of s, γ and α, is well-defined and compact as a map from B to B.
As is common in computer-assisted proofs, we reformulate the fixed point equation for the map φ →φ as a fixed point problem for an associated quasi-Newton map. Since we need three distinct versions of this map, let us first describe a more general setting. The domain of N is defined to be the set of of all h ∈ B with the property that ϕ+Lh ∈ D.
Notice that, if h is a fixed point of N , then ϕ+Lh is a fixed point of F . In our applications, ϕ is an approximate fixed point of F and L is an approximate inverse of I − DF (ϕ).
The following is an immediate consequence of the contraction mapping theorem.
where ε, K are positive real numbers that satisfy ε + Kδ < δ. Then F has a fixed point in ϕ + LB δ . If L is invertible, then this fixed point is unique in ϕ + LB δ .
In our applications below, B is always a subspace of B. The domain parameter ρ and the constant θ that appears in the normalization condition (2.14) are chosen to have the fixed values ρ = 2 5 , θ = 2 −12 . (3.6) The domain parameter ̺ is defined implicitly in our proofs. That is, the lemmas below hold for ̺ > 1 sufficiently close to 1.
Consider first the problem of determining the bifurcation point γ 0 and the associated Let s = 0, and denote by D the set of all functions φ ∈ B with the property that Bφ = 0. Define F : D → B to be the map φ →φ given by (2.13), with γ = γ(φ) and α = α(φ) defined by the equation (2.15). Clearly, F is not only C 1 but real analytic on D. Our proof of this lemma is computer-assisted and will be described in Section 4. By Proposition 3.2, the map F has a unique fixed point φ * ∈ ϕ + L 1 B δ . We define γ 0 = γ(φ * ) and α 0 = α(φ * ).
Our next goal is to construct a branch of periodic solutions for the equation ( Our proof of this lemma is computer-assisted and will be described in Section 4. As a consequence we have the following. Assume that some function ψ ∈ B satisfies DF 0 (φ * )ψ = ψ. We may assume that ψ takes real values for real arguments. A straightforward computation shows that DN 0 (0)L −1 ψ = L −1 ψ. Since DN 0 (0) is a contraction in the real setting, by Lemma 3.4, this implies that ψ = 0. So the operator DF 0 (φ * ) does not have an eigenvalue 1. This operator is compact, since it is the composition of a bounded linear operator with the compact operator L α . Thus, DF 0 (φ * ) has no spectrum at 1. By the implicit function theorem, there exists a complex open ball J , centered at the origin, such that the fixed point equation F s (φ) = φ has a solution φ = φ s for all s ∈ J . Furthermore, the curve s → φ s is analytic, passes through φ * at s = 0, and there is a unique curve with this property. By uniqueness, we also have φs = φ s for all s ∈ J , so φ s is real for real values of s ∈ J . Lemma 3.6. Let F = R. There exists an isomorphism L 0 of B such that the following holds. If N γ 0 denotes the the quasi-Newton map associated with (B, F γ 0 , φ * e , L 0 ), then the derivative DN γ 0 (0) of N γ 0 at the origin is a contraction.
Our proof of this lemma is computer-assisted and will be described in Section 4. As a consequence we have the following.
There exists an open disk I ⊂ C, centered at γ 0 , and an analytic curve γ → φ γ on I with values in B, such that F γ (φ γ ) = φ γ for all γ ∈ I. If γ belongs to the real interval I ∩ R, then φ γ is real. Furthermore, φ γ 0 = φ * e . The proof of this corollary is analogous to the proof of Corollary 3.5. We note that the disk I ∋ γ 0 is disjoint from the disk J ∋ 0 described in Corollary 3.5. So there is no ambiguity in using the notation γ → φ γ and s → φ s for the curve of stationary and periodic solutions, respectively, of the equation (2.11), Based on the results stated in this section, we can now give a Proof of Theorem 1.1. As described in the preceding sections, the curve γ → φ γ for γ ∈ I yields a curve γ → u γ of stationary solutions of the equation (1.4), where u γ = ∂ −1 φ γ . By our choice of function spaces, the function (γ, x, y) → u γ (x, y) is real analytic on I × T 2 , where I = I ∩ R.
Similarly, the curve s → φ s for s ∈ J defines a family of of non-stationary periodic solutions for (1.4), with γ = γ s and α = α s determined via the equation (2.15). To be more precise, the even frequency part φ s,e of φ s determines a vector field u s,e = ∂ −1 φ s,e , and the odd frequency part φ s,o determines a vector field u s,o = ∂ −1 φ s,o . If β is a complex number such that s = β 2 ∈ J , then u = u s,e + βu s,o is a periodic solution of (1.4), with γ = γ s and α = α s . Here, we have used the decomposition (2.10). By our choice of function spaces, the functions (s, t, x, y) → u s,e (t, x, y) and (s, t, x, y) → u s,o (t, x, y) are real analytic on J × T 3 , where J = J ∩ R. Clearly, ∂ t u 0,o (t, . , .) = 0, due to the normalization condition φ −1,1,1 = θ imposed in (2.14). And by construction, we have u = u γ 0 for s = 0.

Remaining estimates
What remains to be proved are Lemmas 3.3, 3.4, and 3.6. Our method used in the proof of Lemma 3.3 can be considered perturbation theory about the approximate fixed point ϕ of F . The function ϕ is a Fourier polynomial with over 20000 nonzero coefficients, so a large number of estimates are involved.
We start by describing bounds on the bilinear function L and on the linear operators L α and L ′ α . These are the basic building blocks for our transformations F , F s , and F γ . The "mechanical" part of these estimates will be described in Subsection 4.4.

The bilinear form L and a proof of Proposition 3.1
Consider the bilinear form L defined by (2.7). Using the identity (2.3), we have (4.1) In order to obtain accurate estimates, it is useful to have explicit expressions for L(Φ)φ in terms of the Fourier coefficients of Φ and φ. Given that L is bilinear, and that the identity (4.1) holds pointwise in t, it suffices to compute L(Φ)φ for the time-independent monomials Φ = sin J × sin K , φ = sin j × sin k , (4.2) with J, K, j, k > 0. A straightforward computation shows that with Θ as defined below. As a result we have (4.5) Proof of Proposition 3.1. Using the Cauchy-Schwarz inequality in R 2 , we find that Since the absolute value of N σ,τ is invariant under an exchange of (j, k) and (J, K), this implies that |N σ,τ | ≤ 1/4 where a ∨ b = max(a, b) for a, b ∈ R. As a result, we obtain the bound Using the nature of the norm (3.2), and the fact that A is a Banach algebra for the pointwise product of functions, this bound extends by bilinearity to arbitrary functions Φ, φ ∈ B.

QED
We note that the bound (4.8) exploits the cancellations that lead to the expression (4.3). A more straightforward estimate loses a factor of 2 with respect to (4.8). But it is not just this factor of 2 that counts for us. The expressions (4.5) for the coefficients N σ,τ and the bounds (4.7) are used in our computations and error estimates. The expression on the right hand side of (4.7) is a decreasing function of the wavenumbers j, k, J, K, so it can be used to estimate L(Φ)φ when Φ and/or φ are "tails" of Fourier series.

Estimating operator norms
Recall that a function φ ∈ B admits a Fourier expansion φ = n∈Z j,k∈N 1 φ n,j,k θ n,j,k , θ n,j,k def = cosi n × sin j × sin k , (4.14) and that the norm of φ is given by Let now n ≥ 0. A linear combination c + θ n,j,k + c − θ −n,j,k will be referred to as a mode with frequency n and wavenumbers (j, k) or as a mode of type (n, j, k). We assume of course that c − = 0 when n = 0. Since (4.15) is a weighted ℓ 1 norm, except for the ℓ 2 norm used for modes, we have a simple expression for the operator norm of a continuous linear where the third supremum is over all nonzero modes u of type (n, j, k). Let now n, j, k ≥ 1 be fixed. In computation where Lθ ±n,j,k is known explicitly, we use the following estimate. Denote by L n,j,k the restriction of L to the subspace spanned by the two functions θ ±n,j,k . For q ≥ 1 define L n,j,k q = sup 0≤p<q Lv p , v p = cos πp q θ n,j,k ρ n ̺ j+k + sin πp q θ −n,j,k ρ n ̺ j+k . Since every unit vector in the span of θ ±n,j,k lies within a distance less than π q of one of the vectors v p or its negative, we have L n,j,k ≤ L n,j,k q + π q L n,j,k . Thus Consider now the operator DF s (φ) described in (2.17), with φ ∈ E {0,1} B fixed. Iḟ φ = u n is a nonzero mode with frequency n ≥ 3, thenφ = 2|∆| −1/2 L 0 (φ)φ belongs to E N B with N = {n − 1, n, n + 1}. Thus, we haveγ =α = 0, and Due to the factor L α in this equation, if u n = c + θ n,j,k + c − θ −n,j,k with (j, k) and c ± fixed, then the ratios DF 0 (φ)u n / u n (4.20) are decreasing in n for n ≥ 3. And the limit as n → ∞ of this ratio is zero.
So for the operator L = DF 0 (φ), the supremum over n ∈ N 0 in (4.16) reduces to a maximum over finitely many terms. The same holds for the operator L = DN 0 (0) = DF 0 (φ * )L + I − L that is described in Lemma 3.4. This is a consequence of the following choice.
Remark 2. The operator L chosen in Lemma 3.4 is a "matrix perturbation" of the identity, in the sense that Lθ n,j,k = θ n,j,k for all but finitely many indices (n, j, k). The same is true for the operators L 1 and L 0 chosen in Lemma 3.3 and Lemma 3.6, respectively.

Computer estimates
Lemmas 3.3, 3.6, and 3.4 assert the existence of certain objects that satisfy a set of strict inequalities. The goal here is to construct these objects, and to verify the necessary inequalities by combining the estimates that have been described so far.
The above-mentioned "objects" are real numbers, real Fourier polynomials, and linear operators that are finite-rank perturbations of the identity. They are obtained via purely numerical computations. Verifying the necessary inequalities is largely an organizational task, once everything else has been set up properly. Roughly speaking, the procedure follows that of a well-designed numerical program, but instead of truncation Fourier series and ignoring rounding errors, we determine rigorous enclosures at every step along the computation. This part of the proof is written in the programming language Ada [25]. The following is meant to be a rough guide for the reader who wishes to check the correctness of our programs. The complete details can be found in [24].
An enclosure for a function φ ∈ B is a set in B that includes φ and is defined in terms of (bounds on) a Fourier polynomial and finitely many error terms. We define such sets hierarchically, by first defining enclosures for elements in simpler spaces. In this context, a "bound" on a map f : X → Y is a function F that assigns to a set X ⊂ X of a given type (Xtype) a set Y ⊂ Y of a given type (Ytype), in such a way that y = f (x) belongs to Y for all x ∈ X. In Ada, such a bound F can be implemented by defining a procedure F(X: in Xtype; Y: out Ytype).
Our most basic enclosures are specified by pairs S=(S.C,S.R), where S.C is a representable real number (Rep) and S.R a nonnegative representable real number (Radius). Given a Banach algebra X with unit 1, such a pair S defines a ball in X which we denote by S, X = {x ∈ X : x − (S.C)1 ≤ S.R}.
When X = R, then the data type described above is called Ball. Bounds on some standard functions involving the type Ball are defined in the package Flts Std Balls. Other basic functions are covered in the packages Vectors and Matrices. Bounds of this type have been used in many computer-assisted proofs; so we focus here on the more problem-specific aspects of our programs.
Consider now the space A for a fixed domain radius ̺ > 1 of type Radius. As mentioned before Remark 2, we only need to consider Fourier polynomials in A. Our enclosures for such polynomials are defined by an array(-I c .. I c ) of Ball. This data type is named NSPoly, and the enclosure associated with data P of this type is where ν is an increasing index function with the property that ν(−i) = −ν(i). The type NSPoly is defined in the package NSP, which also implements bounds on some basic operations for Fourier polynomials in A. Among the arguments to NSP is a nonnegative integer n (named NN). Our proof of Lemma 3.6 and Lemma 3.3 uses I c = n = 0 and I c = n = 1, respectively, and ν(i) = i. Values n ≥ 2 are uses when estimating the norm of Lu for the operator L = DN 0 (0), with u a mode of frequency n. In this case, ν takes values in {−n, n} or {−n − 1, −n, −n + 1, 0, n − 1, n, n + 1}, depending on whether n is odd or even. (The value ν = 0 is being used only for n = 2.) The package NSP also defines a data type NSErr as an array(0 .. I c ) of Radius. This type will be used below. This data type is named Fourier3, and the enclosure associated with F=(F.C,F.E) is

K=1
H J,K (F.E(J, K)) . Here, H J,K (E) denotes the set of all functions φ = where φ i can be any function in B whose coefficients φ i n,j,k vanish unless j ≥ J, k ≥ K, and |n| = ν(i). The type Fourier3 and bounds on some standard functions involving this type are defined in the child package NSP.Fouriers. This package is a modified version of the package Fouriers2 that was used earlier in [11,15,21]. The procedure Prod is now a bound on the bilinear map |∆| −1/2 L 0 . The error estimates used in Prod are based on the inequality (4.7). The package NSP.Fouriers also includes bounds InvLinear and DtInvLinear on the linear operators L α and L ′ α , respectively. These bounds use the estimates described in Subsection 4.3.
As far as the proof of Lemma 3.3 is concerned, it suffices now to compose existing bounds to obtain a bound on the map F and its derivative DF . This is done by the procedures GMap and DGMap in Hopf.Fix. Here we use enclosures of type NN=1.
The type of quasi-Newton map N defined by (3.4) has been used in several computerassisted proof before. So the process of constructing a bound on N from a bound on F has been automated in the generic packages Linear and Linear.Contr. (Changes compared to earlier versions are mentioned in the program text.) This includes the computation of an approximate inverse L 1 for the operator I − DF (ϕ). A bound on N is defined (in essence) by the procedure Linear.Contr.Contr, instantiated with Map => GMap. And a bound on DN is defined by Linear.Contr.Contr, with DMap => DGMap. Bounds on operator norms are obtained via Linear.OpNorm. Another problem-dependent ingredient in these procedures, besides Map and DMap, are data of type Modes. These data are constructed by the procedure Make in the package Hopf. They define a splitting of the given space B into a finite direct sum. For details on how such a splitting is defined and used we refer to [16].
If the parameter NN has the value 0, then the procedures GMap and DGMap define bounds on the map F γ and its derivative, respectively. The operator L 0 used in Lemma 3.6 has the property that M 0 = L 0 − I satisfies M 0 = P 0 M 0 P 0 , where P 0 = E {0} P m 0 for some positive integer m 0 . Here, and in what follows, P m denotes the canonical projection in B with the property that P m φ is obtained from φ by restricting the second sum in (4.14) to wavenumbers j, k ≤ m.
If NN has a value n ≥ 2, then the procedure DGMap defines a bound on the map (φ, ψ) → DF 0 (φ)ψ, restricted to the subspace E {0,1} B × E {n} B. The linear operator L that is used in Lemma 3.4 admits a decomposition L = I + M 1 + M 2 + . . . + M N of the following type. After choosing a suitable sequence n → m n of positive integers, we set M n = P n (L − I)P n , where P 1 = E {0,1} P m 1 and P n = E {n} P m n for n = 2, 3, . . . , N . This structure of L simplifies the use of (4.16) for estimating the norm of L = DN 0 (0). Furthermore, to check that L is invertible, it suffices to verify that I + M n is invertible on the finite-dimensional subspace P n B, for each positive n ≤ N .
The linear operator L 1 that is used in Lemma 3.3 is of the form L 1 = I + M 1 with M 1 as described above.
All the steps required in the proofs of Lemmas 3.3, 3.6, and 3.4 are organized in the main program Check. As n ranges from 0 to N = 305, this program defines the parameters that are used in the proof for NN = n, instantiates the necessary packages, computes the appropriate matrix M n , verifies that I + M n is invertible, reads ϕ from the file BP.approx, and then calls the procedure ContrFix from the (instantiated version of the) package Hopf.Fix to verify the necessary inequalities.
The representable numbers (Rep) used in our programs are standard [27] extended floating-point numbers (type LLFloat). High precision [28] floating-point numbers (type MPFloat) are used as well, but not in any essential way. Both types support controlled rounding. Radius is always a subtype of LLFloat. Our programs were run successfully on a 20-core workstation, using a public version of the gcc/gnat compiler [26]. For further details, including instruction on how to compile and run our programs, we refer to [24].