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

Abstract

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.

1 Introduction and Main Results

Let $\Omega \subset {\mathbb {R}}^2$ be a bounded domain and consider the stationary Navier–Stokes equations

$$\begin{aligned} -\nu \Delta u+(u\cdot \nabla )u+\nabla p=f\ ,\quad \nabla \cdot u=0\quad \text{ in } \Omega \ , \end{aligned}$$

(1)

that model the steady-state motion of an incompressible viscous fluid: u is its velocity, p its pressure, f is an external force, $\nu >0$ is the kinematic viscosity. Equation (1) need to be complemented with some boundary constraint, the most common being the no-slip boundary conditions

$$\begin{aligned} u=0\quad \text{ on } \partial \Omega \, , \end{aligned}$$

(2)

that are physically reasonable if the boundary of $\Omega$ 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 $\partial \Omega$ and the source f. Much more delicate is the uniqueness, which is guaranteed only for large viscosities $\nu$ 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 $\nu$. 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 to slip, and with the tangential component of the strain tensor proportional to the tangential component of the fluid velocity on the boundary. The Navier boundary conditions read

$$\begin{aligned} u\cdot n=(\mathbf{D}u,n)\cdot \tau =0\quad \text{ on } \partial \Omega ,\quad \mathbf{D}u={1\over 2}(\nabla u+\nabla ^Tu), \end{aligned}$$

(3)

where n is the outward normal vector to $\partial \Omega$ while $\tau$ 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\in L^2(\Omega )$. Concerning regularity results, we mention the works by Beir$\tilde{\mathrm{a}}$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 $\nu$, 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:

$$\begin{aligned} \Omega =(0,\pi )^2\, . \end{aligned}$$

(4)

It is known [10] that for flat boundaries, the conditions (3) become of mixed Dirichlet-Neumann type and (1)–(3) read

$$\begin{aligned} \left\{ \begin{array}{l} -\nu \Delta u+(u\cdot \nabla )u+\nabla p=f\quad \text{ in } \Omega ,\\ \nabla \cdot u=0\quad \text{ in } \Omega ,\\ u_1=\partial _xu_2=0\quad \text{ on } \{0,\pi \}\times (0,\pi ),\\ u_2=\partial _yu_1=0\quad \text{ on } (0,\pi )\times \{0,\pi \}. \end{array}\right. \end{aligned}$$

(5)

The mixed conditions are natural, since they guarantee that the boundary term vanishes after an integration by parts:

$$\begin{aligned} \int _\Omega \nabla u:\nabla v=-\int _\Omega v\Delta u+\int _{\partial \Omega }(v_1\partial _nu_1+v_2\partial _nu_2)=-\int _\Omega v\Delta u. \end{aligned}$$

The pressure p is defined up to an additive constant so that one can fix its mean value, for instance

$$\begin{aligned} \int _\Omega p=0\, . \end{aligned}$$

(6)

In Sect. 2 we take advantage of the geometry of the square $\Omega$ 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 $\nu \mapsto f_\nu$ from the positive real line $(0,\infty )$ to $L^2(\Omega )$, with $\Vert f_\nu \Vert _{L^2}=1$ for all $\nu$, such that the problem (5), (6) with $f=f_\nu$ admits a unique solution $(u,p)\in H^2(\Omega )\times H^1(\Omega )$.

For other curves $\nu \mapsto f_\nu$, 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=\pi /2$, that is,

$$\begin{aligned} f_1(x,\pi -y)=f_1(x,y),\quad f_2(x,\pi -y)=-f_2(x,y)\quad \forall (x,y)\in \Omega \, . \end{aligned}$$

(7)

One then expects that (at least) a solution of (5) has the same symmetry, namely

$$\begin{aligned} u_1(x,\pi -y)=u_1(x,y),\quad u_2(x,\pi -y)=-u_2(x,y),\quad p(x,\pi -y)=p(x,y),\quad \forall (x,y)\in \Omega \, . \end{aligned}$$

(8)

We take f analytic and concentrated near the center of the square $\Omega$. Then, as the viscosity decreases from $\nu =+\infty$, we are able to prove nonuniqueness via bifurcations and symmetry breaking. More precisely, we have

Theorem 2

There exists an analytic function f, with $\Vert f\Vert _{L^2}=1$, satisfying (7) and such that:

• for all $\nu >0$, (5) and (6) admits an analytic solution $(u_\nu ,p_\nu )$ satisfying (8), the curve $\nu \mapsto u_\nu$ is analytic and $u_\nu$ is isolated within the subspace of symmetric functions satisfying (8);

• there exists $\nu _1>0$ such that the solution of (5) and (6) is unique for all $\nu >\nu _1$;

• at some positive $\nu _0\le \nu _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 $\nu$ follows from a general statement as well, see Proposition 4 where a lower bound for $\nu _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 $\nu _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 $\nu \mapsto u_\nu$ is known since Foias-Temam [14].

The secondary branch mentioned in Theorem 2 has been determined rigorously only close to the bifurcation value $\nu _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.

Fig. 1

Sketch of the bifurcation branch: in black the parts proved through computer assistance, in gray the parts obtained numerically

Fig. 2

Symmetric and non symmetric solutions at $\nu =1/50$

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 $\nu$: 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.

2 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

$$\begin{aligned} \begin{array}{ccc} G_1=\{v\in L^2(\Omega );\, \nabla \cdot v=0,\, v\cdot n=0 \text{ on } \partial \Omega \}\, ,\qquad U=G_1\cap H^1(\Omega ),\\ G_2=\{v\in L^2(\Omega );\, \exists g\in H^1(\Omega ),\, v=\nabla g\}, \end{array} \end{aligned}$$

(9)

where $v\mapsto n\cdot v$ denotes the boundary-normal trace operator. It is known that

$$\begin{aligned} L^2(\Omega )=G_1\oplus G_2\ ,\qquad G_1\perp G_2\ , \end{aligned}$$

(10)

where orthogonality is with respect to the scalar product in $L^2(\Omega )$.

Assuming that $f\in L^2(\Omega )$, we say that $u\in U$ is a weak solution of (5) if

$$\begin{aligned} \nu \int _\Omega \nabla u:\nabla v+\int _\Omega (u\cdot \nabla )u\cdot v=\int _\Omega fv\qquad \forall v\in U\, . \end{aligned}$$

(11)

Weak solutions may also be found under the less restrictive assumption that $f\in H^{-1}(\Omega )$ but, since our main concern are smooth solutions, we assume here (at least) that $f\in L^2(\Omega )$. In this case, as we shall see, the solution has the additional regularity $u\in H^2(\Omega )$. Then, (11) implies that the first equation in (5), seen as an equality of functions in $L^2(\Omega )$, is satisfied when projected onto $G_1$. Due to (10), $-\nu \Delta u+(u\cdot \nabla )u-f\in G_2$ and, therefore, that there exists $p\in H^1(\Omega )$ 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 $\Omega$ will also be used to determine explicitly the eigenfunctions of $-\Delta$ in U, satisfying the boundary conditions in (5). Let

$$\begin{aligned} \Psi _{1,j,k}(x,y)=\bigl (\sin (jx)\cos (ky)\,\mathbf{,}\,0\bigr ),\quad \Psi _{2,j,k}(x,y)=\bigl (0,\cos (jx)\sin (ky)\bigr ). \end{aligned}$$

(12)

Then the $L^2-$normalized eigenfunctions of $-\Delta$ in U are given by

$$\begin{aligned} \begin{aligned} \Phi _{j,k}(x,y)&={2\over \pi \sqrt{j^2+k^2}}(-k\Psi _{1,j,k}(x,y)+j\Psi _{2,j,k}(x,y))\\&={2\over \pi \sqrt{j^2+k^2}} \left( \begin{array}{c} -k\sin (jx)\cos (ky)\\ j\cos (jx)\sin (ky) \end{array}\right) ,\qquad (j,k=1,2,3,\ldots ). \end{aligned} \end{aligned}$$

(13)

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 $\varphi$ of $-\Delta$ in $H^1_0(\Omega )$, that is,

$$\begin{aligned} -\Delta \varphi =\lambda \varphi \quad \text{ in }\quad \Omega \, ,\qquad \varphi =0\quad \text{ on }\quad \partial \Omega \, . \end{aligned}$$

The eigenfunctions with separated variables are given by $\varphi _{j,k}(x,y)=\sin (jx)\sin (ky)$ with corresponding eigenvalues $\lambda _{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\ne j$), that belong to the same eigenspace. Then

$$\begin{aligned} \Phi _{j,k}={2\over \pi \sqrt{j^2+k^2}} \left( \begin{array}{c} -\partial _y\varphi _{j,k}\\ \partial _x\varphi _{j,k} \end{array}\right) \qquad (j,k=1,2,3,\ldots ) \end{aligned}$$

and the corresponding eigenvalue is again given by $\lambda _{j,k}=j^2+k^2$. $\square$

The least eigenvalue $\lambda _{1,1}=2$ is simple; it is associated to the eigenfunction $\Phi _{1,1}$, and it characterizes the Poincaré embedding constant:

$$\begin{aligned} \min _{u\in U}\ {\Vert \nabla u\Vert _{L^2}^2\over \Vert u\Vert _{L^2}^2}=2\, . \end{aligned}$$

(14)

Other eigenvalues $\lambda _{j,j}$ may be simple, for instance $\lambda _{2,2}=8$, $\lambda _{3,3}=18$, or $\lambda _{4,4}=32$. But there are also multiple eigenvalues: if $k\ne j$, the eigenvalue $\lambda _{j,k}$ is at least double. Note also that $\Phi _{1,7}$, $\Phi _{7,1}$, $\Phi _{5,5}$ are all eigenfunctions associated to the eigenvalue $\lambda =50$, which is then triple. Finally, also higher multiplicities are to be expected. We summarize these results in the following statement.

Proposition 1

The functions $\Phi _{j,k}$ in (13) with $j,k=1,2,3,\ldots$ are a complete $L^2$-orthonormal basis of U, consisting of eigenfunctions of $-\Delta$. The associated eigenvalues are $\lambda _{j,k}=j^2+k^2$.

By exploiting the simple geometry of $\Omega$, we are also able to give a fairly precise picture of the existence result whenever f possesses some symmetries.

Proposition 2

Assume that $f=(f_1,f_2)\in L^2(\Omega )$. Then:

$\bullet$ any $u\in U$ satisfying (11) is a strong solution $u\in U\cap H^2(\Omega )$ of (5);

$\bullet$ there exists a solution $u\in U\cap H^2(\Omega )$ of (5), coupled with a unique $p\in H^1(\Omega )$ satisfying (6).

Furthermore, if f satisfies (7) a.e. in $\Omega$, then:

$\bullet$ there exists (at least) one strong solution $(u_1,u_2,p) \in (U\cap H^2(\Omega ))\times H^1(\Omega )$ of (5) and (6) satisfying the symmetry property (8);

$\bullet$ if $(u_1,u_2,p) \in (U\cap H^2(\Omega ))\times H^1(\Omega )$ is a strong solution of (5), so is $(v_1,v_2,q)$ with

$$\begin{aligned} v_1(x,y)=u_1(x,\pi -y),\quad v_2(x,y)=-u_2(x,\pi -y),\quad q(x,y)=p(x,\pi -y),\quad \forall (x,y)\in \Omega . \end{aligned}$$

(15)

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\in L^2(\Omega )$, one has that $u\in H^2_{\mathrm{loc}}(\Omega )$ by local elliptic regularity. Then one can obtain that $u\in H^2(\Omega )$ 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.$\square$

Note that (8) and (15) hold pointwise (and not just a.e.) because $u\in H^2(\Omega )\subset C^0({\overline{\Omega }})$ since $\Omega$ is a planar domain. Note also that a completely similar statement holds, with obvious changes, if f has the “converse” symmetry:

$$\begin{aligned} f_1(\pi -x,y)=-f_1(x,y),\quad f_2(\pi -x,y)=f_2(x,y)\qquad \text{ for } \text{ a.e. } (x,y)\in \Omega \, . \end{aligned}$$

(16)

In this case, (8) and (15) should be replaced, respectively, by

$$\begin{aligned} u_1(\pi -x,y)=-u_1(x,y),\quad u_2(\pi -x,y)=u_2(x,y),\quad p(\pi -x,y)=p(x,y),\quad \forall (x,y)\in \Omega \, , \end{aligned}$$

(17)

$$\begin{aligned} v_1(x,y)=-u_1(\pi -x,y),\quad v_2(x,y)=u_2(\pi -x,y),\quad q(x,y)=p(\pi -x,y),\quad \forall (x,y)\in \Omega \, . \end{aligned}$$

(18)

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 $\rho >1$, denote by ${\mathcal {C}}_\rho \subset L^2(\Omega )$ the space of functions

$$\begin{aligned} f=\sum _{i,j,k}f_{ijk}\Psi _{i,j,k}\qquad (f_{ijk}\in {\mathbb {C}}) \end{aligned}$$

(19)

such that

$$\begin{aligned} \Vert f\Vert _\rho :=\sum _{i,j,k}|f_{ijk}|\rho ^{j+k}<+\infty . \end{aligned}$$

(20)

The sum in (19) and (20) ranges over $\{1,2\}\times {{\mathbb {N}}}\times {{\mathbb {N}}}$, with the triples (1, 0, k) and (2, j, 0) excluded. Notice that these functions extend analytically to the complex domain $|\mathfrak {I}x|<\log \rho$ and $|\mathfrak {I}y|<\log \rho$. For a forcing in this space, we have

Proposition 3

Let $\rho >1$ and let $f\in {\mathcal {C}}_\rho$. Then any solution (u, p) of (5) is analytic in ${\overline{\Omega }}$.

Remark 2

Proposition 3is not the consequence of elliptic regularity since $\partial \Omega$ 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=(-\pi ,\pi )\times (0,\pi )\supset \Omega$ as follows:

$$\begin{aligned} g_1(x,y)=\left\{ \begin{array}{ll} f_1(x,y) &{} \text{ if } x\ge 0\\ -f_1(-x,y) &{} \text{ if } x\le 0 \end{array}\right. \qquad g_2(x,y)=\left\{ \begin{array}{ll} f_2(x,y) &{} \text{ if } x\ge 0\\ f_2(-x,y) &{} \text{ if } x\le 0 \end{array}\right. \ , \\ v_1(x,y)=\left\{ \begin{array}{ll} u_1(x,y) &{} \text{ if } x\ge 0\\ -u_1(-x,y) &{} \text{ if } x\le 0 \end{array}\right. \qquad v_2(x,y)=\left\{ \begin{array}{ll} u_2(x,y) &{} \text{ if } x\ge 0\\ u_2(-x,y) &{} \text{ if } x\le 0 \end{array}\right. \ , \\ q(x,y)=\left\{ \begin{array}{ll} p(x,y) &{} \text{ if } x\ge 0\\ p(-x,y) &{} \text{ if } x\le 0 \end{array}\right. \ . \end{aligned}$$

Then g is analytic in $R_x$ and (v, q) solves the problem

$$\begin{aligned} \left\{ \begin{array}{l} -\nu \Delta v+(v\cdot \nabla )v+\nabla q=g\quad \text{ in } R_x\\ \nabla \cdot v=0\quad \text{ in } R_x\\ v_1=\partial _xv_2=0\quad \text{ on } \{-\pi ,\pi \}\times (0,\pi )\\ v_2=\partial _yv_1=0\quad \text{ on } (-\pi ,\pi )\times \{0,\pi \}\, , \end{array}\right. \end{aligned}$$

since in the first scalar equation all the terms have maintained their sign whereas in the second scalar equation all the signs have changed; the divergence-free condition is also trivially fulfilled. Therefore, (v, q) is analytic in $\overline{R_x}$ except, at most, in its corners. In particular, it is analytic in a neighborhood of the two points (0, 0) and $(0,\pi )$ and so is the original solution (u, p).

With the same reflection technique, one may also obtain intermediate regularity, that is: if $f\in H^m(\Omega )$ for some integer $m\ge 0$, with f as in (19), then any solution (u, p) of (5) satisfies $u\in H^{m+2}(\Omega )$ and $p\in H^{m+1}(\Omega )$. $\square$

3 Uniqueness of the Solution

3.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.

As a straightforward consequence of Proposition 2, we have a nonuniqueness criterion.

Corollary 1

Assume that $f\in L^2(\Omega )$ satisfies (7) (resp. (16)). If (5) admits an asymmetric solution $(u,p)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$ violating (8) (resp. (17)), then (5) admits at least two more solutions: its reflection (v, q) given by (15) (resp. (18)) and a symmetric solution satisfying (8) (resp. (17)).

Moreover, if $f\in L^2(\Omega )$ satisfies both (7)–(16) and (5) admits an asymmetric solution $(u,p)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$ violating both (8)–(17), then (5) admits at least four more solutions.

When considering the boundary conditions in (5), it is convenient to define the horizontal and vertical edges of $\partial \Omega$, namely

$$\begin{aligned} H=(0,\pi )\times \{0,\pi \}\quad \text{ and }\quad V=\{0,\pi \}\times (0,\pi )\, , \end{aligned}$$

(21)

and then to introduce the spaces of scalar functions

$$\begin{aligned} H^1_{\texttt {H}}=\{w\in H^1(\Omega );\, w=0 \text{ on } H\}\quad \text{ and }\quad H^1_{\texttt {V}}=\{w\in H^1(\Omega );\, w=0 \text{ on } V\}. \end{aligned}$$

To be more precise, $H^1_H$ and $H^1_V$ are defined as the closures in $H^1(\Omega )$ of $C^\infty _c([0,\pi ]\times (0,\pi ))$ and $C^\infty _c((0,\pi )\times [0,\pi ])$) respectively.

A crucial role in uniqueness statements is then is played by the Sobolev constant

$$\begin{aligned} S=\ \min _{w\in H^1_{\texttt {V}}}\ {\Vert \nabla w\Vert _{L^2}^2 \over \Vert w\Vert _{L^4}^2}\, . \end{aligned}$$

(22)

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=(f_1,f_2)\in L^2(\Omega )$. If

$$\begin{aligned} \Vert f\Vert _{L^2}\le \nu ^2\, S\, \sqrt{2} , \end{aligned}$$

(23)

then the solution $(u,p)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$ of (5) and (6) is unique. If f also satisfies (7) (resp. (16)), then the unique solution satisfies (8) (resp. (17)).

Moreover, for any $f\in L^2(\Omega )$, if there exists a solution $(u,p)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$ of (5) and (6) satisfying $\Vert \nabla u\Vert _{L^2}<\nu S$, then (u, p) is the unique solution of (5) and (6).

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

$$\begin{aligned} \int _\Omega (u\cdot \nabla )v\cdot v=0\qquad \forall (u,v)\in U\times U\, . \end{aligned}$$

(24)

Notice that the integral is well-defined due to the embedding $U\subset L^4(\Omega )$ and Hölder’s inequality. If H and V are as in (21), we have (first component)

$$\begin{aligned} \int _\Omega (u_1\partial _xv_1+u_2\partial _yv_1)v_1= & {} {1\over 2}\int _\Omega u_1\partial _xv_1^2+{1\over 2}\int _\Omega u_2\partial _yv_1^2\\ \text{(by } \text{ parts } \text{ with } \text{ respect } \text{ to } \,\, x \,\,\text{ or }\,\, y)= & {} -{1\over 2}\int _\Omega v_1^2(\nabla \cdot u)+{1\over 2}\int _V u_1v_1^2+{1\over 2}\int _H u_2v_1^2=0, \end{aligned}$$

where, in the last step, we used: $\nabla \cdot 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), $\int _\Omega (u\cdot \nabla )u\cdot u=0$ and we obtain

$$\begin{aligned} \nu \int _\Omega |\nabla u|^2=\int _\Omega fu\le \Vert f\Vert _{L^2}\, \Vert u\Vert _{L^2}\le \Vert f\Vert _{L^2}\, {\Vert \nabla u\Vert _{L^2}\over \sqrt{2} }\ \Longrightarrow \ \Vert \nabla u\Vert _{L^2}\le {\Vert f\Vert _{L^2}\, \over \nu \sqrt{2} }\, , \end{aligned}$$

(25)

where we used (14).

Step 2. We show that

$$\begin{aligned} S\Vert u\Vert _{L^4}^2\le \Vert \nabla u\Vert _{L^2}^2\qquad \forall u\in U\, , \end{aligned}$$

(26)

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 $H^1_{\texttt {H}}$ with $H^1_{\texttt {V}}$. Hence, the (scalar) components $u_1\in H^1_{\texttt {V}}$ and $u_2\in H^1_{\texttt {H}}$ of a vector $u=(u_1,u_2)\in U$, both satisfy (26). Then, by the Hölder inequality and these scalar versions of (26), we obtain

$$\begin{aligned} S^2\Vert u\Vert ^4_{L^4}\le S^2\left( \Vert u_1\Vert ^2_{L^4}+\Vert u_2\Vert ^2_{L^4}\right) ^2\le \left( \Vert \nabla u_1\Vert ^2_{L^2}+\Vert \nabla u_2\Vert ^2_{L^2}\right) ^2= \Vert \nabla u\Vert ^4_{L^2}, \end{aligned}$$

which proves (26).

Step 3. Conclusion, arguing (twice) by contradiction. We assume that there exist two solutions $u,v\in U$ of (11), namely

$$\begin{aligned} \nu \int _\Omega \nabla u:\nabla \phi +\int _\Omega (u\cdot \nabla )u\cdot \phi =\int _\Omega f\phi \, ,\qquad \nu \int _\Omega \nabla v:\nabla \phi +\int _\Omega (v\cdot \nabla )v\cdot \phi =\int _\Omega f\phi \qquad \quad \forall \phi \in U\, . \end{aligned}$$

We subtract these two equations so that, if $w=u-v$, we get

$$\begin{aligned} \nu \int _\Omega \nabla w:\nabla \phi =\int _\Omega (v\cdot \nabla )v\cdot \phi -\int _\Omega (u\cdot \nabla )u\cdot \phi \qquad \forall \phi \in U\, . \end{aligned}$$

Take $\phi =w$ and use (24) to obtain

$$\begin{aligned} \nu \int _\Omega |\nabla w|^2=-\int _\Omega (w\cdot \nabla )u\cdot w\, . \end{aligned}$$

Then, by the Hölder inequality, using (25) and (26), we infer that

$$\begin{aligned} \nu \Vert \nabla w\Vert _{L^2}^2\le \Vert w\Vert _{L^4}\Vert \nabla u\Vert _{L^2}\Vert w\Vert _{L^4}\le {\Vert \nabla w\Vert _{L^2}^2\over S} \, {\Vert f\Vert _{L^2} \over \nu \sqrt{2} }\, . \end{aligned}$$

(27)

By going carefully through the above inequalities, we see that at least one of them is strict, if $w\ne 0$. Therefore, the last inequality is also strict if $w\ne 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)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$ of (5) and (6) satisfying $\Vert \nabla u\Vert _{L^2}<\nu S$, we go back to (27) and obtain

$$\begin{aligned} \nu \Vert \nabla w\Vert _{L^2}^2\le \Vert w\Vert _{L^4}\Vert \nabla u\Vert _{L^2}\Vert w\Vert _{L^4}\le {\Vert \nabla w\Vert _{L^2}^2\over S} \, \nu S=\nu \Vert \nabla w\Vert _{L^2}^2\, , \end{aligned}$$

with strict inequality if $w\ne 0$. This shows that $w=0$, and that (u, p) is the unique solution.$\square$

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 $\Vert f\Vert _{L^2}$. This is the reason why Theorem 1 appears important and we prove it in the next section.

3.2 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 $\nu$ and we let $\Vert f\Vert _{L^2}$ vary.

Theorem 3

Let $\nu >0$. There exists a continuous curve $\alpha \mapsto f_\alpha$ from the positive real line to $L^2(\Omega )$, with $\Vert f_\alpha \Vert _{L^2}=\alpha$ for all $\alpha$, such that the problem (5) and (6) with $f=f_\alpha$ admits a unique solution $(u,p)\in H^2(\Omega )\times H^1(\Omega )$.

Proof

Consider the eigenfunctions $\Phi _{j,k}$ in (13) and, for any positive integer n, define the spaces

$$\begin{aligned} V_n=\mathrm{span}\{\Phi _{1,k};\, k\ge n\}\, , \end{aligned}$$

as well as the orthogonal projectors $P_n:U\rightarrow V_n$. Then we obtain the improved Poincaré constants

$$\begin{aligned} \min _{v\in V_n}\ {\Vert \nabla v\Vert _{L^2}^2\over \Vert v\Vert _{L^2}^2}=1+n^2\, , \end{aligned}$$

(28)

to be compared with (14). Next, we claim that

$$\begin{aligned} \text{ if }\, f_\alpha \in V_n \,\text{ and }\, \Vert f_\alpha \Vert _{L^2}=\alpha \le \nu ^2 S\sqrt{1+n^2},\, \text{ then }\, (5)\, \text{ admits } \text{ a } \text{ unique } \text{ solution. } \end{aligned}$$

(29)

In order to prove (29), fix $n\ge 1$ and choose any $f_\alpha$ 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

$$\begin{aligned} \nu \int _\Omega |\nabla u|^2=\int _\Omega f_\alpha u\le \alpha \, \Vert P_nu\Vert _{L^2}\le \alpha \, {\Vert \nabla u\Vert _{L^2}\over \sqrt{1+n^2}}\le \nu ^2 S\Vert \nabla u\Vert _{L^2} \ \Longrightarrow \ \Vert \nabla u\Vert _{L^2}\le \nu S\, , \end{aligned}$$

where we used (28). If $u\ne 0$, then one of the above inequalities is strict, that is, $\Vert \nabla u\Vert _{L^2}<\nu S$. Uniqueness then follows from the last statement in Proposition 4: this proves (29).

Finally, we put $\gamma _n=\nu ^2S\sqrt{1+n^2}$ and $\delta _n=(\gamma _n+\gamma _{n+1})/2$, with $\gamma _0=0$. Then we construct the curve $\alpha \mapsto f_\alpha$ as follows:

– for all $n\ge 0$ and $\alpha \in (\gamma _n,\delta _n]$, we take $f_\alpha =\alpha \Phi _{1,n+1}$;

– for all $n\ge 0$ and $\alpha \in [\delta _n,\gamma _{n+1}]$, we take

$$\begin{aligned} f_\alpha =\alpha \, {(\gamma _{n+1}-\alpha )\Phi _{1,n+1}+(\alpha -\delta _n)\Phi _{1,n+2} \over \sqrt{(\gamma _{n+1}-\alpha )^2+(\alpha -\delta _n)^2} } \in V_{n+1}\, . \end{aligned}$$

A sketch of this curve is given in Fig. 3 below. Note that, for any $\alpha >0$, $\Vert f_\alpha \Vert _{L^2}=\alpha$ and uniqueness for (5) follows from (29). This completes the proof. $\square$

Fig. 3

Continuous path $\alpha \mapsto f_\alpha$ connecting 0 and $\infty$ for which (5) admits a unique solution

Clearly, it is possible to build many different curves $\alpha \mapsto f_\alpha$ 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.

3.3 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(\Omega )$ connecting 0 and $\infty$ 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 $-\Delta$ 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 $\Vert f\Vert _{L^2}$ becomes large.

Theorem 4

Let $j,k\ge 1$ and define

$$\begin{aligned} \Upsilon _{j,k}(x,y)={4\big (k^2\sin ^2(jx)+j^2\sin ^2(ky)\big ) \over \pi ^2(j^2+k^2)}\, . \end{aligned}$$

Let $f_\alpha =\alpha \Psi _{j,k}$ with $\alpha \in {\mathbb {R}}$ and $\Psi _{j,k}$ as defined in (13). If $\Vert f_\alpha \Vert _{L^2}\le \nu ^2S\sqrt{j^2+k^2}$, with S as in (22), then (5) and (6) admits a unique solution $(u,p)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$ which is explicitly given by

$$\begin{aligned} u={\alpha \over \nu (j^2+k^2)}\Phi _{j,k}\, ,\qquad p={\alpha ^2\over \nu ^2(j^2+k^2)^2} \left( {2 \over \pi ^2} - \Upsilon _{j,k}\right) \, . \end{aligned}$$

(30)

Moreover, (30) is the “largest solution”, that is: if $\Vert f_\alpha \Vert _{L^2}=|\alpha |>\nu ^2S\sqrt{j^2+k^2}$ and (5) and (6) admits another solution (v, q), different from (30), then

$$\begin{aligned} \nu S\le \Vert \nabla v\Vert _{L^2}<\Vert \nabla u\Vert _{L^2}\, ,\qquad \int _\Omega (v\cdot \nabla )v\cdot u>0\, ,\qquad v_2\Delta v_1\not \equiv v_1\Delta v_2 \text{ in } \Omega \, . \end{aligned}$$

(31)

Proof

We first remark that

$$\begin{aligned} (\Phi _{j,k}\cdot \nabla )\Phi _{j,k}=\nabla \Upsilon _{j,k} \end{aligned}$$

(32)

by explicit computation, so that $(\Phi _{j,k}\cdot \nabla )\Phi _{j,k}\in G_2$, see (9). Therefore, (30) is indeed a smooth solution of (5) and (6) and the uniqueness statement for small $|\alpha |$ is a straightforward consequence of a generalized version of (29).

Assume now that, for some $|\alpha |>\nu ^2S\sqrt{j^2+k^2}$, there exists another solution (v, q) of (5) and (6). The lower bound for $\Vert \nabla v\Vert _{L^2}$ in (31) follows from Proposition 4 since, otherwise, v would be the unique solution (5).

In view of (30) and (32) we know that $(u\cdot \nabla )u\in G_2$. Therefore, if we test the Eq. (11) satisfied by u with v, we obtain (recall $f_\alpha =\alpha \Phi _{j,k}$)

$$\begin{aligned} \nu \int _\Omega \nabla u:\nabla v=\alpha \int _\Omega \Phi _{j,k}v\, . \end{aligned}$$

Moreover, if we test (11) satisfied by v with v itself and we recall (24), we get

$$\begin{aligned} \nu \int _\Omega |\nabla v|^2=\alpha \int _\Omega \Phi _{j,k}v\, . \end{aligned}$$

By combining these two identities, we infer that

$$\begin{aligned} \int _\Omega |\nabla v|^2=\int _\Omega \nabla u:\nabla v\, . \end{aligned}$$

(33)

In turn, (33) yields

$$\begin{aligned} \int _\Omega |\nabla u|^2-\int _\Omega |\nabla v|^2= \int _\Omega |\nabla u|^2-2\int _\Omega \nabla v\!:\!\nabla u+\int _\Omega |\nabla v|^2=\int _\Omega |\nabla (u\!-\!v)|^2>0 \end{aligned}$$

since $v\ne u$. This proves the upper bound for $\Vert \nabla v\Vert _{L^2}$ in (31).

From (30) we infer that

$$\begin{aligned} \Vert \nabla u\Vert _{L^2}^2={\alpha ^2 \over \nu ^2(j^2+k^2)}\, . \end{aligned}$$

Hence, if we test (11), satisfied by v, with u, then we obtain

$$\begin{aligned} \nu \int _\Omega \nabla v\!:\!\nabla u+\int _\Omega (v\cdot \nabla )v\cdot u=\alpha \int _\Omega \Phi _{j,k}u={\alpha ^2 \over \nu (j^2+k^2)} =\nu \int _\Omega |\nabla u|^2\, , \end{aligned}$$

where we used the explicit form of u in (30) and the fact that $\Vert \Phi _{j,k}\Vert _{L^2}=1$. By using again (33), we then obtain

$$\begin{aligned} {1\over \nu } \int _\Omega (v\cdot \nabla )v\cdot u=\int _\Omega |\nabla u|^2-\int _\Omega \nabla v:\nabla u= \int _\Omega |\nabla (u-v)|^2>0 \end{aligned}$$

since $v\ne u$. This proves the second property in (31).

Finally, assume for contradiction that $v_2\Delta v_1\equiv v_1\Delta v_2$ in $\Omega$. Since $\nabla \cdot v=0$, we then have that

$$\begin{aligned} \begin{array}{c} v_2(\partial _{xx}v_1+\partial _{yy}v_1)=v_1(\partial _{xx}v_2+\partial _{yy}v_2)\ \Longleftrightarrow \ v_1\partial _{xy}v_1+v_2\partial _{yy}v_1=v_1\partial _{xx}v_2+v_2\partial _{xy}v_2\\ \Longleftrightarrow \ \partial _yv_1\partial _xv_1+v_1\partial _{xy}v_1+\partial _yv_1\partial _yv_2+v_2\partial _{yy}v_1= \partial _xv_1\partial _xv_2+v_1\partial _{xx}v_2+\partial _xv_2\partial _yv_2+v_2\partial _{xy}v_2\\ \Longleftrightarrow \ \partial _y\big (v_1\partial _xv_1+v_2\partial _yv_1\big )=\partial _x\big (v_1\partial _xv_2+v_2\partial _yv_2\big )\ \Longleftrightarrow \ (v\cdot \nabla )v\in G_2 \end{array} \end{aligned}$$

so that $\int _\Omega (v\cdot \nabla )v\cdot u=0$ by (10), contradicting the second property in (31).$\square$

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 $\Phi _{j,k}$ satisfies (7) and the solution (u, p) in (30) satisfies (8). If j is even, then $\Phi _{j,k}$ satisfies (16) and the solution (u, p) in (30) satisfies (17). Finally, if both j and k are even, then $\Phi _{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

$$\begin{aligned} -\nu \Delta w+(u\cdot \nabla )w+(w\cdot \nabla )u=0\, ,\qquad \nabla \cdot w=0\qquad \text{ in } \Omega \, . \end{aligned}$$

These equations, and their weak formulation, are naturally associated to the bilinear form

$$\begin{aligned} a_u(w,\varphi ):=\nu \int _\Omega \nabla w:\nabla \varphi +\int _\Omega (u\cdot \nabla )w\cdot \varphi +\int _\Omega (w\cdot \nabla )u\cdot \varphi \end{aligned}$$

over $U\times U$; this form satisfies

$$\begin{aligned} \big |a_u(w,\varphi )\big |\le C_u\Vert \nabla w\Vert _{L^2}\Vert \nabla \varphi \Vert _{L^2}\, ,\qquad a_u(w,w)=\nu \Vert \nabla w\Vert _{L^2}^2+\int _\Omega (w\cdot \nabla )u\cdot w\qquad \forall (w,\varphi )\in U^2\, . \end{aligned}$$

The first inequality shows that $a_u(\cdot ,\cdot )$ 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 $|\alpha |$, the Lax-Milgram Theorem cannot be applied to exclude bifurcations. Clearly, for small $|\alpha |$ 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 $\Vert f_\alpha \Vert _{L^2}=|\alpha |>\nu ^2S\sqrt{j^2+k^2}$ in the case $f_\alpha =\alpha \Phi _{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\in {\mathbb {N}}$ and $\alpha \in {\mathbb {R}}$, (5) with $f_\alpha =\alpha \Phi _{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.

Proposition 5

There exists a dense open set ${\mathcal {O}}\subset {\mathbb {R}}\times L^2(\Omega )$ such that if $(\nu ,f)\in {\mathcal {O}}$ then the problem (5) – (6) admits a finite number of solutions $(u,p)\in (U\cap H^2(\Omega ))\times H^1(\Omega )$.

4 Proof of Theorem 2

4.1 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 $\beta =\nu ^{-2}$ and scale (u, p) to $(\nu u,\nu ^2 u)$, so that (5) becomes

$$\begin{aligned} \left\{ \begin{array}{l} -\Delta u+(u\cdot \nabla )u+\nabla p=\beta f\quad \text{ in } \Omega ,\\ \nabla \cdot u=0\quad \text{ in } \Omega ,\\ u_1=\partial _xu_2=0\quad \text{ on } \{0,\pi \}\times (0,\pi ),\\ u_2=\partial _yu_1=0\quad \text{ on } (0,\pi )\times \{0,\pi \}. \end{array}\right. \end{aligned}$$

(34)

We look for solutions in the space ${\mathcal {C}}_\rho$ defined in (19), with $\rho ={33\over 32}$. For $u=\bigl (u_1,u_2\bigr )\in {\mathcal {C}}_\rho$ we write

$$\begin{aligned} u=\sum _{i,j,k}u_{ijk}\Psi _{i,j,k}, \end{aligned}$$

(35)

where the vector fields $\Psi _{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

$$\begin{aligned} {{\mathbb {P}}}\bigl (c_1\Psi _{1,j,k}+c_2\Psi _{2,j,k}\bigr ) ={(c_1k-c_2j)k\over j^2+k^2}\Psi _{1,j,k} -{(c_1k-c_2j)j\over j^2+k^2}\Psi _{2,j,k}. \end{aligned}$$

(36)

Applying ${{\mathbb {P}}}$ on both sides of the first equation in (34) yields

$$\begin{aligned} {{\mathbb {P}}}(-\Delta )u ={{\mathbb {P}}}\bigl [\beta f-(u\cdot \nabla )u\bigr ], \end{aligned}$$

while the remaining equations in (34) show that $u={{\mathbb {P}}}u$. As it happens for the variational characterization (11), the pressure p does not appear in the equation, when projected onto $G_1$:

Lemma 1

Let $f\in {\mathcal {C}}_\rho$. For $u\in {\mathcal {C}}_\rho$ define

(37)

Equation (34) admits a solution $u\in {\mathcal {C}}_\rho$ if and only if $u={\mathbb {F}}_\beta (u)$.

Proof

Since functions in ${\mathcal {C}}_\rho$ satisfy the boundary conditions in (34), and $\Omega$ is simply connected, then $u\in {\mathcal {C}}_\rho$ is a solution of (34) if and only if both u and $f-(u\cdot \nabla )u$ are divergence free. When restricted to ${\mathcal {C}}_\rho$, the Laplacian commutes with the projection ${{\mathbb {P}}}$, therefore if $u={\mathbb {F}}_\beta (u)$, then u is divergence free and so is $f-(u\cdot \nabla )u$.$\square$

We take $f=(f_1,0)$, where

$$\begin{aligned} f_1(x,y)=\sum _{j,k=1,\ldots ,11}f_{jk}\Psi _{1jk}(x,y), \end{aligned}$$

with the coefficients $f_{jk}$ being approximations of the coefficients ${\hat{f}}_{jk}$ of the expansion of a Gaussian centered in $(\pi /2,\pi /2)$, that is

$$\begin{aligned} \sum _{j,k\ge 1}{\hat{f}}_{jk}\Psi _{1jk}(x,y)=-e^{-30((x-\pi /2)^2+(y-\pi /2)^2)}. \end{aligned}$$

For the exact values of the coefficients $f_{jk}$ we refer to the file data/f in [5], where they are stored in hexadecimal format. For comparison with the uniqueness result presented in Proposition 4, we note that $\Vert f\Vert _{L^2}=0.1578\ldots$

By recalling the change of parameters leading to (34), Theorem 2 will be proved once we prove the following statement.

Theorem 5

For all $\beta \in [0,2500]$ Eq. (34), with f as described above, admits a solution $u_{\beta }\in {\mathcal {C}}_\rho$ symmetric with respect to the reflection $y\mapsto \pi /2-y$, see (8). The function $\beta \mapsto u_\beta$ is analytic. The solution $u_\beta$ is isolated within the subspace of symmetric functions satisfying (8). At some $\beta \in [2065+1485\cdot 2^{-11}, 2065+1487\cdot 2^{-11}]$ there is a pitchfork bifurcation, and the secondary branches are not symmetric. Equation (34) also admits non-symmetric solutions at $\beta =2100$, $\beta =2200$ and $\beta =2500$. The solution at the bifurcation point satisfies $\Vert \nabla u\Vert _{L^2}=66.83\ldots$

Remark 3

As observed in Theorem 4, when the forcing term is proportional to $\Phi _{j,k}$ and sufficiently small, there exists a unique solution also proportional to $\Phi _{j,k}$ and, for such solution, the quadratic term is a gradient, therefore its projection on $G_1$ is zero. Then the eigenvalues of $D{\mathbb {F}}_\beta$ are proportional to $\beta$. We have numerical evidence that such eigenvalues lie all on the imaginary axis for some arbitrary choice of $\beta$, therefore no eigenvalue can reach the value 1 for any value of $\beta$, and since that is a necessary condition to have a bifurcation we conjecture that no bifurcation takes place.

4.2 Branches

Denote by ${\mathcal {D}}_\rho$ the subspace of ${\mathcal {C}}_\rho$ characterized by the symmetry (7). Note that ${\mathbb {F}}_\beta ({\mathcal {D}}_\rho )\subset {\mathcal {D}}_\rho$ and ${{\mathbb {P}}}f\in {\mathcal {D}}_\rho$. To prove that (34) admits a solution $u(\beta )\in {\mathcal {D}}_\rho$ for each $\beta \in [0,2500]$, and that the function $\beta \mapsto u(\beta )$ is analytic, we write all the coefficients in the Fourier expansion of u as Taylor polynomials in $\beta$:

$$\begin{aligned} u(\beta )=\sum _{i,j,k} a_{ijk}(\beta )\Psi _{i,j,k},\quad a_{ijk}(\beta )=\sum _{l=0}^L a_{ijkl}\left( \beta -b\over \delta \right) ^l, \end{aligned}$$

(38)

where $a_{ijkl}=0$ if k is odd and b, $\delta$, L are as in Table 1. As a first step, we choose some Fourier–Taylor polynomial ${\bar{u}}$ that is an approximate fixed point of ${\mathbb {F}}_{\beta }$, and some finite rank operator M such that ${{\mathbb {I}}}-M$ is an approximate inverse of ${{\mathbb {I}}}-D{\mathbb {F}}_{\beta }({\bar{u}})$. Then for $h\in {\mathcal {D}}_\rho$ we define

$$\begin{aligned} {{\mathcal {N}}}_\beta (h)={\mathbb {F}}_\beta ({\bar{u}}+ \Lambda h)-{\bar{u}}+M h, \qquad \Lambda ={{\mathbb {I}}}-M. \end{aligned}$$

(39)

Clearly, if h is a fixed point of ${{\mathcal {N}}}_\beta$, then $u=\bar{u}+\Lambda h$ is a fixed point of ${\mathbb {F}}_\beta$ and, hence, u solves (34) in view of Lemma 1. Given $r>0$ and $w\in {\mathcal {C}}_\rho$, let $B_r(w)=\{v\in {\mathcal {C}}_\rho \,:\,\Vert v-w\Vert _\rho <r\}$. We partition the interval [0, 2500] into six subintervals. The center b and width $\delta$ for each subinterval are shown in Table 1.

Table 1 Branch intervals

Then we prove the following lemma with the aid of a computer, see Sect. 5.

Lemma 2

The following holds for each $i\in \{1,\ldots ,6\}$. Define b, $\delta$, and L as in row i of Table 1. There exist a Fourier–Taylor polynomial ${\bar{u}}(\beta )$ of degree L as described in (38), a bounded linear operator $M(\beta )$ on ${\mathcal {D}}_\rho$, and positive real numbers ${\varepsilon }, r, K$ satisfying ${\varepsilon }+Kr<r$, such that

$$\begin{aligned} \Vert {{\mathcal {N}}}_\beta (0)\Vert \le {\varepsilon },\quad \Vert D{{\mathcal {N}}}_\beta (v)\Vert \le K,\quad \forall v\in B_r(0) \end{aligned}$$

(40)

holds for all $\{\beta \in {\mathbb {C}}: |\beta -b|<\delta \}$. Furthermore, if $i<6$, then either

$$\begin{aligned} B_{r^i}(u(\delta ^i+b^i))\subset B_{r^{i+1}}(u(\delta ^{i+1}-b^{i+1})), \text { or } B_{r^{i+1}}(u(\delta ^{i+1}-b^{i+1}))\subset B_{r^i}(u(\delta ^i+b^i)), \end{aligned}$$

(41)

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 $\beta \in [0,2500]$ there exists a fixed point $u_\beta \in B_r(u_\beta )$ of ${\mathbb {F}}_\beta$, and the curve $\beta \mapsto u_\beta$ is analytic.

4.3 Bifurcations

Let us write the fixed point equation for ${\mathbb {F}}_\beta$ as ${\mathbb {F}}(\beta ,u)=0$, where

(42)

Numerically, we observe that $D{\mathbb {F}}(\beta ,\mathbf{.})$ has a (simple) eigenvalue zero for $\beta$ near $b=2065+743\cdot 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 $\beta$ and a coordinate $\lambda$ for the range of a suitable one-dimensional projection $\ell$. Then we define a two-parameter family of functions $u(\alpha ,\lambda )$ by solving

$$\begin{aligned} ({{\mathbb {I}}}-\ell ){\mathbb {F}}\bigl (\beta ,u(\beta ,\lambda )\bigr )=0,\qquad \ell u(\beta ,\lambda )=\lambda {\hat{u}}, \end{aligned}$$

(43)

where ${\hat{u}}$ is a fixed nonzero function in the range of $\ell$. For ${\hat{u}}$ we choose a Fourier polynomial that approximates the eigenvector of $D{\mathbb {F}}(b,\mathbf{.})$ for the eigenvalue closest to zero. The projection $\ell$ is defined by

$$\begin{aligned} \ell u=\ell _0(u){\hat{u}},\qquad \ell _0(u)=\sum _{i,j,k=1,2}u_{ijk}{\hat{u}}_{ijk}, \end{aligned}$$

(44)

where $u_{ijk},{\hat{u}}_{ijk}$ are the Fourier coefficients of u and ${\hat{u}}$, respectively. Our goal is to show that for a rectangle $I\times J$ in the parameter space, the Eq. (43) has a smooth and locally unique solution $u: I\times J\rightarrow {\mathcal {C}}_{\rho }\,$. Then locally, the solutions of ${\mathbb {F}}(\beta ,u)=0$ are determined by the zeros of the function g,

$$\begin{aligned} g(\beta ,\lambda ){\hat{u}} =\ell {\mathbb {F}}\bigl (\beta ,u(\beta ,\lambda )\bigr ). \end{aligned}$$

(45)

We write all the coefficients in the Fourier expansion of u as Taylor polynomials in $\beta ,\lambda$:

$$\begin{aligned} u(\beta ,\lambda )=\sum _{i,j,k} a_{ijk}\Psi _{1jk},\quad a_{ijk}=\sum _{l+m\le 4} a_{ijklm}\left( \beta -b\over \delta \right) ^l \left( \lambda \over \lambda _1\right) ^m, \end{aligned}$$

(46)

where $\delta =2^{-11}$ and $\lambda _1=2^{-6}$.

Equation (43) for $u=u(\beta ,\lambda )$ is equivalent to the fixed point equation for the map ${\mathbb {F}}_{\beta ,\lambda }\,$, defined by

$$\begin{aligned} {\mathbb {F}}_{\beta ,\lambda }(u)=({{\mathbb {I}}}-\ell ){\mathbb {F}}_\beta (u)+\lambda {\hat{u}}. \end{aligned}$$

(47)

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 ${\mathcal {M}}_{\beta ,\lambda }$ defined by

$$\begin{aligned} {\mathcal {M}}_{\beta ,\lambda }(h)={\mathbb {F}}_{\beta ,\lambda }({\bar{u}}+ \Lambda h)-\bar{u}+M h,\qquad \Lambda ={{\mathbb {I}}}-M. \end{aligned}$$

(48)

Here, ${\bar{u}}$ is an approximate fixed point of ${\mathbb {F}}_{b,0}$, and M is a finite rank operator such that $\Lambda ={{\mathbb {I}}}-M$ is an approximate inverse of ${{\mathbb {I}}}-D{\mathbb {F}}_{b,0}({\bar{u}})$. By the Implicit Function Theorem, the solution then depends analytically on the two parameters $\beta$ and $\lambda$. Denote by $D_r(z)$ a disk in the complex plane of radius r and center z. Let $I=D_{\delta }(b)$ and $J=D_{\lambda _1}(0)$. The following lemma is proved with the aid of a computer; see Sect. 5.

Lemma 3

There exists a Fourier polynomial ${\hat{u}}$, a Fourier–Taylor polynomial ${\bar{u}}(\beta ,\lambda )$ as in (46), and positive constants ${\varepsilon },r,K$ satisfying ${\varepsilon }+Kr<r$, such that

$$\begin{aligned} \Vert {\mathcal {M}}_{\beta ,\lambda }(0)\Vert \le {\varepsilon }, \qquad \Vert D{\mathcal {M}}_{\beta ,\lambda }(v)\Vert \le K. \end{aligned}$$

(49)

for all $v\in B_{r}(0)$ and for all $\beta \in I$ and all $\lambda \in J$.

As we had in the previous subsection, this lemma, together with the contraction mapping theorem and the implicit function theorem, implies the following:

Proposition 7

For every $(\beta ,\lambda )$ in $I\times J$, the Eq. (43) has a unique solution $u(\beta ,\lambda )$ in $B_r({\bar{u}}(\beta ,\lambda ))$, and the map $(\beta ,\lambda )\mapsto u=u(\beta ,\lambda )$ is analytic. For any given real $\beta \in I$, a function u in $B\cap \ell ^{-1}(J{\hat{u}})$ is a fixed point of $F_\beta$ if and only if $u=u(\beta ,\lambda )$ for some real $\lambda \in J$, and $g(\beta ,\lambda )=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 $\dot{f}$ and $f'$ the partial derivatives of f with respect to the first and second argument, respectively. Let $I=[\beta _1,\beta _2]$ and $J=[-b,b]$.

Lemma 4

(pitchfork bifurcation) Let g be a real-valued $C^3$ function on an open neighborhood of $I\times J$, such that $g(\beta ,0)=0$ for all $\beta \in I$, and

Then $g(\beta ,\lambda )=\lambda G(\beta ,\lambda )$ for some $C^2$ function G, and the solution set of $G(\beta ,\lambda )=0$ in $I\times J$ is the graph of a $C^2$ function $\beta =a(\lambda )$, defined on a proper subinterval $J_0$ of J. This function takes the value $\beta _2$ at the endpoints of $J_0\,$, and satisfies $\beta _1<a(z_2)<\beta _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.

Lemma 5

For any $u(\beta ,\lambda )\in B_r({\bar{u}}(\beta ,\lambda ))$ the function satisfies the assumptions of Lemma 4.

The proof of Theorem 2 is concluded by joining the results of Propositions 6 and 7 , together with Lemmas 4 and 5 .

5 Computer Estimates

5.1 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,\delta )$ from a fixed but arbitrary row in Table 1. As a first step, we determine a Fourier–Taylor polynomial ${\bar{u}}_\beta$ which is an approximate fixed point of ${\mathbb {F}}_\beta$ and an approximate inverse of ${{\mathbb {I}}}-D{\mathbb {F}}_{\beta }({\bar{u}}_\beta )$ of the form $\Lambda (\beta )={{\mathbb {I}}}-M(\beta )$, where $M(\beta )$ has finite rank and is a polynomial in $\beta$. The remaining steps are rigorous: we compute an upper bound ${\varepsilon }$ on the norm of ${\mathbb {F}}_\beta (0)$, and an upper bound K on the operator norm of $D{\mathbb {F}}_\beta (h)$ that holds for all h of norm ${\varepsilon }$ or less. This is done simultaneously for all values of $\{\beta \in {\mathbb {C}}: |\beta -b|<\delta \}$. After verifying that $K<\frac{7}{8}$, we choose a positive $r<8{\varepsilon }$ in such a way that ${\varepsilon }+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+\beta _1z_1,\gamma _0+\gamma _1z_2)$, and note that h is analytic in $\{z\in {\mathbb {C}}^2\,:\,|z_1|<1,|z_2|<1\}$ by Proposition 6. Given the Fourier–Taylor polynomial ${\bar{u}}(\beta ,\lambda )$ of Lemma 3, we compute a Taylor polynomial $P(z_1,z_2)$ of degree 4 and a bound E, such that that

$$\begin{aligned} \sup _{|z_1|<1,|z_2|<1}|P(z_1,z_2)-h(z_1,z_2)|<E. \end{aligned}$$

(50)

Derivatives of analytic functions such as $f=(P-h)(.,z_2)$ or $f=(P-h)(z_1,.)$ can be estimated by using the standard Cauchy bound: if f(z) is an analytic function in the unit disk such that $\sup _{|z|\le 1}|f(z)|<E$, then for $0<\varrho <1$ and $n>0$

$$\begin{aligned} \sup _{|z|\le \varrho }|f^{(n)}(z)|\le {n!\over (1-\varrho )^n}E. \end{aligned}$$

(51)

Using $\varrho =205/256$, this bound is applied (separately in each variable) to verify the conditions in Lemma 4.

5.2 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:{{\mathcal {X}}}\rightarrow {{\mathcal {Y}}}$ is a function F that assigns to a set $X\subset {{\mathcal {X}}}$ of a given type (Xtype) a set $Y\subset {{\mathcal {Y}}}$ of a given type (Ytype), in such a way that $y=f(x)$ belongs to Y for all $x\in 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 ${{\mathcal {X}}}$ with unit $\mathbf{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 ${{\mathcal {X}}}$ is $\langle \mathtt{S},{{\mathcal {X}}}\rangle =\{x\in {{\mathcal {X}}}:\Vert x-(\mathtt{S.C})\mathbf{1}\Vert \le \mathtt{S.R}\}$.

When ${{\mathcal {X}}}={\mathbb {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 $\rho$, consider the space ${\mathbb {T}}_\rho$ of all real analytic functions $g(t)=\sum _n g_nt^n$ on the interval $|t|<\rho$, obtained by completing the space of polynomials with respect to the norm $\Vert g\Vert _\rho =\sum _n|g_n|\rho ^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, $\mathtt{P.R}=\rho$, and P.C is an array(0..D) of Ball. The corresponding set in $\langle \mathtt{Taylor1},{\mathbb {T}}_\rho \rangle$ is defined as

$$\begin{aligned} \langle \mathtt{P},{\mathbb {T}}_\rho \rangle =\sum _{n=0}^{m-1}\bigl \langle \mathtt{P.C(n)},{\mathbb {R}}\bigr \rangle p_n +\sum _{n=m}^D\bigl \langle \mathtt{P.C(n)},{\mathbb {T}}_\rho \bigr \rangle p_n,\qquad p_n(t)=t^n, \end{aligned}$$

(52)

where $m=\min (\mathtt{P.F},\mathtt{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 ${{\mathcal {X}}}$. Consider the space ${\mathbb {F}}_\rho$ of all real analytic functions in two variables periodic of period $2\pi$ in both variables, obtained by completing the space of Fourier polynomials $p(w,z)=\sum _{jk}p_{jk}e^{i(jw+kz)}$, $p_{jk}\in {{\mathcal {X}}}$, with respect to the norm $\Vert p\Vert _\rho =\sum _{jk}\Vert p_{jk}\Vert \rho ^{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 $\langle \mathtt{F},{\mathbb {F}}_\rho \rangle$ is the set of all function $u=p+h\in {\mathbb {F}}_\rho$, where

$$\begin{aligned} p(w,z)=\!\!\sum _{j,k=1}^K \langle \mathtt{F.C(J,K)},{{\mathcal {X}}}\rangle \,e^{i(jw+kz)},\quad h=\sum _{j,k=1}^{2K}h^{j,k},\quad h^{j,k}(w,z)=\!\!\!\sum _{m\ge j, n\ge k}\!\! h^{j,k}_{m,n}\,e^{i(mw+nz)}, \end{aligned}$$

(53)

with $p_{J,K}\in \langle \mathtt{F.C(J,K)},{{\mathcal {X}}}\rangle$ and $\Vert h^{J,K}\Vert \le \mathtt{F.E(J,K)}$, for all $J,K\ge 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 ${\mathcal {C}}_\rho$ and ${\mathcal {D}}_\rho$. 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 ${\mathcal {C}}_\rho$, a bound NegInvLap on $-\Delta ^{-1}$, a bound ProjDivFree on the projection onto divergence-free vector fields, a bound NonLinearAdv on the map $u\mapsto f-u\cdot \nabla 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 ${\mathbb {F}}_\beta$ defined in (37) is implemented by the procedure GMap in the package Taylors1.Foor.Fix. Defining and estimating a contraction like ${{\mathcal {N}}}_\beta$ 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 ${\mathcal {C}}_{\beta ,\lambda }$ defined in (48). The bound GMap on the map ${\mathbb {F}}_{\beta ,\lambda }$ defined in (47) is implemented in the package Taylors2.Foor.Fix.

6 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 $\Omega$ 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.

Problem 1

Classical arguments in critical point theory show that a (scalar) minimizer $w\in H^1_{\texttt {V}}$ of the ratio $\Vert \nabla w\Vert _{L^2}^2/\Vert w\Vert _{L^4}^2$ satisfies the Euler–Lagrange equation

$$\begin{aligned} -\Delta w=w^3 \text{ in } \Omega \, ,\quad w=0 \text{ on } V\, ,\quad \partial _yw=0 \text{ on } H\, . \end{aligned}$$

If we take $w_0\in H^1_{\texttt {V}}$ independent of y we find a boundary value problem for an ODE, that is,

$$\begin{aligned} w_0''(x)+w_0(x)^3=0 \text{ in } (0,\pi )\, ,\quad w_0(0)=w_0(\pi )=0\, , \end{aligned}$$

(54)

whose solution is given by

$$\begin{aligned} w_0(x)={T \over 2\pi }\, \text{ cn }\left( {T \over 2\pi }\left( x-{\pi \over 2}\right) ,{1\over \sqrt{2}}\right) \, ,\qquad T=4\int _0^{\pi /2}{dt \over \sqrt{1-{\sin ^2t \over 2}}}\approx 7.416\, , \end{aligned}$$

where cn denotes the Jacobi cosine, see [1]. Multiplying (54) by $w_0$ and integrating by parts yields

$$\begin{aligned} \int _0^\pi w_0'(x)^2\, dx=\int _0^\pi w_0(x)^4\, dx={T^4 \over 16\pi ^4}\int _0^\pi \text{ cn }\left( {T \over 2\pi }\left( x\!-\!{\pi \over 2}\right) ,{1\over \sqrt{2}}\right) ^4\, dx ={16 \over 3\pi ^3 }\left[ \int _0^{\pi /2}{dt \over \sqrt{1\!-\!{\sin ^2t \over 2}}}\right] ^4\, . \end{aligned}$$

Therefore, if we view $w_0$ as a function of two variables, we obtain the following bound for S in (22):

$$\begin{aligned} S=\ \min _{w\in H^1_{\texttt {V}}}\ {\Vert \nabla w\Vert _{L^2}^2 \over \Vert w\Vert _{L^4}^2}\le {\Vert \nabla w_0\Vert _{L^2}^2 \over \Vert w_0\Vert _{L^4}^2}= {\int _0^\pi \int _0^\pi w_0'(x)^2\, dx\, dy \over (\int _0^\pi \int _0^\pi w_0(x)^4\, dx\, dy)^{1/2}}=\sqrt{\pi }\left[ \int _0^\pi \!\! w_0(x)^4 dx\right] ^{1/2}\approx 1.4257\, . \end{aligned}$$

(55)

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 $\nu _0=\nu _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 $\Omega$, and to more general forces f? Can the same technique be extended to 3D domains?

• Is it possible to find $g\in L^2(\Omega )$ such that Eq. (5) with $f=\alpha g$ admits a unique solution for all $\alpha \in [0,+\infty )$? 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 $\Omega$, 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.