Uniqueness and Bifurcation Branches for Planar Steady Navier–Stokes Equations Under Navier Boundary Conditions

The stationary Navier–Stokes equations under Navier boundary conditions are considered in a square. The uniqueness of solutions is studied in dependence of the Reynolds number and of the strength of the external force. For some particular forcing, it is shown that uniqueness persists on some continuous branch of solutions, when these quantities become arbitrarily large. On the other hand, for a different forcing, a branch of symmetric solutions is shown to bifurcate, giving rise to a secondary branch of nonsymmetric solutions. This proof is computer-assisted, based on a local representation of branches as analytic arcs.


Introduction and Main Results
Let Ω ⊂ R 2 be a bounded domain and consider the stationary Navier-Stokes equations − νΔu + (u · ∇)u + ∇p = f , ∇ · u = 0 in Ω , that model the steady-state motion of an incompressible viscous fluid: u is its velocity, p its pressure, f is an external force, ν > 0 is the kinematic viscosity. Equation (1) need to be complemented with some boundary constraint, the most common being the no-slip boundary conditions that are physically reasonable if the boundary of Ω is solid and the flow is regular in a suitable sense. Existence and regularity results for (1) and (2) are classical topics [15], the latter being strongly influenced by the regularity of both the boundary ∂Ω and the source f . Much more delicate is the uniqueness, which is guaranteed only for large viscosities ν or small sources f . Through an application of the Sard-Smale Lemma, Foias-Temam [13,14] were able to prove that (1) and (2) admits a finite number of solutions, generically with respect to f and ν. Non-uniqueness has been obtained in very particular situations such as the Bénard problem [25], see also [20] for the same problem tackled through computer assistance, or the so-called Taylor problem, where one has multiplicity of solutions if f is large, see [32] and also [15, Theorem IX.2.2] for a slightly more general statement. There exists no multiplicity result valid in any situation, nor any detailed description of how the bifurcation from uniqueness might occur; see however [26]. Therefore, a complete comprehension of these phenomena is a challenging task, see [21,Problem 67].
In some situations, such as in 2D geophysical models [24], (2) is no longer suitable to describe the behavior of the fluid at the boundary and a slip boundary condition appears more realistic. In 1827, Navier [22] proposed boundary conditions with friction, with a stagnant layer of fluid close to the wall allowing a fluid where n is the outward normal vector to ∂Ω while τ is tangential. The boundary conditions (3) turn out to be appropriate in many physically relevant cases [29], in particular in presence of permeable walls [8] or of turbulent boundary layers [16,23]. The Navier-Stokes equations (1) under the Navier boundary conditions (3) (with and without friction) have been studied by many authors. The first contribution (in 1973) is due to Solonnikov-Scadilov [30], with external forces f ∈ L 2 (Ω). Concerning regularity results, we mention the works by Beirão da Veiga [9], Amrouche-Rejaiba [4], Acevedo et al. [2] while Clopeau et al. [12] and Iftimie-Sueur [18] studied the inviscid limit of (1) under conditions (3). Regularity results can also be found in the 3D work by Berselli [10] which appears relevant for our purposes, since he considers flat boundaries. Under the Navier boundary conditions (3) as well, uniqueness results for (1) are available only for small f or large ν, see e.g. Theorem 1.5 in [23], while multiplicity results are not known. In this paper we give results on this problem in the simple 2D case of a square: It is known [10] that for flat boundaries, the conditions (3) become of mixed Dirichlet-Neumann type and The mixed conditions are natural, since they guarantee that the boundary term vanishes after an integration by parts: The pressure p is defined up to an additive constant so that one can fix its mean value, for instance Ω p = 0 .
In Sect. 2 we take advantage of the geometry of the square Ω in order to obtain symmetry and regularity results for the solutions of (5). Our particular choice of domain allows us to obtain more precise information and to highlight some phenomena in a simple geometric situation. Still, it seems plausible that similar phenomena might occur for more general situations as well, including the case of no-slip boundary conditions (2). Then we focus our attention on uniqueness and multiplicity results for (5). In Sect. 3.2 we prove the following statement, which shows that small viscosities do not necessarily imply multiplicity of solutions for (5).

Theorem 1. There exists a continuous curve
For other curves ν → f ν , nonuniqueness and bifurcations may occur. Although we believe these to happen in general, we choose a particular forcing term. We take f symmetric with respect to a reflection about the line y = π/2, that is, One then expects that (at least) a solution of (5) has the same symmetry, namely We take f analytic and concentrated near the center of the square Ω. Then, as the viscosity decreases from ν = +∞, we are able to prove nonuniqueness via bifurcations and symmetry breaking. More precisely, we have Theorem 2. There exists an analytic function f , with f L 2 = 1, satisfying (7) and such that: • for all ν > 0, (5) and (6) admits an analytic solution (u ν , p ν ) satisfying (8), the curve ν → u ν is analytic and u ν is isolated within the subspace of symmetric functions satisfying (8); • there exists ν 1 > 0 such that the solution of (5) and (6) is unique for all ν > ν 1 ; • at some positive ν 0 ≤ ν 1 a pitchfork bifurcation occurs and the secondary branches of solutions that arise do not contain solutions satisfying (8).
The existence of a symmetric solution is standard, see Sect. 2, and it does not depend on the specific choice of f . The uniqueness statement for large ν follows from a general statement as well, see Proposition 4 where a lower bound for ν 1 is given. The remaining part of the proof of Theorem 2 is obtained through computer assistance and is described in Sect. 4, where we provide the explicit definition of the function f and the value ν 0 of the bifurcation point. Our computer assisted proof, which takes strong advantage from the Navier boundary conditions, is based on a well-established technique that, in the last few years, has been applied to many different kind of differential equations, see e.g. [6,7], and see [31] for a very recent result on the existence of a periodic solution for the full Navier-Stokes equations. We emphasize that the proof not only provides the existence of the solutions listed in Theorem 2, but also rigorous and tight bounds for their Fourier coefficients. Moreover, the proof of the pitchfork bifurcation and of the existence of the analytic branch of solutions is based on the Taylor expansion of the Fourier coefficients and is, to the best of our knowledge, completely new. We recall that the analyticity of the curve ν → u ν is known since Foias-Temam [14]. The secondary branch mentioned in Theorem 2 has been determined rigorously only close to the bifurcation value ν 0 , except for a few isolated points along a numerically determined continuation of the branch; see Figs. 1 and 2. For further details on our bifurcation results, both rigorous and numerical, we refer to Sect. 4.
Theorems 1 and 2 are complemented with a number of further results. In Sect. 2 we determine explicitly the eigenfunctions of the associated Stokes problem and we make clear how the symmetry property (7) for f influences the symmetry of some solutions of (5) and (6). In fact, the symmetric framework also yields a nonuniqueness criterion, see Corollary 1. Besides the expected uniqueness result for small f (Proposition 4), we prove a statement yielding uniqueness for some special f with arbitrarily large norm (Theorem 4), whose proof is based on the knowledge of explicit solutions of (1)-(3): we take advantage of the fact that the eigenfunctions of the Stokes problem are transformed into conservative vector fields by the nonlinearity. In fact, when f coincides with one of these eigenfunctions, we expect uniqueness of the solution of (1)-(3) for any value of ν: in Theorem 4 we obtain a "strange" necessary condition for uniqueness to fail, see (31). The remaining part of this paper is organized as follows. In Sect. 3 we prove Theorem 1 and provide some additional results concerning the uniqueness of the solution. In Sect. 4 we describe the functional setting of the computer assisted proof of Theorem 2. In Sect. 5 we provide the technical details of the computer assisted proof. Section 6 contains the conclusions and a list of open problems, some of which appear to be quite challenging and of wide interest from several points of view.

Existence and Symmetry of Solutions
Most of the results presented in this section are well-known, but it is useful to emphasise some peculiarities of the problem (5) in the square (4). For this reason, we sketch some steps in the proofs.
We first introduce the usual spaces of the Helmholtz-Weyl decomposition where v → n · v denotes the boundary-normal trace operator. It is known that where orthogonality is with respect to the scalar product in L 2 (Ω). Assuming that f ∈ L 2 (Ω), we say that u ∈ U is a weak solution of (5) if Weak solutions may also be found under the less restrictive assumption that f ∈ H −1 (Ω) but, since our main concern are smooth solutions, we assume here (at least) that f ∈ L 2 (Ω). In this case, as we shall see, the solution has the additional regularity u ∈ H 2 (Ω). Then, (11) implies that the first equation in (5), seen as an equality of functions in L 2 (Ω), is satisfied when projected onto G 1 . Due to (10), −νΔu + (u · ∇)u − f ∈ G 2 and, therefore, that there exists p ∈ H 1 (Ω) such that the full equation is satisfied in a strong form. This is the reason why the pressure p does not show up explicitly in the weak formulation (11), see also Lemma 1. The geometry of Ω will also be used to determine explicitly the eigenfunctions of −Δ in U , satisfying the boundary conditions in (5). Let Ψ 1,j,k (x, y) = sin(jx) cos(ky) , 0 , Ψ 2,j,k (x, y) = 0, cos(jx) sin(ky) .
Then the L 2 −normalized eigenfunctions of −Δ in U are given by In fact, the eigenfunctions of the Stokes operator are defined up to the addition of the gradient of an harmonic function, see [27, (1)-(2)-(3)]; here we take such function to be zero. Remark 1. The eigenfunctions in (13), that will also be used to find explicit solutions of (5) (see the proof of Theorem 4 below), may be obtained as follows. Consider the (scalar!) eigenfunctions ϕ of −Δ in The eigenfunctions with separated variables are given by ϕ j,k (x, y) = sin(jx) sin(ky) with corresponding eigenvalues λ j,k = j 2 + k 2 . The separated-variables-feature excludes, for instance, eigenfunctions that are linear combinations of sin(jx) sin(ky) and sin(kx) sin(jy) (when k = j), that belong to the same eigenspace. Then and the corresponding eigenvalue is again given by λ j,k = j 2 + k 2 .
The least eigenvalue λ 1,1 = 2 is simple; it is associated to the eigenfunction Φ 1,1 , and it characterizes the Poincaré embedding constant: Other eigenvalues λ j,j may be simple, for instance λ 2,2 = 8, λ 3,3 = 18, or λ 4,4 = 32. But there are also multiple eigenvalues: if k = j, the eigenvalue λ j,k is at least double. Note also that Φ 1,7 , Φ 7,1 , Φ 5,5 are all eigenfunctions associated to the eigenvalue λ = 50, which is then triple. Finally, also higher multiplicities are to be expected. We summarize these results in the following statement. By exploiting the simple geometry of Ω, we are also able to give a fairly precise picture of the existence result whenever f possesses some symmetries.
Furthermore, if f satisfies (7) a.e. in Ω, then: • there exists (at least) one strong solution (u 1 , u 2 , p) ∈ (U ∩ H 2 (Ω)) × H 1 (Ω) of (5) and (6) satisfying the symmetry property (8); Proof. The existence of a weak solution follows from the classical Galerkin method, see e.g. the proof of Theorem 1.4 in [23]. Since f ∈ L 2 (Ω), one has that u ∈ H 2 loc (Ω) by local elliptic regularity. Then one can obtain that u ∈ H 2 (Ω) either by arguing as in the proof of [19,Theorem p.403] or by repeating the reflection argument in the proof of Proposition 3 below, see Remark 2. If f satisfies (7), we proceed as in [17,Theorem 3.4] in order to obtain the existence of at least a symmetric solution satisfying (8). We then notice that (15) solves (5) since in the first scalar equation all the terms maintain their sign whereas in the second scalar equation all the signs change, see again the proof of Proposition 3 below.
Note that (8) and (15) hold pointwise (and not just a.e.) because u ∈ H 2 (Ω) ⊂ C 0 (Ω) since Ω is a planar domain. Note also that a completely similar statement holds, with obvious changes, if f has the "converse" symmetry: In this case, (8) and (15) should be replaced, respectively, by Proposition 2 and this variant will be used to obtain multiplicity results, see Corollary 1 below. For our computer assisted proofs we need much more regularity of the forcing term f in (5). Recalling (12) and given ρ > 1, denote by C ρ ⊂ L 2 (Ω) the space of functions such that The sum in (19) and (20) ranges over {1, 2} × N × N, with the triples (1, 0, k) and (2, j, 0) excluded. Notice that these functions extend analytically to the complex domain | x| < log ρ and | y| < log ρ. For a forcing in this space, we have Proposition 3. Let ρ > 1 and let f ∈ C ρ . Then any solution (u, p) of (5) is analytic in Ω.

Remark 2.
Proposition 3 is not the consequence of elliptic regularity since ∂Ω is merely Lipschitz. Instead, it follows directly from the explicit form of f in (19) which enables one to extend the problem and the solutions by periodicity to the whole plane. In particular, one can extend all the functions to the rectangle R x = (−π, π) × (0, π) ⊃ Ω as follows: (v, q) is analytic in R x except, at most, in its corners. In particular, it is analytic in a neighborhood of the two points (0, 0) and (0, π) and so is the original solution (u, p).
With the same reflection technique, one may also obtain intermediate regularity, that is: if f ∈ H m (Ω) for some integer m ≥ 0, with f as in (19), then any solution (u, p) of (5) satisfies u ∈ H m+2 (Ω) and p ∈ H m+1 (Ω).

Uniqueness for Small Forcing f
In the previous section we have analyzed the main properties of the solutions of (5), but we have not discussed uniqueness of the solution. This is a delicate matter and is the subject of the present section.
When considering the boundary conditions in (5), it is convenient to define the horizontal and vertical edges of ∂Ω, namely and then to introduce the spaces of scalar functions A crucial role in uniqueness statements is then is played by the Sobolev constant Uniqueness for the Navier-Stokes equations (1) is expected under some smallness assumption on the data. We are interested in conditions, as precise as possible, with explicit bounds on f . This is why we give here a sketch of the following uniqueness criterion: Proposition 4. Let S be as in (22) and assume that f (5) and (6) is unique. If f also satisfies (7) (resp. (16)), then the unique solution satisfies (8) (resp. (17)). (5) and (6) satisfying ∇u L 2 < νS, then (u, p) is the unique solution of (5) and (6).

Moreover, for any
Proof. We divide the proofs in three steps.
Step 1. We determine an a priori bound for the solution of (5). To this end, we claim that Notice that the integral is well-defined due to the embedding U ⊂ L 4 (Ω) and Hölder's inequality. If H and V are as in (21), we have (first component) where, in the last step, we used: ∇ · u = 0 (first integral), u 1 = 0 on V (second integral), u 2 = 0 on H (third integral). Similarly, one proceeds for the second component thereby proving (24). Assume that u satisfies (11) and take v = u as test function therein. Then, by (24), Ω (u · ∇)u · u = 0 and we obtain where we used (14).
Step 2. We show that where S is the constant defined in (22). By switching the roles of x and y, we see that S in (22) may be defined equivalently by replacing (26). Then, by the Hölder inequality and these scalar versions of (26), we obtain Step 3. Conclusion, arguing (twice) by contradiction. We assume that there exist two solutions u, v ∈ U of (11), namely We subtract these two equations so that, Take φ = w and use (24) to obtain Then, by the Hölder inequality, using (25) and (26), we infer that By going carefully through the above inequalities, we see that at least one of them is strict, if w = 0. Therefore, the last inequality is also strict if w = 0, which contradicts (23). This shows that w = 0, thereby proving uniqueness whenever (23) holds. Moreover, we notice that if f satisfies (7) as well, then the assertion follows from Proposition 2. Finally, if there exists a solution (u, p) ∈ (U ∩ H 2 (Ω)) × H 1 (Ω) of (5) and (6) satisfying ∇u L 2 < νS, we go back to (27) and obtain with strict inequality if w = 0. This shows that w = 0, and that (u, p) is the unique solution.
Proposition 4 ensures uniqueness of the solution (u, p) of (5), provided that f is small. But it says nothing about uniqueness (or multiplicity) of solutions for large values of f L 2 . This is the reason why Theorem 1 appears important and we prove it in the next section.

Proof of Theorem 1
In this section we prove the following statement, equivalent to Theorem 1, up to a scaling. We maintain here a constant viscosity ν and we let f L 2 vary. Theorem 3. Let ν > 0. There exists a continuous curve α → f α from the positive real line to L 2 (Ω), with f α L 2 = α for all α, such that the problem (5) and (6) with Proof. Consider the eigenfunctions Φ j,k in (13) and, for any positive integer n, define the spaces as well as the orthogonal projectors P n : U → V n . Then we obtain the improved Poincaré constants min v∈Vn ∇v 2 to be compared with (14). Next, we claim that if f α ∈ V n and f α L 2 = α ≤ ν 2 S 1 + n 2 , then (5) admits a unique solution.
In order to prove (29), fix n ≥ 1 and choose any f α satisfying the assumptions therein. The existence of a solution (u, p) of (5) follows from Proposition 2. Then take v = u as test function in (5) and proceed as in the proof of (25) in order to obtain where we used (28). If u = 0, then one of the above inequalities is strict, that is, ∇u L 2 < νS. Uniqueness then follows from the last statement in Proposition 4: this proves (29). Finally, we put γ n = ν 2 S √ 1 + n 2 and δ n = (γ n + γ n+1 )/2, with γ 0 = 0. Then we construct the curve α → f α as follows: -for all n ≥ 0 and α ∈ (γ n , δ n ], we take f α = αΦ 1,n+1 ; -for all n ≥ 0 and α ∈ [δ n , γ n+1 ], we take A sketch of this curve is given in Fig. 3 below. Note that, for any α > 0, f α L 2 = α and uniqueness for (5) follows from (29). This completes the proof. Clearly, it is possible to build many different curves α → f α with the same properties, mostly displaying a zig-zag path as in Fig. 3. In the next section, we describe the behavior along a specific straight line, which involves the bifurcation described in Theorem 2.

Larger Bounds for Uniqueness for Some Special f
Theorem 3 (and Fig. 3) show that there exist a continuous piecewise linear curve in L 2 (Ω) connecting 0 and ∞ such that if f belongs to this curve, then (5) admits a unique solution. One can then wonder whether this holds true for some straight lines as well. In this respect, an interesting case occurs whenever f is an eigenfunction of −Δ in U . In the next statement we improve the uniqueness threshold of Proposition 4 for some special f and we describe a property of possible multiple solutions when f L 2 becomes large.
In view of (30) and (32) we know that (u · ∇)u ∈ G 2 . Therefore, if we test the Eq. (11) satisfied by u with v, we obtain (recall f α = αΦ j,k ) Moreover, if we test (11) satisfied by v with v itself and we recall (24), we get By combining these two identities, we infer that In turn, (33) yields  From (30) we infer that .
Hence, if we test (11), satisfied by v, with u, then we obtain where we used the explicit form of u in (30) and the fact that Φ j,k L 2 = 1. By using again (33), we then This proves the second property in (31). Finally, assume for contradiction that v 2 Δv 1 ≡ v 1 Δv 2 in Ω. Since ∇ · v = 0, we then have that (10), contradicting the second property in (31).
Several comments are in order. Theorem 4 may be complemented with some information about the symmetry of the force and of the solution. If k is even, then Φ j,k satisfies (7) and the solution (u, p) in (30) satisfies (8). If j is even, then Φ j,k satisfies (16) and the solution (u, p) in (30) satisfies (17). Finally, if both j and k are even, then Φ j,k satisfies both (7)-(16) and the solution (u, p) in (30) satisfies both (8)- (17).
Although (31) may apply in more general settings, Theorem 4 gives no practical condition ensuring multiplicity of solutions. It appears out of reach to detect a bifurcation through classical arguments since the linearized equations, around the solution u in (30), read −νΔw + (u · ∇)w + (w · ∇)u = 0 , ∇ · w = 0 in Ω .
These equations, and their weak formulation, are naturally associated to the bilinear form The first inequality shows that a u (·, ·) is continuous but the second formula shows that it fails to be coercive if u is large, especially because of the sign arising from (31) (recall that the sign changes if we switch u and w). Since u in (32) increases linearly with |α|, the Lax-Milgram Theorem cannot be applied to exclude bifurcations. Clearly, for small |α| and u the bilinear form a u is coercive but, in this case, we are in the uniqueness regime where bifurcation is automatically excluded. We searched for bifurcations when f α L 2 = |α| > ν 2 S j 2 + k 2 in the case f α = αΦ j,k , and we found numerical evidence that there are none, see Remark 3. Thus we propose the following conjecture, with a weak and a strong version.

Conjecture 1.
For any j, k ∈ N and α ∈ R, (5) with f α = αΦ j,k admits a unique solution. Or, at least, the branch of solutions (u, p) given in (30) is isolated from other possible solutions of (5).
Finally, let us mention that, by repeating the arguments by Foias-Temam [13, Theorem 2.1], one obtains the following statement.

Scaling and Equivalent Formulation
To study (5) numerically, and to make a computer assisted proof of the existence of solutions, it is convenient to define β = ν −2 and scale (u, p) to (νu, ν 2 u), so that (5) becomes We look for solutions in the space C ρ defined in (19), with ρ = 33 32 . For u = u 1 , u 2 ∈ C ρ we write where the vector fields Ψ i,j,k are defined in (12) and satisfy the boundary conditions in (34). Furthermore, the orthogonal projection onto the space G 1 of solenoidal vector fields is given by Applying P on both sides of the first equation in (34) yields while the remaining equations in (34) show that u = Pu. As it happens for the variational characterization (11), the pressure p does not appear in the equation, when projected onto G 1 : Proof. Since functions in C ρ satisfy the boundary conditions in (34), and Ω is simply connected, then u ∈ C ρ is a solution of (34) if and only if both u and f − (u · ∇)u are divergence free. When restricted to C ρ , the Laplacian commutes with the projection P, therefore if u = F β (u), then u is divergence free and so is f − (u · ∇)u.

Remark 3.
As observed in Theorem 4, when the forcing term is proportional to Φ j,k and sufficiently small, there exists a unique solution also proportional to Φ j,k and, for such solution, the quadratic term is a gradient, therefore its projection on G 1 is zero. Then the eigenvalues of DF β are proportional to β. We have numerical evidence that such eigenvalues lie all on the imaginary axis for some arbitrary choice of β, therefore no eigenvalue can reach the value 1 for any value of β, and since that is a necessary condition to have a bifurcation we conjecture that no bifurcation takes place.

Branches
Denote by D ρ the subspace of C ρ characterized by the symmetry (7). Note that F β (D ρ ) ⊂ D ρ and Pf ∈ D ρ . To prove that (34) admits a solution u(β) ∈ D ρ for each β ∈ [0, 2500], and that the function β → u(β) is analytic, we write all the coefficients in the Fourier expansion of u as Taylor polynomials in β: where a ijkl = 0 if k is odd and b, δ, L are as in Table 1. As a first step, we choose some Fourier-Taylor polynomialū that is an approximate fixed point of F β , and some finite rank operator M such that I − M is an approximate inverse of I − DF β (ū). Then for h ∈ D ρ we define Clearly, if h is a fixed point of N β , then u =ū + Λh is a fixed point of F β and, hence, u solves (34) in view of Lemma 1. Given r > 0 and w ∈ C ρ , let B r (w) = {v ∈ C ρ : v − w ρ < r}. We partition the interval [0, 2500] into six subintervals. The center b and width δ for each subinterval are shown in Table 1.
Then we prove the following lemma with the aid of a computer, see Sect. 5.

Lemma 2.
The following holds for each i ∈ {1, . . . , 6}. Define b, δ, and L as in row i of Table 1. There exist a Fourier-Taylor polynomialū(β) of degree L as described in (38), a bounded linear operator M (β) on D ρ , and positive real numbers ε, r, K satisfying ε + Kr < r, such that holds for all {β ∈ C : |β − b| < δ}. Furthermore, if i < 6, then either where a superscript i refers to the data in row i of Table 1. This lemma, together with the Contraction Mapping Theorem and the Implicit Function Theorem, implies the following: Proposition 6. For each β ∈ [0, 2500] there exists a fixed point u β ∈ B r (u β ) of F β , and the curve β → u β is analytic.

Bifurcations
Let us write the fixed point equation for F β as F(β, u) = 0, where Numerically, we observe that DF(β, .) has a (simple) eigenvalue zero for β near b = 2065 + 743 · 2 −10 . This suggests the possibility of a bifurcation which takes place in a two dimensional submanifold. We parametrize this surface by using the parameter β and a coordinate λ for the range of a suitable onedimensional projection . Then we define a two-parameter family of functions u(α, λ) by solving whereû is a fixed nonzero function in the range of . Forû we choose a Fourier polynomial that approximates the eigenvector of DF(b, .) for the eigenvalue closest to zero. The projection is defined by where u ijk ,û ijk are the Fourier coefficients of u andû, respectively. Our goal is to show that for a rectangle I × J in the parameter space, the Eq. (43) has a smooth and locally unique solution u : Then locally, the solutions of F(β, u) = 0 are determined by the zeros of the function g, We write all the coefficients in the Fourier expansion of u as Taylor polynomials in β, λ: where δ = 2 −11 and λ 1 = 2 −6 . Equation (43) for u = u(β, λ) is equivalent to the fixed point equation for the map F β,λ , defined by As in the last subsection, we use the contraction mapping principle to solve this fixed point problem. In place of the map defined in (39), we use the map M β,λ defined by Here,ū is an approximate fixed point of F b,0 , and M is a finite rank operator such that Λ = I − M is an approximate inverse of I − DF b,0 (ū). By the Implicit Function Theorem, the solution then depends analytically on the two parameters β and λ. Denote by D r (z) a disk in the complex plane of radius r and center z. Let I = D δ (b) and J = D λ1 (0). The following lemma is proved with the aid of a computer; see Sect. 5.
As we had in the previous subsection, this lemma, together with the contraction mapping theorem and the implicit function theorem, implies the following: For every (β, λ) in I × J, the Eq. (43) has a unique solution u(β, λ) in B r (ū(β, λ)), and the map (β, λ) → u = u(β, λ) is analytic. For any given real β ∈ I, a function u in B ∩ −1 (Jû) is a fixed point of F β if and only if u = u(β, λ) for some real λ ∈ J, and g(β, λ) = 0.
This leaves the problem of verifying that the zeros of g correspond to a pitchfork bifurcations. A sufficient set of conditions for the existence of such a bifurcation is given below, see [6,Lemma 3.4] for a proof.
If f is any differentiable function of two variables, denote byḟ and f the partial derivatives of f with respect to the first and second argument, respectively. Let g (β 2 , 0) < 0. Then g(β, λ) = λG(β, λ) for some C 2 function G, and the solution set of G(β, λ) = 0 in I × J is the graph of a C 2 function β = a(λ), defined on a proper subinterval J 0 of J. This function takes the value β 2 at the endpoints of J 0 , and satisfies β 1 < a(z 2 ) < β 2 at all interior points of J 0 , which includes the origin.
The following lemma is proved with the aid of a computer; see Sect. 5.
The proof of Theorem 2 is concluded by joining the results of Propositions 6 and 7 , together with Lemmas 4 and 5 .

Verifying Lemmas 2, 3 and 5
In this section we provide a rough description of the proofs of Lemmas 2, 3 and 5 . A complete and detailed description can be found in [5].
The proofs are carried out with the aid of a computer. Our proof of Lemma 2 involves the following steps. Consider the parameter values (b, δ) from a fixed but arbitrary row in Table 1. As a first step, we determine a Fourier-Taylor polynomialū β which is an approximate fixed point of F β and an approximate inverse of I − DF β (ū β ) of the form Λ(β) = I − M (β), where M (β) has finite rank and is a polynomial in β. The remaining steps are rigorous: we compute an upper bound ε on the norm of F β (0), and an upper bound K on the operator norm of DF β (h) that holds for all h of norm ε or less. This is done simultaneously for all values of {β ∈ C : |β − b| < δ}. After verifying that K < 7 8 , we choose a positive r < 8ε in such a way that ε + Kr < r.
The same approach is used to prove Lemma 3. Concerning the computation of bounds, we refer to Sect. 5.2 below. The proof of Lemma 5 consists in verifying some bounds on the derivatives of the function g. Let h(z 1 , z 2 ) = g(b + β 1 z 1 , γ 0 + γ 1 z 2 ), and note that h is analytic in {z ∈ C 2 : |z 1 | < 1, |z 2 | < 1} by Proposition 6. Given the Fourier-Taylor polynomialū(β, λ) of Lemma 3, we compute a Taylor polynomial P (z 1 , z 2 ) of degree 4 and a bound E, such that that sup |z1|<1,|z2|<1 (51) Using = 205/256, this bound is applied (separately in each variable) to verify the conditions in Lemma 4.

Technicalities
The methods used here can be considered perturbation theory: given an approximate solution, prove bounds that guarantee the existence of a true solution nearby. But the approximate solutions needed here are too complex to be described without the aid of a computer, and the number of estimates involved is far too large. The first part (finding approximate solutions) is a strictly numerical computation. The rigorous part is still numerical, but instead of truncating series and ignoring rounding errors, it produces guaranteed enclosures at every step along the computation. This part of the proof is written in the programming language Ada [3]. 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 [5].
In the present 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).
To represent balls in a real Banach algebra X with unit 1 we use pairs S=(S.C,S.R), where S.C is a representable number (Rep) and S.R a nonnegative representable number (Radius). The corresponding ball in X is S, When X = R the data type described above is called Ball. Our bounds on some standard functions involving the type Ball are defined in the packages 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.
The computation and validation of branches involves Taylor series in one variable, which are represented by the type Taylor1 with coefficients of type Ball. The definition of the type and its basic procedures are in the package Taylors1. Given a Radius ρ, consider the space T ρ of all real analytic functions g(t) = n g n t n on the interval |t| < ρ, obtained by completing the space of polynomials with respect to the norm g ρ = n |g n |ρ n . Given a positive integer D, a Taylor1 is a triple P=(P.C,P.F,P.R), where P.F is a nonnegative integer, P.R = ρ, and P.C is an array(0..D) of Ball. The corresponding set in Taylor1, T ρ is defined as where m = min(P.F, D + 1). For the operations that we need in our proof, this type of enclosure allows for simple and efficient bounds. We also need a type Taylor2, representing Taylor series of two variables, for the proof of Lemmas 3 and 5, which is defined similarly in the package Taylors2. Finally, solutions of the Eq. (34) are given as Fourier series in two variables, represented by the type Fourier2 with coefficients in some Banach algebra with unit X . Consider the space F ρ of all real analytic functions in two variables periodic of period 2π in both variables, obtained by completing the space of Fourier polynomials p(w, z) = jk p jk e i(jw+kz) , p jk ∈ X , with respect to the norm p ρ = jk p jk ρ j+k . The type Fourier2 consists of a triple F=(F.T,F.C,F.E), where F.T is a record identifying the parity (see below), F.C is an array(0..K,0..K) of Ball, and F.E is an array(0..2*K,0..2*K) of Radius.
Let f j,k = The corresponding set F, F ρ is the set of all function u = p + h ∈ F ρ , where with p J,K ∈ F.C(J, K), X and h J,K ≤ F.E(J, K), for all J, K ≥ 1.
To be more precise, instead of (53) we use cosine series (parity 0) and sine series (parity 1). The parities are specified by the component F.T. The definition of the type and its basic procedures are in the package Fouriers2. For the operations that we need in our proof, this type of enclosure allows for simple and efficient bounds. We note that Fourier2 (just like Taylor1 and Taylor2) allows a generic type Scalar for its coefficients; and this Scalar can be again a Taylor (or Fourier) series. This feature makes it easy to represent Fourier series whose coefficients depend on parameters.
In particular, our enclosures for the functions described in (38) use an instantiation of Fourier2 with Scalar = TScalar, where TScalar is an istantiation of Taylor1 with Scalar = Ball. The corresponding type VF describes (sets of) vector field in the spaces C ρ and D ρ . These types are defined in the package Taylors1.Foor. Types and bounds that are specific to the Navier-Stokes equation (34) are defined in the child package Fouriers2.VFs. This includes the type VF (a pair of Fouriers2) which is used to specify enclosures on vector fields in the spaces C ρ , a bound NegInvLap on −Δ −1 , a bound ProjDivFree on the projection onto divergence-free vector fields, a bound NonLinearAdv on the map u → f − u · ∇u, and a bound DNonLinearAdv on the derivative of NonLinearAdv.
Analogous types, with Taylor1 replaced by Taylor2, are defined in the package Taylors2.Foor. These types are used as enclosures for for the functions described in (46). A bound on the map F β defined in (37) is implemented by the procedure GMap in the package Taylors1.Foor.Fix. Defining and estimating a contraction like N β is a common task in many of our computer-assisted proofs. An implementation is done via two generic packages, Linear and Linear.Contr. For a description of this process we refer to [7]. The only problem-dependent parts here are the bound GMap, and the type VFMode defined in Taylors1.Foor. The same approach is used for the estimates (49) on the contraction C β,λ defined in (48). The bound GMap on the map F β,λ defined in (47) is implemented in the package Taylors2.Foor.Fix.

Concluding Remarks and Open Problems
In this paper we were able to exhibit a bifurcation for the Navier-Stokes equations under Navier boundary conditions (5) for a particular force f . We exploited the geometry of the square Ω and we highlighted possible symmetry breaking of the solutions. We were also able to follow some branches of arbitrarily large forces for which uniqueness for (5) holds. Although our results shed some light on these fairly difficult topics, many problems remain open. Among them are the following. If we take w 0 ∈ H 1 V independent of y we find a boundary value problem for an ODE, that is, w 0 (x) + w 0 (x) 3 = 0 in (0, π) , w 0 (0) = w 0 (π) = 0 , whose solution is given by where cn denotes the Jacobi cosine, see [1]. Multiplying (54) by w 0 and integrating by parts yields Therefore, if we view w 0 as a function of two variables, we obtain the following bound for S in (22): We do not know if the minimizers for (22) depend on x only. Is the inequality (55) an equality? The precise knowledge of S would allow to obtain more precise statements for uniqueness, see Proposition 4 and Theorem 4.
Several further problems were left open. We summarize them in the following list.
• Is it ν 0 = ν 1 in Theorem 2 and for any f ? In other words, for general f is uniqueness of the solution of (5) lost due to a bifurcation or to the appearance of isolated solutions, far away from the symmetric branch? This uncertainty also occurs for (possibly inhomogeneous) Dirichlet problems, see Théorème 3.3 in [14] and the subsequent picture therein. • Do the results obtained in the present paper hold in more general situations? That is, can Theorems 1 and 2 be extended to any planar bounded domain Ω, and to more general forces f ? Can the same technique be extended to 3D domains? • Is it possible to find g ∈ L 2 (Ω) such that Eq. (5) with f = αg admits a unique solution for all α ∈ [0, +∞)? This problem is related to Conjecture 1. • The bifurcation branches obtained in Theorem 5 only contain non-symmetric solutions. Is it possible to obtain multiplicity of symmetric solutions satisfying (8)? Alternatively, if f satisfies (7), is there always a unique symmetric solution? This problem appears quite challenging also from a numerical point of view: in a symmetric channel containing a circular cylinder, Sahin-Owens [28, Fig.6] (see branches 1, 3, 5 therein) numerically found different symmetric solutions for suitable Reynolds numbers, but for different ratio between the width of the channel and the diameter of the cylinder. Form a theoretical point of view, the answer to this question would enable to understand whether the bifurcation from the symmetric solution is both a necessary and a sufficient condition for the appearance of a lift force on a symmetric obstacle immersed in the fluid, see [11,17]. • Can the construction of eigenvectors of the Stokes problem, as described in Remark 1, be extended to more general domains? Here we exploited the geometry of the square Ω, while in general domains this issue appears to be much more involved. Let us mention that also under the no-slip boundary conditions (2), explicit eigenvalues and eigenfunctions are known only in special domains, see [27] and references therein.