Approximations of energy minimization in cell-induced phase transitions of fibrous biomaterials: $\Gamma$-convergence analysis

TL;DR

This paper develops a new mathematical framework for approximating energy minimization in fibrous biomaterials undergoing phase transitions, ensuring convergence of numerical solutions to the true continuous model using $ ext{Γ}$-convergence and Orlicz spaces.

Contribution

It introduces a novel approach to $ ext{Γ}$-convergence for discontinuous Galerkin methods applied to complex energy minimization problems in nonlinear elasticity.

Findings

Discrete minimizers converge to continuous minimizers as discretization refines

A new $ ext{Γ}$-convergence approach for lower regularity spaces

Development of an Orlicz space framework for penalized interpenetration terms

Abstract

We consider a model of energy minimization arising in the study of the mechanical behavior caused by cell contraction within a fibrous biological medium. The macroscopic model is based on the theory of non rank-one convex nonlinear elasticity for phase transitions. We study appropriate numerical approximations based on the discontinuous Galerkin treatment of higher gradients and used succesfully in numerical simulations of experiments. We show that the discrete minimizers converge in the limit to minimizers of the continuous problem. This is achieved by employing the theory of $\Gamma$-convergence of the approximate energy functionals to the continuous model when the discretization parameter tends to zero. The analysis is involved due to the structure of numerical approximations which are defined in spaces with lower regularity than the space where the minimizers of the continuous variational problem are sought. This fact leads to the development of a new approach to $\Gamma$-convergence, appropriate for discontinuous finite element discretizations, which can be applied to quite general energy minimization problems. Furthermore, the adoption of exponential terms penalising the interpenetration of matter requires a new framework based on Orlicz spaces for discontinuous Galerkin methods which is developed in this paper as well.

Equations

$$\Psi [\mathbf{u}]={\int }_{\Omega }\left[W(\nabla \mathbf{u}(\mathbf{x}))+\Phi (\nabla \mathbf{u}(\mathbf{x}))+\frac{{\epsilon }^2}{2}|\nabla \nabla \mathbf{u}(\mathbf{x}){|}^2\right]d\mathbf{x},$$

$${\Psi }_h[{\mathbf{u}}_h]=\underset{{\mathbf{w}}_h\in {\mathbb{A}}_h^q(\Omega )}{\inf }{\Psi }_h[{\mathbf{w}}_h].$$

$$c_0\left(|1+\nabla \mathbf{u}{|}^2-c_1\right)\le W(\nabla \mathbf{u})\le c_2\left(|1+\nabla \mathbf{u}{|}^m+c_3\right),$$

$$\displaystyle\mathbb{A}(\Omega)=\{\bm{u}\in H^{2}(\Omega)^{2}:{\left.\kern-1.2pt\bm{u}\vphantom{\big{|}}\right|_{\partial\Omega}}={\left.\kern-1.2pt\bm{g}\vphantom{\big{|}}\right|_{\partial\Omega}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\left.\kern-1.2pt\nabla\bm{u}\vphantom{\big{|}}\right|_{\partial\Omega}}={\left.\kern-1.2pt\nabla\bm{g}\vphantom{\big{|}}\right|_{\partial\Omega}}}\}.$$

$$\inf \{\Psi [\mathbf{u}]:\mathbf{u}\in \mathbb{A}(\Omega )\}.$$

$${\rho }_K\ge \frac{h_K}{c},\text{for all }K\in T_h,$$

$$C{h}_K\le h_e\le c{h}_K,\text{for }e\in \partial K\text{ and all }K\in T_h.$$

$$\displaystyle V^{q}_{h}(\Omega)=\{v\in C^{0}(\Omega):{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{K}}\in\mathbb{P}_{q}(K),K\in T_{h}\},\quad q\in\mathbb{N}.$$

$$\displaystyle\tilde{V}^{k}_{h}(\Omega)=\{v\in L^{2}(\Omega):{\left.\kern-1.2ptv\vphantom{\big{|}}\right|_{K}}\in\mathbb{P}_{k}(K),K\in T_{h}\},\quad k\in\mathbb{N}.$$

$$T(E_h):={\Pi }_{e\in E_h}{L}^2(e),$$

$$\{\{\cdot \}\}:$$

$$\displaystyle\{\kern-3.125pt\{\bm{w}\}\kern-3.125pt\}:=\left\{\begin{array}[]{ll}\frac{1}{2}({\left.\kern-1.2pt\bm{w}\vphantom{\big{|}}\right|_{K_{e^{+}}}}+{\left.\kern-1.2pt\bm{w}\vphantom{\big{|}}\right|_{K_{e^{-}}}}),&\quad\text{for }e\in E_{h}^{i}\\ \bm{w},&\quad\text{for }e\in E_{h}^{b},\end{array}\right.$$

$$[[\cdot ]]:$$

$$\displaystyle\llbracket\bm{v}\rrbracket:={\left.\kern-1.2pt\bm{v}\vphantom{\big{|}}\right|_{K_{e^{+}}}}-{\left.\kern-1.2pt\bm{v}\vphantom{\big{|}}\right|_{K_{e^{-}}}},$$

$$[[\cdot ]]:$$

$$\displaystyle\llbracket\bm{v}\otimes\bm{n}_{e}\rrbracket:={\left.\kern-1.2pt\bm{v}\vphantom{\big{|}}\right|_{K_{e^{+}}}}\otimes\bm{n}_{e^{+}}+{\left.\kern-1.2pt\bm{v}\vphantom{\big{|}}\right|_{K_{e^{-}}}}\otimes\bm{n}_{e^{-}},$$

$${\Psi }_h[{\mathbf{u}}_h]$$

$$\displaystyle+\varepsilon^{2}\Bigg{(}\frac{1}{2}\sum_{K\in T_{h}}\int_{K}|\nabla\nabla\bm{u}_{h}|^{2}-\sum_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e\in E_{h}}}\Big{[}\int_{e}\{\kern-3.125pt\{\nabla\nabla\bm{u}_{h}\}\kern-3.125pt\}\cdot\llbracket\nabla\bm{u}_{h}\otimes\bm{n}_{e}\rrbracket+\frac{\alpha}{h_{e}}\int_{e}|\llbracket\nabla\bm{u}_{h}\rrbracket|^{2}\Big{]}\Bigg{)}$$

$$={\int }_{\Omega }[W(\nabla {\mathbf{u}}_h)+\Phi (\nabla {\mathbf{u}}_h)]+{\epsilon }^2{\Psi }_h^{ho}[{\mathbf{u}}_h]$$

$$[[\nabla {\mathbf{u}}_h]]:=\nabla {\mathbf{u}}_h^{-}-\nabla I_h^q\mathbf{g},\text{ and }[[\nabla {\mathbf{u}}_h\otimes {\mathbf{n}}_e]]:=(\nabla {\mathbf{u}}_h^{-}-\nabla I_h^q\mathbf{g})\otimes \mathbf{n},\text{ for }e\in E_h^b=E_h\cap \partial \Omega ,$$

$$\displaystyle\mathbb{A}_{h}^{q}(\Omega)=\{\bm{u}_{h}\in V^{q}_{h}(\Omega)^{2}:{\left.\kern-1.2pt\bm{u}_{h}\vphantom{\big{|}}\right|_{\partial\Omega}}={\left.\kern-1.2pt\bm{g}_{h}\vphantom{\big{|}}\right|_{\partial\Omega}}\},$$

$$\inf \{{\Psi }_h[{\mathbf{u}}_h]:{\mathbf{u}}_h\in {\mathbb{A}}_h^q(\Omega )\}.$$

$$\frac{k}{2}{\int }_{{\Gamma }_i}|\mathbf{u}-\mathbf{g}{|}^{2}ds,\text{ where }{\Gamma }_i\subset \partial \Omega \text{ and }|{\Gamma }_i|>0,$$

$$|\hat{w}{|}_{W^{p,r}(\hat{K})}=h_K^{p-n/r}|w{|}_{W^{p,r}(K)}.$$

$$\Vert w{\Vert }_{L^p(\partial \Omega )}\le C\Vert w{\Vert }_{L^p(\Omega )}^{1-1/p}\Vert w{\Vert }_{W^{1,p}(\Omega )}^{1/p}.$$

$$\Vert u_h{\Vert }_{L^2(e)}^2\le \frac{c_q}{h_e}\Vert u_h{\Vert }_K^2,\forall u_h\in {\mathbb{P}}_q(K),\forall K\in T_h.$$

$$I_h^q:H^s(\Omega {)}^2\to V_h^q(\Omega {)}^2,s\ge 2.$$

$$|\mathbf{u}-I_h^q\mathbf{u}{|}_{H^m(K)}\le c{h}_K^{s-m}|\mathbf{u}{|}_{H^s(K)},0\le m\le s$$

$$|\mathbf{u}-I_h^q\mathbf{u}{|}_{H^m(e)}\le c{h}_K^{s-m-1/2}|\mathbf{u}{|}_{H^s(K)},$$

$$|\mathbf{w}{|}_{H^2(\Omega ,T_h)}^2:=\sum _{K\in T_h}{\int }_K|\nabla \nabla \mathbf{w}{|}^2+\sum _{e\in E_h}\frac{1}{h_e}{\int }_e|[[\nabla \mathbf{w}]]{|}^2.$$

$$\displaystyle\lVert\nabla\bm{w}\rVert_{L^{r}(\Omega)}^{2}\leq c_{\theta_{min}}\Big{(}|\bm{w}|_{H^{2}(\Omega,T_{h})}^{2}+|\widehat{\nabla\bm{w}}|^{2}\Big{)}$$

$$\displaystyle\widehat{\nabla\bm{w}}=\left\{\begin{array}[]{ll}\frac{1}{|\Omega|}\int_{\Omega}\nabla\bm{w}\\ \frac{1}{|\partial\Omega|}\int_{\partial\Omega}\nabla\bm{w}\\ \frac{1}{|\partial\Gamma_{i}|}\int_{\partial\Gamma_{i}}\nabla\bm{w},\quad\Gamma_{i}\subset\partial\Omega\text{ and }|\Gamma_{i}|>0.\end{array}\right.$$

$$\displaystyle{\left.\kern-1.2pt(\nabla_{h}w_{i})\vphantom{\big{|}}\right|_{K}}:=\nabla({\left.\kern-1.2ptw_{i}\vphantom{\big{|}}\right|_{K}}),\quad i=1,...,m,\forall K\in T_{h}.$$

Full text

1 Introduction

Increasingly sophisticated mathematical techniques are needed in order to describe biological phenomena. Biological cells severely deform the fibrous collagen extracellular material (ECM) around them in remarkable ways, causing the formation of complex microstructures of highly localized deformation [33, 52]. In [32], this phenomenon, known as mechanical remodelling of the ECM, is modelled at a macroscopic level, as a phase transition in a nonlinear elastic material with a multi-well energy. The nonconvexity is due to buckling instability of ECM fibers at the microsopic level. Subsequently, computational predictions based on the mathematical model, combined with targeted experiments, lead for the first time to understanding the mechanisms of the observed ECM remodelling. The latter facilitates intercellular communication through the formation of tethers, regions where a densification phase transition takes place, joining distant contracting cells. The associated variational problem involves a non rank-one convex strain-energy function, regularized by a higher gradient term.

Here, our objective is to mathematically justify the above procedure, by showing that appropriate numerical approximations, based on the discontinuous Galerkin treatment of higher gradients, and used very successfully in computational experiments, indeed converge in the limit to minimizers of the continuous problem. This is done by employing the theory of $\Gamma$-convergence of the approximate energy functionals to the continuous model when the discretization parameter tends to zero. This is a rather involved task due to the structure of numerical approximations, which are defined in spaces with lower regularity than the space where the minimizers of the continuous variational problem are sought.

Our work has a number of methodological advances which go beyond the scope of the particular application. In fact, within the field of nonlinear PDEs with possibly singular solutions, calculus of variations and energy minimization has received a lot of attention from the analysis point of view. Although quite interesting and challenging, the numerical analysis of these problems is much less developed, a variety of approaches are discussed in, e.g., [6, 7, 20, 41, 46] and their references. Key issues, related to the design of computational algorithms, their analysis, selection criteria among possible solutions and minimizers remain largely unexplored. Our work contributes to scheme design and analysis of such problems. The problem considered is a typical nonlinear energy minimization problem which admits solutions exhibiting phase transitions. We show that a good scheme design strategy is to regularise at the discrete level the continuum energies by higher order gradients. The regularisation at the discrete level considered is of similar nature to the artificial diffusion in conservation laws and it appears that it enjoys remarkable properties in terms of computational robustness and analytical consistency. Compared to other approaches, such as relaxation, [7], the regularisation by higher gradients has certain distinct desirable characteristics. Computing with the relaxed energy, assuming that we can obtain it, smears over the multiphase mixture (microstructure) or oscillations, which in certain applications (including ours) are very physical as they are observed in experiments. Furthermore, when microstructures appear, quite often a mathematical quantity of interest is an underlying parametrized Young measure. The regularisation parameter $\varepsilon$ sets an upper bound to the frequency of the oscillations and thus allows micro-phenomena to be present in a computationally accessible way; thus such oscillations (not necessarily extremely fine) may provide approximations to the underlying Young measure. The discretization of higher gradients is done using the discontinuous Galerkin framework, thus retaining the regularity of only $C^{0}$ conforming elements. This natural approach on the other hand poses new challenges in the analysis of schemes. New ideas are needed in the $\Gamma$-convergence analysis due to the lack of conformity at the higher-gradient spaces. The adoption of exponential terms penalising the interpenetration of matter requires a new framework based on Orlicz spaces for discontinuous Galerkin methods, developed in this paper as well.

Soon after a preprint version of this work was published [31], a related preprint [10] appeared, which also treats $\Gamma$-convergence of the discetized energy from a nonconvex mechanics problem, using totally discontinuous finite elements (these two works were developed independently of each other). Apart from this similarity, there are substantial differences in both the model and the discretization that set these papers apart.

The model. We consider the problem of minimizing the total potential energy

[TABLE]

where the displacement $\bm{u}\in H^{2}(\Omega)^{2}$ and satisfies some appropriate boundary conditions, $W+\Phi$ is the strain energy function, $\Phi$ is a function that penalizes the interpenetration of matter and is allowed to grow faster than $W$ as the volume ratio approaches zero, and $\varepsilon>0$ is a fixed real parameter (higher gradient coefficient). The energy involves a non rank-one convex strain-energy function, regularized by a higher order term. The penalty term $\Phi$ is important since, although it permits the appearance of phase transitions, it prevents interpenetration of matter from taking place. The strain energy function models the bulk response of the collagen ECM, while the higher order term represents a length scale for the thickness of phase transition layers and the emerging two-phase microstructures. Specifically, the strain energy function models the mechanical response of the ECM which is a collagen material in the form of a random network of fibers at the microscopic level. Biological cells, such as fibroblasts, are embedded in the ECM. They are attached onto the ECM fibers through proteins known as focal adhesions. Through these molecules, cells can detect mechanical alterations to their microenvironment and can deform the surrounding ECM. Cells typically deform the ECM by actively contracting. These tractions are observed to create distinct spatial patterns of localized, severe densification in the ECM between cells, forming a tether connecting them, and around the periphery of the cell in the form of hair-like microstructures, [33, 52, 45, 32]. These microstructures and the associated strain oscillations are strongly reminiscent of fine phase mixtures in the theory of nonlinear elasticity for phase transitions [5, 4]. This similarity was explained in [32] where a model for the strain energy density was obtained via multiscale modelling from the stress-strain behavior of single fibers comprising the bulk ECM material, as summarized below in section 8. The central physical observation related to instability and phase change [38] is that individual collagen fibers can sustain tension, but buckle and collapse under compression. At the microscopic level, this is accounted for as an effective softening behavior in uniaxial compression of fibers in our model. The macroscopic behavior is obtained by averaging over a uniform distribution of fiber orientations. This results in a strain energy function that loses rank-one convexity, and is essentially equivalent to a multi-well potential corresponding to a densification phase transition at the continuum scale. The phases correspond to low- and high-density states and their simulation via energy minimization leads to a remarkable agreement with experimental observations of densification microstructures (an example of an experiment is shown in Figure 5a; the corresponding simulation is shown in Figure 5b). These microstructures are composed of tethers, or relatively straight bands joining different cells, and hairs, thinner multiple bands emanating radially from each cell and tapering off into the ECM. The ECM density within both tethers and hairs can be 3-5 times larger than outside them.

Failure of rank-one convexity implies a loss of ellipticity of the Euler-Lagrange PDEs of the corresponding energy functional. For a wide class of similar problems it is known [5, 4] that there exist oscillatory minimizing sequences with finer and finer microstructures involving increasing numbers of strain jumps. Similar behavior is observed in our model; the numerical approximations—obtained by mesh refinement—of terms of increasing fineness in the minimizing sequences, involve more and thinner hairs (see Figures 1a-1d). They also bear strong similarity with experiments [32].

The higher gradient term in (1.1) regularises the corresponding total potential energy, keeping the aforementioned minimizing sequences from having arbitrarily fine structure (see Figure 3 for numerical examples).

To the best of our knowledge, the deformations observed in [32] and in the present study are the first examples of minimizing sequences in a multi-well compressible isotropic material.

Approximations and results. The approximation of minimizers of $(\ref{equ:total_potential})$ is quite subtle. A straightforward approach would be to seek approximate minimizers in the space of conforming finite elements, i.e. of discrete function spaces which are finite dimensional subspaces of $H^{2}(\Omega)^{2}$. Such spaces are based on elements which require $C^{1}$ regularity across element interfaces, e.g. Argyris elements [14]. However, the conformity in regularity has a very high computational cost under the minimization process and in addition results in much more complicated algorithms as far as the implementation is concerned. Our choice is to use the framework of the discontinuous Galerkin method. In effect this weakens the regularity of the approximating spaces, and counterbalances the resulting nonconformity, by amending appropriately the discrete energy functional. Motivated by the analysis in [43], we introduce an approximate energy which is compatible with $C^{0}$ finite element spaces, and thus requires only $H^{1}$ regularity. Corresponding finite element methods, known as $C^{0}$-interior penalty methods, have been introduced previously for the approximation of the biharmonic equation in [16, 27]; see also [2] for fully discontinuous finite element methods.

Here we study the convergence of discrete absolute minimizers. Specifically, let $(\bm{u}_{h})$ be a sequence of absolute minimizers for the discretized energy functional $\Psi_{h}$, namely

[TABLE]

Equation (1.2) indicates that, for a fixed $h$, $\bm{u}_{h}$ is an absolute minimizer of $\Psi_{h}$. Therefore, as $h\rightarrow 0$, it is natural to ask whether $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$, where $\bm{u}$ is an absolute minimizer of the continuous problem (1.1). Note that $\bm{u}_{h}\in H^{1}(\Omega)^{2}$ and $\bm{u}\in H^{2}(\Omega)^{2}$. To answer this question we assume first that the penalty function $\Phi$ has polynomial growth. Then the convergence result is given in Theorem 6.1, where we have employed the theory of $\Gamma-$convergence and discrete compactness results. The analysis is rather involved due to the lack of regularity of the approximate spaces. A $\Gamma$-convergence result for discrete surface functionals involving high gradients using conforming finite element spaces can be found in [6]. Assuming that the penalty term $\Phi$ has exponential growth, extra embedding results are needed to show that $\Psi_{h}$ $\Gamma-$converges to $\Psi$. For this purpose, it is crucial to use an adaptation of Trudinger’s embedding theorem for Orlicz spaces [53], to the piecewise polynomial spaces admitting discontinuities in the gradients, Theorem 7.1. The analysis in this case is carried out in Section 7. It is to be noted that in our approach we are interested in the limit $h\to 0$ for any fixed $\varepsilon.$ The tools to address the very interesting case $\varepsilon\to 0$ for rather general non rank-one convex functionals as $W$ considered herein, are currently not available. In addition, we remark that the case where $\Phi(\nabla\bm{u})\rightarrow\infty$ as $\det(\bm{1}+\nabla\bm{u})\to 0$ is currently beyond our reach.

The present work addresses key technical issues related to the analysis of approximations of energy minimization problems involving higher gradients. A main technical obstacle in proving $\Gamma$-convergence is the fact that the approximating discrete energy functionals are defined in spaces of lower regularity compared to the limiting functional as a result of the discontinuous Galerkin formulation. Notice here that from a computational perspective the use of $C^{0}$ elements permits direct comparisons with the approximations obtained even in the limit case $\varepsilon=0.$ In our analysis we use certain recovery operators for the higher gradients, well known in the analysis of discontinuous Galerkin methods. As such, previous results from [19, 25] are useful. Notice that compared to the method in [19] where recovery operators were used in the definition of the discretization method as well, our method leads to the natural discontinuous Galerkin method, in the sense that the part corresponding to the higher gradients in the energy introduced herein has as first variation the bi-linear form used in [16]. As mentioned, to treat exponential penalty terms in the energy functional we need to develop an appropriate discontinuous Galerkin framework in Orlicz spaces (Section 7). To this end, results of [35] for averaging operators have proven useful.

This paper is organized as follows. In section 2 we discuss some properties of the continuum problem, a lower bound is proved and the minimization problem is stated. In section 3 the necessary notation, some standard finite elements results, the discrete total potential energy $\Psi_{h}$ and lifting operators used in the next sections are introduced. Equi-coercivity, the $\liminf$ and the $\limsup$ inequalities are proved for the discrete energy functional in section 5, which imply $\Psi_{h}\xrightarrow{\enskip\Gamma\enskip}\Psi$. From the $\Gamma$-convergence result and a discrete compactness property we deduce the convergence of the discrete absolute minimizers, section 6. In section 7 the same convergence result is established when the penalty function has exponential growth. In this section we derive key embeddings of broken polynomial spaces into an appropriate Orlicz space. We conclude with section 8 illustrating some computational results which demonstrate the robustness of the approximating scheme and the model when both the mesh discretization parameter and $\varepsilon$ vary.

2 The Continuum Problem

We assume the following bounds for the terms in equation (1.1)

[TABLE]

for some $m\geq 2$ and positive constants $c_{0},c_{1},c_{2},c_{3}$. Also, we will assume that the penalty term satisfies the conditions

[TABLE]

again for some $m_{0}\geq 1$ and positive constants $C_{0},C_{1},C_{2},C_{3}$.

To define a minimization problem we should declare appropriate boundary conditions. We specify a globally injective, orientation preserving $\bm{g}\in{H^{3}(\Omega)^{2}}$. We encode boundary conditions in the following set:

[TABLE]

Notice that the boundary conditions were chosen for analytical convenience only. A variety of other boundary conditions can be treated with appropriate modifications in the finite element analysis, see Remark 3.1. Now, the minimization problem can be defined as:

[TABLE]

To ensure that the total potential energy has a minimizer $\bm{u}\in\mathbb{A}(\Omega)$, one can prove that $\Psi[\cdot]$ is coercive and lower semicontinuous. The former can be derived from the properties of the strain energy function, i.e. (2.1), and the Poincaré inequality. One way to prove the latter is to show convergence of the lower order term using $Vitali^{\prime}s$ Theorem and the convexity of the higher order term. To avoid repeating similar proofs, these ideas will be used to show that an appropriate discretization of the energy functional $\Gamma-$converge to the continuous total potential energy. Then, using a discrete compactness result we deduce that $\Psi[\cdot]$ admits a minimizer using the Fundamental Theorem of $\Gamma-$convergence [12, 11].

3 Discretization

3.1 Notation

Here we assume for simplicity that the domain $\Omega$ is polygonal and, henceforth, $T_{h}$ denotes the triangulation of the domain $\Omega$ with mesh size $h$. For $K\in T_{h}$, $K$ a triangle, $h_{K}$ is the diameter of $K$ and the mesh size is then defined as $h:=\max_{K\in T_{h}}h_{K}$. The space of polynomials defined on $K$ with total degree less than or equal to $q$ is denoted by $\mathbb{P}_{q}(K)$. Next we will require the partitions of the domain to be shape regular [14], i.e., there exists $c>0$ such that

[TABLE]

where $\rho_{K}$ is the diameter of the largest ball inscribed in $K$.

The boundaries of the elements comprise the set of mesh edges $E_{h}$. The set $E_{h}$ is partitioned into $E^{b}_{h}$, the boundary, and $E^{i}_{h}$, the internal edges, such that $E^{b}_{h}=E_{h}\cap\partial\Omega$ and $E^{i}_{h}=E_{h}\setminus E^{b}_{h}$. For all $e\in E^{i}_{h}$ there exist two distinct elements, we denote them $K_{e^{+}}$ and $K_{e^{-}}$, such that $e\in\partial K_{e^{+}}\cap\partial K_{e^{-}}$. Similarly, if $e\in E^{b}_{h}$, there exists one element $K_{e}$ such that $e\in\partial K_{e}$.

For an edge $e\in\partial K$ we denote by $h_{e}$ its length. Assuming shape regularity it can be shown that there exist constants $C,c$ independent of the mesh size $h$, such that

[TABLE]

To discretize the continuous functional we need to first define our finite element spaces. We use continuous and discontinuous families of Lagrange elements. Consider the space of continuous piecewise polynomial functions $V_{h}^{q}$, viz.

[TABLE]

Also consider the discontinuous finite dimensional space

[TABLE]

We know that $V_{h}^{q}(\Omega)\subset H^{1}(\Omega)$. However, $V_{h}^{q}(\Omega)\not\subset H^{2}(\Omega)$ and thus describing (1.1) over $V_{h}^{q}(\Omega)$ will require the introduction of penalty and jump terms in the discrete functional. Notice that for $u_{h}\in V^{q}_{h}(\Omega)$ we have $\nabla u_{h}\in\tilde{V}^{q-1}_{h}(\Omega)$. In the sequel, we shall use the following notation: The trace of functions in $\tilde{V}^{q}_{h}(\Omega)$ belong to the space

[TABLE]

where we recall that $E_{h}$ is the set of mesh edges. The average and jump operators over $T(E_{h})$ for $\bm{w}\in T(E_{h})^{2\times 2\times 2}$ and $\bm{v}\in T(E_{h})^{2\times 2}$ are defined by:

[TABLE]

[TABLE]

[TABLE]

where $K_{e^{+}}$, $K_{e^{-}}$ are the elements that share the internal edge $e$; $\bm{n}_{e^{+}},\bm{n}_{e^{-}}$ are the corresponding outward normal to the edge and $\bm{v}\otimes\bm{n}_{e}$ is a third order tensor with $\left(\bm{v}\otimes\bm{n}_{e}\right)_{ijk}=\bm{v}_{ij}\bm{n}_{e_{k}}$. Since the boundary terms encode implicitly the boundary conditions we prefer to define them below, see (2.2j). We employ the usual summation convention, e.g. $\bm{v}\cdot\bm{v}=v_{ij}v_{ij}$; also subscripts preceded by a comma indicate partial differentiation with the respect to the corresponding coordinate, e.g. $f_{i,j}=\partial f_{i}=\partial x_{j}$.

3.2 Discretization of the Energy functional

A direct discretization of the minimization problem (1.1) would require an approximation space, a subspace of $H^{2}(\Omega)\times H^{2}(\Omega)$. This means that, for conforming finite elements, we would require $C^{1}$ continuity at the interfaces, i.e. across element internal boundaries. It is well known that the construction of elements that ensure $C^{1}$ continuity is quite complex. Here we adopt to our problem an alternative approach based on the discontinuous Galerkin formulation. Our approximations will be sought on $V_{h}^{q}(\Omega)^{2};$ however the energy functional should be modified to account for possible discontinuities of normal derivatives at the element faces. The appropriate modification of the energy functional proposed below is motivated by the analysis in [43]; the resulting bilinear form of the biharmonic operator obtained via the first variation, will be the form of the $C^{0}$ discontinuous Galerkin method for the linear biharmonic problem, introduced in [27].

The discretized functional for $\bm{u}_{h}\in V^{q}_{h}(\Omega)^{2}$, $q\geq 2$, has the form:

[TABLE]

where the functional $\Psi^{ho}_{h}[\bm{u}_{h}]$ contains the higher order terms. For the boundary faces we use the notation:

[TABLE]

where $\bm{g}$ is given in (2.2c), $I^{q}_{h}$ is the standard nodal interpolation operator in (2.2d) and $\bm{n}$ the outward normal on $\partial\Omega.$ Note that the stabilization term is independent of $\bm{n}_{e}$ because $|\nabla\bm{u}_{h}\otimes\bm{n}_{e}|^{2}=|\nabla\bm{u}_{h}|^{2}$. Although the boundary condition on $\nabla\bm{g}$ is encoded implicitly in the discrete functional through (2.2j), the boundary condition on $\bm{g}$ is enforced explicitly in the discrete space:

[TABLE]

where $\bm{g}_{h}=I^{q}_{h}\bm{g}.$ From now on we shall use the convention that for elements of $\mathbb{A}_{h}^{q}(\Omega)$ the jumps on the boundary faces are given through (2.2j). So, we have to solve the corresponding discrete minimization problem

[TABLE]

Clearly $\nabla\bm{u}_{h}\in\tilde{V}^{q-1}_{h}(\Omega)^{2\times 2}$. On the other hand, $\nabla\nabla\bm{u}_{h}$ does not exist as a function in $L^{2}(\Omega)$ and it can be defined only in the piecewise sense at the element level, i.e. {\left.\kern-1.2pt\nabla\nabla\bm{u}_{h}\vphantom{\big{|}}\right|_{K}}\in\mathbb{P}_{q-2}(K).

Remark 3.1** (Alternative boundary conditions).**

Modeling and simulations of experimental results do not require strict boundary conditions of the function $\bm{u}$ and its derivative. It is observed in the experiments, that when cells contract and pull the ECM fibers, they deform inhomogeneously and change shape under inhomogeneous forces. A minimal model for this deformation has been employed in [32], where the boundary of the cell is connected to the ECM with linear springs. This contributes to the continuous model energy the term

[TABLE]

where $k$ is the stiffness constant and $\Gamma_{i}$ are the boundaries of the cells. See [32, 5.2.2. Model for active particles] for more details. For large enough values of $k$, the above term can model Dirichlet type boundary conditions, in the sense that the displacement is imposed by a stiff linear spring, which is perhaps closer to the experimental approach than the explicit enforcement of the boundary conditions. Then, our analysis remains valid, where now $\bm{u}\in H^{2}(\Omega)^{2}$, $\bm{u}_{h}\in V^{q}_{h}(\Omega)^{2}$ and the terms involving boundary edges in (2.2i) and elsewhere are excluded.

4 Preliminary results

4.1 Preliminary results for finite element spaces

For convenience we briefly state some preliminary results on the finite element spaces which will be useful in the sequel. Following partially the notation of Brenner & Scott, [14], let $\hat{K}=\{(1/h_{K})\bm{x}:\bm{x}\in K\}$ and, for $w\in\mathbb{P}_{q}(K)$, define the function $\hat{w}\in\mathbb{P}_{q}(\hat{K})$ by $\hat{w}(\hat{\bm{x}})=w(h_{K}\hat{\bm{x}})$. Then $w\in W^{p,r}(K)$ is equivalent to $\hat{w}\in W^{p,r}(\hat{K})$ and

[TABLE]

Next we state the well known trace inequality:

Lemma 4.1**.**

Assume that $\Omega$ is bounded and has a Lipschitz boundary. Let $w\in W^{1,p}(\Omega)$, $p\in[1,+\infty]$. Then there exists a constant $C$ depending only on $p$ and $\Omega$ such that

[TABLE]

In general one can have an estimate of the above constant, for instance if $\Omega$ is the unit disk in $\mathbb{R}^{2}$ and $p=2$ then $C\leq 8^{1/4}$, see [14]. We state the discrete trace inequality which is a consequence of Lemma 4.1 and of (2.2a).

Lemma 4.2** (Discrete Trace Inequality).**

Let $T_{h}$ be a shape regular triangulation. Then, there exists a constant $c_{q}$ independent of $h$, but depending on $q$, such that

[TABLE]

The standard nodal interpolation operator will be denoted by $I_{h}^{q}$, where

[TABLE]

Next some well known interpolation error estimates are presented.

Lemma 4.3**.**

Let $\bm{u}\in H^{s}(\Omega)^{2}$, with $s\geq 2$, $T_{h}$ be a shape-regular triangulation of the domain $\Omega$. If $q\geq\left\lceil{s}\right\rceil-1$, where $\left\lceil{s}\right\rceil$ is the smallest integer value greater than or equal to $s$, then there exists a constant $c$ depending only on the domain $\Omega$, a shape parameter of the triangulation and $s$ such that

[TABLE]

where $h_{K}$ denotes the diameter of the element $K$.

Using the trace inequality (2.2b) we obtain the error estimates for norms defined on the mesh edges.

Corollary 4.1**.**

If the assumptions of Lemma 4.3 hold, then for a face $e$ of an element $K$, i.e. $e\in\partial K$ we have the following error estimate of the integral over $e$:

[TABLE]

4.2 Poincaré Inequalities for Broken Sobolev Spaces

To bound $\bm{v}_{h}\in V^{q}_{h}(\Omega)^{2}$ in higher order norms we will need Poincaré inequalities for broken Sobolev spaces. For this purpose we define the broken Sobolev seminorm for $\bm{w}\in\mathbb{A}_{h}^{q}(\Omega)$:

[TABLE]

Since $\Omega\subset\mathbb{R}^{2}$ is an open, connected and bounded set with Lipschitz boundary, the main result in [39] and in [15] implies that

[TABLE]

for all $r\in[1,+\infty)$, where $\widehat{\nabla\bm{w}}$ can take any of the following forms

[TABLE]

The constant $c_{\theta_{min}}$ depends on the minimum angle of the triangles, $\theta_{min}$, and $r$. We have assumed that the family of partitions $\{T_{h}\}$ is shape regular, consequently $\theta_{min}$ is independent of $h$, see [39] for details.

4.3 Lifting and the Discrete Gradient Operators

We adopt in our case of the discontinuous gradient the results of [8, 18] regarding Lifting and Discrete gradients.

Definition 4.1** (The vectorial piecewise gradient).**

Let $\bm{w}\in\tilde{V}^{k}_{h}(\Omega)^{m}$, $\Omega\subset\mathbb{R}^{n}$. The vectorial piecewise gradient $\nabla_{h}:\tilde{V}^{k}_{h}(\Omega)^{m}\rightarrow\tilde{V}^{k-1}_{h}(\Omega)^{m\times n}$ is defined componentwise as

[TABLE]

Consequently for all $\bm{u}_{h}\in V_{h}^{q}(\Omega)^{2}$, $\nabla_{h}\nabla\bm{u}_{h}\in L^{2}(\Omega)^{2\times 2\times 2}$ and

[TABLE]

For all $e\in E_{h}$ we define the linear operator, known as lifting operator [8, 18], $\bm{r}_{e}:L^{2}(e)^{2\times 2}\rightarrow\tilde{V}^{q-2}_{h}(\Omega)^{2\times 2\times 2}$, for all $\bm{\phi}\in L^{2}(e)^{2\times 2}$ as follows

[TABLE]

It can be shown that $\bm{r}_{e}(\phi)$ is non zero only on the elements that contain $e$ on their boundary, i.e. $supp(\bm{r}_{e})=\{K\in T_{h}:e\in\partial K\}$. We also define the global lifting operator

[TABLE]

With the help of this operator we can represent the second term of the functional $\Phi^{ho}_{h}$, see (2.2i), as an integral over $\Omega$; namely,

[TABLE]

Since $\bm{u}_{h}\in\mathbb{A}_{h}^{q}(\Omega)$, we can substitute $\bm{w}_{h}$ by $\nabla_{h}\nabla\bm{u}_{h}$ to obtain

[TABLE]

where recall that at the boundary faces $e\in E^{b}_{h}$ we use the convention $\llbracket\nabla\bm{u}_{h}\rrbracket=\nabla\bm{u}_{h}^{-}-\nabla I_{h}^{q}\bm{g}\,,$ and $\{\kern-3.125pt\{\nabla\nabla\bm{u}_{h}\}\kern-3.125pt\}=\nabla\nabla\bm{u}_{h}^{-}\,.$

We next define the discrete gradient $G_{h}$, which is a combination of the vectorial piecewise gradient and the global lifting operator. In particular, the discrete gradient of $\nabla\bm{u}_{h}$ is defined as

[TABLE]

Later we will show under which conditions $G_{h}(\nabla\bm{u}_{h})\rightharpoonup\nabla\nabla\bm{u}$ in $L^{2}(\Omega)^{2\times 2\times 2}$, when $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$ as $h\rightarrow 0$. For this reason it will be convenient to write the higher order derivatives of the discrete total potential energy in terms of the discrete gradient and the global lifting operator. We do this in Lemma 4.4 below which will be useful in the sequel.

Lemma 4.4**.**

The higher order terms of the discretized total potential energy can be written in terms of the discrete gradient and the global lifting operator, as follows

[TABLE]

Proof.

The proof follows from Definition 4.1 and equation (2.2n). ∎

From Lemma 4.4 the functional $\Psi^{ho}_{h}$, displayed in (2.2i), becomes

[TABLE]

Using inverse inequalities (2.2c), one can show, see [18], that there exists a positive constant $C$, independent of h, such that

[TABLE]

This bound finally implies the next Lemma, see [26, Lemma 4.34] and [19, Lemma 7] for details.

Lemma 4.5** (Bound on global lifting operator).**

For all $\bm{u}_{h}\in\mathbb{A}_{h}^{q}(\Omega)$ there holds

[TABLE]

where the constant $C_{R}$ depends on the constant of (2.2r).

Corollary 4.2** (Bound on Discrete Gradient).**

There exists a constant $C>0$ such that for all $\bm{u}_{h}\in\mathbb{A}_{h}^{q}(\Omega)$ it holds that

[TABLE]

4.4 Analytical preliminaries

In the subsequent sections we examine the convergence of the lower order terms in $L^{1}(\Omega)$, i.e. the terms $W(\cdot)$ and $\Phi(\cdot)$ of (2.2i), when $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$. For this purpose we employ Vitali’s convergence theorem [49].

Theorem 4.1** (Vitali convergence theorem).**

Let $\Omega$ be a set of finite measure and $(f_{n})$ be a sequence of functions in $L^{1}(\Omega)$. Assume $(f_{n})$ is uniformly integrable over $\Omega$, i.e., for each $\epsilon>0$, there exist $\delta>0$ independent of $n$ such that $\text{if }A\subset\Omega\text{ measurable and }|A|<\delta,\text{ then }\int_{A}f_{n}d\bm{x}<\epsilon.$ If $(f_{n})\rightarrow f$ pointwise a.e. on $\Omega$, then $f\in L^{1}(\Omega)$ and

[TABLE]

We now state a useful criterion of uniform integrability, known as the de la Vallée Poussin criterion, [9, Theorem 4.5.9].

Theorem 4.2** (de la Vallée Poussin criterion).**

A family $(f_{n})\subset L^{1}(\Omega)$, is uniformly integrable if, and only if, there exists a non-negative increasing function $G$ on $[0,+\infty)$ such that

[TABLE]

5 $\Gamma-$Convergence of the discretization

In this section we establish the $\Gamma$-convergence of the discretized functionals $\Psi_{h}$ to the continuum energy $\Psi$. The proof consists of three parts: we first prove equi-coercivity, a necessary property to bound the sequence $(\bm{u}_{h})$ when the family of discrete energies $(\Psi_{h}[\bm{u}_{h}])$ is bounded. Then, we show the $\liminf$ inequality which provides a lower bound of the discrete energies by the continuum counterpart. We conclude with the $\limsup$ inequality which, as we will see in Section 6 ensures the attainment of the limit.

5.1 Equi-Coercivity and Convergence of the Discrete Gradient

Proposition 5.1** (Equi-coercivity).**

Assume that ${\alpha}>C_{R}$, i.e. the penalty parameter in (2.2i) is greater than the constant of Proposition 4.5. Let $(\bm{u}_{h})_{h>0}$ be a sequence of displacements in $\mathbb{A}^{q}_{h}(\Omega)$ such that for a constant $C>0$ independent of $h$ it holds that

[TABLE]

Then the there exists a constant $C_{1}>0$ such that

[TABLE]

In addition,

[TABLE]

for a positive constant $C_{2}$, where $C_{1},C_{2}$ are independent of $h$.

Proof.

We have shown, see (2.2i) and (2.2q), that $\Psi^{ho}[\bm{u}_{h}]$ can be written in terms of the discrete gradient $G_{h}$ and the lifting operator $\bm{R}_{h}$. From the assumption $a>C_{R}$ and the bound of the global lifting operator in (2.2s), we see that $\Psi^{ho}[\bm{u}_{h}]$ is nonnegative:

[TABLE]

In particular,

[TABLE]

By (2.1) it holds $W(\nabla\bm{u}_{h})\geq-c$, $c$ is a positive constant, and $\Phi(\nabla\bm{u}_{h})$ is a nonnegative penalty parameter, consequently

[TABLE]

From the assumption that $\Psi_{h}[\bm{u}_{h}]$ is uniformly bounded over $h$, i.e. $\Psi_{h}[\bm{u}_{h}]\leq C$ for all $h>0$, we obtain that all terms appearing in the right hand sides of (2.2d) and (2.2e) are uniformly bounded. Therefore, using the bound of the global lifting operator, (2.2s), we conclude that

[TABLE]

and

[TABLE]

It remains to show that $\lVert\bm{u}_{h}\rVert_{H^{1}(\Omega)^{2}}$ is uniformly bounded, using that $\bm{u}_{h}\in\mathbb{A}_{h}^{q}(\Omega)$. By the coercivity condition on $W$ given in (2.1), equation (2.2e) gives

[TABLE]

where we have used the Cauchy-Schwarz and Young’s inequality. Choosing, for example, $\delta=1/2$, we infer that $C\geq\frac{c_{0}}{2}\int_{\Omega}\left(|\nabla\bm{u}_{h}|^{2}-c\right)d\bm{x}$ and the proof is concluded by Poincaré’s inequality. ∎

5.2 The $\liminf$ inequality

Lemma 5.1** (Convergence of the lower order terms).**

Let $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$, with $\bm{u}_{h}\in H^{1}(\Omega)^{2}$ and $\bm{u}\in H^{2}(\Omega)^{2}$. Suppose further that

[TABLE]

where $C$ is independent of $h$. Then, up to a subsequence,

[TABLE]

Proof.

From the assumption $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$ there exists a subsequence, not relabeled, such that $\nabla\bm{u}_{h}\rightarrow\nabla\bm{u}$ a.e. Note that $W$, see (2.1), is bounded from above. Specifically

[TABLE]

for some $m\in\mathbb{N}$. Since $(\nabla\bm{u}_{h})$ is bounded in all $L^{r}$ norms, by (2.2k), this implies that $(W(\nabla u_{h}))$ is uniformly integrable and from Vitali’s Theorem 4.1 must converge to $W(\nabla\bm{u})$ in $L^{1}(\Omega)$. ∎

Remark 5.1** (Convergence of the penalty term).**

From inequality (2) the penalty term is bounded from above, i.e.

[TABLE]

Then as in Lemma 5.1 $\Phi(\nabla\bm{u}_{h})\rightarrow\Phi(\nabla\bm{u})$, as $h\rightarrow 0$ up to a subsequence.

Lemma 5.2**.**

Let $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$, with $\bm{u}_{h}\in\mathbb{A}^{q}_{h}(\Omega)$. If $|\bm{u}_{h}|_{H^{2}(\Omega,T_{h})}$ is uniformly bounded with respect to $h$, then

[TABLE]

Proof.

This proof is an adaptation of the proof in [25, Theorem 2.2]. Let $\bm{\phi}\in C_{c}^{\infty}(\Omega)^{2\times 2\times 2}$, then from the definition of the piecewise gradient, equation (2.2j), and the divergence theorem, we obtain

[TABLE]

where we recall that $(\nabla\bm{u}_{h}\otimes\bm{n})_{ijk}=(\bm{u}_{h})_{i,j}\bm{n}_{k}$.

First, we show that the last two terms convergence to [math] as $h\rightarrow 0$. To this end, let $\bar{I}^{0}_{h}$ be the piecewise average operator onto $\tilde{V}^{0}_{h}(\Omega)^{2\times 2\times 2}$, i.e. {\left.\kern-1.2pt\bar{I}^{0}_{h}\bm{\phi}\vphantom{\big{|}}\right|_{K}}=1/|K|\int_{K}\bm{\phi}, and define $\bm{\phi}_{h}=\bar{I}^{0}_{h}\bm{\phi}$. Then from standard error estimates, see e.g. [26, Lemma 1.58],

[TABLE]

Now, using the definition of $\bm{R}_{h}$, see (2.2m), for the last two terms of (2.2o) we obtain

[TABLE]

The assumption that $|\bm{u}_{h}|_{H^{2}(\Omega,T_{h})}$ is bounded for all $h>0$, implies that $\sum_{e\in E_{h}}h^{-1}_{e}\lVert\llbracket\nabla\bm{u}_{h}\rrbracket\rVert^{2}_{L^{2}(e)^{2\times 2}}$ is also uniformly bounded. Also, from Proposition 4.5 the global lifting operator is bounded from the jump terms, inequality (2.2s). As a result,

[TABLE]

and $(\bm{R}_{h}(\nabla\bm{u}_{h}))$ is uniformly bounded in the $L^{2}(\Omega)^{2\times 2\times 2}$ norm. From the last relation and the Cauchy-Schwarz inequality we bound $I_{2}$:

[TABLE]

Hence, $I_{2}\rightarrow 0$, as $h\rightarrow 0$ by the error estimate given in (2.2p). Similarly for the term $I_{1}$, using the previous uniform bound for the jump terms, the error estimate for integrals that are defined over an edge $e$, (2.2f), and letting $K_{e}=\{K\in T_{h}:e\in\partial K\}$, we infer that

[TABLE]

Since $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$, taking the limit as $h\rightarrow 0$, employing (2.2s) and (2.2t), the relation (2.2o) implies

[TABLE]

∎

Corollary 5.1**.**

*(Weak Convergence of the Discrete Gradients).

Suppose the assumptions of Lemma 5.2 hold. In addition assume that $\bm{u}\in H^{2}(\Omega)^{2}$, then*

[TABLE]

Proof.

Lemma 5.2 ensures the limit of equation (2.2n) for all $\bm{\phi}\in C^{\infty}_{c}(\Omega)$, hence

[TABLE]

∎

Theorem 5.1** (The $\liminf$ inequality.).**

Assume that ${\alpha}>C_{R}$, i.e. that the parameter of the discrete energy function is larger than the constant of (2.2s). Also, let the penalty term satisfy condition (2). Then for all $\bm{u}\in\mathbb{A}(\Omega)$ and all sequences $(\bm{u}_{h})\subset\mathbb{A}^{q}_{h}(\Omega)$ such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$ it holds that

[TABLE]

Proof.

We assume there is a subsequence, still denoted by $\bm{u}_{h}$, such that $\Psi_{h}[\bm{u}_{h}]\leq C$ uniformly in $h$, otherwise $\Psi[\bm{u}]\leq\liminf\limits_{h\rightarrow 0}\Psi_{h}[\bm{u}_{h}]=+\infty$. The following steps conclude the proof:

From Proposition 5.1, the uniform bound $\Psi_{h}[\bm{u}_{h}]\leq C$ implies that $|\bm{u}|_{H^{2}(\Omega,T_{h})}$ and $\lVert\bm{u}_{h}\rVert_{H^{1}(\Omega)}$ are uniformly bounded. 2. 2.

Corollary 5.1 implies $G_{h}(\nabla\bm{u}_{h})\rightharpoonup\nabla\nabla\bm{u}$ in $L^{2}(\Omega)^{2\times 2\times 2}$. 3. 3.

The term $\int_{\Omega}|G_{h}|^{2}$ is convex which implies weak lower semicontinuity [21]: $\liminf\limits_{h\rightarrow 0}\int_{\Omega}|G_{h}|^{2}\geq\int_{\Omega}|\nabla\nabla\bm{u}|^{2}$. 4. 4.

From the Poincaré inequality for broken Sobolev spaces, (2.2h), there exist a constant $c$ independent of $h$ such that for all $m\in[1,+\infty)$ it holds

[TABLE]

where the last bound holds from step 1. 5. 5.

Finally, the assumed convergence of $\bm{u}_{h}$ to $\bm{u}$ in $H^{1}(\Omega)^{2}$, Lemma 5.1 and Remark 5.1 ensure the convergence of the remaining terms, i.e. $\int_{\Omega}W(\nabla\bm{u}_{h})+\Phi(\nabla\bm{u}_{h})\rightarrow\int_{\Omega}W(\nabla\bm{u})+\Phi(\nabla\bm{u})$.

As a result we deduce that $\Psi[\bm{u}]\leq\liminf\limits_{h\rightarrow 0}\Psi_{h}[\bm{u}_{h}]$. ∎

5.3 The $\limsup$ inequality

In this section we focus on the $\limsup$ inequality. Given $\bm{u}\in\mathbb{A}(\Omega)$, we would like to prove the existence of a sequence $(\bm{u}_{h})\subset\mathbb{A}_{h}^{q}(\Omega)$, such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$ and $\Psi[\bm{u}]\geq\limsup\limits\Psi_{h}[\bm{u}_{h}]$ as $h\rightarrow 0$. To this end, choose a sequence of smooth functions $\bm{u}_{\delta}$ such that $\bm{u}_{\delta}\rightarrow\bm{u}$ in $H^{2}(\Omega)$ and we define $\bm{u}_{h}$ to be an appropriate interpolant of $\bm{u}_{\delta}$, $I^{q}_{h}\bm{u}_{\delta}$, for a chosen $\delta$. A similar strategy was used in [6]. We start with the following Proposition:

Proposition 5.2**.**

For $\bm{u}\in H^{s}(\Omega)^{2}$, if $s>2$ and $q\geq 2$, there exist constants $C_{1},C_{2}>0$ such that

[TABLE]

Proof.

To simplify the notation let $\bm{v}_{h}:=I^{q}_{h}\bm{u}$. First we study (2.2y). Adding and subtracting the term $\nabla\bm{u}$, using the error estimates of Corollary 4.1, yields

[TABLE]

where $K^{+}_{e},K^{-}_{e}$ denote the distinct elements that share the edge $e$, i.e. $e=K^{+}_{e}\cap K^{-}_{e}$. Summing over all edges $e\in E^{i}_{h}$ inequality (2.2aa) gives

[TABLE]

To prove (2.2z) notice that

[TABLE]

The Cauchy Schwarz inequality and the discrete trace inequality (2.2c) imply

[TABLE]

Therefore

[TABLE]

However,

[TABLE]

where we have used the error estimates of (2.2e). Consequently (2.2ae) becomes

[TABLE]

Similarly we show that

[TABLE]

∎

Theorem 5.2** (The $\limsup$ inequality.).**

Let the penalty function $\Phi$ satisfy condition (2). The following property holds: For all $\bm{u}\in\mathbb{A}(\Omega)$, there exists a sequence $(\bm{u}_{h})_{h>0}$ with $\bm{u}_{h}\in\mathbb{A}^{q}_{h}(\Omega)$, such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$ and

[TABLE]

Proof.

Given $\bm{u}\in\mathbb{A}(\Omega)$ we construct an appropriate sequence $(\bm{u}_{h})_{h>0}\subset\mathbb{A}^{q}_{h}(\Omega)$, the recovery sequence, such that $\lVert\bm{u}_{h}-\bm{u}\rVert_{H^{1}(\Omega)^{2}}\rightarrow 0$ and $\lim_{h\rightarrow 0}\Psi_{h}[\bm{u}_{h}]=\Psi[\bm{u}]$.

In particular, we approximate $\bm{u}$ via mollification by a sequence of smooth functions in $\overline{\Omega}$ to find that, for all $\delta>0$, there exists $\bm{u}_{\delta}\in H^{3}(\Omega)^{2}\cap\mathbb{A}(\Omega)$ such that

[TABLE]

Recall, $\bm{u}_{\delta}=\bm{u}=\bm{g}$ and $\nabla\bm{u}_{\delta}=\nabla\bm{u}=\nabla\bm{g},$ on $\partial\Omega$, where $\bm{g}\in H^{3}(\Omega)$ by the definition of $\mathbb{A}(\Omega)$. Here $c$ is independent of $\delta$.

Next, define $\bm{u}_{h,\delta}:=I^{q}_{h}\bm{u}_{\delta}\in H^{1}(\Omega)^{2}$, noting that $\bm{u}_{h,\delta}\in\mathbb{A}^{q}_{h}(\Omega)$. From the error estimates in (2.2e) and the fact that $q\geq 2$ we find that

[TABLE]

Inequality (2.2aj) and the error estimate of (2.2ak) imply that

[TABLE]

Choosing

[TABLE]

we deduce that the sequence $(\bm{u}_{h,\delta_{h}})\subset V^{q}_{h}(\Omega)$ and $\bm{u}_{h,\delta_{h}}\rightarrow\bm{u}$ in $H^{1}(\Omega)$, as $h\rightarrow 0$. It remains to prove

[TABLE]

For the convergence of the lower order terms we use Lemma 5.1 and Remark 5.1. Hence, it suffices to show that $\lVert\nabla\bm{u}_{h,\delta_{h}}\rVert_{L^{q}(\Omega)^{2\times 2}}$, is uniformly bounded with respect to $h$. Indeed, Sobolev’s embedding theorem and classical error estimates, see inequality (2.2al), imply that $\bm{u}_{h,\delta_{h}}$ is uniformly bounded in ${W^{1,q}(\Omega)}$ with respect to $h$, for all $q\in[1,\infty)$:

[TABLE]

To extend the bound over $\Omega$, assume that $q$ is integer and $q\geq 2$; then the multinomial formula implies

[TABLE]

The result can be extended for $q=1$ because $|\Omega|<+\infty$ and for $q\in[1,+\infty)$ using the interpolation bound, [17]:

[TABLE]

where $\lfloor\cdot\rfloor,\lceil\cdot\rceil$ are the floor and ceiling integer functions respectively and $\xi=\frac{(\lceil q\rceil-q)\lfloor q\rfloor}{q}$. Notice that Lemma 5.1 and Remark 5.1 give

[TABLE]

in $L^{1}(\Omega)$, as $h\rightarrow 0$. Also, since we have choosen $\delta_{h}=\sqrt{h}$, Proposition 5.2, (2.2aj), implies that as $h\rightarrow 0$

[TABLE]

For the boundary terms we have:

[TABLE]

Similarly,

[TABLE]

For the remaining higher order terms, we work similarly as before. We begin, showing an inequality for a given element $K$ using the error estimates of (2.2al) and (2.2am):

[TABLE]

Then summing over $K$ in $T_{h}$ and using the Cauchy-Schwarz inequality, gives:

[TABLE]

$\text{as }h\rightarrow 0$, where for the last inequality we have used (2.2aj) and (2.2ao). Furthermore, note that

[TABLE]

Finally (2.2aw) and (2.2ax) give

[TABLE]

which concludes (2.2ap). ∎

6 Compactness and Convergence of Discrete Minimizers

In this section our main task is to use the results of the previous section to show that under some boundedness hypotheses on $\bm{u}_{h}$, a sequence of discrete minimizers $(\bm{u}_{h})$ converges in $H^{1}(\Omega)$ to a global minimizer $\bm{u}$ of the continuous functional, Theorem 6.1. Such results are standard in the $\Gamma-$convergence literature, [11, 22], but the application in our setting is not straightforward. The main reason is that in our case we need certain intermediate results, such as a discrete DG version of the Rellich-Kondrachov theorem, which we show in the sequel. We use related discrete bounds derived previously in [40, 29, 25, 19]

We will use the total variation of a function $v\in L^{1}(\Omega)$, see [28], defined as

[TABLE]

The space of functions of bounded variation in $\Omega$, denoted by $BV(\Omega)$, contains all $L^{1}(\Omega)$ functions with bounded total variation, i.e.,

[TABLE]

The space $BV(\Omega)$ is endowed with the norm $\lVert v\rVert_{BV}=\lVert v\rVert_{L^{1}(\Omega)}+|Dv|(\Omega)$. The following inequality plays a key role in the desired compactness.

Lemma 6.1** (Bounds for the total variation).**

Let $\bm{w}\in V^{q}_{h}(\Omega)^{2}$. Then, there exist a constant $C$ independent of $h$ such that

[TABLE]

A proof can be found in [40, Theorem 3.26] and a generalization in ([19, Lemma 2]. It is based on the observation

[TABLE]

where $\bm{\phi}\in C_{c}^{1}(\Omega)^{2\times 2\times 2}$ and appropriate bounds on the right-hand side.

Proposition 6.1** (Discrete Rellich-Kondrachov).**

Let a sequence $(\bm{u}_{h})\subset\mathbb{A}^{q}_{h}(\Omega)$ be bounded, for $C>0$, as

[TABLE]

Then $(\bm{u}_{h})$ is relatively compact in $W^{1,p}(\Omega)^{2}$ for $1\leq p<+\infty$, i.e. there exists a $\bm{u}\in W^{1,p}(\Omega)^{2}$ such that

[TABLE]

up to a subsequence.

Proof.

Using the Sobolev embedding theorem, the discrete Poincaré inequality (2.2h) and inequality (2.2e) we conclude that

[TABLE]

uniformly with respect to $h$, for all $r\in[1,+\infty)$, where $C$ is a positive constant independent of $h$. The space $W^{1,r}(\Omega)$ is reflexive for $r\in(1,+\infty)$. Therefore, for every bounded sequence there exists a subsequence, not relabeled, and a function $\bm{u}\in W^{1,r}(\Omega)$ such that

[TABLE]

From the Rellich-Kondrachov theorem, [17, Theorem 9.16], and since $\dim{\Omega}=2$, it is known that ${H^{1}(\Omega)\subset\subset L^{p}(\Omega)}$, i.e. it is compactly embedded for $1\leq p<+\infty$. In particular, $\bm{u}_{h}\rightarrow\bm{u}$ in $L^{p}(\Omega)$ for $p\in[1,+\infty)$.

It remains to prove that $(\nabla\bm{u}_{h})$ is relatively compact in $L^{p}(\Omega)$. Similar results have been proved in [25, 19, 29], all of them are based on a bound of the $BV$ norm, $\lVert\cdot\rVert_{BV(\Omega)}$, from the higher order terms. To this end, from inequalities (2.2c) and (2.2e) it follows that $(\nabla\bm{u}_{h})$ is uniformly bounded in $BV(\Omega)^{2\times 2}$. By standard embedding theorems, see for example [28, Theorem 5.5], $\left(\nabla\bm{u}_{h}\right)$ in $BV(\Omega)^{2\times 2}$ is relatively compact in $L^{1}(\Omega)^{2\times 2}$ and, up to a subsequence, there exists $\bm{w}\in BV(\Omega)^{2\times 2}$ such that

[TABLE]

Therefore, from the interpolation inequality, [17], and (2.2g), we obtain that

[TABLE]

where $p\in(1,r)$ and $\theta=\frac{r-p}{p(r-1)}\in(0,1)$. From (2.2h) we conclude that $\bm{u}_{h}\rightarrow\bm{u}$ in $W^{1,p}(\Omega)$ for $p\in[1,+\infty)$. ∎

The discrete Rellich-Kondrachov Theorem ensures that there exist $\bm{u}$ in $H^{1}(\Omega)^{2\times 2}$ such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$ under some boundedness hypotheses on $\bm{u}_{h}$. However, the proof that minimizers of the discrete problem converge to a minimizer of the continuous problem would require higher regularity on $\bm{u}$, i.e. $\bm{u}\in H^{2}(\Omega)^{2}$. Similar arguments were used previously in [25, 19].

Proposition 6.2** **(Regularity of the limit and Weak Convergence of the Discrete

Gradient).

Let a sequence $(\bm{u}_{h})\subset\mathbb{A}^{q}_{h}(\Omega)$ with $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$. If the sequence is bounded in the $H^{2}(\Omega,T_{h})$ seminorm, then

[TABLE]

*up to a subsequence. In addition, $\bm{u}\in\mathbb{A}(\Omega)$. *

Proof.

Here we adopt partially the proof of [25]. From inequality (2.2t), the discrete gradient $G_{h}(\nabla\bm{u}_{h})$ is bounded by $|\bm{u}_{h}|_{H^{2}(\Omega,T_{h})}$ in the $L^{2}(\Omega)^{2\times 2\times 2}$ norm. Hence, there exists $\bm{w}\in L^{2}(\Omega)^{2\times 2\times 2}$ such that, up to a subsequence,

[TABLE]

To prove that $\bm{w}=\nabla\nabla\bm{u}$, let $\bm{\phi}\in C^{\infty}_{c}(\Omega)^{2\times 2\times 2}$. By Lemma 5.2

[TABLE]

which means that $\bm{w}=\nabla\nabla\bm{u}$.

It remains to prove that $\bm{u}\in\mathbb{A}(\Omega)$. If $\bm{u}_{h}\in\mathbb{A}^{q}_{h}(\Omega)$ then {\left.\kern-1.2pt\bm{u}_{h}\vphantom{\big{|}}\right|_{\partial\Omega}}={\left.\kern-1.2pt(I^{q}_{h}\bm{g})\vphantom{\big{|}}\right|_{\partial\Omega}} and $\bm{g}\in{H^{3}(\Omega)^{2}}$. Classical interpolation error estimates, see inequality (2.2f), give

[TABLE]

Combining the last result with the trace inequality yields

[TABLE]

which means $\bm{u}=\bm{g}$ a.e on $\partial\Omega$. For the $\nabla\bm{g}$ term we first notice that

[TABLE]

On the other hand, Lemma 4.1 implies ($K_{e}$ is the element with boundary face $e$)

[TABLE]

The proof is thus complete. ∎

Theorem 6.1** (Convergence of discrete absolute minimizers).**

Assume that ${\alpha}>C_{R}$, i.e. that the stabilization parameter is greater than the constant of (2.2s). Let $(\bm{u}_{h})\subset\mathbb{A}^{q}_{h}(\Omega)$ be a sequence of absolute minimizers of $\Psi_{h}$, i.e.,

[TABLE]

If $\Psi_{h}[\bm{u}_{h}]$ is uniformly bounded then, up to a subsequence, there exists $\bm{u}\in\mathbb{A}(\Omega)$ such that

[TABLE]

Remark 6.1**.**

We note that for each fixed $h>0$ the functional $\Psi_{h}$ admits a minimizer. Indeed, given an infimising sequence in the finite-dimensional space $\mathbb{A}^{q}_{h}(\Omega)$, by coercivity (Proposition 5.1) and the Poincaré inequality (2.2h), it must be bounded in the norm

[TABLE]

Finite-dimensionality of $\mathbb{A}^{q}_{h}(\Omega)$ implies strong convergence of the infimising sequence in the above norm and it is easy to check that the functional $\Psi_{h}$ is continuous with respect to this convergence.

Proof of Theorem 6.1.

The uniform bound for the discrete energies implies from the equi-coercivity property, Proposition 5.1, that

[TABLE]

uniformly with respect to $h$. The discrete Rellich-Kondrachov, Proposition 6.1, ensures that there exists $\bm{u}\in H^{1}(\Omega)^{2}$ such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$ up to a subsequence not relabeled here. Since $|\bm{u}_{h}|_{H^{2}(\Omega,T_{h})}$ is uniformly bounded, Proposition 6.2 implies that $\bm{u}\in H^{2}(\Omega)^{2}$ and also $\bm{u}\in\mathbb{A}(\Omega)$.

To prove that $\bm{u}$ is a global minimizer of $\Psi$ we use the $\liminf$ and $\limsup$ inequalities, Theorems 5.1 and 5.2 respectively. Let $\bm{w}\in\mathbb{A}(\Omega)$, then the $\limsup$ inequality implies that there exist $\bm{w}_{h}\in\mathbb{A}^{q}_{h}(\Omega)$ such that

[TABLE]

Therefore, since $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$ the $\liminf$ inequality and the fact that $\bm{u}_{h}$ are absolute minimizers of the discrete problems imply that

[TABLE]

for all $\bm{w}\in\mathbb{A}(\Omega)$. Therefore $\bm{u}$ is an absolute minimizer of $\Psi$. ∎

7 Incorporating Penalty terms with Exponential Growth

So far we have shown the convergence of discrete absolute minimizers to a global minimizer of the continuous problem. The proof is based on the $\Gamma-$convergence of $\Psi_{h}$ to $\Psi$ and on discrete compactness results. To penalize interpenetration of matter we have added a penalty function $\Phi$ in the total potential energy (1.1), assuming that $\Phi$ has polynomial growth, see (2). One can notice that the polynomial growth penalizes also deformations where $J>1$, $J=\det(\bm{1}+\nabla\bm{u})$ is the Jacobian determinant of the mapping, and as a consequence can affect the material properties.

Using a penalty of the form

[TABLE]

for large enough ${A}>0$, the penalty parameter contributes to the total potential energy when $J<1$, thus an assumption as (2.2a) seems preferable. In addition, computational results related especially to densified phase and its comparison to experimental data indicates that (2.2a) is a better choice, for more details see [30].

In this section we will assume that instead of polynomial growth, (2.2a) holds. However, employing a penalty function $\Phi$ with exponential growth, affects the proofs of the inequalities $\liminf$, Theorem 5.1, and $\limsup$, Theorem 5.2. The main technical difficulty addressed in this section is the proof of the $\limsup$ inequality when (2.2a) is assumed. To show the analog of Theorem 5.2, and in particular that $\Phi$ is uniformly integrable one has to use appropriate Orlicz spaces and corresponding embedding results. To do this in discrete DG spaces requires new ideas, which we describe in this section. As far as we know these bounds are the first embedding estimates for DG spaces using the Orlicz framework.

Below we recall some definitions and basic properties for Orlicz spaces which we require. The reader is referred to [47] for a thorough review of Orlicz spaces. Let $\varphi:{\mathbb{R}}\to[0,+\infty]$ be a continuous, convex and even function satisfying

[TABLE]

The Orlicz class $L_{\varphi}(\Omega)^{2}$ consists of all measurable functions $\bm{u}:\Omega\to{\mathbb{R}}^{2}$ such that

[TABLE]

The Orlicz space $L_{\varphi^{*}}(\Omega)^{2}$ is the linear span of functions in $L_{\varphi}(\Omega)^{2}$ and it becomes a Banach space when equipped with the Luxemburg norm

[TABLE]

Remark 7.1**.**

Let $E_{\varphi}(\Omega)^{2}$ denote the closure of $L^{\infty}(\Omega)^{2}$ in $L_{\varphi^{*}}(\Omega)^{2}$. The convexity of $\varphi$ implies $E_{\varphi}(\Omega)^{2}\subset L_{\varphi}(\Omega)^{2}$, see [47]. Moreover, given a sequence $(\bm{u}_{k})\subset L_{\varphi}(\Omega)^{2}$ and $\bm{u}\in L_{\varphi^{*}}(\Omega)^{2}$, we say that $\bm{u}_{k}$ is mean convergent to $\bm{u}$ if

[TABLE]

Norm convergence in $L_{\varphi^{*}}(\Omega)^{2}$ is stronger than mean convergence.

Next, we equip the space $\tilde{V}_{h}^{r}(\Omega)^{2}$ with the norn

[TABLE]

where

[TABLE]

for all $\bm{w}\in\tilde{V}_{h}^{r}(\Omega)^{2}$. Our strategy relies on proving an embedding theorem of the space $\tilde{V}_{h}^{r}(\Omega)^{2}$, for all $r\geq 1$, into the Orlicz space $L_{\varphi^{*}}(\Omega)$ where

[TABLE]

This embedding is proved in Theorem 7.1 below, which extends Trudinger’s embedding theorem for Orlicz spaces, [53, Theorem 2], to the DG finite element setting.

It will be useful to use the reconstruction operator of Karakashian and Pascal, [35]. The next result is a local version of its approximation properties established in [35].

Lemma 7.1**.**

Let $\omega_{e}(\varv)$ denote the set of edges that contain the node $\varv$, i.e. $\omega_{e}(\varv)=\{e\in E^{i}_{h}|\varv\in e\}$. Then for $\bm{u}\in\tilde{V}_{h}^{r}(\Omega)^{2}$ there exists a reconstruction operator $Q:\tilde{V}_{h}^{r}(\Omega)^{2}\rightarrow V^{r}_{h}(\Omega)^{2}$ such that

[TABLE]

where $\omega_{e}(K)=\cup_{\varv\in K}\omega_{e}(\varv)$.

The proof requires an appropriate adaptation of [36, Theorem 2.1], see also [24]. Let $\omega_{\varv}=\{K_{1},...,K_{|\omega_{\varv}|}\}$, where every consequence pair share the edge $e_{i}^{\varv}=K_{i}\cap K_{i+1}$, when $i=1,...,|\omega_{\varv}|-1.$ Using the explicit definition of $Q,$ as an averaging operator, [36], appropriate scaling arguments and the inverse inequality

[TABLE]

the proof relies on establishing

[TABLE]

For details see [30] and previous versions of this manuscript. The next result will be useful.

Lemma 7.2**.**

Let $\bm{u}\in\tilde{V}_{h}^{r}(\Omega)^{2}$, for all $\bm{x}\in\Omega$ it holds

[TABLE]

where $\hat{\bm{u}}=Q\bm{u}$ and $Q$ the reconstruction operator of Lemma 7.1.

Proof.

For all $\bm{x}\in\Omega$ we have

[TABLE]

$T_{h}$ is a triangulation of the domain $\Omega$, consequently there exists $K^{*}\in T_{h}$ such that $\lVert\bm{u}-\hat{\bm{u}}\rVert_{L^{\infty}(\Omega)}=\lVert\bm{u}-\hat{\bm{u}}\rVert_{L^{\infty}(K^{*})}$. Using inequality (2.2e) we deduce that

[TABLE]

∎

Next we state a crucial lemma, stated in [53, Lemma 1], for a complete proof of the embedding theorem.

Lemma 7.3**.**

Let $\Omega\subset{\mathbb{R}}^{2}$ satisfy the cone condition and $\bm{w}\in H^{1}(\Omega)^{2}$. Then for a.e. $\bm{x}\in\Omega$ it holds that

[TABLE]

We may now prove the embedding.

Theorem 7.1**.**

Let $\Omega\subset R^{2}$ satisfy the cone condition. For each $h>0$, the space $\tilde{V}^{r}_{h}(\Omega)^{2}$ is continuously embedded into the Orlicz space $L_{\varphi^{*}}(\Omega)$ where

[TABLE]

Furthermore, for any $\psi$ such that $\psi(t)\leq\varphi(\lambda t)$ for some $\lambda>0$, the space $\tilde{V}^{r}_{h}(\Omega)^{2}$ is continuously embedded, in the sense of mean convergence, into the Orlicz class $L_{\psi}(\Omega)$; i.e., whenever $\lVert\bm{u}_{k}-\bm{u}\rVert_{H^{1}(\Omega,T_{h})}\to 0$ then

[TABLE]

Proof.

We exhibit the existence of constants $b=b(\bm{u})>0$ and $C=C(\Omega)>0$ such that

[TABLE]

To estimate the exponential of $|\bm{u}|^{2}$, we require to estimate all $L^{q}$ norms of $\bm{u}$ for $1\leq q<\infty$. Since

[TABLE]

using Lemma 7.2, we conclude

[TABLE]

where $\hat{\bm{u}}=Q\bm{u}$ and $\hat{\bm{u}}\in H^{1}(\Omega)^{2}$. As in the proof of [53, Theorem 2] one can show that

[TABLE]

For completeness we highlight some details in the proof of (2.2j) following [53]. Using Lemma 7.3, we conclude

[TABLE]

Regarding the double integral above, applying the Cauchy-Schwarz inequality gives

[TABLE]

We estimate the two double integrals separately. Denoting by $d$ the diameter of $\Omega$, we have

[TABLE]

Hence,

[TABLE]

On the other hand, we find that

[TABLE]

and, therefore as before,

[TABLE]

Then, the second term in (2.2l) becomes

[TABLE]

Combining (2.2m)-(2.2n) and (2.2l), we deduce that

[TABLE]

Replacing the above bound in (2.2k) we obtain (2.2j).

Our aim is to bound (2.2i) with respect to $\lVert\bm{u}\rVert_{H^{1}(\Omega,T_{h})}$. From (2.2d) the following bound holds

[TABLE]

Therefore equation (2.2j) becomes

[TABLE]

Returning to (2.2i) we infer that

[TABLE]

leading to the estimate

[TABLE]

In particular, note that

[TABLE]

so that, choosing $b>0$ such that $bC\lVert\bm{u}\rVert^{2}_{H^{1}(\Omega,T_{h})}<1/e$ we reach (2.2h).

Regarding the embedding in the sense of mean convergence, as in [53], we note that bounded functions are dense in the space $\tilde{V}^{r}_{h}(\Omega)^{2}$ and hence $\tilde{V}^{r}_{h}(\Omega)^{2}\subset E_{\varphi}(\Omega)$ due to the above embedding. In particular, if $\psi(t)\leq\varphi(\lambda t)$, for some $\lambda>0$, $\tilde{V}^{r}_{h}(\Omega)\subset L_{\psi}(\Omega)$ and the embedding is continuous with respect to mean convergence as, by Remark 7.1, norm convergence implies convergence in the mean. ∎

Proposition 7.1**.**

Let $\Phi:{\mathbb{R}}^{2\times 2}\to{\mathbb{R}}$ a continuous function satisfying

[TABLE]

and suppose that

[TABLE]

for $\bm{u}_{h}\in V^{q}_{h}(\Omega)^{2}$ and $\bm{u}\in H^{2}(\Omega)^{2}$. Then, up to extracting a subsequence,

[TABLE]

Proof.

Let $\bm{w}_{h}=\nabla\bm{u}_{h}$, $\bm{w}=\nabla\bm{u}$, then $\|\bm{w}_{h}-\bm{w}\|_{L^{2}(\Omega)}+|\bm{w}_{h}-\bm{w}|_{H^{1}(\Omega,T_{h})}\to 0,\,\,h\to 0.$ Up to extracting a subsequence, $\bm{w}_{h}\to\bm{w}$ pointwise a.e. and, by the continuity of $\Phi$, also $\Phi(\bm{w}_{h})\to\Phi(\bm{w})$ a.e. in $\Omega$.

Next, note that

[TABLE]

where $\psi(t)=e^{4c_{2}|t|^{2}}-1$, is a convex function such that $\psi(t)=\varphi(\lambda t)$ for $\lambda=4c_{2}$, where $\phi$ is given in Theorem 7.1. The convexity of $\psi$ implies that

[TABLE]

By Theorem 7.1, we have that $(\psi(\bm{w}_{h}-\bm{w}))$ converges in $L^{1}(\Omega)$, as $h\rightarrow 0$. Therefore $\psi(\bm{w}_{h}/2)$ is uniformly integrable. But then (2.2q) implies that $(\Phi(\bm{w}_{h}))$ is a uniformly integrable sequence and thus $\Phi(\bm{w}_{h})\to\Phi(\bm{w})$ in $L^{1}(\Omega)$. ∎

Remark 7.2**.**

Proposition 7.1 and the classical embedding of Trudinger for Orlicz spaces [53], also implies that if $(\bm{u}_{\delta})\subset H^{2}(\Omega)$ and $\lVert\bm{u}_{\delta}-\bm{u}\rVert_{H^{2}(\Omega)}\rightarrow 0$, as $\delta\rightarrow 0$ then

[TABLE]

Finally, we establish the $\Gamma$-convergence and the convergence of discrete absolute minimizers when the penalty term $\Phi$ has exponential growth.

Theorem 7.2**.**

Let the penalty function $\Phi$ have the exponential growth (2.2b). Then, the following properties hold:

1. (i)

for all $\bm{u}\in\mathbb{A}(\Omega)$, there exists a sequence $(\bm{u}_{h})_{h>0}$ with $\bm{u}_{h}\in\mathbb{A}^{q}_{h}(\Omega)$, such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)^{2}$ and

[TABLE] 2. (ii)

for all $\bm{u}\in\mathbb{A}(\Omega)$ and all sequences $(\bm{u}_{h})\subset\mathbb{A}^{q}_{h}(\Omega)$ such that $\bm{u}_{h}\rightarrow\bm{u}$ in $H^{1}(\Omega)$ it holds

[TABLE] 3. (iii)

Theorem 6.1 holds, i.e., the discrete absolute minimizers converge to an absolute minimizer of the continuous problem.

Proof.

(i) Following the proof of Theorem 5.2 it remains only to verify $\Phi(\nabla\bm{u}_{h,\delta_{h}})\rightarrow\Phi(\nabla\bm{u})$ in $L^{1}(\Omega)$, from (2.2as), where $(\bm{u}_{h,\delta_{h}})$ is the sequence defined after (2.2aj). It is enough to show

[TABLE]

in $L^{1}(\Omega)$, as $h\rightarrow 0$. By Proposition 7.1, it suffices to show that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\nabla\bm{u}_{h,\delta_{h}}-\nabla\bm{u}_{\delta_{h}}|_{H^{1}(\Omega,T_{h})}}\rightarrow 0,\text{ as }h\rightarrow 0, i.e.,

[TABLE]

as $h\rightarrow 0$. Using (2.2am), (2.2aj) (2.2ao) we obtain the following bounds for the first term

[TABLE]

Note that $\bm{u}_{\delta}\in H^{3}(\Omega)^{2}$ implies $\bm{u}_{\delta}\in C^{1}(\Omega)^{2}$ from the Sobolev embedding. Then working as (2.2at)

[TABLE]

From (2.2v) and (2.2w) we deduce (2.2u). Also Remark 7.2 implies $\Phi(\nabla\bm{u}_{\delta_{h}})\rightarrow\Phi(\nabla\bm{u})$, in $L^{1}(\Omega)$, which concludes the proof of (i).

(ii) This is immediate as $\Phi$ is continuous and thus by Fatou’s Lemma,

[TABLE]

(iii) From the $\Gamma-$convergence result, (i) and (ii) , the proof of Theorem 6.1 can be adopted. ∎

8 Model and Computational Experiments

Numerical Energy Minimization. The potential energy of the continuous model has been discretized using the finite element discretization provided by the FEniCS project [1]. Specifically, we have used quadratic Lagrange elements, i.e., $\bm{u}_{h}\in V^{2}_{h}(\Omega)^{2}$. For the energy minimization procedure we have employed the nonlinear conjugate gradient method, see [44]. A parallelization of the nonlinear conjugate gradient algorithm has been developed using petsc4py data structures [23, 3]. The resulting deformations are visualized with paraview [51].

The nonlinear conjugate gradient method is a popular method for large-scale nonlinear optimization problems. It is a variant of the conventional conjugate gradient method, where for a quadratic function the conjugate directions and the minimum along a given direction can be computed explicitly. For the nonlinear variant, in the $k-$th iteration, given the direction $p_{k}$, the next search direction $p_{k+1}$ (conjugate direction in the conventional method) is computed through the Polar-Ribière method. The minimum value over a given direction is approximated employing line search algorithms that satisfy the strong Wolfe condition, which guarantees in certain circumstances that the computed direction $p_{k}$ is a descent direction. For details see [44].

Strain Energy Density. Here we describe our model for the strain energy density $W+\Phi$ in (1.1) by following [32, Section 5.2.1]. Starting with a (microscopic) fiber, we let the effective stretch $\lambda$ equal the distance between its endpoints divided by its undeformed (relaxed) length. The energy of a single fiber can be expressed as a function $w(\lambda)$ of the effective stretch $\lambda$. When the fiber is in tension, it is straight and $\lambda$ equals the actual stretch (strain $+1$), while $w(\lambda)$ equals the elastic energy due to stretching of the fiber. While in compression, it may be buckled; in which case the elastic (mostly bending) energy of the fiber can still be expressed as a function $w(\lambda)$ of the distance between its endpoints. In order to model a 1D fiber energy $w(\lambda)$ for a single fiber that buckles in compression, we start with the derivative $S(\lambda)=w^{\prime}(\lambda)$, which represents force as a function of stretch. We choose a polynomial that has increasing slope so that $S^{\prime}(\lambda)<S^{\prime}(1)$ for $0<\lambda<1$ (softening in compression due to buckling, and $S^{\prime}(\lambda)>S^{\prime}(1)$ for $\lambda>1$ (stiffening in tension). An example is $S(\lambda)=\mu(\lambda^{5}-\lambda^{3})$ for $\lambda>0$ with $\mu=$const.$>0$. Integrating this with respect to $\lambda$ gives an energy

[TABLE]

We model the ECM as a 2D nonlinear elastic continuum undergoing deformations $\bm{y}(\bm{x})$ where a particle with position vector $\bm{x}$ in the undeformed state is mapped to deformed position $\bm{y}=\bm{y}(\bm{x})=\bm{x}+\bm{u}(\bm{x})$, with $\bm{u}$ the displacement. The strain energy density of the material can be written as a function $\tilde{W}(\bm{F})$ of the deformation gradient $\bm{F}=\nabla\bm{y}$. We model the ECM as an isotropic material, which means that $W$ depends on $\bm{F}$ only through the principal stretches $\lambda_{1}$, $\lambda_{2}$, the eigenvalues of the right stretch tensor $(\bm{F}^{T}\bm{F})^{1/2}$. To connect the single fiber energy with the 2D strain energy density $W$ we follow e.g., [54]. We suppose the ECM consists of uniformly distributed fibers at the microscopic level. A macroscopically affine deformation is equivalent to a biaxial stretch in two orthogonal directions with stretches $\lambda_{1}$, $\lambda_{2}$, modulo rotation. A fiber of undeformed length $l$ making an angle $\theta$ with the stretch axes in the undeformed state, will have endpoints at $(0,0)$ and $(l\cos\theta,l\sin\theta)$. After deformation the latter is $(\lambda_{1}l\cos\theta,\lambda_{2}l\sin\theta)$. As a result, the effective stretch of the fiber is $\lambda=\lambda(\theta)=\sqrt{(\lambda_{1}\cos\theta)^{2}+(\lambda_{2}\sin\theta)^{2}}$, and its energy is $w(\lambda(\theta))$. Summing over all fiber orientation angles $\theta$, we obtain the macroscopic elastic energy density of the fibrous ECM:

[TABLE]

For $w(\lambda)$ given e.g. by (2.2a), this integral can be evaluated explicitly:

[TABLE]

Here the deformation invariants are $I_{1}=I_{1}(\bm{F})=\hbox{tr}(\bm{F}^{T}\bm{F})=\lambda_{1}^{2}+\lambda_{2}^{2}$ and $J=J(\bm{F})=\det\bm{F}=\lambda_{1}\lambda_{2}$, the Jacobian determinant of the deformation. The ratio of deformed to undeformed density equals $1/J$. From [30, Lemma A.1.1], $W$ satisfies the lower bound of (2.1). For an upper bound one has to remove the negative terms and use the inequality $2J\leq I_{1}$. We add a fiber volume penalty term to the energy to account for resistance of densified fibers to complete crushing by virtue of their nonzero volume, and to penalize intepenetration of matter. This term increases the energy abruptly when the Jacobian $J$ becomes less than a small positive constant $b<<1$, while it becomes negligible as $J$ increases from $b$. An example is

[TABLE]

with $\Phi(\nabla\bm{u})=\tilde{\Phi}(\det({\bm{1}}+\nabla\bm{u}))$, where $A$ is a large positive constant. The strain energy function $W+\Phi$ is non rank-one convex. In fact the total energy density expressed as a function of the principal stretch pair $U(\lambda_{1},\lambda_{2})$ is a double-well potential, modulo a null Lagrangian [32], which can be chosen so that the two minima of $U$ are of zero energy which is positive elsewhere. For a typical choice of parameters, the two minima (wells) of $U$ are the reference state $(\lambda_{1},\lambda_{2})=(1,1)$ and the state $(\lambda_{1},\lambda_{2})=(0.2,1.06)$. The latter is a severe compression in one direction ($\lambda_{1}=0.2<1$) combined with a moderate stretch ($\lambda_{2}=1.06>1$) in an orthogonal direction. This compressed well involves an almost fivefold density increase and corresponds to the densified phase consisting of buckled, collapsed fibers. In our simulations, the deformation in tethers and hairs is in the neighborhood of the compressed, densified well, while outside them it is close to the undeformed well.

Simulations. In the experiments of [32], the rather unpredictable live cells were replaced by round active particles which were embedded in the ECM. These contract on demand by a 50% decrease in radius, thereby exerting tractions onto the surrounding ECM that trigger the phase microstructure formation. In our simulations, the ECM is modelled as a homogeneous material with the strain energy function just described. In the undeformed (reference) configuration, the ECM occupies the part $\Omega$ of a rectangular or round domain exterior to one or more circles $C_{i}$ of radius $r_{c}$, which represent active-particle boundaries. Contraction of these particles is modelled by imposing Dirichlet boundary conditions on the displacement $\bm{u}$. The displacement field $\bm{u}$ is required to vanish on the outside boundary, which is thus assumed fixed. At the inner boundaries $C_{i}$, where $|\bm{x}-\bm{z}_{i}|=r_{c}$ (circles of radius $r_{c}$ and center $\bm{z}_{i}$) $\bm{u}$ is specified as

[TABLE]

This represents a radial contraction which maps each particle boundary to a smaller circle of radius $r^{\prime}_{c}=r_{c}-u_{0}$. The constant $u_{0}$ is obtained from experimental data and has a typical value of 0.5$r_{c}$. For a more sophisticated model of active particles that allows shape deviations and motion of centers due to deformation see Remark 3.1 and [32, Section 2.4, Section 5.2.2]. Computations involving a single active particle contracted by $50\%$ at the center of the ECM domain are shown in Figure 8, with a color map of the ratio of deformed over reference density. The densified phase (red) occurs in radial hairlike bands that taper off into the undensified phase (blue). Figures 1c and 1d are examples of fine phase mixtures. The mixture of low and high strain phases is energetically preferable because the average strain it produces is compatible with the Dirichlet boundary conditions, whereas the strains in the densified phase, despite their low energy, are not compatible with the boundary deformation.

In the case of two contracting active particles depicted in Figure 2, the material between the particles is stretched along the axis passing through the active particles’ centers and is compressed in the transverse direction. This renders the strain state of the densified energy well energetically favorable. As depicted in Figure 2, above some critical value of active particle contraction a tract in the densified phase, namely a tether between particles emerges, while radial bands emanate from each particle boundary as in the case of a single particle.

Examples of the agreement between simulations and experiments are shown in Figure 8. In the experiments of [32],a tether is sometimes observed to split into thinner parallel bands (Figure 5a). This is similar to the phenomenon of twinned martensite in crystals, namely splitting and tapering of twin bands in a crystal near an incompatible boundary [34]. Here as well, energy minimization forces strains to stay close to energy-density minima. The active particles in Figure 5a contracted by $u_{0}/r_{c}$ = 38%. The azimuthal stretch $\lambda_{\theta}=1-u_{0}/r_{c}=62\%$ imposed at the particle boundary by contraction is incompatible with the stretch $\lambda_{1}=20\%$ corresponding to the densified-phase energy well. To avoid this mismatch while maintaining displacement continuity, the tethers splits into narrow bands to minimize contact with the particle boundary (experiment:Figure 5a, simulation: Figure 5b).

Can our model predict how close particles should be and how much should they contract in order for a tether to form betwewen them? Particle contraction and distance between particles have been varied in multiple simulations in [32]. This provided a separatrix curve of average particle radial strain versus distance between particles (blue curve in Figure 5c). Above this curve, our model predicts that a tether forms joining the two particles; below the curve no tether will form. Data from our experiments agreed with this prediction: blue points in Figure 5c are data points from experimental particle pairs with a tether observed joining them, red points correspond to pairs without a tether between them. Furthermore, our model predicts that the correct displacement decay rate with distance from a single contracting particle, Figure 5d, which is much slower than in materials that do not suffer phase transition. This is an important feature for long-range mechanosensing. Experiments with fibroblasts [45], reveal the same displacement decay with distance as the one predicted by our computations, of the form $r^{-0.5}$. This means that the displacement fields propagate over a longer range compared to the range predicted by linear elasticity, where decay is proportional to $r^{-2}$. For more details see [30].

If the regularization (higher gradient) term is omitted, i.e. $\varepsilon=0$, then the computed solutions depend on the mesh size; similar results can be found in [42]. As depicted in Figures 1a-1d, mesh resolution must be fine enough to capture localized deformations, but further increases of resolution result in more and thinner bands around the active particle. The appearance of phase boundaries is due to the ellipticity failure [37, 48] and the rank-one connected minima. The appearance of finer and finer phase mixtures as resolution is increased is related to incompatibility of the wells with the boundary conditions[34]. As a result, the minimum is not attained, but minimizing sequences develop more and finer oscillations in order to create less incompatible deformations of lower energy [5].This is responsible for the mesh dependence, as increasing mesh resolution simply captures terms further along such minimising sequences as in Figures 1a-1d. The higher gradient term restores the ellipticity of the Euler-Lagrange equation, consequently regularizing the solution, ensuring the existence of a minimizer of limited fineness, and eliminating mesh dependence. In addition, $\varepsilon$ can be considered an internal length scale, controlling the thickness of transition layers that replace gradient discontinuities. This means that smaller values of $\varepsilon$ permit finer microstructures with more and thinner hairs; see an example with varying $\varepsilon$ in Figure 3. Figures 1e-1h illustrate that for fixed $\varepsilon>0$ numerical solutions converge as mesh size tends to zero, as expected from the analysis in this work, although it is possible that the limit state is merely an isolated local minimum of the energy in this numerical example. The analysis in this work could conceivably apply to the case of convergence to an isolated local minimum as mesh size approaches zero in the presence of fixed $\varepsilon>0$ in the sense of [13]; see also [6]. When $\varepsilon=0$, the energy functional in not lower semicontinuous, but for any $\varepsilon>0$ the presented theory holds. In the numerical minimization, the local mesh size should be smaller than the length scale imposed by $\varepsilon$ in order to resolve fine structure as is illustrated by Figures 1e-1h. Therefore, for very small values of $\varepsilon$ the degrees of freedom increase substantially due to the necessity of very high resolution, which may raise practical issues with the discrete minimization process.

More complex cases include the contraction of multiple active particles. Then the densified tehters connecting two cells can appear or disappear, influenced by other neighboring cells. An example can be seen in Figure 4.

In [32] extensive simulations of the model are performed which exhibit excellent agreement with experimental results. In particular the elastic phase transition model, combined with the numerical scheme analysed here is capable of predicting and explaining intricate details of the geometry of observed multiphase microstructures in fibrous collagen biomaterials.

Φ≥0,  and  Φ(∇u ) ≤&

Φ≥0,  and  Φ(∇u ) ≤&

Figure 5: Comparison with experimental data. (a)-(c) reproduced from [32] with permission of the authors. LABEL:sub@fig:fn3a experiment and LABEL:sub@fig:fn3b simulation of active particle pair contracted by 38%. LABEL:sub@fig:fn3a A tether fully split into multiple thin bands. Insert shows high-contrast version of area within dotted rectangle. LABEL:sub@fig:fn3b Simulation with initial radii, distances and contractions matched with LABEL:sub@fig:fn3a also results in split tether. Colorbar: ratio of deformed to undeformed density. LABEL:sub@fig:fn3c Predicting whether a tether forms between two particles. Blue curve: separatrix obtained from multiple simulations. Axes: % decrease in particle radius vs deformed distance (in deformed particle radii). Above blue curve, tethers are predicted to form between particle pairs. No tether is predicted to form below blue curve. Our experimental data (each particle pair is one point) agreed with the prediction: blue points: tether has formed. Red points: no tether has formed. LABEL:sub@fig:fn3d Displacement norm $|\bm{u}(\bm{x})$ over radial distance ($r=|\bm{x}|$). Red dots: computed displacement norm at mesh points, blue curve: graph of $Ar^{-n}$ fitted to the red dots. Least squares fit yielded parameters $A=1.43$ and $n=0.5$.

Acknowledgment

Partially supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no: 642768 ModCompShock. The research of GG was also supported by a Vannevar Bush Postdoctoral Fellowship.