Higher Sobolev Regularity of Convex Integration Solutions in Elasticity

In this article we discuss quantitative properties of convex integration solutions arising in problems modeling shape-memory materials. For a two-dimensional, geometrically linearized model case, the hexagonal-to-rhombic phase transformation, we prove the existence of convex integration solutions $u$ with higher Sobolev regularity, i.e. there exists $\theta_0>0$ such that $\nabla u \in W^{s,p}_{loc}(\mathbb{R}^2)\cap L^{\infty}(\mathbb{R}^2)$ for $s\in(0,1)$, $p\in(1,\infty)$ with $0<sp<\theta_0$. We also recall a construction, which shows that in situations with additional symmetry much better regularity properties hold.


Introduction
In this article we are concerned with the detailed analysis of certain convex integration solutions, which arise in the modeling of solid-solid, diffusionless phase transformations in shape-memory materials. We seek to precisely analyze the regularity properties of these constructions in a simple, two-dimensional, geometrically linear model case. Shape-memory materials undergo a solid-solid, diffusionless phase transition upon temperature change (see e.g. [Bha03] and the references given there): In the high temperature phase, the austenite phase, the materials form very symmetric lattices. Upon cooling down the material, the symmetry of the lattice is reduced, the material transforms into the martensitic phase. Due to the loss of symmetry, there are different variants of martensite, which make these materials very flexible at low temperature and give rise to a variety of different microstructures. Mathematically, it has proven very successful to model this behavior variationally in a continuum framework as the following minimization problem [BJ89]: minΩ W (∇y, θ)dx. (1) Here Ω ⊂ R n is the reference configuration of the undeformed material. The mapping y : Ω → R n describes the deformation of the material with respect to the reference configuration. It is assumed to be of a suitable Sobolev regularity. The function W : R n×n × R → [0, ∞) denotes the energy density of a given deformation gradient M ∈ R n×n at a certain temperature θ ∈ R. Due to frame indifference, W is required to be invariant with respect to rotations, i.e.
Modeling the behavior of shape-memory materials, the energy density further reflects the physical properties of these materials. In particular, it is assumed that at high temperatures θ > θ c the energy density W has a single minimum (modulo SO(n) symmetry), which (upon normalization) we may assume to be given by the SO(n) orbit of α(θ)Id, where α : R → (0, ∞) with α(θ c ) = 1 (c.f. [Bal04]). This is the (austenite) energy well at temperature θ. Upon lowering the temperature below a critical temperature θ c , the function W displays a (discrete) multi-well behavior (modulo SO(n)): There exist finitely many matrices U 1 (θ), . . . , U m (θ) ∈ R n×n + , m ∈ N, such that The matrices U j (θ) represent the variants of martensite at temperature θ < θ c and are referred to as the (martensite) energy wells. At the critical temperature θ = θ c both the austenite and the martensite wells are energy minimizers. In the sequel, we assume that θ < θ c is fixed, so that only the variants of martensite are energy minimizers. We seek to study the quantitative behavior of minimizers for energies of the type (1). Here we make the following simplifications: (i) Reduction to the m-well problem. Instead of studying the full variational problem (1), we only focus on exact minimizers. Restricting to the low temperature regime, this implies that we seek solutions to the differential inclusion ∇y ∈ m j=1 SO(n)U j (θ), for some θ < θ c . (ii) Small deformation gradient case, geometric linearization. We further modify (2) and assume that ∇y is close to the identity. This allows us to linearize the problem around this constant value (c.f. Chapter 11 in [Bha03]). Instead of considering (2), we are thus lead to the inclusion problem e(∇u) := ∇u + (∇u) T 2 ∈ {e 1 , . . . , e m }.
The symmetrized gradient e(∇u) represents the infinitesimal deformation strain associated with the displacement u, which is defined as u(x) := y(x) − x (with slight abuse of physical convention in the sequel we do not distinguish between the deformation and the displacement in our use of language and will simply refer to both as a "deformation"). The symmetric matrices e 1 , . . . , e m ∈ R n×n are the exactly stress-free strains representing the variants of martensite. While this linearizes the geometry of the problem (by replacing the symmetry group SO(n) by an invariance with respect to the linear space Skew(n)), the differential inclusion (3) preserves the inherent physical nonlinearity, which arises from the multi-well structure of the problem. (iii) Reduction to two dimensions and the hexagonal-to-rhombic phase transformation. In the sequel studying an as simple as possible model case, we restrict to two dimensions and a specific two-dimensional phase transformation, the hexagonal-to-rhombic phase transformation (this is for instance used in studying materials such as Mg 2 Al 4 Si 5 O 18 or Mg-Cd alloys undergoing a (three-dimensional) hexagonal-to-orthorhombic transformation, c.f. [CPL14], [KK91], and also for closely related materials such as Pb 3 (VO 4 ) 2 , which undergo a (three-dimensional) hexagonal-to-monoclinic transformation, c.f. [MA80a], [MA80b], [CPL14]). From a microscopic point of view, the hexagonal-to-rhombic phase transformation occurs, if a hexagonal atomic lattice is transformed into a rhombic atomic lattice. From a continuum point of view, we model it as solutions to the differential inclusion u : R 2 → R 2 , 1 2 (∇u + (∇u) T ) ∈ K a.e. in Ω, where Ω ⊂ R 2 is a bounded Lipschitz domain and We note that all the matrices in K are trace-free, which corresponds to the (infinitesimal) volume preservation of the transformation. We note that the set K is "large" (its convex hull is a two-dimensional set in the threedimensional ambient space of two-by-two, symmetric matrices, c.f. Lemma 2.7). In the sequel, we study the problem (4), (5) and investigate regularity properties of its solutions.
1.1. Main result. The geometrically linearized hexagonal-to-rhombic phase transformation is a very flexible transformation, which allows for numerous exact solutions to the associated three-well problem (4) with different types of boundary data. Here the simplest possible solutions are so-called simple laminates, for which the strain is a one-dimensional function e(∇u)(x) = f (x · n) for some vector n ∈ S 1 and for which f (x · n) ∈ {e (i1) , e (i2) } a.e. in Ω, i 1 , i 2 ∈ {1, 2, 3} and i 1 = i 2 , i.e. e(∇u) only attains two values. The possible directions of these laminates, as given by the vector n ∈ S 1 are (up to sign reversal) six discrete values, which arise as the symmetrized rank-one directions between the energy wells: For each i 1 , i 2 ∈ {1, 2, 3} with i 1 = i 2 there exists (up to sign reversal and exchange of the roles of a i1,i2 and n i1,i2 and renormalization) exactly one pair (a i1,i2 , n i1,i2 ) ∈ R 2 \ {0} × S 1 with the property that e (i1) − e (i2) = a i1,i2 n i1,i2 := 1 2 (a i1,i2 ⊗ n i1,i2 + n i1,i2 ⊗ a i1,i2 ).
The possible vectors are collected in Lemma 2.13. In addition to these "simple" constructions, there are further exact solutions to the three-well problem associated with the hexagonal-to-rhombic phase transformation, e.g. there are patterns involving all three variants as depicted in Figures 24 and 25 in the Appendix (Section 7).
In the sequel, we study solutions to the hexagonal-to-rhombic phase transformation with affine boundary conditions, i.e. we consider u ∈ W 1,p loc (R 2 ) with p ∈ (2, ∞] such that u : R 2 → R 2 , ∇u = M a.e. in R 2 \ Ω, 1 2 (∇u + (∇u) T ) ∈ K a.e. in Ω. (6) Here we investigate the rigidity/ non-rigidity of the problem by asking whether it has non-affine solutions: (Q1) Are there (non-affine) solutions to (6) with M ∈ R 2×2 ? Clearly, a necessary condition for this is that e(M ) ∈ conv(K). Using the method of convex integration, Müller andŠverák [MŠ99] (c.f. also the Baire category arguments of [Dac07], [DM12]) constructed multiple solutions to related differential inclusions, displaying the existence of a variety of solutions to the problem. Noting that these techniques are applicable to our set-up of the three-well problem, ensures that for any M with e(M ) ∈ intconv(K) there exists a non-affine solution to (6). In general these convex integration solutions are however very "wild" in the sense that they do not possess very strong regularity properties (c.f. [DM95b]). As our inclusion (6) is motivated by a physical problem, a natural question addresses the relevance of this multitude of solutions: (Q2) Are all the convex integration solutions physically relevant? Or are they only mathematical artifacts? Is there a mechanism distinguishing between the "only mathematical" and the "really physical" solutions? Guided by the physical problem at hand and the literature on these problems, natural criteria to consider are surface energy constraints and surface energy regularizations. For our differential inclusion these translate into regularity constraints and lead to the question, whether unphysical convex integration solutions have a natural regularity threshold. Here an immediate regularity property of solutions to (6) is that e(∇u) ∈ L ∞ (R 2 ). With slightly more care, it is also possible to obtain solutions with the property that u ∈ W 1,∞ loc (R 2 ). However, prior to this work it was not known whether these solutions can enjoy more regularity, i.e. whether for instance there are convex integration solutions with ∇u ∈ W s,p (Ω) for some s > 0, p ≥ 1. Motivated by these questions, in this article, we study the regularity of a specific convex integration construction and obtain higher Sobolev regularity properties for the resulting solutions: Theorem 1. Let Ω ⊂ R 2 be a bounded Lipschitz domain. Let K be as in (5) and let M ∈ R 2×2 be such that e(M ) := M +M T it is natural to expect that convex integration constructions deteriorate for matrices M with e(M ) approaching the boundary of conv(K). The precise dependence on the behavior towards the boundary however is less intuitive. In this context, it is interesting to note that the regularity threshold θ 0 > 0 does not depend on the distance to the boundary of K, but rather on the angle, which is formed between the initial matrix e(M ) and the boundary of K. This is in agreement with the intuition that the larger the angle is, the better the convex integration algorithm becomes, as it moves the values of the iterations, which are used to construct the displacement u, further into the interior of K. In the interior of K it is possible to use larger length scales, which increases the regularity of solutions. Whether this dependence is necessary in the value of the product of sp or whether the product sp should be independent of this and only the value of the corresponding norm should deteriorate with a smaller angle, is an interesting open question. We remark that in the special case of additional symmetries it is possible to construct much better solutions. An example is given in the appendix for the case M = 0 (c.f. also [Pom10] and [CPL14]). It is an important and challenging open question, whether it is possible to exploit further symmetries and thus to construct further solutions with these much better regularity properties.
1.2. Literature and context. A fascinating problem in studying solid-solid, diffusionless phase transformations modeling shape-memory materials is the dichotomy between rigidity and non-rigidity. Since the work of Müller andŠverák [MŠ99], who adapted the convex integration method of Gromov [Gro73], [EM02] and Nash-Kuiper [Nas54], [Kui55] to the situation of solid-solid phase transformations, and the work of Dacorogna and Marcellini [Dac07], [DM12], it is known that under suitable conditions on the convex hulls of the energy wells, there is a very large set of possible minimizers to (1) (c.f. also [SJ12] and [Kir03] for a comparison of these two methods). More precisely, the set of minimizers forms a residual set (in the Baire sense) in the associated function spaces. However, in general convex integration solutions are "wild"; they do not enjoy very good regularity properties. This has rigorously been proven for the case of the geometrically nonlinear two-well problem [DM95b], [DM95a], the geometrically nonlinear three-well problem in three dimensions (the "cubic-to-tetragonal phase transformation") [Kir98], [CDK07] and (under additional assumptions) for the geometrically linear six-well problem (the "cubic-to-orthorhombic phase transformation") [Rül16]. In these works it has been shown that on the one hand convex integration solutions exist, if the deformation gradient is only assumed to be L ∞ regular. If on the other hand, the deformation gradient is BV regular (or a replacement of this), then solutions are very rigid and for most constant matrices M the analogue of (6) does not possess a solution. Thus, convex integration solutions cannot exist at BV regularity for the deformation gradient; at this regularity solutions are rigid. At L ∞ regularity they are however flexible and a multitude of solutions exist. Similarly as in the related (though much more complicated) situation of the Onsager conjecture for Euler's equations [SJ12], [DLSJ16] or the situation of isometric embeddings [CDLSJ12], it is hence natural to ask whether there is a regularity threshold, which distinguishes between the rigid and the flexible regime. It is the purpose of this article to make a first, very modest step into the understanding of this dichotomy by analyzing the W s,p regularity of a (known) convex integration scheme in an as simple as possible model case.
1.3. Main ideas. In our construction of solutions to the differential inclusion (6) we follow the ideas of Müller andŠverák [MŠ99] (in the version of [Ott12]) and argue by an iterative convex integration algorithm. For the hexagonal-to-rhombic transformation this is particularly simple, since the laminar convex hull equals the convex hull of the wells and since all matrices in the convex hull are symmetrized rank-one-connected with the wells (c.f. Lemma 2.7). As a consequence it is possible to construct piecewise affine solutions (in the language of [Kir03], Chapter 4). This simplifies the convergence of the iterative construction drastically. It is one of the reasons for studying the hexagonal-to-rhombic phase transformation as a model problem. Yet, in spite of the (relative) simplicity of obtaining convergence of the iterative construction to a solution of (6) and hence of showing existence, substantially more care is needed in addressing regularity. In this context we argue by an interpolation result (c.f. Theorem 2 and Proposition 5.5): While our approximating deformations u k : R 2 → R 2 are such that the BV norms of the iterations increase (exponentially), the L 1 norm of their difference decreases exponentially. If the threshold θ 0 > 0 is chosen appropriately, the W s,p norm for 0 < sp < θ 0 is controlled by an interpolation of the BV and the L 1 norms, which can be balanced to be uniformly bounded. To ensure this, we have to make the iterative algorithm quantitative in several ways: (i) Tracking the error in strain space. In order to iterate the convex integration construction, it is crucial not to leave the interior of the convex hull of K in the iterative modification steps. In qualitative convex integration algorithms, it suffices to use errors, which become arbitrarily small and to invoke the openness of intconv(K). As the admissible error in strain space is however coupled to the length scales of the convex integration constructions (c.f. Lemma 3.3) and as these in turn are directly reflected in the solutions' regularity properties, in our quantitative algorithm we have to keep track of the errors in strain space very carefully. Here we seek to maximize the possible length scales (and hence the error) without leaving intconv(K) in each iteration step. This leads to the distinction of various possible cases (the "stagnant", the "push-out", the "parallel" and the "rotated" case, c.f. Notation 3.6, Definition 3.10 and Algorithm 3.8).
In these we quantitatively prescribe the admissible error according to the given geometry in strain space. (ii) Controlling the skew part without destroying the structure of (i). Seeking to construct W 1,∞ solutions, we have to control the skew part of our construction. Due to the results of Kirchheim, it is known that this is generically possible (c.f. [Kir03], Chapter 3). However, in our quantitative construction, we cannot afford to arbitrarily change the direction of the rank-one connection, which is chosen in the convex integration algorithms, at an arbitrary iteration step. This would entail BV bounds, which could not be compensated by the exponentially decreasing L 1 bounds in the interpolation argument. Hence we have to devise a detailed description of controlling the skew part (c.f. Algorithm 3.11). (iii) Precise covering construction. In order to carry out our convex integration scheme we have to prescribe an iterative covering of our domain by constructions, which successively modify a given gradient. As our construction in Lemma 3.3 relies on triangles, we have to ensure that there is a class of triangles, which can be used for these purposes (c.f. Section 4). In particular, we have to quantitatively control the overall perimeter (which can be viewed as a measure of the BV norm of ∇u k ) of the covering at a given iteration step of the convex integration algorithm. This crucially depends on the specific case ("rotated" or "parallel"), in which we are in.
1.4. Organization of the article. The remainder of the article is organized as follows: After briefly collecting preliminary results in the next section (interpolation results, results on the convex hull of the hexagonal-to-rhombic phase transition), in Section 3 we begin by describing the convex integration scheme, which we employ.
Here we first recall the main ingredients of the qualitative scheme (Section 3.1) and then introduce our more quantitative algorithms in Sections 3.2-3.3.2. As this algorithm crucially relies on the existence of an appropriate covering, we present an explicit construction of this in Section 4. Here we also address quantitative covering estimates for the perimeter and the volume. The ingredients from Sections 3 and 4 are then combined in Section 5, where we prove Theorem 1 for a specific class of domains. In Section 6 we explain how this can be generalized to arbitrary Lipschitz domains. Finally, in the Appendix, Section 7, we recall a symmetry based construction for a solution to (6) with M = 0 with much better regularity properties.

Preliminaries
In this section we collect preliminary results, which will be relevant in the sequel. We begin by stating the interpolation results of [CDDD03], on which our W s,p bounds rely. Next, in Section 2.2 we recall general facts on matrix space geometry and in particular apply this to the hexagonal-to-rhombic phase transformation and its convex hulls.
2.1. An interpolation inequality and Sickel's result. Seeking to show higher Sobolev regularity for convex integration solutions, we rely on the characterization of W s,p Sobolev functions. Here we recall the following two results on an interpolation characterization [CDDD03] and on a geometric characterization of the regularity of characteristic functions [Sic99]: Theorem 2 (Interpolation with BV, [CDDD03]). We have the following interpolation results: (i) Let p ∈ [2, ∞) and assume that 1 q = 1−θ p + θ for some θ ∈ (0, 1). Then (ii) Let p ∈ (1, 2] and let 1 q = 1−θ p + θ for some θ ∈ (0, 1). Let further (θ 1 , q 1 ) ∈ (0, 1) × (1, ∞) be such that , for some τ ∈ (0, 1), where 1− denotes an arbitrary positive number slightly less than 1. Then, Before proceeding to the proof of Theorem 2, we present an immediate corollary of it: For functions, which are "essentially" characteristic functions, we obtain the following unified result: Corollary 2.1. Let u : R n → R n be a function, such that u L ∞ (R n ) < ∞ and |u(x)| ≥ c 0 > 0 for a.e. x ∈ supp(u).  Figure 1. For functions, which are "essentially" characteristic functions in the sense that condition (9) of Corollary 2.1 holds, we obtain the interpolation inequality (10), which is valid in the whole coloured region in the figure (green and blue). Here the blue region is already covered in Theorem 2 (i). In order to also obtain the green region, we have to be able to simplify the statement of (8), which in general is only valid for functions, which are "essentially" characteristic functions. In our application of Corollary 2.1 (c.f. Proposition 5.5, 5.8), we will restrict to the region to the left of the dashed line. Remark 2.2 shows that having a bound for the product of the right hand side of (10) for a specific value (θ 0 , 1) already allows to deduce a bound for all exponents (θ, q) on the associated dashed line connecting (θ 0 , 1) with (0, ∞).
Proof of Corollary 2.1. By virtue of Theorem 2 (i) and equation (7), it suffices to consider the regime, in which p ∈ (1, 2). In this case the statement follows from a combination of equation (8) and the fact that for functions satisfying (9) we have for 1 < p 1 ≤ r ≤ p 2 , r −1 = σp −1 1 + (1 − σ)p −1 2 and σ ∈ (0, 1). We postpone a proof of (11) to the end of this proof, and observe first that it indeed suffices to show (11) to conclude the claim of (10). To this end, we note that the exponents in (8) obey the relation This in turn is a consequence of the three identitites Here we note that τ 1−θ ∈ (0, 1). Hence (11) (applied to r = p, p 1 = 1+, p 2 = 2 and σ = τ 1−θ ) together with (8) yields the claim of (10). It thus remains to prove (11). To this end, we observe that for any r ∈ [1, ∞] where for abbreviation we have set C 1 := u L ∞ (R n ) . With this we infer This concludes the argument.
After this discussion, we come to the proof of Theorem 2: Proof of Theorem 2. If p ≥ 2, the interpolation result is a special case of Theorem 1.4 in [CDDD03] (where in the notation of [CDDD03] we have chosen s = 0, t = θ): Indeed, for γ < 1 − 1 n and (s, p) satisfying (s − 1)p * 1 n = γ − 1 with p * being the dual exponent of p, the estimate in Theorem 1.4 from [CDDD03] reads We note that in the setting of Theorem 2 the estimate (12) is applicable, as in the notation of [CDDD03] and with dimension n we have that γ := − p p−1 1 n + 1 = 1 − 1 n − 1 p−1 1 n < 1 − 1 n , which implies the validity of (7). The simplification from (12) to (7) is then a consequence of the facts that [BM01]), • and for p ≥ 2 the embedding L p (R n ) → B 0 p,p (R n ) is valid (Theorem 2.41 in [BCD11]).
This concludes the argument for (i). To obtain (ii), we combine (i) with an additional interpolation inequality, which becomes necessary, as the inclusion L p (R n ) → B 0 p,p (R n ) is no longer valid for p ∈ (1, 2). Hence, we rely on the following interpolation estimate (c.f. Lemma 3 in [BM01]) which is valid for −∞ < s 0 < s 1 < ∞, 0 < q 0 , q 1 ≤ ∞, 0 < p 0 , p 1 ≤ ∞, 0 < τ < 1 with Here the spacesF s r,l denote the (modified) Triebel-Lizorkin spaces from [BM01]. The main advantage of the estimate (13), which goes back to Oru [Oru98], is that there are no conditions on the relations between l, l 1 , l 2 in this estimate. In particular, we can choose l 1 = 2, l 2 = p 0 and l = r. Using that •F s r,2 (R n ) = L s,r (R n ) for s ∈ R, 1 < r < ∞ and that for this range we can simplify (13) to yield which is valid for 0 < s 1 < ∞, 1 < p 0 < ∞, 1 < p 1 < ∞ with We apply (14) with p 0 = 1+, s = θ, r = q and (s 1 , p 1 ) = (θ 1 , q 1 ) lying on the boundary of the interpolation region from (i) (c.f. the blue region in Figure 1), i.e.
Although this theorem provides good geometric intuition and could have been used as an alternative means of proving Theorem 1, we do not pursue this further in the sequel, but postpone its discussion to future work.
Remark 2.2. We note that the estimate (15) in Theorem 3 yields a condition on the product θq > 0, while, at first sight, Theorem 2 and Corollary 2.1 pose a restriction on θ, q individually. As we are dealing with bounded (or even characteristic) functions, we however observe that it is also possible to obtain an analogous condition on the product θq in Theorem 2 and Corollary 2.1: Indeed, assume that u ∈ L ∞ (R 2 ) is such that for some θ 0 ∈ (0, 1) the product is bounded. Then, we claim that for q ∈ (1, ∞),θ := θ 0 q −1 and for p = 1 −θ θ also the product is bounded. To derive this, we first observe that the L ∞ bound for u allows us to infer that for any p ∈ (1, ∞) As a consequence, we deduce that Here we have made use of the specific choices of exponents from (17) and the boundedness of u, which allowed us to invoke (18). This concludes the argument for the claim. Thus, relying on the bound (19), we infer that given a bound on (16), we obtain that for all exponents q,θ, p from (17) Here we applied Theorem 2 (or Corollary 2.1), for which we noted that the respective exponents are admissible. On the one hand, this is the desired analogue of the condition from Theorem 3 and allows us to obtain a whole family of W θ,q bounds for u, where θq < θ 0 . On the other hand, it shows that although p = 1 is not admissible in Theorem 2 and Corollary 2.1, for our purposes, it still suffices to consider the case p = 1 and to prove a control for (16), which then gives the full range of expected exponents in the form of the estimate (20).
Remark 2.3 (Fractal packing dimension). Following Sickel [Sic99], Proposition 3.3 (c.f. also [JM96], Theorem 2.2) we remark that for a characteristic function its W s,p regularity has direct consequences on the packing dimension (c.f. [JM96], [Mat99]), which we denote by dim P , of its boundary: If for some set E ⊂ R n its characteristic function χ E satisfies χ E ∈ W s,p (R n ) for some s > 0 and 1 ≤ p < ∞, then dim P (S δ (∂E)) ≤ min{n, n − sp + δ}. 2.2. Matrix space geometry. Before discussing our convex integration scheme, we recall some basic notions and properties of the hexagonal-to-rhombic phase transformation, which we will use in the sequel. We begin by introducing notation for the symmetric and antisymmetric part of two matrices. 2.2.1. Lamination convexity notions. Relying on the notation from Definition 2.4, in the sequel we discuss the different notions of lamination convexity. Here we distinguish between the usual lamination convex hull (defined by successive rankone iterations) and the symmetrized lamination convex hull (defined by successive symmetrized rank-one iterations): Definition 2.5 (Lamination convex hull, symmetrized lamination convex hull).
We define the following notions of lamination convex hulls: We refer to U lc as the laminar convex hull of U and to L k (U ) as the laminates of order at most k.
Here a b := 1 2 (a⊗b+b⊗a). We refer to U lc sym as the symmetrized laminar convex hull of U and to L k sym (U ) as the symmetrized laminates of order at most k. (iii) We denote the convex hull of a set U ⊂ R m by conv(U ).
Remark 2.6. We note that if U ⊂ R n×n or U ⊂ R n×n sym is (relatively) open, then also U lc or U lc,sym is (relatively) open.
Proof. The first point follows from an observation of Bhattacharya (c.f. [Bha03] and also Lemma 4 in [Rül16]). The second point either follows from a direct calculation or by an application of Lemma 2.8 below.
The following lemma establishes a relation between rank-one connectedness and symmetrized rank-one connectedness. It in particular shows that in two dimensions all symmetric trace-free matrices are pairwise symmetrized rank-one connected.
Proof. We refer to [Rül16], Lemma 9 for a proof of this statement.
This lemma allows us to view symmetrized rank-one connectedness essentially as equivalent to rank-one connectedness.
2.2.2. Skew parts. We discuss some properties of the associated skew symmetric parts of rank-one connections, which occur between points in the interior of intconv(K). To this end, we introduce the following identification: Notation 2.9 (Skew symmetric matrices). As the two dimensional skew symmetric matrices are all of the form we use the mapping S →ω to identify Skew(2) with R. We define an ordering on Skew(2) by the corresponding ordering on R, i.e.
We begin by estimating the symmetric and skew-symmetric parts of a symmetrized rank-one connection: Lemma 2.10. Let a ∈ R 2 , n ∈ S 1 with a · n = 0. Then a n = |a|/2 = ω(a ⊗ n) .
As n, 1 |a| a forms an orthonormal basis, this shows that a n = |a|/2.
Using the previous result, we can control the size of the skew part which occurs in rank-one connections with K: Lemma 2.11. For all matrices N with e(N ) ∈ intconv(K) and with N being rankone connected with a matrix e (j) ∈ K it holds ω(N ) ≤ 10.
2.2.3. Geometry of the hexagonal-to-rhombic phase transformation. In this subsection, we discuss the specific matrix space geometry of the hexagonal-to-rhombic phase transformation. To this end we decompose each matrix of the form α β β −α into a component in v 1 = 1 0 0 −1 and a component in v 2 = 0 1 1 0 direction.
With this notation we make the following observations: Lemma 2.12. Let v 1 , v 2 ∈ R 2×2 sym be as above. Then, Furthermore, we have that .
Proof. Using the trigonometric identities an immediate computation shows the claim. the respective pairs (a ij , n ij ) are marked in the same color. We note that by symmetry it is also possible to pass from (a, n) to (−a, −n). After normalizing appropriately, symmetry further allows to exchange the roles of a, n. The figure on the right depicts the normals, which arise in the decomposition of differences e−e (i) with e (i) ∈ K and e ∈ C d (c.f. Lemma 2.15). The colors represent the well, with which e ∈ C d is connected: Red corresponds to the cone at e (1) , blue to the cone at e (2) and black to the cone at e (3) . As C d does not contain the full convex hull of K (c.f. Figure 3), only vectors, which lie within the colored zones arise as possible normals (in particular, there is a gap between these vectors and the vectors, which arise as decomposition of differences of the wells).
In other words, Lemma 2.12 allows us to identify all lines in matrix space (through the origin) by their rotation angle. In particular, this gives a simple description of the possible rank-one connections between the energy wells (c.f. also Figure 2 (a)). In our application to the hexagonal-to-rhombic phase transformation we have to take into account that non-trivial differences of the matrices e (1) , e (2) , e (3) lie on the sphere of radius √ 3 in matrix space (with respect to the spectral norm), which yields slightly different normalization factors for a: Lemma 2.13. Let K be as in (5). Then we have that  Proof. This is a consequence of Lemma 2.12 and the form of the matrices in (5).
d Figure 3. The colored domain represents the interior of the star C d from Lemma 2.15. As not the full cone of rank-one directions between the wells are possible (c.f. 22), the difference in the angle between different cones is bounded strictly away from zero and π.
Proof. This is a direct consequence of Lemma 2.12 and of the fact that the set K forms an equilateral triangle in strain space.
As an immediate consequence of Lemma 2.13 and Lemma 2.14, we infer the following result, which is graphically illustrated in Figure 2  Lemma 2.15. Let K be as in (5) and let with Q being a rotation by π 6 and with e (i1) , e (i2) ∈ K and e (i1) = e (i2) . Define α m1,m2 ∈ (0, π) as Then there exists a constant C = C(d) ∈ (0, π/4) such that Proof. Arguing as for (21) by using the definition of the set C d , we infer that the angles ϕ that occur in the representation from Lemma 2.12 for e−e (i) e−e (i) with e ∈ C d and e (i) ∈ K satisfy Here C(d) ∈ (0, π/3) is a constant, which depends only on d. The associated symmetrized rank-one connection is determined by ϕ as stated in Lemma 2.14. Applied to the situation in Lemma 2.15 this implies that a i1 , n i1 , a i2 , n i2 are expressed in terms of ϕ (as in Lemma 2.14). Since d > 0, the sectors parametrized by ϕ however do not overlap for i 1 = i 2 . As a consequence of this and of the options in (22), only the claimed angles α m1,m2 occur.

The Convex Integration Algorithm
In this section we present and analyze our convex integration algorithm (c.f. Algorithms 3.8, 3.11). Our discussion of this consists of four parts: First in Section 3.1 we introduce a replacement construction, in which a deformation gradient can be modified (c.f. Lemma 3.1-3.5). Here we follow Otto's Minneapolis lecture notes [Ott12] and refer to this construction as a version of Conti's construction (c.f. [Con08], [CT05] (Appendix) but also [Kir03]). Next in Section 3.2 we explain how the Conti construction can be exploited to formulate the convex integration algorithm (c.f. Algorithms 3.8, 3.11). Here we deviate from the more common qualitative algorithms by precisely prescribing error estimates in strain space, by specifying a covering construction and by controlling the skew part quantitatively. In Section 3.3 we analyze our algorithms and show that they are well-defined (c.f. Proposition 3.12). We further provide a control on the skew part of the resulting construction (c.f. Proposition 3.15). Finally, in Section 3.4 we use Algorithms 3.8, 3.11 to deduce the existence of solutions to the inclusion problem (6), c.f. Proposition 3.16.
We remark that our version of the convex integration scheme is based on particular properties of our set of strains: For the hexagonal-to-rhombic phase transition the laminar convex hull equals the convex hull (c.f. Lemma 2.7). Moreover, we can connect any matrix in intconv(K) with the wells K (c.f. Lemma 2.8). For a general inclusion problem this is no longer possible and hence more sophisticated arguments are necessary. In spite of the restricted applicability of the scheme, we have decided to focus on the hexagonal-to-rhombic phase transformation, as it yields one of the simplest instances of convex integration and illustrates the difficulties and ingredients, which have to be dealt with in proving higher Sobolev regularity in the simplest possible set-up. We begin by recalling Otto's variant of Conti's construction [Ott12], see also the video at [Ott]: Lemma 3.1 (Undeformed Conti construction). Let Ω = (−1, 1) 2 and let Then there exists u : R 2 → R 2 Lipschitz such that {x : ∇u(x) = M i } is a finite union of rectangles and triangles, Following Otto [Ott12], we generalize this construction slightly by allowing variable volume fractions: Then there exists u : R 2 → R 2 Lipschitz such that Multiplying u with 1 2λ and setting λ = 1 − λ , we obtain the matrices on page 56 of [Ott12]: and volumes (total volume 4) Proof of Lemma 3.2 and 3.1. Lemma 3.1 is a special case of Lemma 3.2 and (24) with λ = 1/2. We thus describe the general construction for arbitrary λ ∈ (0, 1) and then show that this reduces to the matrices from Lemma 3.1 in the case λ = 1 2 . In order to construct the function u, we prescribe the value of u at the points (0, ±λ), (±λ, 0) for λ ∈ (0, 1) to be determined and then consider linear interpolations. This then yields a piecewise affine Lipschitz map. It remains to verify that all matrices are divergence-free and are as claimed in the lemmata (c.f. also Figure  4).
We start with the ansatz given in Figure 5. The value of u at the points (±λ, 0), (0, ±λ) is chosen in such a way that linear interpolation in the triangles on the sides of the square in Figure 4 yields ∇u ∈ {M 0 , M 4 }, i.e.

Figure 5. A symmetric ansatz
By linear interpolation on the inner diamond we hence obtain that Choosing λ = 1 2 , we thus obtain ∇u = M 1 . It remains to check the value of ∇u on the triangles, which interpolate between the sides of the inner diamond and the corners of the outer square. By symmetry it suffices to consider the lower left triangle. Using again linear interpolation, there ∇u has to satisfy Setting λ = 1/2 this equals which is the matrix M 2 from Lemma 3.1.
Using this construction as a basic building block, the following lemma allows us to replace a general matrix M ∈ R 2×2 and to restrict the replacement matrices to an -neighborhood of a rank-one line passing through M . Figure 6. The horizontal axis corresponds to the upper right component of the matrix: Lemma 3.3 (Deformed Conti construction, page 57 of [Ott12]). Let M, M 0 , M 1 be given matrices such that Then, for every > 0 there exist matricesM 1 ,M 2 ,M 3 ,M 4 with a rectangle Ω ⊂ R 2 of aspect ratio δ = 20|a| and a Lipschitz map u : Furthermore, the level sets of ∇u are given by the union of at most 16 triangles.
Remark 3.4. We remark that both the side ratio δ as well as the error remain unchanged under rescalings of the form µu( x µ ) (as this leaves the gradient invariant). We now show how to apply Lemma 3.3 to the setting of symmetric matrices in our three-well problem (5): Lemma 3.5 (Application to the three-well-problem, pages 60 ff. of [Ott12]). Sup- and let e (i) with i ∈ {1, 2, 3} be such that

100
. Then, for every 0 < < 0 there exist a Lipschitz function u : R 2 → R 2 , a rectangular domain Ω (with ratio 1 : Proof. Let M and e (i) be given. Since e (1) , e (2) , e (3) are arranged in an equilateral triangle with side lengths √ 3 (with respect to the spectral norm) and as (27) holds, there existsẽ 1 ∈ intconv(K) such that Next let S := ω(M ) ∈ Skew(2) and letS ∈ Skew(2) to be determined. Then we obtain Since we are in two dimensions, any two symmetric, trace-free matrices are symmetrized rank-one connected (c.f. Lemma 2.8). Thus, there exist vectors a ∈ R 2 \ {0}, n ∈ S 1 such that Furthermore, as tr(e (i) ) = tr(ẽ 1 ), a and n are orthogonal. Choosing we thus obtain that the matrices are rank-one connected (with difference a ⊗ n or n ⊗ a, respectively) and We may hence apply the construction of Lemma 3.3 with M, M 0 , M 1 as defined above. Noting that e(M ) − e (i) = |a| 2 (c.f. Lemma 2.10), Lemma 3.3 implies the statement on the side ratio for Ω. Finally, we note that the -closeness of the matricesM 1 ,· · · ,M 4 also implies that their symmetric parts are -close.
Notation 3.6. In the preceding Lemma 3.5 the matricesM 0 , · · · ,M 4 obey the same (convexity) relations as the ones in Lemma 3.3, where for the matrices M 0 and M 1 we insert the ones from (29), c.f. Figures 6, 7. The error estimates in (25) thus • motivate us to refer to the matrixM 4 as stagnant (with respect to the replaced matrix M ). • The matricesM 1 ,M 2 ,M 3 will also be called pushed-out matrices (with the factors 4 3 and 16 15 respectively), since by construction and similarly for the other matrices. In order to emphasize the dependence on M , we also use the notatioñ Although the matricesM 0 , . . . ,M 4 also depend on the choice of e (i) , in the sequel we will often suppress this additional dependence for convenience as the reference well will be clear in most of our applications. We refer to the construction of Lemma 3.5 as the ( , δ) Conti construction with respect to M, e (i) . If some of the parameters of this are self-evident from the context, we also occasionally omit them in the sequel.
We emphasize that in our construction in Lemma 3.5, we have the choice between two different solutions, which differ in the sign of their skew symmetric component and thus in the choice of the corresponding rank-one connection (c.f. (28)). This freedom of choice is a central ingredient in the control over the skew symmetric part of the iterated constructions. We summarize this observation in the following corollary.
Corollary 3.7. Let M, e (i) , 0 , be as in Lemma 3.5. Then there exist two Lipschitz functions u + , u − : R 2 → R 2 such that on the set where e(∇u ± ) = e (i) Furthermore, up to an error of size the skew parts on the other level sets are given by Proof. From (29) we read off the skew symmetric parts of M 0 , M 1 . The skew symmetric part of M 2 := 1 5 M 0 + 4 5 M 1 is a consequence of that. The result then follows from Lemma 3.3.
3.2. The convex integration algorithm. In this subsection we formulate our convex integration algorithm. It consists of two parts, Algorithms 3.8 and 3.11. The first part (Algorithm 3.8) determines the symmetric part of the iterated deformation vector field, while the second part (Algorithm 3.11) deals with the choice of the "correct" skew component. After formulating the algorithms, we prove their well-definedness (i.e. show that it is indeed possible to iterate this construction as claimed).
In the whole section we assume that the domain Ω and the matrix M in (6) fit together in the sense that Ω = Q β [0, 1] 2 , where Q β is the rotation of the Conti construction from Lemma 3.3 for M (and the closest energy well e (i) ). These "special" domains will play the role of the essential building blocks of the construction of convex integration solutions in general Lipschitz domains (c.f. Section 6).
We define our convex integration scheme: Algorithm 3.8 (Quantitative convex integration algorithm, I). We consider the following construction: Step 0: State space and data.
(a) State space. Our state space is given by Here j ∈ N and u j : Ω → R 2 is a piecewise affine function. The sets is constant on each of the sets Ω j,k . It essentially keeps track of the well closest to e(∇u j | Ω j,k ) for each j, k. The functions are constant on each set Ω j,k and vanish in Ω \ Ω j . They correspond to the error and side ratio in the Conti construction, which is to be applied in Ω j,k . The functions j , δ j are coupled by the relation Hence, in the following (update) steps, we will mainly focus on j and assume that δ j is modified accordingly.
We explain how to update u j and j , δ j on Ω j,k . We seek to apply the construction of Lemma 3.5 with j,k : in a part of Ω j,k . To this end, we cover the domain Ω j,k by a union of finitely many (up to null sets) disjoint triangles and rectangles. The rectangles are chosen as translated and rescaled versions of the domains in the ( j,k , δ j,k ) Conti construction with respect to the matrices from (31). We denote these rectangles by R k j,l , l ∈ {1, . . . , K j,k }, for some K j,k ∈ N and require that they cover at least a fixed volume fraction v 0 > 0 of the overall volume of Ω j,k (which is always possible, c.f. Section 4 for our precise covering algorithm).
We define new setsΩ k j+1,l , l ∈ {1, . . . ,K j,k }: These are given by the trian- R k j,l and by the triangles which form the level sets of the deformed Conti rectangles R k l .
Further we set Ω k j+1,l :=Ω k j+1,l . Carrying this out for all k ∈ {1, . . . , J j } hence yields a collection of triangles In the sets R k j,l we apply the Conti construction with the matrices from (31). In this application we choose the skew part according to Algo- Step 2 (a), we define u j+1 |Ωk j+1,l as the function from the corresponding Conti construction. More precisely, in each of the rectangles R k j,l the matrix M j,k has been replaced by the matrices j,k . For each x ∈Ω k j+1,l withΩ k j+1,l as above, we define . For the definition of δ j+1 we recall its coupling with j+1 . We further set for ∇u j+1 |Ω j+1,k =M 0 (M j,k ). Here we choose an arbitrary possible minimizer if there is non-uniqueness.
While this algorithm prescribes the symmetric part of the iteration, we complement it with an algorithm, which defines the choice of the skew part. Here the main objectives are to keep the resulting skew parts uniformly bounded (which is necessary, if we seek to obtain bounded solutions to (6)) and simultaneously to ensure the choice of the "right" rank-one direction (c.f. Section 5, Lemma 5.2). Here the rank-one direction has to be chosen "correctly" in the sense that the successive Conti constructions are not rotated too much with respect to one another (which corresponds to the "parallel" case, c.f. Definition 3.10). In order to make this precise, we introduce two definitions: The first (Definition 3.9) allows us to introduce an "ordering" on the triangles in {Ω j,k } k∈{1,...,Jj } for different values of j ∈ N. With this at hand, we then define the notions of being parallel or rotated (c.f. Definition 3.10).
..,J j+l } is (part of ) a level set of ∇u j+l and is obtained from D by an l-fold application of the update step of Algorithm 3.8 (where we specify the covering to be the one, which is described in Section 4). The set of descendants of D of order l is denoted by D l (D). We define . We then writeD ∈ P l (D) and also use the notation P(D) for the set of all predecessors of D.
With this we define the parallel and the rotated cases: j,k be as in Algorithm 3.8. Let D ∈ {Ω j,k } k∈{1,...,Jj } for j ≥ 1. Let j 0 = 0 be the smallest index, for which P j0 (D) D = D (i.e. P j0 (D) was the last triangle in Algorithm 3.8, to which Step 2 (b) was applied instead of Step 2 (a)). Then, if for a.e. x ∈ D e (p) we say that in step j the triangle D is in the parallel case. If there is no possible confusion, we also just refer to D as in the parallel case. If for a.e. x ∈ D e (p) we say that in step j the triangle D is in the rotated case. If there is no possible confusion, we also just refer to D as in the rotated case.
Let us comment on this definition: Intuitively, its objective is to describe whether successive Conti constructions can be chosen as essentially parallel or whether they are necessarily substantially rotated with respect to each other (hence, these notions will also play a crucial role in Section 4, where we construct our precise covering). More precisely, let SP j be as in Algorithm 3.8 and let j, j 0 , D,D be as in Definition 3.10. Then, at the iteration step j 0 the triangleD was a subset of one of the Conti rectangles R k j−j0,l . Thus, u j−j0 is modified according to the Conti construction with respect to ∇u j−j0 |D, e Step 2(b) of Algorithm 3.8, we can hence ensure that the successive Conti constructions are essentially parallel, if (32) is satisfied. We remark that for this argument to hold and for it to yield new, significant information, it was necessary in Definition 3.10 to mod out the cases, in which Step 2(a) was active, i.e. P l (D) = {D}, as during these there are no changes. If (33) holds, then the directions of the successive Conti constructions are necessarily substantially rotated with respect to each other (c.f. Lemma 4.3 for the precise bounds). In this case we cannot substantially improve the situation to being more parallel by choosing the skew part appropriately in Corollary 3.7. Thus, in the sequel, we will exploit these instances as possibilities to control the size of the skew part and to use this, if necessary, to change the sign of the skew direction. The precise formulation of this is the content of Algorithm 3.11.
Algorithm 3.11 (Quantitative convex integration algorithm, II). Let Ω, u j : Ω → R 2 and SP j for j ≥ 1 be as in Algorithm 3.8. We further consider This function will be defined to be piecewise constant on Ω and to be constant on each triangle Ω j,k . It will define the skew part of ∇u j on Ω j,k , i.e.
In the initialization step of Algorithm 3.8 we choose ω 1 arbitrarily.
Step 2: Update. Let j ∈ N, j ≥ 1. Let ω j and Ω j,k be given. Suppose thatΩ k j+1,l withΩ k j+1,l ∈ D 1 (Ω j,k ) is constructed from Ω j,k by our covering argument (c.f. Step 2 in Algorithm 3.8). Then we define ω j+1 as follows: (a) IfΩ k j+1,l is not part of one of the Conti constructions in the covering, then we set ,l is part of one of the Conti constructions in the covering, then by Algorithm 3.8 we seek to apply the construction of Lemma 3.5 with scale j | Ω j,k and e (p) j | Ω j,k , ∇u j | Ω j,k . Thus, by Corollary 3.7 we have two possible choices for the skew part of ∇u j+1 . These are determined by their sign. To define the sign, let j 0 ∈ N be the smallest integer such that D := P j0 (Ω j,k ) = Ω j,k . We then choose the sign of the new skew direction ω j+1 |Ωk j+1,l (and hence determine the whole corresponding skew part) according to After having carried out the relabeling step, in which we pass fromΩ k j+1,l to Ω j+1,l , the function ω j+1 is constant on each of the triangles in Ω j+1,l . Together with Algorithm 3.8 this completes the construction of ∇u j+1 .
Let us comment on these algorithms: Due to the structure of the convex hulls (Lemma 2.7), our convex integration algorithm produces a (countably) piecewise affine solution (in contrast to the solutions obtained by means of an in-approximation scheme). This is reflected in the fact that the deformation u j is not further modified in the piecewise polygonal domains in Ω \ Ω j . The preceding algorithm differs from a non-quantitative version of a convex integration scheme in several aspects: • We consider finite coverings of Ω \ Ω j instead of directly covering the whole domain. • We prescribe the choice of j quantitatively.
• We prescribe the skew part quantitatively.
These points are central in our higher regularity argument: As we seek to prove higher regularity by means of the interpolation result from Theorem 2 or Corollary 2.1, we have to control the BV norm of the resulting deformation gradients. However, by a countably infinite (self-similar) covering of the whole domain, this is in general not possible (the total perimeter of the covering triangles is not bounded in general). Hence we only consider finite coverings, which produce a controlled (but growing) BV norm and simultaneously allows us to cover a sufficiently large volume fraction v 0 of our domain Ω j . That it is possible to satisfy these two competing aims is content of the covering results of the next sections (c.f. Propositions 4.16, 4.19). This finite covering of Ω j,k is the cause for the splitting of Step 2 into two parts. Part (a) deals with the triangles which are not covered by Conti constructions and are in this sense "errors" (in the sense that u j is not modified here), while part (b) deals with the part of the domain that is covered by Conti constructions, on which u j is modified. The specification of j is of key relevance as well. It distinguishes in a quantitative way whether a new rank-one connection is rotated or not with respect to the corresponding last rank-one connection. In our BV estimate this leads to different bounds (c.f. the perimeter estimates in Propositions 4.16, 4.19). In particular we cannot afford substantial rotations, as long as j 0 is very small, since this would yield superexponential growth for the BV norms, which cannot be compensated in our estimates (c.f. Figure 8 and the corresponding explanations for the intuition behind this). Due to the relation between the size of the scales δ j (which itself is directly coupled to the admissible error j ) and our regularity estimates, we in general seek to choose the value of j as large as possible without leaving intconv(K). By the intercept theorem, it is always possible to choose j to be "relatively large" in the push-out steps (c.f. Notation 3.6). However, for stagnant matrices, this is no longer possible.
Here we have to ensure a choice of j , which is summable in j ∈ N (in Algorithm 3.8 we choose it geometrically decaying), in order to avoid leaving intconv(K). These considerations lead to the case distinction in the definition of j+1 in Step 2 (b) of Algorithm 3.8. Finally, the quantitative prescription of the skew part is central to deduce the quantitative BV bound of Lemma 5.2, as we have to take care that, as long as we remain "parallel" in strain space (c.f. Definition 3.10), we approximately preserve the same skew direction. This is necessary to prevent the Conti constructions from being substantially rotated with respect to each other if j is very small and constitutes a crucial ingredient in the derivation of our perimeter and BV estimates in Sections 4 and 5 (c.f. Figure 8 for the intuition behind this). The normalization of the initial skew part is convenient (though not necessary).
3.3. Well-definedness of Algorithms 3.8, 3.11. We now proceed to prove that Algorithms 3.8 and 3.11 are well-defined. Here in particular, it is crucial to show that with our choice of the admissible error j , we do not leave intconv(K) in the iteration except to attain one of the energy wells in K (c.f. Proposition 3.12). Moreover, we seek to construct solutions to (6), which are Lipschitz regular. These points are the content of the following two Propositions 3.12, 3.15, which deal with the symmetric and anti-symmetric parts respectively. To show these we will rely on several auxiliary observations. 3.3.1. Symmetric part. We begin by discussing the symmetric part and by showing that in our construction it does not leave intconv(K), except to reach K.
Proposition 3.12 (Symmetric part). Let and let SP j and M be as in Algorithm 3.8. Then for every j, k ∈ N and every domain Ω j,k ∈ {Ω j,k } k∈{1,...,Jj } there holds In particular, for all j ≥ 1 it holds that ∇u j (x) ∈ intconv(K) for almost all x ∈ Ω j .
Proof. We prove the statement inductively. For j = 0, we note that this holds since j = 0 and ∇u 0 = M . Let thus ∇u j | Ω j,k =: M j,k be given. We only show that the result remains true for j + 1 in the regions, in which the Conti construction is applied, as in the other regions it holds by the induction hypothesis (as j+1 = j for these regions). LetM 0 (M j,k ), . . . ,M 4 (M j,k ) be the matrices, by which M j,k is replaced in the application of the Conti construction of Lemma 3.5. We consider first the pushed out matrices (see also Notation 3.6). If the edge of ∂ conv(K) closest tõ M l (M j,k ), l = 1, 2, 3, is different from the edge closest to M j,k , then by construction d(M l (M j,k )) ≥ 1/16. It thus remains to discuss the situation, in which this is not the case. In this situation the intercept theorem and the induction hypothesis, for l = 1, 2, 3, (for which j+1 |Ωk Here we used the definition of 0 (c.f. Step 0 (b)) and the induction hypothesis combined with the bound Finally, forM 4 (M j,k ) we estimate This concludes the proof.
3.3.2. Skew symmetric part. In order to deal with the skew part and to show its boundedness, we need several auxiliary results. These are targeted at controlling the maximal number of push-out steps in the parallel case (c.f. Lemma 3.14), where the notions "parallel" and "rotated" are used as in Definition 3.10. With the control of the maximal number of push-out steps at hand, we can then present a bound on the skew part of the gradients from Algorithms 3.8 and 3.11 (c.f. Proposition 3.15). Together with the boundedness of the symmetrized gradient this yields the uniform L ∞ bounds on ∇u j .
We begin by estimating the distance to the wells.
Lemma 3.13. Let SP j be the j-th step of the convex integration construction obtained in Proposition 3.16. Then for every level set Ω j,k it holds where d 0 , d K , 0 are as in Step 0(b) in Algorithm 3.8.
The statement of this lemma is very similar to the result of Proposition 3.12. However, instead of controlling the distance to the boundary, we here estimate the distance to the wells. This can be substantially larger than the distance to the boundary.
Proof. The proof follows along the same lines of the one of Proposition 3.12. We note that the statement is true for j = 0 (by the definition of 0 ) and proceed by induction. Let thus M j,k := ∇u j | Ω j,k be given. With slight abuse of notation we set j := j | Ω j,k . It suffices to show that the values of ∇u j+1 , which were obtained from M j,k by an application of the Conti construction, still satisfy the desired estimates (in the domains, in which u j is unchanged the estimate holds by the induction assumption). The application of the Conti construction yields matrices M 0 (M j,k ), · · · ,M 4 (M j,k ). As e(M 0 (M j,k )) ∈ K, we only consider the other matrices. We consider the matricesM 1 (M j,k ), · · · ,M 3 (M j,k ), which are constructed by "pushing-out" (c.f. Notation 3.6). Without loss of generality (c.f. the argument in Proposition 3.12), we only discuss the case that the closest well for e(M i (M j,k )) is the same as for e(M j,k ). For i ∈ {1, 2, 3} we have Here we used the induction assumption for M j,k as well as the estimate for d(M j,k ) from Proposition 3.12 and the definition of 0 .
Here we used the definition of j+1 := j /2 on the subset of the Conti construction, on whichM 4 (M j,k ) is attained.
Using Lemma 3.13 and recalling Definitions 3.9, 3.10, we bound the maximal number of possible push-out steps in the parallel situation: Lemma 3.14. Let SP j and ω j be as in Algorithms 3.8 and 3.11. Assume that D ∈ {Ω j0+n,k } k∈{1,...,Jj 0 +n } and suppose that the construction of D fromD ∈ P n (D) involves k with k ∈ N ∪ {0} push-out steps (c.f. Notation 3.6). Further assume that for a.e. x ∈ D and for all r ∈ {1, . . . , n} e (p) Then there exists a number Proof. The proof relies on the definition of 0 and the control on the distance to the wells, which was obtained in Lemma 3.13. Indeed, let Ω j0+l,m ∈ {Ω j0+l,m }m ∈{1,...,J j 0 +l } with Ω j0+l,m ⊂D be arbitrary but fixed. Without loss of generality, we assume that in all the iteration steps j 0 , . . . , j 0 + l Step 2(b) occurs on our respective domain (as there is no change, if Step 2(a) occurs, and as we are only interested in the maximal number of steps, in which a specific change, i.e. a push-out, occurs). Let M j := ∇u j | Ω j 0 +l,m be given. Suppose that a matrix M j0+n+1 is obtained from M j0+n for some n ∈ {1, . . . , l} by push-out and that M j0+n is obtained from M j0 by stagnating n-times. Then, dist(e(M j0+n+1 ), e where d K is defined as in Step 0 (b) in Algorithm 3.8. Therefore, by Step 2 of Algorithm 3.8 (i.e. the update for e (p) j ) after at most push-out steps, we are no longer in the parallel case. This yields the desired upper bound.
Relying on the previous lemma, we obtain a uniform bound on the skew part: Proposition 3.15 (Skew symmetric part). Let SP j and ω j be as in Algorithms 3.8 and 3.11. Suppose that N 0 > 0 is the number from Lemma 3.14. Definē C := max{100, 20(N 0 + 1)(1 + 0 )}. Then, and Proof. We prove the claims inductively and note that ω 0 = 0 satisfies them. We first discuss (35) and show that it remains true for ω j with j ∈ N. To this end, let l ∈ N and D ⊂ {Ω j+l,k } k∈{1,...,J j+l } . For abbreviation we set M j := ∇u j | D , ω j := ω j | D (and recall that ω j | D = ω(∇u j | D )) and first assume thatω j ≤ 0 (see Notation 2.9). We begin by making the following additional assumption: Assumption 1. We suppose that the skew matrixω j+l is derived fromω j by an l-fold application of Algorithms 3.8 and 3.11, where in the Conti construction of Corollary 3.7 we always choose the positive skew direction.
We point out that this assumption can occur both in the parallel and in the rotated case, but ensures that the skew direction was not changed in this process. In other words, Assumption 1 implies that the sign of the skew direction, which is chosen in Corollary 3.7 remains fixed. We hence refer to this situation as the "fixed sign case". We further introduce the auxiliary functions N l,1 , N l,2 , N l : Ω → N ∪ {0}, with N l := N l,1 + N l,2 . Here for given l ∈ N andω j+l , we define N 1,l as the number of 4/3 push-out steps in the process of obtainingω j+l fromω j , and N 2,l as the number of 16/15 push-out steps. By Lemma 3.14 we know that 0 ≤ N l ≤ N 0 .
Step 1: Upper bound in the fixed sign case. We first deal with the upper bound forω j+l . To this end we note that We iterate this estimate: Here we used the estimate dist(e(∇u j )| D , K) ≤ 3, the fact that in each push-out step an error 0 is possible, while in each stagnant step the error is decreased by a factor two. Recalling the definition ofC and the fact that N l ≤ N 0 hence implies ω j+l ≤C/2 +ω j .
Using thatω j ≤ 0, therefore allows us to conclude that ω j+l ≤C/2.
Step 2: Lower bound in the fixed sign case. Still working under the assumptions from above, we now bound the negative part ofω j+l . Here we show that We first consider the push-out steps. LetM 1 (M j+l−1 ),M 2 (M j+l−1 ),M 3 (M j+l−1 ) be the push-out matrices in the corresponding Conti construction of Algorithms 3.8, 3.11. Their skew parts are contained in 0 neighborhoods of ω(M j+l−1 ) + 1 3Ŝ , ω(M j+l−1 ) + 1 15Ŝ .
By Lemma 2.10, Lemma 2.11 and Proposition 3.12, we obtain that Here d : R 2×2 → R denotes the function from Proposition 3.12, and a ⊗ n is the rank-one connection, which appears in the Conti construction. In the last estimate we have used the estimate from Proposition 3.12. Hence the skew parts ofM 1 (M j+l−1 ),M 2 (M j+l−1 ),M 3 (M j+l−1 ) are respectively bounded by which shows the claimed estimate (35) with j = 0 . ForM 4 (M j+l−1 ) we have that which also proves the desired result. This concludes the proof of (35) in the fixed sign case.
Step 3: Sign change. Let j + l + 1 be the first index, in which the sign of the difference of the skew parts changes according to Algorithm 3.11. By Assumption 1 and by the definition of our Algorithms 3.8, 3.11, this can only be the case if ω j+l ≥ 0. The definition ofC ensures that theM 4 (M j+l+1 ) obeys the upper bound ω(M 4 (M j+l+1 )) ≤C 2 + 0 ≤C.
For the pushed out parts,M 1 (M j+l+1 ),M 2 (M j+l+1 ),M 3 (M j+l+1 ), we argue similarly as we did in Step 1, but now with a change of signs: By the intercept theorem, the resulting skew parts become strictly smaller than the one ofω j+l (potentially they even become negative). This then improves the upper bound (37). For the lower bound we argue as in Step 1 but with reversed sign in Assumption 1. This concludes the proof of (35).
Step 4: Proof of (36). In order to obtain the estimate (36), we notice that the skew parts associated with values of e(∇u j ) ∈ K may on the one hand be strictly larger than the bound given in (35). But on the other hand, they are derived as añ M 0 matrix in one of the Conti constructions, in which matrices satisfying (35) are modified. This implies that at most a gain of 5 in the modulus of the corresponding skew part is possible, which yields the bound (36). As these domains are not further modified in the convex integration algorithm this bound cannot deteriorate in the course of the application of Algorithms 3.8 and 3.11.
3.4. Existence of convex integration solutions. Finally, in this last subsection, we show that Algorithms 3.8, 3.11 can be used to deduce the existence of solutions to our problem (6). Proof. We apply Algorithm 3.8 withM := M − ω(M ). By the results of Propositions 3.12 and 3.15 this algorithm is well-defined and can be iterated with j → ∞. This yields a sequence of functions u j : R 2 → R 2 with bounded gradient (with ∇u j L ∞ (R 2 ) depending on d K , c.f. Lemma 3.14). We prove the convergence of this sequence and show that the limiting function u 0 solves (5) with boundary datā M . We note that for k ≥ j Moreover, ∇u k =M a.e. in R 2 \ Ω.
By construction ∇u j is bounded, hence ∇u j ∇u 0 in the L ∞ loc (R 2 ) weak- * and the L 2 loc (R 2 ) weak topologies. By Poincaré's inequality u j → u 0 in L 2 loc (R 2 ). We observe that Step 2 in Algorithm 3.8 decreases the total volume of the Ω j , i.e. of the part of Ω, on which e(∇u j ) does not yet attain one of the wells: Combined with (38) and the L ∞ bound, this implies the desired convergence ∇u j → ∇u 0 with respect to the L 2 loc (R 2 ) topology, where ∇u 0 ∈ L ∞ (R 2 ) is a solution to the problem (5) with boundary dataM . Defining u(x) := u 0 (x) + ω(M ) hence concludes the proof of Proposition 3.16.
In Sections 4 and 5 we present a more refined analysis of this construction algorithm. In particular, we give an explicit quantitative construction for the covering procedure from Step 2 in Algorithm 3.8.

Covering Constructions
In the following section we present the details of the coverings, which we use in the Algorithms 3.8, 3.11. Here we pursue two (partially) competing objectives: Given a triangle D, • we seek to cover an as large as possible volume fraction of it, but at least a given fixed volume fraction, v 0 > 0. • We have to control the perimeters of the triangles in the resulting new covering. In the context of these considerations, it turns out that the parallel and the rotated cases (c.f. Definition 3.10) differ quantitatively and hence have to be discussed separately. This can be understood when considering possible coverings of rectangles by parallel or rotated rectangles.
We illustrate this in two extreme situations (c.f. Figure 8): Given a rectangle R 1,δ0 with sides of length 1 and δ 0 , we seek to cover it with rectangles, which have a fixed side ratio r and whose long sides are either parallel or orthogonal to the long side of the original rectangle R 1,δ0 . In order to illustrate the differences between these situations, we for instance assume that r = δ 0 /2. In the situation, in which the original rectangle R 1,δ0 is covered by rectangles, whose long side is parallel to the long side of R 1,δ0 , the covering can be achieved by splitting R 1,δ0 1 δ 1 δ 0 δ j δ 0 Figure 8. Covering a rectangle R 1,δ0 of side lengths 1 and δ 0 by (a) a parallel rectangle of half its aspect ratio, (b) an orthogonal rectangle of aspect ratio r = δ j (which could for instance be r = δ 0 /2). along its central line as illustrated in Figure 8 (a). Thus, the resulting perimeter (we view it as a measure of the BV energy of the characteristic functions in the Conti covering), which is necessary to cover the volume of R 1,δ0 is bounded by twice the perimeter of R 1,δ0 . If the long sides of the covering rectangles of ratio δ 0 /2 are however orthogonal to the long side of R 1,δ0 , the covering of R 1,δ0 can only be achieved by 2δ −2 0 small rectangles of side lengths δ 0 and δ 2 0 /2 (c.f. Figure 8 (b)). The necessary perimeter for this covering is thus proportional to δ −1 0 Per(R 1,δ0 ). For a small value of δ 0 this makes a substantial difference and accounts for the losses in the estimates for the rotated situation. The difference of the parallel and the rotated situation become even more apparent, if we consider a sequence of coverings: Here we start with the rectangle R 1,δ0 and first consider an iterative covering of it by parallel rectangles, which in the j-th iteration step are of side ratio δ j := 2 −j δ 0 (and such that the long side is parallel to the long side of R 1,δ0 ). The desired covering of R 1,δ0 in the iteration step k can be achieved by splitting the rectangles from the covering at the iteration step k − 1 along their central lines. In each iteration step the overall perimeter increases at most by a factor two, so that after j iteration steps the overall perimeter can be estimated by 2 j Per(R 1,δ0 ).
If in comparison, we consider the case, in which the covering rectangles are rotated in every step by π/2 with respect to the preceding rectangles and again choose a ratio δ j := 2 −j δ 0 in the j-th iteration step, we inductively obtain a bound of the form for the overall perimeter after the j-th step. In contrast to the parallel situation this has superexponential behavior in j. If we consider the π/2 rotated situation with fixed ratio δ j = δ 0 , this bound improves to an exponential bound of the form δ −j 0 Per(R 1,δ0 ). Hence, the estimates in the rotated situation are substantially worse than the ones in the parallel situation. In order to avoid superexponential behavior, we have to take care that the rotated case can only occur, if the value of δ j is controlled from below. These heuristics a posteriori justify our careful choice of j and ω j in Algorithms 3.8, 3.11.
Although the level sets of the Conti construction consist of triangles and hence our coverings {Ω j,k } k∈{1,...,Jj } will be coverings of triangles by triangles (instead of the previously described rectangular coverings), the heuristics from above still persist.
Motivated by these heuristic considerations, in the sequel we seek to provide covering results and associated BV bounds, which can be applied in Algorithms 3.8, 3.11. We organize the discussion of this as follows: In Section 4.1, we introduce some of the fundamental objects (c.f. Definitions 4.4, 4.7) and formulate the main covering result (Proposition 4.9). Here we consider a similar distinction into a parallel and a rotated situation as described in the above heuristics (c.f. Definition 4.7). With the class of triangles from Definition 4.4 at hand we distinguish several different cases and discuss different covering scenarios. The respective coverings are tailored to the specific situation and are made such that we do not leave our class of triangles during the iteration. Their discussion is the content of Sections 4.2-4.5. Finally, the various different cases are combined in Section 4.6 to provide the proof of Proposition 4.9.
4.1. Preliminaries. In this section we introduce the central objects of our covering (c.f. Definition 4.4, 4.7) and state our main covering result (Proposition 4.9).
As a preparation for the main part of this section, we begin by discussing auxiliary results on matrix space geometry. We first estimate the angle formed in strain space between two matrices: Remark 4.2. Applied to a triangle D ∈ {Ω j,k } k∈{1,...,Jj } in the parallel case (c.f. Definition 3.10), Lemma 4.1 implies that the rotation angle α 1 , with which the consecutive Conti constructions are rotated with respect to each other (and which is defined as in Lemma 4.1), is bounded by Here j , δ j are the functions from the convex integration Algorithm 3.8.
Proof of Lemma 4.1. As sketched in Figure 9, we may estimate where j is the error in matrix space. Here the first estimate follows by a Taylor approximation and by noting that |α 1 | is small, so that in particular |α 1 | ≤ π 6 (in which range the tangent is invertible and for which the Taylor expansion is valid).
Next we observe the following bounds on the rotation angles: α ∼ d Figure 9. The angle between rank-one connections of nearby matrices is small. In particular, the associated rectangle constructions from Lemma 3.3 are close to being parallel.

Lemma 4.3 (Angles).
Let D ∈ {Ω j,k } k∈{1,...,Jj } for j ≥ 1. LetD ∈ D 1 (D). Assume that the triangleD is in the rotated case (c.f. Definition 3.10). Let α 1 denote the angle between the long sides of the current and the following Conti constructions. Then, we have that Proof. This is an immediate consequence of Lemma 2.15.
With these auxiliary results at hand, we proceed to the discussion of our central covering objects. In order to define our set of covering triangles, {Ω j,k } k∈{1,...,Jj } , we consider a subclass of triangles with, for our purposes, suitable properties. To this end, we can not ensure that all domains appearing in our covering argument are right angle triangles (due to the presence of the green triangles in the Conti construction in Figure 10), for which one of the other angles is approximately of size δ j . However, the following definition provides a family of sets with similar properties. This will allow us to formulate a precise, iterative covering result.
Definition 4.4. Let δ 0 be as in Algorithm 3.8. The triangle D is said to be δ-good with respect to a reference direction n ∈ S 1 , for δ ∈ (0, δ 0 ], if (1) One angle, α, satisfies α ∈ δ[ 1 10 , 1000], (2) The other two angles are contained in π 2 + 2δ[−1000, 1000], (3) One of the long sides encloses an angle in δ[−1000, 1000] with n. We refer to the long side of the triangle, which satisfies the requirement of 3.) as the direction of D. We also say that the triangle D is oriented parallel to n or that D is aligned to n. If a triangle D satisfies 1.) and 2.) but not necessarily 3.) we call it δ-good (which allows for a possible change of orientation).
Remark 4.5. We note that if α is small, both long sides could satisfy condition 3.) at the same time. In this case both directions are valid as directions of D.
In our construction one prominent reference direction is obtained from Conti's construction, as detailed in the following definition.
Definition 4.6. Let D be a level set of ∇u j and let e = e (p) j | D be the reference well. Let further n ∈ S 1 be the direction of the long side of the Conti rectangle from Step 2 in Algorithm 3.11. We say that n is the direction of the relevant Conti construction (at step j). A δ-good triangle is parallel to Conti's construction if one of its long sides is parallel to n.
In the sequel, we will give a precise covering result, which shows that in the j-th step of our convex integration Algorithms 3.8 and 3.11, we may assume that only very specific triangles are present in the collection {Ω j,k } k∈{1,...,Jj } as (parts of) level sets of ∇u j . To this purpose we define the following classes of triangles: Definition 4.7. Let SP j be as in the Algorithms 3.8, 3.11. Then, a triangle D j ∈ {Ω j,k } k∈{1,...,Jj } is in the case: Remark 4.8. The cases above, as stated, are not distinct since we allow for a factor in our definition of being δ-good (c.f. Definition 4.4). For instance, there might be triangles which are in both case (P 1) and (P 2). However, in such situations also the constructions and perimeter estimates are comparable. In situations, in which the estimates would differ significantly and where δ j−1 δ j = δ 0 , the above definitions yield distinct cases.
Let us comment on this classification: The basic distinction criterion separating the triangles into the different cases is given by checking whether the corresponding triangles are roughly aligned (as in the cases (P1), (P2)) or whether they are substantially rotated (as in the cases (R1), (R2)) with respect to the direction of the relevant Conti construction (the case (R3) is a special "error situation", which does not entirely fit into this heuristic consideration). Roughly speaking, this determines whether we are in a situation analogous to the first or to the second picture in Figure 8. This distinction is necessary, as else a control of the arising surface energy is not possible in a, for our purposes, sufficiently strong form. This distinction (essentially) coincides with our definition of the parallel and the rotated cases (c.f. Definition 3.10): If a triangle D j ∈ {Ω j,k } k∈{1,...,Jj } is in the parallel case in step j, then the directions of the Conti construction, which gave rise to D j , and of the relevant Conti construction at step j (i.e. the construction, by which D j is (in part) covered) are essentially parallel (c.f. Lemma 4.1 and Remark 4.2). Letting j 0 ∈ N denote the index from Definition 3.10 and assuming that j 0 = 1, there are three possible scenarios for the relation of δ j | Dj and δ j−1 | Dj : (i) δ j | Dj = δ j−1 | Dj /2. In this case ∇u j | Dj was produced as the stagnant matrix (c.f. Notation 3.6) in the iteration step j − 1. In this case, our covering construction will ensure that D j is in case (P1) (not exclusively, c.f. Remark 4.8, but as one option). (ii) δ j | Dj = δ 0 but δ j−1 | Dj = δ 0 . This can for instance occur in a parallel push-out step. In this case, our covering construction will ensure that D j is in case (P2) (not exclusively (depending on the value of δ j−1 ), c.f. Remark 4.8, but as one option). (iii) δ j | Dj = δ 0 = δ j−1 | Dj . This case can for instance occur in two successive push-out steps. In this case, our covering ensures that D j is in the case (P1). If a triangle D j ∈ {Ω j,k } k∈{1,...,Jj } is in the rotated case in step j, then the direction of the Conti construction, which gave rise to D j , and the direction of the relevant Conti construction at step j are necessarily substantially rotated with respect to each other (c.f. Lemma 4.3). Again assuming that the index j 0 = 1 (where j 0 denotes the index from Definition 3.10), we now distinguish two cases for the relation between δ j | Dj and δ j−1 | Dj : Here we first note that necessarily (by definition of the rotated case, as occurring only after a push-out step) we have δ j | Dj = δ 0 . Then there are two options for δ j−1 | Dj : (i) δ j−1 | Dj = δ 0 . This case can for instance occur in the situation of two successive push-out steps. In this case our covering ensures that D j can be taken to be in the case (R1). (ii) δ j−1 | Dj = δ 0 . This case can for instance occur in the case, in which ∇u j−1 | Dj is produced in a stagnant and ∇u j | Dj in a push-out step. In this situation our covering ensures that D j can be taken to be in the case (R2). The case (R3) only occurs as an error case as a consequence of our specific covering procedure for the triangles of the types (R1) and (R2). We relate the different cases to the heuristics given at the beginning of Section 4 (c.f. Figure 8). We view the cases (P1) and (R1) as the "model cases" without and with substantial rotation and corresponding to the parallel and orthogonal (triangular) situation depicted in Figure 8. In both cases (P1) and (R1) the aspect ratio of the given triangle D j is roughly of order δ j | Dj (i.e. the quotient of its shortest and of its longest sides are roughly of that order) and we seek to cover it with a Conti construction of comparable ratio δ j | Dj . The cases (P2) and (R2) are situations, in which the underlying triangle D j is roughly of side ratio δ j−1 | Dj (i.e. the quotient of its shortest and of its longest sides are roughly of that order), where we however seek to cover the triangle with Conti constructions with ratio δ 0 . This mismatch is a consequences of our construction of the function δ j | Dj in Algorithm 3.8: Here we prescribe that the matrices, which are pushed out (c.f. Notation 3.6), are allowed to have an error tolerance of δ 0 . In particular, it may occur that δ j−1 | Dj δ j | Dj = δ 0 , which is the situation described in (P2), (R2) either without or with substantial rotation. The case (R3) is a consequence of how we deal with "remainders" in our covering constructions for the cases (R1), (R2).
Our main result of the present section states that it is possible to find a covering of the level sets, which respects Algorithms 3.8, 3.11, such that only the specific triangles from Definition 4.7 occur. Moreover, we provide bounds for the remaining uncovered "bad" volume and the resulting perimeters.  Here v 0 ∈ (0, 1) is a small constant, which is independent of 0 , d 0 and d K .
In the remainder of this section we seek to prove this result and to construct the associated covering. To this end, in Section 4.2 we first explain that the "natural covering" of the Conti construction, which is achieved by splitting it into its level sets, satisfies the requirements of Proposition 4.9. In particular, this implies that the covering, which is obtained in Step 1 of Algorithm 3.8, satisfies the properties of Proposition 4.9 (the resulting triangles are of the types (P1), (P2) or (R1), (R2)). Hence, in the remaining part of the section, it suffices to prove that given a triangle of the type (P1)-(R3), we can construct a covering for it, which obeys the claims of Proposition 4.9. To this end, in Section 4.3, we first describe a general construction, on which we heavily rely in the sequel. With this construction at hand, in Section 4.4 and its subsections we then deal with the cases (P1), (P2), in which there is no substantial rotation involved. Subsequently, we discuss the cases (R1)-(R3) with non-negligible rotations in Section 4.5. Finally, in Section 4.6 we provide the proof of Proposition 4.9.
The generalization to more generic domains is detailed in Section 6. 4.2. Covering the Conti construction by triangles. We begin with our covering construction by explaining that a Conti construction of ratio δ j | Dj can be divided into a finite number of triangles, which are all of the types (P1), (P2) and (R1)-(R3).
Lemma 4.10. Let SP j be as in Algorithm 3.8 and let D j ∈ {Ω j,k } k∈{1,...,Jj } . Suppose that R ⊂ D j is a Conti rectangle of ratio δ j | Dj . Let M 1 , . . . , M 4 denote the gradients occurring in the Conti construction with the same convention as in Notation 3.6. Then all level sets in R, on which M 4 is attained, can be decomposed into (at most two) triangles, which are of the type (P1). The level sets with M 1 , M 2 , M 3 can be decomposed into triangles of the type (P1), (P2) or (R1), (R2).
Proof. We recall that (after a suitable splitting into in total 16 triangles as depicted in Figure 10) all except for four triangles in the undeformed Conti construction (c.f. Lemma 3.2) are axis-parallel and have aspect ratio approximately 1 : δ j | Dj (with a factor depending on λ; for λ = 1 4 a factor in the interval (1/4, 4) is more than sufficient). After rescaling the x 2 -axis by δ j | Dj (as in Lemma 3.3), these aspect ratios are then comparable to 1 : It remains to discuss the remaining triangles contained in F (green in Figure  10). These are again δ j | Dj -good by a similar estimate on the aspect ratios, and by an estimate on the angle of rotation with respect to the x 1 -axis. As by our convex integration Algorithm 3.8, Step 2 (b), δ j+1 | F = δ 0 , they are in general of the type (P2) or (R2), if δ j | Dj = δ 0 , but could also be of the type (P1) or (R1), if δ j | Dj = δ 0 .
Remark 4.11. We observe that Lemma 4.10 in particular implies that the triangles, which are obtained in Step 1 in Algorithm 3.8, all satisfy the claim of Proposition 4.9. Hence, in the sequel it suffices to provide a covering algorithm, which preserves this property.

4.3.
A basic building block. We begin our iterative covering statements by presenting a general building block, which we will frequently use in the sequel. Given a triangle D we seek to reduce the discussion to that of a rectangle R 2 , whose long side is aligned with the direction of D and which is of similar volume as the original triangle. Only in the covering of this rectangle will the situations (P1), (P2) and (R1), (R2) differ. For the case (R3) we argue differently.
Proposition 4.12. Let D be a δ-good triangle with 0 < δ ≤ δ 0 . Let furtherR 2 be a rectangle of aspect ratio r ∈ [ 1 10 , 1000]δ, whose long side is aligned to the direction of D. Then there exists a rescaled and translated copy R 2 ⊂ D of the rectangleR 2 (of aspect ratio 1 : r), for which three of its corners lie on ∂D and such that: (i) |R 2 | ≥ 10 −6 |D|.
(ii) One corner divides a side of the triangle in the ratio 2 3 : 1 3 . (iii) The set D \ R 2 consists of at most 100 δ-good triangles, which are aligned with the direction of D.
3 tan(α) Figure 12. As the angle at the bottom right is very close to π 2 , the box is well approximated by a rectangle of side-lengths 1 3 : 2 3 tan(α). Partitioning the box into N slices of the same height, an estimate on the tan of the opening angles γ j of the corresponding rectangles shows that tan(γ j ) ≈ 2 N tan(α). Choosing N ∈ {1, 2, 3} appropriately, we thus obtain δ-good triangles.
Proof. Let D be a given δ-good triangle and let α denote the corresponding angle from Definition 3.9 (1). Without loss of generality we may assume that the triangle D is aligned with the x 1 -axis, that the tip of the triangle lies at the origin and that (after rescaling) the x 1 -axis-parallel side is given by the interval [0, 1] × {0} (c.f. Figure 11).
Let P 1 = ( 2 3 , 0) and let g be the line of slope −r through P 1 . Then g intersects ∂D in exactly one other point P 2 . Being aligned along the x 1 -axis, the rectangle R 2 is then uniquely determined by requiring that P 1 and P 2 are two of its corners. By construction it has aspect ratio 1 : r. In order to infer the bound on the volume, we compute the coordinates of P 2 = ( 2r tan(α)+r , 2r tan(α) 3(tan(α)+r) ). Hence the volume of R 2 is given by .
Since the volume of D is comparable to tan(α) 2 , this results in a volume fraction of approximately min(r, tan(α)) max(r, tan(α)) .
Using the fact that r ≥ δ 10 and α, arctan(r) ≤ 1000δ, we infer the desired estimate on the volume fraction. In particular, we note that it is independent of δ.
Adding a vertical line through P 1 and a horizontal line through P 1 +(0, 2 3 tan(α)) ∈ ∂D, we obtain an axis parallel triangle of opening angle α to the left of R 2 , another axis parallel triangle of opening angle α above R 2 , a four-sided box B to the right of R 2 and triangle self-similar to D above the box (c.f. Figure 11).
As D is δ-good, so are the above mentioned three triangles and it hence remains to discuss the box B on the right of R 2 (c.f. Figure 12). By construction the bottom side of B is axis parallel and of length 1 3 and the left side is also axis-parallel and of length 2 3 tan(α). Furthermore, since D is δ-good, the angle on the bottom right of B is given by π 2 − γ for some γ ∈ δ [−2000, 2000] and in particular |γ| ≤ 1 10 . Hence, the length of the axis-parallel top side differs from 1 3 by | tan(γ) 2 3 tan(α)| ≤ 1 10 . Introducing further horizontal lines, we may partition B into N boxes with three axis-parallel sides of height 2 3N tan(α) and length close to 1 3 (c.f. Figure 12). Bisecting along the diagonals, we hence obtain opening angles γ j with Choosing N ∈ {1, 2, 3} appropriately and noting that the remaining angle is either Figure 13. Constructing an axis parallel rectangle starting from R. We begin with the inner rectangle R =R (rotated with an angle β with respect to the x 1 -axis) and successively add the eight outer rectangles to obtain the boxR 2 from Lemma 4.14. The explicit construction of the four white boxes, which together with R form the inner "cross" are described in detail in Figures 15 and 16. All the outer rectangles are of aspect ratio approximately 1 : δ j and can hence be decomposed into δ j -good triangles. After a translation, rescaling and reflection with respect to the x 2 -axis, we may assume that β ≥ 0 and that the corners ofR are given by P 1 = (0, 0), P 2 = (1, tan(β)), P 3 = (1, tan(β)) + r 0 1 + tan(β) 2 (− sin(β), cos(β)) = (1, tan(β)) + r 0 (− tan(β), 1), (1, tan(β)) r0(− tan(β), 1) Figure 15. Adding δ-good triangles on the left and right, we can achieve axis-parallel boundaries. β γ2 γ2 + β Figure 16. Adding δ j -good triangles on the top and bottom, we can achieve axis-parallel boundaries. Since β might include a large or small factor in the definition, we may either allow and opening angle γ 2 + β or introduce a horizontal line to obtain two triangles with opening angle γ 2 and β.
In the following, we add quadrilaterals with three axis-parallel sides with aspect ratios r 0 and r 2 = tan(γ 2 ) for to be chosen later.
We begin by adding quadrilaterals on the left and right by inserting the following four points (c.f. Figure 14) (1, 0), (1, 0), We, in particular, note that the lines Q 1 Q 2 and Q 3 Q 4 are parallel to the x 2 -axis. Furthermore, the triangles Q 2 P 1 P 4 and P 4 Q 1 Q 2 have opening angles arctan(r 0 ), are parallel to the x 1 -axis and are either right-angled or have an angle π 2 −β. Hence all of these triangles are δ-good with direction e 1 . Similar observations hold for the triangles, which are constructed from P 2 , P 3 , Q 2 , Q 3 .
Here, the aspect ratio r 2 is chosen flexibly to account for the facts that our construction is horizontally of a length, which is slightly larger than 3, and that the rotated rectangle has height r 0 + tan(β) ≥ r 0 . By symmetry, we may restrict ourselves to discussing the rectangle P 1 Q 5 Q 6 P 2 . The axis-parallel right angled triangle P 1 Q 5 Q 6 has opening angle γ 2 (as defined in (39)) and is thus δ-good. For the remaining triangle P 1 Q 6 P 2 , we distinguish two cases: • If β ∈ δ[ 1 10 , 1000], we additionally introduce the point Q 9 = (1, 0) and note that P 1 Q 9 P 2 is δ-good and axis-parallel, as are P 1 Q 9 Q 6 and Q 6 Q 5 P 1 (which both have an opening angle γ 2 ).
• If 0 ≤ β ≤ δ 1 10 , we note that by our restriction on γ 2 , , 300], which ensures that P 1 Q 6 P 2 is δ-good and parallel to the long side n ofR.
Finally, we complete our thus far roughly cross-shaped construction to the desired axis-parallel rectangle R 2 by adding four rectangles as in Figure 13 (the green rectangles there). These rectangles have side lengths (1 + r 0 tan(β)) : r 2 , 1 : (r 2 + tan(β)).
For the second rectangle with ratio as in (41), we again distinguish two cases: • If β ∈ δ[ 1 10 , 1000], we divide the rectangle by a horizontal line through P 1 , which yields two rectangles of lengths 1 : tan(β), 1 : r 2 , which are δ-good.
Again this is satisfied due to our restrictions on γ 2 , r 0 and β.
We thus obtain a large family of admissible values r 2 . We note that the aspect ratio of R 2 is comparable to 2r 2 + tan(β) + r 0 , as is the area of R 2 . Hence, as a particular choice, we may take r 2 comparable to r0+tan(β) 2 (within a factor 3 to ensure that γ 2 satisfies the above restriction). Then the aspect ratio is comparable to 1 : and the volume ratio is comparable to |R| |R 2 | ≥ r 0 3 · 2(r 0 + tan(β)) = 1 6 r 0 r 0 + tan(β) .
Remark 4.15. One should think of the rectangle R in Lemma 4.14 as a Conti construction, which we seek to fit into a triangle D as in Proposition 4.12. In doing so, we however have to be careful, since in the cases (P1), (P2) we have to avoid creating new triangles, which are substantially rotated with respect to the original one (c.f. Figure 17 and the explanations at the beginning of the next section). The box construction of Lemma 4.14 ensures this.

4.4.
Covering in the cases (P1), (P2). In this section we explain how, given a triangle D j ∈ {Ω j,k } k∈{1,...,Jj } , which is of type (P1) or (P2), we can cover it by a combination of the relevant Conti constructions and some remaining triangles, which are again of the types (P1), (P2) and (R1), (R2). Moreover, we seek to achieve two partially competing objectives: On the one hand, we have to control the volume of D j , which is covered by Conti constructions, from below. On the other hand, we aim at keeping the resulting overall perimeter of the new covering geometry as small as possible. The construction of a covering, which balances these two objectives, is the content of Proposition 4.16, which is the main result of this section.
Motivated by the heuristic considerations at the beginning of Section 4 (c.f. Figure 8 (a)), we expect that in the cases (P1) and (P2), in which there is no substantial rotation with respect to the relevant Conti construction, the two competing objectives of sufficient volume coverage (Proposition 4.16 (1)) and of a good perimeter bound (Proposition 4.16 (3)), can be satisfied with a surface energy, which is independent of δ j and δ 0 . Indeed, it is possible to show that in the situation without substantial rotation, in each iteration step the overall perimeter of the covering of a triangle is comparable to the perimeter of the original triangle up to a loss of a controlled universal factor.
Proposition 4.16. Let D j be as in (P1)-(P2) with j ≥ 1. Then there exists a covering and a constant C > 1 (independent of δ 0 ) such that: (1) A volume fraction of at least 10 −12 |D j | is covered by finitely many rescaled and translated Conti constructions from Lemma 3.3. The Conti constructions can again be covered by finitely many triangles of the types occurring in the cases (P1)-(P2) and (R1)-(R2), where j is replaced by j + 1. Figure 17. Problems which could arise in the covering algorithm: In our parallel covering result we have to avoid rotated triangles as the aspect ratios are very small for these.
(2) The complement of the Conti constructions is covered by finitely many triangles occurring in the cases (P1)-(P2), where j is replaced by j + 1.
(3) The overall surface energy of the new triangles D j+1,l ∈ D 1 (D k ), is controlled by In the proof of Proposition 4.16 we have to be careful in the choice of the covering, in order to keep all the resulting triangles parallel to the direction of D j or parallel to the relevant Conti construction (c.f. Definition 4.6). This is necessary to ensure a covering such that the sum of the resulting perimeters is comparable to the original perimeter; in particular no factor of δ 0 occurs here. We emphasize that this alignment with the directions of the original triangle or the relevant Conti construction is a central point, since if a (substantial) rotation angle with respect to these directions were to be obtained (e.g. as illustrated in Figure 17, where the covering gives rise to triangles which are rotated by an angle of π 2 ), we would inevitably fall into cases similar as the situations described in (R1), (R2), however with a ratio δ j , which might be substantially smaller than δ 0 . As explained at the beginning of Section 4, this would entail a growth of the perimeters of the covering by a factor δ j . As a consequence our BV estimate from Section 5 would become a superexponential bound, which could no longer be compensated by the only exponential L 1 decay. This would hence destroy all hopes of deducing good higher regularity estimates for the convex integration solutions.
The remainder of this section is organized into three parts: We first discuss the covering constructions for the cases (P1) and (P2) separately in Sections 4.4.1, 4.4.2. Then in Section 4.4.3 we combine these cases, in order to provide the proof of Proposition 4.16. 4.4.1. The case (P1). We begin by explaining the covering in the case (P1).
Lemma 4.17. Let D be a δ j -good triangle oriented along the x 1 -axis. Let R be a rectangle of aspect ratio 1 : δ j /2 such that its long axis is rotated by an angle of β ∈ δ j [−10, 10] with respect to the x 1 -axis. Then there exists a covering by (i) K j , with K j ∈ [1, 100], δ j -good, up to null-sets disjoint triangles, D l , l ∈ {1, . . . , K j }, which are either oriented along the x 1 -axis or the long side of R, (ii) a rescaled and translated copyR of R, such that |R| ≥ 10 −6 |D|.

Moreover,
where C > 1 is a universal constant (in particular independent of δ j and δ 0 ).
Proof. We first invoke Lemma 4.14 with R and δ = δ j . This yields an axis-parallel boxR 2 of side ratio r ∈ [ 1 10 , 10]δ j . This boxR 2 is admissible in Proposition 4.12. An application of this proposition with D,R 2 and δ = δ j hence yields a covering of D by δ j -good triangles, which all have e 1 as their direction, and a box R 2 , which is covered as described in Lemma 4.14. We note that the triangles within R 2 are thus also δ j -good and have as their directions either e 1 or the long side of R. The estimate on the perimeter follows, since all the covering triangles have perimeter controlled by Per(D) and as K j ≤ 100. The estimate on the volume fraction is a consequence of Proposition 4.12 (i) and Lemma 4.14 (i).

4.4.2.
The case (P2). As in the case (P1) we have the following main covering result: Lemma 4.18. Let D be a δ j−1 -good, axis-parallel triangle. Let R be a rectangle of aspect ratio 1 : δ 0 , whose long side is rotated with respect to the axis by an angle β ∈ δ 0 [−1000, 1000]. Then there exists a covering of D into where C > 1 is a universal constant (in particular independent of δ j and δ 0 ).
Proof. We apply Lemma 4.14 with δ = δ 0 and the box R. This yields a boxR 2 of ratio approximately δ 0 and a boxR ⊂R 2 , which is a translated and rescaled copy of R of volume comparable to the volume ofR 2 . Stacking K j := δ0 δj−1 translated copies of the boxesR 2 along the x 1 -axis next to each other and denoting the individual boxes byR 2,k (each containing a translated copyR k ofR), yields as their union a new boxR 2 of aspect ratio 1 : δ j−1 (c.f. Figure 18). With respect to this rectangleR 2 and with δ = δ j−1 we now apply Proposition 4.12, which yields a rectangle R 2 of the same aspect ratio as that ofR 2 . As the volume of eachR k is comparable to the volume ofR 2,k , the claim (ii) of Lemma 4.18 follows from Proposition 4.12, since this ensures that R 2 has volume comparable to D. It remains to bound the perimeters. Here we only estimate the sum of the perimeters of the rectanglesR 2,k , as the remaining parts of the covering are controlled by a multiple of this. We note that 1 δj δj δ0 Figure 18. The stacking construction of the boxes. Each of the smaller boxes is a suitably translated copy of the rectangleR 2 , which is roughly of aspect ratio 1 : δ 0 . In each of these we insert a rescaled and translated version of the construction from Lemma 4.14, c.f. Figure 13.
The main difference of Lemma 4.18 with respect to Proposition 4.12 is the step, in which we bridge the mismatch in the ratios of the triangle D (ratio δ j−1 ) and the given box R (ratio δ 0 ). Here we pass from a box of ratio approximately δ 0 (which is prescribed for R and hence forR 2 ) to a box with ratio approximately δ j (forR 2 ) by stacking translates of the boxesR 2,k next to each other. Proof. The first property of the Proposition follows from Lemma 4.10 in combination with Lemma 4.17 (ii) (in the case (P1)) or Lemma 4.18 (ii) (in the case (P2)). In particular, by Lemma 4.10 all the triangles, which are used to cover the Conti constructions, are δ j+1 -good with respect to the relevant Conti construction. The second property is a consequence of Lemma 4.17 (i) combined with Lemma 4.14 (in the case (P1)) or Lemma 4.18 (i), (iii) (in the case (P2)). We emphasize that all these triangles are either parallel to the original triangle D or to the relevant Conti construction, implying that both the angles and the orientations are within the admissible margins. Finally, the bound on the perimeters follows from the corresponding claims in Lemmata 4.17 and 4.18. 4.5. Covering in the cases (R1)-(R3). In this section we deal with the covering in the cases (R1)-(R3). As in Section 4.4 we seek to simultaneously control the perimeter of the resulting covering and the volume of the domain, which is covered by Conti constructions. Motivated by the discussion from the beginning of Section 4, we however expect that it is unavoidable to produce estimates, in which the ratio δ 0 appears.
With this expectation, we are less careful in our covering constructions and for instance do not seek to preserve the direction n, in which the corresponding δ jgood triangles are oriented. Yet, we still heavily rely on Proposition 4.12 and only modify the construction within the block R 2 . This will give rise to certain new "error triangles", which are of the type (R3). In analogy to Proposition 4.16 we have: Proposition 4.19. Let D j be as in (R1)-(R3) with j ≥ 1. Then there exists a covering and a constant C > 0 independent of δ 0 such that: (1) A volume fraction of 10 −6 |D j | is covered by finitely many rescaled and translated Conti constructions. The Conti-constructions can again be covered by finitely many triangles of the types occurring in the cases (P1)-(P2) and (R1)-(R3), where j is replaced by j + 1. Lemma 4.20. Let β ∈ (Cδ 0 , π 2 −Cδ 0 ). Assume that D is a δ 0 -good triangle oriented parallel to the x 1 -axis. LetR be a rectangle of aspect ratio 1 : δ 0 , which encloses an angle β with respect to the orientation of D. Then D can be covered by the union of (i) M j , with M j ∈ [1, 100], δ 0 -good triangles, D 1,k , which are aligned with the direction of D, (ii) 0 < K j ≤ Cδ −2 0 many translated, up to null-sets disjoint and rescaled copies R 2,k of the rectangleR, whose union covers a volume of size at least 10 −6 |D|, (iii) 0 ≤ L j ≤ Cδ −2 0 many triangles D 2,k , which are of the type (R3).
Here C > 1 is a universal constant. The overall perimeter of the resulting triangles and rectangles is controlled by Proof. Let c 1 ∈ (1/4, 1), c 2 ∈ (1, 4). We begin by stacking K j ∈ [c 1 , c 2 ]δ −2 0 ∩ N many rectanglesR 2,k , which are translated copies ofR, next to each other in such a way that their lowest corners lie on the x 1 -axis (c.f. Figure 19). LetR 2 denote the enveloping axis-parallel rectangle. By adapting the constants c 1 , c 2 we can arrange thatR 2 has an aspect ratio r allowing for an application of Proposition 4.12 with δ = δ 0 ,R 2 and r. This yields a rectangle R 2 of aspect ratio r. Thus, the set D \ R 2 consists of the triangles described in (i). Moreover R 2 is covered as in Figure 19 by K j -many rectangles R 2,k with aspect ratio δ 0 and by a comparable number of Figure 19. The covering of the box R 2 of ratio 1 : δ 0 . The dashed rectangles correspond to the K j stacked (rescaled and translated) copies ofR, which we denote by R 2,k . Their envelope is a (rescaled) copy ofR 2 , which we denote by R 2 . It is the rectangle, which is returned as the output of Proposition 4.12. The parts of R 2 , which are not covered by the rectangles R 2,k , consist of the triangles D 2,k .
"error" triangles D 2,k . By definition, the constant K j satisfies the bounds in (ii). Using elementary geometry, we calculate that This implies the claim of (ii). We note that the error triangles D 2,k are all right angle triangles. Moreover, one of the other angles coincides with the rotation angle β. At least one of triangles' sides is parallel to the orientation of the rectangles R 2,k . Thus the triangles D 2,k are of the type (R3). The bound on the perimeters of the rectangles R 2,k and of the triangles D 2,k results from the following observations: 4.5.2. The case (R2). The case (R2) is the rotated analogue of the case (P2). As we are in a rotated case, we have to be less careful about preserving orientations, and proceed similarly as in the case (R1). Again the main issue is the covering of the rectangle R 2 . However, in contrast to the case (R1) we now have to deal with a mismatch between the ratio of the triangle D (with ratio δ j = δ 0 ) and the ratio of the Conti construction (with ratio δ 0 ). Similarly as in the case (P2) we overcome this issue by a "stacking construction", which compensates the mismatch.
Lemma 4.22. Let β ∈ (Cδ 0 , π 2 − Cδ 0 ). Assume that D is a δ j -good triangle with direction parallel to the x 1 -axis. LetR be a rectangle of side ratio δ 0 , which encloses an angle β with respect to the long side of D. Then D can be covered by the union of 1 δ j ∼ δ 0 δ j Figure 20. The figure shows a (rescaled) copy of the enveloping rectangle R 2 and the stacked (and rescaled) copies of the rectanglē R, which we denote by R 2,k . The triangles correspond to the ones, which we denote by D 2,k in Lemma 4.22 (iii).
There exists a universal constant C > 1 such that the perimeter of the resulting triangles and rectangles is bounded by Proof. We construct a boxR 2 as in the case (R1) but now by stacking K j ∈ [c 1 , c 2 ](δ 0 δ j ) −1 ∩ N many of the boxesR next to each other, where c 1 ∈ (1/4, 1) and c 2 ∈ (1, 4). We denote these stacked boxes byR 2,k and defineR 2 as the enveloping axis-parallel rectangle. By adapting the values of c 1 , c 2 it is possible to obtain a ratio r forR 2 (c.f. Figure 20), which is admissible in applying Proposition 4.12 with δ = δ j ,R 2 and r. This yields a box R 2 , which is covered by rescaled copies R 2,k of the rectanglesR 2,k and by "error" triangles D 2,k . By construction and by elementary geometry (as in (43)) these satisfy the requirements in (ii), (iii). By Proposition 4.12 also (i) holds true. It remains to estimate the perimeter of the union of the rectangles R 2,k and the triangles D 2,k . To this end, we note that: • Each rectangle R 2,k has perimeter controlled by Cδ j Per(D).
• There are Cδ −1 j δ −1 0 many such rectangles R 2,k . • The perimeters of the error triangles D 2,k are up to a factor controlled by the perimeters of the rectangles R 2,k .
Thus, the resulting perimeter is up to a constant bounded by Per(D 2,k ) . This concludes the argument of the lemma. 4.5.3. The case (R3). We deal with the error triangles from the previous step. All of them are right angle triangles, in which the other two angles are bounded from below and above by Cδ 0 and π 2 − Cδ 0 . We show that in this situation we can reduce to two model cases, which we discuss below. This allows us to obtain the following result: Lemma 4.23. Let D be a triangle of type (R3). Let R be a rectangle of side ratio δ 0 , which is parallel to one of the sides of D. Then it is possible to cover D by finitely many scaled and translated copies of itself and by finitely many translated and scaled copies R k of the rectangle R such that Our main ingredient in proving this is the following lemma: Lemma 4.24 (Covering of triangles by rectangles). Let D 1,m denote a right angle triangle, in which the sides enclosing the right angle are of side lengths 1 and m. Assume that m ∈ (0, 50δ −1 ). Let R be a rectangle, which has side ratio δ ∈ (0, 1). Assume that the longer side of R is parallel to the side of the triangle D 1,m , which is of length m.
Then there exist a number L = L(m) and disjoint, rescaled and translated copies R k of R with the properties that: (ii) The sum of the perimeters of the rectangles R k satisfies Remark 4.25. We remark that for our application, the bound on m does not impose an additional requirement. Indeed, the triangles of type (R3) only occur as artifacts of the coverings in Lemmas 4.20, 4.22. Here we may estimate m by tan(β) for the error triangles in . For β = π 2 − δ, a Taylor expansion of sin( π 2 −δ) cos( π 2 −δ) entails the desired estimate.
Before explaining Lemma 4.24, we show how our main covering result, Lemma 4.23, can be reduced to the situation of Lemma 4.24.
Proof of Lemma 4.23. We first claim that without loss of generality D can be assumed to be of type (R3) with R being parallel to one of the short sides of the triangle D. Indeed, if R is parallel to the long side of the triangle D, then this side is opposite of the right angle of the triangle. In this case, we split the triangle D into two smaller triangles D (1) , D (2) by connecting the corner, at which D has its right angle, by the shortest line to the long side. The resulting triangles D (1) , D (2) have the same angles as the original triangle (and in particular satisfy the non-degeneracy conditions for the angles, which are required in condition (R3)), but are now such that R is parallel to one of their short sides. After this reduction, we seek to apply Lemma 4.24 with δ = δ 0 for each of the triangles D (1) , D (2) . To this end we note that as β 0 ≥ Cδ 0 , we have that m ≤ Cδ −1 0 . As a consequence, Lemma 4.24 yields the desired result (by observing that Per(D (1) ) + Per(D (2) ) ≤ 2 Per(D)).
Proof of Lemma 4.24. We construct the desired covering by a "greedy" type algorithm. We begin by fitting in the largest possible copy R 1 of R, which touches the side of the triangle D 1,m , which is of length m (c.f. Figure 21). As the rectangle R 1 has to have a side ratio of 1 : δ, its side lengths can be computed explicitly to be l 1 = m 1+mδ , l 2 = δl 1 . Choosing R 1 as the first rectangle in the desired covering, we have created a decomposition of the triangle D 1,m into three parts, the rectangle R 1 and two triangles, which are self-similar to the original triangle: Here a 1 := l 2 and hence 1 − a 1 = 1 1+δm denote the similarity factors with respect to the original triangle D 1,m . We iterate this procedure in the new triangle D (1−a1),(1−a1)m , while ignoring the (smaller) triangle D a1,a1m . After L steps of this algorithm we have obtained L (up to null sets) disjoint rectangles R 1 , ..., R L .
We claim that if L is chosen sufficiently large, the covering L k=1 R k has the desired properties. Indeed, we choose L such that 1 1+δm L ∈ ( 1 4 , 1 2 ) and first note that the construction of the rectangles R k is based on a self-similar iterative process with similarity factor λ = 1 1+δm . Thus, we infer that Here we used the disjoint construction of the covering, the choice of L, the value of λ and δm ≤ 50. As |D 1,m | = 1 2 m, this yields the first claim. Similarly, Here we used that Per(R 1 ) ≤ 2l 1 and that Per(D m,1 ) ≥ 1.
4.6. Proof of Proposition 4.9. We initialize the construction by applying Step 1 in the Algorithms 3.8, 3.11. As these initial triangles are obtained as the level sets of a Conti construction with ratio δ 0 , they all form δ 0 good triangles (c.f Lemma 4.10) and hence satisfy the properties of the theorem. It therefore remains to argue that this is preserved in our constructions from Sections 4.4 and 4.5. Given one of the triangles D j as in the theorem, the results of Propositions 4.16 and 4.19 ensure this, once the rotation angle of the successive Conti constructions is controlled. This however is the achieved by virtue of Remark 4.2 and Lemma 4.3.

Quantitative analysis
After having recalled the qualitative construction of convex integration solutions in Section 3, we now focus on controlling the scheme quantitatively. Here we rely on the quantitative covering results from Section 4 (c. f. Propositions 4.16 and 4.19), which allow us to obtain bounds on the BV norm of the iterates u k and the corresponding characteristic functions associated with the well e (i) ∈ K (Lemma 5.2). Combined with an L 1 estimate and the interpolation inequality from Theorem 2 or from Corollary 2.1, this then yields the desired W s,q regularity of the characteristic function of the phases.
As in Section 3.2, given a matrix M with e(M ) ∈ intconv(K), we here assume that Ω := Q β [0, 1] 2 , where Q β is the rotation, which describes how the Conti construction with respect to M and e (p) 0 is rotated with respect to the x 1 -axis. This special case will play the role of a crucial building block in the situation of more general domains (c.f. Section 6).
We begin by defining the characteristic functions associated with the corresponding wells: Definition 5.1 (Characteristic functions). We define the characteristic functions, χ k associated with e (1) , e (2) , e (3) in the k-th step of the Conti construction as We denote their point-wise a.e. limits as k → ∞ by χ (i) , i ∈ {1, 2, 3}.
We emphasize that these point-wise limits exists, since for a.e. point x ∈ Ω there exists an index k x ∈ N such x ∈ Ω \ Ω kx . By our convex integration algorithm and by Definition 5.1, the value of χ l (x) remains fixed for l ≥ k x .
Using the covering results from Section 4, we can address the BV bounds for the characteristic functions χ (i) k , i ∈ {1, 2, 3}: denote the deformation and characteristic functions, which are obtained in the j-th step of the convex integration scheme from Proposition 3.16. Let δ 0 be as in Algorithm 3.8, Step 0 (b). Then, there exist constants C 0 , C 1 > 0 (independent of δ 0 ) such that , 2, 3}. Proof. We deduce the following iterative bound for the size of the BV norm: To this end, let Ω j,k be a triangle from the covering {Ω j,k } j∈{1,...,Jj } at the j-th iteration step. In particular, e(∇u j ) is constant on Ω j,k . We apply Algorithms 3.8, 3.11 and in these specify the choice of our covering to be the one of Proposition 4.16 or the one of Proposition 4.19. In order to bound the resulting BV norm, we distinguish two cases: (a) The parallel case. Assume that Ω j,k is of the type (P1) or (P2) (which, by the explanations below Definition 4.7 holds in the parallel case). Thus, the Conti construction in the j-th and (j + 1)-th step are nearly aligned in the direction of their degeneracy. In this case, Proposition 4.16 is applicable and implies that for some absolute constant C > 0. (b) The rotated case. Assume that Ω j,k is of the type (R1)-(R3) (which, by the explanations below Definition 4.7 holds in the rotated case). Thus, the Conti construction in the j-th and (j + 1)-th step are not aligned. By Lemma 4.3 there are even lower bounds on the degree of alignment. In this case, Proposition 4.19 is applicable and yields that Combining both estimates (45), (46) and summing over all domains Ω j,k for fixed j implies (44). From this we infer that Jj k=1 Per(Ω j,k ) ≤ Cδ −j 0 .

As by construction
we therefore obtain the statement of the Lemma.
Using the explicit construction of our convex integration scheme, we further estimate the difference of two successive iterates in the L 1 norm: denote the deformation and characteristic functions, which are obtained in the k-th step of the convex integration scheme from Lemma 3.5. Then, Remark 5.4. In our realization of the covering argument, which is described in Section 4, we have chosen v 0 = 10 −6 . In particular, it is independent of the boundary condition M in (6).
Proof. The proof follows immediately from the Conti construction and the observations that and that χ j (x) for a.e. x ∈ Ω \ Ω j , if k ≥ j. Combining Lemma 5.2 and 5.3 with Theorem 2 or Corollary 2.1 yields the following regularity result: Proposition 5.5 (Regularity of convex integration solutions). Let M with e(M ) ∈ intconv(K) and assume that Ω = Q β [0, 1] 2 . Let δ 0 > 0 be as in Step 0(b) in Algorithm 3.8. Let u : Ω → R 2 be a convex integration solution obtained according to the Algorithms 3.8, 3.11 and described in Proposition 3.16. Then it is possible to obtain BV L q1 Figure 22. By interpolation, our decay and growth bounds for the L 1 and BV norms of the differences χ j yield Cauchy sequences in the W θ,q spaces inside the triangular region. Here, we use that the functions under consideration are characteristic functions and hence all L p norms with 1 ≤ p < ∞ can be compared.
Proof. The proof is along the lines of the proof of Proposition 5.5. However, instead of estimating χ j , we bound ∇u j+1 − ∇u j . Here the BV bound follows from the bound for χ (i) j+1 by noting the uniform boundedness of ∇u j (c.f. Proposition 3.15) and the fact that for the estimate for χ (i) j+1 we used the whole resulting perimeter. Hence For the L 1 estimate we use the L ∞ bound for ∇u j (which follows from Proposition 3.15) and the fact that e(∇u j ) ∈ conv(K)) in combination with the fact that in the j-th iteration step ∇u j is only changed on a volume fraction of 1 − 7 8 v 0 j . Thus, Hence the same interpolation as above yields the W s,p (R 2 ) regularity of ∇u − M , which implies the desired result.

General Domains
In this section we explain how to construct the desired "regular" convex integration solutions in arbitrary Lipschitz domains by using the bounds from the special cases, which were discussed in Section 5. In this context our main result is the following: Proposition 6.1. Assume that Ω ⊂ R 2 is a bounded Lipschitz domain and suppose that M ∈ R 2×2 with e(M ) ∈ intconv(K). Let β ∈ [0, 2π) be the angle, with which the Conti construction for M is rotated with respect to the x 1 -axis and let χ (i) k be defined as in Definition 5.1. Let θ 0 > 0 be the W s,p exponent for the regularity of χ (i) with respect to the rotated unit square Q β [0, 1] 2 adapted to M , i.e. let θ 0 be such that for all s ∈ (0, 1), p ∈ (1, ∞] with 0 < sp < θ 0 and for some µ(s, p) ∈ (0, 1) it holds Then there exists a constant C(Ω, M, s, p) and a family of subsetsΩ k ⊂ Ω such that where Q m l := [0, λ l ] 2 + x l,m are (up to null-sets) disjoint cubes with x l,m ∈ Ω and λ l := 2 −l such that (in the sense of the convergences of their characteristic functions), (49) remains valid for all s, p with s ∈ (0, 1), p ∈ (1, ∞] and 0 < sp < θ 0 . In the dependences the constant C(s, p) however is replaced by C(Ω, M, s, p) and µ(s, p) replaced by C(θ)µ(s, p). Here θ = θ(s, p) is the interpolation exponent associated with s, p > 0.
As an immediate consequence we infer the following corollary: Corollary 6.2. Suppose that Ω ⊂ R 2 is a bounded Lipschitz domain and assume that M ∈ intconv(K). Let β ∈ [0, 2π] be the angle, with which the Conti construction for M is rotated with respect to the x 1 -axis. Let θ 0 > 0 be the limiting W s,p exponent for the regularity of χ (i) with respect to the rotated unit square Q β [0, 1] 2 adapted to M . Then for all s, p > 0 with 0 < sp < θ 0 (i) the point-wise limitχ (i) of the functionsχ Proof of Corollary 6.2. We first note that the point-wise limitχ (i) exists, sincē Ω k → Ω and since χ (i) k → χ k in a point-wise sense as k → ∞. With this at hand, the proof of Corollary 6.2 follows from Proposition 6.1 by interpolation in an analogous way as explained in Proposition 5.5. We therefore omit the details of the proof of the corollary.
We proceed to the proof of Proposition 6.1. Here we argue by covering our general domain Ω by the special domains from Section 5 (Steps 1 and 2). On each of the special domains, we apply the construction from Section 5 (c.f. also Algorithms 3.8, 3.11). In order to obtain a sequence with bounded W s,p norm, we however do not refine to arbitrarily fine scales immediately, but proceed iteratively (c.f. Step 3). A central point here is to control the necessary number of cubes at each scale (Claim 1), since this has to be balanced with the corresponding energy contribution (c.f. Step 4). To this end, we use a "volume argument", which by the Lipschitz regularity of the domain allows us to infer information on the number of cubes on each scale (c.f. Proof of Claim 1).  1)]. By symmetry we may further assume that f (x 1 ) ≥ 0 for all x 1 ∈ [0, 1]. For general domains Ω, by the compactness and Lipschitz regularity, we may locally reduce to a similar case, where f is a Lipschitz curve, but not necessarily a graph. However, all arguments in the following extend to that case as well.
Step 2: Counting cubes. LetΩ l := K l k=1 Q k l , where Q k l ⊂ Ω are (up to zero sets) disjoint, grid cubes of an axis-parallel grid of grid size λ l := 2 −l . We choose K l ∈ N maximal. Thus, by definition we have thatΩ l ⊂Ω l+1 ⊂ Ω for all l ∈ N. In the limit l → ∞ the setsΩ l eventually cover the whole set Ω (which we assume to be as in Step 1). We estimate the number of the cubes, which are contained in the setsΩ l+1 \Ω l . For these we claim: Claim 1. The setΩ l+1 \Ω l contains at most C f λ −1 l+1 of the grid cubes Q k l+1 ⊂ Ω. Proof of Claim 1. Indeed, we first observe that for a sufficiently large constant C f (depending on f , c.f. Remark 6.3) and for a sufficiently large value of l ∈ N every point x ∈ Ω in the subgraph S(f, l) of f − C f λ l is contained in a cube Q k l ⊂Ω l . Hence at least a volume of size |S(f, l) ∩ Ω| : is completely covered by cubes of size λ l . There may be additional cubes of size λ l contained inΩ l . As however only a volume of size C f λ l is left and as each cube has volume λ −2 l , the number of these additional cubes is controlled by #{grid cubes of size λ l inΩ l \ S(f, l)} ≤ 2C f λ l λ −2 l ≤ 2C f λ −1 l . Combining this with the observation that |Ω ∩ S(f, l + 1)| − |Ω ∩ S(f, l)| = C f (λ l − λ l+1 ) = C f λ l+1 , which implies that S(f, l + 1) has at most C f λ −1 l+1 more cubes of size λ l+1 than S(f, l), we infer thatΩ l+1 \Ω l contains at most 4C f λ −1 l+1 cubes of size λ l+1 . Figure 24. The zero-homogeneous exactly stress-free configurations for the hexagonal-to-rhombic phase transition. The white sectors correspond to the variant e (1) , the gray ones to e (2) and the black ones to e (3) : There is a twelve-fold corner, up to symmetry (i.e. rotations by π 2 ) a single variant of a four-fold corner and up to symmetry two different variants of a six-fold corner.

Appendix: A Construction Using Symmetry
In this final section we recall a special solution to (6) with M = 0, which enjoys BV regularity (c.f. [Con08], [Pom10], [CPL14]). This construction crucially relies on symmetry. It hence gives rise to the question whether it is possible to exploit symmetry in a more systematic way in constructing "regular" convex integration solutions.
In describing this particular solution, we first present all possible zero homogeneous solutions to the differential inclusion (4) in Section 7.1. In Section 7.2 we then rely on this to construct the desired solution with BV regularity. 7.1. Exactly stress-free configurations for the hexagonal-to-rhombic phase transformation. The high degree of non-rigidity of the hexagonal-to-rhombic phase transition is reflected in a comparably large number of possible solutions to (4). There is already a large number of solutions with homogeneous strain, i.e. solutions u of (4) such that all the phases intersect in a single point and the strains are zero-homogeneous functions e(∇u)(λx) = e(∇u)(x) (c.f. Figure 24). To verify this, we recall, c.f. Lemma 17 in [Rül16], that the following conditions are necessary and sufficient for the presence of such a corner.
Lemma 7.1 (Compatibility condition at a zero-homogeneous corner). Let e : R 2 → R 2×2 sym be a zero-homogeneous tensor field. Let A 1 , . . . , A m ∈ R 2×2 sym be symmetric matrices such that A j = A j+1 , where j, j + 1 are considered modulo m. Assume that along a closed circle surrounding the origin, e successively attains the values A 1 , . . . , A m (as for instance in Figure 24). Then e is a strain tensor, i.e. there exists a function u ∈ W 1,∞ loc (R 2 ) such that e = e(∇u), if and only if the following two conditions are satisfied: (1) There exist vectors a i ∈ R 2 \ {0}, n i ∈ S 1 such that A i − A i+1 = 1 2 (a i ⊗ n i + n i ⊗ a i ) for i ∈ {1, ..., m}.
(2) m i=1 a i ⊗ n i = 0. Figure 25. The homogeneous corners can be combined to yield the above patterns (up to symmetry they correspond to a single pattern). As before white corresponds to the variant e (1) , gray to the variant e (2) and black to the variant e (3) .
Here the first condition corresponds to tangential continuity along the interfaces of the jumps. The second requirement ensures the compatibility of the skew symmetric part of the gradient.
Keeping this in mind, a symbolic Mathematica computation allows to determine all possible zero-homogeneous corners. This leads to the following classification result: Observation 1. Apart from simple laminates (which trivially exist due to the symmetrized rank-one connectedness of the strains) there are the following compatible zero-homogeneous configurations of strains (c.f. Figure 24): • a single configuration involving twelve strains, • (up to symmetry) two types of configurations involving six strains, • (up to symmetry) one configuration involving four strains.
These homogeneous configurations can further be combined to yield compatible "zig-zag" configurations, c.f. Figure 25. The zero-homogeneous configurations will in the sequel serve as the building blocks of our constructions. 7.2. A BV construction for zero boundary data. In this section we present an explicit construction of a solution to (6) with M = 0 and e(∇u) ∈ BV (and with ∇u ∈ L p (R 2 ) for all p ∈ (1, ∞) but ∇u / ∈ L ∞ (R 2 )). In our construction (c.f. Proposition 7.2) we crucially rely on symmetry properties of the strains. As the wells form an equilateral triangle in strain space, it appears plausible to expect the best regularity properties arise in the center of the convex hull of K, i.e. for solutions with zero boundary data.
For an arbitrary domain this construction (which is motivated by the constructions of Conti [Con08] and which was similarly already used in [Pom10] and [CPL14]) can be used to obtain a solution to (6) with e(∇u) ∈ W s,p (R 2 ) for all s ∈ (0, 1), p ∈ (1, ∞) with sp < 1.
As we are in two-dimensions, symmetrized rank-one connections exist between any pair of symmetric matrices with vanishing trace (c.f. Lemma 2.8). In particular, all the strains e (1) , e (2) , e (3) are compatible with any (constant) boundary condition. Yet, a priori it is not clear, whether this can be turned into a global configuration involving as few strains as possible. In the case of zero boundary data, this is indeed possible, while preserving very good regularity properties: Proposition 7.2 (Zero boundary data construction). There exists a bounded domain Ω ⊂ R 2 , Ω = ∅, and a solution u : R 2 → R 2 of (6) with M = 0, such Figure 26. The first iteration step in the deformation corresponding to zero boundary data. It is possible to interpret the construction as (a linearization at the identity of) the deformation depicted in the first two pictures: The undeformed reference configuration (left) is deformed into the configuration on the right. that e(∇u) ∈ BV (R 2 ).
We emphasize that this construction is not new and that related constructions, using the symmetry of the domain, already appeared in [Con08], [Pom10], [CPL14].
Proof. In order to obtain the desired construction, we consider an equilateral triangle rotated by π 12 with respect to the x 1 -axis and a self-similar copy of it which is homothetically positioned in the larger one at a length ratio 4 − 2 √ 3 and then rotated by π 3 with respect to the outer triangle (c.f. Figure 26, left). Now we rotate the inner triangle by π 3 , so that it turns into a homothetically scaled version of the outer triangle and stretch it by a factor 2 while preserving the boundary of the larger triangle (c.f. Figure 26, right). This leads to the gradient distribution depicted in Figure 27. In particular, the only strains, which are used consist of e (1) , e (2) , e (3) as well as the zero strain. This can be iterated in the respectively smaller triangles. As only the skew symmetric part of the strains grow, while the symmetric strains are fixed in the set of our wells and the zero matrix, and as the new interior triangle is a self-similar copy of the outer original triangle with a constant ratio of 2(4 − 2 √ 3), this yields the claimed energy contributions: In each step the construction leads to a bound of the form where c ∈ (0, 1) and C > 1 is a universal constant. Thus, passing to the limit n → ∞ implies the desired regularity result for the symmetrized part.
Remark 7.3. We remark that this can easily be turned into a scaling result for associated elastic and surface energies on the domain Ω. It gives an energy scaling for minimizers subject to ∇u = 0 on the boundary of Ω, which corresponds to a "surface energy contribution".
A second way of producing the construction of Proposition 7.2 relies on the homogeneous building blocks, which were described in the previous section. We • compute all the possible four-fold corners consisting of the strains e (1) , e (2) , e (3) , 0. There is an admissible four-fold corners with a large portion (more precisely involving an angle of 5π 3 ) of the zero phase (this is one of the corner depicted in the left picture in Figure 26), • act on this configuration via a rotation of 2π 3 : This yields two new compatible four-fold corners (these are the other two corners in the left picture in Figure 26), • combine the corners in an equilateral triangle as in Figure 26, • check the compatibility (by means of Lemma 7.1) of the resulting central corners with the strain e = 0.
Similar constructions (but with different symmetries) can be applied in other two-dimensional configurations in matrix space (e.g. with the symmetries of a square), if the strains are arranged in a symmetric polygon in strain space (e.g. in square, in which case one possible solution would be the one given in Lemma 3.1, c.f. also Figure 4).
In a general domain the construction from above can be applied with a (logarithmic) loss: Proposition 7.4 (General zero boundary data construction). Let Ω ⊂ R 2 be a non-empty Lipschitz domain. Then there exists a configuration such for all s ∈ (0, 1), p ∈ (1, ∞) with sp < 1 we have e(∇u) ∈ W s,p (R 2 ).
As the passage from the special domain to an arbitrary domain follows by a covering argument analogous to the one presented in Section 6, we omit the proof here.
Remark 7.5. Similarly as explained in Remark 7.3, Proposition 7.4 can also be transformed into a scaling result. In comparison to the "surface energy scaling", which is obtained from Proposition 7.2, we however lose a logarithmic factor in the corresponding construction. It is an interesting and challenging open problem to decide whether this logarithmic loss is necessary in a general domain.