---
abstract: 'In this paper we consider the branched transportation problem in 2D associated with a cost per unit length of the form $1+\alpha m$ where $m$ denotes the amount of transported mass and $\alpha>0$ is a fixed parameter (notice that the limit case $\alpha=0$ corresponds to the classical Steiner problem). Motivated by the numerical approximation of this problem, we introduce a family of functionals ($\{{\mathcal{F}}_{\varepsilon}\}_{{\varepsilon}>0}$) which approximate the above branched transport energy. We justify rigorously the approximation by establishing the equicoercivity and the $\Gamma$-convergence of $\{{\mathcal{F}}_{\varepsilon}\}$ as ${\varepsilon}\downarrow0$. Our functionals are modeled on the Ambrosio-Tortorelli functional and are easy to optimize in practice. We present numerical evidences of the efficiency of the method.'
author:
- 'A. Chambolle'
- 'B. Merlet'
- 'L. Ferrari'
bibliography:
- 'bib.bib'
title: 'A simple phase-field approximation of the Steiner problem in dimension two'
---

Introduction
============

In this paper, we introduce a phase-field approximation of a branched transportation energy for lines in the plane. Our main goal is to derive a computationally tractable approximation of the Steiner problem (of minimizing the length of lines connecting a given set of points) in a phase-field setting. Similar results have recently be obtained by [@Bon_Lem_San], however we believe our approach is slightly simpler and numerically easier to implement. We show that we can modify classical approximations for free discontinuity problems [@Mod_Mort; @Amb_Tort1; @Iur; @ContiFocardiIurlano] to address our specific problem, where the limiting energy is concentrated only on a singular one-dimensional network (and roughly measures its length). Numerical results illustrate the behaviour of these elliptic approximations. In this first study, we limit ourselves to the two-dimensional case, as in that case lines can locally be seen as discontinuities of piecewise constant functions, so that our construction is deriving in a quite simple ways from the above mentioned previous works on free discontinuity problems. Higher dimension is more challenging, from the topological point of view; an extension of this approach is currently in preparation.

We now introduce precisely our mathematical framework. Let ${\Omega}\subset {\mathbf{R}}^2$ be a convex, bounded open set. We consider measures $\s\in \M({\overline}{{\Omega}},{\mathbf{R}}^2)$ that write $$\s=\t\xi\cdot\H^1{ \!\!\mbox{{ \Large$\llcorner$}} }M,$$ where $M$ is a $1$-dimensional rectifiable set orientated by a Borel measurable mapping $\xi:M{\rightarrow}\S^1$ and $\t:M{\rightarrow}{\mathbf{R}}_+$ is a Borel measurable function representing the multiplicity. Such measure is called a rectifiable measure. We follow the notation of [@OS2011] and write $\s=U(M,\t,\xi)$. Given a cost function $f\in C({\mathbf{R}}_+,{\mathbf{R}}_+)$, we introduce the functional defined on $\M({\overline}{\Omega},{\mathbf{R}}^2)\to{\mathbf{R}}_+\cup\{+\infty\}$ as $${\mathcal{E}}_f(\s):=
\begin{dcases}
\int_Mf(\t)\dH^1&\text{if }\s=U(\t,\xi,M),\\
\quad +{\infty}&\text{in the other cases.}
\end{dcases}$$ Given a sequence of $N+1$ distinct points $S=(x_0,\dots,x_N)\in {\Omega}^{N+1}$, we consider the minimization of ${\mathcal{E}}_f(\s)$ for $\s\in \M({\overline}{\Omega},{\mathbf{R}}^2)$ satisfying the constraint $$\label{eq:divconstraint1}
{\nabla}\cdot\s=N\d_{x_0}-\sum_{i=1}^N\d_{x_i}\quad\mbox{ in }{\mathcal{D}}'({\mathbf{R}}^2).$$ The distributional support of such $\sigma$ connects the source in $x_0$ to the sinks in $x_1,\cdots,x_N$. In general, a model for branched transport connecting a set of sources to a set of sinks (represented by two discrete probabilistic measures supported on a set of points in ${\Omega}$) is obtained by choosing $f(\t)=|\t|^\a$ with $0<\a<1$ and minimizing the associated functional under a divergence constraint similar to equation . The direct numerical optimization of the functional ${\mathcal{E}}_f$ is not easy because we do not know [*a priori*]{} the topological properties of the tree $M$. For this reason it is interesting to optimize an “approximate" functional defined on more flexible objects such as functions. Such approximate model has been introduced in [@OS2011] where the authors study the $\Gamma$-convergence (see [@Bra2]) of a family of functionals inspired by the well known work of Modica and Mortola [@Mod_Mort]. Another effort in this direction can be found in the work [@Bon_Lem_San] where is studied an approximation to the Steiner Minimal Tree problem ([@Gilb_Poll], [@Ambr_Tilli] and [@Paol_Step]) by means of analogous techniques.

Here, we consider variational approximations of some energies of the form ${\mathcal{E}}_f$ through a family of functionals modeled on the Ambrosio-Tortorelli functional [@Amb_Tort1]. For being more precise, we need to introduce some material. Let $\rho\in C^\infty_c({\mathbf{R}}^2,{\mathbf{R}}_+)$ be a classical radial mollifier with $\operatorname{supp}\rho\subset B_1(0)$ and $\int \rho=1$. For ${\varepsilon}\in (0,1]$, we set $\rho_{\varepsilon}(x)={\varepsilon}^{-2}\rho({\varepsilon}^{-1}x)$ and we define the space $V_{\varepsilon}({\Omega})$ of square integrable vector fields whose weak divergence satisfy the constraint $$\label{eq:divconstraint2}
{\nabla}\cdot\s_{\varepsilon}\, =\, {\left}(N\delta_{x_0} - \sum_{j=1}^N \delta_{x_j} {\right})\ast \rho_{\varepsilon}.$$ For $\eta=\eta({\varepsilon})>0$, we note $$\label{eq:Weps}
W_{\varepsilon}({\Omega})\,=\,{\left}\{\phi\in H^1({\Omega})\,:\, \eta \leq\phi\leq 1\mbox{ in }{\Omega},\ \phi\equiv 1\mbox{ on }\de{\Omega}{\right}\}.$$ Then we define the energy ${\mathcal{F}}_{\varepsilon}:\M({\Omega},{\mathbf{R}}^2)\times L^1({\Omega}){\rightarrow}[0,+{\infty}]$ as $$\label{eq:functional}
{\mathcal{F}}_{\varepsilon}(\s,\phi):=
  \begin{dcases}
    \;\int_{{\Omega}} \frac{1}{2{\varepsilon}}\phi^2 |\s|^2 \dx+ \int_{{\Omega}} \frac{{\varepsilon}}{2}|{\nabla}\phi|^2  + \frac{(1-\phi)^2}{2{\varepsilon}} \dx, & \text{if } (\s,\phi)\in V_{\varepsilon}({\Omega})\times W_{\varepsilon}({\Omega}),\\
    \quad+{\infty}& \text{in the other cases.}
  \end{dcases}$$ The first integral in the definition of the energy will be refered to as the “constraint component” while the second integral will be regarded as the “Modica-Mortola component”. Let us briefly describe the qualitative properties of the associated minimization problem. First notice that the constraint  enforces $\s$ to be non zero on a set connecting $S$. Next, the constraint component of the energy strongly penalizes $\phi^2 |\s|^2$ so that $\phi$ should be small in the region where $|\sigma|$ is large. On the other hand the behavior of $\phi$ is controlled by the Modica Mortola component that forces $\phi$ to be close to $1$ away from a one-dimensional set as ${\varepsilon}$ converges to 0. As a consequence, we expect the support of $\sigma$ and the energy to concentrate on a one dimensional set connecting $S$. The main part of the paper consists in making rigorous and quantitative this analysis.

From now on, we assume that there exists some $\a\geq 0$ such that $$\label{eq:limalpha}
\dfrac\eta{{\varepsilon}}\ \stackrel{{\varepsilon}\downarrow0}\longrightarrow\ \a.$$ We note $\M_S({\overline}{\Omega})$ the set of ${\mathbf{R}}^2$-valued measures $\s\in \M({\mathbf{R}}^2,{\mathbf{R}}^2)$ with support in ${\overline}{{\Omega}}$ such that the constaint  holds. We define the limit energy ${\mathcal{E}}_\a:\M({\overline}{\Omega},{\mathbf{R}}^2)\times L^1({\Omega}){\rightarrow}[0,+{\infty}]$ as $$\label{eq:limit}
   {\mathcal{E}}_\a(\s,\phi)=
   \begin{dcases}
    \;\int_{M}(1+\a\,\t)\dH^1 & \text{if } \phi\equiv1,\ \s\in \M_S({\overline}{\Omega})\mbox{ and }\s=U(M,\t,\xi),\\
  ~  \quad+{\infty}&\text{in the other cases.}
   \end{dcases}$$ We prove the $\Gamma$-convergence of the sequence $({\mathcal{F}}_{\varepsilon})$ to the energy ${\mathcal{E}}_\a$ as ${\varepsilon}{\downarrow}0$. More precisely the convergence holds in $\M({\overline}{\Omega},{\mathbf{R}}^2)\times L^1({\Omega})$ where $\M({\overline}{\Omega},{\mathbf{R}}^2)$ is endowed with the weak star topology and $L^1({\Omega})$ is endowed with its classical strong topology.\
We establish the following lower bound.

\[teo:sigma\_liminf\] For any sequence $(\s_{\varepsilon},\phi_{\varepsilon})\subset\M({\Omega},{\mathbf{R}}^2)\times L^1({\Omega})$ such that $\s_{\varepsilon}\rwstar\s$ and $\phi_{\varepsilon}{\rightarrow}\phi$ in the $L^1({\Omega})$ topology, with $(\s,\phi)\in\M({\Omega},{\mathbf{R}}^2)\times L^1({\Omega})$ $$\liminf_{k{\rightarrow}+{\infty}}{\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon})\geq {\mathcal{E}}_\a(\s,\phi).$$

In this statement and throughout the paper, we make a small abuse of language by noting $(a_{\varepsilon})_{{\varepsilon}\in(0,1]}$ and calling sequence a family $\{a_{\varepsilon}\}$ labeled by a continuous parameter ${\varepsilon}\in(0,1]$. In the same spirit, we call subsequence of $(a_{\varepsilon})$, any sequence $(a_{{\varepsilon}_j})$ with ${\varepsilon}_j{\rightarrow}0$ as $j{\rightarrow}+{\infty}$.\
To complete the $\Gamma$-convergence analysis, we establish the matching upper bound.

\[teo:sigma\_limsup\] For any $(\s,\phi)\subset\M({\Omega},{\mathbf{R}}^2)\times L^1({\Omega})$ there exists a sequence $(\s_{\varepsilon},\phi_{\varepsilon})$ such that $\s_{\varepsilon}\rwstar\s$ and $\phi_{\varepsilon}{\rightarrow}\phi$ in the $L^1({\Omega})$ topology and $$\limsup_{k{\rightarrow}+{\infty}}{\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon})\leq {\mathcal{E}}_\a(\s,\phi).$$

Moreover, under the assumption $\a>0$ we prove the equicoercivity of the sequence $({\mathcal{F}}_{\varepsilon})$.

\[teo:sigma\_equicoercive\] Assume $\a>0$. For any sequence $(\s_{\varepsilon},\phi_{\varepsilon})_{{\varepsilon}\in(0,1]}\subset\M({\Omega},{\mathbf{R}}^2)\times L^1({\Omega})$ with uniformly bounded energies, [*i.e.*]{} $$\sup_{{\varepsilon}} {\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon})\ <\ +{\infty},$$ there exist a subsequence ${\varepsilon}_j{\downarrow}0$ and a measure $\s\in\M_S({\overline}{\Omega},{\mathbf{R}}^2)$ such that $\s_{{\varepsilon}_j}\rightarrow \s$ with respect to the weak-$*$ convergence of measures and $\phi_{{\varepsilon}_j} \rightarrow 1$ in $L^1({\Omega})$. Moreover, $\s$ is a rectifiable measure ([*i.e.*]{} it is of the form $\s=U(M,\t,\xi)$).

Observe that letting $\a=0$ in equation  we obtain ${\mathcal{E}}_0(\sigma)=\H^1(\{x\in M\,:\, \t(x)>0\})$ where $\s=U(M,\t,\xi)$. This is exactly the functional associated with the Steiner Minimal Tree problem. Unfortunately, the hypothesis $\a>0$ is necessary in the compactness Theorem \[teo:sigma\_equicoercive\].

The fact that we are working in dimension $2$ is fundamental for the proof of Theorem \[teo:sigma\_liminf\] as it allows to locally rewrite the vector field $\s_{\varepsilon}$ as the rotated gradient of a function.

**Structure of the paper:** In Section \[sec:Not\] we introduce some notation and several tools and notions on $SBV$ functions and vector field measures. In Section \[sec:LocRes\] we study a first family of energies obtained by substituting ${\nabla}u$ for $\s$ in the definition of ${\mathcal{F}}\e$. In Section \[sec:Comp\] we prove the equicoercivity result, Theorem \[teo:sigma\_equicoercive\] and we establish the lower bound stated in Theorem \[teo:sigma\_liminf\]. In Section \[sec:UppBound\] we prove the upper bound of Theorem \[teo:sigma\_limsup\]. Finally, in the last section, we present and discuss various numerical simulations.

Notation and Preliminary Results {#sec:Not}
================================

In the following ${\Omega}\subset\subset\hat{\Omega}\subset{\mathbf{R}}^d$ are bounded open convex sets. Given $X\subset{\mathbf{R}}^d$ (in practice $X={\Omega}$ or $X=\hat{\Omega})$, we denote by $\A(X)$ the class of all open subsets of $X$ and by $\A_S(X)$ the subclass of all simply connected open sets $O\subset X$ such that $\overline{O}\cap S=\varnothing$. We denote by $(e_1,\dots,e_d)$ the canonical orthonormal basis of ${\mathbf{R}}^d$, by $|\cdot|$ the euclidean norm and by $\langle\cdot,\cdot\rangle$) the euclidean scalar product in ${\mathbf{R}}^d$. The open ball of radius $r$ centered at $x\in{\mathbf{R}}^d$ is denoted by $B_r(x)$. The $(d-1)$-dimensional Hausdorff measure in ${\mathbf{R}}^d$ is denoted by $\H^{d-1}$. We write $|E|$ to denote the Lebesgue measure of a measurable set $E\subset {\mathbf{R}}^d$. When $\mu$ is a Borel meaure and $E\subset {\mathbf{R}}^d$ is a Borel set, we denote by $\mu{ \!\!\mbox{{ \Large$\llcorner$}} }E$ the measure defined as $\mu{ \!\!\mbox{{ \Large$\llcorner$}} }E(F)=\mu(E\cap F)$.\
Let us remark that from Section \[sec:Comp\] onwards, we work in dimension $d=2$.\
For $O\in\A({\Omega})$, the functional ${\mathcal{F}}_{\varepsilon}(\cdot,\cdot;O)$ is the functional obtained by substuting $O$ for ${\Omega}$ in the definition of ${\mathcal{F}}_{\varepsilon}$ (see ). Similarly we define the local version ${\mathcal{E}}_\a(\cdot,\cdot;O)$ of ${\mathcal{E}}_\a$.

BV(${\Omega}$) functions and Slicing
------------------------------------

BV(${\Omega}$) is the space of functions $u\in L^1({\Omega})$ having as distributional derivative $Du$ a measure with finite total variation. For $u\in BV({\Omega})$, we denote by $S_u$ the complement of the Lebesgue set of $u$, that is $x\not \in S_u$ if and only if $$\lim_{\rho{\rightarrow}0^+}\frac{1}{|B_\rho(x)|}\int_{B_\rho(x)}|u(y)-z|\dy=0$$ for some $z\in{\mathbf{R}}$. We say that $x$ is an approximate jump point of $u$ if there exist $\xi \in \S^{d-1}$ and distinct points $a,b \in {\mathbf{R}}$ such that $$\lim_{\rho{\downarrow}0}\frac{1}{|B^+_\rho(x,\xi)|}\int_{B^+_\rho(x,\xi)}|u(y)-a|\dy=0\qquad\text{and}\qquad\lim_{\rho{\downarrow}0}\frac{1}{|B^-_\rho(x,\xi)|}\int_{B^+_\rho(x,\xi)}|u(y)-b|\dy=0,$$ where $B^\pm_\rho(x,\xi):=\{y\in B_\rho(x)\;:\;\pm\langle y-x,\xi\rangle\geq0\}$. Up to a permutation of $a$ and $b$ and a change of sign of $\xi$, this characterizes the triplet $(a,b,\xi)$ which is then denoted by $(u^+,u^-,\nu_u)$. The set of approximate jump points is denoted by $J_u$. The following theorem holds [@Am_Fu_Pal].

\[teo:federer\] The set $S_u$ is countably $\H^{d-1}$-rectifiable and $\H^{d-1}(S_u\setminus J_u)=0$. Moreover $Du{ \!\!\mbox{{ \Large$\llcorner$}} }J_u=(u^+-u^-)\nu_u\H^{d-1}{ \!\!\mbox{{ \Large$\llcorner$}} }J_u$ and $$Tan^{d-1}(J_u,x)=\nu_u(x)^\perp$$ for $\H^{d-1}$-a.e. $x\in J_u$.

We write the Radon-Nikodym decomposition of $Du$ as $Du= {\nabla}u \dx+D^s u$. Setting $D^c u=D^s u { \!\!\mbox{{ \Large$\llcorner$}} }({\Omega}\setminus S_u)$ we get the decomposition $$Du={\nabla}u \dx+ (u^+-u^-)\nu_u\H^{d-1}{ \!\!\mbox{{ \Large$\llcorner$}} }J_u+D^c u.$$ Moreover the Cantor part is such that if $\H^{d-1}(E)<\infty$, we have $|D^cu|(E)=0$ [@Am_Fu_Pal Thm. 3.92]. In particular, we have the following useful consequence $$\label{BVproperty}
\H^{d-1}(E)=0\quad\Longrightarrow\quad |Du|(E)=0.$$ We frequently use the notation $[u]$ for the jump function $(u^+-u^-):J_u{\rightarrow}{\mathbf{R}}$.\
When $d=1$ we use the symbol $u'$ in place of ${\nabla}u$ and $u(x^\pm)$ to indicate the right and left limits of $u$ at $x$.

Let us introduce the space of special functions of bounded variation and a variant: $$SBV({\Omega}):=\{u\in BV({\Omega}):D^cu=0\},$$ $$GSBV({\Omega}):=\{u\in L^1({\Omega}):\max(-T,\min(u,T))\in SBV({\Omega})\;\forall T\in{\mathbf{R}}\}.$$ Eventually, in Section \[sec:LocRes\], the following space of piecewise constant functions will be useful. $$\P({\Omega})=\{u\in GSBV({\Omega}):{\nabla}u=0\}.$$ To conclude this section we recall the slicing method for functions with bounded variation. Let $\xi\in\S^{d-1}$ and let $$\Pi_\xi:=\{y\in{\mathbf{R}}^d:\langle y,\xi\rangle=0\}.$$ If $y\in\Pi_\xi$ and $E\in{\mathbf{R}}^d$, we define the one dimensional slice $$E_{\xi,y}:=\{t\in{\mathbf{R}}: y+t\xi\in E\}.$$ For $u:{\Omega}{\rightarrow}{\mathbf{R}}$, we define $u_{\xi,y}:{\Omega}_{\xi,y}{\rightarrow}{\mathbf{R}}$ as $$u_{\xi,y}(t):=u(y+t\xi), \quad t\in {\Omega}_{\xi,y}.$$ Functions in $GSBV({\Omega})$ can be characterized by one-dimensional slices (see [@Bra1 Thm. 4.1])

\[teo:braides\] Let $u\in GSBV({\Omega})$. Then for all $\xi \in\S^{d-1}$ we have $$u_{\xi,y}\in GSBV({\Omega}_{\xi,y})\quad \text{for } \H^{d-1}\text{-a.e. }y\in\Pi_\xi.$$ Moreover for such $y$, we have $$\begin{gathered}
 u'_{\xi,y}(t)=\langle{\nabla}u(y+t\xi),\xi\rangle\quad \text{for a.e. } t\in{\Omega}_{\xi,y},\\
 J_{u_{\xi,y}}=\{t\in {\mathbf{R}}:y+t\xi\in J_u\},\end{gathered}$$ and $$u_{\xi,y}(t^\pm)=u^\pm(y+t\xi)\quad \text{or}\quad u_{\xi,y}(t^\pm)=u^\mp(y+t\xi)$$ according to whether $\langle\nu_u,\xi\rangle>0$ or $\langle\nu_u,\xi\rangle<0$. Finally, for every Borel function $g:{\Omega}{\rightarrow}{\mathbf{R}}$, $$\label{eq:braidesricostruzione}
 \int_{\Pi_\xi}\sum_{t\in J_{u_{\xi,y}}} g_{\xi,y}(t)\dH^{d-1}(y)=\int_{J_u}g|\langle\nu_u,\xi\rangle|\dH^{d-1}.$$ Conversely if $u\in L^1({\Omega})$ and if for all $\xi \in\{e_1,\dots,e_d\}$ and almost every $y\in\Pi_\xi$ we have $u_{\xi,y}\in SBV({\Omega}_{\xi,y})$ and $$\int_{\Pi_\xi}|Du_{\xi,y}|({\Omega}_{\xi,y})\dH^{d-1}(y)<+{\infty}$$ then $u\in SBV({\Omega})$.

Rectifiable vector Measures
---------------------------

Let us introduce the linear operator $\perp$ that associates to each vector $v=(v_1,v_2)\in {\mathbf{R}}^2$ the vector $v^\perp=(-v_2,v_1)$ obtained via a $90^\circ$ counterclockwise rotation of $v$. Notice that the $\perp$ operator maps divergence free ${\mathbf{R}}^2$-valued measures onto curl free ${\mathbf{R}}^2$-valued measures. Let $O\subset{\mathbf{R}}^2$ be a simply connected and bounded open set. By Stokes Theorem, for any divergence free measure $\s\in\M(O,{\mathbf{R}}^2)$ there exists a function $u\in BV({\Omega})$ with zero mean value such that $\s= Du^\perp$. On the other hand for $u\in \P({\Omega})$ $\sigma:=Du^\perp$ is divergence free and by Theorem \[teo:federer\], $\sigma =(u^+-u^-)\nu_u^\perp \H^1=U(J_u,[u],\nu_u^\perp)$.

Let us now produce an elementary example of measure $\g$ of the form $U(M,\t,\xi)$.

\[ex:div2punti\] Given two points $x,y\in{\Omega}$ we consider the smooth path from $x$ to $y$ defined as $$r(t):=x+t\frac{y-x}{|x-y|}\quad \mbox{ for } t\in [0,|x-y|].$$ We define the measure $\g\in \M({\Omega},{\mathbf{R}}^2)$ by $$(\phi,\g)\,:=\,\int\langle\phi(r(t)),\dot{r}(t)\rangle\dt\qquad\mbox{for any } \phi\in \C(\overline{{\Omega}},{\mathbf{R}}^2).$$ We then have $\g=U([x,y],1,\xi)$ with $\xi=\frac{y-x}{|y-x|}$. Notice that ${\nabla}\cdot\g=\d_{x}-\d_y$.\
Similarly we obtain a measure satisfying this property by substituting for $r$ any Lipschitz path from $x$ to $y$.

We will make use of the following construction.

\[lemma:completediv\] Given $S=(x_0,\cdots,x_N)\in {\Omega}^{N+1}$ a sequence of $N+1$ distinct points, there exist a vector measure $\g= U(M_\g,\t_\g,\xi_\g)$ and a finite partition $({\Omega}_i)\subset \A({\Omega})$ of ${\Omega}$ such that

1.  ${\nabla}\cdot\g=-N\d_{x_0}+\sum_{i=1}^{N} \d_{x_i}$,

2.  $\t_\g:M_\g{\rightarrow}\{1,N\}$,

3.  each ${\Omega}_i$ is a polyhedron,

4.  $ M_\g\subset \bigcup_i \de{\Omega}_i$,

5.  ${\Omega}_i$ is of finite perimeter for each $i$ and ${\Omega}_i\cap{\Omega}_j=\emptyset$ for $i\neq j$,

6.  $\L^2({\Omega}\sm\cup_i{\Omega}_i)=0$.

Moreover if $M$ is a $1$ dimensional countably rectifiable set, we can choose $\gamma$ and $({\Omega}_i)$ such that $\H^1(M\cap \bigcup_i \de{\Omega}_i)=0$.

Let us fix a point $p\in{\Omega}\setminus S$. By Example \[ex:div2punti\] we can construct a measure $\g_i$ with ${\nabla}\cdot\g_i=\d_{x_i}-\d_{p}$ for $i\in\{0,\cdots,N\}$. We define $$\g=-N \g_0 +\sum_{i=1}^N\g_i.$$ By construction $(a)$ holds true. Moreover, up to a small displacement of $p$ we may assume that $[p,x_i]\cap[p,x_j]=\{p\}$ for $i\neq j$ so that $(b)$ holds.\
Next, let $D_j$ be the straight line supporting $[p,x_j]$. We define the sets $({\Omega}_i)$ as the connected components of ${\Omega}\setminus \left(D_0 \cup\cdots\cup D_N\right)$. We see that $(c,d,e,f)$ hold true.\
For the last statement, we observe that by the coarea formula, we have $\H^1(M\cap \bigcup_i \de{\Omega}_i)=0$ for a.e. choice of $p$.

![Example of the construction of the $\H^1$-rectifiable measure $\g$ (red) and of the partition $\{{\Omega}_i\}$ (gray) in the case $\s$ (green) is being an $\H^1$-rectifiable vector measure.[]{data-label="fig:partizione"}](partizione){width=".4\textwidth"}

Local Result {#sec:LocRes}
============

In this section we introduce a localization of the family of functionals $({\mathcal{F}}_{\varepsilon})$ (see ). We establish a lower bound and a compactness property for these local energies.

\[procedure\] Let $O\in A_S({\Omega})$, for $u{\varepsilon}\in H^1(O)$ and $\phi_{\varepsilon}\in H^1(O)$, we define $$\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O)\ :=\ {\mathcal{F}}_{\varepsilon}({\nabla}u_{\varepsilon},\phi_{\varepsilon};O).$$ Notice that for ${\varepsilon}<d(O,S)$, we have ${\nabla}\cdot \s_{\varepsilon}\equiv 0$ in $O$ for any $\s_{\varepsilon}\in V_{\varepsilon}({\Omega})$. By the Stokes theorem we have $Du_{\varepsilon}=\s_{\varepsilon}^\perp$ for some $u\in H^1(O)$ and we have $${\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon};O)=\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O).$$

The remaining of the section is devoted to the proof of

\[teo:localResult\] Let $(u_{\varepsilon})_{{\varepsilon}\in(0,1]}\subset H^1(O)$ be a family of functions with zero mean value and let $(\phi_{\varepsilon})\subset H^1(O)$ such that $\phi_{\varepsilon}\in H^1(O,[\eta({\varepsilon}),1])$. Assume that $c_0:=\sup_{\varepsilon}\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O)$ is finite, then there exist a subsequence ${\varepsilon}_j$ and a function $u\in BV({\Omega})$ such that

a)  $\phi_{{\varepsilon}_j} {\rightarrow}1$ in $L^2(O)$,

b)  $u_{{\varepsilon}_j}{\rightarrow}u$ with respect to the weak-$*$ convergence in $BV$,

c)  $u\in \P(O)$.

Furthermore for any $u\in\P(O)$ and any sequence $(u_{\varepsilon},\phi_{\varepsilon})$ as above such that $u_{\varepsilon}\rwstar u$, we have the following lower bound of the energy: $$\liminf_{{\varepsilon}{\rightarrow}0}\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O)\geq \int_{J_u\cap O}{\left}[1+\alpha|[u]|{\right}]\dH^{d-1}.$$

The proof is achieved in several steps and mostly follows ideas from [@Iur] (see also [@ContiFocardiIurlano]). In the first step we obtain *(a)* and *(b)*. In steps 2. we prove *(c)* and the lower bound for one dimensional slices. Finally in step 3. we prove *(c)* and the lower bound in dimension $d$. The construction of a recovery sequence that would complete the $\Gamma$-limit analysis is postponed to the global model in Section \[sec:UppBound\].

*Step 1.* Item *(a)* is a straightforward consequence of the definition of the functional. Indeed, we have $$\int_O(1-\phi_{\varepsilon})^2\, \dx\, \leq {\varepsilon}\, \LF_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon})\, \leq\, c_0\, {\varepsilon}\,  \stackrel{{\varepsilon}{\downarrow}0}\longrightarrow \ 0.$$ For *(b)*, since $(u_{\varepsilon})$ has zero mean value, we only need to show that $ \sup_{\varepsilon}\{|Du_{\varepsilon}|(O): k\in{\mathbf{N}}\}<+\infty$. Using Cauchy-Schwarz inequality we get $$\label{eq:compactness1}
{\left}[|Du_{\varepsilon}|(O){\right}]^2={\left}(\int_{O}|{\nabla}u_{\varepsilon}|{\right})^2\leq {\left}({\varepsilon}\int_{O}\frac{1}{\phi^2\e}{\right}){\left}(\frac{1}{{\varepsilon}}\int_{O}\phi^2\e|{\nabla}u_{\varepsilon}|^2{\right}).$$ By assumption, the second therm in the right hand side of is bounded by $2c_0$. In order to estimate the first term we split $O$ in the two sets $\{\phi_{\varepsilon}<1/2\}$ and $\{\phi_{\varepsilon}\geq1/2\}$. We have,$${\varepsilon}\int_{O}\frac{1}{\phi\e}\,=\, {\varepsilon}\int_{\{\phi_{\varepsilon}<1/2\}}\frac{1}{\phi^2\e}+{\varepsilon}\int_{\{\phi\e\geq1/2\}}\frac{1}{\phi^2\e}.$$ Since $\phi_{\varepsilon}\geq \eta$ and $(1-t)^2$ is decreasing in the interval $O$ the following inequalities hold $$\begin{aligned}
\int_{\{\phi\e<1/2\}}\frac{1}{\phi^2\e}&\leq& \frac{2 {\varepsilon}}{\eta^2(1-1/2)^2}\int_{O}\frac{(1-\phi\e)^2}{2{\varepsilon}}\,\leq\,\frac {8{\varepsilon}}\eta^2 c_0,\\
\nonumber\int_{\{\phi_{\varepsilon}\geq1/2\}}\frac{1}{\phi^2\e}&\leq &\int_{O}\frac{1}{(1/2)^2} \, =\, 4|O|.\end{aligned}$$ Combining these estimates with  we obtain $$\label{eq:limitatezzagradu}
{\left}[|Du_{\varepsilon}|(O){\right}]^2\, \leq\, \frac{{\varepsilon}^2}{\eta^2}\:16c_0^2 + 8{\varepsilon}|O| c_0\,\stackrel{{\varepsilon}{\downarrow}0}\longrightarrow\,  \frac{16c_0^2}{\a^2}\,<\,\infty.$$ This establishes *(b)*.

*Step 2.* In this step we suppose $O$ to be an interval of ${\mathbf{R}}$. We first prove that $u$ is piecewise constant. The idea is that in view of the constraint component of the energy, variations of $u_{\varepsilon}$ are balanced by low values of $\phi_{\varepsilon}$. On the other hand the Modica-Mortola component of the energy implies that $\phi_{\varepsilon}\simeq 1$ in most of the domain and that transitions from $\phi_{\varepsilon}\simeq 1$ to $\phi_{\varepsilon}\simeq 0$ have a constant positive cost (and therefore can occur only finitely many times).

*Step 2.1. proof of $u\in \P(O)$.* Let us define $$\label{def:insiemi}
B_{\varepsilon}:= {\left}\{x\in O:\phi_{\varepsilon}(x)<\frac{3}{4}{\right}\}\,\supset\,  A_{\varepsilon}:= {\left}\{x\in O:\phi_{\varepsilon}(x)<\frac{1}{2}{\right}\},$$ and let $$\label{def:C}
C_{\varepsilon}= \{I \text{ connected component of } B_{\varepsilon}: I\cap A_{\varepsilon}\neq \varnothing \}.$$ Let us show that the cardinality of $C_{\varepsilon}$ is bounded by a constant independent of ${\varepsilon}$. Let ${\varepsilon}$ be fixed and consider an interval $I\in C_{\varepsilon}$. Let $a,b\in \bar I$ such that $\{\phi_{\varepsilon}(a),\phi_{\varepsilon}(b)\}=\{1/2,3/4\}$. Using the usual Modica-Mortola trick, we have $$\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};I)\,\geq\,\int_I {\varepsilon}|\phi'_{\varepsilon}|^2  + \frac{(1-\phi_{\varepsilon})^2}{4{\varepsilon}} \dx \geq \int_{(a,b)} |\phi'_{\varepsilon}|(1-\phi_{\varepsilon})\dx\,\geq\, \int_{1/2}^{3/4} (1-t)\,dt \, =\, \dfrac3{2^5}.$$ Since all the elements of $C_{\varepsilon}$ are disjoint, we deduce from the energy bound that $$\# C_{\varepsilon}\leq2^5c_0/3,$$ where we note $\# C_{\varepsilon}$ the cardinality of $C_{\varepsilon}$. Next, up to extraction we can assume that $\# C_{\varepsilon}=N$ is fixed. The elements of $ C_{\varepsilon}$ are written on the form $I^{\varepsilon}_i=(m^{\varepsilon}_i-w_i^{\varepsilon},m^{\varepsilon}_i+w_i^{\varepsilon})$ for $i=1,\cdots,N$ with $m^{\varepsilon}_i<m^{\varepsilon}_{i+1}$. Since $\phi_{\varepsilon}{\rightarrow}1$ in $L^1(O)$ we have $$\label{misura:zero}
\sum_{I^{\varepsilon}_i\in C_{\varepsilon}}|I^{\varepsilon}_i|= \sum_i2w^{\varepsilon}_i {\rightarrow}0.$$ Up to further extraction, we can assume that each sequence $(m_i^{\varepsilon})$ converges in ${\overline}O$. We call $m_1\leq m_2 \leq \dots \leq m_N$ their limits. We now prove that $|D u|(O\setminus \{m_i\}_{i=1}^N) =0$. For this, we fix $x\in O\setminus\{m_i\}_{i=0}^N$ and establish the existence of a neighborhood $B_\delta(x)$ of $x$ for which $|Du|(B_\delta(x))=0$. Let $0<\delta\leq 1/2\min_i|x-{m}_i|$. Equation  ensures that for ${\varepsilon}$ small enough $B_\delta(x)\cap C_{\varepsilon}=\varnothing$. Notice that from the definitions in  and  we have that $\phi_{\varepsilon}\geq 1/2$ outside $C_{\varepsilon}$. Hence, using Cauchy-Scwarz inequality, we have for ${\varepsilon}$ small enough, $${\left}(\int_{B_\delta(x)}| u'_{\varepsilon}|\dx{\right})^2\leq 2\delta \int_{B_\delta(x)}| u'_{\varepsilon}|^2\dx\leq (2\delta)(2{\varepsilon}) 4{\left}(\frac{1}{2{\varepsilon}}\int_{B_\delta(x)}\phi_{\varepsilon}|^2 u'_{\varepsilon}|^2\dx{\right})\, \leq\, 16c_0{\varepsilon}\delta \,\stackrel{{\varepsilon}{\downarrow}0}\longrightarrow\,0.$$ By lower semicontinuity of the total variation on open sets we conclude that $|Du|(B_\delta(x))=0$. This establishes $u\in \P(O)$ with $J_u\subset \{m_1,\cdots,m_N\}$.

*Step 2.2. Proof of the lower bound.* Without loss of generality, we prove the lower bound in the case $J_u=\{0\}$ and $ D:=u(0^+)=-u(0^-)>0$. Using the same argument as in [@Iur Pag. 7] for any $0<d<D$ there exist six points $y_1<x^1_{\varepsilon}\leq\tilde{x}^1_{\varepsilon}<\tilde{x}^2_{\varepsilon}\leq x^2_{\varepsilon}<y_2$ such that $$\begin{gathered}
\lim_{{\varepsilon}{\rightarrow}0}\phi_{\varepsilon}(y_1)= \lim_{{\varepsilon}{\rightarrow}0} \phi_{\varepsilon}(y_2) = 1,\\
\lim_{{\varepsilon}{\rightarrow}0}\phi_{\varepsilon}(x^1_{\varepsilon})= \lim_{{\varepsilon}{\rightarrow}0} \phi_{\varepsilon}(x^2_{\varepsilon}) =0,\\
u_{\varepsilon}(\tilde{x}^1_{\varepsilon})= -D+d,\quad\quad \; u_{\varepsilon}(\tilde{x}^2_{\varepsilon})= D-d.\end{gathered}$$ Using the Modica-Mortola trick in the intervals $(y_1,x^1_{\varepsilon})$ and $(x^2_{\varepsilon},y_2)$ as above, we compute: $$\label{dis1:liminf} 
\liminf_{{\varepsilon}{\downarrow}0}\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};(y_1,x^1_{\varepsilon})\cup(x^2_{\varepsilon},y_2))\geq\liminf_{{\varepsilon}{\downarrow}0}\int^{x^1_{\varepsilon}}_{y_1}(1-\phi_{\varepsilon})|\phi'_{\varepsilon}| \dx+\int_{x^2_{\varepsilon}}^{y_2}(1-\phi_{\varepsilon})|\phi'_{\varepsilon}| \dx\,\geq\,1.$$ For the estimation on the interval $I_{\varepsilon}=(\txk^1,\txk^2)$ let us introduce: $$\begin{gathered}
G_{\varepsilon}:= {\left}\{ w\in H^1(I_{\varepsilon}) : w(\txk^1)=-D+d,\, w(\txk^2)=D-d {\right}\},\\
Z_{\varepsilon}:= {\left}\{ z\in H^1(I_{\varepsilon}) :\eta\leq z\leq1 \text{ a.e. on } I_{\varepsilon}{\right}\},\\
H_{\varepsilon}(w,z):=\int_{I_{\varepsilon}}{\left}(\frac{1}{2{\varepsilon}}z^2|w'|^2+\frac{(1-z)^2}{2{\varepsilon}}{\right})\dx,\\
h_{\varepsilon}(z)= \inf_{w\in W_{\varepsilon}}H_{\varepsilon}(w,z) \text{ for } z\in Z_{\varepsilon}.\end{gathered}$$ Note that by Cauchy-Schwarz inequality, we have for $w\in G_{\varepsilon}$ and $z\in Z_{\varepsilon}$, $$\int_{I_{\varepsilon}}z^2|w'|^2\, \geq \, \left(\int_{I_{\varepsilon}}|w'|\dx\right)^2{\left}(\int_{I_{\varepsilon}}\frac{1}{z^2}{\right})^{-1}\, \geq\, 4(D-d)^2{\left}(\int_{I_{\varepsilon}}\frac{1}{z^2}{\right})^{-1}.$$ We deduce the lower bound $$\label{eq:minimo}
h_{\varepsilon}(z) \, \geq  \, 4(D-d)^2 {\left}(2{\varepsilon}\int_{I_{\varepsilon}}\frac{1}{z^2}\dx{\right})^{-1}+ \int_{I_{\varepsilon}}{\left}(\frac{(1-z)^2}{2{\varepsilon}}{\right})\dx.$$ Let us remark that optimizing $H_{\varepsilon}(w,z)$ with respect to $w\in G_{\varepsilon}$ we see that this inequality is actually an equality.\
Consider for $0 <\lambda<1$ the inequalities: $$\int_{\{x\in I_{\varepsilon}: \phi_{\varepsilon}\geq \lambda\}}\frac{1}{\phi^2_{\varepsilon}}\leq \frac{\L^1(I_{\varepsilon})}{\lambda^2}\qquad\mbox{ and }\qquad\int_{\{x\in I_{\varepsilon}: \phi_{\varepsilon}< \lambda\}}\frac{1}{\phi^2_{\varepsilon}}\leq \frac{1}{(1-\lambda)^2}\frac{2{\varepsilon}}{\eta^2}{\left}(\int_{I_{\varepsilon}}\frac{(1-\phi_{\varepsilon})^2}{2{\varepsilon}}\dx{\right}).$$ Applying both of them in  we obtain $$\begin{aligned}
\nonumber
\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon},I_{\varepsilon})&\geq  h_{\varepsilon}(\phi_{\varepsilon})\\ \nonumber
&\geq\frac{2(D-d)^2}{\frac{{\varepsilon}\L^1(I_{\varepsilon})}{\lambda^2}+\frac{1}{(1-\lambda)^2}\frac{2{\varepsilon}^2}{\eta^2}{\left}(\int_{I\e}\frac{(1-\phi\e)^2}{2{\varepsilon}}\dx{\right})}+ \int_{I\e}{\left}(\frac{(1-\phi\e)^2}{2{\varepsilon}}{\right})\dx \\
&\geq 2(1-\lambda)\frac{\eta}{{\varepsilon}}(D-d)-(1-\lambda)^2\frac{\eta^2}{2{\varepsilon}}\frac{\L^1(I\e)}{\lambda^2} \label{dis2:liminf}\end{aligned}$$ where the latter is obtained by minimizing the function: $$t \mapsto \frac{2(D-d)^2}{\frac{{\varepsilon}\L^1(I_{\varepsilon})}{\lambda^2}+\frac{1}{(1-\lambda)^2}\frac{2{\varepsilon}^2}{\eta^2}t}+ t.$$ Therefore we can pass to the limit in  and obtain: $$\liminf_{{\varepsilon}{\downarrow}0}\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon},I_{\varepsilon})\geq(1-\lambda)\a\,2(D-d).$$ Sending $\lambda$ and $d$ to $0$ and recalling the estimation in  we get $$\label{eq:liminf1d}
\liminf_{{\varepsilon}{\downarrow}0}\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon},(y_1,y_2))\, \geq\, 1+\a\,2D \, =\, 1+\a|u(0^+)-u(0^-)|.$$ *Step 3.* Using Fubini’s decomposition we can rewrite the energy obtaining for every $\xi\in \S^{d-1}$, $$ \LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O)\,\geq\,\int_{\Pi_\xi}{\left}( \int_{O^\xi_y}\frac{1}{2{\varepsilon}}(\phi^2_{\varepsilon})^\xi_y |(u'_{\varepsilon})^\xi_y|^2  + \frac{{\varepsilon}}{2}|(\phi'_{\varepsilon})^\xi_y |^2  + \frac{(1-(\phi_{\varepsilon})^\xi_y )^2}{2{\varepsilon}} \dt{\right})\dH^{d-1}(y).$$ Since $\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O)$ is bounded, for $\H^{d-1}$ almost every $y\in \Pi_\xi$, the inner integral in the latter is also bounded, furthermore it corresponds to the functional on the one dimensional slice $O^\xi_y$ studied in the previous step evaluated on the couple $((u_{\varepsilon})^\xi_y,(\phi_{\varepsilon})^\xi_y)$. Taking the $\liminf$ of the above quantity, using , we get by Fatou’s lemma that for any $\xi\in \S^{d-1}$ and $\H^{d-1}$ almost every $y\in {\Omega}_\xi$ $$\int_{\Pi_\xi}\sum_{m_i\in(J_{u})^\xi_y}{\left}[1+\a|u^\xi_y(m_i^+)-u^\xi_y(m_i^-)| {\right}]\mathrm{d}\H^{d-1}(y)\leq \liminf_{{\varepsilon}{\downarrow}0} \LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O).$$ Therefore in force of Theorem \[teo:braides\] we have $u\in SBV(O)$. Moreover, since $(u')^\xi_y=0$ on each slice, we have $u\in \P(O)$. Applying identity  we get $$\label{eq:liminfxi} 
\liminf_{{\varepsilon}{\rightarrow}0}\LF_{\varepsilon}(u_{\varepsilon}, \phi_{\varepsilon};O) \geq \int_{J_u\cap O}|\nu_u\cdot \xi|{\left}[1+\a|[u]|{\right}]\dH^{d-1}.$$ In order to conclude, we use the following localization method stated by Braides in [@Bra1 Prop. 1.16].

\[lem:braides\] Let $\mu:\A(X){\rightarrow}[0,+{\infty})$ be a superadditive set function and let $\lambda$ be a positive measure on $X$. For any $i\in{\mathbf{N}}$ let $\psi_i$ be a Borel function on $X$ such that $\mu(A)\geq\int_A \psi_i\dif \lambda$ for all $A\in\A(X)$. Then $$\mu(A)\geq\int_A\psi\dif\lambda$$ where $\psi:=\sup_i\psi_i$.

We introduce the superadditive increasing set function $\mu$ defined on $\A(O)$ by $$\mu(A):=\Gamma-\liminf_{{\varepsilon}{\rightarrow}0} \LF_{\varepsilon}(u;)),\quad\text{for any }A\in\A(O)$$ and we let $\lambda$ be a Radon measure defined as $$\lambda:=[1+\a|u(x^+)-u(x^-)|]\H^{d-1}{ \!\!\mbox{{ \Large$\llcorner$}} }J_u.$$ Fix a sequence $(\xi_i)_{i\in{\mathbf{N}}}$ dense in $\S^{d-1}$. By  we have $$\mu(O)\geq \int_O \psi_i\dif \lambda, \quad i\in {\mathbf{N}},$$ where $$\psi_i(x):=
\begin{dcases}
   |\langle\nu_u(x),\xi_i\rangle|&\text{if } x\in J_u,\\
   0&\text{if } x\in O\setminus J_u.
  \end{dcases}$$ Hence by Lemma \[lem:braides\] we finally obtain $$\liminf_{{\varepsilon}{\rightarrow}0}\LF_{\varepsilon}(u_{\varepsilon},\phi_{\varepsilon};O)\geq \int_O \sup_i\psi_i(x)\dif \mu =\int_{J_u\cap O}{\left}[1+\a|[u]|{\right}]\dH^{d-1}.$$

Equicoercivity and $\Gamma$-liminf {#sec:Comp}
==================================

We first prove the compactness property stated in the introduction. Let us consider a sequence $(\s_{\varepsilon},\phi_{\varepsilon})\in\M({\Omega},{\mathbf{R}}^2)$ uniformly bounded in energy by $c_0<+{\infty}$, $$\label{eq:unifbound}
0\leq{\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon})\leq c_0\qquad\mbox{for } {\varepsilon}\in(0,1].$$

First observe that by definition  and equation , we have $\s_{\varepsilon}\in V_{\varepsilon}({\Omega})$ and $\phi_{\varepsilon}\in W_{\varepsilon}({\Omega})$.

Next, substituting $|\s_{\varepsilon}|$ for $|{\nabla}u_{\varepsilon}|$ in the argument of *Step 1.* of the proof of Theorem \[teo:localResult\], inequality  reads $$ |\s_{\varepsilon}|({\Omega})\leq\sqrt{16\frac{{\varepsilon}^2}{\eta^2}\:c_0^2 + 8{\varepsilon}|{\Omega}| c_0} \,\stackrel{{\varepsilon}{\downarrow}0}\longrightarrow\,  \frac{4 c_0}{\a}\,<\,\infty.$$ Thus the total variation of $(\s_{\varepsilon})$ is uniformly bounded and there exists $\s\in \M_S({\overline}{{\Omega}})$ such that up to extraction $\s_{\varepsilon}\to\s$ weakly-$*$ in $\M({\overline}{\Omega})$.

Now, considering the last term in the energy  we have $$\int_{\Omega}(1-\phi_{\varepsilon})^2 \dx\leq 2{\varepsilon}\;{\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon})\leq 2{\varepsilon}\;c_0 \rightarrow 0.$$ Hence, $\phi_{\varepsilon}\to1$ in $L^2({\Omega})$.

Let us now study the structure of the limit measure $\s$. Let us recall that $\hat{\Omega}$ is a bounded convex open set such that ${\overline}{\Omega}\subset\hat{\Omega}$ and let us extend $\s_{\varepsilon}$ by 0 and $\phi_{\varepsilon}$ by 1 in $\hat{\Omega}\setminus{\overline}{\Omega}$. Obviously we have ${\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon};\hat{\Omega})={\mathcal{F}}_{\varepsilon}(\s_{\varepsilon},\phi_{\varepsilon};{\Omega})$, therefore for any $O\in\A_s(\hat{\Omega})$ applying the localization described in Section \[sec:LocRes\] we can associate to each $\s_{\varepsilon}$ a function $ u_{\varepsilon}\in H^1(O)$ with mean value $0$ such that $\s_{\varepsilon}={\nabla}^\perp u_{\varepsilon}$ in $O$. By Theorem \[teo:localResult\] there exists $u\in\P(O)$ such that up to extraction $u_{\varepsilon}\rwstar u$. Eventually, by uniqueness of the limit, we get $$\s{ \!\!\mbox{{ \Large$\llcorner$}} }O=-[u]\nu_{J_u}^\perp\H^1{ \!\!\mbox{{ \Large$\llcorner$}} }(J_u\cap O).$$ Since we can cover ${\overline}{\Omega}\setminus S$ by finitely many sets $O\in \A_s(\hat{\Omega})$, this shows that $\s$ decomposes as $$\s\, =\, U(M_\s,\t_\s,\xi_\s)+\underbrace{\sum_{j=0}^N c_j \delta_{x_j}}_{\mu}.$$ By Lemma \[lemma:completediv\] there exists a rectifiable measure $\g=U(M_\g,\t_\g,\xi_\g)$ such that ${\nabla}\cdot(\s+\g)=0$ and $\H^1(M_\g\cap M_\s )=0$. Then there exists $u\in BV({\Omega})$ such that $Du =\s^{\perp}+\g^{\perp}$. From , we deduce $|Du|(S)=0$ which implies $|\mu|(S)=\sum|c_j|=0$. Hence $c_j=0$ for $j=0,\dots,N$ and $\s$ writes in the form $U(M_\s,\t_\s,\xi_\s)$.

Let us now use the local results of Section \[sec:LocRes\] to prove the lower bound.

Let $(\s_{{\varepsilon}},\phi_{{\varepsilon}})$ as in the statement of the theorem. Without loss of generality, we can suppose that ${\mathcal{F}}_{{\varepsilon}}(\s_{{\varepsilon}},\phi_{{\varepsilon}})<+{\infty}$. Theorem \[teo:sigma\_equicoercive\] then ensures the existence of a rectifiable measure $\s=U(M_\s,\t_\s,\xi_\s)$ with $\operatorname{supp}(\s)\subset\overline{\Omega}$ such that $\s_{{\varepsilon}}\rwstar\s$.\
Let $\hat{\Omega}$ be as in the previous proof and let us define $\mu=\Gamma-\liminf_{{\varepsilon}}{\mathcal{F}}_{{\varepsilon}}(\s_{{\varepsilon}},\phi_{{\varepsilon}})$ and $\lambda= \a|\s|+\H^{1}{ \!\!\mbox{{ \Large$\llcorner$}} }M_\s$. Consider the countable family of sets $\{O_{i}\}\subset\A_{S}(\hat{\Omega})$ made of the open rectangles $O_i\subset \hat{\Omega}\setminus S$ with vertices in ${\mathbf{Q}}^{2}$ and let $\psi_{i}:=1_{O_{i}}$. The local result stated in Theorem \[teo:localResult\] gives for any $i\in{\mathbf{N}}$ $$\mu(A)\geq\mu(O_{i}\cap A)\geq\lambda(O_{i}\cap A)=\int_{A}\psi_{i}\dif \lambda.$$ Therefore Lemma \[lem:braides\] gives $$\begin{aligned}
\Gamma-\liminf_{{\varepsilon}{\downarrow}0}{\mathcal{F}}_{{\varepsilon}}(\s_{{\varepsilon}},\phi_{{\varepsilon}})&=\mu(\hat{\Omega})\geq\lambda(\hat{\Omega})=\a|\s|(\overline{\Omega})+\H^{1}(M_\s) \end{aligned}$$ since $\sup_{i}\psi_{i}$ is the constant function $1$.

Upper bound {#sec:UppBound}
===========

A density result
----------------

In order to obtain the upper bound we first provide a density lemma. We show that measures which have support contained in a finite union of segments, are dense in energy.\
Without loss of generality let us assume that $\s\in \M_S({\overline}{\Omega})$ is such that ${\mathcal{E}}_\alpha(\s)<\infty$. In particular $\s=U(M_\s,\t_\s,\xi_\s)$ is a $\H^1$-rectifiable measure. Applying Lemma \[lemma:completediv\] we obtain a $\H^1$-rectifiable measure $\g=U(M_\g,\t_\g,\xi_\g)$ and a partition of ${\Omega}$ made of polyhedrons $\{{\Omega}_i\}$ such that $M_\g\subset \cup_i\de{\Omega}_i$, $\H^1(M_\s \cap\cup_i\de{\Omega}_i)=0$ and $\s+\g$ is divergence free.\
From the above properties, we can write $$\s^\perp+\g^\perp=Du$$ for some $u\in\P({\Omega})$. Our strategy is the following, using existing results [@Bel_Cham_Gold], we build an approximating sequence for $u$ on each ${\Omega}_j$ whose gradient is supported on a finite union of segments. We then glue these approximations together to obtain a sequence $(w_j)$ approximating $u$ in $\hat{\Omega}$. The main difficulty is to establish that $Dw_j { \!\!\mbox{{ \Large$\llcorner$}} }[\cup_i\de{\Omega}_i]$ is close to $Du{ \!\!\mbox{{ \Large$\llcorner$}} }[\cup_i\de{\Omega}_i]=\g^\perp$.

\[lem:uapproximation\] There exists a sequence $(w_j)\subset\P(\hat{\Omega})$ with the following properties:

1.  $w_j\to u$ weakly in $BV(\hat{\Omega})$,

2.  $\operatorname{supp}w_j\subset {\overline}{\Omega}$,

3.  $\limsup_{j{\rightarrow}{\infty}}{\mathcal{E}}_\a(w_j,1)\leq{\mathcal{E}}_\a(u,1)$,

4.  $J_{w_j}$ is contained in a finite union of segments for any $j\in {\mathbf{N}}$,

5.  $|Dw_j-Du|(\cup \de{\Omega}_i)\to0$.

*Step 1.* In order to apply the results of [@Bel_Cham_Gold], we first need to modify $u$ and the energy. Let us note the energy density function $f(t)=1+\a t$ and for $k\geq0$ and $t\geq 0$ let us introduce the approximation $$\label{eq:definizionefk}
 f_k(t):=\min\{(2^{k/2}+\a 2^{-k/2})\sqrt{t},f(t)\}.$$

![Graph of f and two of its approximations $f_{k_1} $ and $f_{k_2}$ with $k_1<k_2$.[]{data-label="fig:approssimazionef"}](approssimazionef){width=".3\textwidth"}

We have $0\leq f_k\leq f$ and $f_k\equiv f$ on $[2^{-k},+\infty)$. Notice that $f_k$ is continuous, sub-additive and increasing on $[0,+{\infty})$ and that $f_k(0)=0$ with $\lim_{t{\rightarrow}0} \frac{f_k(t)}{t}=+{\infty}$. We define the associated energy for functions $v\in\P(\hat{\Omega})$ as ${\mathcal{E}}_{f_k}(v,\hat {\Omega}):=\int_{J_{v}\cap\hat{\Omega}}f_k([v])\dH^1$.\
Now we note $\P_k(\hat{\Omega})$ the set of functions $v\in \P(\hat{\Omega})$ such that $v(\hat{\Omega})\subset 2^{-k}{\mathbf{Z}}$. For these functions we have $|v^+(x)-v^-(x)|\geq 2^{-k}$ for $\H^1$-almost every $x\in J_v$. Consequently, there holds $${\mathcal{E}}_{f_k}(v)\, =\, {\mathcal{E}}_f(v).$$ For each fixed $k\geq0 $, let us introduce the function $$u_k=2^{-k}\lfloor2^k u\rfloor$$ where $\lfloor t\rfloor$ denotes the integer part of the real $t$. Note that $u_k\in \P_k(\hat{\Omega})$ with $J_{u_k}\subset J_u$ and $\|u-u_k\|_\infty\leq 2^{-k}$. Notice also since ${\left}|(u_k^+-u_k^-)-(u^+-u^-){\right}|\leq 2^{-k}$ we have $$\label{eq:diffuuk}
|Du_k-Du|(\hat{\Omega})\, \leq\, 2^{-k}\H^1(J_u)$$ In particular $u_k\to u$ strongly in $BV(\hat{\Omega})$. Moreover, we see that $$\label{eq:Ef=Efk}
{\mathcal{E}}_{f_k}(u_k) \, =\, {\mathcal{E}}_f(u_k) \, \leq\ \, {\mathcal{E}}_f(u) + \a 2^{-k}\H^1(J_u).$$

*Step 2.* Let us approximate the function $u_k$. Let us fix $k\geq 0$ and ${\Omega}_i$. We can apply Lemma 4.1 of [@Bel_Cham_Gold] to the function $u_k{ \!\!\mbox{{ \Large$\llcorner$}} }{{\Omega}_i}$ and to the energy ${\mathcal{E}}_{f_k}(\cdot,{\Omega}_i)$. We obtain a sequence $(w^i_j)$ which enjoys the following properties: $$\label{eq:properties}
  \begin{aligned}
        & w^i_j({\Omega}_i)\subset u_k({\Omega}_i)\subset 2^{-k}{\mathbf{Z}},\quad \forall j\in{\mathbf{N}}, \mbox{ hence }w^i_j\in \P_k(\hat{\Omega}), \\
        & w^j_i\to u_k \mbox{ in } L^1({\Omega}_i)\mbox{ as }j{\rightarrow}+{\infty}, \\
        & \lim_{j{\rightarrow}+{\infty}} {\mathcal{E}}_{f_k}(w^i_j,{\Omega}_i)= \lim_{j{\rightarrow}+{\infty}} {\mathcal{E}}_{f}(w^i_j,{\Omega}_i)={\mathcal{E}}_{f}(u_k,{\Omega}_i),\\  
        & J_{w^j_i} \mbox{ is contained in a finite union of segments for any } j\in {\mathbf{N}},\\
        & \int_{\de {\Omega}_i}|Tw^i_j-Tu_k|\dH^1{\rightarrow}0\mbox{ where } T:BV({\Omega}_i){\rightarrow}L^1(\de {\Omega}_i) \mbox{ denotes the trace operator}.\\
     \end{aligned}$$ Let us now define globally $$w_j:=\sum_j w_j^i 1_{{\Omega}_i}.$$ From the above properties, we have $w_j\rwstar u_k$, $$\label{eq:convEfw}
\lim{\mathcal{E}}_{f}(w^i_j,\hat {\Omega})={\mathcal{E}}_{f}(u_k,{\Omega}_i)$$ and $$\label{eq:diffukwj}
|D w_j - D u_k|(\cup_i\de{\Omega}_i)\, \to\, 0\quad\mbox{as $j\to\infty$}.$$ Eventually, using a diagonal argument, we have proved the existence of a sequence $(w_j)\subset \P(\hat{\Omega})$ complying to items *(a)*,  *(b)* and *(d)* of the lemma. Moreover, item *(c)* is the consequence of  and  and item *(e)* follows from  and .

Going back to the $\H^1$-rectifiable measures $\s=U(M_\s,\t_\s,\xi_\s)$, we define the sequence $$\s_j:=-Dw^\perp_i-\g.$$ We recall that $\g=U(M_\g,\t_\g,\xi_\g)$ with $M_\g\subset \cup \de {\Omega}_i$. In particular $\g=-Du^\perp{ \!\!\mbox{{ \Large$\llcorner$}} }(\cup_i \de{\Omega}_i)$. We deduce from the previous lemma:

\[lem:density\] There exists a sequence $(\s_j)\in\M_S({\overline}{\Omega})$ with the properties:

1.  $\s_j{\rightarrow}\s$ with respect to weak-$*$ convergence of measures,

2.  $\s_j= U(M_{\s_j},\t_{\s_j},\xi_{\s_j})$ with $M_{\s_j}$ contained in a finite union of segments,

3.  $\limsup_{j{\rightarrow}{\infty}}{\mathcal{E}}_\a(\s_j,1)\leq{\mathcal{E}}_\a(\s,1)$.

Construction of a recovery sequence
-----------------------------------

Let us prove the $\Gamma$-limsup inequality stated in Theorem \[teo:sigma\_limsup\]. Recall that the latter consists in finding a sequence $(\s\e,\phi\e)$ for any given couple $(\s,\phi)\in\M(\overline{\Omega},{\mathbf{R}}^2)\times L^1({\Omega})$ such that $\s\e\rwstar\s$, $\phi\e{\rightarrow}\phi$ in $L^1({\Omega})$ and $$\label{eq:limsup}
\limsup_{{\varepsilon}{\downarrow}0}{\mathcal{F}}\e(\s\e,\phi\e)\leq{\mathcal{E}}_\a(\s,\phi).$$ When ${\mathcal{E}}_\a(\s,\phi)=+\infty$ the inequality is valid for any sequence therefore by definition  we can assume $\s=U(M,\t,\xi)$ and $\phi=1$. Furthermore by density Lemma \[lem:density\] is sufficient to consider measures of the form $$\label{eq:strutturadense}
\s=\sum_{i=1}^n U(M_i,\t_i,\xi_i),$$ where $M_i$ is a segment, $\t_i\in{\mathbf{R}}_+$ is $\H^1$-a.e. constant and $\xi_i$ is an orientation of $M_i$ for each $i$. Without loss of generality we can suppose that for each couple of segments $M_i$, $M_j$, for $i \neq j$, the intersection $M_i\cap M_j$ is at most a point (called branching point) not belonging to the relative interior of $M_i$ and $M_j$. We firstly produce the estimate  for $\s$ composed by a single segment thus let us assume $\s=\t e_1\cdot \H^1{ \!\!\mbox{{ \Large$\llcorner$}} }(0,l)\times\{0\}$.

*Notation:* Let us fix the values $$\label{eq:infinitesimi}
a_{\varepsilon}:=\begin{dcases}
         \frac{\t\a\,{\varepsilon}}{2} &\mbox{ if } \a>0\\
         \quad{\varepsilon}& \mbox{ if } \a=0
        \end{dcases},
\qquad\qquad
b_{\varepsilon}:={\varepsilon}\ln{\left}(\frac{1-\eta}{{\varepsilon}}{\right})
\qquad\mbox{ and }\qquad
r\e=\max\{{\varepsilon},a\e\}.$$ Let $d_\infty(x,S)$ be the distance function from $x$ to the set $S\subset {\Omega}$ relative to the infinity norm on ${\mathbf{R}}^2$ and $Q_r(P)=\{x\in{\mathbf{R}}^2: d_\infty(x,P)\leq r\}$ the square centered in $P$ of size $2r$ and sides parallel to the axes. Introduce the sets $$\begin{aligned}
\label{def:rectangular}
I_{a\e}&:=\{x\in{\mathbf{R}}^2: d_\infty(x,[0,l]\times\{0\})\leq a\e\}\cup Q_{r\e}(0,0)\cup Q_{r\e}(l,0),\\
I_{b\e}&:=\{x\in{\mathbf{R}}^2: d_\infty(x,I_{a\e})\leq b\e\},\\
I_{c\e}&:=\{x\in{\mathbf{R}}^2: d_\infty(x,(I_{a\e}\cup I_{b\e}))\leq {\varepsilon}\},\\
I_{d\e}&:={\Omega}\setminus(I_{a\e}\cup I_{b\e}\cup I_{c\e}),\end{aligned}$$ and define $R\e= I_{a\e}\setminus(Q_{r\e}(0,0)\cup Q_{r\e}(l,0))$.

![Example of the neighborhoods of the segment $[0,l]\times\{0\}$. On the left the case $r\e={\varepsilon}$ on the right the case in which $r\e=a\e>{\varepsilon}$. The stripped region is $R\e$ and $I_{a\e}=R\e\cup(Q_{r\e}(0,0)\cup Q_{r\e}(l,0))$. Remark that $\operatorname{supp}(\rho\e)=B(0,{\varepsilon})$.[]{data-label="fig:intorninuovo"}](nuovointorni){width="\textwidth"}

*Costruction of $\s\e$:* We build $\s\e$ as a vector field supported on $I_{a\e}$. In particular we add together three different constructions performed respectively on $R\e$, $Q_{r\e}(0,0)$ and $ Q_{r\e}(l,0)$. Let $r=r_{\varepsilon}/{\varepsilon}$ and consider the problem

$$\begin{dcases}
\Delta u=\pm\t\delta_{x_0}*\rho &\mbox{ on }  Q_{r}(0,0),\\
\frac{\de u}{\de \nu}=\frac{\pm\t}{\H^1(\Sigma)}  &\mbox{ on }  \Sigma^\pm=\{x\in{\mathbf{R}}^2: x_1=\pm1,\;|x_2|\leq\frac{\t\a}{2}\}.
\end{dcases}$$

![image](problema00){width=".6\textwidth"}

Let $u^+$ be the solution relative to the problem in which every occurrence of $\pm$ is replaced by $+$ and let $u^-$ defined accordingly. Then set $$\label{def:se}
 \s\e=
 \begin{dcases}
       \frac{{\nabla}u^+(x/{\varepsilon})}{{\varepsilon}} &\mbox{ on } Q_{r\e}(0,0),\\
       \quad\frac{\t}{2a_{\varepsilon}}\cdot e_1 &\mbox{ on } R\e,\\
       \frac{{\nabla}u^-((x-(l,0))/{\varepsilon})}{{\varepsilon}} &\mbox{ on } Q_{r\e}(l,0).
  \end{dcases}$$ By construction we have that ${\nabla}\cdot\s\e=\t(\d_{(0,0)}-\d_{(l,0)})\ast\rho\e$ and $\s\e\rwstar\s$. Let us point out as well that there exists a constant $c(\a,\t)$ such that $$\label{eq:energiapunti}
c(\a,\t):=\int_{Q_{r\e}(l,0)}|\s\e|^2\dx=\int_{Q_{r\e}(0,0)}|\s\e|^2\dx=\int_{Q_{r}(0,0)}{\left}|{\nabla}u^+(x){\right}|^2\dx=\int_{Q_{r}(0,0)}{\left}|{\nabla}u^-(x){\right}|^2\dx.$$

*Costruction of $\phi\e$:* Most of the properties of $\phi\e$ are a consequence of the inequalities obtained in Theorem \[teo:localResult\] and the structure of $\s\e$. On one hand we need $\phi\e$ to attain the lowest value possible on $I_{a\e}$ in order to compensate the concentration of $\s\e$ in this set, on the other, as shown in inequality , we need to provide the optimal profile for the transition from this low value to $1$. For this reasons we are led to consider the following ordinary differential equation associated with the optimal transition $$\label{eq:ode}
\begin{dcases}
w'_{\varepsilon}= \frac{1}{{\varepsilon}}(1-w_{\varepsilon}), \\
w_{\varepsilon}(0)= \eta. \\
\end{dcases}$$ Observe that $w_{\varepsilon}=1-(1-\eta)\exp{\left}(\frac{-t}{{\varepsilon}}{\right})$ is the explicit solution of equation  and set $$\label{def:phie}
  \phi_{\varepsilon}(x):=\begin{dcases}
        \eta&\mbox{ if } x\in I_{a_{\varepsilon}},\\
           w_{\varepsilon}(d_\infty(x, I_{a\e})) &\mbox{ if } x\in I_{b_{\varepsilon}},\\
         d_\infty(x, I_{b\e})-{\varepsilon}+1&\mbox{ if } x\in I_{c_{\varepsilon}},\\
            1&\mbox{ otherwise.}
                \end{dcases}$$ The choice of the behavior in the region $I_{c\e}$ is given by the fact that following the optimal profile we will reach the value $1$ only at $+\infty$ thus a linear correction on a small set ensures that this value is achieved with a small cost in energy.

*Evaluation of ${\mathcal{F}}\e(\s\e,\phi\e)$:* We prove inequality  for the sequence we have produced. Since the sets $I_{a\e}$, $I_{b\e}$, $ I_{c\e}$ and $I_{d\e}$ are disjoint we can split the energy as follows $$\label{dislimsup0}
{\mathcal{F}}\e(\s_{\varepsilon},\phi\e)={\mathcal{F}}\e(\s_{\varepsilon},\phi\e;I_{a\e})+{\mathcal{F}}\e(\s_{\varepsilon},\phi\e;I_{b\e})+{\mathcal{F}}\e(\s_{\varepsilon},\phi\e;I_{c\e})+{\mathcal{F}}\e(\s_{\varepsilon},\phi\e;I_{d\e})$$ and evaluate each component individually. Since $\s\e$ is null and $\phi\e$ is constant and equal to $1$ in $I_{d\e}$ we have that ${\mathcal{F}}\e(\s,\phi\e;I_{d\e})=0$. For the other components we strongly use the definitions in  and . Firstly we split again the energy on the set $I_{a\e}$ as following $${\mathcal{F}}\e(\s\e,\phi\e;I_{a\e})={\mathcal{F}}\e(\s\e,\phi\e;R\e)+{\mathcal{F}}\e(\s\e,\phi\e;Q_{r\e}(0,0))+{\mathcal{F}}\e(\s\e,\phi\e;Q_{r\e}(l,0)).$$ Now identity  leads to the estimate $${\mathcal{F}}\e(\s\e,\phi\e;Q_{r\e}(0,0))={\mathcal{F}}\e(\s\e,\phi\e;Q_{r\e}(l,0))=\frac{\eta^2}{2{\varepsilon}}c(\a,\t)+\frac{(1-\eta)^2}{2{\varepsilon}}\;r\e^2$$ and $${\mathcal{F}}\e(\s\e,\phi\e;R\e)={\left}[\frac{1}{2{\varepsilon}}\eta^2 {\left}|\frac{\t}{2a\e}{\right}|^2+\frac{(1-\eta)^2}{2{\varepsilon}}{\right}]|R\e|\leq {\left}[\frac{(\t\eta)^2}{8{\varepsilon}a\e^2}+\frac{1}{2{\varepsilon}}{\right}]2a\e l.$$ Then passing to the limsup we obtain $$\label{dislimsup1}
 \limsup_{{\varepsilon}{\downarrow}0}{\mathcal{F}}\e(\s\e,\phi\e;I_{a\e})\leq  \t\a l= \t_{}\a\H^1([0,l]\times\{0\}).$$ To obtain the inequality on the sets $I_{b\e}$ and $I_{c\e}$ we are going to apply the Coarea formula therefore let us observe that for both $d_{\infty}(x, I_{a\e})$ and $d_{\infty}(x, I_{b\e})$ there holds $|{\nabla}d_{\infty}(x,\cdot)|=1$ for a.e. $x\in{\Omega}$ and that there exist a constant $k=k(\a,\t)$ such that the level lines $\{d_{\infty}(x,\cdot)=t\}$ have $\H^1$ length controlled by $2l+kt$. In force of these remarks we obtain $$\begin{aligned}
\label{dislimsup2}
 {\mathcal{F}}\e(\s\e,\phi\e;I_{b\e})&=\int_{I_{b\e}} {\left}[\frac{{\varepsilon}}{2}|{\nabla}\phi_{\varepsilon}|^2  + \frac{(1-\phi_{\varepsilon})^2}{2{\varepsilon}} {\right}]|{\nabla}d_{\infty}(x,I_{a\e})|\dx\nonumber\\
 &=\int_{0}^{b_{\varepsilon}}{\left}[\frac{(1-w_{\varepsilon}(t))^2}{2{\varepsilon}}+\frac{{\varepsilon}}{2}|w'_{\varepsilon}(t)|^2 {\right}] \H^1(\{d_\infty(\cdot, I_{a\e})=t\})\;\dt\nonumber\\
 &\leq(2l+k{\varepsilon}){\left}[\frac{1}{2}(1-w_{\varepsilon}(t))^2{\right}]_{0}^{b_{\varepsilon}}\nonumber\\
 &={\left}(l-\frac{k{\varepsilon}}{2}{\right}){\left}[(1-\eta)^2-{\varepsilon}^2{\right}]\xrightarrow [{\varepsilon}{\downarrow}0]{} l= \H^1([0,l]\times\{0\})\end{aligned}$$ and $$\begin{aligned}
\label{dislimsup3}
 {\mathcal{F}}\e(\s\e,\phi\e;I_{c\e})&=\int_{I_{c\e}} {\left}[\frac{{\varepsilon}}{2}|{\nabla}\phi_{\varepsilon}|^2  + \frac{(1-\phi_{\varepsilon})^2}{2{\varepsilon}} {\right}]|{\nabla}d_{\infty}(x,I_{b\e})|\dx\nonumber\\
 &=\int_{0}^{{\varepsilon}}{\left}[\frac{(1-t+{\varepsilon}-1)^2}{2{\varepsilon}}+\frac{{\varepsilon}}{2} {\right}] \H^1(\{d_\infty(\cdot, I_{b\e}\cup I_{a\e})=t\})\;\dt\nonumber\\
 &\leq(2l+k{\varepsilon}){\left}[\frac{(t-{\varepsilon})^3}{6{\varepsilon}}+\frac{{\varepsilon}}{2} t{\right}]_{0}^{{\varepsilon}}\nonumber\\
 &=(2l+k{\varepsilon})\,{\varepsilon}^2\,\frac{2}{3}\xrightarrow [{\varepsilon}{\downarrow}0]{} 0.\end{aligned}$$ Finally adding up equations , ,  and  we obtain $$\limsup_{{\varepsilon}{\downarrow}0}{\mathcal{F}}\e(\s_{\varepsilon},\phi\e)\leq (1+\a\;\t)\; \H^1([0,l]\times\{0\}).$$

*Case $\s$ of the form :*

Let us call $\s\e^i$, $\phi\e^i$ the functions obtained above for each $\s_i=\t_i\xi_i \H^1{ \!\!\mbox{{ \Large$\llcorner$}} }M_i$ and set $$\s\e=\sum_{i=1}^n \s\e^i,\qquad\qquad \phi\e=\min_{i}\;\phi\e^i.$$ Let us remark that in force of the local construction we have made at the ending points of each segment and since $\s$ satisfies equation  for each ${\varepsilon}$ there holds $${\nabla}\cdot\s_{\varepsilon}\, =\, {\left}(N\delta_{x_0} - \sum_{j=1}^N \delta_{x_j} {\right})\ast \rho_{\varepsilon}.$$ We now prove inequality . The following inequality holds true $$\begin{aligned}
{\mathcal{F}}\e(\s\e,\phi\e)&=\int_{\Omega}\frac{1}{2{\varepsilon}}|\min_i \phi\e^i|^2|\sum_{i=1}^n \s^i\e|^2 +\frac{{\varepsilon}}{2}|{\nabla}(\min_i \phi\e^i)|^2 + \frac{(1-\min_i \phi\e^i)^2}{2{\varepsilon}}\dx\nonumber\\
&\leq\int_{\Omega}\frac{1}{2{\varepsilon}}|\min_i \phi\e^i|^2|\sum_{i=1}^n \s^i\e|^2 \dx+\sum_{i=1}^n\int_{\Omega}\frac{{\varepsilon}}{2}|{\nabla}\phi\e^i|^2 + \frac{(1-\phi\e^i)^2}{2{\varepsilon}}\dx,\label{dis:decomposizione}\end{aligned}$$ therefore we look into an estimation of the first integral in the latter. Observe that for ${\varepsilon}$ sufficiently small we can assume that all the $R\e^i$ are pairwise disjoint thus we study the behavior in the squares. Let $M_{i_1},\dots,M_{i_{m_P}}$ be the segments meeting at a branching point $P$. For $j=i_1,\dots,i_{m_P}$ let us call $Q_{r\e^j}(P)$ the squared neighborhood of $P$ relative to the segment $M_j$ as constructed previously. Let us recall that by definition $\phi\e$ is constant and equal to $\eta$ on $\cup_{j=i_1}^{m_P} Q_{r\e^j}(P)$ then we have the estimation $$\begin{aligned}
\int\limits_{\cup_{j=i_1}^{m_P} (R\e^j\cup Q_{r\e^j}(P))}\frac{\phi^2\e}{2{\varepsilon}}\,|\s\e|^2\dx&=\sum_{j=i_1}^{m_P}\int_{R\e^j}\frac{\phi^2\e}{2{\varepsilon}}\,|\s\e|^2\dx+\int_{\cup_{i=i_1}^{m_P} Q_{r\e^j}(P)}\frac{\phi^2\e}{2{\varepsilon}}\,|\sum_{j=i_1}^{m_P}\s^j\e|^2\dx\nonumber\\
&\leq \sum_{j=i_1}^{m_P}\int_{R\e^j}\frac{\phi^2\e}{2{\varepsilon}}\,|\s\e|^2\dx+m_P\,\frac{\eta^2}{2{\varepsilon}}\sum_{j=i_1}^{m_P}\int_{Q_{r\e^j}(P)}|\s^j\e|^2\dx\nonumber\\
&\leq \sum_{j=i_1}^{m_P}\int\limits_{(R\e^j\cup Q_{r\e^j}(P))}\frac{1}{2{\varepsilon}}|\phi\e^j|^2|\s^j\e|^2\dx+\underbrace{(m_P-1)\,{\left}(\sum_{j=i_1}^{i_{m_P}}{c(\a,\t_j)}{\right})\frac{\eta^2}{2{\varepsilon}}}_{c(m_P,\a,\t_{i_1},\dots,\t_{i_{m_P}}){\varepsilon}}.\label{dis:punti}\end{aligned}$$ Applying inequality  on each branching point in equation  and recomposing the integral gives $$\begin{aligned}
\limsup_{{\varepsilon}{\downarrow}0}{\mathcal{F}}\e(\s\e,\phi\e)&\leq\limsup_{{\varepsilon}{\downarrow}0}\sum_{i=1}^n{\mathcal{F}}\e(\s^i\e,\phi^i\e) +n\; c(n,\a,\t_i,\dots,\t_n){\varepsilon}\\
&\leq \sum_{i=1}^n(1+\a\;\t_i)\; \H^1(M_i)\\
&=\int_{\operatorname{supp}(\s)}(1+\a\,\t)\dH^1 ={\mathcal{E}}_\a(\s,1)\end{aligned}$$ which ends the proof.

Numerical Approximation {#sec:NumApprox}
=======================

Equations
---------

In this section we present some numerical simulations of the $\Gamma$-convergence result we have shown. The first issue we address is how to impose the divergence constraint. To this aim is convenient to introduce the following notation $$\begin{aligned}
 f\e&= {\left}(N\delta_{x_0} - \sum_{j=1}^N \delta_{x_j} {\right})\ast \rho_{\varepsilon},\\
 G\e(\s,\phi)&=\begin{dcases}
                   \int_{{\Omega}} {\left}[\frac{1}{2{\varepsilon}}|\phi|^2 |\s|^2 {\right}]\dx& \text{if } \s\in V\e,\\
                   +{\infty}&\mbox{otherwise in } L^2({\Omega},{\mathbf{R}}^2),
                 \end{dcases}\\
 \Lambda\e(\phi)&=\begin{dcases}
        \int_{{\Omega}}{\left}[ \frac{{\varepsilon}}{2}|{\nabla}\phi|^2 +\frac{(1-\phi^2)}{2{\varepsilon}} {\right}]\dx& \text{if } \phi\in W\e,\\
          +{\infty}&\mbox{otherwise in } L^1({\Omega}).
        \end{dcases}\end{aligned}$$ Then let us observe that the following equality holds $$\min_{\s\in L^2({\Omega},{\mathbf{R}}^2)} G\e(\s,\phi)= \inf_{\s\in L^2({\Omega},{\mathbf{R}}^2)}{\left}\{ \sup_{u\in H^1({\Omega})}\int_{{\Omega}} \frac{1}{2{\varepsilon}}|\phi|^2 |\s|^2 +u({\nabla}\cdot\s-f\e)\dx{\right}\}.$$ By Von Neumann’s min-max Theorem [@Att_Hed_Butt Thm. 9.7.1] we can exchange inf and sup obtaining for each ${\varepsilon}>0$ and $\phi \in W\e$ $$\begin{aligned}
\min_{\s} G\e(\s,\phi)&=\sup_{u}\inf_{\s}\int_{{\Omega}} \frac{1}{2{\varepsilon}}|\phi|^2 |\s|^2 -({\nabla}u\s+u f\e)\dx\\
&=-\min_{u}\int_{{\Omega}} \frac{{\varepsilon}|{\nabla}u|^2}{2|\phi|^2}  +uf\e\dx=-\min_u G'\e(u,\phi).\end{aligned}$$ With the relation $\s=\frac{{\varepsilon}{\nabla}u}{\phi^2}$. This naturally leads to the following alternate minimization problem: given an initial guess $\phi_0$ we define $$\begin{aligned}
\s_j&:=\frac{{\varepsilon}{\nabla}u_j}{\phi_j^2} \quad\mbox{ where } \quad u_j:=\operatorname{argmin}G'\e(u,\phi_j), \\
\phi_{j+1}&:=\operatorname{argmin}G\e(\s_j,\phi)+\Lambda\e(\phi).\end{aligned}$$ We supplement the alternate minimization with a third step where we optimize the component $\Lambda\e$ with respect to a deformation of the domain. Let us describe this step. For $T:{\Omega}{\rightarrow}{\Omega}$ a smooth map we define $$\phi_T=\phi\circ T(x) \quad\mbox{and}\quad\Lambda\e(T)=\Lambda\e(\phi_T).$$ By a change of variables we get $$\Lambda\e(T)=\int_{\Omega}{\left}[\frac{{\varepsilon}}{2}|({\nabla}T\circ T^{-1}){\nabla}\phi|^2 +\frac{(1-\phi^2)}{2{\varepsilon}}{\right}]\det({\nabla}T^{-1})\dy$$ In particular we choose $T$ to be of the form $x+V(x)$ and evaluate the gradient obtaining $$\langle d\Lambda\e(T),W\rangle=\int_{{\Omega}}{\left}[{\varepsilon}({\nabla}\phi_T;{\nabla}W{\nabla}\phi_T)-\frac{{\varepsilon}}{2}|{\nabla}\phi_T|^2{\nabla}\cdot W-\frac{1}{2{\varepsilon}}(1-\phi_T)^2{\nabla}\cdot W{\right}]\dx$$ Representing in $H^1({\Omega},{\Omega})$ the gradient of the functional $\Lambda\e$ evaluated for $T(x)=x$ obtains the elliptic problem $$\int_{\Omega}{\left}({\nabla}V,{\nabla}W{\right})\dx+\int_{{\Omega}}{\left}[{\varepsilon}({\nabla}\phi;{\nabla}W{\nabla}\phi)-\frac{{\varepsilon}}{2}|{\nabla}\phi|^2{\nabla}\cdot W-\frac{1}{2{\varepsilon}}(1-\phi)^2{\nabla}\cdot W{\right}]\dx=0$$ This method enhances the length minimization process since, as we already pointed out, $\Lambda\e$ is a variation of Modica-Mortola’s functional.

Discretization
--------------

We define a circular domain ${\Omega}$ containing the points in $S$ endowed with a uniform mesh and four values $\a$, ${\varepsilon}_{in}, {\varepsilon}_{end}$ and $N_{iter}$ and a gaussian convolution kernel $\rho_{{\varepsilon}_{end}}$ in order to define $f\e$. For the discrete spaces we have chosen for $u$, $\phi$ and the vector field $V$ to be piecewise polynomials of order $1$. This leads to the following algorithm

$S=\{x_0,\ldots,x_N\}$, ${\varepsilon}_{in},\quad{\varepsilon}_{end}$,$N_{iter}$,$\a$, index. Set $f\e=(N\d_{x_0}-\sum_{i=1}^N \d_{x_i})*\rho_{{\varepsilon}_{end}}$ and $\phi_0=1$ ${\varepsilon}_j= {\left}(\frac{j-N_{iter}}{N_{iter}}{\right}){\varepsilon}_{in}-{\left}(\frac{j}{N_{iter}}{\right}){\varepsilon}_{end}$ $\tilde\phi\leftarrow L^1$-projection of $\phi_{j-1}^2$ Set $u_j$ as the minimizer of $G'_{{\varepsilon}_j}(\cdot,\phi_{j-1})$ Set $\s_j=\frac{{\varepsilon}_j {\nabla}u_j}{\tilde\phi_{j-1}}$ Set $\phi_j$ as the minimizer of $G_{{\varepsilon}_j}(\s_j,\cdot)+\Lambda\e(\cdot)$ Solve $\langle d\Lambda_{{\varepsilon}_j}(T),W\rangle=0$ Set $\phi_j=\phi_j(x+T)$ Set $\phi_j=\max\{\eta,\phi_j\}$ $\phi_{N_{iter}}$, $\s_{N_{iter}}$.

We have implemented the algorithm in FREEFEM++. In the next figures we show the graphs obtained for the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ via the approximation algorithm with the choices $\a=0.05$, ${\varepsilon}_{in}=0.5$, ${\varepsilon}_{end}=0.05$, $\a=0.05$, $N_{iter} = 500$ and $index = 300$. We have chosen to make simulations for points located on the vertices of regular polygons of respectively 3, 4, 5 and 6 vertices. This choice allows a direct visual perception of the results.

![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](3puntiFA.png "fig:"){width="3cm"}![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](4puntiFA.png "fig:"){width="3cm"}![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](5puntiFA "fig:"){width="3cm"}![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](6puntiFA "fig:"){width="3cm"} ![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](3puntiA "fig:"){width="3cm"}![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](4puntiA "fig:"){width="3cm"}![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](5puntiA "fig:"){width="3cm"}![Graph of the couple $(\s_{N_{iter}},\phi_{N_{iter}})$ obtained via Algorithm \[alg:1\] in the case of 3, 4, 5 and 6 points located on the vertices of a regular polygon.[]{data-label="fig:simulazione1"}](6puntiA "fig:"){width="3cm"}

Finally let us point out the need of the third minimization step. In the following figure we have the graph of the solution obtained for a simulation in which the third step is omitted. Even from visual perception is possible to recognize that the solution differs both from the solution of the Steiner Tree and the minimizer of the ${\mathcal{E}}_\a$ energy as evident from the figure. Furthermore we do not obtain the classical straight segments we would expect in studying geodesic in the euclidean metric. We suppose that these alterations are a consequence of the alternate minimization method that could not lead to a global minimum and therefore we introduced the third step in the algorithm in order to perturbate local solutions.

![On the left: Graph of $\phi$ obtained via Algorithm \[alg:1\] in which the gradient descend method is omitted. On the right: in red, one of the solutions to the Steiner problem for four points on the vertices of a square, while in blue, a minimizer of the energy ${\mathcal{E}}_\a$ associated to the same constraint.[]{data-label="fig:4sbagliato"}](4sbagliato "fig:"){width="3.5cm"}![On the left: Graph of $\phi$ obtained via Algorithm \[alg:1\] in which the gradient descend method is omitted. On the right: in red, one of the solutions to the Steiner problem for four points on the vertices of a square, while in blue, a minimizer of the energy ${\mathcal{E}}_\a$ associated to the same constraint.[]{data-label="fig:4sbagliato"}](soluzionisteinerapprossimate "fig:"){width="3cm"}

To ensure that this step is reasonable we have studied several experiments and plotted the numerical energy of each experiment and observed that we are always led to a lower energy. The following plot shows the behavior of the energy for the iterations concerning the third step for the first two solutions in figure \[fig:simulazione1\]. Is possible to observe that although there are increments

![Behavior of the estimated energy of the last 200 iterations of Algorithm \[alg:1\] referring to the first two figures in figure \[fig:simulazione1\] .[]{data-label="fig:grafici_energia"}](grafici_energia){width=".6\textwidth"}

Acknowledgments {#acknowledgments .unnumbered}
===============

The authors have been supported by the ANR project Geometrya, Grant No. ANR-12-BS01-0014-01. A.C. also acknowledges the hospitality of Churchill College and DAMTP, U. Cambridge, with a support of the French Embassy in the UK, and a support of the Cantab Capital Institute for Mathematics of Information.
