Convex Integration Arising in the Modelling of Shape-Memory Alloys: Some Remarks on Rigidity, Flexibility and Some Numerical Implementations
Angkana Rüland Jamie M. Taylor Christian Zillinger
math.AP


Flexibility, Rigidity and Some Numerical Implementations]Convex Integration Arising in the Modelling of Shape-Memory Alloys: Some Remarks on Rigidity, Flexibility and Some Numerical Implementations







Max-Planck-Institute for Mathematics in the Sciences, Inselstraße 22, 04103 Leipzig, Germany
rueland@mis.mpg.de



Department of Mathematics,
University of Southern California,
3620 S. Vermont Avenue,
Los Angeles, CA 90089-2532, US
zillinge@usc.edu


Kent State University - Kent Campus, Department of Mathematical Sciences, Mathematics and Computer Science Building 233, Summit Street, Kent OH 44242, US
jtayl139@kent.edu




We study convex integration solutions in the context of the modelling of shape-memory alloys. The purpose of the article is two-fold, treating both rigidity and flexibility properties: Firstly, we relate the maximal regularity of convex integration solutions to the presence of lower bounds in variational models with surface energy. Hence, variational models with surface energy could be viewed as a selection mechanism allowing for or excluding convex integration solutions. Secondly, we present the first numerical implementations of convex integration schemes for the model problem of the geometrically linearised two-dimensional hexagonal-to-rhombic phase transformation. We discuss and compare the two algorithms from <cit.> and <cit.>.



[
    Christian Zillinger
    December 30, 2023
=======================






§ INTRODUCTION



Shape-memory alloys are materials undergoing a first order, diffusionless solid-solid phase transformation in which symmetry is lost. Here typically the high temperature phase, the austenite, is highly symmetric (e.g. having atomistically a highly symmetric unit cell structure), while the low temperature phase, the martensite, loses some of this symmetry. This loss of symmetry amounts to the presence of multiple, energetically equivalent variants of martensite. This renders these materials rather complex and gives rise to many complicated, interesting microstructures in these materials <cit.>.





 §.§ The Ball-James variational model


Mathematically, various features of the microstructures in shape-memory alloys have been quite successfully modelled by (static) energy minimization <cit.>, <cit.>:

    min_y(x)=M(x)  on ∂Ω∫_Ω W(∇ y, θ) dx .

Here Ω⊂^3 denotes the reference configuration (which is often chosen to be the material configuration in the austenite phase at fixed temperature), y:Ω→^3 denotes the deformation of the material which is subjected to the boundary conditions y(x)=M(x) on ∂Ω.
The function θ:Ω→_+ represents temperature (in normalized units) and W: ^3× 3_+ →_+ the stored energy function. Physical requirements on W are:


  (i) frame indifference: W(QM,θ)=W(M,θ) for each M∈^3× 3_+, Q ∈ SO(3) and θ∈_+.

  (ii) material symmetry: W(MH, θ)= W(M,θ) for each M∈^3× 3_+, θ∈_+ and H ∈𝒫, where 𝒫⊂ SO(3) denotes the point group of the material, which encodes the symmetries of the material. 

The symmetry of the phase transformation in conjunction with these two conditions then allows one to describe the energy wells of W for a martensitic phase transformation

    W(M,θ)=0 
    ⇔ M ∈ K(θ)= {[           α(θ)SO(3) θ>θ_c,; ⋃_j=1^m SO(3)U_j(θ) θ≤θ_c, ].
     U_j(θ)^t= U_j(θ).

Here θ_c ∈_+ denotes the transformation temperature between the austenite and the martensite, α(θ)SO(3) with α: _+ →_+ models the austenite phase at temperature θ with its thermal expansion α(θ), and the matrices SO(3)U_j(θ) model the (energetically) equivalent variants of martensite at temperature θ∈_+ <cit.>.
The two conditions (i), (ii) render the minimization problem (<ref>) highly non-(quasi)convex and hard to analyse mathematically (even at a fixed temperature, which is the usual mathematical set-up) <cit.>. 



 §.§ The m-well problem


As a consequence of the difficulty of dealing with the full problem (<ref>), if one is interested in studying low energy microstructures, it is often useful to investigate exactly stress-free deformations and microstructures. These are the solutions to the differential inclusion

    ∇ y ∈ K(θ) Ω ,  y(x) = M(x) ∂Ω,

where K(θ) is the set from (<ref>). If θ<θ_c, (<ref>) is also called an m-well problem (due to the presence of m variants of martensite). A natural question here is the existence of solutions to (<ref>), respectively, the existence of energy zero solutions to (<ref>). While highly oscillatory minimizing sequences with limiting energy zero had been shown to exist for a large class of boundary data (related to the quasiconvex hull of K(θ); in general these sequences however only have measure valued limits), <cit.>, the existence of exact energy zero solutions was only resolved later by relying on the technique of convex integration <cit.>. These quite surprising solutions however are rather “wild" and oscillatory in that they use a cascade of different scales (similar to convex integration solutions in fluid mechanics and geometry, c.f. <cit.> and the references therein). This “wildness" can be quantified in terms of the regularity of the solutions to (<ref>). In particular, it has been shown in different models <cit.> that convex integration solutions can only exist at rough regularities, giving rise to a dichotomy. For instance for the two-well problem <cit.> (where there are exactly two rank-one connections between the wells) the following is known: 


  (a) If a solution is BV regular, i.e. if ∇ y ∈ BV, then ∇ y is a simple laminate, i.e. it is (up to boundary effects) “one-dimensional" with a directional dependence which is determined by the structure of the set K(θ) from (<ref>).

  (b) For any fixed θ<θ_c and M ∈ (K(θ)^lc), where K(θ)^lc denotes the lamination convex hull of K(θ), there exist solutions to (<ref>). By default these solutions are L^∞ regular, i.e. ∇ y ∈ L^∞.

Thus, a natural question is whether there is a threshold between these two regimes in terms of a critical regularity threshold (similar to the case in fluid mechanics and geometry <cit.>). In contrast to the situation in fluid mechanics, the problem (<ref>) does not have a natural scaling, which would help to indicate a critical threshold. In previous work together with B. Zwicknagl <cit.>, we have hence studied the possibility of deriving higher regularity for the solutions from (b) for simplified “linearised" models. It is one of the purposes of this article to show that scaling laws yield lower bounds on the possible regularity of these solutions (and could be viewed as “viscosity approximations" of the limiting model). Moreover, we present numerical simulations of convex integration solutions for the simplest possible model setting. We explain this in more detail in the following.



 §.§ Surface energy regularisation and scaling laws

While various important qualitative properties, e.g. the emergence of preferred twinning directions as well as the presence of clear volume fractions of certain variants of martensite, have been successfully explained by (<ref>) and (<ref>), more quantitative properties such as the emergence of length scales are not captured by this. In the literature (c.f. for instance <cit.>) this has been remedied by adding a higher order surface energy contribution to (<ref>), which is weighted with a small (material dependent) prefactor ϵ∈_+:

    min_y(x)=M(x)  on ∂Ω( ∫_Ω W(∇ y, θ) dx + ϵ^2 ∫_Ω|∇^2 y|^2 dx ).

From a mathematical point of view the presence of the higher order term in (<ref>) regularises the problem and restores compactness to it. In general, the exact physical form of the higher order term in (<ref>) is not known and various types of different regularisations have been used in the literature, including the diffuse, L^2 variant stated in (<ref>), but also the sharp, BV variant (c.f. <cit.> for an overview on different models for this). In a sense (<ref>) can be viewed as a selection criterion for certain microstructures that arise in (<ref>) (and (<ref>)); the additional surface term in (<ref>) distinguishing between wild, rough and highly oscillatory (approximate) solutions in (<ref>) and (approximate) solutions to (<ref>), which are more regular (and hence use less surface energy). The properties of minimisers are often analysed by studying scaling limits of the functional (<ref>), c.f. for instance <cit.>. While not predicting the exact form of the selected minimisers in (<ref>), scaling does yield important information on these.
Recently, a more precise limiting analysis has been carried out in <cit.>. In spite of the various impressive advances in this direction, many points remain open in this context.



 §.§ The main results

The purpose of this note is two-fold and addresses both rigidity and flexibility properties of the underlying differential inclusions: Firstly, we show that scaling laws complement the analysis of convex integration solutions in a precise sense, in that the presence of a scaling law (or rather a lower bound on the energy (<ref>)) yields an upper bound on the maximal regularity of convex integration solutions (which might be well-known to experts in the field, but which we could not find in the literature). Here we exploit that, while (<ref>) itself does not have a prescribed scale predicting a possibly critical regularity, the model (<ref>) has.



Let Ω⊂^n be an open bounded C^1,1 domain. Let E_ϵ be the energy

    E_ϵ
    = min_∇ u = M  a.e. in ^n∖Ω{∫_Ω^2(∇ u, K) dx 
    + ϵ^2 ∫_Ω|∇^2 u|^2 dx  },

where K=K(θ) for some fixed θ∈_+ is as in (<ref>).
Assume that there exist constants C>1 and μ∈ (0,1/2) such that for all ϵ∈ (0,ϵ_0)

    ϵ^2μ≤ C E_ϵ.
 
Suppose that u is a solution to 

    ∇ u    ∈ K Ω,
    ∇ u    = M ^n ∖Ω.

If v(x):=u(x)-Mx-b ∈ H^1+s(^n) for some b∈^n, s∈ with (v)⊂Ω and ∇ v ∈ L^∞(^n), then s≤μ.

Loosely speaking this implies that for ∇ u on an H^s_loc(^n) scale, any convex integration solution ∇ u ∈ L^∞(Ω) can be at most H^μ_loc(^n) regular.

We remark that the bound μ∈ (0,1/2) is related to the presence of trace estimates (c.f. the proof of Theorem <ref> below).


Similar results can be obtained by using different forms of surface energies (c.f. Propositions <ref> and <ref> in Section <ref>) and the corresponding Sobolev and Besov spaces. 


Secondly, in addition to providing this upper bound on the maximal regularity of convex integration solutions in terms of the behaviour of scaling laws, we present numerical implementations of the convex integration schemes discussed in <cit.> and <cit.> in the case of the model setting of the geometrically linearised hexagonal-to-rhombic phase transformation (c.f. Section <ref> for a more detailed description of this phase transformation). 



 §.§ Outline of the article

The remainder of the article is organised as follows: First, in Section <ref>, we address the relation between scaling laws and rigidity results. In particular, we present the proof of Theorem <ref> (and discuss variants of it). Next, in Section <ref>, we recall the convex integration algorithms from <cit.> and <cit.> for the model case of the hexagonal-to-rhombic phase transformation. Numerical implementations of these are presented in Section <ref>, where we also compare the main features of the two algorithms. Finally, in Section <ref> we summarise our findings and discuss interesting questions and models related to these. In the Appendix, we provide a proof of a fractional Poincaré inequality in the form in which we use it in Section <ref>.



 §.§ Acknowledgements
 
The research of Jamie M Taylor leading
to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no 291053. Christian Zillinger acknowledges the support of an AMS-Simons Travel Grant. 




§ RIGIDITY AND SCALING LAWS



In this section we present the derivation of the bound on the maximal possible regularity of convex integration solutions. Here we do not yet specify the set-up to that of a particular phase transformation but remain in the general framework, which was layed out in Section <ref>.


We explain the close relation between scaling laws and the possible regularity of convex integration solutions (which is reminiscent of viscosity approximations, c.f. <cit.> for these ideas in the context of fluid mechanics):
Assume that for a bounded Lipschitz domain Ω⊂^n a function u ∈ W^1,∞(Ω), which in addition satisfies ∇ u ∈ H^γ_loc(^n), is a solution to the differential inclusion

    ∇ u    ∈ K Ω,
    ∇ u    = M ^n ∖Ω,

where K:=K(θ) denotes the wells from (<ref>) for some fixed temperature θ>0 (this corresponds to the case of affine boundary data in (<ref>), which is already an interesting case). 
Then, any sequence ∇ u_k in H^γ_loc(^n) with ∇ u_k →∇ u in L^2_loc(^n) and with ∇ u_k = M a.e. in ^n ∖Ω is an admissible competitor in the minimization problem for E_ϵ, where E_ϵ is the energy from (<ref>). The lower bound in the scaling law from Theorem <ref>, then yields a bound on the maximal regularity of ∇ u.

The special role of the space H^μ (and the W^s,p Sobolev  exponent ps= 2 μ) can be guessed if we suppose that the scaling law E_ϵ(∇ u) ∼ϵ^2μ holds. Assuming an equipartition of energy and treating the first term in (<ref>) heuristically as an L^2 energy, we obtain that for a minimiser u_ϵ with fixed ϵ>0

    ∇ u_ϵ_L^2^2    ∼ϵ^2μ,  ∇^2 u_ϵ_L^2^2 ∼ϵ^2μ-2.

Hence if there is a critical function space for ∇ u, it should be determined by the requirement that its scaling behaviour is independent of ϵ>0. By interpolation, we thus guess that it is given by ∇ u ∈ [ L^2, H^1]_μ = H^μ.

Similar arguments indicate “critical" function spaces if the energies in (<ref>) are not diffuse, but for instance sharp interface models.

This intuition can be made rigorous and directly yields the proof of Theorem <ref> (c.f. Section <ref>).



 §.§ Function spaces



Before proceeding with the discussion of the relation between regularity and scaling, for reference we first collect the function spaces, which we will be using in the sequel.
Denoting the Fourier transform by 
w (k) = ∫_^n e^i k · xw(x) dx, for s∈ we consider the L^2 based fractional Sobolev spaces

    H^s(^n):={u∈𝒟'(^n):  u_H^s(^n):= ^-1[(|ξ|^2 +1)^s/2 u]_L^2(^n)<∞}.


Moreover, we will also use Besov spaces. To this end we work with a Littlewood-Paley decomposition. Following <cit.>, we let φ denote a non-negative, radially symmetric bump function supported in {k ∈^d:  |k|≤ 2}, which is equal to one on {k ∈^d:  |k|≤ 1}. Let N ∈ 2^, i.e. assume that there exists j∈ such that N=2^j. We then define the Littlewood-Paley projectors P_<N, P_≥ N, P_N to be 

    (P_< N f)(k)   := φ(k/N) f(k),
    (P_≥ N f)(k)   := (1-φ(k/N)) f(k),
    (P_N f)(k)   := (φ(k/N)- φ(2k/N)) f(k).

With this at hand, we
recall that for s∈ and p∈ [1,∞) the space of Besov functions f ∈ B^s_p,p(^n) is defined as

    B^s_p,p(^n)={f∈𝒟'(^n):  f_B^s_p,p(^n)<∞},

where

    f_B^s_p,p(^n):= ( ∑_N∈ 2^ N^spP_N f_L^p(^n)^p)^1/p,

c.f. for instance <cit.> or <cit.>.
For later use, we recall that B_1,1^0(^n) ⊂ L^1(^n) (c.f. Theorem 2.41 in <cit.>).
For s> 0,  s∉ this space coincides with the fractional Sobolev space W^s,p(^n), which is defined by the norm

    u_W^s,p(^n)^p = u_L^p(^n)^p + ∇^k u_L^p(^n)^p + ∫_^n∫_^n|u(x)-u(y)|^p/|x-y|^n+σ p dy dx,

where s=k+σ, k∈, σ∈ (0,1) (c.f. <cit.> for a comparison of several fractional Sobolev type function spaces).






 §.§ Proof of Theorem <ref>



We begin the discussion of the relation between scaling and regularity by presenting the proof of Theorem <ref>.


We argue by regularization (through a frequency cut-off in the corresponding Littlewood-Paley decomposition) and a spacial cut-off (in order to ensure the validity of the boundary conditions).


Let u be a solution to (<ref>). As a preliminary step, to simplify the set-up, we shift the problem (<ref>) by introducing v(x)=u(x)-Mx-b for b∈^n, which solves

    ∇ v    ∈K̃Ω,
    ∇ v    = 0 ^n ∖Ω,

where K̃=K-M:={ N-M:  N ∈ K}. The constant b∈^n is chosen such that v ∈ H^s+1(^n) and (v) ⊂Ω. We seek to show that s≤μ.
Rewritten in terms of the shifted functions v, the energy (<ref>) turns into

    Ẽ_ϵ= min_∇ v = 0  a.e. in ^n ∖Ω∫_Ω^2(∇ v, K̃) dx + ϵ^2 ∫_Ω |∇^2 v|^2 dx.
 
For P_<N, P_≥ N denoting the standard Littlewood-Paley projector (c.f. the definition in (<ref>)), in the sequel, we use the following notation:

    (∇ u)_<N:= P_<N(∇ u), (∇ u)_≥ N:= P_≥ N(∇ u).

Given the function v, we cut-off its high frequencies and correct the frequency cut-off of v so that it satisfies the desired boundary conditions, i.e. we consider the function

    ṽ_N(x):= η_δ(x) v_<N(x) ,

where η_δ is a smooth, positive cut-off function which is equal to one in Ω_2δ:={x∈Ω: (x,∂Ω)>2δ}, which vanishes outside of Ω_δ:={x∈Ω: (x,∂Ω)>δ} and which satisfies

    |∇η_δ| ≤C/δ,  |∇^2 η_δ| ≤C/δ^2.

The value of the parameter δ will be determined in the sequel. Moreover, here and in the following arguments, the constant C>0 will be generic and might change from line to line, but will always be independent of δ, N, ϵ.
We observe that for v with ∇ v ∈ L^∞(^n)∩ H^s(^n) Bernstein estimates (c.f. the Appendix in <cit.>, and Section 2.1.1. in <cit.>) yield

    ∇^2 v_<N_L^2(Ω)   ≤∇^2 v_<N_L^2(^n)≤  N^1-s∇ v_H^s(^n),
    (∇ v)_≥ N_L^2(Ω)   ≤(∇ v)_≥ N_L^2(^n)≤ N^-s∇ v_H^s(^n).



We now argue in two steps and first estimate the size of (∇ṽ_N, K̃)_L^2(Ω)^2. We show that this quantity is controlled by Cϵ^2s for some uniform constant C>1. Afterwards we prove a similar bound for ϵ^2∇^2 ṽ_N_L^2(Ω)^2.  From this we will derive a bound on the possible regularity of ∇ v (and hence for ∇ u).


Step 1: Estimate for (∇ṽ_N,K̃)_L^2(Ω)^2.
Using the triangle inequality, we estimate:

    ( ∇ṽ_N, K̃) _L^2(Ω)^2
       ≤ C((η_δ∇ v_<N, K̃) _L^2(Ω)^2 + v_<N∇η_δ_L^2(Ω)^2) 
       ≤ C((∇ v, K̃) _L^2(Ω)^2 
    + ((1-η_δ)∇ v, K̃)_L^2(Ω)^2 .
         . + ∇ v - ∇ v_<N_L^2(Ω)^2 
    + v_<N∇η_δ_L^2(Ω)^2 
    )
       ≤ C(δ+ v_<N∇η_δ_L^2(^n)^2 
    + ∇ v - ∇ v_<N_L^2(^n)^2
    ).

Here we used that v is a solution to (<ref>) and that since ∇ v ∈ L^∞(Ω) and since K̃ is bounded, the support condition for η_δ implies

    ∫_Ω^2((1-η_δ)∇ v, K̃) dx ≤ C δ.

We estimate the two terms on the right hand side of (<ref>) separately: First, we apply the Bernstein estimates from (<ref>) to infer

    ∇ v - ∇ v_<N_L^2(^n)^2
     = ∇ v_≥ N_L^2(^n)^2
     ≤ C N^-2s∇ v_H^s(^n)^2.

Next, we estimate

    v_<N∇η_δ_L^2(^n)^2≤ C (v ∇η_δ_L^2(^n)^2 + v_≥ N∇η_δ_L^2(^n)^2).

In order to bound the second term on the right hand side of (<ref>), we note that by the bounds for η_δ (c.f. (<ref>)) and Bernstein estimates

    v_≥ N∇η_δ_L^2(^n)^2
       ≤ Cδ^-2v_≥ N_L^2(^n)^2
    ≤ C δ^-2 N^-2-2s∇ v_H^s(^n)^2.

For the first term on the right hand side of (<ref>), using that (∇η_δ) ⊂Ω_δ∖Ω_2δ, we obtain 

    v ∇η_δ_L^2(^n)^2 
    ≤ C δ^-2v_L^2(Ω_δ∖Ω_2δ)^2
    ≤ C  ∇ v_L^2(Ω∖Ω_2δ)^2
    ≤ C δ^2s∇ v_H^s(^n)^2.

Indeed, the last estimate in this chain of inequalities can be inferred by using the support assumption for ∇ v and a fractional Poincaré inequality (c.f. Lemma <ref> in Section <ref>).

Combining the bounds from (<ref>)-(<ref>), we thus obtain as the estimate for the elastic energy

    (∇ṽ_N, K̃)_L^2(Ω)^2
    ≤ C (N^-2s + δ^-2N^-2-2s + δ^2s)∇ v_H^s(Ω)^2 + C δ.

Choosing N = δ^-1= ϵ^-1, hence implies the bound

    (∇ṽ_N, K̃)_L^2(Ω)^2
    ≤ C (ϵ^2s + ϵ).




Step 2: Conclusion.
Next we estimate the surface energy:

    ∇^2 ṽ_N_L^2(Ω)^2
       ≤ C (∇^2 v_<N_L^2(^n)^2 + v_<N∇^2 η_δ_L^2(^n)^2 + ∇η_δ∇ v_<N_L^2(^n)^2)
       ≤ C(N^2-2s∇ v_H^s(^n)^2 + v_<N∇^2 η_δ_L^2(^n)^2 + ∇η_δ∇ v_<N_L^2(^n)^2).

We again estimate the two last terms separately: Using a fractional Poincaré inequality once more in combination with Bernstein estimates, we obtain for the third term in (<ref>)

    ∇η_δ∇ v_<N_L^2(^n)^2
       ≤ C (∇η_δ∇ v_L^2(Ω)^2 + ∇η_δ∇ v_≥ N_L^2(Ω)^2)
       ≤ C δ^-2δ^2s∇ v_H^s(^n)^2 + C δ^-2N^-2s∇ v_H^s(^n)^2.

For the second term in (<ref>), we similarly note by Poincaré and Bernstein

    v_<N∇^2 η_δ_L^2(Ω)^2
       ≤ C v∇^2 η_δ_L^2(^n)^2 
    + C v_≥ N∇^2 η_δ_L^2(^n)^2
       ≤ C ∇ v_H^s(^n)^2 (δ^-2+2s + δ^-4N^-2-2s).

Combining (<ref>)-(<ref>) we therefore obtain

    ϵ^2∇^2 ṽ_N_L^2(Ω)^2
    ≤ Cϵ^2(N^2-2s+δ^-2+2s+δ^-4N^-2-2s + δ^-2N^-2s)∇ v_H^s(^n)^2.

Finally, combining (<ref>) with (<ref>) and choosing N^-1=δ = ϵ, we conclude

    ϵ^2μ≤ C (∇ṽ_N, K̃)^2_L^2(Ω) + C ϵ^2 ∇^2 ṽ_N_L^2(Ω)^2 ≤ C ϵ^2s∇ v_H^s(^n)^2 + C ϵ.

As ∇ v_H^s(^n)^2 is bounded independently of ϵ, we observe that, for ϵ∈ (0,1) sufficiently small and μ∈ (0,1/2), (<ref>) can only hold as long as s ≤μ. This concludes the argument.






The restriction μ∈ (0,1/2) enters rather naturally and is connected to the presence/absence of a trace. For most (though not for all, c.f. Section 7 in <cit.>) convex integration solutions, one would expect that no trace estimates are available, the prescribed boundary data enforcing too strong oscillations in ∇ u close to the boundary ∂Ω. For the convex integration solutions in <cit.> the presence of a trace for instance provides the sharp obstructions for the existence of convex integration solutions.






 §.§ Variations on Theorem <ref>



Instead of working in L^2 based function spaces, we can also work with a Littlewood-Paley decomposition and the associated Bernstein estimates in other function spaces, i.e. for different energy functionals.




Consider for instance the (sharp interface) energy functional

    E_ϵ,1
    = min_∇ u = M  a.e. in ^n ∖Ω{∫_Ω^2(∇ u, K) dx 
    + ϵ |∇^2 u|(Ω)},

where |∇^2 u|(Ω) denotes the variation measure of the distributional derivative of ∇ u.

As in Theorem <ref> we obtain the following bound for possible convex integration solutions:



Let E_ϵ,1 be the energy from (<ref>). Assume that there exist constants C>1 and μ∈ (0,1/2) such that for all ϵ∈ (0,ϵ_0)

    ϵ^2μ≤ C E_ϵ,1.
 
Suppose that u is a solution to (<ref>).
If for some b ∈^n and s∈ we have
v(x)=u(x)-Mx-b ∈ B_1,1^1+s(^n), (v)⊂Ω and ∇ v ∈ L^∞(^n), then s≤ 2μ.


Here B^1+s_1,1(^n) denotes the standard Besov space (defined for instance through a suitable Littlewood-Paley decomposition, c.f. Section <ref> and <cit.>). Loosely speaking the proposition states, that on a Besov scale, we can at most have ∇ u ∈ (B^2μ_1,1)_loc(^n).


As a first preliminary step, we again switch to the situation where

    v(x) = u(x)-Mx-b,

and b∈^n is chosen such that v∈ B_1,1^1+s(^n) and (v)⊂Ω.
The modified energy wells turn into K̃:= K-M, and the energy (<ref>) becomes 

    Ẽ_ϵ,1:= min_∇ v = 0  a.e. in ^n ∖Ω{(∇ v, K̃)_L^2(Ω)^2 + ϵ |∇^2 v|(Ω) }.


Given a convex integration solution u of (<ref>) or v of (<ref>), we again consider the functions

    ṽ_N(x):= η_δ(x) v_<N(x) ,

where η_δ is a smooth, positive cut-off function which is equal to one in Ω_2δ:={x∈Ω: (x,∂Ω)>2δ}, which vanishes outside of Ω_δ:={x∈Ω: (x,∂Ω)>δ} and which satisfies

    |∇η_δ| ≤C/δ,  |∇^2 η_δ| ≤C/δ^2.

The value of the parameter δ will be determined in the sequel. 
Using Bernstein type inequalities <cit.> and the embeddings B_1,1^0(^n)⊂ L^1(^n) (c.f. Theorem 2.41 in <cit.>) and W^1,1(Ω) ⊂ BV(Ω), 
for v with ∇ v ∈ L^∞(^n)∩ B^s_1,1(^n) we deduce

    |∇ (∇ v)_<N|(Ω)
       ≤(∇ v)_<N_Ẇ^1,1(^n)≤  C(∇^2 v)_<N_B^0_1,1(^n)≤  CN^1-s∇ v_B^s_1,1(^n),
    (∇ v)_≥ N_L^1(Ω)   ≤(∇ v)_≥ N_L^1(^n)≤ C(∇ v)_≥ N_B^0_1,1(^n)≤ CN^-s∇ v_B^s_1,1(^n).

With this at hand, using that

    ∇ṽ_N = η_δ∇ v_<N + v_<N∇η_δ,

and
using the properties of the cut-off η_δ (c.f. (<ref>)), we estimate as follows:

    (∇ṽ_N, K̃)_L^2(Ω)^2
    ≤ 2(η_δ∇ v_<N, K̃)_L^2(Ω)^2 
    + 2v_<N∇η_δ_L^2(Ω)^2
       ≤ 2(η_δ∇ v_<N, K̃)_L^2(Ω)^2 
    + 2v_<N∇η_δ_L^∞(^n)v_<N∇η_δ_L^1(^n)
       ≤ 2(η_δ∇ v_<N, K̃)_L^2(Ω)^2 
    + 2C_ν( 1/δ N^2-ν+δ)∇ v_L^∞(^n)^2,

where ν >0 is an arbitrarily small constant.
In order to infer this, we used the bounds

    v_<N∇η_δ_L^q(^n)   ≤ C( v ∇η_δ_L^q(^n) 
    +  v_≥ N∇η_δ_L^q(^n)
    )

for q∈{1,∞}. By Bernstein estimates, we control the second term in (<ref>) for q=∞

    v_≥ N∇η_δ_L^∞(^n)   ≤ C δ^-1 N^n/pv_≥ N_L^p(^n)≤ C δ^-1 N^n/p N^-1|∇| v_≥ N_L^p(^n)
       ≤ C δ^-1 N^n/p N^-1∇ v_≥ N_L^p(^n)≤ C δ^-1 N^n/p N^-1∇ v_L^p(^n)
       ≤ C δ^-1 N^-1+ν/2∇ v_L^∞(^n),

with ν/2:= n/p (and p>1 can be chosen arbitrarily large),
and for q=1

    v_≥ N∇η_δ_L^1(^n)≤∇η_δ_L^∞(^n)v_≥ N_L^∞(^n)|Ω_δ∖Ω_2δ|
           ≤ C 1/δ1/N^1-ν/2∇ v_L^∞(^n)δ = C/N^1-ν/2∇ v_L^∞(^n).

The respective first terms on the right hand side of (<ref>) are controlled by invoking the fundamental theorem

    v ∇η_δ_L^∞(^n)   ≤ C δ^-1v _L^∞(Ω_δ∖Ω_2δ)≤ C δ^-1δ∇ v_L^∞(Ω∖Ω_2δ) ,
    v ∇η_δ_L^1(^n)   ≤ C δ^-1v _L^1(Ω_δ∖Ω_2δ)≤ C δ^-1δ∇ v_L^1(Ω∖Ω_2δ)≤ C  δ∇ v_L^∞(Ω∖Ω_2δ) .

Further, we observe that

    (η_δ∇ v_<N, K̃)_L^2(Ω)^2
    ≤ C(η_δ∇ v_<N, K̃)_L^1(Ω) (∇ v_<N_L^∞(Ω)+C_K̃)
       ≤ C((η_δ∇ v, K̃)_L^1(Ω) 
    + ∇ v - ∇ v_<N_L^1(Ω))
    (∇ v_<N_L^∞(^n)+C_K̃)
       ≤ C(( ∇ v, K̃)_L^1(Ω) + ((1-η_δ) ∇ v, K̃)_L^1(Ω) 
    + ∇ v - ∇ v_<N_L^1(^n))
           × (∇ v_<N_L^∞(^n)+C_K̃)
       ≤ C(δ
    + N^-s∇ v_≥ N_B^s_1,1(^n))(N^ν/2∇ v_L^∞(^n)+C_K̃)
       ≤ C(δ
    + N^-s∇ v_B^s_1,1(^n))(N^ν/2∇ v_L^∞(^n)+C_K̃).

Combining (<ref>)-(<ref>), we infer for the elastic energy

    (∇ṽ_N, K̃) _L^2(Ω)^2 
       ≤  C(δ N^ν/2
    + N^-s+ν/2∇ v_B^s_1,1(^n)+Cδ^-1N^-2+ν)
         × 
    (∇ v_L^∞(^n)+C_K̃) .

For the surface energy we estimate 

    |∇^2 ṽ_N|(Ω)
       ≤∇^2 v_<N_L^1(Ω) + 2 ∇η_δ∇ v_<N_L^1(Ω) +  v_<N∇^2 η_δ_L^1(Ω)
       ≤∇^2 v_<N_L^1(Ω) + 2 ∇η_δ∇ v_≥ N_L^1(^n)
    + 2 ∇η_δ∇ v _L^1(^n)
           + v_≥ N∇^2 η_δ_L^1(^n)
    + 2 v ∇^2 η_δ_L^1(^n)
       ≤ C(∇^2 v_<N_L^1(Ω) + 1/δ∇ v_≥ N_L^1(^n) + 1/δ^2 v_≥ N_L^1(^n)
          + 1/δ∇ v_L^1(Ω∖Ω_2δ) + 1/δ^2 v_L^1(Ω∖Ω_2δ))
    
       ≤ C( N^1-s∇ v_B^s_1,1(^n) + 2/δ N^s∇ v_B^s_1,1(^n) + 1/δ^2 N^1+s∇ v _B^s_1,1(^n)
          + 2 δ^s/δ∇ v_B^s_1,1(^n) + δ^1+s/δ^2∇ v_B^s_1,1(^n)).

Here we have used some of the Bernstein type estimates from (<ref>) in combination with a fractional Poincaré inequality (c.f. Section <ref>).
Combining (<ref>), (<ref>) and the assumed lower bound for E_ϵ,1 and choosing N=δ^-1 yields

    ϵ^2μ   ≤Ẽ_ϵ,1≤(∇ṽ_N,K)_L^2(Ω)^2 + ϵ |∇^2 ṽ_N|(Ω)
       ≤ C(ν,∇ v_L^∞(^n), ∇ v_B^s_1,1(^n)) (N^-1+ N^-s+ν/2 + N^-1+ν + ϵ N^1-s).


Choosing N ∼ϵ^-1, we obtain

    ϵ^2μ   ≤Ẽ_ϵ,1≤  C(ν,∇ v_L^∞(^n), ∇ v_B^s_1,1(^n))( ϵ^s - ν/2 + ϵ^1-ν).

As a consequence, by choosing ν depending on μ sufficiently small and considering the limit ϵ→ 0, we infer that on a Besov scale no (L^∞-) convex integration solution can be more than B^2μ_1,1(^n)-regular for μ∈ (0,1/2).



More generally, the previous situation can be studied with regularisations in the Sobolev spaces W^s,p. 
In this setting we have the following result:



Let 1< p ≤ 2 and denote by E_ϵ,p the energy 

    E_ϵ,p
    = min_∇ u = M  a.e. in ^n ∖Ω{∫_Ω^2(∇ u, K) dx 
    + ϵ^p ∇^2 u_L^p(Ω)^p }.

Assume that there exists constants C_p>1 and μ∈ (0,1/2) such that for all ϵ∈ (0,ϵ_0)

    ϵ^2μ≤ C_p E_ϵ,p.
 
Suppose that u is a solution to (<ref>). If for some b∈^n and s∈ it holds that v(x)=u(x)-Mx-b ∈ W^s,p(^n) with (v)⊂Ω and ∇ v ∈ L^∞(^n), then s≤2μ/p.


In other words, the scaling behaviour yields an upper bound on the regularity of convex integration solutions, i.e. on a Sobolev W^s,p_loc scale, we can at most have ∇ u ∈ W^2μ/p,p_loc(^n)∩ L^∞(^n).

The proof proceeds as in the BV setting by using Bernstein type estimates.


As in the previous proofs, we reduce to the shifted problem by setting v(x)=u(x)-Mx-b for which we assume that for some b∈^n we have v∈ W^s,p(^n) and (v) ⊂Ω. Hence, if u is a solution to (<ref>), then v is a solution to (<ref>). The correspondingly modified energy reads

    Ẽ_ϵ,p= min_∇ v = 0  a.e. in ^n ∖Ω{∫_Ω^2(∇ v, K̃)dx + ϵ^p ∇^2 v_L^p(Ω)^p }.

As in the proof of Proposition <ref>, we use the standard Littlewood-Paley decomposition to define

    ṽ_N(x):= η_δ(x) v_<N(x),

and estimate for N=δ^-1

    (∇ṽ_N, K̃)_L^2(Ω)^2
    ≤ C((∇ v_<N,K̃)_L^2(Ω)^2 + v_<N∇η_δ_L^2(Ω)^2 + δ)
       ≤ C((∇ v_<N,K̃)_L^2(Ω)^2 + v_<N∇η_δ_L^1(^n)v_<N∇η_δ_L^∞(^n) + δ)
       ≤ C_ν((∇ v_<N,K̃)_L^2(Ω)^2 + δ^1-ν∇ v_L^∞(Ω)^2 + δ),

where we used similar estimates as in (<ref>), (<ref>) in the proof of Proposition <ref>.
Furthermore,

    (∇ v_<N,K̃)_L^2(Ω)^2
       ≤(∇ v_<N, K̃)_L^p(Ω)^p (∇ v_<N, K̃)_L^∞(Ω)^2-p
       ≤ C((∇ v, K̃)_L^p(Ω) + ∇ v- ∇ v_<N_L^p(^n))^p 
         × (N^ν∇ u_L^∞(^n)^2-p+C_K̃)  
       ≤  C ∇ v_≥ N_L^p(^n)^p (N^ν∇ u_L^∞(^n)^2-p+C_K̃)
       ≤ C N^-ps∇ v_W^s,p(^n)^p (N^ν∇ u_L^∞(^n)^2-p+C_K̃).

Here we used that by Bernstein estimates, Calderón-Zygmund estimates, and the compact support of ∇ v,

    ∇ v_<N_L^∞(^n)^2-p≤ C_ν N^ν∇ v_<N_L^p(^n)^2-p≤ C_ν N^ν∇ v_L^p(^n)^2-p≤ C_ν N^ν∇ v_L^∞(^n)^2-p

for ν = ν(p)>0 being arbitrarily small.
Thus, for N^-1=δ the elastic energy is controlled by

    (∇ṽ_N, K̃)_L^2(Ω)^2
    ≤ C( δ^ps - ν∇ v_W^s,p(^n)^p +  δ^1-ν∇ v_L^∞(^n)^2 + δ)

with C=C(K̃, ν, n, ∇ v_L^∞(^n)).
For the surface energy we estimate (by using Poincaré and Bernstein)

    ∇^2 ṽ_N_L^p(Ω)^p
       ≤ C (∇^2 v_<N_L^p(Ω)^p + ∇η_δ∇ v_<N_L^p(Ω)^p + v_<N∇^2 η_δ_L^p(Ω)^p)
       ≤ C (∇^2 v_< N_L^p(^n)^p + ∇η_δ∇ v_≥ N_L^p(^n)^p + v_≥ N∇^2 η_δ_L^p(^n)^p
          +  ∇η_δ∇ v_L^p(^n)^p + v∇^2 η_δ_L^p(^n)^p)
       ≤ C (N^p-sp+δ^-pN^-sp + δ^-2pN^-sp-p
            + δ^-pδ^sp+ δ^-2pδ^p+sp)∇ v_W^s,p(^n)^p
       ≤ C δ^-p+sp∇ v_W^s,p(^n)^p,

where we have again chosen δ=N^-1.
Thus, by recalling the assumed lower bound for the energy, we infer

    ϵ^2μ≤Ẽ_ϵ,p≤ C(ν, ∇ v_L^∞(Ω), ∇ v_W^s,p(Ω))( ϵ^p N^p-sp+ C N^-1+ν + C N^-1 + N^-ps+ν),

where C≥ 1.
Further setting N=ϵ^-1, we hence deduce that

    ϵ^2μ   ≤Ẽ_ϵ,p≤ C (ϵ^sp - ν  + ϵ^1-ν).

As a consequence, for μ∈ (0,1/2) on a W^s,p Sobolev scale no convex integration solution can be more than W^2μ/p,p-regular (which can be observed by choosing ν depending on μ sufficiently small and by passing to the limit ϵ→ 0).




Similarly as in <cit.> and in <cit.>, we have thus obtained that if a family of energy functionals is controlled, then the possible W^s,p Sobolev regularity is determined in terms of the product sp and not by s,p individually. Here the L^∞ bound was crucial, as it allowed us to “break scaling".





The work of Dacorogna, Marcellini and Paolini on origami constructions
<cit.> illustrates that in some cases it is possible to saturate the sharp regularity threshold originating from scaling laws by convex integration solutions: Considering the energy E_ϵ,1, in their case, it is possible to approach μ→1/2 and self-similar convex integration solutions of BV regularity exist. This shows that in their very flexible setting of the inclusion ∇ u ∈ O(2) (or also ∇ u ∈ O(n) and zero boundary conditions) the only obstruction is the presence or absence of trace estimates.





§ FLEXIBILITY



In this section, we recall the convex integration results from <cit.>, <cit.> in the model case of the geometrically linearised hexagonal-to-rhombic phase transformation (which is discussed in Section <ref>). In particular, after recalling the general outline of the convex integration algorithm in Section <ref>, we discuss the two different convex integration schemes, which were used in <cit.> and <cit.> to derive the higher regularity results (c.f. Sections <ref>-<ref>). In Section <ref> we will then discuss (non-quantitative) numerical implementations of the schemes and compare these.




 §.§ The geometrically linearised hexagonal-to-rhombic phase transformation



We recall the model and some properties of the geometrically linearised hexagonal-to-rhombic phase transformation. 



A difficulty of the variational model from Section <ref> is the two-fold “nonlinearity" of the problem, manifested in the two physical requirements (i), (ii). This is also still reflected in the m-well problem from Section <ref>, in that the set K(θ) still displays SO(3) invariance, and has a multi-well structure at temperatures below θ_c. Hence, in order to simplify this, it is often convenient to “linearise" the frame indifference assumption and hence to pass from SO(3) to Skew(3) invariance (we recall that Skew(3) is the linearisation of SO(3) at the identity). Mathematically, this has the advantage of dealing with an invariance which is given by a vector space structure in contrast to a nonlinear group structure, while at the same time preserving the “material nonlinearity", i.e. the multi-well structure. For one-well problems this has been rigorously justified in <cit.>. In the small rotation regime, it is expected that this linearisation still captures important features of the nonlinear original problem (<ref>) (although care is required, in particular in the large rotation regime, c.f. <cit.>). 

In the sequel, we study such a geometrically linearised problem in two-dimensions. Passing formally from the deformation y(x) to the displacement u(x)=y(x)-x and dropping all “higher order terms" (c.f. the discussion in <cit.>), this turns the nonlinear problem (<ref>) (at fixed temperature θ<θ_c) into a linearised m-well problem of the type

    e(∇ u):= 1/2(∇ u + (∇ u)^t) ∈{e^(1),…,e^(m)} + Skew(3).

Here e(∇ u) denotes the (infinitesimal) strain tensor and e^(1), …, e^(m)∈^3× 3_sym represent the
variants of martensite. In the sequel, we specify the transformation to be the two-dimensional hexagonal-to-rhombic phase transformation (which is related to transformations occurring in materials such as Mg_2Al_4Si_18, Mg-Cd alloys or Pb_3(VO_4)_2, see <cit.>, <cit.>), c.f. Figure <ref> for a microscopic two-dimensional derivation of the deformation matrices (the third direction can be ignored as the material only undergoes an affine change there). Assuming that the material undergoes at most an affine deformation in the third direction, allows us to reduce the three-dimensional problem to a two-dimensional one. Hence, in the following we study the differential inclusion

    e(∇ u):= 1/2(∇ u + (∇ u)^t) ∈{e^(1),e^(2),e^(3)} + Skew(2),

where now u:Ω⊂^2 →^2 and 

    e^(1) = [  1  0;  0 -1 ],  
    e^(2) = 1/2[   -1 √(3); √(3)    1 ],  
    e^(3) = 1/2[    -1 -√(3); -√(3)     1 ].

Note that the traces of all three matrices in (<ref>) vanish, which corresponds to the modelling assumption that the transformation is (infinitesimally) volume preserving.

We recall two important features of the model (<ref>), (<ref>), which were discussed in more detail in Section 2 in <cit.> (c.f. also Figure <ref>):



Let e^(1), e^(2), e^(3) be as in (<ref>). Then, for each pair i,j ∈{1,2,3}, with i≠ j, there exist (up to a change of sign) exactly two vectors n_ij^1, n_ij^2∈ S^1 and two vectors a_ij^1, a_ij^2∈^2 ∖{0} such that 

    e^(i)-e^(j) = 1/2(a_ij^k⊗ n_ij^k + n_ij^k ⊗ a_ij^k)=: a_ij^k⊙ n_ij^k,  k ∈{1,2}.





This lemma in particular implies that for every pair e^(i), e^(j) with i ≠ j there exist simple laminates. These are one-dimensional microstructures in which the strain only attains the values e^(i), e^(j), i.e. there exists ν∈ S^1 such that

    e(∇ u)(x) = f(x ·ν)  f(x·ν) ∈{e^(i), e^(j)}.

Moreover, (up to a change of sign) ν∈{n_ij^1, n_ij^2}. These are however not the only possible piecewise affine microstructures, which can occur in the hexagonal-to-rhombic phase transformation. For a complete list of “homogeneous" deformations we refer to Section 7 in <cit.>. Moreover, we remark that these “homogeneous" deformations can also be concatenated into more complex microstructures.

As a consequence of the presence of the rank-one connections between the wells, the transformation displays a very flexible behaviour. This is manifested in the size of its convex hulls (c.f. Figure <ref>):



Let e^(1), e^(2), e^(3) be as in (<ref>). Then,

    {e^(1), e^(2), e^(3)}^lc = 
    ({e^(1), e^(2), e^(3)}).

In particular, 

    ({e^(1), e^(2), e^(3)}^lc) = 
    (({e^(1), e^(2), e^(3)}))=2.





For a discussion of the argument leading to this result we refer to Section 2 in <cit.> and <cit.>. In particular, in full matrix space, the inclusion (<ref>) is a co-dimension one problem (the only constraint coming from the trace condition). Hence, there is enough room to construct convex integration solutions:



Let Ω⊂^2 be an open Lipschitz domain.
Let e^(1),e^(2),e^(3) be as in (<ref>). Then for any M∈^2× 2 with

    e(M):=1/2(M+M^t) ∈({e^(1),e^(2),e^(3)})
 
there exists a solution u of (<ref>). Moreover, it is possible to ensure that for any ϵ>0 there is a solution u_ϵ to (<ref>) with

    u_ϵ-Mx_L^∞(Ω)< ϵ.



This (non-quantitative) result follows for instance from the arguments of Müller-Šverák <cit.> or from the results of Dacorogna-Marcellini <cit.>. As seen from the approximation result (<ref>) these solutions are highly non-unique; in <cit.> this is made more precise in a Baire category sense.
We remark that in invoking the arguments from <cit.> and <cit.>, there is a slight subtlety here, in that the inclusion is formulated on the level of the symmetrised gradient and not on the level of the full gradient. However, there are various ways of overcoming this, one being to simply “pull up" the inclusion problem to an inclusion problem for the full gradient.

In <cit.> and <cit.> together with B. Zwicknagl we analysed the underlying construction schemes more precisely and showed that it is possible to improve the regularity of solutions to (<ref>) on a (W^s,p) Sobolev scale by choosing the underlying lengths scales carefully. The iterative algorithms in <cit.> and <cit.> however differed quantitatively, in that they produced different dependences: While the regularity of solutions constructed by the scheme from <cit.> depended on the position of the boundary data M in matrix space, the scheme from <cit.> could produce solutions with a regularity which did not depend on this (or more precisely, where we could give bounds which were independent of this). 

In the sequel, we discuss and compare explicit numerical implementations of these schemes, which illustrate the differences which occur here. While convex integration solutions in the context of elasticity have been known for quite a while, the present note seems to contain the first numerical implementation of these. In our numerics, we focus on the qualitative convex integration result from Theorem <ref> and do not seek to optimize the underlying partitions, which was necessary in the quantitative analysis in <cit.>, <cit.>. This is due to the introduction of highly fractal structures which would exceed the capability of our computers (and which from a certain level onwards would also not be seen as major changes without zooming into the structures). Also, we believe that, while being analytically convenient, the exact coverings from <cit.>, <cit.> are far from optimal. Since the numerical implementation without this already provides interesting insights, we opted to focus on the problems without the additional layer of (analytical and numerical) difficulty originating from the complicated covering structures from <cit.>, <cit.>. 




 §.§ The outline of the convex integration scheme for the hexagonal-to-rhombic phase transformation



We recall the convex integration schemes from <cit.>, <cit.> applied to the hexagonal-to-rhombic phase transformation. In the following section, we then discuss and compare some numerical implementations of these. 

We begin by formulating the rough outline of the convex integration algorithms which are used in <cit.> and <cit.> in terms of “pseudo-code". In their main structure the two algorithms are similar, however there are important differences which are mainly encoded in the functions, which are used in Step (1b).



Let M ∈ℳ:={N ∈^2× 2:  e(N)∈({e^(1), e^(2), e^(3)})}, where e^(1), e^(2), e^(3) are as in (<ref>). Let Ω⊂^2 be a triangle.


  (1a) Variables. We consider 


  * the displacement u_k: Ω→^2 at step k, 

  * a collection of (up to null-sets disjoint) triangles Ω_k = {Ω_1^k,…,Ω_j_k^k}, which cover Ω, 

  * and the error in matrix space ϵ_k: Ω_k →  (0,1) at step k, which is constant on each subset of Ω_k.


  (1b) Functions. 
We consider a “covering function", which covers a given triangle by “good sets", on which the deformation u_k will be improved, and a “remainder". More precisely,

    _v:    𝒯×ℳ→ℛ×𝒯,
       (T,M) ↦ ({Q_1,…,Q_j(M,T)} , {T_1,…, T_l(M,T)}).

Here 𝒯 is the set of all triangles, ℛ is a set of certain quadrilaterals (these differ in <cit.> and <cit.>) and ℳ:={M∈^2× 2:  e(M) ∈({e^(1), e^(2), e^(3)}) }. For each T ∈𝒯 and M ∈ℳ we have an (up to null sets) disjoint covering

    T = ⋃_Q ∈_v(T,M)[1] Q ∪⋃_T̃∈_v(T,M)[2]T̃,

with

    |⋃_Q ∈_v(T,M)[1] Q |/|T|≥ v.

The sets Q ∈_v(T,M)[1] are the “good sets", on which the current displacement gradient ∇ u_k will be modified and pushed towards the energy wells. The sets T̃∈_v(T,M)[2] are the “remainders" on which the displacement u_k is not changed.

We further consider a “replacement function", which improves the current displacement gradient in the sense that the replaced deformation gradient is closer to the wells (or at least closer to the wells on a large portion of the domain). This depends on the current displacement gradient, the underlying domain and the error in matrix space. As an output it yields 


  * the level sets (in the form of a finite collection of triangles) of the new improved deformation, 

  * a piecewise affine function whose symmetric gradient attains values in ({e^(1), e^(2), e^(3)}),

  * and an updated error in matrix space. 

More precisely,

    :     ℛ×𝒜_({e^(1), e^(2), e^(3)})× (0,1) 
         →𝒯×𝒜_({e^(1), e^(2), e^(3)})× (0,1),
          (Q,w_old,ϵ) ↦ ({T_1,…,T_j_0} , w, ϵ̃)

Here 𝒜_({e^(1), e^(2), e^(3)}) denotes the set of piecewise affine deformations with symmetric gradients in {e^(1),e^(2),e^(3)}. Furthermore, w|_∂ Q=w_old|_∂ Q.

  (2) Initialization. We begin by setting u_0(x)=Mx, ϵ_0 = ϵ_0(M), Ω_0 = {Ω}. 

  (3) Iteration step. The algorithm proceeds iteratively: Assume that u_k, ϵ_k and Ω_k with k≥ 0 are already given.
Then on each Ω_j^k ⊂Ω_k for which  e(∇ u_k)|_Ω_j^k∉{e^(1), e^(2), e^(3)}, 
apply the function _v(Ω_j^k, ∇ u_k|_Ω_j^k). Let _v(Ω_j^k,∇ u_k|_Ω_j^k)[1]= {Ω_j,1^k,…,Ω_j,m(j,k)^k}.

For each Ω_j,m^k ∈_v(Ω_j^k,∇ u_k|_Ω_j^k)[1] apply the function 

(Ω_j,m^k, u_k |_Ω_j^k,ϵ_k). This yields 


  (i) an up to null-sets disjoint covering of Ω_j,m^k into triangles 

    {Ω_j,m,1^k,…,Ω_j,m,l(j,m,k)^k}:=(Ω_j,m^k, u_k|_Ω_j^k,ϵ_k)[1];


  (ii) a function v_j,m,k:=(Ω_j,m^k, u_k|_Ω_j^k,ϵ_k)[2]: Ω_j,m^k →^2 whose gradient is constant on each of the sets Ω_j,m,l^k with l∈{1,…, l(j,m,k)} and for which v_j,m,k(x)=u_k(x)  x ∈∂Ω_j,m^k;

  (iii) a parameter ϵ̃_j,m,k:=(Ω_j,m^k,  u_k |_Ω_j^k,ϵ_k)[3].

We then set

    u_k+1|_Ω_j,m^k=v_j,m,k u_k+1|__v(Ω_j^k,∇ u_k|_Ω_j^k)[2]=u_k|__v(Ω_j^k,∇ u_k|_Ω_j^k)[2],
       Ω_k+1 = ⋃_j,m(Ω_j,m^k,  u|_Ω_j^k,ϵ_k)[1] ∪⋃_j_v(Ω_j^k,∇ u_k|_Ω_j^k)[2],
       ϵ_k+1|_T = {[                                               ϵ̃_j,m,k T ⊂Ω_j,m,l^k ∈_v(Ω_j^k,∇ u_k|_Ω_j^k)[1]  j ∈{1,…,j_k},;                                                 ϵ_k|_T            T ∈_v(Ω_j^k,∇ u_k|_Ω_j^k)[2]  j ∈{1,…,j_k}. ].

If on Ω_j^k we already have e(∇ u_k)|_Ω_j^k∈{e^(1), e^(2), e^(3)}, we set u_k+1|_Ω_j^k = u_k|_Ω_j^k and

    Ω_k+1=Ω_k+1∪{Ω_j^k∈Ω_k: ∇ u_k|_Ω_j^k∈{e^(1), e^(2),e^(3)}}.




It has been shown in <cit.> and <cit.> that this procedure is well-defined, if the function  is chosen appropriately (there is a slight subtlety in that in those articles we approximate domains by rectangles, but this does not matter for the non-quantitative algorithm which is used here).
In spite of their similar overall structure, the algorithms from <cit.>, <cit.> however differ substantially in their underlying replacement constructions encoded by the function . We discuss the different constructions for the function  in Sections <ref> and <ref>. In Section <ref> we also present the details on the function _v.




 §.§ The replacement constructions



As they form the core of the convex integration algorithm, we discuss the replacement constructions, which are encoded in the function  in the Algorithm <ref>, in more detail. Here we used different options in <cit.> and <cit.>. In spite of this, the main aim of the replacement construction is the same: Given a prescribed gradient distribution, we seek to modify and improve it in the sense that it becomes “closer" to being an exact solution of the differential inclusion (<ref>), (<ref>).

On the one hand, in the article <cit.>, we achieved this by a “piecewise affine" construction in which for a given point x∈Ω the gradient distribution ended up being exactly in the well after a finite (but x-dependent) number of steps. However, on parts of the domain, the gradient distribution deteriorated, in the sense that on a part of the domain, the gradient was pushed away from the wells.
On the other hand, in <cit.> in general, we had to modify the gradient distribution for each point countably many times, but improved the overall gradient distribution in each step essentially in a “uniform way". We describe this in more detail in the following two sections. 



  §.§.§ The replacement construction from <cit.>



In <cit.> the replacement construction (i.e. the function ) is based on the following lemma, which is a modification of a construction due to Conti <cit.> (c.f. Figure <ref>):



Let M∈^2× 2 with 

    e(M):= 1/2(M+M^t) ∈(K)

and let e^(i) with i∈{1,2,3} (as in (<ref>)) be such that

    |e(M)-e^(i)| ≤(e(M),K) + 4 ϵ_0,

where ϵ_0∈ (0,(e(M), ∂(K)))/100. 
Then there exists a rectangular domain Ω_δ(M) (which in general is
rotated with respect to the coordinate axes) and a Lipschitz function u:^2
→^2, with symmetric gradients e^(i)∈ K, ẽ_1, …, ẽ_4 ∈(K) such that 


  (i) e(∇ u) ∈{e^(i), ẽ_1, …, ẽ_4}⊂(K) Ω_δ(M) and

    |{x∈Ω_δ(M): e(∇ u)(x)=e^(i)}|= 1/4|Ω_δ(M)|.


  (ii) u(x)=Mx ^2 ∖Ω_δ(M).





The distribution of the symmetric gradients ẽ_1,…, ẽ_4 in strain space is illustrated in Figure <ref>: On a fixed volume fraction of the domain the modified gradient is pushed exactly into one of the wells, but on the other parts of the domain it either remains close to the original boundary datum, or is pushed even further into the interior of the equilateral triangle spanned by e^(1), e^(2), e^(3).





In the context of Algorithm <ref> this implies the following: Given
the current displacement u_k, we modify it on each of its level sets
Ω_j^k by first covering Ω_j^k by finitely many (translations and
rescalings) of the domains Ω_δ(∇ u_k|_Ω_j^k) from
the Lemma <ref> and a remainder, which consists of finitely many
triangles (this can always be achieved, c.f. the discussion of the function
_v in Section <ref>). The number of domains in the covering is
chosen such that a sufficiently large volume fraction v∈ (0,1) is covered
(v=1 would be ideal, but in general, this would require countably many sets, which is not
numerically feasible; in the quantitative convex integration algorithm from <cit.> we also excluded this, as it would yield infinite surface energy in general). On the (translations and rescalings of the) domains Ω_δ(∇ u_k) we invoke Lemma <ref>. 

In particular, the construction from Lemma <ref> is such that the displacement u_k is replaced by a displacement u_k+1 whose symmetric gradient is already exactly contained in one of the wells on a subset of Ω_δ(∇ u_k) of a given volume fraction (which is in particular uniformly bounded from below, e.g. by 1/8). This set is no longer touched in the remainder of the Algorithm <ref>. Hence, the volume fraction of the original domain, on which the displacement gradient is modified in the k-th step, is bounded by (1-v/8)^k. In particular, for almost any point x∈Ω there exists an integer k(x) such that the gradient ∇ u_k(x) is already exactly in the wells after k(x) iteration steps. 
On the complement of this set, i.e. on the set, where the displacements are not yet in the energy wells, we push a part of the gradient distribution further into the interior of {e^(1), e^(2),e^(3)} and on another part we remain in an ϵ_k neighbourhood of the original e(∇ u_k)|_Ω_j^k. In particular, on these parts of the domain, we either do not improve or even deteriorate the underlying displacement construction (c.f. Figure <ref>).
This is iterated.

Due to the “pushing into the interior", after a finite number of steps, the closest well changes (i.e. the choice of e^(i) in (<ref>) changes). This corresponds to a change in the rank-one direction in matrix space, which in turn is reflected in the orientation of the underlying domains Ω_δ(∇ u_k). In particular, here any possible orientation (in the union of certain cones, c.f. the blue cones in Figure 3 in <cit.>) may arise in general, they are not prescribed by the rank-one directions (c.f. Lemma <ref>) of the energy wells {e^(1), e^(2), e^(3)}. 


The key features of the replacement construction in <cit.> can hence be summarized as follows:


  * It is “finite", in the sense that for almost every point x∈Ω there exists k(x)∈ such that the iteration creates a gradient distribution for which e(∇ u_k)(x) already attains one of the values e^(1), e^(2) or e^(3) after a finite number k(x) of steps. After this has been achieved the gradient distribution at that point will not be altered anymore.

  * It is “non-uniform", in that by using the “pushing out" construction the gradient distribution is deteriorated on part of the domain (i.e. moved further away from the wells) in each iteration step.

  * It possibly involves a continuum of possible orientations and normals.







  §.§.§ The replacement construction from <cit.>



In the algorithm from <cit.> the replacement construction is based on the following observation:



Let M ∈^2× 2 be such that e(M)∈(K). Assume that
(M,K)=d_0. Then, for a constant γ∈ (0,1), there exists a domain Ω_δ(M), which is diamond-shaped (but possibly rotated), and a piecewise affine map w:Ω_δ(M) →^2 such that


  (i) (∇ w, K)≤γ d_0,

  (ii) w(x)=Mx on ∂Ω_δ(M).





Regarding the constant γ∈ (0,1), we remark that it can be chosen in a uniform way in the sense that there exists σ∈ (0,1) (independent of M; almost 1/2 up to ϵ errors) such that the distance between ∇ w and K improves by a factor σ, where the distance is measured in
l^∞ with respect to the convex coefficients of M= λ_ie^(i).





The construction of Lemma <ref> is illustrated in Figure <ref>. We
remark that the directions along which ∇ w is modified and along which Ω_δ(M) is oriented correspond to the
rank-one directions from Lemma <ref>. In particular, only
finitely many of these directions occur. 



As before, in the application of Lemma <ref> in Algorithm <ref>, we first
cover a sufficiently large volume fraction of the domains Ω_j^k by
finitely many translations and rescalings of the sets Ω_δ(∇
u_k|_Ω_j^k) ∈_v(Ω_j^k, ∇ u_k|_Ω_j^k)[1]. On these we apply the
replacement construction from Lemma <ref>. We note that in contrast
to the scheme from Section <ref>, the new piecewise affine
deformation is in general not exactly contained in the energy wells on
any of its level sets – not even on a part of the underlying domain. In
particular, the volume of the domain on which the displacement gradients have to
be improved in the next step, does not decrease; in the above sense this
algorithm is not “finite" but “countable". However, in contrast to the
previous scheme, we now improve the gradient distribution “uniformly" and (essentially)
dyadically on all of its level sets (on which it is modified), which yields geometric convergence of the
gradient distribution. On the modified domain the gradient distribution is
always improved by a factor of (nearly) two with respect to an l^∞-based metric on the barycentric coordinates (c.f. property (i) in Lemma <ref>).


Hence the properties of this scheme can be summarized as follows:


  * It is “countable", in the sense that in general the construction does not push the gradient distribution exactly into one of the wells within a finite iteration time.

  * It is “uniform", in that the gradient distribution is always improved by a factor two in the whole domain on which the function  is acting on. 

  * It involves a finite set of possible normals and orientations.


We emphasize that in both replacement constructions we obtain good control on the skew part. While this plays a role in the quantitative argument in <cit.>, it is not of major importance in our qualitative scheme. Hence, we do not discuss this further here.



 §.§ The covering constructions



As the form of the replacement function 
has already been discussed in the previous section, it remains to describe our implementation
of the function _v. Given a triangle T ∈𝒯 and a matrix M ∈ℳ, i.e. e(M)∈({e^(1), e^(2), e^(3)}), we essentially proceed by a greedy algorithm. We cover T by a grid of dyadic scales
consisting of dyadic rescalings and translations of the domains Ω_δ(M) from Lemmas <ref> and <ref>. If a grid parallelogram is a strict subset of T, we keep it and add it to the
list _v(T, M)[1]. If a parallelogram is a strict subset of T^c, we discard the parallelogram. For all parallelograms intersecting the boundary of T, we dyadically refine the parallelogram, and repeat the process on these refined parallelograms. 
We iterate this until we have covered a volume fraction of size at least v ∈ (0, 1). At each step, we have a list _v(T, M)[1] of “good" parallelograms (green in Figure <ref>), the discarded parallelograms which do not need to be stored any longer (black in Figure <ref>), and a remaining list of parallelograms to be further refined (white in Figure <ref>). As our shapes are convex polytopes, efficient intersection algorithms can be implemented.
The remainder of the domain, which is not covered by this procedure is split into
triangles, which are inserted into _v(T, M)[2]. 









§ NUMERICAL IMPLEMENTATION OF THE CONVEX INTEGRATION SCHEMES FOR THE HEXAGONAL-TO-RHOMBIC PHASE TRANSFORMATION



In this section, we present some of the output of Algorithm <ref> in both of the variants described above. Here we use a colour coding based on the CMYK colour model, where the colours cyan, magenta and yellow (we do not use key) correspond to one of the wells respectively. More precisely, the well e^(1) corresponds to cyan, the well e^(2) to magenta, and e^(3) to yellow.
A matrix which is included in the interior of the convex hull is correspondingly
depicted as a convex combination of these colourings (c.f. Figure
<ref>).

We remark that in commercial printing the CMYK colour model is used to blend colours as mixtures of the
primary colours cyan, magenta, yellow (and black). In particular,
for our purposes this is very convenient, as it allows us to work with “barycentric coordinates" based on cyan, magenta and yellow. 
Moreover, too fine
structures in our construction, which are hard to see in the picture due to limited
resolution and rasterisation, are then “homogenised" automatically in a way which agrees with the colours of the averaged matrix values. 



 §.§ Output of the implementation of Algorithm <ref>


We compare the output of the two variants of Algorithm <ref> for three different scenarios. These are chosen such that the overall boundary data M ∈({e^(1), e^(2),e^(3)}) are such that either (c.f. Figure <ref>)


  (i) M is close to the barycenter 1/3e^(1)+1/3e^(2)+1/3e^(3)=0 of the equilateral triangle spanned by e^(1), e^(2), e^(3).

  (ii) M is close to one of the wells.

  (iii) M is close to the boundary of (e^(1),e^(2), e^(3)) and essentially in between two wells.

We remark that, due to the Skew(2) symmetry, our construction in Algorithm 3.3
essentially only depends on the symmetric part e(M) of the matrix M (the
skew part only amounts to orthogonal translations of the matrices). We discuss the model cases (i)-(iii) separately in the sequel.









  §.§.§ An example of case (i)

As a first case, we discuss boundary data which are chosen to lie very close to the barycenter of the equilateral triangle spanned by the matrices e^(1), e^(2), e^(3). 





In the implementations of the two algorithms from <cit.> and <cit.> we see major differences, which are illustrated in Figure <ref>. As discussed in Section <ref> we observe the following differences:


  * finiteness: The construction of <cit.> is finite, which is
  manifested in the fact that already after three iteration steps there are
  large patches (approximately 20 percent), in which the final colouring is attained (in Figure <ref> for instance there are large yellow patches). In contrast to this, the algorithm from <cit.> is countably infinite and although in the fine structure one already sees colours close to the wells, these are not exactly in the wells in general. In particular this yields a very fine-scale structure, which  “homogenises" in the “eye norm". 
  

  * uniformity: In the implementation of the algorithm <cit.> the replacement constructions always contain patches which are “pushed away" into the interior (c.f. Figure <ref>). In the implementation of the algorithm from <cit.>, we in contrast see that the replacement constructions use constructions whose colours are increasingly close to the wells (c.f. Figure <ref>, right panel).

  * normal directions: While the normals which are used in the algorithm from <cit.> may vary in a continuum, the ones in the algorithm from <cit.> are discrete. Although this is hard to see by the “eye norm", it is manifested in the replacement constructions shown in Figure <ref>.  

Both constructions are of self-similar, fractal structure (however leading to quite different overall structures).







  §.§.§ An example of case (ii)




In addition to the differences which had already been discussed in the case (i),
the most striking point which can is reflected in the implementations of the case
(ii) is the “inertia" which is underlying the algorithm from <cit.>
before the well is changed. This can best be understood by recalling the
replacement algorithm from Section <ref>, which is also illustrated
in Figure <ref>. Here the original boundary data are replaced
by five new matrices, of which some are “pushed further into the interior" of
the convex hull of the strains e^(1), e^(2), e^(3). This implies that
the closest well eventually changes. However, since the push-out is achieved
with a fixed factor, this is slower, the closer the original boundary datum had
been to one of the wells. In particular, if one starts with a matrix very close to one of the wells, it might take a large number of iteration steps before any other well is used (as the algorithm will have reached very small scales by then, it is also hard to see this, even if this happens).
   
In contrast to this, the algorithm from <cit.> uses values close to two of the wells in each iteration step: Hence even if one starts close to a single well, the other wells will emerge within the next two iteration steps (c.f. the second panel of Figure <ref>).


Again both of the constructions (the one from <cit.> and the one from <cit.>) already display fractal behaviour.




  §.§.§ An example of case (iii)


Last but not least, we discuss the setting in which the boundary data are very close to the rank-one line between two wells. Here both constructions mainly consist of shades of cyan and magenta. The angle between the closest well and the boundary is very flat. Hence, we expect that the construction from <cit.> (left panel in Figure <ref>) is forced to use very fine scales depending on the boundary distance. In particular, the regularity in the finely twinned areas should be bad, eventually yielding regularities which are significantly rougher than the ones from cases (i), (ii) (c.f. Figures <ref>, <ref>). In contrast, the construction from <cit.> (right panel in Figure <ref>) is not significantly more complex than the constructions from Figures <ref>, <ref>.







§ DISCUSSION AND SUMMARY



We conclude the article by comparing our findings with some related models. Here we in particular return to their relation to scaling laws and comment on their relation to the Ball-Planes model <cit.>.



 §.§ Relation to scaling limits

As already discussed in the introduction, scaling laws for energies such as (<ref>) play an important role in the mathematical analysis of shape-memory materials. It is not known whether energies which include surface terms and which hence penalise too high oscillations are compatible with convex integration solutions. In our Theorem <ref> we could however relate the presence of lower bounds to regularity thresholds for convex integration solutions. In this sense, models involving surface energies could be interpreted as a selection mechanism.



 §.§ Probabilistic models on nucleation


An interesting, probabilistic model for nucleation displaying “wild" structures has been proposed by the groups around J. Ball and A. Planes and E. Vives <cit.>. This model takes the dynamic process of nucleation into account and assumes that in the nucleation process the different phases are nucleated according to the following rules (there are different subschemes of this): Let Ω⊂^2 be a bounded Lipschitz domain representing the sample. The model is then based on the following algorithm:


  * Choose a random point x∈Ω.

  * According to a prescribed probability distribution, choose from a set of finitely many possible orientations, which also each have an associated colour. 
Physically, these essentially model the rank-one lines between the wells and hence also the variants of martensite associated with the given orientation. 

  * “Nucleate" the chosen variant of martensite with the prescribed orientation at the point x∈Ω. This is achieved by inserting a line segment (with or without prescribed width) with a prescribed associated colour (for each orientation) through the point x∈Ω of the prescribed orientation. This line segment is maximally extended through the point x∈Ω with the chosen orientation until it hits the boundary of Ω or another already present line, in which case it terminates.

It is assumed that the transformation process is irreversible, i.e. an already inserted line cannot be erased in the explained algorithm.

This seemingly “simple" process is used to explain the lack of a fixed length scale, the experimentally observed fractal behaviour and the intermittency in the nucleation process of martensite <cit.>. A certain “universality" of the distribution of the lengths of line segments is identified <cit.>. However, compatibility considerations do not enter the model.


Our geometrically more complicated setting is not viewed as a dynamic process (although it might be possible to give a dynamic meaning to the iteration parameter k ∈ in the convex integration algorithm). However, in contrast to the Ball-Planes-Vives model, we include compatibility in our considerations. The lack of a fixed length scale and the fractal behaviour are also reflected in our convex integration solutions. We also expect fat tailed distributions for the Fourier transforms of our length scale distribution functions, c.f. <cit.>. Mathematically this and the universality can be viewed as an (ir)regularity result. It would be interesting to see whether there is a closer relation between the two models and whether the two models could be combined.

We note that both models produce highly non-unique solutions. Yet, it might be possible to identify statistical properties characterising these (possibly uniquely).



 §.§ Experiments

It is not known whether convex integration solutions really appear in nature. Although, there are indications that in for instance the cubic-to-monoclinic phase transformation, the nucleation of martensite in austenite creates quite wild microstructures <cit.>, it would be interesting to further compare the experimental results with numerical/analytical predictions.



 §.§ Summary

In the present note we have presented two main results: 

On the one hand, we have related the maximal possible regularity of convex integration solutions to the existence of lower bounds in the associated scaling laws. This could present a possibility of deducing bounds on the maximal regularity of convex integration solutions (at least in “sufficiently simple" models for which one can hope to obtain scaling laws).


On the other hand, we have presented a first implementation of convex integration solutions in the context of the hexagonal-to-rhombic phase transformation, which we view as a model problem. We have compared the numerical outputs of the two algorithms from <cit.> and <cit.>.

Of course this is only a first modest step towards an improved understanding of convex integration solutions, which we seek to further study in future work. Although experimental results on nucleation suggest quite wild structures in certain phase transformations <cit.>, we have no direct comparison between these solutions and numerical predictions yet (which in particular is expected to be difficult to obtain due to the intrinsic non-uniqueness of convex integration solutions; a possible remedy could be the statistical analysis of these solutions). We seek to pursue this comparison in future work. In particular, we emphasise that the question whether convex integration solutions arise in nature in the context of phase transformations in shape-memory alloys remains an exciting open problem.





§ APPENDIX



In this last section, for self-containedness, we recall a possible proof of the fractional Poincaré inequality (for fractional Besov type spaces), which is used in different forms in Section <ref>. We start with the L^2 based version:



Let s∈ (0,1) and n≥ 1.
Let Ω⊂^n be an open, bounded C^1,1 domain.
Let u∈ H^s(^n) with  (u) ⊂Ω and denote Ω_δ={x∈Ω: (x,∂Ω)≥δ} for δ>0.
Then, there exists C=C(s,Ω,n)>1 such that

    u_L^2(Ω_δ∖Ω_2δ)≤ C δ^su_H^s(^n).





We use the Caffarelli-Silvestre extension u of u <cit.>, which is the H^1(^n+1_+, x_n+1^1-2s) solution to the equation

    ∇· x_n+1^1-2s∇u    = 0 ^n+1_+,
    u    = u ^n ×{0}.

For x_1 ∈Ω_δ∖Ω_2δ,
x_2=p(x_1) ∈^n ∖Ω with |x_1-x_2|≤ 4 δ and δ̃∈ [δ, 2δ], the fundamental theorem and Hölder's inequality yield

    |u(x_1)|^2    = |u(x_1,0)|^2
    ≤ 2|u(x_1,δ̃)|^2 + 
    2 ( ∫_0^δ̃|_t u(x_1,t)| dt )^2
       ≤ 2|u(x_1,δ̃)|^2 + 
    C δ̃^2s∫_0^δ̃ t^1-2s|_t u(x_1,t)|^2 dt
       ≤ 2 |u(p(x_1),δ̃)|^2 + C δ̃^2 ∫_0^1|∇' u(τ x_1 + (1-τ)p(x_1), δ̃)|^2 dτ
          + C δ̃^2s∫_0^δ̃ t^1-2s|_t u(x_1,t)|^2 dt
       ≤  2δ̃^2 ∫_0^1|∇' u(τ x_1 + (1-τ)p(x_1), δ̃)|^2 dτ + C δ̃^2s∫_0^δ̃ t^1-2s|_t u(x_1,t)|^2 dt
         
    + C δ̃^2s∫_0^δ̃ t^1-2s|_t u(p(x_1),t)|^2 dt.

Here the constant C>1 changes from line to line and depends on s but not on δ. In deriving the above estimate, we applied the fundamental theorem thrice: First in the normal direction (where we then used Hölder's inequality to insert the weight t^1-2s), then in the tangential directions and finally once more in the normal direction (where we used the vanishing Dirichlet data for u(p(x_1))=u(p(x_1),0)).
Estimating the integrals involving the normal derivative by using δ <δ̃≤ 2δ, we thus infer

    |u(x_1)|^2 
       ≤  2δ̃^2 ∫_0^1|∇' u(τ x_1 + (1-τ)p(x_1), δ̃)|^2 dτ + C δ^2s∫_0^2δ t^1-2s|_t u(x_1,t)|^2 dt
         
    + C δ^2s∫_0^2δ t^1-2s|_t u(p(x_1),t)|^2 dt.

Averaging over δ̃∈ [δ, 2δ] yields

    |u(x_1)|^2 
       ≤ 2δ^-1∫_δ^2δδ̃^2 ∫_0^1|∇' u(τ x_1 + (1-τ)p(x_1), δ̃)|^2 dτ d δ̃
    + C δ^2s∫_0^2δ t^1-2s|_t u(x_1,t)|^2 dt
         
    + C δ^2s∫_0^2δ t^1-2s|_t u(p(x_1),t)|^2 dt
        =  2δ^-1∫_δ^2δδ̃^1+2sδ̃^1-2s∫_0^1|∇' u(τ x_1 + (1-τ)p(x_1), δ̃)|^2 dτ d δ̃
    + C δ^2s∫_0^2δ t^1-2s|_t u(x_1,t)|^2 dt
         
    + C δ^2s∫_0^2δ t^1-2s|_t u(p(x_1),t)|^2 dt
       ≤  C δ^-1δ^1+2s∫_δ^2δδ̃^1-2s∫_0^1|∇' u(τ x_1 + (1-τ)p(x_1), δ̃)|^2 dτ d δ̃
    + C δ^2s∫_0^2δ t^1-2s|_t u(x_1,t)|^2 dt
         
    + C δ^2s∫_0^2δ t^1-2s|_t u(p(x_1),t)|^2 dt
       ≤  C δ^2s∫_0^2δ∫_0^1 t^1-2s |∇' u(τ x_1 + (1-τ)p(x_1), t)|^2 dτ dt
    + C δ^2s∫_0^2δ t^1-2s|_t u(x_1,t)|^2 dt
         
    + C δ^2s∫_0^2δ t^1-2s|_t u(p(x_1),t)|^2 dt.

Finally, integrating over x_1 ∈Ω_δ∖Ω_2δ leads to

    u_L^2(Ω_δ∖Ω_2δ)^2
        ≤ C δ^2s∫_0^1 ∫_0^2δ t^1-2s∇' u(τ· + (1-τ)p(·),t)_L^2(Ω_δ∖Ω_2δ)^2 dt dτ
          + 2C δ^2st^1-2s/2_t u_L^2((Ω_δ∖Ω_2δ) × [0,2δ])^2
       ≤ C δ^2st^1-2s/2∇u_L^2(Ω× [0,2δ])^2
       ≤ C δ^2st^1-2s/2∇u_L^2(^n+1_+)^2.

Hence, using that (c.f. <cit.>)

    t^1-2s/2∇u_L^2(^n+1_+)^2
    ≤ C u_Ḣ^s(^n),

we obtain the desired estimate (<ref>).


Similarly, we have the analogue of this in Sobolev and Besov spaces:



Let s∈ (0,1), p∈ (1,∞) and n≥ 1.
Let Ω⊂^n be an open, bounded C^1,1 domain.
Let u∈ W^s,p(^n) with  (u) ⊂Ω and denote Ω_δ={x∈Ω: (x,∂Ω)≥δ} for δ>0.
Then, there exists C=C(s,p,Ω,n)>1 such that

    u_L^p(Ω_δ∖Ω_2δ)≤ C δ^su_W^s,p(^n).





The argument follows essentially as above. However for the estimate in the normal direction, we use

    |u(x_1,0)|^p
       ≤ C |u(x_1,δ)|^p + C(∫_0^δ_t u dt)^p
       ≤ C|u(x_1,δ)|^p + C∫_0^δ |t^1-1/p-s _t u|^p dt (∫_0^δ t^-1+sp/p-1dt )^p-1
       ≤ C|u(x_1,δ)|^p +  Cδ^sp∫_0^δ t^p-1-s p|_t u|^p dt .

Combining this with estimates in the tangential directions similarly as in the proof of Lemma <ref>, we obtain 

    u_L^p(Ω_δ∖Ω_2δ)   ≤ C δ^st^1-1/p-s∇u_L^p(^n+1_+)≤ C δ^su_W^s,p(^n).

Here we used the trace characterisation of W^s,p(^n), c.f. <cit.>, where the authors rely on the characterisation from <cit.>.


A similar argument using the characterisation of Besov spaces from 
<cit.> (where the authors again rely on the characterisation from <cit.>) also yields a similar Poincaré estimate in Besov spaces.

alpha


