Numerical approaches for investigating quasiconvexity in the context of Morrey's conjecture

Deciding whether a given function is quasiconvex is generally a difficult task. Here, we discuss a number of numerical approaches that can be used in the search for a counterexample to the quasiconvexity of a given function $W$. We will demonstrate these methods using the planar isotropic rank-one convex function \[ W_{\rm magic}^+(F)=\frac{\lambda_{\rm max}}{\lambda_{\rm min}}-\log\frac{\lambda_{\rm max}}{\lambda_{\rm min}}+\log\det F=\frac{\lambda_{\rm max}}{\lambda_{\rm min}}+2\log\lambda_{\rm min}\,, \] where $\lambda_{\rm max}\geq\lambda_{\rm min}$ are the singular values of $F$, as our main example. In a previous contribution, we have shown that quasiconvexity of this function would imply quasiconvexity for all rank-one convex isotropic planar energies $W:\operatorname{GL}^+(2)\rightarrow\mathbb{R}$ with an additive volumetric-isochoric split of the form \[ W(F)=W_{\rm iso}(F)+W_{\rm vol}(\det F)=\widetilde W_{\rm iso}\bigg(\frac{F}{\sqrt{\det F}}\bigg)+W_{\rm vol}(\det F) \] with a concave volumetric part. This example is therefore of particular interest with regard to Morrey's open question whether or not rank-one convexity implies quasiconvexity in the planar case.


Introduction
The question whether or not rank-one convexity implies quasiconvexity in the planar case is considered one of the major open problems in the calculus of variations [58,14,54,35]. Morrey [47,48] conjectured that this is not the case in general. Of course, in order to demonstrate that these convexity properties are indeed distinct, it is sufficient to identify a single rank-one convex function which is not quasiconvex, as was done by Sverak in the non-planar case [66]. While numerous viable criteria are known for rank-one convexity, it remains highly difficult to explicitly show that a function is not quasiconvex. It is therefore common to apply numerics to the problem [22,57,10].
In this article, we discuss different numerical approaches for demonstrating the non-quasiconvexity of a given function. We will primarily apply our methods to a single example, which has been the subject of a previous article [70]: The planar isotropic energy function W + magic , which is defined via with λ max ≥ λ min > 0 as the ordered singular values of the deformation gradient F = ∇ϕ ∈ GL + (2). While it has already been shown that this function is rank-one convex but nowhere strictly elliptic and not polyconvex [70], it remains open whether or not W + magic is quasiconvex. The energy candidate (1.1) is very compelling: It emerged from the investigation of planar isotropic elastic energies with a so-called additive volumetric-isochoric split, i.e. of energies that can be written as the sum of an isochoric part depending only on the product λmax λmin and a volumetric part depending only on λ max λ min . Details of this structure and its characteristics are discussed in Section 2.1. It can be shown [70] that W + magic (F ) is a limit case in the investigation of "least rank-one convex" candidates in the family of energy functions with an additive volumetric-isochoric split. In this specific setting, it was demonstrated that the question of quasiconvexity for W + magic (F ) also determines whether or not any other rank-one convex function with an additive volumetric-isochoric split and a concave volumetric part is quasiconvex.
This work is primarily meant to serve as a guideline towards numerical investigations in the context of Morrey's conjecture. In particular, while the general methods described here are mostly applicable to a large class of energy functions, 1 we utilize a number of specific invariance properties exhibited by the particular example W + magic , as discussed in Section 2. Exploiting such properties can vastly improve the efficiency of numerical approaches and, when investigating other types of functions, it should be kept in mind that similar invariances could be identified and used to simplify the more general numerical algorithms for finding counterexamples to quasiconvexity.
The methods we describe reliably find such counterexamples for functions which are known to be nonelliptic and therefore non-quasiconvex. For the rank-one convex energy candidate W + magic , a number of microstructures with the same energy level as the homogeneous deformation were (re-)discovered. However, we were unable to demonstrate the non-quasiconvexity of W + magic with any numerical approach. Although inconclusive, these numerical results certainly suggest that the function W + magic is in fact quasiconvex and therefore not suited for answering Morrey's conjecture.

Convexity properties of energy functions
We start by recalling the classical definitions of generalized convexity properties [19,65,63].
for any domain Ω ⊂ R n with Lebesgue measure |Ω|. The energy function is strictly quasiconvex if the inequality in (1.2) is strict for all ϑ = 0.
While quasiconvexity, together with suitable growth conditions, is sufficient to ensure weak lower semicontinuity of the energy functional, it has the disadvantage of being notoriously difficult to prove or disprove directly. This led to the introduction of various sufficient and necessary conditions for quasiconvexity. Whereas polyconvexity provides a sufficient criterion for quasiconvexity, a necessary condition is given by the rank-one convexity of a function.
Definition 1.3. The energy function W : If the energy function is twice differentiable, rank-one convexity is equivalent to the Legendre-Hadamard ellipticity condition which expresses the ellipticity of the Euler-Lagrange equation Div DW (∇ϕ) = 0 corresponding to the variational problem The energy function is strictly rank-one convex if inequality (1.4) is strict.
Overall, for any W : R n×n → R ∪ {+∞} we have the well known hierarchy [4,7,19] polyconvexity =⇒ quasiconvexity =⇒ rank-one convexity , (1.7) with none of the reverse implications holding in general for n ≥ 3. The remaining question whether rank-one convexity implies quasiconvexity for n = 2 is known as Morrey's problem.
In continuum mechanics and related applications, energy functions are often more naturally defined on the group GL + (n) instead of R n×n since a deformation gradient F with non-positive determinant would imply local self-intersection. For such functions we introduce: is quasi-/poly-/rank-one convex.

Periodic boundary conditions
The definition of quasiconvexity (1.2) can be reformulated in terms of deformations with periodic boundary conditions. As will be demonstrated in Section 4, this can be helpful to find periodic microstructures numerically. In the context of planar elasticity, periodic boundary conditions have the form with a homogeneous deformation F 0 ∈ GL + (2) and a periodic superposition ϑ # ∈ W 1,∞ per (Ω). The domain Ω has to be in such a geometrical shape that R 2 can be covered by periodic replications of Ω and the superposition ϑ # must be Ω-periodic, cf. Figure 1. Seemingly, periodic boundary conditions are a more general concept than Dirichlet boundary conditions with test functions ϑ ⊂ W 1,∞ 0 (Ω; R n ) because they also allow for changes on the boundary itself. 2 However, for the notion of quasiconvexity, the two can be used interchangeably.
2 While a periodic superposition ϑ # ∈ W 1,∞ per ([a, b]) allows for more modifications of the homogeneous deformation F 0 x than a test function ϑ ∈ W 1,∞ 0 ([a, b]) , the average rate of change of ϕ(x) for any domain Ω ⊂ R n with Lebesgue measure |Ω| such tht R 2 can be covered by periodic replications of Ω. The energy function is strictly quasiconvex if the inequality in (1.9) is strict for all ϑ = 0.
In the context of planar elasticity, it is sufficient to consider Ω = [0, 1] 2 , and without loss of generality we can assume ϑ # (0, 0) = (0, 0). The periodic replication of Ω to cover R 2 implies that the values of ϑ # coincide on opposite edges, i.e. every point at the boundary belongs to two separate unit squares (four at the corners of the square). Thus we can write Ω-periodicity as and

Previous results related to Morrey's conjecture
Morrey's conjecture has long been considered one of the most important open questions in the calculus of variations, and the remaining problem of the planar case has been the subject of extensive research.
The  function W : is a homogeneous polynomial of degree four; here, λ i denote the signed singular values 3 of F . It has been shown [23] that the function is rank-one convex if γ ∈ [0, 4 √ 3 ≈ 2.309], but polyconvex only if γ ≤ 2. This result can be extended to the more general class At present, it is not known whether this expression is quasiconvex for γ ∈ (2, 4 √ 3 ]. However, extensive numerical calculations suggest that it is quasiconvex [22,21]. In the same works, Dacorogna et al. study several rank-one convex functions, including an example by Ball [6,20]: Together with an example by Aubert [2], the Dacorogna-Marcellini function has been the first example given in the literature of a planar function which is rank-one convex but not polyconvex. Many planar functions used in the context of Morrey's conjecture have the structure W (F ) = g( F 2 , det F ) for which additional numerical optimization is available [27,26] or are composed of polynomials up to the degree four [31,9].
Guerra and da Costa [29] recently employed a systematic numerical approach towards the question of Morrey's conjecture. Their findings suggest that in the planar case, rank-one convexity implies N -wave quasiconvexity for N ≤ 5. Since any function which is N -wave quasiconvex for all N ∈ N is quasiconvex [29, Proposition 3.6] (cf. [64]), these results provide some evidence for the conjecture that rank-one convexity indeed implies quasiconvexity for planar energies.
Previous applications of machine learning in nonlinear elasticity, as we consider in Section 4.1, have mostly focused on the energy function itself [25,41]. The application to deformation functions presented here is based on the concept of physics-informed neural networks, which have recently been employed for finding approximate solutions to various partial differential equations [59,34].

Exploitable properties of functions
Before employing numerical methods to investigate whether a function is quasiconvex, it is generally useful to identify invariances and similar properties of the specific function that may allow to improve the efficiency of the numerical approach. In the following, we will focus on the particular energy function W + magic , which exhibits three important properties that can be exploited numerically: isotropy, scaling invariance and the specific form of a volumetric-isochoric split.

The volumetric-isochoric split
The energy W + magic emerged from the investigation of the family of planar isotropic energies W : GL + (2) → R with an additive volumetric-isochoric split . (2.1) We will motivate both the additional structure one can achieve with this type of energy functions as well as the candidate W + magic itself. By [46,Lemma 3.1], energies of the type (2.1) can be written as where λ 1 , λ 2 > 0 denote the singular values of F and h, f are real-valued functions. In nonlinear elasticity theory, energy functions with an additive volumetric-isochoric split are widely used, primarily to model the behaviour of slightly compressible materials [15,33,53,50].
For a further representation of W , we introduce the (nonlinear) distortion function or outer distortion where . denotes the Frobenius matrix norm with for all a > 0 , R ∈ SO(2) .
We can then express every conformally invariant energy W on GL + (2) as W (F ) = h(K(F )) = Ψ(K(F )) with h , Ψ : [1, ∞) → R [46]. Note that in general, In a previous paper [70] we motivated the reduction to the case f (det F ) = c log det F for arbitrary additive volumetric-isochoric split energies with a newly developed rank-one convexity criterion [69]. More specifically, we showed that if there exists a rank-one convex energy function with an additive volumetric-isochoric split that is not quasiconvex, then we can find such a function in the set as well. It is therefore sufficient to consider M * instead of the general class of volumetric-isochoric split energies when discussing Morrey's conjecture.

Scaling invariance
In general an arbitrary energy W ∈ M * is neither simply scaling invariant, i.e. W (αF ) = W (F ) for all α ∈ R and F ∈ GL + (2), nor tension-compression symmetric, i.e. W (F ) = W (F −1 ) for all F ∈ GL + (2). However, M * provides additional invariance properties that hold for rank-one convexity and quasiconvexity which we will use to simplify numerical calculations in the following sections by reducing the number of deformation gradients F ∈ GL + (2) that we must consider.
Lemma 2.1. Let W ∈ M * be twice differentiable. Then the ellipticity domain of W is scaling invariant, i.e. a cone: if W is elliptic at F 0 ∈ GL + (2), then W is elliptic at αF 0 for every α > 0.
Proof. Let α > 0. The isochoric part W iso (F ) = h λ1 λ2 of W is conformally invariant by definition, which implies W iso (αF ) = W iso (F ), and therefore for all H ∈ R 2×2 . For the volumetric part W vol (F ) = c log det F we calculate Together with (2.8) this implies for all F, H ∈ R 2×2 with rank(H) = 1 and all α > 0, which implies the scaling invariance of the ellipticity domain of W ∈ M * .
Proof. Due to the isotropy of every W ∈ M * , and since the singular values of F and F T are identical, W (F ) = W (F T ) and therefore thus ellipticity at F and F T are equivalent, cf. [43]. Moreover, Lemma 2.1 states that the ellipticity domain of W is scaling invariant, i.e. ellipticity at F implies W is elliptic for all αF with α > 0. Combining both lemmas and choosing α = 1 det F > 0, we find that ellipticity at Cof F would imply ellipticity at 1 det F (Cof F ) T = F −1 . Therefore, it remains to show that W is elliptic at Cof F .
In the planar case, the singular values of F and Cof F are identical 4 and thus W (Cof F ) = W (F ) for all F ∈ GL + (2) due to the isotropy of the energy. Furthermore, because of the simple shape of the cofactor matrix in the planar case Cof Note carefully that these properties do not hold for dimension n > 2. Thus which completes the proof because rank(H) = 1 implies 5 rank(Cof H) = 1.
for any domain Ω ⊂ R 2 , then the energy is quasiconvex at αF 0 for all α > 0.
Proof. We show that the so-called energy gap i.e. the difference between the energy of ϕ(x) = F 0 x + ϑ(x) and the homogeneous solution ϕ 0 (x) = F 0 x, is scaling invariant. Note that quasiconvexity at F 0 , Hence the eigenvalues of B and Cof B = (Cof F )(Cof F ) T and thus the singular values of F and Cof F are identical. 5 The rank of a matrix M ∈ R n×n is equal to the maximal size of nonzero minors. The cofactor Cof M is formed by implies that the energy gap is always non-negative. We write F = ∇ϕ and compute Remark 2.4. With Lemmas 2.1 and 2.3 we showed that both rank-one convexity and quasiconvexity are scaling invariant for any energy W ∈ M * . Regarding Morrey's question for planar isotropic energies with volumetric-isochoric split we can therefore assume det F 0 = 1 without loss of generality. More specifically, for arbitrary W ∈ M * we just need to prove the rank-one convexity for all F 0 ∈ GL + (2) with det F 0 = 1 to obtain rank-one convexity at all F ∈ GL + (2). Likewise, it is sufficient to only check for quasiconvexity starting with a homogeneous deformation x → F 0 x with det F 0 = 1 in place of arbitrary F ∈ GL + (2), cf. Figure 2. [ellipticity/quasiconvexity] on the curve det F = λ 1 λ 2 = c (black) implies [ellipticity/quasiconvexity] for the corresponding cone. Right: For an arbitrary volumetric-isochoric split energy this invariance is lost.

The least rank-one convex candidate W + magic (F )
Continuing with the class M * , i.e. with energy functions with c log det F as volumetric part, we focus on positive c > 0. As multiplying a function by a scalar does not change its convexity behavior, we can consider c = 1 without loss of generality and assume that the isochoric part is convex. 6 Thus we consider the class Using sharp rank-one convexity conditions [69], it is possible to identify "least" rank-one convex candidates by searching for functions that satisfy those conditions by equality. Surprisingly [70], it is possible to show that the question of quasiconvexity reduces to the single energy function Hence, if W + magic were quasiconvex, then every function in the class M * + and thus every rank-one convex planar isotropic energy function with an additive volumetric-isochoric split whose isochoric part h(t) is convex would be quasiconvex as well. A first analytic observation (see Appendix A) could not conclusively answer this question, but opens the possibility to interesting microstructures having the same energy value as the homogeneous deformation. This motivates to proceed numerically and try to falsify the quasiconvexity of W + magic .

The classical finite element approach
A first way to show that a given energy density is not quasiconvex is the finite element method. For this, we discretize the displacement field ϑ by Lagrange finite elements [16]. We use triangle grids and first-order finite elements only. That way, the deformation gradient and hence the integrand are piecewise constant and the hyperelastic energy Ω W (∇ϕ) dx can be computed without quadrature error which is important in view of exactly calculating the energy gap (2.14). Our implementation is based on the Dune libraries [11,61].

Testing for quasiconvexity
We perform tests on two domains: The square [−1, 1] 2 and the unit disc B 1 (0). Both are filled with coarse triangle grids as shown in Figure 3. For ease of implementation, we approximate the boundary of the disc by six quadratic arcs (dashed lines). The finite elements grids are then constructed by uniform refinement of the coarse grids. For the disc grid, new boundary vertices are placed not at edge midpoints but on the curved arcs approximation the boundary. The final grids consist of 16641 vertices and 32768 elements for the square and 12481 vertices and 24576 elements for the disc. Due to the scaling invariance (the results of Section 2.2) it is sufficient to test with an F 0 such that det F 0 = 1, thus due to isotropy we consider only with arbitrary a > 0. For the result shown here we select a = 2. We minimize the hyperelastic energy Ω W + magic (F ) dx with a trust-region algorithm. These algorithms have been thoroughly studied in the literature, and they can be shown to always converge to stationary points of the energy [17]. As trust-region methods are descent methods, maximizers and saddle-points are not attractive points, and convergence is therefore typically towards (local) minimizers only. Trust-region methods perform sequences of quadratic minimization problems with convex inequality constraints. We use a trust-region defined in terms of the maximum norm, and therefore the convex constraints are a set of independent bound constraints. For the quadratic bound-constrained minimization problem we then use a monotone multigrid method [42] as suggested in [60]. Such methods achieve multigrid convergence rates even for bound-constrained problems. We solve each inner problem until the maximum-norm of the correction drops below 10 −5 . The large but sparse tangent matrices are computed using the ADOL-C algorithmic differentiation system [71].
When looking for global minimizers with a descent method, the question of initial iterates is of central importance. As shown exemplary in the next section, when the energy is not quasiconvex, imperfections caused by the finite-precision arithmetic are sufficient to drive the system towards microstructures even starting from the homogeneous configuration. For the particular energy W + magic , however, this did not lead to any energy decrease. We obtained the same negative results for some other "obvious" initial iterates, such as random perturbations of the homogeneous state. More involved constructions of non-homogeneous initial iterates are described in Sections 3.2 and 5.

Non-elliptic microstructure
In the following, we introduce several additional ideas to search for microstructures with energy levels below the homogeneous state with more adept methods.
In order to better understand the shape of a possible microstructure, we ensured its existence by considering slightly modified problems. We start with the weakened energy which is not rank-one convex anymore but satisfies lim . We are interested in the resulting microstructures especially if they do not appear to be simple laminations caused by the loss of ellipticity [5,24,45]. Any local minimizer found for W c can then be used as an initial deformation for minimizing W + magic again with the hope to maintain the non-homogeneous structure and thereby disprove quasiconvexity of the energy W + magic itself. While our material W c indeed shows microstructures that contain a lamination structure, we observe radially symmetric contracting regions as well, cf. Figure 4. However, we already know that deformations of this kind cannot lower the energy value of W + magic (see Appendix A) and thus they are not suited for finding a new microstructure by using them as an initial configuration for minimizing W + magic (F ) again. This is confirmed by direct numerical experiments: when letting the trust-region solver of the previous section minimize W + magic starting from the configurations shown in Figure 4, all we obtain is convergence to the homogeneous state. Therefore, we continue with an alternative numerical experiment to produce more convoluted microstructures.
For values slightly larger than c = 1 we observe that all three circles contract in a radially symmetric fashion, cf. Figure 5. Again, we note that Appendix A shows that radial symmetric contracting deformations have the same energy value for the limit case c = 1. Increasing the weighting of log det F by raising c results in lower energy values compared to the homogeneous deformation. For higher values of c, the microstructure becomes more convoluted but keeps its contraction structure, cf. Figure 5.
We note that both microstructures are primarily located inside these circles (where ellipticity is lost), while for the rest of the material with c = 1, the deformation remains mostly homogeneous. In particular, the borders of the inner circles maintain their shape to a certain extent, even though we do not impose additional internal boundary conditions to ensure this. We interpret these observations as a first indicator towards our energy candidate W + magic being quasiconvex since the response of W + magic is "stable" toward the assumed non-elliptic perturbation in the circle.

Deep neural networks
One may conjecture that minimizing W + magic in a finite element space fails to find microstructure because each finite element coefficient influences only a very local part of the deformation function. Additionally, when the strain energy density lacks quasiconvexity (which we consider possible for W + magic ), FEM-methods generally fail [44] due to the non-uniqueness of the solution, i.e. the microstructures. In this chapter, we experiment with an alternative approach where this relationship is more global.

Physics-informed neural networks
We use a numerical scheme which is based on deep neural networks as an ansatz for solving partial differential equations, an idea also referred to as physics-informed neural networks [59,34]. In principle, this is similar to the ansatz constructed classically using finite element functions. However, deep neural networks generally lead to highly nonlinear and more efficient (in terms of number of parameters) approximation spaces with considerable approximation properties even for low numbers of coefficients.
In the following, we consider periodic deformations only, cf. Section 1.2. Consider the following ansatz for the periodic superposition: where F ω f , G ωg : [0, 1] → R 2 and H ω h : [0, 1] 2 → R 2 are deep neural networks [62] parameterized by parameters ω f , ω g , and ω h , respectively. The ansatz is intentionally constructed this way to identically satisfy the periodic boundary conditions (1.10). The neural network architectures are given by repeated composition of successive high-dimensional linear and nonlinear transformations as Here, L i→j ω ,k , k = 1, · · · , 5 ( is a placeholder for 'f ', 'g', and 'h') denotes the k th -linear transformation parameterized by the set of weights and biases ω , The linear transformations are interleaved with element-wise nonlinear transformations R(·) = tanh(·).
Since each layer of a neural network is a differentiable operation, the gradient of the superposition field ϑ # can be computed using the chain rule. This is efficiently implemented using an automatic differentiation engine [12]. Note that, unlike numerical differentiation by, e.g., finite differences, the gradients computed via automatic differentiation are analytically exact.
For numerical integration of the strain energy density over the domain Ω we discretize the domain with a uniform grid {(x α 1 , x β 2 ) | α, β = 1, . . . , N } of N × N points. For a given ∇ϑ #,ω , the total strain energy is approximated via the trapezoidal integration rule as where the integration weights ξ(α) are Whether the trapezoidal integration rule over-or underestimates the integral depends on whether the integrand is convex or concave, respectively, in the interval of the integration. However, if the integrand exhibits an inflection point (which is also observed here), over-/underestimation of the integral cannot be guaranteed with this rule.
The optimal parameters ω = {ω f , ω g , ω h } of the neural networks are then obtained as minimizers of the total energy The minimization problem is solved iteratively using Adam [39], an efficient first-order gradient-based stochastic optimization method. The derivatives of the objective function with respect to the parameters ω are computed via automatic differentiation again. Following the minimization, the periodic superposition field ϑ # = ϑ #,ω is given by (4.1). The Adam optimizer was used for 2000 iterations with an initial learning rate of 10 −3 which was decayed by a factor of 0.1 after the 700 th , 1400 th , and 1800 th iterations. The numerical scheme was implemented in PyTorch [55]. Figure 6 illustrates the representative microstructures obtained via the numerical scheme for F 0 equal to diag 3, 1 3 and diag 10, 1 10 on a grid of resolution N = 128. For both values of F 0 the microstructure has the form of a smooth laminate and its energy seems to equal the one of the corresponding homogeneous deformation up to machine precision 7 which motivates the search for a precise form of the analytical solution.

Smooth laminates
Guided by the numerical findings of the previous section, it turned out that the smooth laminates can be explained analytically as well. Let Ω = [0, 1] 2 be the unit square and consider W + magic (F ) = λmax λmin + 2 log λ min with the ordered singular values λ max ≥ λ min of F ∈ GL + (2). For any homogeneous deformation ϕ 0 (x) = F 0 x with F 0 = diag(a 1 , a 2 ) and a 1 ≥ a 2 > 0, the elastic energy Proof. Let ϑ # be as in (4.7). We find that This implies λ max = a 1 + f # (x 1 ) and λ min = a 2 . Periodic boundary conditions as defined by (1.10) imply f (0) = f (1) as well as f (0) = f (1). Thus We also show explicitly that any ϕ as defined in Lemma 4.1 is indeed an equilibrium point of I by direct computation. For the corresponding Euler-Lagrange equations we must compute the first Piola-Kirchhoff stress S 1 (F ) = DW + magic (F ). We start by considering the restriction to deformations of the type (4.7) as an a priori constraint and verify that the single resulting reduced Euler-Lagrange equation holds: Since in the class of deformations of the type (4.7), all functions are stationary points of this constrained problem. This is a necessary condition for the general problem [67,68]. It remains to show that any such function is also a stationary point of the full Euler-Lagrange equations S 1 (F ) = DW + magic (F ) as well. For this we identify and start with the first derivative of the nonlinear distortion function K(F ): Furthermore, using the notation from (4.10), we make use of (4.14) Altogether, we find Next, we insert deformations of the type (4.7) as an a posteriori constraint [67]: Thus we arrive at (4.17) The full Euler-Lagrange equations are therefore satisfied since Note that the case f # ≡ 0 corresponds to the homogeneous solution and that this is the only superposition which also complies with Dirichlet boundary conditions given by F 0 . As deformations of the type (4.7) only allow for a smooth displacement in one coordinate direction, we refer to them as smooth laminates, cf. Figure  6. Remark 4.2. We can construct a continuous map remaining in the class of smooth laminates (i.e. which only contains stationary points of W + magic ) from one candidate f # to the identity f # ≡ 0 and from there to another solution f # at constant energy value I(ϕ). Therefore, neither the homogeneous solution nor any of these smooth laminates are stable. A similar argument holds for the radially symmetric deformations considered in [70]. Remark 4.3. If we allow ϑ # to be non-smooth (cf. Figure 7) we can obtain various (classical) lamination patterns as well. For example, as a limit of a smooth laminate microstructure of the type (4.7), we may consider which yields the simple two-phase laminate Recall that we can only construct laminations that satisfy the constraint f # (x) + a 1 ≥ −a 2 , cf. Lemma 4.1.

Relaxation to non-gradient fields with a Curl-based penalty
While the smooth laminates discussed in Section 4 as well as the contracting deformations in [70] (cf. Appendix A) allow for non-homogeneous deformations whose energy level for W + magic is as low as the homogeneous one, we have yet failed to find strictly lower energy values.
For a new attempt, we expand our set of possible deformations ϕ. To this end, we extend the energy functional from gradient fields F = ∇ϕ to more general matrix fields P , but control the distance of P to the set of compatible mappings (i.e. gradient fields) with the Curl 2D operator. This approach is commonly used when dealing with local dislocations in gradient plasticity with plastic spin [52,49] and with relaxed micromorphic continua [51]. We expect to obtain new microstructures numerically and, in order to gain additional insight into our original variational problem, we will observe the material behavior as a function of the weight parameter on the penalty term Curl 2D P . The hope is that when such new microstructures are used as initial iterates, the optimization algorithm of Section 3 will converge to an energy level below the homogeneous one.
More specifically, instead of the classical minimization problem in the space W 1,2 (Ω, R 2 ), i.e.
with W + magic as in (2.17), we consider in the larger space H(Curl). Here, the vector field τ is the unit tangent vector to ∂Ω 8 and L c ∈ R + is a penalty parameter. The planar operator Curl 2D is discussed in Appendix B. Note that the mapping P : Ω ⊂ R 2 → R 2×2 does not need to be a gradient field, i.e. P may be incompatible, but that inf I 2 (P ) ≤ inf I 2 (∇ϕ) = inf I 1 (ϕ) and lim Lc→∞ I 2 (P ) = +∞ if P is not a gradient field , since we are considering contractible domains only, and on such domains Curl 2D P = 0 if and only if P = ∇ϕ for some weakly differentiable ϕ.
We minimize the functional I 2 of (5.2) numerically by approximating the matrix field P by piecewise polynomial functions P n such that each row of P n is a first-order Nédélec finite element of the first kind [40] which are elements of the space H(Curl) by construction. As the domain, we use the unit ball Ω = B 1 (0) and the same grid as in Section 3.2. The algorithm used to minimize the discretized functional is the same trust-region multigrid algorithm as in Section 3.2 as well.
We consider the minimization problem (5.2) for different parameters L c and choose the boundary deformation We must note that for the range of parameters L c considered here, the trust-region algorithm of Section 3 would not converge to a stationary point. Rather, due to the severe degeneracy of the energy function (5.2) for low L c , the solver would get stuck eventually. In such situations, the trust-region control would decrease the trust-region radius further and further, without ever finding an acceptable new iterate. In view of the finite-precision arithmetic used in actual simulation, this is not a contradiction to the trust-region theory, which claims that such a new iterate will always be found.
All figures show the last step computed. For L c = 0.5, 1, 2 we observe different non-trivial microstructures with significantly lower energy value than the homogeneous solution, cf. Figure 8. The larger L c , the closer the energy approaches the energy of the homogeneous state.
After calculating P from I 2 (P ) → min by the above numerical method, we construct a deformation field ϕ : Ω → R 2 with a deformation gradient close to P by computing in the space of first-order Lagrange finite elements. While the compatible deformations ϕ all look similar to their corresponding incompatible matrix fields P , cf. Figure 9, the energy value I 1 ( ϕ) is always higher then the homogeneous one I 1 (ϕ 0 ). Starting the minimization algorithm for I 1 from ϕ always results in the homogeneous configuration.  As a further attempt to reach low values of I 2 (P ), we also started the optimization from a non-compatible matrix field instead of the homogeneous deformation gradient P = F 0 we used for the numerical methods before. For this we chose a checkerboard pattern with squares of size 1 b × 1 b and alternately set the values F 1 = (1 − δ)F 0 and F 2 = (1 + δ)F 0 with δ = 0.5, cf. Figure 10.
This pattern has a lower energy value without considering the regularization term Curl P 2 as it only activates the concave volumetric part log det F of W + magic (F ). As shown in Figure 11, we indeed arrive at different microstructures. However, the compatible deformations obtained from minimizing ∇ϕ − P 2 L 2 (Ω) Figure 10: Checkerboard pattern used as initial iterate. The pattern is not a gradient and combines an inhomogeneous determinant with a constant distortion K.
again have higher energy values than the homogeneous state. Furthermore, minimizing the original energy I 1 from there lead back to the homogeneous configuration once again.

Gradient Young measures and laminates
More recently, the problem of Morrey's conjecture has been approached from the point of view of gradient Young measures and laminates [36,37,38,29,28,32]. In the following, we give a brief overview of this alternative approach and its relation to the optimization methods used above.
Throughout this section, let Ω = [0, 1] 2 denote the unit square. 9 In order to positively answer Morrey's conjecture, we would need to find a rank-one convex function W : R 2×2 → R that is not quasiconvex. More explicitly, this task can be stated as follows: (I) Find a rank-one convex function W : R 2×2 → R, a Lipschitz mapping ϕ : Ω → R 2 and F 0 ∈ R 2×2 with ϕ(x) = F 0 .x for all x ∈ ∂Ω such that In particular, such a full solution to Morrey's problem requires both an energy function W : R 2×2 → R and a mapping ϕ : Ω → R 2 .
The classical approach to this problem, as followed in the previous sections, consists of choosing a plausible candidate W first before trying to find a mapping ϕ that satisfies (6.1) in order to prove the non-quasiconvexity of W . However, recent investigations of Morrey's conjecture have often taken a slightly different point of view, emphasizing the search for a suitable mapping ϕ instead. More specifically, these approaches try to establish the existence of a Lipschitz mapping ϕ : Ω → R 2 such that the pushforward measure induced by ∇ϕ violates a Jensen-type inequality. Figure 11: Two compatible deformations (for L c = 0.5 and L c = 2, respectively) starting from a checkerboard pattern for F 0 = diag √ 2, 1 √ 2 . This pattern remains visible for the first calculation because as explained in the text, the trust-region algorithm did not find a local minimizer for I 2 (P ) → min.
This change in perspective is mainly based on the seminal results by Kinderlehrer and Pedregal [36,37,38] on the relation between quasiconvexity and gradient Young measures. In particular, these results allow for Morrey's conjecture (I) to be rephrased in terms of properties of probability measures.

Equivalent formulations of Morrey's conjecture
First, we consider the following rephrasing, which is obtained by a simple substitution on the left-hand side of (6.1).
This phrasing of Morrey's problem in terms of probability measures 10 is closely related to its formulation in terms of homogeneous gradient Young measures. The exact relation between pushforward measures of gradients and homogeneous gradient Young measures can be established by the Averaging Theorem [37, Theorem 2.1], which directly implies the following lemma as a corollary.
Then ν ϕ is a homogeneous gradient Young measure on Ω with barycenter ν ϕ = F 0 .
By virtue of Lemma 6.1, if (II) can be solved, so can the following problem: (III) Find a rank-one convex function W : R 2×2 → R and a homogeneous gradient Young measure ν such that In order to see that (III) is also sufficient for (and thus equivalent to) solving (II), recall that for any homogeneous gradient Young measure ν there exists a sequence (ϕ k ) k∈N ⊂ W 1,∞ (Ω; R 2 ) with affine linear boundary values induced by F 0 = ν such that for any continuous function W : R 2×2 → R. In particular, this implies R 2×2 W (A) dν ϕ k (A) < W (ν) = W (F 0 ) for sufficiently large k if ν satisfies (6.2).
Finally, we can rephrase (III) by employing the notion of a laminate. In the planar case, a laminate can be characterized as a probability measure ν on R 2×2 such that the Jensen-type inequality holds for any rank-one convex function W : R 2×2 → R [56]. Since ν is a homogeneous gradient Young measure if and only if (6.3) holds for any quasiconvex energy W according to the Kinderlehrer-Pedregal Theorem [37,38], the following problem is equivalent to (III): (IV) Find a homogeneous gradient Young measure ν which is not a laminate.
Note that (IV) seems to make no reference to any energy function W . In practice, however, (approximately) rank-one convex energy functions need to be applied to a given measure ν in order to numerically establish whether it is a laminate [29]. On the other hand, it is often obvious by construction that ν is a homogeneous gradient Young measure. In particular, due to Lemma 6.1, this is the case if ν is obtained as the pushforward measure with respect to the gradient of a mapping ϕ : Ω → R 2 .

The numerical search for non-laminate gradient Young measures
A specific numerical method for finding non-laminate homogeneous gradient Young measures numerically has been suggested by Guerra et al. [29]. This approach is based on selecting a dense subset K(Ω) ⊂ W 1,∞ (Ω; R 2 ) such that each ϕ ∈ K(Ω) induces a discrete-valued gradient field ∇ϕ. Then the following problem needs to be solved numerically: Both the barycenter ν ϕ and the energy value R 2×2 W (A) dν(A) can be easily computed numerically for the measure ν ϕ corresponding to any ϕ ∈ K(Ω) and a given energy function W on R 2×2 . However, in order to demonstrate that a given measure ν ϕ is not a laminate, it is necessary to find a rank-one convex energy W for which (6.3) is violated. The numerical approach suggested by Guerra et al. [29] consists of generating a large number of (approximately) rank-one convex energy functions by computing the rank-one convex envelope for a specific class of functions based on earlier considerations byŠverák [66].

Application to W + magic
In contrast, based on our results in [70], we conjecture that the functional W + magic defined in 2.17 is a good candidate energy to detect non-laminate measures. Therefore, we do not need to perform the computationally expensive calculations of multiple rank-one convex envelopes. Instead, we directly include the class K(Ω) of deformations in our search for a counterexample to the quasiconvexity inequality (1.2). Since the energy values of the inhomogeneous deformations ϕ ∈ K(Ω) and of the corresponding homogeneous deformations given as the energy at the barycenter W (F 0 ) = W (ν ϕ ) are easily computed numerically (without the need of FEM), we can thereby test W + magic with a large number of additional deformations with periodic boundary conditions. We call these inhomogeneous deformations ϕ ∈ K(Ω) combined laminates.
For a given homogeneous deformation gradient F 0 and given number N of laminates to combine, cf. (6.4), our implementation selects random parameter values for η i , ξ i , c i with i ∈ {1, · · · , N }, validating that det F > 0 for all such combinations, and computes the energy value of the resulting combined laminate. Figures 12  and 13 show examples of such deformations with N = 3 and F 0 = 1; the "phases" of the superimposed deformations, i.e. the local values of h ( x, η i + c i ), are indicated by (+/−).
Our numerical tests focused on the case N ≥ 5. 11 After testing slightly more than one million combinations with a routine written in Python, where we started with different 1,10] and tried N = 4, 5, 6, 7, we were once more unable to obtain an energy level below W + magic (F 0 ). As with our previous approaches, we always find non-trivial microstructures when we change our energy function to be non-rank-one convex. Figure 14, for example, shows such a microstructure for the energy W c defined in (3.1) with a modified volumetric part c log det F in place of W + magic ; for c > 1, this method once more finds random configurations with lower energy than the homogeneous state.

Discussion
We have presented several different numerical approaches to check for quasiconvexity of a given function W . • In Section 3, we demonstrated a classical finite element approach that can find easily microstructures if we perturb the energy candidate to be slightly non-rank-one convex. In addition, we showed a method of disturbing the homogeneous structure of the solution by modifying the energy values on a subdomain and computing the microstructure resulting from this material inhomogeneity.
• In Section 4, under the assumption of periodic boundary conditions, we introduced a numerical scheme that is based on deep neural networks and thereby discovered a new non-trivial microstructure (smooth laminates) with the same energy value as the homogeneous deformation for the considered energy function, which we then investigated analytically.
• For the relaxation technique considered in Section 5, we extended our numerical calculations from gradient fields F to general matrix field P by introducing a penalty term based on Curl P . This approach resulted in various non-compatible fields with lower energy value than the energy of the homogeneous deformation, even for a quasiconvex energy candidate. The fields found this way were then used as new starting configurations for the finite elements approach. This, unfortunately, did never lead to energies below the one of the homogeneous state.
• In Section 6, we discussed an alternative numerically straightforward way of checking for quasiconvexity connected to the theory of gradient Young measures rather than to optimization. Several rank-one laminations were combined so that the resulting deformation remained piece-wise homogeneous and their energy value were compared to the homogeneous deformation. Again, we only found lower values for a non-rank-one convex energy density.
We tested all these methods with the energy W + magic from [70]. If we change this energy to be non-elliptic and thus non-quasiconvex, the presented methods were all able to produce non-trivial microstructures, i.e. deformations that are neither homogeneous nor a simple first-order laminate. This demonstrates their viability as numerical tests for quasiconvexity. On the other hand, while non-quasiconvexity of an energy W can conceivably be proven by numerical methods (since identifying a single suitable deformation would be sufficient), no amount of numerical testing can be considered an actual proof that a given function is quasiconvex. However, if all the described methods fail to yield a deformation energetically more optimal than the homogeneous one -as was the case for the energy W + magic (F ) = λmax λmin + 2 log λ min considered here -then this should be interpreted as a strong indication that the function is indeed quasiconvex and therefore not a viable candidate for answering Morrey's conjecture. where Div    P11 · · · P1n . . . . . .
Since in the planar case, curl can be seen as a rotated div, i.e.
Because of the involved nonlinear form of W + magic (P ) it is, however, not clear whether or not P ∈ L 2 (Ω) holds. 15 15 For the modified energy it is possible to achieve W 0 (F ) ≥ c + F q .