<figure><img src="image_1.png" /></figure>

HAL
open science

A CONVEX RELAXATION METHOD FOR FREE BOUNDARY PROBLEMS

Guy Bouchitté, Minh Phan

<figure><img src="image_2.png" /></figure>

To cite this version:

Guy Bouchitté, Minh Phan. A CONVEX RELAXATION METHOD FOR FREE BOUNDARY PROBLEMS. 2024. hal-04820298

HAL Id: hal-04820298

https://univ-tln.hal.science/hal-04820298v1

Preprint submitted on 5 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A CONVEX RELAXATION METHOD FOR FREE BOUNDARY PROBLEMS

GUY BOUCHITTÉ AND MINH PHAN

**ABSTRACT.** We derive a convex relaxation principle for a large class of non convex variational problems where the functional to be minimized involves a one homogeneous gradient energy. This applies directly to free boundary or multiphase problems in the case of the classical total variation or of some anisotropic variants. The underlying argument is an exclusion principle which states that any global minimizer avoids taking values in the intervals where the lower order potential is nonconvex. This allows using duality methods and deriving a saddle point characterization of the global minimizers. A numerical validation of our principle is presented in the case of several free boundary and multiphase problems that we treat through a primal-dual algorithm. The accuracy of the interfaces and the convergence of the algorithm benefit in a large way of a new epigraphical projection method that we introduced to tackle the non differentiability of the convexified Lagrangian.

**Keywords:** duality, relaxation, non convex problems, free boundaries

**Mathematics Subject Classification:** 49J45, 28A25

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. Exclusion principle and main results</td><td>2</td></tr><tr><td>3. Proof of the exclusion principle</td><td>5</td></tr><tr><td>4. Counterexamples</td><td>8</td></tr><tr><td>5. Application to free boundary and multiphase problems</td><td>12</td></tr><tr><td>5.1. A three phases problem</td><td>12</td></tr><tr><td>5.2. A four phases problem</td><td>16</td></tr><tr><td>6. Min-max approach. A primal-dual algorithm for non differentiable Lagrangians</td><td>17</td></tr><tr><td>6.1. Explicit scheme with projection on epi(g**)</td><td>19</td></tr><tr><td>6.2. Semi-implicit scheme with projection on epi(g*)</td><td>20</td></tr><tr><td>Appendix A. Proof of the duality Theorem</td><td>22</td></tr><tr><td>A.1. Background around the space BV</td><td>22</td></tr><tr><td>A.2. Convexified problem in N + 1 dimensions</td><td>23</td></tr><tr><td>A.3. Proof of Theorem 3.1</td><td>26</td></tr><tr><td>References</td><td>28</td></tr></table>

## 1. INTRODUCTION

We are interested in a new global approach for non convex variational problems of the kind

$$ (P) \qquad \inf \left\{ \int_{\Omega} (h(\nabla u) + g(u) - f(x)u) dx : u \in W^{1,1}(\Omega), u = u_0 \text{ on } \Gamma \right\} $$

where $g: \mathbb{R} \to (-\infty, +\infty]$ is a non convex function, $\Omega$ a bounded domain of $\mathbb{R}^N$ and $\Gamma$ a Borel subset of $\partial\Omega$. In view of applications to free boundary problems; this function $g$ is assumed to be merely lower semicontinuous with possibly several discontinuity points. Recently a general duality theory
has been introduced for such problems [7] [9] in the case where the integrand $h: \mathbb{R}^N \to [0, +\infty)$ is a convex function satisfying growth conditions of order $p > 1$ (superlinear case). This theory allows a direct approach to global minimizers by means of min-max formulations, and it seems to be powerful in the case of free boundary problems; let us mention, for example, the numerical methods developed in computer vision [19, 20] which are based on the same kind of convexification principle. The main drawback of this method is that we have to increase the dimension $N$ of the initial problem to the dimension $N + 1$, which is costly in terms of numerical simulations.

Surprisingly, when the convex integrand growths linearly at infinity ($p = 1$), a major simplification of the method occurs when the two following conditions are satisfied:

a) The integrand $h$ is one homogeneous

$$h(tz) = t h(z) \quad \text{for all } z \in \mathbb{R}^N \text{ and } t \ge 0.$$

b) The boundary data $u_0$ does not range in the subset $Z := \{t \in \mathbb{R} : g^{**}(t) < g(t)\}$.

It turns out that, in this case, the duality scheme developed in [7] can be reduced to the dimension $N$ with the consequence that the initial problem ($\mathcal{P}$) can be relaxed directly (in a sense to be precised later) by substituting $g$ with its convexified $g^{**}$. This is the main contribution of the paper (Theorem 2.4). It rests upon an *exclusion principle* (see Theorem 2.2) that we deduce from the original duality theory in dimension $N + 1$. In fact this exclusion principle has been discovered by A. Visintin in the 1990's (see [18]) in the different context of Cahn-Hilliard type models where $g$ is a double well potential and $u$ represents the density of a two-phase fluid. In this work where the Dirichlet boundary condition is replaced by a prescribed total mass constraint $\int_{\Omega} u = m$, the arguments to prove such an exclusion principle rely on rearrangements tricks which apply also when the term $\int_{\Omega} h(\nabla u)dx$ is replaced by a general convex 1-homogeneous functional satisfying a co-area formula.

The method we use in this paper is based on the calibration theory developed in [7] that we adapt to the case $p = 1$ and assuming a non vanishing source term $f \in L^{\infty}(\Omega)$. In a preliminary step, the primal problem ($\mathcal{P}$) has to be relaxed in the space $BV(\Omega)$ (see (2.4)). In this formulation the energy to be minimized (see (2.5)) includes, besides the Cantor part, an additional term concentrated on the jump set of $u$ and a contribution on the part of boundary $\Gamma$ where the condition $u = u_0$ is not satisfied. Then, assuming that minimizers range in a finite interval $I := [m, M]$, we construct a dual problem ($\mathcal{P}^*$) on the cylinder $\Omega \times I$ in the same line as in [7], where the competing vector fields $\sigma \in L^{\infty}(\Omega \times I; \mathbb{R}^{N+1})$ are divergence free and satisfy a vanishing normal trace condition on the lateral boundary subset $\Gamma \times I$. Our main result is then obtained in a straightforward way by comparing this dual problem to the one obtained when substituting $g$ with $g^{**}$.

Next we investigate the necessity of conditions a) and b). In Section 4, we provide two counterexamples. In particular the necessity of condition b) is illustrated by considering a minimal surface problem with free boundary.

In the remaining of the paper, we apply our result to a class of free boundary and multiphase problems. To that aim we choose $g$ so that the coincidence set $\{t \in \mathbb{R} : g(t) = g^{**}(t) < +\infty\}$ consists of a finite number of values $\{t_1, ..., t_n\}$ (sorted in increasing order) while the boundary data $u_0$ avoids all the intermediate values. Doing so we can reduce some $n$-multiple phases problem to a classical convex optimization problem for which many efficient algorithms are available. Note however that, by construction, $g^{**}$ is piecewise affine and not differentiable at the values $t_i$ corresponding to the searched phases. To overcome this difficulty we propose two numerical methods: in the first we use an explicit primal-dual scheme with a geometric (polyhedral) projection on the epigraph of $g^{**}$ whereas the second one combines a semi-implicit algorithms with a projection on the epigraph of $g^*$.
The efficiency of the method is evidenced in many numerical simulations including the case where the homogeneous integrand $h$ is associated with a crystalline norm.

## 2. EXCLUSION PRINCIPLE AND MAIN RESULTS

Let $\Omega$ be an open bounded domain of $\mathbb{R}^N$ with Lipschitz boundary and let $\Gamma$ be a Borel subset of $\partial\Omega$. We denote by $\nu_\Omega$ the unit exterior normal on $\partial\Omega$. For a given $u_0 \in L^1(\partial\Omega)$, we consider the infimum problem

$$ \inf \left\{ \int_{\Omega} [h(\nabla u) + g(u) - f(x)u] dx : u \in W^{1,1}(\Omega), u = u_0 \text{ on } \Gamma \right\} \quad (2.1) $$

We recall that, by Gagliardo's Theorem, the trace map: $u \in W^{1,1}(\Omega) \to \text{Tr}(u) \in L^1(\partial\Omega)$ is continuous and *surjective*.

The standing assumptions that we make on $h, g$ and $f$ are listed below:

(H1) The function $h: \mathbb{R}^N \to (-\infty, +\infty]$ is convex lower semicontinuous, positively 1-homogeneous and satisfies the standard growth condition

$$ \forall z \in \mathbb{R}^N, \quad C_1|z| \le h(z) \le C_2(1 + |z|). \quad (2.2) $$

(H2) The function $g: \mathbb{R} \to (-\infty, +\infty]$ is lower semicontinuous with non empty domain $\text{dom } g := \{t \in \mathbb{R} : g(t) < +\infty\}$ and such that $\lim_{|t| \to \infty} \frac{g(t)}{|t|} = +\infty$. Furthermore, we assume that $g$ is upper semicontinuous on the complement of a *Lebesgue negligible* subset $D \subset \mathbb{R}$, namely

$$ \limsup_{s \to t} g(s) \le g(t), \quad \forall t \in \mathbb{R} \setminus D. $$

(H3) The source term $f$ belongs to $L^\infty(\Omega)$.

The assumption (H1) implies that $h$ is the support function of some convex compact subset of $K \subset \mathbb{R}^N$ whose indicator function coincides with $h^*$ (the Fenchel conjugate of $h$). More precisely it holds:

$$ h(z) = \sup \{z \cdot z^* : z^* \in K\} \quad , \quad K := \text{dom}(h^*). \quad (2.3) $$

Under the assumptions (H1)(H2)(H3), minimizing sequences for (2.1) turn out to be bounded in $W^{1,1}(\Omega)$, thus relatively compact in $L^1(\Omega)$. However, as it is usual when handling problems with linear growth conditions (condition (H1)), we cannot expect to obtain in the limit a solution in $W^{1,1}(\Omega)$ but merely a relaxed solution in the larger space $BV(\Omega)$ where discontinuities are allowed (see for instance [13] and also the beginning of Section A for a short background). A natural relaxed formulation of problem (2.1) is the following :

$$ (P_g) \qquad \inf \{F_g(u) : u \in BV(\Omega)\} \quad (2.4) $$

being $F_g$ defined on $BV(\Omega)$ by

$$ F_g(u) = \int_{\Omega} h(Du) + \int_{\Omega} (g(u) - f(x)u) dx + \int_{\Gamma} h((u_0 - u) \nu_{\Omega}) d\mathcal{H}^{N-1}, \quad (2.5) $$

where $\int_{\Omega} h(Du) := \int_{\Omega} h(\frac{dDu}{|Du|})d|Du|$. Note that in the notation above, we underlined the dependence with respect to $g$ since we will need to consider different choices for $g$. We emphasize that $g$ is a priori a *nonconvex* function with possibly many jumps and infinite values. This is crucial in order to link $(P_g)$ with a free boundary or a multi-phase problem. The existence of a solution for $(P_g)$ follows classically from the so called direct method of the calculus of variations.
Remark 2.1. Obviously one has inf(2.1) ≥ min($\mathcal{P}_g$). Unfortunately, the inverse inequality may fail, in particular inf(2.1) = +∞, if the range of the boundary datum $u_0$ meets an interval where $g$ is infinite. For this reason, we omit the initial formulation (2.1) and consider only the relaxed formulation ($\mathcal{P}_g$). Note, however, that the equality of the two infima can be established under the additional condition that $g$ is continuous and satisfies the growth condition $|g(t)| \le C(1 + |t|^{\frac{N}{N-1}})$ for suitable $C > 0$.

The main result of this paper (Theorem 2.4) states that, under suitable conditions, the problem ($\mathcal{P}_g$) shares the same minimum as its convexified counterpart ($\mathcal{P}_{g^{**}}$), being $g^{**}$ the convex envelope of $g$. The key argument is the following fundamental exclusion principle:

**Theorem 2.2.** Let $-\infty < a < b < +\infty$. Assume that $g_{a,b} < g < +\infty$ in $(a, b)$ where $g_{a,b}$ is the affine interpolant given by:

$$g_{a,b}(t) := \begin{cases} g(t) & \text{if } t \notin (a,b) \\ m_{a,b}(t-a) + g(a) & \text{otherwise,} \end{cases} \quad \text{and} \quad m_{a,b} := \frac{g(b) - g(a)}{b-a}.$$

Then if $u_0 \notin (a, b)$ a.e. on $\Gamma$, it holds $\inf(\mathcal{P}_g) = \inf(\mathcal{P}_{g_{a,b}})$. Moreover any solution $u$ to $(\mathcal{P}_g)$ is such that $u(x) \notin (a, b)$ for a.e. $x \in \Omega$.

The exclusion principle stated in Theorem 2.2 whose proof is postponed to Section 3 can be applied in fact to a countable union of intervals $(a, b)$ on which $g > g_{a,b}$, namely to the following set

$$Z := \{g > g^{**}\} \qquad (2.6)$$

Indeed, we can apply to $g$ the following result:

**Lemma 2.3.** Let $g: \mathbb{R} \to (-\infty, +\infty]$ be a l.s.c. non convex function whose Fenchel biconjugate $g^{**}$ is proper. Let $Z$ be the non empty set $Z := \{t \in \mathbb{R} : g(t) > g^{**}(t)\}$ and let $[m, M]$ be the closure of the interval $\{t \in \mathbb{R} : g^{**}(t) < +\infty\}$ ($-\infty \le m < M \le +\infty$). Then, we have:

(i) $Z$ is an open subset of $(m, M)$;

(ii) Assume that $Z$ is bounded. Then there exists an at most countably disjoint family of finite intervals $\{(a_i, b_i)\}_{i \in I}$ so that $Z = \cup_{i \in I}(a_i, b_i)$ and $g^{**} = g_{a_i, b_i}$ in each interval.

In view of Lemma 2.3 and of Theorem 2.2, we consider the following additional assumptions:

(H4) $Z := \{g > g^{**}\}$ is a bounded subset of $\mathbb{R}$;

(H5) $u_0 \in \mathbb{R} \setminus Z$ a.e. on $\Gamma$,

We are now in position to state our main “convex relaxation result”:

**Theorem 2.4.** Under (H1-H5), it holds:

i) $\inf(\mathcal{P}_g) = \inf(\mathcal{P}_{g^{**}})$

ii) If $u$ is a solution to $(\mathcal{P}_g)$, then $u(x) \in \mathbb{R} \setminus Z$ for a.e. $x \in \Omega$. Therefore $u$ is also a solution to $(\mathcal{P}_{g^{**}})$.

iii) Assume that $\text{dom}(g) \setminus Z$ consists of finitely many real numbers $t_1 < \dots < t_k$. Then solutions to $(\mathcal{P}_g)$ are all of the form:

$$u = \sum_{i=1}^{k} t_i \mathbb{1}_{A_i}, \quad \text{where } \{A_i, 1 \le i \le k\} \text{ is a Borel partition of } \Omega. \qquad (2.7)$$
Proof of Theorem 2.4. Let $u$ be a solution to $(\mathcal{P}_g)$ (the existence of $u$ is provided by Lemma ??). Let $\{(a_i, b_i)\}_{i \in I}$ be the family of intervals given by Lemma 2.3. Then, for every $i \in I$ it holds $u_0 \notin (a_i, b_i)$ a.e. whereas $g > g_{a_i, b_i}$. Thus, by applying Theorem 2.2, we infer that $u \notin (a_i, b_i)$ a.e. Since $I$ is at most countable, we deduce that $u \notin Z$ a.e.

i) and ii) is a straightforward consequence of the exclusion principle and Lemma 2.3 where we apply Theorem 2.2 on each subinterval $(a_i, b_i)$ of $Z = \{g^{**} < g\} = \cup_i (a_i, b_i)$. iii) Assume that $\text{dom}(g) \setminus Z$ consists of finitely many reals numbers $t_1 < \dots < t_k$, thus $Z = \cup_{i=1}^{k-1} (t_i, t_{i+1})$. According to the exclusion principle, if $u$ is a solution to $(\mathcal{P}_g)$ then $u$ only ranges in $\{t_1, \dots, t_k\}$. Let $A_i := \{u = t_i\}$ for $i = 1, \dots, k$. Then, it holds

$$u = \sum_{i=1}^{k} t_i \mathbf{1}_{A_i}.$$

□

Remark 2.5. As shown in the assertion ii), we have merely the inclusion $\text{Argmin}_{\mathcal{P}_g} \subset \text{Argmin}_{\mathcal{P}_{g^{**}}}$. Clearly, since $\text{Argmin}_{\mathcal{P}_{g^{**}}}$ is a convex subset, the equality $\text{Argmin}_{\mathcal{P}_g} = \text{Argmin}_{\mathcal{P}_{g^{**}}}$ is untrue in the case where $(\mathcal{P}_g)$ has several solutions. However the equality holds in the case where $\mathcal{P}_{g^{**}}$ admits a unique solution. This uniqueness issue for $(\mathcal{P}_{g^{**}})$ seems to be a very difficult task since $g^{**}$ is not strictly convex (at all on $Z$) and we refer for that to the recent paper for a uniqueness result in a similar context [4]. We believe that uniqueness for $(\mathcal{P}_g)$ implies uniqueness for $(\mathcal{P}_{g^{**}})$ although we are not yet able to prove it. Actually in many examples of multiphase problems as presented in Section 5, we observed a very good convergence of the numerical approximation of $(\mathcal{P}_{g^{**}})$ to a solution of the form (2.7). Therefore we believe that uniqueness holds generically with respect to the boundary data and the parameters entering in the definition of $g$, $\Omega$.

Remark 2.6. We stress the fact that the assumption (H1) (homogeneity of integrand $h$) and the assumption (H5) ($u_0 \in Z$) are necessary in order to ensure the validity of Theorem 2.4. Related to this, two illuminating counterexamples are detailed in Section 4.

Proof of Lemma 2.3. (i) Obviously one has $Z \subset [m, M]$ since $g(t) = g^{**}(t) = +\infty$ for $t \notin [m, M]$. On the other hand, as we are on the real line, the convex envelope $\text{co} g$ of $g$ is given by:

$$g^{**}(t) = (\text{co} g)(t) = \inf \{\theta g(t') + (1-\theta)g(t'') : \theta \in [0,1], t', t'' \in [m, M], \theta t' + (1-\theta)t'' = t\}$$

In particular, $g^{**}$ agrees on the extreme points of $[m, M]$, hence $m, M \notin Z$ and $Z \subset (m, M)$.

We show now that $Z^c := \mathbb{R} \setminus Z$ is a closed subset of $\mathbb{R}$. Let $t_n \in Z^c$ be such that $t_n \to t$. We need to check that $t \in Z^c$ which is already known if $t \notin (m, M)$. If $t \in (m, M)$, we exploit the continuity of the convex function $g^{**}$ on the open interval $(m, M)$ together with the lower semicontinuity of $g$:

$$g^{**}(t) = \lim_{n} g^{**}(t_n) \geq \inf_{n} g(t_n) \geq g(t).$$

Thus $g^{**}(t) = g(t)$ and $t \in Z^c$.

(ii) The bounded open subset of $Z \subset (m, M)$ can be written as $Z = \cup_{i \in I} (a_i, b_i)$ where $I$ is finite or countable and $(a_i, b_i) \subset (m, M)$ (thus $g(a_i), g(b_i)$ are finite). It remains to show that $g^{**}$ is affine in each interval $(a_i, b_i)$. Let us denote by $\hat{g}$ the convex function which coincides with $g^{**}$ on $Z^c$ and with the affine interpolant of $g^{**}$ on each interval $[a_i, b_i]$. Clearly one has:

$$g^{**} \leq \hat{g}, \quad g^{**} = \hat{g} = g \text{ in } Z^c, \quad \partial \hat{g}(t) = \{(\hat{g})'(t)\} = \{m_i\} \quad \forall t \in (a_i, b_i),$$

where $m_i := \frac{g(b_i) - g(a_i)}{b_i - a_i}$ denotes the slope of $\hat{g}$ on $[a_i, b_i]$. We are done if we can show the inequality $g \geq \hat{g}$ since then $g^{**} \geq \hat{g}$ by taking the convex envelope, thus $g^{**} = \hat{g}$. We observe that the function $\varphi := g - \hat{g}$ vanishes on $Z^c$ and is l.s.c. (since $g$ is l.s.c. and $\hat{g}$ is continuous on $(m, M)$). Assume that $\inf \varphi < 0$. Then, recalling that $Z$ is bounded, $\varphi$ achieves its global minimum at some $\bar{t} \in Z$. Let
$i \in I$ such that $\bar{t} \in (a_i, b_i)$. Then $\hat{g}(t) \ge \hat{g}(\bar{t}) + m_i(t - \bar{t})$ whereas $\varphi(t) \ge \varphi(\bar{t})$ for all $t \in \mathbb{R}$. So far we have obtained that $g(t) \ge g(\bar{t}) + m_i(t - \bar{t})$ thus $m_i \in \partial g(\bar{t})$. This is impossible since $g^{**}(\bar{t}) < g(\bar{t})$. $\square$

### 3. PROOF OF THE EXCLUSION PRINCIPLE

This Section is devoted to the proof of Theorem 2.2. The key tool is a duality argument introduced in [7] for super linear growth functionals that we adapt to our case. To shorten the presentation, we will treat only the case where $\text{dom}(g) \subset I := [m, M]$ where $-\infty < m < M < +\infty$, thus reducing the class of competitors to a subset of $BV(\Omega; I)$, that is

$$ (\mathcal{P}_g) \qquad \min \{F_g(u) : u \in BV(\Omega; I)\}. \qquad (3.1) $$

Notice that if $f = 0$ and $u_0 \in (m, M)$ a.e. on $\Gamma$, such a reduction can be obtained by using a truncation argument provided $g$ satisfies the milder condition $g(t) > g(m)$ on $(-\infty, m)$ and $g(t) > g(M)$ on $(M, +\infty)$.

**Duality result** Following the strategy developed in [9], [7], we are led to consider a dual problem of $(\mathcal{P}_g)$ in the form of a constrained optimal flow problem in higher dimension. This problem involves divergence free vector fields $\sigma := (\sigma^x, \sigma^t) \in L^\infty(\Omega \times I; \mathbb{R}^{N+1})$ and reads as follows:

$$ (\mathcal{P}_g^*) \qquad \sup \{L(\sigma) : \sigma \in \mathcal{B}_g\} \qquad (3.2) $$

Here the affine function $L(\sigma)$ is defined by

$$ L(\sigma) := - \int_{\Omega} \sigma^t(x, m) dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_{\Omega} \mathbf{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt \qquad (3.3) $$

where $\sigma^x \cdot \nu_{\Omega}$ denotes the normal trace of $\sigma$ on $\partial\Omega \times I$, and $\mathbf{1}_{u_0}(x, t) = \begin{cases} 1 & \text{if } t \le u_0(x) \\ 0 & \text{otherwise} \end{cases}$.

The vector field $\sigma$ runs over the class $\mathcal{B}_g \subset L^\infty(\Omega \times I; \mathbb{R}^{N+1})$ defined by:

* $\sigma$ satisfies the equations:

$$ -\operatorname{div} \sigma = 0 \quad \text{in } \Omega \times I, \qquad (3.4) $$

$$ \sigma^x \cdot \nu_{\Omega} = 0 \quad \text{on } (\partial\Omega \setminus \Gamma) \times I \qquad (3.5) $$

* $\sigma$ belongs to the convex subset $\mathcal{K}_g$ associated with the pointwise constraints:

$$ \sigma^x(x, t) \in K \quad \text{a.e. in } \Omega \times I \quad (K \text{ defined in (2.3)}) \qquad (3.6) $$

$$ \sigma^t(x, t) + g(t) - t f(x) \ge 0 \quad \forall t \in I \text{ and for } \mathcal{L}^N\text{-a.e. } x \in \Omega, \qquad (3.7) $$

As the vector field $\sigma$ is bounded and divergence free, its normal trace $\sigma^t(\cdot, a)$ is well defined $\mathcal{L}^N$-a.e. $x \in \Omega$ on every horizontal interface $t = a$. Its normal trace is defined as well on the lateral boundary $\partial\Omega \times I$. We denote it by $\sigma^x \cdot \nu_{\Omega}$ and we may apply the generalized Green's formula given in (A.7). As a consequence the affine function $L(\sigma)$ appearing (3.3) can be alternatively be written as

$$ L(\sigma) := - \int_{\Omega} \sigma^t(x, a) dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_{\Omega} (\mathbf{1}_{u_0} - \mathbf{1}_a) d\mathcal{H}^{N-1} \otimes dt \qquad (3.8) $$

for any real $a \in [m, M]$, being $\mathbf{1}_a$ the counterpart of $\mathbf{1}_{u_0}$ for a constant function.

Note that in [7], the interface condition (3.7) was imposed only for $t \in D$ where $D$ is the discontinuity subset of $g$ appearing in assumption (H3). A key point is the following no duality gap statement:

**Theorem 3.1.** It holds $\inf(\mathcal{P}_g) = \sup(\mathcal{P}_g^*)$.
The proof of Theorem 3.1, postponed to the Appendix, is quite technical and relies on an adaptation of the calibration method developed in [7]. Thanks to the duality result of Theorem 3.1, we are able to prove the exclusion principle stated in Theorem 2.2. This is the purpose of the rest of this section.

## Proof of Theorem 2.2

In order to prove the theorem, we need only to establish the following inequality:

$$ \inf(\mathcal{P}_{g_{a,b}}) \leq \inf(\mathcal{P}_g). \qquad (3.9) $$

Indeed, the converse inequality is trivial since $g \geq g_{a,b}$. Then, if $u$ is a minimizer for $(\mathcal{P}_g)$, we will deduce that

$$ F_{g_{a,b}}(u) \leq F_g(u) = \inf(\mathcal{P}_g) = \inf(\mathcal{P}_{g_{a,b}}), $$

thus $u \in \text{Argmin}(\mathcal{P}_{g_{a,b}})$ and $F_{g_{a,b}}(u) = F_g(u)$. Therefore, since $g > g_{a,b}$ on $(a, b)$ (while $g = g_{a,b}$ outside), we will conclude that $u \notin (a, b)$ a.e. in $\Omega$.

For proving Claim (3.9), we exploit the duality result presented above considering the dual problem $(\mathcal{P}_g^*)$ in the form given by (3.2) where the affine function $L(\sigma)$ is written in the form (3.8), that is

$$ L(\sigma) = - \int_{\Omega} \sigma^t(x, a) dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_{\Omega} (\mathbf{1}_{u_0} - \mathbf{1}_a) d\mathcal{H}^{N-1} \otimes dt . $$

In a similar way, we consider the dual problem associated with $(\mathcal{P}_{g_{a,b}})$ namely

$$ (\mathcal{P}_{g_{a,b}}^*) \qquad \sup \left\{ L(\sigma) : \sigma \in \mathcal{B}_{g_{a,b}} \right\}, $$

where $\mathcal{B}_{g_{a,b}}$ is the counterpart of the subset $\mathcal{B}_g$ where we keep conditions (3.4), (3.5), (3.6) while we substitute (3.7) with the condition

$$ \sigma^t(x, t) + g_{a,b}(t) - t f(x) \geq 0 \quad \forall t \in I \text{ and for } \mathcal{L}^N\text{-a.e. } x \in \Omega. \qquad (3.10) $$

By applying Theorem 3.1 to $(\mathcal{P}_g)$ and $(\mathcal{P}_{g_{a,b}})$ respectively, we infer that proving (3.9) amounts to show that $\sup(\mathcal{P}_g^*) \leq \sup(\mathcal{P}_{g_{a,b}})$, that is:

$$ L(\sigma) \leq \sup(\mathcal{P}_{g_{a,b}}^*) , \quad \text{for all } \sigma \in \mathcal{B}_g \qquad (3.11) $$

Let us take such a $\sigma \in \mathcal{B}_g$ and modify it in the strip $\Omega \times (a, b)$ so that the new vector field $\sigma_{a,b}$ is admissible for $(\mathcal{P}_{g_{a,b}}^*)$ while $L(\sigma) = L(\sigma_{a,b})$. This is done in a straightforward way by averaging the horizontal component $\sigma^x(x, \cdot)$ with respect to $t \in (a, b)$ while the vertical component $\sigma^t(x, \cdot)$ is changed into an affine function between the two interfaces $t = a$ and $t = b$ (see Figure 1). Precisely, let us set

$$ \theta(t) := \inf \left\{ \left( \frac{t-a}{b-a} \right)_+ , 1 \right\} , \quad q(x) := \frac{1}{b-a} \int_a^b \sigma^x(x, s) ds , \qquad (3.12) $$

and define

$$ \sigma_{a,b}(x,t) = \begin{cases} (q(x), (1-\theta(t))\sigma^t(x,a) + \theta(t)\sigma^t(x,b)) & \text{if } t \in [a,b], \\ \sigma(x,t) & \text{if } t \notin (a,b). \end{cases} \qquad (3.13) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 1. The field $\sigma_{a,b}$ defined by (3.13).</figcaption></figure>

By our assumption (H5), it is straightforward that

$$L(\sigma) = L(\sigma_{a,b}) \qquad (3.14)$$

Indeed, as $u_0(x) \notin (a, b)$ for a.e. $x \in \Gamma$, the function $(\mathbf{1}_{u_0} - \mathbf{1}_a)(x, \cdot)$ is constant in $(a, b)$, while by (3.13) the average of $\sigma^x - \sigma_{a,b}^x$ vanishes (see (3.13)). Therefore, as $\sigma = \sigma_{a,b}$ in $\Omega \times (I \setminus (a, b))$, we have:

$$L(\sigma) - L(\sigma_{a,b}) = \int_{\Gamma \times (a,b)} (\sigma^x - \sigma_{a,b}^x) \cdot \nu_{\Omega} (\mathbf{1}_{u_0} - \mathbf{1}_a) d\mathcal{H}^{N-1} \otimes dt = 0$$

Now we show that $\sigma_{a,b}$ belongs to $\mathcal{B}_{g_{a,b}}$ by checking the conditions (3.4), (3.5), (3.6) and (3.10).

• The relation (3.10) for $\sigma_{a,b}$ is deduced directly by taking the affine interpolant (between $t = a$ and $t = b$) of the left hand side of (3.7)).

• Since $K$ is a compact convex subset, it is clear that $\sigma^x(x, t) \in K$ a.e. implies that the average $q(x)$ given in (3.12) still satisfies $q(x) \in K$ a.e. in $\Omega$. Thus $\sigma_{a,b}$ satisfies (3.6).

• Let us check the condition (3.4). By applying generalized Green's formula on $\Omega \times (a, b)$ to $\sigma$ and a smooth function $\varphi(x, t) = \psi(x)$ with $\psi \in \mathcal{D}(\Omega)$, we obtain:

$$\int_{\Omega} (\sigma^t(x, b) - \sigma^t(x, a))\psi(x) dx = \int_{\Omega \times (a,b)} \sigma \cdot \nabla \varphi dx dt = (b-a) \int_{\Omega} q(x) \cdot \nabla \psi(x) dx,$$

where $q$ is defined by (3.12). As $\psi$ is arbitrary, we deduce that $q$ satisfies the equation

$$-\mathrm{div}_x q = \frac{\sigma^t(x, b) - \sigma^t(x, a)}{b-a} \quad \text{in the distributional sense on } \Omega.$$

In view of (3.13), it follows that $\sigma_{a,b}$ is divergence free on the open set $\Omega \times (a, b)$, thus on all $\Omega \times \mathbb{R}$ since the normal trace $\sigma^t(x, \cdot)$ does not jump on the interfaces $t = a$ and $t = b$.

• Let us check now that $\sigma_{a,b}$ still satisfies the normal trace condition (3.5). This is clearly the case on subset $(\partial\Omega \setminus \Gamma) \times (I \setminus (a, b))$. And for $\mathcal{H}^N$-a.e. $(x, t) \in (\partial\Omega \setminus \Gamma) \times (a, b)$, we have

$$(\sigma_{a,b})^x(x, t) \cdot \nu_{\Omega}(x) = q(x) \cdot \nu_{\Omega}(x) = \frac{1}{b-a} \int_a^b \sigma^x(x, s) \cdot \nu_{\Omega}(x) ds = 0.$$

Summarizing we have shown that $\sigma_{a,b}$ is admissible for $(\mathcal{P}_{g_{a,b}}^*)$. Thus, recalling (3.14), we conclude that $L(\sigma) = L(\sigma_{a,b}) \le \sup(\mathcal{P}_{g_{a,b}}^*)$. This is the claim (3.11) and the proof of Theorem 2.2 is finished. $\square$
4. COUNTEREXAMPLES

In this section, we discuss about the necessity of the one homogeneity assumption on h ( i.e. (H1)) and of the condition (H5) on the boundary data $u_0$. We are going to show that the conclusions of Theorem 2.4 fail, in particular that $\inf(\mathcal{P}_{g_{\lambda}^{**}}) < \inf(\mathcal{P}_{g_{\lambda}})$ can occur, if one of the two assumptions is violated. In the two following examples taken in dimension $N = 2$, we choose the non convex $g$ to be $g = g_{\lambda}$ where

$$g_{\lambda}(t) = \begin{cases} 0 & \text{if } t \in [0, 1) \\ -\lambda & \text{if } t = 1 \\ +\infty & \text{if } t \notin [0, 1] \end{cases}$$

and $\lambda$ is a positive real parameter.

In that case, the subset $Z$ defined in (2.3) is given by $Z = (0, 1)$ and the convex envelop reads

$$g_{\lambda}^{**}(t) = \begin{cases} -\lambda t & \text{if } t \in [0, 1] \\ +\infty & \text{if } t \notin [0, 1] \end{cases} \qquad (4.1)$$

Since $\text{dom}(g_{\lambda}) \setminus Z := \{0, 1\}$, the application of Theorem 2.4 should lead to binary solutions $u$ i. e. ranging in $\{0, 1\}$.

**Example 4.1.** In this example the assumption (H5) is not satisfied. Let $\Omega$ be the unit square in $\mathbb{R}^2$, i.e. $\Omega = (0, 1)^2$. We take the boundary data $u_0 = \frac{1}{2}$ and consider the variational problem

$$(\mathcal{P}_{g_{\lambda}}) \qquad \inf \left\{ \int_{\Omega} |Du| + \int_{\partial\Omega} |u - \frac{1}{2}| d\mathcal{H}^1 + \int_{\Omega} g_{\lambda}(u) dx : u \in \text{BV}(\Omega) \right\}.$$

Let $u$ be optimal for problem $(\mathcal{P}_{g_{\lambda}})$. Then it satisfies $\frac{1}{2} \le u \le 1$. Indeed, if we set $\tilde{u} = u \lor \frac{1}{2}$, then

$$\int_{\Omega} |\nabla u| dx - \lambda |u| = 1| \ge \int_{\Omega} |\nabla \tilde{u}| dx - \lambda |\tilde{u}| = 1|$$

so that equality holds and $\tilde{u} = u$ a.e. We can then apply to $u$ the following identity

$$\int_{\Omega} |Dw| + \int_{\partial\Omega} |w - \frac{1}{2}| d\mathcal{H}^1 = \int_{\frac{1}{2}}^{1} P(w > t) dt \quad \text{for all } w \in \text{BV}\left(\Omega; \left[\frac{1}{2}, 1\right]\right) \quad (4.2)$$

(which is derived from the coarea formula applied to the function $v$ defined by $v = w - \frac{1}{2}$ in $\Omega$ and $v = 0$ in $\mathbb{R}^2 \setminus \Omega$).

Now by taking into account $\{u = 1\} \subset \{u > t\}$ for all $t \in [\frac{1}{2}, 1]$, we infer that

$$\min(\mathcal{P}_{g_{\lambda}}) = \int_{\Omega} |Du| + \int_{\partial\Omega} |u - \frac{1}{2}| d\mathcal{H}^1 - \lambda |u| = 1| \ge \int_{\frac{1}{2}}^{1} (P(u > t) - 2\lambda|u > t|) dt. \quad (4.3)$$

Let $h_{\Omega}$ denote the Cheeger constant of $\Omega$, that is

$$h_{\Omega} := \inf_{A \subset \Omega} \frac{P(A)}{|A|}.$$

Then, if $\lambda \le \frac{h_{\Omega}}{2}$, it holds

$$\int_{\frac{1}{2}}^{1} (P(u > t) - 2\lambda|u > t|) dt \ge \int_{\frac{1}{2}}^{1} (P(u > t) - h_{\Omega}|u > t|) dt \ge 0.$$

Therefore, from (4.3), we conclude that

$$\inf(\mathcal{P}_{g_{\lambda}}) \ge 0 \quad \text{whenever } \lambda \le \frac{h_{\Omega}}{2}.$$
Now if we consider the convex problem ($P_{g_{\lambda}^{**}}$) obtained by substituting $g_{\lambda}$ with $g_{\lambda}^{**}$, by taking $u \equiv \frac{1}{2}$ as a competitor, we get

$$ \inf(P_{g_{\lambda}^{**}}) \leq -\frac{\lambda}{2}|\Omega| = -\frac{\lambda}{2} $$

Therefore we conclude that $\inf(P_{g_{\lambda}^{**}}) < \inf(P_{g_{\lambda}})$ holds whenever $\lambda \leq \frac{h_{\Omega}}{2}$.

Let us notice that $u \equiv \frac{1}{2}$ is the unique solution to problem ($P_{g_{\lambda}^{**}}$) for all $\lambda < h_{\Omega}$. Indeed, as for ($P_g$), any solution $u$ to ($P_{g_{\lambda}^{**}}$) belongs to $BV(\Omega; [\frac{1}{2}, 1])$ so that by applying (4.2), we get

$$ \int_{\Omega} |Du| + \int_{\partial\Omega} |u - \frac{1}{2}| d\mathcal{H}^1 - \lambda \int_{\Omega} u \geq \int_{\frac{1}{2}}^{1} (P(u > t) - \lambda|u > t|) dt - \frac{\lambda}{2} |\Omega|. $$

Since $\lambda < h_{\Omega}$, the integral in the right hand side is strictly positive unless $u = \frac{1}{2}$ a.e.

**Example 4.2.** In the following, taking the same $g_{\lambda}$ as before, we consider now the case where $h$ is not homogeneous. More precisely we consider two different integrands

$$ h(z) = \sqrt{1 + |z|^2}, \qquad h_0(z) = 1 + |z|. $$

Correspondingly, we consider the infimum problems:

$$ \beta(\lambda) = \inf \left\{ \int_{\Omega} \sqrt{1 + |\nabla u|^2} dx + \int_{\Omega} |D_s u| + \int_{\partial\Omega} |u| d\mathcal{H}^1 - \lambda \left| \{u = 1\} \right| : u \in BV(\Omega) \right\}, $$

$$ \beta_0(\lambda) = \inf \left\{ \int_{\Omega} |Du| + \int_{\partial\Omega} |u| d\mathcal{H}^1 - \lambda \left| \{u = 1\} \right| : u \in BV(\Omega) \right\}. $$

Let us observe that $\beta(\lambda)$ is the relaxed version of a quite challenging problem minimal surface problem with free boundary (see for instance [11, 10, 3, 13, 12, 14, 15]) in which the free boundary consists of the boundary of the unknown set $\{u = 1\}$, namely

$$ \inf_{u \in W_0^{1,1}(\Omega)} \left\{ \int_{\Omega} \sqrt{1 + |\nabla u|^2} dx - \lambda |u = 1| \right\}. $$

The infimum problem corresponding to $\beta_0(\lambda)$ is much simpler. In particular, as it fits to the assumptions (H1)-(H5) of our convexification result (Theorem 2.4), its solutions only take value in $\{0, 1\}$ and by substituting $u$ with $\mathbb{1}_A$, we obtain the following Cheeger's type problem:

$$ \beta_0(\lambda) = |\Omega| + \inf \{P(A) - \lambda|A| : A \subset \Omega\}. $$

Suppose that Theorem 2.4 holds also for problem $\beta(\lambda)$, then its solutions would range in $\{0, 1\}$ as well. Then observing that the two functionals to be minimized in $\beta(\lambda)$ and $\beta_0(\lambda)$ coincide for binary functions $u = \mathbb{1}_A$, we infer that the equality $\beta(\lambda) = \beta_0(\lambda)$ would hold for every $\lambda \geq 0$. We are going to show that such an equality is untrue in general, even in simple case of a radial geometry.

To that aim we consider $\Omega = \{x \in \mathbb{R}^2 : |x| < R\}$. In this case, we obtain that the unique solution of problem $\beta_0(\lambda)$ is $u_0 \equiv 0$ for $\lambda < h_{\Omega}$ where $h_{\Omega} = \frac{2}{R}$ is the Cheeger constant of $\Omega$. In contrast the unique solution is $u_1 \equiv 1$ for $\lambda > h_{\Omega}$. For the precise value $\lambda = h_{\Omega}$, we obtain the coexistence of the two solutions $u_0, u_1$. Accordingly,

$$ \beta_0(\lambda) = \begin{cases} |\Omega| & \text{if } \lambda \leq h_{\Omega} \\ P(\Omega) + (1 - \lambda)|\Omega| & \text{if } \lambda > h_{\Omega}. \end{cases} $$

Let us now turn to the determination of the infimum $\beta(\lambda)$. By a rearrangement argument, we can prove that solutions are all radial of the form $u(x) = \varphi(\frac{|x|}{R})$ being $\varphi(t) : [0, 1] \to [0, 1]$ monotone
non increasing. The plateau {$u=1$} is associated with an interval $t \in [0, \rho]$ for a suitable value of $\rho \in [0, 1]$ to be determined. For such a plateau, the minimal surface problem reads

$$J(\rho) := \inf_{\substack{\varphi(\rho)=1 \\ \varphi(1)=0}} I(\varphi), \quad I(\varphi) := R \int_{\rho}^{1} \sqrt{R^2 + \varphi'^2} t dt. \qquad (4.4)$$

The first integral of Euler equation for this minimization problem reads

$$\frac{t\varphi'}{\sqrt{R^2 + \varphi'^2}} = \mu \qquad (4.5)$$

for some constant $\mu$. As $\varphi(1) = 0$, we are led to the explicit form

$$\varphi(t) = K(\mu, t), \quad K(\mu, t) := \mu R \log \left( \frac{1 + \sqrt{1 - \mu^2}}{t + \sqrt{t^2 - \mu^2}} \right), \qquad (4.6)$$

provided we can find $\mu \in [0, \rho]$ such that $\varphi(\rho) = K(\mu, \rho) = 1$. In fact, the function $\mu \mapsto K(\mu, \rho)$ is strictly increasing on $[0, \rho]$ so that such a $\mu$ is unique and exists if and only if

$$1 \le K(\rho, \rho) = \rho R \log \left( \frac{1 + \sqrt{1 - \rho^2}}{\rho} \right). \qquad (4.7)$$

If the inequality above is strict, then we obtain that the solution to (4.4) is unique and smooth. In the limit case where (4.7) is an equality, one has $\mu = \rho$ and $\varphi'(\rho) = +\infty$. If $K(\rho, \rho) < 1$, then no solution to (4.4) exists but merely a (unique) relaxed solution which exhibits a jump at $t = \rho$ of amplitude $1 - K(\rho, \rho)$. In all the cases, we have determined, in term of parameter $\rho \in [0, 1]$, an optimal radial configuration whose plateau {$u=1$} agrees with the disk $B(0, \rho)$. Its total energy is given by

$$E(\rho) := 2\pi J(\rho) + (1 - \lambda)\pi\rho^2 R^2 .$$

In order to minimize $E(\rho)$ on interval $[0, 1]$, we introduce

$$\bar{\mu}(\rho) := \sup_{0 \le \mu \le \rho} \{\mu : K(\mu, \rho) \le 1\}.$$

It is easy to check that for every $\rho \in [0, 1]$ such that if $K(\rho, \rho) < 1$ then $\bar{\mu}(\rho) = \rho$. Otherwise, $\bar{\mu}(\rho)$ is the unique solution of equation $K(\mu, \rho) = 1$. After a straightforward computation and exploiting (4.5), we obtain

$$J(\rho) = R^2 \int_{\rho}^{1} \frac{t^2}{\sqrt{t^2 - \bar{\mu}(\rho)^2}} dt + \rho R (1 - K(\bar{\mu}(\rho), \rho)).$$

Thus, noticing that $K(\bar{\mu}(\rho), \rho) \le 1$, we are led to:

$$E(\rho) = \pi R^2 \left( \sqrt{1 - \bar{\mu}(\rho)^2} - \rho \sqrt{\rho^2 - \bar{\mu}(\rho)^2} + \bar{\mu}(\rho)^2 \log \frac{1 + \sqrt{1 - \bar{\mu}(\rho)^2}}{\rho + \sqrt{\rho - \bar{\mu}(\rho)^2}} + \frac{2\rho(1 - K(\bar{\mu}(\rho), \rho))}{R} + (1 - \lambda)\rho^2 \right)$$

Summarizing, an optimal $\bar{\rho}$ for $E(\rho)$ will give a radial function $\bar{u}(x) = \bar{\varphi}(\frac{|x|}{R}) = K(\bar{\mu}(\bar{\rho}), \frac{|x|}{R})$ (defined in (4.6)) which minimizes the relaxed problem associated with $\beta(\lambda)$. This solution is continuous if $K(\bar{\rho}, \bar{\rho}) = 1$ and otherwise exhibits a jump of amplitude $1 - K(\bar{\rho}, \bar{\rho})$ before reaching the value 1 on the plateau. This is illustrated in Figure 2.
<figure><img src="image_5.png" /><figcaption>FIGURE 2. Illustration for an optimal $u$ and optimality conditions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 3. Numerical computation of functions $\beta_0(\lambda)$, $\beta(\lambda)$, critical values of $\lambda$: $\lambda_0(R)$, $\lambda_1(R)$, and Cheeger constant $h_\Omega = \frac{2}{R}$.</figcaption></figure>

The minimization of $E(\rho)$ is performed by using Matlab for different values of $R$ and $\lambda$. It turns out that:

* For $R > 1$, the exists constants $\lambda_0(R) < h_\Omega < \lambda_1(R)$ where $h_\Omega = \frac{2}{R}$ (Cheeger constant of $\Omega$) and such that $\beta(\lambda) < \beta_0(\lambda)$ for every $\lambda \in (\lambda_0, \lambda_1)$. Solutions to $\beta(\lambda)$ do or do not have jump depending on the position of $\lambda$ with respect to 1: $[u] = 0$ if $\lambda \le 1$ and $[u] > 0$ if $\lambda > 1$ (see Figure 2 and the middle subfigure of Figure 3). In contrast it holds $\beta(\lambda) = \beta_0(\lambda)$ for $\lambda \le \lambda_0$ (in this case $u \equiv 0$ solves both problems) or for $\lambda \ge \lambda_1$ (in that case $u \equiv 1$ is a solution).

* For small balls $R \le 1$, we have the equality $\beta(\lambda) = \beta_0(\lambda)$ (and $\lambda_0 = \lambda_1 = \frac{2}{R}$ which is shown in the right hand side of Figure 3).

## 5. APPLICATION TO FREE BOUNDARY AND MULTIPHASE PROBLEMS

In this Section we consider several multiphase problems on the unit square $\Omega = (0, 1)^2$ of $\mathbb{R}^2$. The boundary data $u_0$ will take two values $u_0 = 0$ on $\Gamma_0$ and $u_0 = 1$ on $\Gamma_1$, with $\{\Gamma_0, \Gamma_1\}$ being a partition of $\Gamma := \partial\Omega$. We first treat a three phases problem (in Subsection 5.1) and then a four phases problem (in Subsection 5.2).

### 5.1. A three phases problem. We consider the variational problem

$$ \min \left\{ \int_{\Omega} |Du| + \int_{\Gamma_0} |u| d\mathcal{H}^1 + \int_{\Gamma_1} |1-u| d\mathcal{H}^1 + \lambda \left| \left\{ u \neq \frac{1}{2} \right\} \right| : u \in BV(\Omega; [0, 1]) \right\}. \quad (5.1) $$
Here the term $\lambda|\{u \neq \frac{1}{2}\}|$ appearing in the energy can be written as $\int_{\Omega} g(u)dx$ where the non convex function $g(t)$ is given by

$$g(t) = \begin{cases} 0 & \text{if } t = \frac{1}{2} \\ \lambda & \text{if } t \in [0, \frac{1}{2}) \cup (\frac{1}{2}, 1] \\ +\infty & \text{otherwise.} \end{cases} \qquad (5.2)$$

Obviously, the convex envelope of $g(t)$ reads

$$g^{**}(t) = \begin{cases} 2\lambda|t - \frac{1}{2}| & \text{if } t \in [0, 1] \\ +\infty & \text{otherwise.} \end{cases}$$

In Figure 4 below, we describe the two function $g$, $g^{**}$ and $\text{epi}(g) := \{(t, \alpha) \in \mathbb{R}^2 : g(t) \leq \alpha\}$ the epigraph of $g$.

<figure><img src="image_7.png" /><figcaption>FIGURE 4. (a) Representation of $g$ and of $g^{**}$ (in blue). (b) Representation of $g^*$. The epigraphs $\text{epi}(g)$ and $\text{epi}(g^*)$ are filled in dashed.</figcaption></figure>

Notice that in this case the set $Z = \{g^{**} < g\}$ is the complement of $\{0, \frac{1}{2}, 1\}$ in the interval $[0, 1]$. In view of Theorem 2.4, a solution to the relaxed problem in $BV(\Omega)$ associated with (5.1) takes only three values represented by the phases $A_0 := \{u = 0\}$, $A_{\frac{1}{2}} := \{u = \frac{1}{2}\}$, $A_1 := \{u = 1\}$. The convexified problem that we are going to solve in order to recover these three phases reads

$$\min \left\{ \int_{\Omega} |Du| + \int_{\Gamma_0} |u| d\mathcal{H}^1 + \int_{\Gamma_1} |1 - u| d\mathcal{H}^1 + 2\lambda \int_{\Omega} |u - \frac{1}{2}| dx : u \in BV(\Omega; [0, 1]) \right\}. \quad (5.3)$$

If the solution $u$ to (5.3) is unique, then (see Theorem 2.4) it is of the form $u = \frac{1}{2} \mathbf{1}_{A_{\frac{1}{2}}} + \mathbf{1}_{A_1}$ which allows to determine the searched optimal partition $\{A_0, A_{\frac{1}{2}}, A_1\}$ of $\Omega$. This scenario is confirmed by numerical simulations, namely in the particular case described in Figure 6 (see also the 4 phases variant of (5.3) in Figures 11).

## Dual problem and interface conditions.

The classical dual problem to $(\mathcal{P}_{g^{**}})$ reads:

$$(Q) \quad \sup \left\{ \int_{\Gamma_1} q \cdot \nu_{\Omega} d\mathcal{H}^{N-1} - \int_{\Omega} g^*(\text{div } q)dx : q \in L^{\infty}(\Omega; \mathbb{R}^N), |q| \le 1 \right\}$$

Let $\bar{q}$ solve $(Q)$. Then a function $\bar{u}(x)$ taking values in $\{0, \frac{1}{2}, 1\}$ for a.e. $x \in \Omega$ solves problem $(\mathcal{P}_g)$ (then also $(\mathcal{P}_{g^{**}})$) if and only if

$$\bar{q} = \nu_{\bar{u}} \text{ on } S_{\bar{u}}, \quad \text{div} \bar{q} \in \partial g^{**}(\bar{u}) \text{ a.e. in } \Omega$$
where $S_{\bar{u}}$ denotes the jump set of $\bar{u}$.

We now give optimality conditions for minimal partitions, which are called recovering interface conditions. Euler equation for shape variations of problem $(\mathcal{P}_{g^{**}})$ is written in the form of conservation law (see [8])

$$ \mathrm{div} \mathcal{A} = 0 \quad \text{in } \Omega $$

where $\mathcal{A}$ is a tensor defined on the product space $\mathrm{Argmin}(\mathcal{P}_{g^{**}}) \times \mathrm{Argmax}(\mathcal{Q})$:

$$ \mathcal{A}(u, q) := (|Du| + g^{**}(u))I_2 - Du \otimes q. $$

Let $(\bar{u}, \bar{q})$ be an element of $\mathrm{Argmin}(\mathcal{P}_{g^{**}}) \times \mathrm{Argmax}(\mathcal{Q})$. The vector field $\bar{q}$ identifies to $D\bar{u}/|D\bar{u}|$, and the function $\bar{u}$ only takes values in $\{0, 1/2, 1\}$ so that $g^{**}(\bar{u}) = g(\bar{u})$. Thus, we have

$$ \mathcal{A}(\bar{u}, \bar{q}) = (I_2 - \nu_{\bar{u}} \otimes \nu_{\bar{u}}) [\bar{u}] \delta_{S_{\bar{u}}} + g(\bar{u}) \mathbf{1}_{\Omega} I_2. $$

In particular, on the interface $S_{\bar{u}} \cap \Omega$ the conservation law reduces to

$$ [\bar{u}]\kappa_{S_{\bar{u}}} + [g(\bar{u})] = 0 \qquad (5.4) $$

provided $(I_2 - \nu_{\bar{u}} \otimes \nu_{\bar{u}})\delta_{S_{\bar{u}}}$ being the tangential projection of the curve vector measure $\nu_{\bar{u}}\delta_{S_{\bar{u}}}$ and the distribution $\kappa_{S_{\bar{u}}}\nu_{\bar{u}} = \mathrm{div}((I_2 - \nu_{\bar{u}} \otimes \nu_{\bar{u}})\delta_{S_{\bar{u}}})$ being the algebraic curvature vector of $S_{\bar{u}}$, see [5, 6]. Notice that the curvature form (5.4) is achieved thanks to the fact that $[\bar{u}]$ is constant along the interface $S_{\bar{u}} \cap \Omega$. In our case, $[\bar{u}] = \frac{1}{2}$, $[g(\bar{u})] = \lambda$, and hence $|\kappa_{S_{\bar{u}}}| = 2\lambda$. It is to say that the interfaces of partitions are curves with constant curvature (arcs of radius $\frac{1}{2\lambda}$), see the configuration in Figure 5.

<figure><img src="image_8.png" /><figcaption>FIGURE 5. (a) Configuration of minimal partitions. (b) Equilibrium of junction.</figcaption></figure>

On the other hand, the vector measure $F := \mathrm{div} \mathcal{A}$ defines a force field on $\bar{\Omega}$ which, by the stabilization of shape variations, belongs to the normal cone to the domain $\Omega$. From this point of view, if a junction $J$ of three phases occurs then it must be balanced, therefore the force must vanish at that point, i.e. $F(J) = 0$. Inversely, $F$ is singular at the triple junction and $F(J) = \vec{\tau}_{0,1} + \frac{1}{2}(\vec{\tau}_{0,\frac{1}{2}} + \vec{\tau}_{\frac{1}{2},1}) \neq \vec{0}$, see Figure 5.

## Numerical results.

We present here numerical simulations in $\mathbb{R}^2$ where we solve the convexified problem (5.3) by using the algorithm described in Section 6. They are done over the unit square $\Omega = (0, 1)^2$ and $\lambda = 3/4$. The Dirichlet conditions are varied on the boundary:

$$ u = 0 \text{ on } \Gamma_0^\alpha \quad \text{and} \quad u = 1 \text{ on } \Gamma_1^\alpha \qquad (5.5) $$
where $\Gamma_1^\alpha = \partial\Omega \setminus \Gamma_0^\alpha$ being $\Gamma_1^\alpha = ((0,1] \times \{0\}) \cup (\{1\} \times [0,\alpha))$ for $\alpha \in \{0,0.5,1\}$.

<figure><img src="image_9.png" /><figcaption>FIGURE 6. Problem of three phases or minimal partitions (in (5.3)). Respectively, $\{u=0\}$ is blue, $\{u=1/2\}$ is green, $\{u=1\}$ is red.</figcaption></figure>

## Anisotropic variant.

According to hypothesis (H1), we may substitute the classical total variation $\int |Du|$ with $\int h(Du)$ being $h$ any crystalline norm, that is

$$h(x) = \|x\|_B := \inf\{r \ge 0 : x \in rB\},$$

where $B$ is a convex compact polyhedron containing 0 in its interior. The related anisotropic perimeter $P_B$ in $\mathbb{R}^N$ is defined as

$$P_B(A) := \int_{\partial A} \|\nu_A\|_B d\mathcal{H}^{N-1} \quad \text{for } A \subset \mathbb{R}^N,$$

where $\nu_A$ denotes the outward unit normal to the boundary $\partial A$ of $A$. Let $K$ be the polar of $B$ that is

$$K = \{z^* \in \mathbb{R}^N : z \cdot z^* \le 1, \forall z \in B\}.$$

As the crystalline norm $h(z)$ coincides with the support function of $K$, optimal multiphase configurations will favour interfaces whose normal are orthogonal to the faces of $K$. This is confirmed by the numerical simulations below (see Figures 8 and 9) which have been performed in $\mathbb{R}^2$ in the case where $K$ is a square or a regular hexagon that we rotate with different angles $\theta$, see Figure 7.

<figure><img src="image_10.png" /><figcaption>FIGURE 7. Rotation is applied to convex sets in order to generate variants of orientation on which we define anisotropic seminorms.</figcaption></figure>

The numerical simulations are enforced with the primal-dual algorithms stated in Section 6.
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Here $h$ is the support function of $K_\theta$ deduced from $(-1, 1)^2$ after a rotation of angle $\theta$. From the left to the right the solution for $\theta = 0, \frac{\pi}{6}, \frac{\pi}{4}, \frac{\pi}{3}$, all with $\lambda = 1$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 9. Here $h$ is the support of $K = H_\theta$ where $H_\theta$ is a regular hexagon deduced from $H_0$ by a rotation of angle $\theta$. From the left to the right, the solution for $\theta = 0, \frac{\pi}{6}, \frac{\pi}{4}, \frac{\pi}{3}$, all with $\lambda = 4/5$.</figcaption></figure>

5.2. **A four phases problem.** We exhibit here various numerical simulations of a 4 phases problem which highlight the validity of Theorem 2.4. The non convex term $g$ and its convex envelop are set up as below

$$g(t) = \begin{cases} \lambda & \text{if } t \in [0, a) \\ -\frac{\kappa}{b-a}(t-a) & \text{if } t \in [a, b] \\ \mu & \text{if } t \in (b, 1] \\ +\infty & \text{otherwise} \end{cases}$$

$$g^{**}(t) = \begin{cases} -\frac{\lambda}{a}t + \lambda & \text{if } t \in [0, a) \\ -\frac{\kappa}{b-a}(t-a) & \text{if } t \in [a, b] \\ \frac{\mu+\kappa}{1-b}t - \frac{b\mu+\kappa}{1-b} & \text{if } t \in (b, 1] \\ +\infty & \text{otherwise} \end{cases}$$

where the positive parameters $a, b, \lambda, \mu$ and $\kappa$ satisfy

$$0 < a < b < 1, \quad 0 < \kappa < \frac{\lambda(b-a)}{a}.$$

The second condition ensures that the convex envelop of $g$ shows three different slopes (see Figure 10). If its is not satisfied then $g^{**}$ is affine in the interval $[0, b]$ and, by the exclusion principle, optimal solutions will range in $\{0, b, 1\}$ (avoiding the value $a$ so that $(P_g)$ degenerates to a 3 phases problem).
<figure><img src="image_13.png" /><figcaption>FIGURE 10. (a) Graphs of functions g and g** (in blue). (b) Illustration of the Fenchel conjugate g*. Epigraphs epi g, epi g* are displayed in dashed.</figcaption></figure>

We shall use the algorithms proposed in Section 6 for the simulation. In our numerical experiment, we maintain the prescribed boundary conditions on the partition {$\Gamma_0^{1/2}, \Gamma_1^{1/2}$} of $\partial\Omega$ as the settings (5.5) for the 3-phase problem studied in Section 5.1. This arrangement permits having boundary junctions which are one singular point and one regular point of the domain's boundary $\partial\Omega$. It will serve on further discussions about interface conditions and phase-junctions.

<figure><img src="image_14.png" /><figcaption>FIGURE 11. Problem of 4-phases. Each row corresponds to a choice of certain $(\lambda, \theta)$. In a left-to-right order, the norm $h$ is generated by a ball (Euclidian norm), by a square and by a regular hexagon (crystalline norms).</figcaption></figure>

The numerical simulations shown in Figure 11 represent the solutions of the 4 phases problem for each choice of the norm $h$ and for different values of the $\theta$-rotation (as discussed in Section 5.1) and of the weight factor $\lambda$. The parameters $a, b, \mu, \kappa$ are set up as below:

$$a = 0.25, \quad b = 0.75, \quad \mu = 0.3, \quad \kappa = 0.7,$$

We observe that, in the case of the total variation, the geometric optimality condition (5.4) still holds on $S_{\bar{u}} \cap \Omega$ (which consists of all interfaces appearing in Figure 11). The determination of
the counterpart of condition (5.4) in case of crystalline norms is a delicate issue worth of further investigations.

## 6. MIN-MAX APPROACH. A PRIMAL-DUAL ALGORITHM FOR NON DIFFERENTIABLE LAGRANGIANS

We present several numerical approaches for the (non convex) multiphase problems introduced in the previous sections. There are all based on primal-dual formulations of the convexified problem obtained in Theorem 2.4 which bring us to the search of saddle points for a suitable convex-concave Lagrangian. Recall

(P) is a non convex problem. The convexification principle fundamentally based on Theorem 2.4 is applied to bring us a convex problem under the form

$$ \inf_{u \in C} \int_{\Omega} [h(\nabla u) + g^{**}(u) - f(x)u]dx $$

where $C := \{u \in W^{1,1}(\Omega) : u = u_0$ on $\Gamma$} $ is a closed convex subset in $W^{1,1}(\Omega)$. By the classical duality argument, we can rewrite this infimum problem as a saddle point problem

$$ \inf_{u \in C} \sup_{q \in K} \int_{\Omega} [\nabla u \cdot q + g^{**}(u) - f(x)u]dx \quad (6.1) $$

where $K = \text{dom}(h^*)$. Remark that under the assumption on positive homogeneity of $h$ (see (H1)), one has $h^*(q) = 0$ for all $q \in K$, $h^*(q) = +\infty$ otherwise. In particular, if we replace the density $h$ by a gauge (Minkowski functional) of a closed convex set then the constraint $K$ is indeed the polar of that closed convex set (see Anisotropic variant in Section 5.1). Let us denote by $L(u, q)$ the Lagrangian associated to the inf-sup problem (6.1):

$$ L(u, q) = \int_{\Omega} [\nabla u \cdot q + g^{**}(u) - f(x)u]dx. $$

We can use classical primal-dual algorithms based on Arrow-Hurwicz' method to seek a saddle point of $L(u, q)$ in $C \times K$. But, to that aim, we shall need a regularization of function $g^{**}$ since this function is usually piecewise affine and hence non differentiable. Take for instance the 3 phases problem introduced in Section 5.1, $g^{**}(t)$ is not differentiable at $t = 1/2$. A regularization for $g^{**}$ in that case can be done with

$$ g_{\varepsilon}(t) = \begin{cases} 2\lambda \sqrt{(t - \frac{1}{2})^2 + \varepsilon^2} & \text{if } t \in (0, 1) \\ +\infty & \text{otherwise.} \end{cases} $$

Correspondingly, we denote $L_{\varepsilon}$ the regularized Lagrangian which is obtained by substituting $g^{**}$ with its regularization $g_{\varepsilon}$. Then, using the explicit algorithm described in [16], we produce an iterative process with initialization $\eta, \theta > 0$, $(u_0^h, q_0^h) \in C \times K$, $\bar{u}_0^h = u_0^h$,

$$ q_{n+1}^h = \Pi_K^h \left( q_n^h + \eta \frac{\partial L_{\varepsilon}}{\partial q} \left( \bar{u}_n^h, q_n^h \right) \right) $$

$$ u_{n+1}^h = \Pi_C^h \left( u_n^h - \theta \frac{\partial L_{\varepsilon}}{\partial u} \left( u_n^h, q_{n+1}^h \right) \right) $$

$$ \bar{u}_{n+1}^h = 2u_{n+1}^h - u_n^h $$

where $\Pi_K^h, \Pi_C^h$ are respectively the orthogonal projections on closed convex sets $K, C$. The superscript $h$ indicates that a discretization with mesh size $h$ was done. This algorithm converges to a saddle point (proved in [16]) under the following constraint on the step sizes

$$ 0 < \eta, \quad 0 < \theta < \frac{2}{L_{g_{\varepsilon}'}}, \quad \eta\theta \|\nabla^h\|^2 + \frac{\theta L_{g_{\varepsilon}'}^2}{2} < 1 $$
where $L_{g'_\epsilon}$ denotes the Lipschitz constant of $g'_\epsilon$. In [16], we also proposed a semi-implicit scheme which considerably accelerates the convergence of the iterative algorithm. This implicit scheme is realized with inspiration from the splitting method for Navier-Stokes system (in this circumstance, the couple $(u, q)$ stands for the pressure and the speed of the fluid):

$$
\begin{align*}
q_{n+1}^h &= \Pi_K^h (q_n^h + \eta \nabla^h \bar{u}_n^h) \\
u_{n+1}^h &= u_n^h - \theta(I - \Delta^h)^{-1}(-\operatorname{div}^h q_{n+1} + g'_\epsilon(u_n) - f^h) \\
\bar{u}_{n+1}^h &= 2u_{n+1}^h - u_n^h
\end{align*}
$$

with $\Delta^h$ being a discretized Laplacian and $I$ standing for the identity operator. Note that the projection on $C$ is just to maintain the boundary condition $u = u_0$ on $\Gamma$ and this can be simply implemented within Laplacian operator, thus it is hidden. This scheme gives an iterative process converging under conditions which are more flexible

$$0 < \eta, \quad 0 < \theta < \frac{2}{L_{g'_\epsilon}}, \quad \eta\theta < 1.$$

However, the step size $\theta$ is still constrained by the Lipschitz constant of the derivative of the regularization $g_\epsilon$ which is of order $O(1/\epsilon)$. This is a limited aspect of regularization method since handling the inverse of Laplacian within an iterative algorithm, as we know, is costly. To overcome this challenge, we introduce a new approach using a geometric projection on epigraph. And we shall show that the step sizes are free from Lipschitz constants.

**6.1. Explicit scheme with projection on epi(g**$^{**}$).** Let us rewrite problem (6.1) by introducing a new variable $\alpha \in L^1(\Omega)$:

$$\inf_{g^{**}(u) \le \alpha} \sup_{q \in K} \int_{\Omega} (\nabla u \cdot q - f(x)u + \alpha)dx \quad (6.2)$$

and denote by $\tilde{L}(u, q)$ the corresponding Lagrangian

$$\tilde{L}(u, q) := \int_{\Omega} (\nabla u \cdot q - f(x)u + \alpha)dx.$$

We introduce some notations

$$
\begin{align}
\tilde{u} &:= (u, \alpha), & \tilde{q} &:= (q, -f, \beta), \nonumber \\
\tilde{C} &:= \{(u, \alpha) : u \in C, g^{**}(u) \le \alpha\}, & \tilde{K} &:= \{(q, -f, 1) : q \in K\}, \tag{6.3} \\
A\tilde{u} &:= (\nabla u, u, \alpha), & A^*\tilde{q} &:= (-\operatorname{div} q - f, \beta). \nonumber
\end{align}
$$

Then, $A$ defines a linear operator from $W^{1,1}(\Omega) \times L^1(\Omega)$ to $L^1(\Omega; \mathbb{R}^N) \times L^1(\Omega) \times L^1(\Omega)$, and its adjoint is $A^*$. Problem (6.2) now falls into the form

$$\inf_{\tilde{u} \in \tilde{C}} \sup_{\tilde{q} \in \tilde{K}} \langle \langle A \tilde{u}, \tilde{q} \rangle \rangle \quad (6.4)$$

where $\langle\langle\cdot\rangle\rangle$ is a scalar product defined by

$$\langle\langle A \tilde{u}, \tilde{q} \rangle\rangle := \int_{\Omega} (\nabla u \cdot q - fu + \alpha\beta)dx.$$

We remark that the presence of $\beta$ is just to define the scalar product $\langle\langle\cdot\rangle\rangle$, it will be kept equal to 1 by the convex constraint $\tilde{K}$. And the convex $\tilde{C}$ is indeed (a part of) the epigraph of $g^{**}$ ($u \in C$ is
the boundary condition). Problem (6.4) is a standard problem which has the simplest form. We can easily apply the explicit algorithm based on Arrow-Hurwicz' method described in [16] to obtain

$$
\begin{align*}
\tilde{q}_{n+1}^h &= \Pi_K^h (\tilde{q}_n^h + \eta A \bar{u}_n^h) \\
\tilde{u}_{n+1}^h &= \Pi_C^h (\tilde{u}_n^h - \theta A^* \tilde{q}_{n+1}^h) \\
\bar{u}_{n+1}^h &= 2\tilde{u}_{n+1}^h - \tilde{u}_n^h
\end{align*}
$$

The choice of the step sizes is now of order $O(h)$ and it is no longer squeezed by Lipschitz constants

$0 < \eta, \quad 0 < \theta, \quad \eta\theta\|\nabla^h\|^2 < 1.$

It is evident that transforming problem (6.1) to the standard one (6.4) with the projection on epigraph allows to face the difficulties resulting from the non-differentiability of $g^{**}$. At the moment, a semi-implicit scheme should be generated to help speed up the convergence of the algorithm. The splitting technique may be applied to the Lagrangian $\tilde{L}(u, q)$ to create such a semi-implicit scheme. More precisely, we rewrite the inf-sup problem (6.2) into the form

$$
\inf_{(\phi, u, \alpha) \in \ddot{C}} \sup_{q \in K} \int_{\Omega} (\phi \cdot q - f(x)u + \alpha) dx \quad (6.5)
$$

with $\ddot{C}$ given by

$$
\ddot{C} := \{ (\nabla u, u, \alpha) : u \in C, g^{**}(u) \le \alpha \}.
$$

And the implicitness of the scheme in variable $u$ is effectively enforced thanks to the following projection

$$
\inf_{\substack{v=u_0 \text{ on } \Gamma \\ g^{**}(v) \le a}} \frac{1}{2} \int_{\Omega} \left( |\nabla v - \phi|^2 + |v-u|^2 + |a-\alpha|^2 \right) dx.
$$

This is a slightly heavy projection. But pay attention that if we remove the epigraph component ($g^{**}(v) \le a$) then that projection just becomes a proximal operator of a quadratic form with easily
solvable resolvent which appears as $(I - \Delta^h)^{-1}$ in discretized schemes. For more details in this topic,
we refer to the discussion in [16]. Fortunately, removing the epigraph part in variable $u$ is doable
and even simple by exploiting the conjugate $g^*$ instead of $g^{**}$. In the next subsection, we shall see
that the semi-implicit scheme is still well adapted for the epigraph projection method. This will be
a quite complete setting in order to reach a better acceleration in convergence of our algorithms.

**6.2. Semi-implicit scheme with projection on epi(g").** We start with recasting the saddle
point problem (6.1) as

$$
\inf_{u \in C} \sup_{\tau \in L^\infty(\Omega)} \int_{\Omega} [\nabla u \cdot q + \tau u - g^*(\tau) - f(x)u] dx = \inf_{u \in C} \sup_{q \in K} \int_{\Omega} [\nabla u \cdot q + \tau u - \beta - f(x)u] dx \quad (6.6)
$$

We use here the epigraph of the conjugate $g^*$ (instead of $g^{**}$). Let us set some notations to be
suitable for the current context

$$
\begin{align}
& \hat{u} := (u, \alpha), \quad \hat{q} := (q, -f, \tau, \beta), \nonumber \\
& \hat{C} := \{(u, -1) : u \in C\}, \quad \hat{K} := \{(q, -f, \tau, \beta) : q \in K, g^*(\tau) \le \beta\}, \tag{6.7} \\
& \hat{A}\hat{u} := (\nabla u, u, u, \alpha), \quad (\hat{A})^*\hat{q} = (-\operatorname{div} q - f + \tau, \beta). \nonumber
\end{align}
$$
We observe that the convex set $\hat{C}$ is merely the boundary condition whilst $\hat{K}$ is independently combined by the convex constraint $K$ and the epigraph of $g^*$. With these configuration, problem (6.6) becomes

$$ \inf_{\hat{u} \in \hat{C}} \sup_{\hat{q} \in \hat{K}} \langle \langle \hat{A} \hat{u}, \hat{q} \rangle \rangle \qquad (6.8) $$

where the corresponding scalar product $\langle\langle\cdot, \cdot\rangle\rangle$ is given by

$$ \langle\langle \hat{A} \hat{u}, \hat{q} \rangle\rangle := \int_{\Omega} (\nabla u \cdot q - fu + \tau u + \alpha\beta) dx. $$

The semi-implicit scheme (described in [16]) is available and perfectly fits the standard problem (6.8). It provides a convergent iterative algorithm

$$ \begin{aligned} \hat{q}_{n+1}^h &= \Pi_K^h (\hat{q}_n^h + \eta \hat{A} \bar{u}_n^h) \\ \hat{u}_{n+1}^h &= \Pi_C^h (\hat{u}_n^h - \theta (\hat{A}^* \hat{A})^{-1} (\hat{A}^* (\hat{q}_{n+1}^h))) \\ \bar{u}_{n+1}^h &= 2\hat{u}_{n+1}^h - \hat{u}_n^h \end{aligned} \qquad (6.9) $$

under a rather comfortable choice of step sizes

$$ 0 < \eta, \quad 0 < \theta, \quad \eta\theta < 1. $$

For conveniences in practice, let us unfold the algorithm (6.9) in the expression of the primal variables $u$, $q$ and the extra variables $\tau$, $\beta$:

$$ \begin{aligned} q_{n+1}^h &= \Pi_K^h (q_n^h + \eta \nabla^h \bar{u}_n^h) \\ (\tau_{n+1}^h, \beta_{n+1}^h) &= \Pi_{\text{epi}(g^*)}^h ((\tau_n^h, \beta_n^h) + \eta(\bar{u}_n^h, -1)) \\ u_{n+1}^h &= u_n^h - \theta(2I - \Delta^h)^{-1} \left( -\text{div}^h(q_{n+1}^h) + \tau_{n+1}^h - f^h \right) \\ \bar{u}_{n+1}^h &= 2u_{n+1}^h - u_n^h \end{aligned} \qquad (6.10) $$

To simplify our discussion, we denote in the following

* EReg: explicit scheme combined with regularization,

* EEpi: explicit scheme combined with epigraph projection,

* IEpi: semi-implicit scheme combined with epigraph projection.

As shown in Figures 12, 13 and Table 1, the combination of the semi-implicit scheme and the epigraph projection method in IEpi is more efficient in term of precision of the interfaces and reduces many necessary iterations. We conclude that the algorithm (6.10) is an effective combination between the epigraph projection method and the semi-implicit scheme. This coupling provides a quite powerful numerical approach to treat non-differentiable Lagrangians of saddle point problems of the same type, and contemporaneously brings considerable profit on acceleration of convergence of the algorithm.
<figure><img src="image_15.png" /><figcaption>FIGURE 12. Convergence of primal-dual gap (in $N = 1500$) and comparison in terms of iteration and computational time for the 3-phase problem with the classical total variation and $\lambda = 0.75$.</figcaption></figure>

<table><caption>TABLE 1. Comparison of numerical methods in terms of iteration and computational time for the 3-phase problem with the classical total variation and $\lambda = 0.75$.</caption><caption>primal-dual gap &lt; 10<sup>-3</sup> with MPI in 6 processes</caption><thead><tr><th rowspan="2">N</th><th colspan="3">iteration</th><th colspan="3">time (second)</th></tr><tr><th>EReg</th><th>EEpi</th><th>IEpi</th><th>EReg</th><th>EEpi</th><th>IEpi</th></tr></thead><tbody><tr><td>201</td><td>6329</td><td>3706</td><td>1566</td><td>2.361</td><td>1.599</td><td>49.047</td></tr><tr><td>401</td><td>14164</td><td>7600</td><td>2127</td><td>20.414</td><td>13.269</td><td>110.211</td></tr><tr><td>601</td><td>23393</td><td>11379</td><td>3045</td><td>87.592</td><td>55.751</td><td>316.240</td></tr><tr><td>801</td><td>45089</td><td>14939</td><td>3854</td><td>356.529</td><td>158.577</td><td>760.536</td></tr><tr><td>1001</td><td>63759</td><td>18264</td><td>4790</td><td>897.927</td><td>342.014</td><td>1551.827</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>FIGURE 13. Approximation of interfaces of discontinuity. Here are contour lines of numerical solution $u$ for 3-phase problem ($N = 1500$, $\lambda = 0.75$, primal-dual gap $< 10^{-3}$).</figcaption></figure>

## APPENDIX A. PROOF OF THE DUALITY THEOREM

This section is devoted to the proof of Theorem 3.1. Before proceeding, we need some background about BV functions (Subsection A.1) and then, as a premilary step, we present a recipe allowing to embed the non convex problem ($P_g$) into a convex problem in higher dimension (Subsection A.2).
A.1. **Background around the space BV.** Recall that $BV(\Omega)$ is the set of functions $u$ in $L^1(\Omega)$ whose distributional gradient $Du$ is an element of $\mathcal{M}(\Omega; \mathbb{R}^N)$ the set of Radon vector measure from $\Omega$ to $\mathbb{R}^N$. Equipped with the norm $\|u\| = \int_{\Omega} |u| + |Du|(\Omega)$, $BV(\Omega)$ is a Banach space. If $\Omega$ is bounded, every bounded sequence $(u_n)$ in $BV(\Omega)$ admits a subsequence $(u_{n'})$ converging weakly to some $u \in BV(\Omega)$ in the following sense

$$u_{n'} \to u \text{ in } L^1(\Omega), \quad Du_{n'} \to Du \text{ weakly star in } \mathcal{M}(\Omega; \mathbb{R}^N).$$

By considering only elements $u \in BV(\Omega)$ ranging into a given closed interval $I \subset \mathbb{R}$, we obtain a closed subspace denoted $BV(\Omega; I)$.

For every function $u \in BV(\Omega)$, we denote by $u^{\pm}(x)$ the upper and lower approximate limits of $u$, $S_u := \{x \in \Omega : u^{-}(x) < u^{+}(x)\}$ is the set of all approximate jump points of $u$. $Du$ is a bounded Radon measure, which can be decomposed into

$$Du = \nabla u dx + D^c u + (u^+ - u^-) \nu_u d(\mathcal{H}^{N-1} \llcorner S_u) \quad (A.1)$$

where $\nu_u$ denotes the Radon-Nikodým density of $Du$ with respect to its total variation $|Du|$, i.e. $\nu_u := dDu/d|Du|$. Note that $D^c u$ is the Cantor part of the measure $Du$. The quantity $[u] := u^+ - u^-$ is called the jump of $u$ across the interface $S_u$ and the direction of the jump is given by $\nu_u$ along $S_u$. Accordingly, the complete graph of function $u$, denoted by $\bar{G}_u$, is defined by

$$\bar{G}_u := \bigcup_{x \in \Omega} \left( \{x\} \times [u^{-}(x), u^{+}(x)] \right).$$

It is a $N$-rectifiable subset of $\Omega \times \mathbb{R}$ with an oriented unit normal denoted by $\hat{\nu}_u$. This oriented normal $\hat{\nu}_u$ is $\mathcal{H}^N \llcorner \bar{G}_u$ a.e. determined by

$$\hat{\nu}_u(x, t) = (\nu_u(x), 0) \quad \text{for } x \in S_u \text{ and } t \in [u^{-}(x), u^{+}(x)], \quad (A.2)$$

on the vertical part of $\bar{G}_u$, whereas, on the approximately continuous part $G_u := \{(x, u^{+}(x)) : u^{-}(x) = u^{+}(x)\}$, it is identified as

$$\hat{\nu}_u(x, u^+(x)) = \frac{(\nabla u(x), -1)}{\sqrt{1 + |\nabla u(x)|^2}} \quad (A.3)$$

if $u$ is approximately differentiable at $x$ (with its approximate gradient $\nabla u(x)$), and it is horizontal, i.e. $\hat{\nu}_u(x, u^+(x)) = (\nu_u(x), 0)$, at points in the support of the Cantor part of $Du$. Notice that $\nu_u = dD^c u/d|D^c u|$ $|D^c u|$-a.e in $\Omega$. We remark also that the complete graph $\bar{G}_u$ of functions $u$ belonging to $W^{1,1}(\Omega)$ agrees with the continuous graph $G_u$ (on which $\hat{\nu}_u(x, u(x))$ is given by (A.3)).

**Subgraph functions.** Any element $u$ of $BV(\Omega; I)$ ($I = [m, M]$), can be identified with a function in one more dimension, given by the characteristic function $1_u$ of its subgraph, defined on $Q := \Omega \times (m, M)$ by

$$1_u(x, t) := \begin{cases} 1 & \text{if } t \le u(x) \\ 0 & \text{if } t > u(x). \end{cases} \quad (A.4)$$

It turns out that the distributional gradient of $1_u$ is a bounded measure on $\Omega \times I$ supported on $\bar{G}_u$ namely

$$D1_u = \hat{\nu}_u \mathcal{H}^N \llcorner \bar{G}_u \quad (A.5)$$

Accordingly if $u \in BV(\Omega; I)$ where $I = [m, M]$ a compact interval, then $1_u$ can be identified as an element of $BV(\Omega \times (m, M); [0, 1])$ such that for a.e. $x \in \Omega$, the function $1_u(x, \cdot)$ is monotone non increasing with $1_u(x, m+0) = 0$ and $1_u(x, M-0) = 1$ (if $I = \mathbb{R}$, one has merely that $1_u \in L_{\text{loc}}^1(\Omega \times \mathbb{R}))$).
**Duality pairing with divergence regular vector fields.** Let us assume that $I = [m, M]$ is a bounded interval and set $Q := \Omega \times (m, M)$. Let $v \in BV(Q)$ and let $\sigma \in X_\infty(Q)$ where

$$X_\infty(Q) := \{\sigma := (\sigma^x, \sigma^t) \in L^\infty(Q; \mathbb{R}^{N+1}) : \operatorname{div}\sigma \in L^\infty(Q)\}, \quad (A.6)$$

where the divergence is taken in the distributional sense in the open subset $Q$. Then we obtain a duality pairing between $\sigma$ and $v$, by integrating over $Q$ the scalar measure $\sigma \cdot Dv$ defined by:

$$\langle (\sigma \cdot Dv), \varphi \rangle := - \int_{\Omega \times \mathbb{R}} v (\sigma \cdot \nabla \varphi + \varphi \operatorname{div} \sigma) dx \quad \forall \varphi \in C_0^\infty(Q).$$

It turns out (see [2, Thm 1.5 and Corollary 1.6]) that $(\sigma \cdot Dv)$ is a Radon measure on $Q$ which is absolutely continuous with respect to $|Dv|$ and satisfies

$$\int_Q |(\sigma \cdot Dv)| \le \|\sigma\|_{\infty} \int_Q |Dv|.$$

In addition, as $\Omega$ is Lipschitz, $v$ admits a trace in $L^1(\partial Q)$ and $\sigma$ admits a normal trace in $L^\infty(\partial Q)$. Notice that the boundary $\partial Q$ consists in the union of the lateral part $\partial \Omega \times I$ with exterior normal $(\nu_\Omega, 0)$ and of the lower and upper interfaces $\bar{\Omega} \times \{m, M\}$. The following generalized integration by parts formula holds ([1]):

$$\int_Q (\sigma \cdot Dv) = - \int_Q v \operatorname{div} \sigma \, dx dt + \int_{\partial \Omega \times I} (\sigma^x \cdot \nu_\Omega) v \, d\mathcal{H}^N + \int_\Omega [\sigma^t(x, M)v(x, M) - \sigma^t(x, m)v(x, m)] \, dx. \quad (A.7)$$

Note that the normal trace of $\sigma$ is well defined on every interface $t = a$ for $a \in [m, M]$ so that the vertical component $\sigma^t(\cdot, a)$ is defined a.e. in $\Omega$.

**A.2. Convexified problem in $N+1$ dimensions.** We assume that $\mathrm{dom}(g) \subset I := [m, M]$ where $-\infty < m < M < +\infty$, thus reducing the class of competitors to a subset of $BV(\Omega; I)$. As above we set $Q = \Omega \times (m, M)$. Then following the general procedure studied in [9], we are led to a convexified problem of $(\mathcal{P}_g)$ in dimension $N+1$ of the form

$$(Q) \qquad \inf \{ \hat{F}_g(v) : v \in \mathcal{A} \} \qquad (A.8)$$

where the admissible set is

$$\mathcal{A} := \{ v \in BV(Q; [0, 1]) : v = 1 \text{ on } \Omega \times \{m\}, v = 0 \text{ on } \Omega \times \{M\} \} \quad (A.9)$$

and the convex energy to be minimized is given by

$$\hat{F}_g(v) := J(v) + \ell(v) \quad \text{for } v \in BV(Q; [0, 1]), \quad (A.10)$$

where $J(v)$ and $\ell(v)$ are defined as follows:

$$J(v) := \int_Q \tilde{h}(t, Dv) \quad \text{with} \quad \tilde{h}(t, z^x, z^t) := \begin{cases} h(z^x) - z^t(g(t) - tf) & \text{if } z^t \le 0 \\ +\infty & \text{if } z^t > 0, \end{cases} \quad (A.11)$$

$$\ell(v) := \int_{\Gamma \times I} h((\mathbf{1}_{u_0} - v)\nu_{\Omega}) \, d\mathcal{H}^{N-1} \otimes dt. \quad (A.12)$$

The key property which relates this convex minimization problem $(Q)$ to the original problem $(\mathcal{P})$ is the following:

**Lemma A.1.** For every $u \in BV(\Omega; I)$, it holds:

$$\hat{F}_g(\mathbf{1}_u) = F_g(u) \qquad (A.13)$$
Proof. By (A.5),(A.3), we have:

$$
\begin{align*}
J(\mathbf{1}_u) &= \int_Q \tilde{h}(t, D\mathbf{1}_u) = \int_{\bar{G}_u} \tilde{h}(t, \hat{\nu}_u) \, d\mathcal{H}^N \\
&= \int_{G_u} \tilde{h}(t, \hat{\nu}_u) \, d\mathcal{H}^N + \int_{S_u} \int_{u^-(x)}^{u^+(x)} \tilde{h}(s, (\nu_u, 0)) \, ds \, d\mathcal{H}^{N-1},
\end{align*}
$$

where in the second line we split $\bar{G}_u$ into its approximately continuous part $G_u$ and its jump part.
Recalling that $\tilde{h}$ defined in (A.11) is positively one-homogeneous and that $\hat{\nu}_u(x, u(x))$ is given
respectively by (A.3) for Lebesgue a.e. and by (A.2) for $|Cu|$ a.e. $x \in \Omega$, we have:

$$
\int_{G_u} \tilde{h}(t, \hat{\nu}_u) d\mathcal{H}^N = \int_{\Omega} (h(\nabla u) + g(u) - f u) dx + \int_{\Omega} h\left(\frac{dD^c u}{d|D^c u|}\right) d|D^c u| \quad (A.14)
$$

On the other hand, as for $x \in S_u$, il holds $\hat{\nu}_u(x, t) = (\nu(x), 0)$ for every $t \in [u^-(x), u^+(x)]$ and in view of (A.1), we infer that

$$
\int_{S_u} \int_{u^-(x)}^{u^+(x)} \tilde{h}(s, (\nu_u, 0)) \, ds \, d\mathcal{H}^{N-1} = \int_{S_u} h(\nu_u(x))(u^+ - u^-)(x) \, \mathcal{H}^{N-1}(dx) = \int_{S_u} h(Du) \quad (A.15)
$$

Thus by adding the equalities (A.14) and (A.15), we are led to the identity:

$$
J(\mathbf{1}_u) = \int_{\Omega} h(Du) + \int_{\Omega} (g(u(x)) - u(x)f(x))dx \quad \forall u \in BV(\Omega; I). \quad (A.16)
$$

Thus, in order to show (A.13), we are reduced to check that $\ell(v)$ defined in (A.12) satisfies:

$$
\ell(\mathbf{1}_u) = \int_{\Gamma} h((u_0 - u)\nu_{\Omega})d\mathcal{H}^{N-1} \quad \forall u \in BV(\Omega; I). \quad (A.17)
$$

This follows from the fact that for $x \in \Gamma$, the function $\mathbf{1}_{u_0}(x, \cdot) - \mathbf{1}_u(x, \cdot)$ remains nonnegative on $I$
if $u_0(x) \ge u(x)$ while it remains nonpositive if $u_0(x) \le u(x)$. Therefore, by the homogeneity of $h$
and since $\int_I |\mathbf{1}_{u_0}(x, t) - \mathbf{1}_u(x, t)|dt = |u_0(x) - u(x)|$, we deduce that:

$$
\int_I h((\mathbf{1}_{u_0} - \mathbf{1}_u)(x, t) \nu_\Omega(x)) \, dt = h((u_0(x) - u(x)) \nu_\Omega(x)).
$$

Integrating the last equality with respect to $\mathcal{H}^{N-1} \llcorner \Gamma$ leads to (A.17). $\square$

<figure><img src="image_17.png" /></figure>

From Lemma A.1, it is obvious that $\inf(Q) \le \inf(P)$. We are now going to show that the converse inequality holds by using a slicing argument. Firstly we observe that $\tilde{h}(z^x, z^t)$ given in (A.11) is infinite for $z^t > 0$. Therefore any $v$ with finite energy $J(v) < +\infty$ is such that $v(x, \cdot)$ is monotone non increasing. For such a non increasing $v \in A$ and for all $s \in [0, 1]$, we define

$$
v_s(x) := \inf \{ \tau \in [m, M] : v(x, \tau) \le s \}. \quad (A.18)
$$

Clearly the subgraph function of $v_s$ agrees with the characteristic function of the set $\{(x, \tau) \in \Omega \times I : v(x, \tau) > s\}$. More precisely we have

$$
\mathbf{1}_{v_s}(x, t) = \mathbf{1}_{\{v>s\}}(x, t) \quad \text{for a.e. } (x, t) \in Q. \qquad (A.19)
$$

In view of (A.18), we can therefore associate with every element $v \in A$ such that $J(v) < +\infty$ a parametrized family $\{v_s : 0 \le s \le 1\}$. A key result is the following coarea type formula which allows to recover solutions of the primal problem ($P_g$) from any solution of its convexified problem ($Q$).

**Proposition A.2.** For every $v \in A$ such that $\hat{F}_g(v) < +\infty$, the function $v_s$ defined in (A.18) belongs to $BV(I)$ for a.e. $s \in [0, 1]$ and it holds

$$
\hat{F}_g(v) = \int_0^1 F_g(v_s) \, ds. \qquad (A.20)
$$
As a consequence, we have $\inf(P_g) = \inf(Q)$ and any solution $u$ to $(P)$ provides a solution $\mathbf{1}_u$ to $(Q)$. Conversely if $v$ is a solution to $(Q)$, then $v_s$ is a solution to $(P)$ for a.e. $s \in (0, 1)$.

*Proof.* Firstly, $J(v)$ satisfies the generalized co-area formula, that is

$$J(v) = \int_{-\infty}^{+\infty} J(\mathbf{1}_{\{v>s\}})ds = \int_{-\infty}^{+\infty} J(\mathbf{1}_{v_s})ds \quad \forall v \in L_{\text{loc}}^{1}(\Omega \times \mathbb{R}) \quad (A.21)$$

where the second equality is derived from the identity (A.19). Let $v \in A$ such that $\hat{F}_g(v) < +\infty$. Since $v \in A$ and $v(x, \cdot)$ is non increasing (because of the finiteness of $J(v)$), $v$ only takes values into $[0, 1]$. We recall that the layer cake representation formula for any function $w$ ranging in $[0, 1]$ is given by $w(y) = \int_0^1 \mathbf{1}_{\{w>s\}}(y)ds$. Thus, noticing that $\mathbf{1}_{\{v>s\}} = \mathbf{1}_{v_s}$ (see (A.19)), we have

$$v(x,t) = \int_{0}^{1} \mathbf{1}_{\{v>s\}}(x,t)ds = \int_{0}^{1} \mathbf{1}_{v_s}(x,t)ds.$$

As we knew that $J(v)$ (defined by (A.11)) satisfies the coarea formula and the condition (A.16), we derive that

$$J(v) = \int_{-\infty}^{+\infty} J(\mathbf{1}_{\{v>s\}})ds = \int_{0}^{1} J(\mathbf{1}_{v_s})ds = \int_{0}^{1} \left( \int_{\Omega} h(Dv_s) + \int_{\Omega} (g(v_s) - v_s f(x))dx \right) ds. \quad (A.22)$$

Recall that $\hat{F}_g(v) = J(v) + \ell(v)$ (see (A.10)), proving (A.20) remains to verify that

$$\begin{aligned}
\ell(v) &= \int_{\Gamma \times I} h((\mathbf{1}_{u_0} - v)\nu_{\Omega}) \, d\mathcal{H}^{N-1} \otimes dt \\
&= \int_0^1 \int_{\Gamma \times I} h((\mathbf{1}_{u_0} - \mathbf{1}_{v_s})\nu_{\Omega}) \, d\mathcal{H}^{N-1} \otimes dt \, ds = \int_0^1 \int_{\Gamma} h((u_0 - v_s)\nu_{\Omega}) \, d\mathcal{H}^{N-1} \, ds.
\end{aligned}
\quad (A.23)$$

Taking into account the homogeneity of $h$, we need only to check (it is straightforward) that for every $s \in [0, 1]$ it holds

$$\operatorname{sgn}(\mathbf{1}_{u_0}(x,t) - v(x,t)) = \operatorname{sgn}(\mathbf{1}_{u_0}(x,t) - \mathbf{1}_{v_s}(x,t)) = \operatorname{sgn}(u_0(x) - v_s(x)) \quad \forall t \in I.$$

Therefore, after integrating with respect to $t$, we get:

$$\begin{aligned}
\int_I h((\mathbf{1}_{u_0} - v)\nu_{\Omega})dt &= \int_I \int_0^1 |\mathbf{1}_{u_0} - \mathbf{1}_{v_s}| h(\operatorname{sgn}(\mathbf{1}_{u_0} - v)\nu_{\Omega}) \, ds \, dt \\
&= \int_I \int_0^1 |\mathbf{1}_{u_0} - \mathbf{1}_{v_s}| h(\operatorname{sgn}(\mathbf{1}_{u_0} - \mathbf{1}_{v_s})\nu_{\Omega}) \, ds \, dt \\
&= \int_0^1 |u_0 - v_s| h(\operatorname{sgn}(u_0 - v_s)\nu_{\Omega}) \, ds
\end{aligned}$$

whence the equalities (A.23) by ultimately integrating in $x$ with respect to $\mathcal{H}^{N-1} \llcorner \Gamma$. In summary, we have proved the slicing formula (A.20). It is now easy to show that $\inf(P_g) = \inf(Q)$. Indeed, if $v \in A$, then for $\mathcal{L}^1$-a.e. $s \in (0, 1)$, the function $v_s$ belongs to $BV(\Omega; I)$ and therefore satisfies the inequality $F_g(v_s) \ge \inf(P_g)$. Then, by applying (A.20), we deduce that $\hat{F}_g(v) \ge \inf(P_g)$. Since $v \in A$ was arbitrary, we conclude that $\inf(P_g) \ge \inf(Q)$. The converse inequality is trivial since, for every $u \in BV(\Omega)$, one has $\mathbf{1}_u \in A$ while $F_g(u) = \hat{F}_g(\mathbf{1}_u)$ by (A.13). .

□
### A.3. Proof of Theorem 3.1.
As, by Proposition A.2, we have $inf(Q) = inf(P_g)$, we are reduced to establish the equality $inf(Q) = sup(P_g^*)$. For this purpose, we follow the same strategy as in [7]: we introduce a perturbation of the dual problem $(P_g^*)$, by setting for every $\eta \in C_0(\bar{\Omega} \times I; \mathbb{R}^{N+1})$:

$$ \Phi(\eta) := \inf \left\{ \int_{\Omega} \sigma^t(x, m)dx - \int_{\Gamma \times I} \sigma^x \cdot \nu_{\Omega} \mathbb{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt : \sigma \in X_{\infty}(Q), \right. \\ \left. \text{div } \sigma = 0, \sigma^x \cdot \nu_{\Omega} = 0 \text{ on } (\partial\Omega \setminus \Gamma) \times I, \sigma + \eta \in \mathcal{K}_g \right\} $$

We can easily verify the convexity of the map $\eta \mapsto \Phi(\eta)$, and evidently, for $\eta = 0$, $-\Phi(0)$ is nothing else but $(P_g^*)$, i.e.

$$ \sup(P_g^*) = -\Phi(0). \qquad (A.24) $$

Moreover, the convex function $\Phi$ is continuous at 0 as it is majorized on a small ball of $C_0(\bar{\Omega} \times I; \mathbb{R}^{N+1})$ centered at the origin, namely there exists $\delta > 0$ such that, for any $\eta$ with $\|\eta\|_{\infty} \le \delta$, there is $\sigma_0 \in \mathcal{K}_g$ such that

$$ \Phi(\eta) \le \int_{\Omega} \sigma_0^t(x, m)dx - \int_{\Gamma \times I} \sigma^x \cdot \nu_{\Omega} \mathbb{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt. $$

In fact, we can take $\sigma_0$ of the kind

$$ \sigma_0(x, t) := (0, b(x)) $$

where $b(x) \in L^{\infty}(\Omega)$ is determined later. Let us show that $b(x)$ can be chosen so that $\sigma_0 + \eta \in \mathcal{K}_g$ for every $\eta$ with $\|\eta\|_{\infty} \le \delta$. For $\delta$ being small enough, since $K (= \text{dom}(h^*))$, see (2.3)) is a convex compact containing the origin, one has $\sigma_0^x + \eta^x \in K$. On the other hand, we should chose $b(x)$ so that $\sigma_0 + \eta$ satisfies the conditions (3.7). To this end, the following implication must be ensured

$$ |\eta^t| < \delta \implies \sigma_0^t(x, t) + \eta^t(x, t) + g(t) - tf(x) \ge b(x) - \delta + g(t) - tf(x) \ge 0 $$

Clearly, $b(x)$ can be taken as large as possible so that $b(x) \ge \delta - \inf g + M \sup f_+$.

From the continuity of $\Phi$ at 0 and thanks to a classical duality result due to J.J.Moreau, we deduce that

$$ -\Phi(0) = -\Phi^{**}(0) = \min(\Phi^{*}). \qquad (A.25) $$

where $\Phi^*$ denotes the Moreau-Fenchel conjugate of $\Phi$. Notice that the conjugate is taken in the duality between continuous functions and bounded measures. No we claim that

$$ \Phi^*(\lambda) = \begin{cases} \hat{F}_g(v) & \text{if } \lambda = Dv + (\mathbb{1}_{u_0} - v)(\nu_{\Omega}, 0) \mathcal{H}^N \llcorner (\partial\Omega \times I) \text{ with } v \in \mathcal{A}, \\ +\infty & \text{otherwise.} \end{cases} \qquad (A.26) $$

If (A.26) is true, then we have the equality $\min(\Phi^*) = \min\{\hat{F}_g(v) : v \in \mathcal{A}\}$. Then, by (A.24)-(A.25), we can conclude that $inf(Q) = sup(P^*)$ and the proof of Theorem 3.1 is finished.

Let us sketch the proof of the claim (A.26). This very technical proof can be achieved in the same line azs in [7, Section 7] for which we refer for further details. Let $\lambda$ be a bounded vector measure such that $\Phi^*(\lambda) < +\infty$. The proof of claim (A.26) is supported by the following intermediate statements:

(i) Let $\hat{h} : \mathbb{R} \times \mathbb{R}^N \times \mathbb{R} \to (-\infty, +\infty]$ be defined by

$$ \hat{h}(t, z^x, z^t) := \begin{cases} h(z^x) - z^t g(t) & \text{if } z^t \le 0 \\ +\infty & \text{if } z^t > 0. \end{cases} \qquad (A.27) $$
and $\mathcal{K}_g^0$ denote the subclass of vector fields in $\sigma \in X_\infty(Q)$ (see (A.6)) such that:

$$ \sigma^x(x, t) \in K \quad \text{a.e. in } \Omega \times I \quad (K \text{ defined in (2.3)}) \qquad (A.28) $$

$$ \sigma^t(x, t) + g(t) \ge 0 \quad \forall t \in I \quad \text{and for } \mathcal{L}^N\text{-a.e. } x \in \Omega, \qquad (A.29) $$

Then, for every $\sigma \in X_\infty(Q)$ and every $v \in BV \cap L^\infty(Q)$, it holds

$$ \int_Q \hat{h}(t, Dv) = \sup \left\{ \int_Q (\sigma + \eta) \cdot Dv : \eta \in \mathcal{D}(Q; \mathbb{R}^{N+1}), \sigma + \eta \in \mathcal{K}_g \right\}. \qquad (A.30) $$

(ii) For every compact neighborhood $U$ of the boundary $\Gamma$ and every bounded continuous $\psi: \bar{\Omega} \times I \to \mathbb{R}^{N+1}$, there hold

$$ \langle \lambda - D \mathbf{1}_{u_0}, \psi \rangle = 0 \quad \text{whenever } \operatorname{div} \psi = 0 \text{ in } Q \text{ and } \psi = 0 \text{ on } U \times I, \\ \int_{(\Omega \setminus U) \times I} \tilde{h}(t, \lambda) < +\infty. $$

(iii) There exists a scalar function $v \in L_{\text{loc}}^1(Q)$ with $v(x, \cdot)$ monotone non-increasing such that

$$ \lambda = Dv + (\mathbf{1}_{u_0} - v)(\nu_\Omega, 0)\mathcal{H}^N \llcorner (\partial\Omega \times I). \qquad (A.31) $$

Furthermore, up to adding a constant to $v$, we have $v \in \mathcal{A}$, i.e.

$$ v \in BV(Q; [0, 1]), \quad v(x, m) = 1, \quad v(x, M) = 0. \qquad (A.32) $$

Let $\sigma$ be an element of $X_\infty(Q)$. By the assertion (iii) above, the duality pair $(\sigma \cdot \lambda) = (\sigma \cdot Dv) + \sigma^x \cdot \nu_\Omega(\mathbf{1}_{u_0} - v)\mathcal{H}^N \llcorner (\partial\Omega \times I)$ is well defined. On the other hand, by (i), it holds

$$ \int_Q \tilde{h}(t, Dv) = \sup \left\{ \int_Q (\sigma + \eta) \cdot Dv : \eta \in \mathcal{D}(\bar{Q}; \mathbb{R}^{N+1}), \sigma + \eta \in \mathcal{K}_g \right\}. $$

Let us compute the Moreau-Fenchel conjugate of $\Phi$:

$$ \begin{aligned} \Phi^*(\lambda) &= \sup_{\eta} \left\{ \int_Q \eta \cdot \lambda - \Phi(\eta) : \eta \in C_0(\bar{Q}; \mathbb{R}^{N+1}) \right\} \\ &= \sup_{\eta} \left\{ \int_Q \eta \cdot \lambda - \Phi(\eta) : \eta \in \mathcal{D}(\bar{\Omega} \times I; \mathbb{R}^{N+1}) \right\} \\ &= \sup_{\eta, \sigma} \left\{ \int_Q \eta \cdot \lambda - \int_{\Omega} \sigma^t(x, m) dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_{\Omega} \mathbf{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt : \eta \in \mathcal{D}(\bar{\Omega} \times I; \mathbb{R}^{N+1}), \right. \\ & \qquad \left. \sigma \in X_\infty(Q), \operatorname{div} \sigma = 0, \sigma^x \cdot \nu_{\Omega} = 0 \text{ on } (\partial\Omega \setminus \Gamma) \times I, \sigma + \eta \in \mathcal{K}_g \right\} \end{aligned} \qquad (A.33) $$

where the second line follows from the density of $\mathcal{D}(\bar{Q}; \mathbb{R}^{N+1})$ in $C^0(\bar{Q}; \mathbb{R}^{N+1})$ and the third one from the definition of $\Phi(\eta)$.
Let us evaluate the energy to be maximized in the last line of (A.33) taking an arbitrary competitor $\sigma \in X_\infty(Q)$ such that $\text{div}\,\sigma = 0$ in $Q$ and $\sigma^x \cdot \nu_\Omega = 0$ on $(\partial\Omega \setminus \Gamma) \times I$. We have:

$$
\begin{align*}
& \int_Q \eta \cdot \lambda - \int_\Omega \sigma^t(x, m)dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_\Omega \mathbf{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt \\
&= \int_Q (\eta + \sigma) \cdot \lambda - \int_Q \sigma \cdot \lambda - \int_\Omega \sigma^t(x, m)dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_\Omega \mathbf{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt \\
&= \int_Q (\eta + \sigma) \cdot Dv + \int_{\Gamma \times I} (\eta^x + \sigma^x) \cdot \nu_\Omega(\mathbf{1}_{u_0} - v) d\mathcal{H}^{N-1} \otimes dt \\
&\quad - \int_Q \sigma \cdot Dv - \int_{\Gamma \times I} \sigma^x \cdot \nu_\Omega(\mathbf{1}_{u_0} - v) d\mathcal{H}^{N-1} \otimes dt - \int_\Omega \sigma^t(x, m)dx + \int_{\Gamma \times I} \sigma^x \cdot \nu_\Omega \mathbf{1}_{u_0} d\mathcal{H}^{N-1} \otimes dt \\
&= \int_Q (\eta + \sigma) \cdot Dv + \int_{\Gamma \times I} (\eta^x + \sigma^x) \cdot \nu_\Omega(\mathbf{1}_{u_0} - v) d\mathcal{H}^{N-1} \otimes dt, \tag{A.34}
\end{align*}
$$

where we substituted $\lambda$ with its expression in (A.31) to obtain the second equality while the third equality follows from the fact that, by applying the generalised integration by parts formula (A.7) to $\int_Q \sigma \cdot Dv$ and considering that $v(x, M) = 0$, the expression in the fourth line of (A.34) vanishes. Combining (A.33) and (A.34), we deduce that:

$$
\begin{align*}
\Phi^*(\lambda) &= \sup_{\eta, \sigma} \left\{ \int_Q (\eta + \sigma) \cdot Dv + \int_{\Gamma \times I} (\eta^x + \sigma^x) \cdot \nu_\Omega (\mathbf{1}_{u_0} - v) d\mathcal{H}^{N-1} \otimes dt : \eta \in \mathcal{D}(\bar{\Omega} \times I; \mathbb{R}^{N+1}) \right. \\
&\qquad \left. \sigma \in X_\infty(Q), \text{div}\,\sigma = 0, \sigma^x \cdot \nu_\Omega = 0 \text{ on } (\partial\Omega \setminus \Gamma) \times I, \sigma + \eta \in \mathcal{K}_g \right\} \\
&= \int_Q \tilde{h}(t, Dv) + \int_{\Gamma \times I} h((\mathbf{1}_{u_0} - v)\nu_\Omega) d\mathcal{H}^{N-1} \otimes dt \\
&= \hat{F}_g(v).
\end{align*}
$$

Note that the supremum in $(\eta, \sigma)$ is attained by using the assertion (i) and the fact that $(\eta^x + \sigma^x) \in K = \text{dom}(h^*)$. The claim (A.26) is validated, proving the theorem 3.1. $\square$

REFERENCES

[1] G. Anzellotti, *Traces of bounded vector fields and the divergence theorem*, Preprint, Uni. Trento.

[2] G. Anzellotti, *Pairings between measures and bounded functions and compensated compactness*, Ann. Mat. Pura Appl. (4) **135** (1983), 293–318 (1984).

[3] L. AMBROSIO, N. FUSCO, and D. PALLARA, *Functions of Bounded Variation and Free Discontinuity Problems*, Clarendon Press, Oxford, 2000.

[4] G. BOUCHITTÉ and P. BOUSQUET, *On a degenerate problem in the calculus of variations*, Trans. Amer. Math. Soc. **371** (2019), no. 2, 777–807.

[5] G. BOUCHITTÉ, G. BUTTAZZO, and I. FRAGALÀ, *Mean curvature of a measure and related variational problems*, Ann. Scuola Norm. Sup. Pisa Cl. Sci., XXV (4) (1997), p. 179–196.

[6] G. BOUCHITTÉ and I. FRAGALÀ, *Variational theory of weak geometric structures: the measure method and its applications*, in Variational methods for discontinuous structures, Springer, 2002, pp. 19–40.

[7] G. BOUCHITTÉ and I. FRAGALÀ, *A duality theory for non-convex problems in the calculus of variations*, Archive for Rational Mechanics and Analysis, 229 (2018), pp. 361–415.

[8] G. BOUCHITTÉ, I. FRAGALÀ, and I. LUCARDESI, *Shape derivatives for minima of integral functionals*, Mathematical Programming, 148 (2014), pp. 111–142.

[9] G. BOUCHITTÉ and M. PHAN, *A duality recipe for non-convex variational problems*, Volume spécial des Comptes Rendus Mécanique, 346 (3) (2018), pp. 206–221.

[10] L. CAFFARELLI and K. LEE, *Homogenization of oscillating free boundaries: the elliptic case*, Communications in Partial Differential Equations, 32 (2007), pp. 149–162.
[11] G. DAVID AND T. TORO, *Regularity of almost minimizers with free boundary*, Calculus of Variations and Partial Differential Equations, 54 (2015), pp. 455–524.

[12] R. FINN, *Remarks relevant to minimal surfaces, and to surfaces of prescribed mean curvature*, Journal d'Analyse Mathematique, 14 (1965), pp. 139–160.

[13] E. GIUSTI, *Minimal surfaces and functions of bounded variation*, Monogr. Math., 80 (1984).

[14] H. JENKINS AND J. SERRIN, *The dirichlet problem for the minimal surface equation in higher dimensions*, J. reine angew. Math, 229 (1968), pp. 170–187.

[15] M. MIRANDA, *Un principio di massimo forte per le frontiere minimali e una sua applicazione alla risoluzione del problema al contorno per l'equazione delle superfici di area minima*, Rend. Sem. Mat. Univ. Padova, 45 (1971), pp. 355–366.

[16] M. PHAN AND C. GALUSINSKI, *A semi-implicit scheme based on Arrow-Hurwicz method for saddle point problems*, submitted for publication, (2017). https://arxiv.org/abs/1712.03888.

[17] R. T. ROCKAFELLAR, *Convex analysis*, vol. 28, Princeton Mathematical Series, 1970.

[18] A. VISINTIN, *Nonconvex functionals related to multiphase systems*, SIAM journal on mathematical analysis, 21 (1990), pp. 1281–1304.

[19] T. POCK, A. CHAMBOLLE, D. CREMERS, AND H. BISCHOF, *A convex relaxation approach for computing minimal partitions*, IEEE, (2009).

[20] B. GOLDLUECKE AND D. CREMERS, *Convex relaxation for multilabel problems with product label spaces*, in European conference on computer vision, Springer, 2010, pp. 225–238.

IMATH, UNIVERSITÉ DE TOULON, BP 20132, 83957 LA GARDE CEDEX- FRANCE
*Email address: bouchitte@univ-tln.fr*

DEPARTMENT OF MATHEMATICS, THE UNIVERSITY OF DANANG, UNIVERSITY OF SCIENCE AND EDUCATION, 459 TÔN DÚC THANG, DANANG, VIETNAM
*Email address: ptdminh@ued.udn.vn*