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ármán to linearised von Kármán theories. While explicit formulae are available in the linearised regimes, the von Kármán theory turns out to be critical and a phase transition from cylindrical (as in linearised Kirchhoff) to spherical or saddle-shaped (as in linearised von Kármán) configurations is observed there. We analyse this behaviour with the help of a whole family (IvKθ)θ∈(0,∞)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(\mathcal{I}^{\theta }_{\mathrm{vK}})_{\theta \in (0,\infty )}$\end{document} of effective von Kármán functionals which interpolates between the two linearised regimes. We rigorously show convergence to the respective explicit minimisers in the asymptotic regimes θ→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\theta \to 0$\end{document} and θ→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\theta \to \infty $\end{document}. Numerical experiments are performed for general θ∈(0,∞)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\theta \in (0,\infty )$\end{document} which indicate a stark transition in a critical region of parameters θ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\theta $\end{document}.


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 B B. Schmidt bernd.schmidt@math.uni-augsburg.de M. de Benito Delgado m.debenito.d@gmail.com 1 Universität Augsburg, 86135 Augsburg, Germany 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.
There is, by now, a considerable body of work on thin plates with pre-strain. We mention, with no attempt to be exhaustive, [12,24,25,28,29,38,39] and refer to our companion paper [11] for further discussion and references. As our primary interest is to model multilayers, our set-up is quite different from those of the aforementioned contributions where the pre-strain depends only on the in-plane variables. In the bending dominated regime, models in which both pre-strain and material parameters were allowed to vary in the thin film direction were already discussed in [38,39]. The recent paper [29] also considers a variety of scaling regimes for thickness dependent pre-strains, while our companion paper generalises those results by allowing for thickness dependent material properties, in particular, for mismatching equilibria in multilayers with different elastic moduli. Except for finite bending configurations in [38,39], however, it appears that little attention has been devoted so far to the issue that is our main focus here-namely, how the choice of the scaling regime influences the geometry of the energy-minimizing configurations.
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 prestrains in very thin layers tend to cause cylindrical shapes whereas smaller pre-strains in thicker layers lead to spherical caps, [13-15, 22, 30, 37]. 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 [19,32,35,40]. On the macroscopic scale, recent experiments on self-folding of two dimensional elastomer polydimethylsiloxane figures into rather complex three dimensional structures are described in [13].
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, α > 2.
The particular case α = 2 with a misfit of the order h of the aspect ratio has been investigated in [7,38,39]. The appropriate plate theory is the nonlinear Kirchhoff theory (in the finite bending regime) and energy minimisers 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 [16]-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 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 thelimit of the 3d nonlinearly elastic energy, rescaled by h 1−2α , cf. [11]. This energy rescaling is tailored so as to capture the effect of the misfit to leading order: Larger rescalings would lead to infinite energy contributions caused by the pre-strain whereas at smaller rescalings the effect of the misfit becomes negligible. For a minimiser (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, nondegenerate parabolas (infinitesimal parts of an elliptical cap or a saddle) are energy minimisers. Only in the latter case, however, the minimiser 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 minimisers. 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 minimisers 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 minimisers in their dependence on θ . Interestingly, our results seem to indicate a stark change of material response in a critical parameter region, showing a symmetry breaking 'phase transition' from a nearly spherical cap to an approximate cylinder.
The occurrence of the three characteristic types of optimal shapes in dependence of the energy scaling, and thus the magnitude of deflections in relation to the aspect ratio, seems to be in line with experimental results (see the above cited literature). In particular, this is exemplified by the isotropic bilayers with isotropic misfit described in [13]. Here, even though the mismatch of equilibria of the two layers is fixed, different scaling regimes can be explored by varying the aspect ratio h. At constant thickness, this simply amounts to testing samples of varying lateral dimensions. Indeed, 'thick' films with rather large aspect ratio are seen to form shallow spherical caps whereas specimens at very small aspect ratio are observed to prefer cylindrical configurations, cf. [13]. A transition between caps and (parts of) cylinders is found with deflections being of the order of the film thickness (corresponding to the von Kármán regime). In fact, at the experimentally observed 'bifurcation points' the measured curvature is close to the aspect ratio, see [13,Fig. 7].

Outline
We begin by recalling our main results from [11] in order to provide the appropriate plate theories in Sect. 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 minimisers. We then discuss the linearised regimes α ∈ (2, 3) and α > 3 as well as the asymptotic von Kármán regimes θ → 0 and θ → ∞ in Sect. 3. The structure of minimisers for small θ is investigated in Sect. 4. Finally, Sect. 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 byconvergence. We then determine the effective (homogenised) elastic moduli and corresponding quadratic energy densities of the plates in terms of the moments of the pointwise elastic constants of the layers.
The existence of minimisers of (6) (7) and (8) follows by a standard application of the direct method or, in the setting of [11], as a direct consequence of -convergence and compactness.

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 sym → 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 To see this, fix ∈ R 2×2 and note that for all x ∈ R 2 \ {0} . t would imply that (tI − )x = 0 in contradiction to having at most two eigenvalues. Expanding the square we get We now introduce a scalar β 0 and vectors b 1 , b 2 ∈ R 3 by setting and we let B 1 , B 2 ∈ R 2×2 sym be the 2 × 2 symmetric matrices corresponding to b 1 and b 2 . Let Q 2 be given as in (4). 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 . 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 Sect. 5, we will study the energy (15) as a function of θ . 2. As a second example we consider a BGaAs/InGaAs bilayer consisting of an InGaAs layer on top of which a BGaAs film is grown epitaxially and where the thickness of the upper layer is approximately 0.8 times the thickness of the lower layer, [32,39]. The linearised energy within the two individual layers is given by with suitable positive constants C 11 , C 44 , C 12 (only depending on t ), see below. As a consequence we obtain For InGaAs (i.e., − 1 2 < t < 1 18 ), respectively, BGaAs (i.e., 1 18 The lattice constants are 0.58031 nm for InGaAs and 0.56313 nm for BGaAs, which corresponds to a misfit of approximately 3%. The misfit tensor B h = B is isotropic, sopossibly after rescaling the in-plane coordinates-it can be chosen in such a way that the first moment B 1 vanishes. With these properties we get Choosing for definiteness B(t) ≈ −1.49I for −1/2 < t < 1/18 and B(t) ≈ 1.56 · I for 1/18 < t < 1/2 of order 1, we see that our theory is relevant for a layer of aspect ratio h and a scaling parameter α whenever h α−1 is of the order 1%. We also compute B 2 = 44.9I and, using (10), (13) and (14), γ = β 0 = 547, F 0 = 4.59I , E 0 = 0 and thus arrive at the explicit formulae

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 Sect. 1 that we are primarily interested 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,Sect. 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 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,Sect. 6]. We allow for a general bounded Lipschitz domain ω in these theorems.

Theorem 1
The minimisers of I lvK , Eq. (8), 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
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 Sect. 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 If #N ≥ 3, 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 minimiser, or #N = 2 and there are precisely two minimisers, or N is an affine ellipse and to each 'winding direction' Re, e ∈ S 1 , there is a unique curvature λ = λ(e) such that   (17).

Proof of Theorem 1 By (8) and (12)
Proof of Theorem 3 a) is immediate from [11,Theorems 7,10,11]. b) directly follows from [11,Theorems 7,8,9] if ω is convex. For general ω first note that the compactness result in [11,Theorem 7] does not use convexity, (12) and (13). Then forū = u + θ −1/2 u we have by (14) 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 -convergence 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.
As for Theorem 2, it is straightforward to see that v as defined in the theorem is a minimiser 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. [33]):

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.
Finally, we come to: Proof of Theorem 2 To recapitulate, according to (6) and (14) the linearised Kirchhoff energy is given by for v ∈ W 2,2 sh (and ∞ otherwise). We observe first that the set N = argmin{Q * is non-negative and strictly convex, but it also need not consist of just one point. Note next that v is a minimiser of (19) iff ∇ 2 v(x) ∈ N for almost every x ∈ ω: Every minimiser has finite energy and thus det ∇ 2 v = 0 a.e. So x → 1 2 x F x for any F ∈ N is a minimiser. But then a general v ∈ W 2,2 sh (ω) is a minimiser if and only if ∇ 2 v ∈ N a.e.
Next we show that any two elements F , G of N are linearly independent. Indeed, by strict convexity we have for all λ ∈ (0, 1): Hence λF + (1 − λ)G / ∈ N or else F , G would not be minimisers. Because Q * 2 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 − λ) . 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 (18). 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 ω.

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. For the sake of simplicity we restrict to the prototypical model from (15): where we allow for θ ∈ [0, ∞). Here the θ -dependent term θ 1/2 (∇ s u + 1 2 ∇v ⊗ ∇v) decouples from the second contribution ∇ 2 v and moreover the average of the misfit B 1 = 1/2 −1/2 Q 2 (t,B(t)) dt vanishes, which considerably simplifies the following computations. In fact the extension to the general case is not obvious. The results in Sect. 3 show that the transition from spherical to cylindrical shapes occurs in the interpolated von Kármán regime 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.
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 a critical bifurcation value of θ , as well as to consider the full model derived in (7). 2 We recall that the existence of minimisers is guaranteed, cf. Remark 1. Without loss of generality we assume that the barycenter of ω is 0. So with (f ) ω := 1 |ω| ω f (x) dx for a function f we have in particular (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 . (Here and in the sequel we abbreviate · k,p,ω = · W k,p (ω) and also drop ω if the domain is clear.) By the arguments in [11,Remark 2] working with these spaces does not lead to a loss of generality either: For an affine function Theorem 4 There exists an ε > 0, a unique point u 0 ∈ X u and a uniquely determined 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 minimisers for θ 1 lies in the fact that minimisers at θ = 0 are not unique. Indeed, as can be readily checked. This is addressed in Sect. 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 Sect. 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 as a function of θ ∈ [0, ε) with I θ vK given by (15). 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, (20) 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 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 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 ).

Theorem 5
There exists an open set W in X, an ε > 0, a point u 0 ∈ X u such that w 0 = (u 0 , v 0 ) ∈ W and a uniquely determined for all w ∈ W and θ ∈ [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 (21) 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 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 obtain 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., where π : L 2 (ω; R 2×2 sym ) → L 2 (ω; R 2×2 sym ) 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 Finally, we compute d dε | ε=0 f (u 0 + εϕ 2 , v 0 + εψ 2 ; 0)[ϕ 1 , ψ 1 ] to have the derivative of f : Step 3: The map F : X → L(X, R) is an isomorphism. Note first that the map defines a scalar product in X, with positive-definiteness following from Korn-Poincaré's and Poincaré's inequality. Then we can write F as where we definedπ := ∇ −1 s • π , a continuous map from L 2 (ω; R 2×2 sym ) to X u . The Riesz representation for F (ϕ 2 , ψ 2 ) in L(X, R) is then (ϕ 2 +π((∇v 0 ⊗ ∇ψ 2 ) sym ), 1 12 ψ 2 ) and the map is clearly an isomorphism in X, with continuity for ψ 2 →π((∇v 0 ⊗ ∇ψ 2 ) sym ) following from the continuity ofπ and the Sobolev embedding W 1,2 → L 4 .

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 minimisers for small non zero values of the parameter θ . 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 (23) 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 θ minimisẽ 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 inW ofĨ θ 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 (24). 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 that 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 (25): 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 (26) follows if θ c and η are chosen sufficiently small. From now on, we letw θ = (ũ θ ,ṽ θ ) be a critical point ofĨ θ vK with and we prove that it is in fact the unique global minimiser.
As above, the last line holds becausew θ minimisesĨ θ vK in a 2 3 η-neighbourhood of itself.
Proof of Theorem 4 It suffices to observe that, for ε > 0 sufficiently small, the set W in Theorem 5 can be chosen so small that (θ (u − u 0 ), v) ∈W for all (u, v) ∈ W , whereW is the set obtained in Theorem 6, and to recall that (u θ , v θ ) is a minimiser of I θ vK if and only if (θ (u θ − u 0 ), v θ ) is a minimiser ofĨ θ vK .

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 (15). We experimentally evaluate the conjectured existence of critical parameter values θ > 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 [21, §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 tailored ones. To avoid the higher number of degrees of freedom, non-conforming methods can be used instead, 3 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 [9], 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. 4 Crucially, these papers use DKTs for the discretisation of the out-of-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: (7). 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 E, F ∈ R 2×2 . We assume that ω ⊂ R 2 is a bounded simply connected domain with Lipschitz boundary and barycenter 0. We implement (projected) gradient descent in a nonconforming 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.
The goal is to solve: 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 [18, p. 416]. We have not explicitly 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. 6 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.

-convergence of the Discrete Energies
The first step in the proof that J θ ε → J θ is dispensing with the interpolation operators for numerical integration: due to the good properties ofÎ ε , we can assume that we work with the true integrals Q 2 instead of Q ε 2 : Lemma 2 (Numerical integration) Let u ε , z ε ∈ W 1,2 (ω; R 2 ) be uniformly bounded in W 1,2 and let Q ε 2 =Î ε • Q 2 as above. Let A ε := θ 1/2 (∇ s u ε + 1 2 z ε ⊗ z ε ), −∇z ε . Then, as ε → 0: Proof By the local interpolation estimate Lemma 1: . Now, the first term is simply 1 + |A ε | 0,2,ω ≤ |ω| 1/2 + A ε 0,2,ω which is uniformly bounded since z ε ⊗ z ε 0,2 = z ε 2 0,4 z ε 2 1,2 , and for the second we use that both ∇ s u ε and ∇z ε are piecewise constant so that for i = 1, 2, We plug this into the preceding computation to obtain The last two norms being uniformly bounded, we conclude: The second step is, as usual, to ensure that we can focus on smooth functions for simplicity in the construction of the upper bound: 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.
We claim now that Î ε (|z ε | 4 ) − |z ε | 4 0,1 = O(ε). Indeed, by the local interpolation estimate (Lemma 1) and Hölder's inequality for integrals and for sums: Fig. 2 Final configurations after gradient descent starting with a flat disk viewed from the top. From left to right, top to bottom: θ = 1, 81, 91 and 150. Color represents the magnitude of the displacements |w|, from blue at its minimum (near the center for small θ ) to red at the maximum (near the boundary only) (Color figure online) (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 w 0 ε = 0. Note that because the model is pre-strained, the ground state is non-trivial and the plate "wants" to reach a lower energy state. In Fig. 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 Fig. 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", see Fig. 4. Testing this shape will highlight the effect of the initial configuration on the final curvature. We examine its strains and symmetry in Fig. 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.
Focusing in the critical region and comparing the energy values we see that the configurations with flat initial condition have a slightly lower energy, around 0.2% far from the transition and up to 0.8% around it. Yet we believe this to be a lack of precision in the discretisation around the critical region requiring a more detailed and thorough investigation which has to be postponed to future work.   Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.