Energy minimising configurations of pre-strained multilayers

We investigate energetically optimal configurations of thin structures with a pre-strain. Depending on the strength of the pre-strain we consider a whole hierarchy of effective plate theories with a spontaneous curvature term, ranging from linearised Kirchhoff to von K\'arm\'an to linearised von K\'arm\'an theories. While explicit formulae are available in the linearised regimes, the von K\'arm\'an theory turns out to be critical and a phase transition from cylindrical (as in linearised Kirchhoff) to spherical (as in von linearised K\'arm\'an) configurations is observed there. We analyse this behavior with the help of a whole family $(\mathcal{I}^{\theta}_{\rm vK})_{\theta \in (0,\infty)}$ of effective von K\'arm\'an functionals which interpolates between the two linearised regimes. We rigorously show convergence to the respective explicit minimisers in the asymptotic regimes $\theta \to 0$ and $\theta \to \infty$. Numerical experiments are performed for general $\theta \in (0,\infty)$ which indicate a stark transition at a critical value of $\theta$.


Introduction
The topic of this paper is motivated by experimental observations on optimal energy configurations in thin (heterogeneous) structures with a pre-strain. The simplest example of such a structure is the classical bimetallic strip which consists of two strips of different materials with different thermal expansion coefficients joined together throughout their length. If heated or cooled, due to the misfit of equilibria, internal stresses develop. The flat reference configuration is no longer optimal and the strip bends in order to reduce elastic energy. This behavior can effectively be modelled with a 1d energy functional comprising a temperature dependent spontaneous curvature term.
In this paper we will investigate thin layers whose two lateral dimensions are much larger than their very small height and whose flat reference configuration is subject to internal stresses (one speaks of pre-strained or pre-stressed bodies). Examples of such structures are heated materials (with inhomogeneous expansion coefficients as in the bimetallic strip referred to above or homogeneous materials with a temperature gradient), crystallisations on top of a substrate as in epitaxially grown layers, or biological materials whose internal misfit is caused by swelling and growing tissue. Our main focus will be on multilayered heterogeneous plates, for which the effective plate theories have been provided in [11]. Our findings, however, apply equally to different situations as long as they are described by the same effective functionals, cf. Remark 1 below.
As a matter of fact, the situation is much more complicated and interesting for two dimensional plates than for one dimensional strips. It has been found that the assumed shape depends on the strength of the pre-strain and the aspect ratio of the specimen: Large pre-strains in very thin layers tend to cause cylindrical shapes whereas smaller pre-strains in thicker layers lead to spherical caps, [25,31,13,14,21,12]. To explain this observation one argues that locally the energy is best released if a spherical shape is assumed. If, however, the aspect ratio is very small, i.e., the lateral dimensions are very large compared to the thickness, then this leads to geometric incompatibilities: non-zero Gauß curvature introduces a change of the metric which by far has too high elastic energy. In contrast, cylindrical shapes do not lead to such incompatibilities.
A thorough theoretical understanding of this mechanism through which 'misfit' of equilibria is converted into mechanical displacement is not only interesting from a mathematical point of view. In view of applications it has proved to constitute a convenient and feasible method to access and manipulate objects even at the nanoscale. By way of example we mention experiments on the self-organised fabrication of nano-scrolls, as reported in [34,18,27].
The aim of this paper is to shed light on the geometry of energetically optimal configurations of pre-strained heterostructures with the help of two-dimensional plate theories. More precisely, we consider effective plate theories for multilayers with reference configuration Ω h = ω × (−h/2, h/2), 0 < h 1, whose (small) misfit pre-strain is described by a matrix h α−1 B h , scaling with h.
The particular case α = 2 with a misfit of the order h of the aspect ratio has been investigated in [32,33,7]. The appropriate plate theory is the nonlinear Kirchhoff theory (in the finite bending regime) and energy minimizers turned out to be (portions of) cylinders whose possible winding directions and radii are determined explicitly. Therefore, in order to be able to encounter different behavior one has to consider weaker scalings of the misfit.
In [11] -based on the homogeneous case explored in [15] -we have found a whole hierarchy of effective plate theories for the scalings α > 2. Suitably rescaled, one obtains only three different limiting plate theories: the linearised Kirchhoff theory for α ∈ (2, 3), the von Kármán theory for α = 3 and the linearised von Kármán theory for α > 3. With a view to our present investigation, we have moreover derived a fine scale θ in the critical von Kármán scale which interpolates continuously between the the two linearised theories.
For such small misfits one is lead to describe a deformation y h : Ω h → R 3 in terms of the scaled and averaged in-plane, respectively, out-of-plane displacements where θ ≡ 1 unless α = 3 and A limiting plate theory in terms of the limiting quantities (u, v) is then derived as the Γ-limit of the 3d nonlinearly elastic energy, rescaled by h 1−2α , cf. [11]. For a minimizer (u, v) of the limiting theory one obtains the shape of an optimal configuration at finite 0 < h 1: After descaling, its x 3 -averaged displacement is given approximately by Since γ > α − 2, the in-plane components are indeed much smaller than the out-of-plane component. In his sense, the shape is to leading order described by v : ω → R only. In the linearised regimes our results give the following picture: If α < 3, degenerate parabolas (infinitesimal parts of cylinders) are seen to be optimal, whereas for α > 3, non-degenerate parabolas (infinitesimal parts of an elliptical cap) are energy minimizers. Only in the latter case, however, the minimizer is unique (up to affine terms). Yet, even in case α < 3 it turns out the geometric shape is uniquely determined as an infinitesimal part of a cylinder while the winding direction and radius may have several optimal values. In both cases we explicitly determine these minimizers. A basic observation shows that for α = 3 these configurations are still asymptotically optimal in the 'almost linearised' regimes θ 1 and θ 1, respectively. The von Kármán regime is much more subtle. We focus on a prototypical functional in order to understand better the material response if the misfit (and hence θ) is increased from 0 to a finite value. We show that for finite, although small, values of θ there is a unique branch of global minimizers emanating from a spherical cap. For a further study for general values of θ ∈ (0, ∞) we then rely on computer experiments. To this end, we develop a penalised, nonconforming finite element discretisation using P 1 elements and employ projected gradient descent to solve the ensuing nonlinear problems while ensuring constraints are met. We first show Γ-convergence of the discrete problems to the continuous one, then investigate the minimizers in their dependence on θ. Interestingly, our results seem to indicate a stark change of material response at a critical value of θ, showing a symmetry breaking 'phase transition' from a nearly spherical cap to an approximate cylinder.

Outline
We begin by recalling our main results from [11] in order to provide the appropriate plate theories in Section 2. There we also identify the effective elastic moduli and spontaneous curvature terms explicitly so as to transform the problem into a more amenable form to identify minimizers. We then discuss the linearised regimes α ∈ (2, 3) and α > 3 as well as the asymptotic von Kármán regimes θ → 0 and θ → ∞ in Section 3. The structure of minimisers for small θ is investigated in Section 4. Finally, Section 5 contains our numerical findings.

Effective plate theories
We first recall the main results of our contribution [11] on a hierarchy of plate theories for pre-strained multilayers derived from non-linear three dimensional elasticity by Γ-convergence. We then determine the effective (homogenised) elastic moduli and corresponding quadratic energy desnities of the plates in terms of the moments of the pointwise elastic constants of the layers.

Dimension reduction for pre-strained multilayers
Working exactly in the setting of [11] we consider a thin domain where ω ⊂ R 2 is bounded with Lipschitz boundary, 0 < h 1, subject to a deformation w : Ω h → R 3 . Changing variables form x 3 to x 3 /h we obtain a deformation mapping y(x) = w(x 1 , x 2 , hx 3 ) and the energy per unit volume where the elastic energy density W h α depends on a scaling parameter α ∈ (2, ∞) and is given by describing the internal misfit and W 0 the stored energy density of the reference configuration. For α = 3 we include an additional parameter θ > 0 controlling further the amount of misfit in the model: We take W 0 fulfilling the usual assumptions of smoothness around SO(3), frame invariance, boundedness and quadratic growth which are detailed in [11]. After linearising around the identity, one obtains the Hessian for t ∈ (−1/2, 1/2) , F ∈ R 3×3 and defines Q 2 by minimising away the effect of transversal strain on Q 3 : for t ∈ (−1/2, 1/2) , G ∈ R 2×2 , e 3 = (0, 0, 1) ∈ R 3 , andĜ ∈ R 3×3 has G as its upper left 2 × 2 submatrix and zeros in the third column and the third row. The functions Q 2 (t, ·), t ∈ (−1/2, 1/2), are quadratic forms on R 2×2 which are positive definite on R 2×2 sym and vanish on antisymmetric matrices. Moreover, they satisfy the bounds for constants c, C > 0 and a.e. t ∈ (−1/2, 1/2). We also denote byB(t) the 2 × 2 matrix which arises from B(t) ∈ R 3×3 by deleting its last row and last column. ThenB From Q 2 (t, ·) andB(t) we define the effective form: and its relaxation In [11] it is shown that h 2−2α E h α Γ-converges for the convergence of the averaged in-plane and out-of-plane displacements (u h , v h ) (u, v) in W 1,2 (ω; R 3 ) modulo a global rigid motion, cf. (1), to the following effective limiting functionals: For the scaling α ∈ (2, 3) as defined in [11] and convex ω, the linearised Kirchhoff energy is given by For α = 3 we have the von Kármán type energy 3 Finally, in the regime α > 3 we have the linearised von Kármán energy 3 As in [11] we slightly overload the notation in what would be a double definition of I h 3 , using the letter in the subindex to dispel the ambiguity.
Remark 1 The precise assumptions on W h α from [11] are not essential for the results of the present contribution. In what follows we will only need that the Q 2 (t, ·), t ∈ (−1/2, 1/2), are quadratic forms on R 2×2 that vanish on antisymmetric matrices and satisfy (2) and thatB satisfies (3).
The existence of minimizers of (5) (6) and (7) follows by a standard application of the direct method or, in the setting of [11], as a direct consequence of Γ-convergence and compactness.
Example. For a homogeneous material Q 2 (t, A) = Q 2 (A) with linear internal misfit B(t) = tI one has for v ∈ W 2,2 sh (ω), respectively, (u, v) ∈ W 1,2 (ω; R 2 ) × W 2,2 (ω; R). These functionals, where the elastic coefficients do not depend on the out-of-plane component, can model for instance a single-layer material under thermal stress. In Section 5, we will study the energy (8) as a function of θ.

Effective moduli and minimising strains
This subsection serves to give explicit formulae relating the homogenised effective elastic moduli found above to the zeroth, first and second moment in t of the individual Q 2 (t, ·). We also identify their pointwise minimiser so as to rewrite the effective quadratic forms in their most convenient form. The computations are completely elementary, we indicate the main steps.
Because Q 2 vanishes on antisymmetric matrices we may restrict our attention to F ∈ R 2×2 sym . From now on, we identify matrices and analogously F → f ,B → b, A → a. Then, for each t ∈ (−1/2, 1/2) there exists some symmetric, positive definite matrix M (t) such that for all A ∈ R 2×2 sym : We define the moments of M as It is easy to see that (2) implies that M 0 and M 2 are positive definite. We claim that also M * := M 2 − M 1 M −1 0 M 1 is positve definite. To see this, fix Λ ∈ R 2×2 and note that for all x ∈ R 2 \ {0} e. t would imply that (tI − Λ)x = 0 in contradiction to Λ having at most two eigenvalues. Expanding the square we get 0 < Let Q 2 be given as in (2.1). Elementary calculations show that We define the linear mappings L i , L * : R 2×2 sym → R 2×2 sym , i = 1, 2, 3, by and the positive definite quadratic forms Q 0 2 and Q * 2 on R 2×2 sym by In terms of these quantities our computation reads with Minimizing out E yields 3 Optimal configurations in the linearised and the asymptotic critical regimes In this section we develop a characterisation of minimisers for the lower range α ∈ (2, 3) and for the upper range α > 3 of scalings. Recall from the discussion in Section 1 that we are primarily intested in the shape of the out-of-plane component v. The results indicate that the characteristic shapes in the limit h → 0 are (infinitesimal) cylinders and paraboloids respectively. Invoking the Γ-convergence results with respect to the interpolation parameter θ from [11, Section 6] this will also shed light on the optimal shapes in the asymptotic regimes θ → 0 and θ → ∞ for the von Kármán scaling α = 3. We collect our results in the following three theorems, where indeed Theorem 1 is indeed rather an elementary observation based on our preparations form the previous section and Theorem 3 is a direct consequence of [11,Section 6]. We allow for a general bounded Lipschitz domain ω in these theorems.

Theorem 1
The minimisers of I lvK , eq. (7), are of the form with E 0 , F 0 ∈ R 2×2 sym the constants from (13). u is unique up to an infinitesimal rigid motion and v up to the addition of an affine transformation.
Theorem 2 Up to the addition of an affine transformation, the minimisers of I lKi , eq. (5), are of the form where Q * 2 , F 0 are given in (11) and (13), respectively.
Remark 2 Describing symmetric 2 × 2 matrices A by vectors a ∈ R 3 as in Section 2.2, the set N is the set of touching points of the two quadrics {a ∈ R 3 : a 1 a 2 − a 2 3 = 0} (a cone) and {a ∈ R 3 : a M * a = c m } (an ellipsoid), where intersecting with an affine plane P containing three distinct points of N shows that N ∩ P is an ellipse and then even N ⊂ P . This shows that either #N = 1 and there is a unique minimizer, or #N = 2 and there are precisely two minimizers, or N is an affine ellipse and to each 'winding direction' Re, e ∈ S 1 , there is a unique curvature λ = λ(e) such that Theorem 3 Suppose that (u θ , v θ ) are minimisers of I θ vK , eq. (6). a) As θ → 0, up to infinitesimal rigid motions in the in-plane component and up to the addition of affine transformations in the out-of-plane compenent, b) As θ → ∞, up to the addition of affine transformations in the out-of-plane component and up to passing to a subsequence, v θ v in W 2,2 (ω; R) with v as in (16).
With the help of the Vitali covering theorem we can exhaust ω up to a set of negligible measure with disjoint convex subdomains ω 1 , ω 2 , . . .. Denoting the accordingly restricted functionals by I θ vK ( · ; ω n ), I lKi ( · ; ω n ) we have where we have made use of the lower bound in the Γ-convegence of I θ vK (·; ω n ) to I lKi (·; ω n ), see [11,Theorem 8], in the third step and of Theorem 2 in the fourth step. So we must have I lKi (v; ω n ) = I lKi (v; ω n ) for all n and hence ∇ 2 v ∈ N a.e. on ω and so the claim follows from Theorem 2. 2 As for Theorem 2, it is straightforward to see that v as defined in the theorem is a minimisers of I lKi . However, the proof that every minimiser of I lKi is necessarily of this form needs some work. The difficulty lies in excluding the possibility of constructing a minimiser by piecing together functions whose Hessian belongs to the set N , all with minimal energy but lacking a nice global structure. Yet it is possible to obtain a global representation of the Hessian which shows that it must be constant over ω so that minimisers are (up to an affine transformation) indeed cylindrical. In order to do this we require (cf. [28]): Definition 1 Let ω ⊂ R 2 a convex bounded domain and y ∈ W 1,2 (ω , R 3 ) be an isometry. A connected maximal subdomain of ω where ∇y is constant and y is affine whose boundary contains more than two segments inside ω is called a body. A leading curve is a curve orthogonal to the preimages of ∇y on the open regions where ∇y is not constant, parametrised by arc-length. We define an arm to be a maximal subdomain ω(γ) which is covered (parametrised) by some leading curve γ as follows: where ν(t) = γ (t) ⊥ . We also speak of a covered domain.
γ ω ′ Figure 1: The partition of ω into bodies and arms. ∇y is constant in the bodies (colored) and along each of the straight lines making up the arms (white).
The existence of covered domains for isometric immersions y ∈ W 1,2 is shown in [28,Corollary 1.2].
Proof We may without loss of generality assume that ω is convex. Using [15, and v k 1,∞ C. By scaling v k with η > 0 we can extend ηv k to an isometry y ([15, Theorem 7]) with ηv k = y 3 . Then, because y is an isometry: where n = y ,1 ∧ y ,2 is the normal and II (y) = (∇y) ∇n the second fundamental form of the surface y(ω). Since ∇ 2 y = 0 a.e. near x 0 , there is a neighbourhood U of x 0 covered by some leading curve γ, that is: withλ ∈ L 2 . Now, [19, Proposition 1, eq. (12)] shows that ∇y(γ(t) + sν(t)) is independent of s, hence n 3 = (y ,1 ∧ y ,2 ) 3 is also independent of s and we can subsume it into the functionλ. Setting λ(t) = −n 3 (t)λ(t)/η we obtain the representation (17). 2 Finally, we come to: Proof of Theorem 2 To recapitulate, according to (5) and (14) the linearised Kirchhoff energy is given by for v ∈ W 2,2 sh (and ∞ otherwise). We observe first that the set is non-negative and strictly convex, but it also need not consist of just one point. Note next that v is a minimiser of (18) iff ∇ 2 v(x) ∈ N for almost every x ∈ ω: On the one hand, every minimiser has finite energy and thus ∇ 2 v must be pointwise a.e. in the set {F ∈ R 2×2 sym : det F = 0}. On the other, any function F : ω → R 2×2 sym with F (x) ∈ N a.e. minimises the integrand in (18) pointwise and thus the energy.
Next we show that any two elements F, G of N are linearly independent. Indeed, by strict convexity we have for all λ ∈ (0, 1): attains a lower value here we must have det(λF + (1 − λ)G) = 0. But then it cannot be that G = ρF for any scalar ρ ∈ R or else it would hold that det(λF + (1 − λ)G) = det(λF + (1 − λ)ρF ) = C det F = 0, a contradiction. Consequently, we have in particular 0 ∈ N unless N = {0}. But in that case ∇ 2 v ≡ 0 and the proof would be concluded. Let now v ∈ W 2,2 sh be a minimiser for I lKi . Note first that ∇v cannot be constant over open sets: indeed we just saw that w.l.o.g. 0 ∈ N and consequently the condition ∇ 2 v = 0 is excluded for a minimiser on any set of positive measure. Consider then some point x 0 ∈ ω with a neighbourhood U where ∇v is not constant and use the representation (17). We have that, pointwise a.e. and over U : If κ(t) = 0, by varying s we obtain distinct, linearly dependent matrices ∇ 2 v(t, s). Because ∇ 2 v ∈ N a.e., this shows that κ(t) = 0 for a.e. t. As a consequence, γ must be constant. But then λ is also constant or again we would have points at which ∇ 2 v is linearly dependent. Since this holds locally around every x = γ(t) + sγ (t), we deduce that ∇ 2 v is constant on U and because we can cover ω in this manner, there exists F ∈ N such that ∇ 2 v ≡ F a.e. over ω. 2 4 Structure of minimisers for I θ vK for small θ The second main contribution of this work is a first study of the properties of minimisers in the interpolating regime, "close" to the linearised von Kármán model. The results in Section 3 show that the transition from spherical to cylindrical shapes occurs in the interpolated von Kármaán as the strength θ of the misfit increases. We will see that for small θ > 0 indeed there exists a unique stable branch of solutions emanating from a perfect spherical cap at θ = 0. For the sake of clarity we restrict to the prototypical model from (8): Natural subsequent steps along this line of work, which we do not take here, are to consider the regime of large values of θ and to investigate the existence of the conjectured critical value θ c , as well as to consider the full model derived in (6). 4 We recall that the existence of minimizers is guaranteed, cf. Remark 1. Without loss of generality wee assume that the barycenter of ω is 0. So with (f ) ω := 1 |ω| ω f (x) dx for a function f we in particular have (x) ω = 0. In order to avoid ambiguities (and to apply Korn's and Poincaré's inequalities) we restrict the functions w = (u, v) to lie in the Banach space and norm (u, v) X = ( u 2 1,2 + v 2 2,2 ) 1/2 . By the arguments in [11,Remark 2] working with these spaces does not lead to a loss of generality either: For an affine function g, ∇(v + g) ⊗ ∇(v + g) − ∇v ⊗ ∇v is a symmetrised gradient.
For small values of the parameter θ we have the following structural result on the set of minimizers showing the existence of a smooth branch of unique global minimisers. Let v 0 (x) = 1 2 |x| 2 − c 0 with c 0 = 1 2 (|x| 2 ) ω .
Theorem 4 There exists an ε > 0, a unique point u 0 ∈ X u and a uniquely determined C 1 map φ : [0, ε) → X such that φ(0) = (u 0 , v 0 ) and for each θ ∈ [0, ε): The proof is a direct consequence of Theorems 5 and 6 that are proved in the following two subsections. The main difficulty in obtaining a local branch of minimizers for θ 1 lies in the fact that minimisers at θ = 0 are not unique. Indeed, (u, v 0 ) ∈ argmin I 0 vK for u arbitrary, as can be readily checked. This is addressed in Subsection 4.1. The proof that in fact these minimisers are global is achieved by an application of a Taylor expansion for a carefully perturbed functional in Subsection 4.2.

A branch of solutions for θ 1
Notation In this section, the parameter θ will be explicitly included in the arguments of the functional and differentiation is understood to be with respect to the variables w = (u, v), unless otherwise stated, i.e.
We are interested in the existence and uniqueness of solutions w = (u, v) to the equation DI θ vK (u, v; θ) = 0 as a function of θ ∈ [0, ε) with I θ vK given by (8). We will in fact prove the existence of a point (u 0 , v 0 ) ∈ X such that there exists a (locally) unique function φ(θ), starting for θ = 0 at (u 0 , v 0 ), such that every φ(θ) ∈ X is a critical point for I θ vK . However, lack of uniqueness of minimisers at θ = 0, (19) will thwart what would be a natural application of the implicit function theorem. The problem manifests itself as a lack of injectivity of the first derivative at which for θ = 0 is and this vanishes at every u ∈ X u and the unique v(x) = 1 2 |x| 2 + a · x + b, a ∈ R 2 , b ∈ R, such that (v) ω = 0 and (∇v) ω = 0, i.e., v = v 0 . Because of this the equation DI θ vK (u, v; θ) = 0 in L(X, R) cannot be uniquely solvable for (u, v) ∈ X as a function of θ, even locally. Nevertheless, after some computations one can see that the problem is the presence of a leading factor θ which we can dispense with, because we may apply the implicit function theorem to the set of equivalent equations These equations are equivalent to DI θ vK (u, v; θ) = 0 for any θ > 0 and by an application of the implicit function theorem around a specific point (u 0 , v 0 ; 0) we determine the existence of a solution function φ : Then we have DI θ vK (φ(θ); θ) = 0 for θ > 0 because of the equivalence mentioned and DI θ vK (φ(0); 0) = 0 by the choice of (u 0 , v 0 ).
Proof We first define a new set of equations to solve, then show that the second derivative of I θ vK is one to one and then the conclusion is exactly that of the implicit function theorem. For brevity we write These define a scalar product and a norm in L 2 (ω; R 2×2 sym ) since Q 2 is by construction bilinear and symmetric and it is positive definite on this space. Even though Q 2 vanishes on antisymmetric matrices, during the proof we keep track of symmetrised arguments to these functions for the sake of clarity.
Step 1: Equivalent equations. From the computations leading to (20) we have: for all (ϕ, ψ) ∈ X. We observe first that, because 1 θ ∂ u I θ vK is independent of θ the right hand side makes sense even if θ = 0. Now, on the one hand, for any fixed value of θ 0 solving the system implies solving: where f : X × R → L(X, R) is given by On the other hand, solving f (u, v; θ) = 0 for θ > 0 is equivalent to solving the original problem DI θ vK (u, v; θ) = 0 as we desired.
Step 2: A zero and the derivative of f . Since we are interested in the behaviour around θ = 0, we evaluate here and We can compute a zero of f (·, ·; 0) by first considering the last term, which vanishes for all ψ ∈ X v if and only if v = v 0 . We next observe that the first term encodes the orthogonality of ∇ s u+ 1 2 ∇v 0 ⊗∇v 0 to the space of symmetrised gradients SG u := {∇ s ϕ : ϕ ∈ X u } with respect to the scalar product induced by Q 2 . The u ∈ X u realizing this is attained by projecting onto SG u , i.e.
is the orthogonal projection onto SG u given by By the Korn-Poincaré inequality this determines u 0 ∈ X u uniquely. We have then a point w 0 = (u 0 , v 0 ) such that f (u 0 , v 0 ; 0) = 0 in L(X, R).

Uniqueness and globality of minimisers
In addition to the previous local result, we can prove that the critical points found in the previous subsection are the unique global minimizers for small non zero values of the parameter θ. We do this in two steps: close to the origin (u 0 , v 0 ) of the branch of solutions, we would like to perform a Taylor expansion and use that the second differential at (u 0 , v 0 ) is "almost" positive definite.
The key idea is to slightly modify the energy by a shift and a rescaling in order to obtain derivatives as those appearing in the equivalent equations (22) of Theorem 5, thus obtaining a positive definite second derivative. We set I θ vK (ũ,ṽ) := I θ vK u 0 +ũ θ ,ṽ and then (ũ θ ,ṽ θ ) is a minimiser ofĨ θ vK if and only if (u 0 +ũ θ /θ,ṽ θ ) is a minimiser of I θ vK . In other words, if (u θ , v θ ) is a minimiser of I θ vK , thenũ θ = θ(u θ − u 0 ) andṽ θ = v θ minimiseĨ θ vK . We namew 0 the point around which we investigate the modified functional: Theorem 6 There exists θ c > 0 and a neighborhoodW ⊂ X withw 0 ∈W such that for every θ ∈ (0, θ c ), every critical point of DĨ θ vK is the unique global minimiser ofĨ θ vK . Proof We proceed in three steps. First we prove that there is some θ c > 0 such that D 2Ĩ θ vK (w) is positive definite for all θ ∈ (0, θ c ) if w −w 0 < η for some suitable η > 0 andw 0 = (0, v 0 ) as defined in (23). Then we use this to determine a neighbourhood ofw 0 where (local) minimisers ofĨ θ vK will be global by first considering points close to one such minimiser and finally those far away. We will need the first two derivatives ofĨ θ vK . For the first differential we apply the chain rule to obtain D uĨ θ vK (ũ,ṽ) = 1 θ D u I θ vK u 0 +ũ θ ,ṽ and substitute: For the second differential we can compute another directional derivative: Step 1: Local positive definiteness. We show there exist η > 0 and θ c > 0 s.t. D 2Ĩ θ vK (w) is positive definite for all θ < θ c and all w −w 0 X < η. More precisely, we even show that there exists somec > 0 such that for all θ < θ c , w −w 0 X ≤ η and (ϕ, ψ) ∈ X. Let then η > 0 be fixed and to be determined later and letw = (ũ,ṽ) ∈ X with w −w 0 X < η. We start by bringing terms together in (24): Given f, g ∈ W 1,2 (ω; R 2 ) we have, by the bounds (2) for Q 2 and Hölder (with the Sobolev embedding W 1,2 (ω; R 2 ) → L 4 (ω; R 2 )): Using this, the first and last term above can be estimated using Korn-Poincaré and Poincaré's inequality: for constants c 1 , C 1 ,C 1 > 0, where in the last step we used the assumption ṽ − v 0 2,2 < η to bound ṽ 2 2,2 by some constant independent of η 1. For the second term, use Cauchy-Schwarz for Q 2 , and the same ideas as above: Again, we used that by assumption ũ 1,2 < η and ṽ − v 0 2,2 < η.
Finally we estimate the third term in D 2Ĩ θ vK with analogous arguments and obtain (c) c 2 ψ 2 2,2 , for a c 2 > 0. Bringing the previous computations together, with a C 2 > 0 we have: from which (25) follows if θ c and η are chosen sufficiently small. From now on, we letw θ = (ũ θ ,ṽ θ ) be a critical point ofĨ θ vK with w θ −w 0 X η/3 (26) and we prove that it is in fact the unique global minimizer.
As above, the last line holds becausew θ minimisesĨ θ vK in a 2 3 η-neighbourhood of itself. 2

Discretisation of the interpolating theory
Our goal in this section is to study the qualitative behaviour of minimisers in the interpolating regime α = 3. To this end, we develop a simple numerical method to approximate minimisers and prove Γ-convergence to the continuous problem. Numerical computations are then conducted for the prototypical example from (8). We experimentally evaluate the conjectured existence of a critical value θ c > 0 for which the symmetry of minimisers is "strongly" broken. We will not provide a full theoretical analysis, but instead adduce some empirical evidence to support the claim. As can only be expected from a topic originating in structural mechanics, numerical methods for plate models are a vast field with a long history and as such a comprehensive review falls well beyond the scope of this contribution. However, it can be said that a significant portion of finite element approaches focus on the Euler-Lagrange equations. For von Kármán-like theories like our interpolating regime, these are transformed into an equivalent form in terms of the Airy stress function [20, §2.6.2]. The resulting system of equations is of fourth order and can be solved with conforming C 1 elements like Argyris or specifically taylored ones. To avoid the higher number of degrees of freedom, non-conforming methods can be used instead, 5 but a poor choice of the discretisation can suffer from locking, as briefly described in Remark 4. Some successful classical methods employ C 0 Discrete Kirchhoff triangles (DKT), but it is also possible to employ standard Lagrange elements with penalty methods [8], as we will do.
A recent line of work, upon which we heavily build in this section, is that of [4,6], where the author develops discrete gradient flows for the direct computation of (local) minimisers of non-linear Kirchhoff and von Kármán models. Γ-convergence and compactness results are also proved showing the convergence of the discrete energies to the continuous ones, as well as their respective minimisers. 6 Crucially, these papers use DKTs for the discretisation of the outof-plane displacements, allowing for a representation of derivatives at nodes in the mesh which is decoupled from function values. This enables e.g. the imposition of an isometry constraint for the non-linear Kirchhoff model, but also the computation of a discrete gradient ∇ ε projecting the true gradient ∇v ε of a discrete function v ε into a standard piecewise P 2 space. The operator ∇ ε has good interpolation properties circumventing the lack of C 1 smoothness of DKTs which would otherwise make them unsuitable to approximate solutions in H 2 . We refer to the book [5] for a systematic and mostly self-contained introduction to these methods.

Discretisation
We wish to investigate minimal energy configurations of the following functional: where (u, v) ∈ W 1,2 (ω; R 2 ) × W 2,2 (ω; R 2 ), cf. (6). We recall the representation of Q 2 derived in (12), which in particular shows that Q 2 is a strictly convex polynomial of degree 2 on R 2×2 sym × R 2×2 sym . It is extended to a convex quadratic function on R 2×2 × R 2×2 by our setting for F, G ∈ R 2×2 . We assume that ω ⊂ R is a bounded simply connected domain with Lipschitz boundary and barycenter 0. We implement (projected) gradient descent in a non-conforming method using C 0 linear Lagrange elements. The first step is to transform the problem into one of constrained minimisation reducing the order of the elements required.
Note that our assumptions on ω guarantee that Z = {∇v : v ∈ W 2,2 (ω)}. We can now use H 1 -conforming elements but, for simplicity of implementation, instead of adding the constraint into the discrete spaces to obtain a truly conforming discretisation, we add a penalty term µ ε curl z ε 2 to ensure that the solutions z ε are close to gradients.
Assume from now on that ω is a polygonal domain. For fixed ε > 0, introduce a quasi-uniform triangulation T ε of ω with triangles T of uniformly bounded diameter c −1 ε ε T cε for some c > 0 and all ε > 0 and T ∈ T ε . 7 Such a mesh is in particular said to be, in virtue of the uniform upper bound, shaperegular. We denote by N ε the set of all nodes of the triangulation. Define V ε to be the standard piecewise affine, globally continuous Lagrange P 1 finite element space S 1 (T ε ) in two dimensions: Quadrature rules will be chosen to be exact for this polynomial degree and the first integrand in the energy interpolated for this to apply by means of the interpolated quadratic function This is defined (with a slight abuse of notation) component-wise using the element-wise nodal interpolantÎ ε , defined for functions v ∈ L ∞ (ω) such that v |T ∈ C(T ) for all T ∈ T ε aŝ where ϕ z|T is the truncation by zero outside T of the global basis function ϕ z ∈ S 1 . Because this is a linear combination of truncated global basis functions, the range ofÎ ε is the spaceŜ 1 (T ε ) of discontinuous, piecewise affine Lagrange elements. In cases where the function to be interpolated is continuous, the elementwise nodal interpolant coincides with the standard nodal interpolant into the space S 1 of globally continuous, piecewise affine functions, which is defined as Notice that the shape functions ϕ z are not truncated. In order to control the error incurred by the interpolation. When working with discontinuous functions inŜ 1 , we will use the following local result. This follows from standard nodal interpolation estimates (see e.g. [17,Theorem 4.28] or [9, (4.4

Remark 3 (Scaling of the constants)
The penalty µ ε = µ(ε) needs to explode as ε → 0 in order for the functionals to Γ-converge (Theorem 7). However, large penalties negatively affect the condition number of the system, so that an adequate choice for µ ε , dependent on the mesh size ε, is required [17, p.416]. We have not explictly investigated how this requirement interacts with the Γconvergence of the functionals, but in our proof we require only that µ ε → ∞ not faster than ε −2 . In the implementation we use µ ε = ε −1/2 . Analogously, large values of the Lamé constants have a similar effect and therefore hinder convergence, so one needs to scale them to the order of the problem.
Remark 4 (Common issues with FEM for plates) Discretisations for lower dimensional theories can face complications due to the infamous locking phenomena. In a nutshell, these mean that as the thickness of the plate tends to zero, discrete solutions "lock" to stiff states of lower, or even vanishing, bending or shearing than the analytic ones. 8 Another instance of unexpected behaviour is known as the Babuška paradox [2], again a failure to converge as expected, which can happen in e.g. the Kirchhoff model when both vertical and tangential displacements are fixed at the boundaries of a polygonal domain: these so-called "hard" support constraints are not enforced in the same manner as in the continuous model because of the approximated domain.
There are two potential sources of locking in our setting: the penalty term µ , which is akin to the shear strain in Timoshenko beams, and θ. We have not obtained any a priori bounds on the error in this work, but a rigorous treatment of the problem would require estimates which are uniform in these parameters as the mesh diameter goes to zero. For the regimes studied and the geometries considered we have found the issue to be of moderate practical relevance, but it does manifest itself e.g. with more complicated domains or higher values of θ.
Finally, our simulations will not suffer from Babuška's paradox because we do not prescribe boundary conditions.
Proof Because of Lemma 2 we can substitute Q 2 for Q ε 2 in J θ ε . Also, by Lemma 3 it is enough to consider smooth functions for the upper bound. Set Step 1: Upper bound.
Let (u, z) ∈ W 1,2 (ω; R 2 ) × Z be C ∞ up to the boundary and define u ε := I ε (u), z ε := I ε (z), where I ε is the nodal interpolant of (31). Note that because u and z are smooth, we can apply standard interpolation estimates to show strong convergence in W 1,2 of these sequences towards u and z. By the compact Sobolev embedding W 1,2 → L 4 we have z ε → z in L 4 , and z ε ⊗ z ε → z ⊗ z in L 2 , so we have that A ε → A in L 2 . Since Q 2 is a polynomial of degree 2, this implies as ε → 0. By the same interpolation estimate above and the assumption on µ ε we have that µ ε curl(Î ε (z) − z) 2 0,2 = o(1) as ε → 0, and consequently Step 2: Lower bound. Let u ε , z ε ∈ V ε ⊂ W 1,2 with u ε u, and z ε z weakly in W 1,2 to u ∈ W 1,2 (ω; R 2 ), z ∈ Z. Because z ε ⊗ z ε → z ⊗ z in L 2 , we have that A ε A in L 2 . Moreover, curl z ε curl z. If linf ε→0 J θ ε = ∞, the assertion is trivial. If not, then µ ε k ω | curl z ε k | 2 dx C and curl z ε k 0,2 → 0 for a subsequence ε k → 0. But then curl z = 0. Dropping the (non-negative) curl term in J θ ε and by the weak sequential lower semicontinuity of all integrands involved (Q 2 being a convex quadratic function), we then get The final ingredient of this subsection is a proof that sequences with bounded energy are (weakly) precompact. The fundamental theorem of Γ-convergence then shows convergence of global minimisers. In order for this to work, we need to assume conditions in the space which provide Korn and Poincaré inequalities. We can do this using functions with zero mean, zero mean of the gradient or zero mean of the antisymmetric gradient as we do above, but including these conditions in the discrete spaces is not entirely trivial. Because the energies are invariant under the transformations which are factored out by taking quotient spaces as described in the sections mentioned, it is enough for our purposes to claim compactness modulo these transformations and to exclude them in the implementation via projected gradient descent.

sums:
and α j determined with line search is energy decreasing. A line search means computing the maximal α j ∈ {2 −k : k ∈ N} such that where ρ ∈ (0, 1/2) is the proverbial fudge factor.
The existence of α j > 0 is guaranteed as long as J θ ε ∈ C 2 (V 2 ε ) because then we can perform a Taylor expansion and use again (35): Remark 5 (Caveat: local and global minimisers) Even though we now know that the discrete energies correctly approximate the continuous one, as well as any global minimisers, gradient descent on each discrete problem is only guaranteed to converge to some local minimiser w ε . Lacking some means of tracking a particular w ε as ε → 0, there is not much one can do to prove that our method actually approximates the true global minimisers of I θ vK . Unless θ 1, in which case we know local minimisers to be global (cf. Theorem 6).

Experimental results
For the implementation of the discretisation detailed above, we employ the FEniCS library [1] in its version 2017.1.0. The code is available at [10] and includes the model, parallel execution, experiment tracking using Sacred [16] with MongoDB as a backend and exploration of results with Jupyter [22] notebooks, Omniboard [35] and a custom application. Everything is packaged using docker-compose for simple reproduction of the results and one-line deployment.
We set ω =B 1 (0), a (coarse) polygonal approximation of the unit disc and test several initial conditions. The space V ε has ∼7000 dofs. We implement a general Q 2 for isotropic homogeneous material with the two (scaled) Lamé constants set to those of steel at standard conditions. We apply neither body forces nor boundary conditions, but hold one interior cell to fix the value of the free constants. We compute minimisers for increasing values of θ and µ ε ∼ 1/ √ ε via projected gradient descent (onto the space of admissible functions V ε ∩ X u ) and examine the symmetry of the final solution. The choice ε −1/2 has shown to provide the fastest convergence results while keeping the violation of the constraint in the order of 10 −4 (higher penalties have the expected effect of adversely affecting convergence). We track two magnitudes as measures of symmetry: on the one hand we compute the mean bending strain over the domain and on the other, as a second simple proxy we employ the quotient of the lengths of the principal axes.
The first initial configuration is the trivial deformation y 0 ε = 0. Note that because the model is prestrained, the ground state is non-trivial and the plate "wants" to reach a lower energy state. In Figure 2 we depict the results of running the energy minimisation procedure for multiple values of θ. We further highlight the behaviour of the solution as a function of θ in Figure  3. In the first plot we compute the mean bending strains 1 |ω| ω (∇ 2 v) ii dx with i ∈ {1, 2}.
As mentioned, these act as an easy to compute proxy for the (mean) principal curvatures. We observe how as θ increases both strains decrease almost by an equal amount as the body gradually opens up and flattens out, while retaining its radial symmetry. However, around θ ≈ 86 a stark change takes place and one of the principal strains decreases while the other increases. This reflects the abrupt change of the minimiser to a cylindrical shape. We observe the same phenomenon with the quotient of the principal axes of the deformed disk in the right plot of the same Figure. The second initial condition tested is an orthotropically skewed paraboloid. Basically, a spherical cap is pressed from the sides to obtain a "potato chip". Testing this shape will highlight the effect of the initial configuration on the final curvature. We examine its strains and symmetry in Figure  5.
Again there is a critical value of θ ≈ 50 around which the shape of the minimiser drastically changes. Note however how the change is now gradual and we see intermediate shapes.