A continuously perturbed Dirichlet energy with area-preserving stationary points that ‘buckle’ and occur in equal-energy pairs

We exhibit a family of convex functionals with infinitely many equal-energy C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document} stationary points that (i) occur in pairs v±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{\pm }$$\end{document} satisfying det∇v±=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \nabla v_{\pm }=1$$\end{document} on the unit ball B in R2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^2$$\end{document} and (ii) obey the boundary condition v±=id\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{\pm }=\text {id}$$\end{document} on ∂B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \partial B$$\end{document}. When the parameter ϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} upon which the family of functionals depends exceeds 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{2}$$\end{document}, the stationary points appear to ‘buckle’ near the centre of B and their energies increase monotonically with the amount of buckling to which B is subjected. We also find Lagrange multipliers associated with the maps v±(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{\pm }(x)$$\end{document} and prove that they are proportional to (ϵ-1/ϵ)ln|x|\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\epsilon -1/\epsilon )\ln |x|$$\end{document} as x→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x \rightarrow 0$$\end{document} in B. The lowest-energy pairs v±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{\pm }$$\end{document} are energy minimizers within the class of twist maps (see Taheri in Topol Methods Nonlinear Anal 33(1):179–204, 2009 or Sivaloganathan and Spector in Arch Ration Mech Anal 196:363–394, 2010), which, for each 0≤r≤1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\le r\le 1$$\end{document}, take the circle {x∈B:|x|=r}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{x\in B: \ |x|=r\}$$\end{document} to itself; a fortiori, all v±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{\pm }$$\end{document} are stationary in the class of W1,2(B;R2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W^{1,2}(B;{\mathbb {R}}^2)$$\end{document} maps w obeying w=id\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w=\text {id}$$\end{document} on ∂B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial B$$\end{document} and det∇w=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \nabla w=1$$\end{document} in B.


Introduction
Let > 1 be a real parameter, letx = x/|x| for any non-zero x in R 2 and let B be the unit ball in R 2 . We study the functional D given by and our goal is to find a global minimizer of D in A.
In this direction, we prove that A contains stationary points of D which occur in pairs, v ± , say, such that (a) D (v + ) = D (v − ) and (b) v + and v − are distinguished by the way they 'buckle' near the centre 0 of the ball B, each executing a 'turn' there through ±π/2. In the incompressible setting, a stationary point v ± satisfies, for an appropriate pressure function p ± , the weak Euler-Lagrange equation (See Theorem 1.2, Sect. 3.1 for details.) The maps v ± are twist maps, i.e. they are of the form v ± (x) = Rot (±k(|x|))x, where k is a suitably chosen function and Rot (h) stands for the 2 × 2 matrix which rotates vectors by h radians anticlockwise about the origin. Thus twist maps take circles centred at the origin to themselves, only allowing for a change in the polar angle of image points. The main theorem in this paper concerning the existence of the buckling stationary points associated with the functional D is as follows: Theorem 1.1. Let D be given by (1.1) and A by (1.2). Then, for each > √ 2, there are infinitely many pairs of C 1 (B) twist maps v ±,j , j ∈ N, of the form v ±,j (x) = Rot (±k j (|x|))x, where k j ∈ C ∞ (0, 1), such that for each j:

j is a stationary point of D in A;
(iii) ∇v ±,j (0) = ±J, where J is the matrix Rot (π/2).
Before setting this in its context, we first gather together the following properties of the integrand W given in (1. (P4) W (x, F ; 1) = |F | 2 , and in fact W (x, F ; ) can be viewed as a 'perturbation' of W (x, F ; 1).
Properties (P1) and (P2) are easily verified. To see (P3), note that W (x, F ; ) is 2-homogeneous with respect to F , so its Hessian obeys for any 2 × 2 matrix A. Property (P4), by contrast, is less obvious, so we give the details here. Let J be the 2 × 2 matrix representing a rotation of π/2 radians anticlockwise, and note that for any non-zero x in R 2 we may write the identity matrix, 1, as 1 =x ⊗x + Jx ⊗ Jx. Then, by writing the inner product between two matrices F 1 and F 2 as F 1 · F 2 = tr (F T 1 F 2 ), we see that a ⊗ b · c ⊗ Jb = (a · c)(b · Jb) = 0 for any two-vectors a, b, c, and hence To pass from the second to the third line, note that by orthogonality. The fourth line follows from the third by using the fact that J T is an isometric linear map, and in reaching the final line it helps to know the easily-verified identity adj F = J T F T J. Finally, by weighting the first and second of these terms with prefactors 1 and respectively, we obtain W (x, F ; ). Theorem 1.1 points to the possibility that there are distinct but equalenergy global minimizers of the functional D (·). Recall that in the context of energy functionals of the form Ω f (∇u) dx, where f is strictly polyconvex, Ω is a domain homeomorphic to a ball and u| ∂Ω agrees with an affine map, the question of the uniqueness of minimizers, as set out in Problem 8 in [4], was settled by Spadaro in [24]. To relate that work to ours, first note that property (P3) enables us (via (1.5) below) to extend W to a polyconvex function. Thus by seeking maps v ± that agree with the identity on ∂B and which minimize D , our aim is to tackle a question that is analagous to Problem 8 in [4], but with two important distinctions. These are that (1) admissible maps v, say, are incompressible, i.e. det ∇v = 1 a.e., and (2) x−dependent integrands are permitted. For clarity, we stress that while the maps v ± which feature in this paper are certainly stationary points of D , we do not know whether they are minimizers. Nevertheless, we believe that to establish their existence and stationarity is an important intermediate step towards this goal. We remark that since Spadaro's minimizers are such that their Jacobians take different signs on different parts of the domain he considers, they do not seem to be directly relevant to our incompressible setting.
It is perhaps striking that, in the pure displacement problem we consider, the maps v ± appear to emulate the buckling beam examples frequently cited as examples of non-uniqueness of equilibria in nonlinear elasticity when mixed boundary conditions are applied 1 . See Fig. 1. Indeed, by viewing W as the finite part of the polyconvex function W given by can be interpreted as the energy stored by an incompressible material occupying the ball B in a reference configuration. By setting = 1, D 1 (v) models the stored energy of an isotropic, incompressible neo-Hookean material (see [10,Section 4.10], for example). When > 1, the connection between nonlinear elasticity and the functional D (v) is maintained, perhaps tenuously and only as a result of a choice made by the authors, by regarding W (x, F ; ) as a perturbation of |F | 2 , as described in Property (P4) earlier.
Works which tackle non-uniqueness in finite elasticity include but are not limited to [2,3,20]. 2 The problems treated are of mixed displacement-traction or pure traction type. Non-uniqueness in the context of pure displacement problems is considered in [6,16,19], and we discuss this point further below in relation to the pure displacement problem dealt with in this paper. There is also a significant body of work on the question of uniqueness of equilibria in finite elasticity which includes but is not limited to [5,16,17,22,23,29,30].
That there is a global minimizer of D in A can be established as follows. Firstly, by property (P3) and [1,Theorem II.4], it follows that the functional D is sequentially lower semicontinuous with respect to weak convergence in H 1 (B; R 2 ). Then, by applying [11,Theorem 8.20,part 1], which tells us that A is closed with respect to weak convergence in H 1 (B; R 2 ), the direct method of the calculus of variations yields the existence of a global minimizer of D in A. When > √ 2, a careful numerical calculation (given in Sect. 4) shows that the global minimizer is not the identity. 3 Figure 1. The map v + with = 5.0, applied to a sector of angle μ = 0.3. The effect of v + is to 'buckle' the ball B near its centre. Stationary points which exhibit more exotic buckling are possible when > √ 2: see Proposition 2.16 and the ensuing discussion, as well as Fig. 5 It is worth pointing out that pure displacement variational problems with infinitely many equilibria were first proposed by F. John in [16] and later confirmed rigorously in [19], with further analysis in [6]. The domain in these cases was an annulus on which the identity boundary condition was imposed; solutions are distinguished by the integer number of times the outer boundary is rotated relative to the inner boundary. By contrast, the maps v ±,j we construct here are defined on the ball B, all display a similar buckling behaviour near the origin, and it seems that a different mechanism gives rise to the countable number of solutions we find when > √ 2. We note that the role played by topology in establishing non-uniqueness is studied carefully in [27].
The idea of considering energies that are related to the Dirichlet energy can also be found in the recent work of Taheri and Morris [18], who study minimizers of functionals of the form (see [18,Equation 1 where X = X[a, b] is the annulus {x ∈ R 2 : a < |x| < b} and where the function F satisfies certain hypotheses ((H1) and (H2), in their notation). The competing functions u are incompressible and agree with the identity on ∂X, and the minimizers are twist maps (which they call whirl maps). Thus their work has significant features in common with ours, but their emphasis is on uniqueness within different homotopy classes. Moreover, our integrand W (x, F ; ) does not factorize as do theirs, and the fact that buckling takes place at x = 0 would seem to rule out the consideration of annuli X [a, b].
Let us note that maps v in A have a continuous representative (see [14,Theorem 2.3.2] or [25,Theorem 5]) and thus their topological degree, d(v, ·, ·), is classically defined. Using this and the change of variables formula [12,Theorem 5.35], it holds that for any open set G ⊂ B and any function g in L ∞ (B), Since the degree is determined by the trace, and since det ∇v = 1 a.e., we can take g = χ v(G) and, after a short calculation, thereby obtain |G| = |v(G)|. Hence maps in A are measure-preserving. They are also homeomorphisms of the ball B, which fact follows by applying the argument of [18,Theorem 2.1] with B in place of X(a, b), as defined earlier.
We conclude this introduction by describing the system of PDEs that one might expect a stationary point of D in A to satisfy. In view of the incompressibility constraint det ∇v = 1 a.e. in B obeyed by admissible maps v, we seek a so-called pressure term, p, say, such that For sufficiently smooth functions v and p, the weak form stated above is equivalent to the pointwise equation supplemented by the conditions implied by membership of v in A. However, since we do not a priori know that the functions v and p are sufficiently regular, the approach we take has to be more indirect. We postpone these details to the discussion given at the beginning of Sect. 3, and instead state the main result of this paper concerning the Lagrange multipliers, or pressures, associated with the buckling stationary points described in Theorem 1.1 above.

Theorem 1.2.
Let v ±,j be a twist map stationary point of D , as described in Theorem 1.1, of the form v ±,j (x) = Rot (±k j (|x|))x. Then the pressure function is such that In particular,

Plan of the paper
After introducing necessary notation and terminology in Sect. 1.2, the rest of the paper is organised as follows. Section 2 is broadly concerned with establishing the existence of nontrivial stationary points of the functional D in the class of twist maps defined in (2.7). Using the direct method, we confirm in Sect. 2.1 that there are indeed energy-minimizing twist maps v(x) = Rot (k(|x|))x, say, and in Sect. 2.2 we establish the relevant Euler-Lagrange equation which k must satisfy, together with appropriate natural and imposed boundary conditions. We then transform this equation into a planar, autonomous ODE and, in Sect. 2.3, identify the non-trivial solutions we seek as heteroclinics connecting certain rest points in the phase plane, making extensive use of a conservation law (2.27) that we derive along the way. The asymptotic detail of these heteroclinic orbits shows that when 1 < ≤ √ 2 there are no non-trivial solutions, while, for > √ 2, we see in Sect. 2.4 that there are countably many pairs v ±,j , say, with j ∈ N, of solutions, all of which agree with the identity map on ∂B. Their energies can be calculated and we give the numerical results in Sect. 4. Section 3 opens with a description of the strategy we use to find and characterize the pressure p associated with the buckling stationary points, and is then divided into two further sections. We show in the first of these, Sect. 3.1, that the twist maps derived in Sect. 2 are stationary points of D in the sense that We use this to exploit results of [13] and [8] to deduce the existence of an associated Lagrange multiplier belonging to the space ∩ p>1 L p (B). Section 3.2 calculates the pressure in terms of the solution k of the relevant Euler-Lagrange equation derived in Sects. 2.2, and 3 concludes by showing that the two entities-pressure and Lagrange multiplier-are connected by a change of variables involving the stationary point itself.

Notation
The inner product between two matrices A and B is given by A·B = tr (A T B), where, as usual, tr A denotes the trace of A. For any two vectors a and n in R 2 , the notation a⊗n will denote the matrix of rank one whose (i, j) entry in some chosen basis is a i n j . The cofactor of a 2 × 2 matrix A is written cof A and is calculated by cof A = J T AJ, where J is the 2×2 matrix representing a rotation anticlockwise through π/2 radians. For invertible A such that det A = 1, the expression cof A = A −T is also valid. The convenient notation for the unit vectorsx = e R (μ) and Jx = e τ (μ) will feature, and we note that in these terms the gradient of a planar map v, when it exists, can be expressed as where v R := ∂v/∂R, v τ := 1 R ∂v/∂μ and R = |x|. Our notation for Sobolev spaces and the various types of convergence therein is generally standard, and where it is not it is explained at the relevant stage. Unless stated otherwise, sequences are indexed by the natural numbers. At certain points in the paper, especially Sect. 2.3, some terminology from dynamical systems is used. We believe this is standard; if in doubt, the authoritative references [9,15] should be consulted. For any set S we denote by χ S the indicator function of S, i.e. χ S (s) = 1 if s lies in S and χ S (s) = 0 otherwise. We reserve the notation 1 for the 2 × 2 identity matrix. Finally, we note that the terms 'map' and 'function' are used interchangeably at various points in the paper, and we trust that no confusion arises as a result.
Definition 2.1. Let k belong to W 1,1 ((0, 1); R) and let x = (R cos μ, R sin μ). (2.1) When R = 0, a short calculation shows that ∇v(x) = (Rk (R)e τ (μ + k(R)) + e R (μ + k(R))) ⊗ e R (μ) + e τ (μ + k(R)) ⊗ e τ (μ), (2.2) and from this it easily follows that det ∇v(x) = 1 a.e. in B for any choice of k. By choosing k(1) = 0 in the sense of trace, we ensure that v(x) = x whenever x ∈ ∂B. Taking v as above, we have Denoting the prefactor of and noting that W (x, 1; ) = + 1 , it follows, after an integration by parts, that  (2.6) and the class C by (2.7) We seek minimizers of E (k) in C, and to do so it is natural to use the direct method of the calculus of variations.

Minimizing E over C
In order to minimize E over C using the direct method we must confirm both that C is closed and that E is sequentially lower semicontinuous with respect to a notion of weak convergence that suits the class C. This is the point of the next three lemmas, the first of which establishes a useful Poincaré type inequality which applies to members of C. In the following we shall abbreviate the space H 1 ((a, b); R) to H 1 (a, b), where a < b, and we will not relabel subsequences. (2.8) Proof. The proof begins, as usual, by supposing for a contradiction that there is a sequence k j in C such that and for which Fix 0 < δ < 1 and note that by (2.9) and (2.10)k j | (δ,1) is bounded in H 1 (δ, 1), and hence there is somek(·, δ) in H 1 (δ, 1) such that, for a subsequence,k j | (δ,1) k (·, δ). By (2.9) and sequential weak lower semicontinuity of the functional 1 δ R 3k 2 j dR we must havek (·, δ) = 0 a.e. on (δ, 1). By Sobolev embedding, we may assume thatk j →k(·, δ) in L 2 (δ, 1), and hence in particular that a subsequence satisfiesk j →k(·, δ) a.e. in (δ, 1). Now, for all j and all δ ≤ R ≤ 1 we havek j (R) = 1 Rk j (s) ds, so by letting j → ∞ it must be thatk(R, δ) = 1 Rk (s, δ) ds a.e. in (δ, 1). Both sides are continuous and so they agree everywhere, and sincek (·, δ) = 0 by the above,k(·, δ) = 0 on (δ, 1).
Thereforek j → 0 in L 2 (δ, 1), and by letting γ 2 < 1 and applying (2.10), it follows that a.e. in (0, δ) then by squaring and integrating we contradict (2.11). Hence the set {R ∈ (0, δ) : R 1 2 |k j (R)| > γδ − 1 2 } has positive measure, and, sincek j is absolutely continuous, there must be 0 < 1 2 . This forces 1 δ R 3k 2 j dR to be bounded away from zero, as we now show. Claim: Withk j as above and R 1 2 . Set ϕ(R) = R 3 2 (k j (R) − z(R)) and use convexity to deduce that To pass from the first to the second line above we have used the facts that R and let k j be a sequence in C such that ||k j || is bounded. Then there is k in C and a subsequence such that the following convergences hold: Proof. That ||·|| is a norm and not merely a seminorm is clear from Lemma 2.3 and the definition of C. Let k j be a sequence in C such that ||k j || is bounded. Let w j (R) = R 3/2 k j (R), so that w j = 3 2 R 1/2 k j + R 3/2 k j a.e. in (0, 1). Then the sequence w j is bounded in H 1 (0, 1), and so there is w in H 1 (0, 1) such that w j w in H 1 (0, 1). By the Rellich-Kondrachov compactness theorem, we may assume that, for a subsequence, w j → w strongly in L 2 (0, 1). Furthermore, since and since, by (2.8), the right-hand side is bounded uniformly in j, it must be that w j /R converges weakly to some W , say, in L 2 (0, 1). Since w j → w in L 2 (0, 1), we can assume that w j → w a.e. in (0, 1), and hence, in particular, that w j /R → w/R a.e.. Hence W = w/R a.e., and so w j /R w/R in L 2 (0, 1) and w j w in H 1 (0, 1). It is now natural to define k(R) := w(R) R 3 2 and to note that (i) k belongs to W 1,1 loc ((0, 1); R), (ii) the weak derivative k satisfies R 3/2 k = w − 3w 2R , which belongs to L 2 (0, 1), (iii) k(1) = 0, by letting j → ∞ in the representation k j (R) = 1 R k j (s) ds and then, as in Lemma 2.3, employing a continuity argument, and (iv), by (2.8), R 1/2 k belongs to L 2 (0, 1). Thus k lies in C and the stated convergences hold by translating information about w j into information about k j .

Lemma 2.5. Let E and C be as in Definition 2.2. Then there is a minimizer of E in C.
Proof. By inspection, E is bounded below, so we may assume that there is a sequence k j in C such that E (k j ) → inf C E where the infimum is finite. The uniform bound 1/ ≤ Ω(k) ≤ in particular implies that ||k j || is bounded, so we may assume without loss of generality that there is k in C such that the convergences described in (2.13), (2.14) and (2.15) hold. The last of these together with the bounded convergence theorem easily implies that the second term in E (k j ) obeys −2p as j → ∞. To deal with the first term in E (k j ), it is helpful to reuse some of the notation introduced in Lemma 2.4, namely let w j := R 3/2 k j , and recall that a.e., with a similar expression for R 3/2 k . Then, by convexity, The rightmost integrand in the line above is then Ω(k j )fg j , which we write as Now, by (2.13) and (2.14), g j 0 in L 2 (0, 1), so (2.15) and Fatou's lemma, the first term on the right-hand side of (2.17) obeys lim inf

Minimizers and the Euler-Lagrange equation
Recall that (2.20) Let I ⊂ (0, 1) be compact and let ϕ be a smooth, real-valued function whose support is contained in I. Then if k is a minimizer of E in C, a standard NoDEA Area-preserving, buckling stationary points Page 13 of 50 6 argument implies that the weak Euler-Lagrange equation holds, or, in integrated form (via DuBois Reymond's lemma [7, Lemma 1.8]), for some constant c. One can now adapt the proof of [7, Theorem 4.1] to conclude that k is smooth on (0, 1) by (i) noting that the ellipticity condition (2.20) implies that F vv (R, z, v) > 0 for all R ∈ I, (ii) applying the (local) argument of [7, Proposition 4.4] on I and (iii) concluding, by the implicit function theorem argument given in [7,Proposition 4.4], that k is smooth on I. Since the interval I ⊂ (0, 1) is arbitrary, we reach the following conclusion.
Proposition 2.6. Let k in C minimize E given by (2.6). Then k is smooth in (0, 1) and Moreover, the natural boundary condition holds.
Proof. The smoothness of k has already been shown, and equation (2.22) follows easily from this and (2.21). Condition (2.23) is derived by taking variations k + hϕ of k in C, where h is a real parameter and spt ϕ is compactly contained in [0, 1) and is, in particular, allowed to include 0. In such circumstances the first variation of E about k in the direction ϕ, δE (k)[ϕ], both vanishes and obeys where, for brevity, we have reused the notation , the integral term vanishes, and, since ϕ(0) is arbitrary, (2.23) follows.

Phase plane analysis of the transformed Euler-Lagrange equation
Let > 1, let k ∈ C be a solution of the Euler-Lagrange equation (2.22) subject to (2.23) and let Ω(k) be given by (2.3). The change of variables R = e t and x(t) := k(R) transforms (2.22) to the autonomous equation and let X = (x, y) T . Then (2.24) can be written as the autonomous systeṁ subject to x(0) = 0 and e 2t y(t) → 0 as t → −∞. For convenience, let us refer to the variable t as time.
Notice that the change of variables X → −X is a symmetry, and that rest points lie at X k := (kπ/2, 0) T for each k ∈ Z. A standard linearization analysis, the precise results of which we record later, shows that each point X 2k+1 for k ∈ Z is a hyperbolic saddle while X 2k is a hyperbolic sink. The phase portrait of (2.26) reveals that the solutions we seek must, in view of the slightly unusual boundary condition 'at −∞', lie on a heteroclinic orbit connecting X −1 to X 0 and, also owing to the symmetry mentioned before, on a heteroclinic orbit connecting X 1 and X 0 . Before describing these orbits and the phase portrait more generally, we first take note of the following conservation law associated with (2.26).
Proof. Since the function F in the system (2.26) is smooth in X, according to standard theory for autonomous ODEs the solution is also (locally) smooth as a function of t. Suppose that y = 0 and divide the equationẏ = F 2 (X), given by the second row of (2.26), by the equationẋ = y. This gives the expression From this it follows that d dx and hence, on multiplying by y =ẋ, that d dt One can now drop the restriction y = 0 by noting that (2.28) holds by straightforward differentiation and by applying (2.26). According to the Proposition 2.6, smooth trajectory exists, X(t) solves (2.26) subject to x(0) = 0 and e 2t y(t) → 0 as t → −∞. Let Z(t) = (2 − y 2 (t))Ω(x(t)), so that (2.28) implies The bounded convergence theorem easily implies that, for fixed t 1 , the righthand side of (2.31) converges to a finite limit as t 0 → 0 − , and hence so too must as t → 0, it follows that y is monotonic on (−δ, 0) for some small δ > 0 provided y is bounded away from zero on that set, and hence, in such circumstances, y 0 := lim t0→0− y(t 0 ) also exists. We note that y(t) cannot approach 0 as t → 0 − since a hyperbolic sink cannot be reached in finite time other than when x ≡ 0. Now set t 0 = 0 and t 1 = t in (2.31) and rearrange to obtain where Z(0) = (2 − y 2 0 ) has been used. Multiplying through by e 4t Ω(x(t)) and applying the condition e 2t y(t) → 0 as t → −∞ gives Putting this expression for y 0 into (2.32) gives it is clear that (2.34) may also be rewritten as (2.27).
Note that the integrand appearing in (2.34) is positive, which immediately implies that |y(t)| < √ 2 for all times t. The next result, which we shall need in due course, shows that this simple bound can be used to deduce the range of values of y(t) compatible withẏ(t) = 0.
The following series of results establish that the phase trajectory φ (−∞,0] , defined in (2.29), of the solution to (2.26) forms part of a heteroclinic orbit joining the rest points X −1 to X 0 or X 1 to X 0 . We begin by noting that the conservation law (2.27) implies that φ (−∞,0] belongs to a certain 'box' in the phase plane. Lemma 2.9. Suppose that X = (x, y) T obeys (2.26) subject to the conditions x(0) = 0 and e 2t y(t) → 0 as t → −∞, and let the trajectory φ (−∞,0] be given by (2.29). Then Proof. (i) The assumptions of the lemma imply that (2.27) and, equivalently, that (2.34) holds. We have already observed that the latter implies that |y(t)| < √ 2 for all finite t, so it must be that To control the horizontal extent of φ (−∞,0] , suppose for a contradiction that there is t * < 0 such that x(t * ) = ±π/2. Then Ω(x(t * )) ≤ Ω(x(s)) for all s, and it follows from (2.27) that y 2 (t * ) ≤ 0. This implies that φ (−∞,0] emanates from one of the rest points (±π/2, 0) at a finite time t * , which is impossible.
or it does not. If it does then y 2 (t) must be bounded away from the set L := (−π/2, π/2) × {0}. To see this, note that the phase portrait of (2.26) immediately implies that trajectories which approach L either cross it or approach the rest point X 0 . But since X 0 is a hyperbolic sink and φ (−∞,0] is parametrized by t ∈ (−∞, 0], X 0 cannot be an (ω-)limit point. Hence there is c > 0 such that y 2 ≥ c for all t < 0, from which it follows that |x(t)| must be unbounded. In particular, we can find t * < 0 such that x(t * ) = −π/2, contradicting the argument given in the first part of the proof of the lemma. If y crosses L only finitely many times then there is T < 0, say, such that φ (−∞,T ) lies either entirely in the lower half-plane or entirely in the upper half-plane. Either way, we can again argue that y 2 (t) must then be bounded away from 0 for all t < T for some T < T , and the same contradiction can be reached. This proves part (ii) of the lemma.
It is clear from (2.26) that F (x, 0) = −(p /Ω(x)) sin 2x(0, 1) T . Hence, as time t decreases the trajectory φ (−∞,0) must cross L − := (−π/2, 0) × {0} from the upper half-plane, to emerge in the lower half-plane, and that L + := (0, π/2) × {0} is crossed with the opposite orientation. We exploit this observation in the next result. Lemma 2.10. Let X and φ (−∞,0] be as in the statement of Lemma 2.9, suppose that X(0) = (0, y 0 ) T where y 0 > 0, and suppose that neither of X ±1 is an α−limit point of the trajectory. Then either: for all k and such that for all j ∈ N, Proof. Since, by hypothesis, neither of X ±1 is an α−limit point of the trajectory, we can apply Lemma 2.9(ii). Using this and the reasoning immedately preceding the statement of Lemma 2.10, the two possibilities for φ (−∞,0] must be as shown in Figs connecting the points P 0 and P 1 , as indicated. On the part of φ [T,0] lying in the upper half-plane, where y =ẋ, x(t) decreases monotonically as t decreases, so there is Arguing similarly for the part of φ [T,0] lying in the lower half-plane, there is Fig. 3. Since, by Lemma 2.9(b), φ (−∞,0] crosses L infinitely often, there must be a strictly decreasing sequence t k , say, of times such that y(t k ) = 0 for points X = (x, y) T belonging to φ (−∞,0] . Since y 0 > 0 by hypothesis, x(t 1 ) < 0, and, by the observation immediately preceding Lemma 2.10, φ (−∞,0] crosses from the upper half-plane to the lower half-plane and can only cross back into the upper half-plane at x(t 2 ) > 0. Since y| [t2,t1] < 0, x is monotonically decreasing for s in [t 2 , t 1 ]. Continuing in this way, it must be that there are times t 2j < t 2j−1 for all j ∈ N such that . Now consider, for j ≥ 2, the points (x(t 2j+1 ), 0) and (x(t 2j−1 ), 0) in the phase plane, , both of which lie in the upper half-plane, to cross, which is impossible for the smooth trajectory φ (−∞,0] . Hence x(t 2j−1 ) forms a strictly decreasing sequence contained in [−π/2, 0), which proves part (b)(iv) of the lemma. Part (b)(v) follows similarly.
The proof of (2.41) relies on (i) the more refined bound on y proved in Lemma 2.8 and (ii) on the geometry of φ (t2j ,t2j−1) , as depicted in Fig 3. Suppose that the point (x,ȳ) ∈ φ (t2j ,t2j−1) is such thatȳ ≤ y(t) for all t 2j ≤ t ≤ t 2j−1 . Since y(t 2j−1 ) = y(t 2j ), it must be that there is at least one t such that y(t) =ȳ andẏ(t) = 0. By construction, φ (t2j ,t2j−1) lies in the lower half-plane, which in particular means that for sufficiently large j there is c > 0 and times c, and control the length of each of these subintervals.
First consider the interval (t − c , t 2j−1 ). The second equation of (2.26) implies that When (2.46) is integrated over (t − c,j , t 2j−1 ) and the conditions y(t − c,j ) = −c and y(t 2j−1 ) = 0 are applied, we obtain Applying the inequality e 2r ≥ 1 + 2r with r = t 2j−1 − t − c,j and rearranging, we The corresponding argument for the interval (t 2j , t + c ) yields which is (2.41).

Lemma 2.12.
Suppose that X = (x, y) T obeys (2.26) subject to the conditions x(0) = 0 and e 2t y(t) → 0 as t → −∞, and let the trajectory φ (−∞,0] be given by (2.29). Then, provided X(t) is not identically zero, φ (−∞,0] must have α-limit point either X 1 or X −1 . Proof. Note that X solves (2.26) subject to the conditions stated above if and only if −X obeys the same system and conditions. Therefore, either both X −1 and X 1 are α-limit points or neither is. If neither of X ±1 is an α-limit point then, by exchanging X and −X if necessary and assuming that X(t) is not identically zero, we can suppose y(0) > 0. Hence either part (a) or part (b) of Lemma 2.10 must hold. If (a) holds then there are times then, by definition of Ω(x), this means that Ω(x(s)) > Ω(x(t 1 )) almost everywhere in the set (−∞, t 1 ). But, by (2.27), is then strictly negative, which is a contradiction. If |x(t 1 )| < x(t 2 ) then the same contradiction can be reached by replacing t 1 with t 2 in the above. Now assume that part (b) of Lemma 2.10 holds, and let t 2j−1 , t 2j be the times introduced there. Note that Lemma 2.11 also now applies. By (2.27) and (2.39), we have ).

(2.50)
NoDEA Area-preserving, buckling stationary points Page 23 of 50 6 The first term on the right-hand side of the line above converges to 0 as j → ∞. We will argue, using (2.40) and (2.41) that the second term is bounded away from 0. Indeed, from (2.40), it easily follows that , and so to bound e 4(s − j −t2j−1) away from zero it is enough to bound the right-hand side of (2.41) above uniformly in j.
We may regard c as a small but fixed quantity that is given via the geometry of φ (−∞,0] . Therefore the only way in which the right-hand side of (2.41) could become unbounded as j → ∞ is if (I) at least one of b + j − π/2 and b − j + π/2 tends to 0, (II) at least one of b ± j tends to 0, or both, were to hold for a subsequence. We can rule out (II) since it is clear from the construction that for a fixed and sufficiently small c, both b + j and b − j are bounded away from 0. To rule out (I), note that the behaviour of φ (−∞,0] in a neighbourhood of X −1 , say, is constrained by the behaviour of the linearized system associated with (2.26). Carrying out that analysis in a neighbourhood of X −1 reveals that the invariant stable and unstable manifolds at X −1 are tangent to the vectors (1, −1± This bounds b − j away from −π/2. A similar argument bounds b + j away from π/2, thereby ruling out (I). In conclusion, it follows from (2.41) that there is a constant M > 0 depending on l, l , c and such that By taking j large enough and bearing in mind the convergence x(t 2j−1 ) → l, we obtain the desired contradiction.
We now turn our attention to ω-limit points of extensions of φ (−∞,0] .
Let us focus on φ + R which, being constrained to lie in K, must approach one of X −1 , X 0 or X 1 as t → ∞. Suppose for a contradiction that ω-lim φ + R = X −1 . Then φ + R must enter the third quadrant R 2 −− := {(x, y) T : x < 0, y < 0} while, by symmetry, φ − R would enter R 2 ++ , with obvious notation. Bearing in mind that trajectories can only change quadrant in a clockwise manner as t increases, we find that φ + R and φ − R must cross, which is a contradiction. Referring to Fig. 4 may help here. Now suppose for a contradiction that ω-lim φ + R = X 1 . If φ + R crosses the line (0, π/2) × {0} then it must cycle at least once around the origin before an approach to X 1 is possible. But in so doing the trajectory would cross itself, which is impossible. Hence φ + R remains in the upper half-plane and x(t) increases monotonically from −π/2 to π/2. Let η > 0 and take t so large that cos 2 x(t) < η. Using the monotonicity of x(t), define t > 0 by cos 2 x(t ) = 2η, and note that we then have Using this inequality, we see that Since t > 0, the second term is negative and bounded away from zero, while the first term can be made arbitrarily small provided t is chosen larger still, if necessary. This contradiction shows that φ + R cannot approach X 1 , and so it must cross the line (0, π/2) × {0} and thereafter remain in the box K.
According to the Poincaré-Bendixson Theorem [9, Theorem 1.115], the ω−limit set of φ + R is either a fixed point, a regular periodic orbit, or it consists of finitely many fixed points joined by non-closed orbits. The function ϕ(x, y) := Ω(x) is a Dulac function for the system (2.26) in the sense that div (ϕ(x, y)F (x, y)) = −2 Ω(x) ≤ −2, which, in almost exactly the same way as Bendixson's criterion, [9, Proposition 1.133], rules out regular periodic orbits.
Thus X 0 must feature in the ω−limit set (because X 1 and X −1 cannot), and since it is a sink, as the linear analysis of (2.26) readily shows, points X(t) in the orbit φ + R must converge to X 0 as t → ∞. This concludes the proof for φ + R , and a similar argument yields ω-lim φ − R = X 0 .

Properties of the functions v ±
We recall that the twist maps v ± are given by v ± (x) = Re R (μ ± k(R)) The dependence of k(R) on is through this ODE, and, strictly speaking, each k should be written k (R), say. However, for brevity we write k(R) when we mean k (R) solving (2.54). It was established in Proposition 2.6 that k(R) is smooth on the open interval (0, 1), so it follows that v ± are smooth on B \ {0}. We now use properties of the solutions of system (2.26) established in the previous subsection to study the behaviour of ∇v ± , and to prove, in particular, that these maps are continuous on the ball B. Proposition 2.14. Let v ± be given by (2.1) where k lies in C and satisfies (2.22) subject to (2.23). Then ∇v ± is continuous on B.

It follows that
and hence that the Hölder exponent α, say, of ∇v + must obey (2.55) It follows from this and the fact that stationary points are smooth away from 0 that v + ∈ C 1,α( ) (B), where α( ) is given in (2.55) for the stated range of .  (2.22) subject to the conditions k(1) = 0 and lim R→0 R 3 k (R) = 0.
Proof. (i) Our contention here is that, when ∈ (1, √ 2), the heteroclinic φ + R connecting X −1 and X 0 is entirely contained in the region Δ, which we define to be the interior of the triangle formed by the points (−π/2, π/2), X 0 and X −1 .
In particular, it cannot cross the y−axis. Consider the line S, say, joining the points (−π/2, π/2) and X 0 , and let n = (1/ √ 2, 1/ √ 2) be its upward-pointing normal. Then, with F (x, y) given by (2.25), we find that on S 2 −1 and note that g (x) > 0 if and only if For 1 < ≤ √ 2 and −π/2 < x < 0, the infimum of the left-hand side of (2.57) agrees with the supremum of the right-hand side, and in fact for all values apart from = √ 2 and x = 0, (2.57) holds with strict inequality. Hence it must be that g (x) < g (0) = 0, from which it follows from (2.56) that n · F < 0 at all points along S. Thus trajectories in Δ cannot leave by crossing at any point on S. Moreover, since it is also clear that trajectories cannot leave Δ by crossing the lines joining X −1 to (−π/2, π/2) or X −1 to X 0 , and so they must remain in Δ. By Lemma 2.12, if X solving (2.26) is not identically zero, then its trajectory φ + R has α−limit point X −1 or X 1 , and, by symmetry, we may assume without loss of generality that the former is true. Thus φ + R meets Δ non-trivially and hence, by the above, we must have φ + R ⊂ Δ. In particular, φ + R cannot cross the yaxis, contradicting the assumption that X solving (2.26) and its associated boundary conditions is non-zero. It follows that X ≡ 0 is the only possible solution, which concludes the proof of part (i) of the lemma. (ii) By Lemma 2.13, if X solving (2.26) and its associated boundary conditions is non-zero then X(t) → X 0 as t → ∞. Therefore, for sufficiently large T , φ ± (T,∞) lies in a neighbourhood of X 0 in which the behaviour of solutions to (2.26) is governed by the linear systemẊ = ∇F (X 0 )X. When > √ 2, the eigenvalues of DF (0) are a complex conjugate pair To see that there are infinitely many solutions of (2.26) subject to x(0) = 0 and e 2tẋ (t) → 0 as t → −∞ when > √ 2, consider the following. Let X(t) lie on φ + R , the heteroclinic connecting X −1 and X 0 , and let t j > 0 be one of the times such that x(t j ) = 0, noting that t j exists precisely because φ + R describes a spiral when > When > √ 2, each additional solution X j (t) =: (x j (t), y j (t)) described above yields a new stationary point of the form v +,j (x) = Re R (μ + k j (R)), say, in which k j (R) = x j (t) and R = e t . Properties of the trajectory φ + (−∞,0] imply that v +,j introduces 'additional buckling' near the centre of the ball B, and this is borne out by numerical simulations, an example of which is given in Fig. 5. The first buckle occurs when R ∼ 0.2 and the second is just visible when R ∼ 0.01. Another family of maps v −,j can be obtained similarly from φ − R , and by symmetry we have The index j is a useful proxy for the number of buckles a solution produces near the centre of B. In Sect. 4 the energies E j ( ) := E (v −,j ) are calculated and are found to form a sequence that increases with j, thereby confirming our natural intuition that the more a solution buckles, the greater the energetic cost.

Stationary points of D subject to the area-preserving constraint: Lagrange multipliers with logarithmic growth
The purpose of this section is twofold. Firstly, we make clear that the maps v ± are not just stationary points in the class C, but that they are also stationary in the larger class A in the sense that By invoking results of [8,13], a natural corollary of this condition is that a Lagrange multiplier, or pressure, λ ± , say, can be rigorously associated to each function v ± through the equation See Lemmata 3.2, 3.3 and Proposition 3.4 for the details. Next, since each v ± is C 1 and invertible, the change of variables p ± (x) := λ ± (y) Thus the second main purpose of this section is to characterise, as far as possible, the functions p ± appearing in (3.3). To do this we exploit the facts that: (i) each p ± is unique up to an additive constant, and (ii) p ± belongs to ∩ s>1 L s (B), both of which follow from the results of [8,13], as we explain below. The rest of the argument goes as follows, and with U := v ± for brevity: Step 1 find, for any δ ∈ (0, 1) a function P , say, depending only on R = |x|, such that Step 2 verify, by an approximation argument that, as a result of (3.4), which is just (3.3) with P in place of p ± , and Step 3 invoke the uniqueness articulated in (ii) above to conclude that there is a constant c, say, depending on v ± such that p ± = P + c a.e. in B. (3.5) The function P determined in Lemma 3.7 is given explicitly in terms k solving (2.54), and the characterisation so obtained shows, in conjunction with (3.5), that: (a) p ± is a smooth, radial function on B \ {0}, and, to within a constant, (b) p ± ∼ 2p ln |x| as |x| → 0, where p = − 1 .

Equation (3.4) is a form of the equation
Following, for example, [28, Section 2], in which a useful, formal derivation is given of the strong form of (3.6) in the case of the Dirichlet energy, we note that away from the origin, i.e. where U is smooth, (3.6) together with Piola's identity div cof ∇U = 0 gives div (D F W (x, ∇U ; )) = cof ∇U (x) ∇P (x) in B \ B(0, δ).

Stationarity in the sense of (3.1)
We begin with the following straightforward algebraic lemma.
Lemma 3.1. Let W (x, F ; ) be given by (1.4), that is Then Proof. For any fixed a in R 2 , it is the case that as is readily checked. Taking first a =x and then a = Jx, (3.8) follows from this, the definition of W (x, F ; ) and the fact that |adj Fx| 2 = |F Jx| 2 .
Throughout the rest of this section, it will be convenient to make the following standing assumptions.
and V (y) = Re R (μ − k(R)) for y = Re R (μ) (3.10) denote the functions v + and v − respectively; (SA) 2 The function k(R) in (3.9) and (3.10) solves (2.22) for R ∈ (0, 1), subject to k(1) = 0 and lim R→0+ R 3 k (R) = 0; (SA) 3 We suppress the dependence of U , V and k on , but bear in mind that there is always an -dependence through equation (2.22). Having introduced U and V as above, we note that and a direct calculation shows that and (3.14) The aim is now to demonstrate that (3.1) holds, and we do this in two steps. The first, in Lemma 3.2 below, reformulates (3.1) as the equivalent problem of showing that where V is given by (3.14) and Σ(R, V ) is a vector quantity defined in (3.19) below. Lemma 3.3 confirms (3.15), and thereby completes the proof of (3.1).

23)
NoDEA Area-preserving, buckling stationary points Page 33 of 50 6 and Now, a short calculation shows that for any φ in C ∞ c (B, R 2 ) and any function g in C 1 (0, 1), it holds that which is true provided [R 2 k Q 2 φ] 1 0+ = 0. The upper limit vanishes because φ has compact support in B, while the lower limit is zero by appealing to Proposition 2.14. The latter shows that ∇U is continuous on the ball B, which in particular implies that Rk is bounded as R → 0+. Hence it is certainly the case that R 2 k → 0, so that (3.27) holds.
Hence, by integrating (3.26), we must have The integrand of the right-hand side of the last equation agrees with where Q 3 is given by (3.18). But a short calculation confirms that Hence (3.28) holds, and the proof is complete. 3 be in force, and suppose that Σ(R, V ) is given by (3.19). Then condition (3.20) is satisfied, and hence U is a stationary point of D in the sense of (3.1).
Proof. Since any function φ ∈ C ∞ c (B, R 2 ) such that div φ = 0 can be written (3.29) so it remains to calculate e τ · Σ and e R · Σ. To do this, the following facts, which are easily derived from (3.16), (3.17) and (3.18) and (3.12), are needed: It follows that By (2.54), the collection of terms in parentheses vanishes on (0, 1), and so we conclude that It is not necessary to calculate Σ(R, V ) · e R in detail. Instead, it is enough to observe that the coefficients of e R and e τ in Q 1 , Q 2 and Q 3 depend only on R. Hence, which vanishes trivially. Using this and (3.30) in (3.29) concludes the proof of the lemma. We are now in a position to apply the results of [13] or [8] to deduce the existence of a Lagrange multiplier satisfying (3.2).

Proposition 3.4. Let the conditions of Lemma 3.3 be in force. Then there exists
Proof. Proceeding as in the first three lines of the proof of Lemma 3.2, and in the notation introduced there, (3.21) and (3.1) imply that the linear functional vanishes on divergence-free φ. Let α > 1. Since V belongs to C 1 (B), it follows that F can be extended by a density argument to a linear functional F : W 1,α 0 → R which, again, vanishes on divergence-free functions, suitably interpreted. In the notation of [13,Equation 2.1], F is now defined on the class U 1,α (B), and so [13,Theorem 3.1] applies. Thus there is λ ∈ L α (B) such that (3.33) Since α > 1 was arbitrary, it follows that λ ∈ ∩ s>1 L s (B), as claimed.

Properties of the pressure
We now follow Steps 1-3 of our argument as set out in the introduction to Sect. 3. Specifically, we begin by seeking a function P (x) obeying (3.34) below, i.e. for each δ > 0, and with (SA) 1 -(SA) 3 in force, we wish to solve (3.34) To do so we require a short technical lemma whose aim is to calculate efficiently div D F W (x, ∇U ; ) when U is given by (3.11).
Lemma 3.6. Let U = Re R (μ + k(R)) represent one of the functions v ± , and suppose δ > 0. Then (3.36) Proof. Using (3.8) together with (3.11) gives Hence, for any A short calculation shows that the right-hand side of the last equation is − B H · η dx, where H is given by (3.36). Since η is arbitrary, (3.35) follows.
The next result identifies the function P satisfying (3.34) by calculating ∇U T div D F W (x, ∇U ; ), which turns out to be a function of the form g (R)e R for a suitably chosen g, and then setting P (x) = g(R). The details are as follows.
Then P is such that NoDEA Area-preserving, buckling stationary points Page 37 of 50 6 Proof. By (3.13) and Lemma 3.6, we see that and, by letting H = H R e R + H τ e τ and using (3.11), this gives (3.38) Using (3.36) and then (2.54), we see that To deal with the coefficient of e R in (3.38), for brevity we first label each occurrence of R 2 k + 3Rk in H R and H τ with γ, say, and calculate The second and last terms cancel, since (2.54) can easily be rearranged to read Using (3.39), the remaining terms can be written as Using the easily-checked identity cos 2 k + sin 2 k cos 2 k + sin 2 k − p 2 sin 2 k cos 2 k = 1, the term with prefactor R 2 k 2 is while the other term can be handled by using the identity cos 2 k + sin 2 k cos(2k) + 2p sin 2 k cos 2 k = cos 2 k − sin 2 k . Thus By choosing P as in (3.37), we see that (3.34) holds.
This completes Step 1 of the argument outlined at the start of Sect. 3.
Step 2 consists in showing that with P and U as in Lemma 3.7, the equation holds. This can be done using a short cut-off function argument, which we now outline.
Lemma 3.8. Let (SA) 1 -(SA) 3 be in force, and let P be given by (3.37). Then (3.40) Proof. It has already been established in the introduction to Sect. 3 that (3.34), as derived above, is equivalent to the pointwise condition so that it is certainly the case that for any ϕ ∈ C ∞ c (B\B(0, δ), R 2 ) and any δ ∈ (0, 1), say. Let η δ be a smooth cutoff function which satisfies η δ (x) = 0 if x ∈ B(0, δ), η δ = 1 if x ∈ B \ B(0, 2δ), |∇η δ | ≤ c δ , and 0 ≤ η δ ≤ 1 on B. Let φ ∈ C ∞ c (B, R 2 ) be arbitrary, take ϕ := η δ φ in (3.41) and, for brevity, let Y (x) = D F W (x, ∇U ; ) − P (x) cof ∇U . The result is that the first term of which converges to B Y · ∇φ dx as δ → 0, while the second obeys, for any s > 1 and some constant c depending only on η, the estimate by Hölder's inequality. By Lemma 3.9 below, we see that Y belongs to ∩ s>1 L s (B), and so by taking This concludes the proof of the lemma.
The following result is needed to complete Step 2.
Lemma 3.9. Let (SA) 1 -(SA) 3 be in force, and let P be given by (3.37 belongs to L s (B) for all s > 1.
Proof. By (3.37), we have In order to control the term R 2 k 2 , we exploit the conservation law (2.34) derived in Sect. 2.3, which is expressed in terms of the variables x(t) = k(e t ), R = e t and y(t) =ẋ(t) by To convert x · ∇P (x) to these variables, note that y(t) = Rk (R), and let and .
Then it follows that 8G(s, t)e 4(s−t) ds, (3.46) where the fact that t −∞ 4e 4(s−t) ds = 1 has been used. Now, by applying the elementary estimates −1 ≤ Ω 1 ≤ and −2 ≤ Ω 2 ≤ 2 , we see that Substituting into (3.46), we obtain Hence, by writing where A := + 3 3 , part (b) follows too. These bounds can be refined for small R. Specifically, by applying Lemma 2.13, we deduce that cos 2 x(t) → 0 and sin 2 x(t) → 1 as t → −∞, and hence that the limit being taken with respect to s obeying s < t since this is the relevant range for the integration performed in (3.46). Thus x · ∇P (x) ∼ 2p as R → 0.
Arguing as above, we therefore have The final assertion in the statement of the lemma follows by putting together the fact that U belongs to C 1 (B) and estimate (3.43).

NoDEA
Area-preserving, buckling stationary points Page 41 of 50 6 An approximation argument almost identical to that which leads to (3.42) in Lemma 3.8 then implies that It follows that the weak derivative ∇λ + exists and agrees a.e. with ∇(P • V ), from which the conclusion of the lemma easily follows.
Remark 3.11. By inspecting the expression (3.37) defining P , it is clear that and in particular that ∇P is the same in each case. Therefore, by Proposition 3.10, what distinguishes the Lagrange multipliers λ + and λ − is the second map in the composition, namely modulo constants.

Numerics
Although useful in many ways, the analysis of the stationary points given in the preceding sections tells us little about the ordering, if any, of the energies E(v ±,j ) for j ∈ N. It seems natural, therefore, to study this problem numerically, and this is the point of the current section. We begin by considering numerical solutions to the ODE (2.24) in the form of the system (2.26), which we recall is given bẏ subject to the initial conditions x(0) = 0, e 2t y(t) → 0 as t → −∞. We refer to (4.1) with these initial conditions as the 2-d system. When we also need to compute the energy obtained by substituting x(t) = k(R) and R = e t into (2.6), and where the integral takes place along the heteroclinic orbit φ + R , we solve numerically the 2-d system augmented by the additional ODĖ E = e 2t H(x(t), y(t)), (4.3) referring to (4.1) together with (4.3) as the 3-d system. lutions v ± as varies and (b) to compute the sequence given by where the usual initial conditions apply and where > √ 2. When this is the case, the times t * j , j ∈ N, form a strictly increasing sequence for which x(t * j ) = 0 and y(t * j ) has alternating sign, with t * 1 = 0 and y(0) > 0. For reasons indicated in the discussion at the end of Sect. 2.3, we expect that the E j ( ) will form an increasing sequence.
In the rest of this section, we discuss the computation of: (1) E(0, ), which is obtained by solving the 3-d system, with the given value of and the usual initial conditions, from t = −∞ to t = 0; (2) y 0 ( ), ∈ ( √ 2, ∞), the point where the heteroclinic orbit φ − R first crosses the y-axis, this point being (0, y 0 ( )); and (3) E j ( ) as defined in (4.5).
We remark that the computation of E j ( ) is most easily done by first splitting the range of integration (−∞, t * j ] in (4.4) into (−∞, 0] and [0, t * j ] and then using the results of computation (1) to calculate the contribution from (−∞, 0). The integral over [0, t * j ] is found by solving the 3-d system with initial conditions x(0) = 0 and y(0) = y 0 ( ). Hence, it is necessary to calculate y 0 ( ) first in computation (2).
We remark that, from a numerical point of view, both the 2-d and 3-d systems are smooth and non-stiff, and finding accurate numerical solutions is comparatively straightforward. We give some details of how to compute, among other things, the energy at specified times, given a standard ODE solver implemented in arbitrary precision. 4 Such a solver needs, first of all, to have controllable accuracy, and second, to have the capability to solve an ODE until a given function of the solution reaches a specified value (so-called 'event handling'). Note also that the solver must allow the initial conditions to apply at any given finite time, with integration proceeding from that point in either direction in time. Several suitable ODE solvers are provided, for instance, in the algebraic manipulation program Maple, in which the accuracy can be controlled both by setting the number of significant figures for computation, and by specifying the desired absolute and relative error tolerances in the solution. where H is defined in (4.2). The integral is along the heteroclinic orbit φ − R between the points X −1 and the first crossing of the y-axis, which takes place at (x, y, t) = (0, y 0 ( ), 0). In practice, we estimate E(0, ) numerically by approximating the integral in (4.6) as where H num is the same function as H(x(t), y(t)) but with the dependence on the variables of interest made explicit. In particular, we show the initial conditions used when solving the 2-d system, these being specified at t = T i , which is also the lower limit of the integral; T i will end up being negative but finite.
Refer now to Fig. 6. Fix δ small and positive and let C be a circle of radius proportional to δ, say bδ, with centre X −1 . The eigenvector v + is shown intersecting C at point P = (−π/2 + δ, δ(γ − 1)). Hence, b 2 = (γ − 1) 2 + 1. Choosing a small value for δ, we then solve the 2-d system using point P as the initial condition at t = 0, and stopping when the condition x(T ) = ξ(T ) − π/2 = 0 is met -this is where the event handling capability of the ODE solver comes into play. This then gives an estimate of T (δ), with T (δ) increasing as δ decreases (at least, for δ not too small). Since the 2-d system is autonomous, we can then set T i = −T in (4.7) and solve the 3-d system, also starting from P , which we take to be the initial condition but now at t = −T . We set the third initial condition E(−T ) = 0, where E(t) solves (4.3), and our estimate of E(0, ) is then just E(0). In practice, we set δ to values in the range 10 −5 -10 −14 , depending on , computing T in each case. We then compare values of E(0) + c(T ) for different values of δ. Their closeness gives us confidence in our estimate of E(0, ). The results are shown in Fig. 7, from which it is clear that, in particular, E(0, ) increases with , at least for the range of considered.
With the initial conditions x(0) = 0,ẋ = y 0 , this has the closed-form solution x(t) = arcsin y 0 e −t sin t , (4.10) which we shall use presently.
For finite values of > √ 2, however, it appears that only numerical solutions are available. In fact, we compute y 0 ( ), where (0, y 0 ( )) is the point at which the heteroclinic orbit φ − R crosses the y-axis at t = 0. We carry out this computation in two ways. In the first computation, we integrate the 2-d system forwards in time, for given , starting from point P = (−π/2 + δ, δ(γ − 1)), as in the previous section, until the solution first crosses the positive y-axis. The point at which it does this is an estimate of y 0 ( ), the accuracy of the estimate depending critically on the closeness of the initial point to (−π/2, 0) -that is, on how small δ is. We cannot decrease δ indefinitely, however, as solving the ODE numerically becomes difficult if the initial condition is too close to X −1 .
For the second computation, we fix and choose two positive values of y(0), y lo 0 < y hi 0 , where, as t → −∞, the solution X(t), starting from X(0) = (0, y lo 0 ) crosses the x-axis in the interval (−π/2, 0), whereas the solution starting from X(0) = (0, y hi 0 ) does not. The latter case can easily be detected by checking, as the solution proceeds, whether, while y(t) > 0, x(t) < −π/2: since y =ẋ > 0, as t → −∞, x can only decrease, so a solution that crosses the line x = −π/2 with y > 0 cannot then cross the x-axis in the interval (−π/2, 0). A bisection algorithm can then be used to refine the original bounds y lo 0 , y hi 0 until the desired accuracy of the estimate of y 0 ( ) is achieved.
The results of the two computations agree to many significant figures, and are shown in Fig. 8, which gives a plot of y 0 ( ) versus . The first algorithm is significantly faster. Figure 8 implies that (i) y 0 ( ) monotonically increases for ∈ ( √ 2, ∞). (ii) lim → √ 2+ y 0 ( ) = 0.