Luca A D Ferrari

Carolin Rossmanith

Benedikt Wirth

Phase field approximations of branched transportation problems

In branched transportation problems mass has to be transported from a given initial distribution to a given final distribution, where the cost of the transport is proportional to the transport distance, but subadditive in the transported mass. As a consequence, mass transport is cheaper the more mass is transported together, which leads to the emergence of hierarchically branching transport networks. We here consider transport costs that are piecewise affine in the transported mass with N affine segments, in which case the resulting network can be interpreted as a street network composed of N different types of streets. In two spatial dimensions we propose a phase field approximation of this street network using N phase fields and a function approximating the mass flux through the network. We prove the corresponding Γ-convergence and show some numerical simulation results.

Introduction

Branched transportation problems constitute a special class of optimal transport problems that have recently attracted lots of interest (see for instance [26, §4.4.2] and the references therein). Given two probability measures µ + and µ - on some domain Ω ⊂ R n , representing an initial and a final mass distribution, respectively, one seeks the most cost-efficient way to transport the mass from the initial to the final distribution. Unlike in classical optimal transport, the cost of a transportation scheme does not only depend on initial and final position of each mass particle, but also takes into account how many particles travel together. In fact, the transportation cost per transport distance is typically not proportional to the transported mass m, but rather a concave, nondecreasing function τ : [0, ∞) → [0, ∞), m → τ (m). Therefore it is cheaper to transport mass in bulk rather than moving each mass particle separately. This automatically results in transportation paths that exhibit an interesting, hierarchically ramified structure (see Figure 1). Well-known models (with parameters α < 1, a, b > 0) are obtained by the choices τ (m) = m α ("branched transport" [START_REF] Gilbert | Minimum cost communication networks[END_REF]), τ (m) = min(m, am + b) ("urban planning" [START_REF] Buttazzo | A model for the optimal planning of an urban area[END_REF][START_REF] Buttazzo | A mass transportation model for the optimal planning of an urban region[END_REF]), τ (m) = 0 if m = 0, 1 else ("Steiner tree problem" [START_REF] Ambrosio | Topics on analysis in metric spaces[END_REF][START_REF] Gilbert | Steiner minimal trees[END_REF][START_REF] Paolini | Existence and regularity results for the Steiner problem[END_REF]). There exist multiple equivalent formulations for branched transportation problems. One particular formulation models the mass flux F as an unknown vector-valued Radon measure that describes the transport from µ + to µ -. The cost functional of the flux, which one seeks to minimize, is defined as

E µ+,µ-[F] = S τ (|σ(x)|) dH 1 (x) + τ (0)|F ⊥ |(R n ) (1)
if div F = µ + -µ -in the distributional sense (so that indeed F transports mass µ + to µ -) and F = σH 1 S + F ⊥ for some countably 1-rectifiable set S ⊂ R n , an H 1 S-measurable function σ : S → R n , and a diffuse part F ⊥ (composed of a Cantor measure and an absolutely continuous measure). Otherwise, E µ+,µ-[F] = ∞.

In this article we devise phase field approximations of the functional E µ+,µ- in two space dimensions in the case of a piecewise affine transportation cost τ (m) = min{α 0 m, α 1 m + β 1 , . . . , α N m + β N } with positive parameters α i , β i . In fact, for a positive phase field parameter ε we consider (a slightly improved but less intuitive version of) the phase field functional

E µ+,µ- ε [σ, ϕ 1 , . . . , ϕ N ] = Ω min α 0 |σ(x)|, min i=1,...,N {ϕ i (x) 2 + α 2 i ε 2 /β i } |σ(x)| 2 2ε + N i=1 β i 2 ε|∇ϕ i (x)| 2 + (ϕ i (x) -1) 2 ε dx if div σ = µ ε + -µ ε -and E µ+,µ- ε
[σ, ϕ 1 , . . . , ϕ N ] = ∞ otherwise. Our main result Theorem 2.1 shows that E µ+,µε Γ-converges in a certain topology to E µ+,µ- as ε 0. Here the vector field σ approximates the mass flux F, the auxiliary scalar phase fields ϕ 1 , . . . , ϕ N disappear in the limit, and µ ε + , µ ε -are smoothed versions of µ + , µ -. One motivation for our particular choice of τ is that this allows a phase field version of the urban planning model, the other motivation is that any concave cost τ can be approximated by a piecewise affine function. Those phase field approximations then allow to find numerical approximations of optimal mass fluxes.

Related work

Phase field approximations represent a widely used tool to approach solutions to optimization problems depending on lower-dimensional sets. The concept takes advantage of the definition of Γ-convergence [START_REF] Braides | Approximation of free-discontinuity problems[END_REF][START_REF] Braides | Γ-convergence for beginners[END_REF][START_REF] Dal | An introduction to Γ-convergence, volume 8 of Progress in Nonlinear Differential Equations and their Applications[END_REF] to approximate singular energies by smoother elliptic functionals such that the associated minimizers converge as well. The term phase field is due to Modica and Mortola [START_REF] Modica | Un esempio di Γ --convergenza[END_REF] who study a rescaled version of the Cahn and Hilliard functional which models the phase transition between two immiscible liquids and which turns out to approximate the perimeter functional of a set. Subsequently a similar idea has been used by Ambrosio and Tortorelli in [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] to obtain an approximation to the Mumford-Shah functional [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associated variational problems[END_REF]. Later on these techniques have been used as well in fracture theory and optimal partitions problems [START_REF] Iurlano | Fracture and plastic models as Γ-limits of damage models under different regimes[END_REF][START_REF] Conti | Phase field approximation of cohesive fracture models[END_REF]. More recently such phase field approximations have also been applied to branched transportation problems and the Steiner minimal tree problem (to find the graph of smallest length connecting a given set of points) which both feature high combinatorial complexity [START_REF] Karp | Reducibility among combinatorial problems[END_REF]. For instance, in [START_REF] Oudet | A Modica-Mortola approximation for branched transport and applications[END_REF] the authors propose an approximation to the branched transport problem based on the Modica-Mortola functional in which the phase field is replaced by a vector-valued function satisfying a divergence constraint. Similarly, in [START_REF] Bonnivard | Approximation of length minimization problems among compact connected sets[END_REF] Santambrogio et al. study a variational approximation to the Steiner minimal tree problem in which the connectedness constraint of the graph is enforced trough the introduction of a geodesic distance depending on the phase field. Our phase field approximations can be viewed as a generalization of recent work by two of the authors [START_REF] Chambolle | A phase-field approximation of the Steiner problem in dimension two[END_REF][START_REF] Chambolle | Variational approximation of size-mass energies for k-dimensional currents[END_REF], in which essentially [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] for τ (m) = αm + β with α, β > 0 is approximated by an Ambrosio-Tortorelli-type functional defined as

Ẽµ+,µ- ε [σ, ϕ] = Ω 1 2ε ϕ(x) 2 |σ(x)| 2 + β 2 ε|∇ϕ(x)| 2 + 1 ε (ϕ(x) -1) 2 dx
if div σ = ρ ε * (µ + -µ -) for a smoothing kernel ρ ε and if ϕ ≥ α √ β ε almost everywhere, and Ẽµ+,µε [σ, ϕ] = ∞ otherwise. The function 1 -ϕ may be regarded as a smooth version of the characteristic function of a 1-rectifiable set (the Steiner tree) whose total length is approximated by the Ambrosio-Tortorelli phase field term in parentheses. In addition, the first term in the integral forces this set to contain the support of a vector field σ which encodes the mass flux from µ + to µ -.

Notation

Throughout, Ω ⊂ R 2 is an open bounded domain with Lipschitz boundary. The spaces of scalar and R n -valued continuous functions on the closure Ω are denoted C(Ω) and C(Ω; R n ) (for m times continuously differentiable functions we use C m ) and are equipped with the supremum norm • ∞ . Their topological duals are the spaces of scalar and R n -valued Radon measures (regular countably additive measures) M(Ω) and M(Ω; R n ), equipped with the total variation norm • M . Weak-* convergence in these spaces will be indicated by * . The subset P(Ω) ⊂ M(Ω) shall be the space of probability measures (nonnegative Radon measures with total variation 1). The total variation measure of some Radon measure λ will be denoted |λ|, and its Radon-Nikodym derivative with respect to |λ| as λ |λ| . The restriction of λ to a measurable set A is abbreviated λ A. The pushforward of a measure λ under a measurable map T is denoted T # λ. The standard Lebesgue and Sobolev spaces are indicated by L p (Ω) and W k,p (Ω), respectively; if they map into R n we write L p (Ω; R n ) and W k,p (Ω; R n ). The associated norms are indicated by • L p and • W k,p . Finally, the n-dimensional Hausdorff measure is denoted H n .

Structure of the paper:

In Section 2 we introduce the approximating energies and precisely state our results. Section 3 is devoted to the proof of the Γ-convergence result. Its first three subsections deal with the Γ -lim inf inequality which is obtained via slicing, while the remaining two subsections prove an equicoercivity result and the Γ-lim sup inequality. Finally, in Section 4 we introduce a numerical discretization and algorithmic scheme to perform some exemplary simulations of branched transportation networks.

Model summary and Γ-convergence result

Here we state in more detail the considered variational model for transportation networks and its phase field approximation as well as the Γ-convergence result.

Introduction of the energies

Before we state the original energy and its phase field approximation, let us briefly recall the objects representing the transportation networks.

Definition 1 (Divergence measure vector field and mass flux).

1. A divergence measure vector field is a measure F ∈ M(Ω; R 2 ), whose weak divergence is a Radon measure, div F ∈ M(Ω), where the weak divergence is defined as

R 2 ψ d div F = - R 2 ∇ψ•dF for all ψ ∈ C 1 (R 2 ) with compact support.
By [START_REF] Šilhavý | Divergence measure vectorfields: their structure and the divergence theorem[END_REF], any divergence measure vector field F can be decomposed as

F = m F θ F H 1 S F + F ⊥ , where S F ⊂ Ω is countably 1-rectifiable, m F : S F → [0, ∞) is H 1 S F - measurable, θ F : S F → S 1 is H 1 S F -
measurable and orients the approximate tangent space of S F , and F ⊥ is H 1 -diffuse, that is, it is singular with respect to the one-dimensional Hausdorff measure on any countably 1-rectifiable set (it is the sum of a Lebesgue-continuous and a Cantor part).

2. A divergence measure vector field F is polyhedral if it is a finite linear combination

F = M i=1 m i θ i H 1 e i
of vector-valued line measures, where M ∈ N and for i = 1, . . . , M we have m i ∈ R, e i ⊂ R 2 a straight line segment, and θ i its unit tangent vector.

3. Given µ + , µ -∈ P(Ω) with compact support in Ω, a mass flux between µ + and µ -is a divergence measure vector field F with div

F = µ + -µ -.
The set of mass fluxes between µ + and µ -is denoted

X µ+,µ-= {F ∈ M(Ω; R 2 ) | div F = µ + -µ -} .
A mass flux between µ + and µ -can be interpreted as the material flow that transports the material from the initial mass distribution µ + to the final mass distribution µ -. Next we specify the cost functional for branched transportation networks for which we will propose a phase field approximation.

Definition 2 (Cost functional). 1. Given N ∈ N and α 0 > α 1 > • • • > α N > 0, 0 = β 0 < β 1 < • • • < β N < ∞, we define the piecewise affine transport cost τ : [0, ∞) → [0, ∞), τ (m) = min i=0,...,N {α i m + β i } . If α 0 = ∞ we interpret τ as τ (m) = 0 if m = 0 , min i=1,...,N {α i m + β i } else.
2. We call µ + , µ -∈ P(Ω) an admissible source and sink if they have compact support in Ω and in the case of α 0 = ∞ can additionally be written as a finite linear combination of Dirac masses.

3. Given admissible µ + , µ -∈ P(Ω), we define the cost functional E µ+,µ-:

X µ+,µ-→ [0, ∞], E µ+,µ-[F] = S F τ (m F (x)) dH 1 (x) + τ (0)|F ⊥ |(Ω)
for α 0 < ∞ (above, τ (0) = α 0 denotes the right derivative in 0) and otherwise

E µ+,µ-[F] = S F τ (m F (x)) dH 1 (x) if F ⊥ = 0, ∞ else.
The phase field functional approximating E µ+,µ-will depend on a vector field σ approximating the mass flux and N Ambrosio-Tortorelli phase fields ϕ 1 , . . . , ϕ N that indicate which term in the definition of τ is active. Definition 3 (Phase field cost functional). Let ε > 0 and ρ : R 2 → [0, ∞) be a fixed smooth convolution kernel with support in the unit disk and R 2 ρ dx = 1.

1. Given admissible source and sink µ + , µ -∈ P(Ω), the smoothed source and sink are

µ ε ± = ρ ε * µ ± , where ρ ε = 1 ε 2 ρ( • ε ) .
2. The set of admissible functions is

X µ+,µ- ε = (σ, ϕ 1 , . . . , ϕ N ) ∈ L 2 (Ω; R 2 ) × W 1,2 (Ω) N div σ = µ ε + -µ ε -, ϕ 1 = . . . = ϕ N = 1 on ∂Ω .
3. The phase field cost functional is given by E

µ+,µ- ε : X µ+,µ- ε → [0, ∞), E µ+,µ- ε [σ, ϕ 1 , . . . , ϕ N ] = Ω ω ε α 0 , γ ε (x) ε , |σ(x)| dx + N i=1 β i L ε [ϕ i ] ,
where we abbreviated (with some p > 1)

L ε [ϕ] = 1 2 Ω ε|∇ϕ(x)| 2 + (ϕ(x) -1) 2 ε dx , γ ε (x) = min i=1,...,N ϕ i (x) 2 + α 2 i ε 2 /β i , ω ε α 0 , γ ε (x) ε , |σ(x)| = γε ε |σ| 2 2 if |σ| ≤ α0 γε/ε α 0 (|σ| -α0 2γε/ε ) if |σ| > α0 γε/ε + ε p |σ(x)| 2 for α 0 < ∞ , ω ε α 0 , γ ε (x) ε , |σ(x)| = γ ε ε |σ| 2 2 for α 0 = ∞ .
Note that the pointwise minimum inside γ ε is well-defined almost everywhere, since all elements of X µ+,µε are Lebesgue-measurable. Note also that for fixed phase fields ϕ 1 , . . . , ϕ N the phase field cost functional E µ+,µε is convex in σ. This ensures the existence of minimizers for E µ+,µε , which follows by a standard application of the direct method.

. , φN ) ∈ X µ+,µ- ε with E µ+,µ- ε [σ, φ1 , . . . , φN ] < ∞. Now consider a minimizing sequence (σ k , ϕ k 1 , . . . , ϕ k N ) ∈ X µ+,µ- ε , k ∈ N, with E µ+,µ- ε [σ k , ϕ k 1 , . . . , ϕ k N ] → inf E µ+,µ- ε monotonically as k → ∞. Since E µ+,µ- ε is coercive with respect to H = L 2 (Ω; R 2 ) × W 1,2 (Ω) N , (σ k , ϕ k 1 , . . . , ϕ k N
) is uniformly bounded in H so that we can extract a weakly converging subsequence, still indexed by k for simplicity,

(σ k , ϕ k 1 , . . . , ϕ k N ) (σ, ϕ 1 , . . . , ϕ N ) .
Due to the closedness of X µ+,µε with respect to weak convergence in H we have

(σ, ϕ 1 , . . . , ϕ N ) ∈ X µ+,µ- ε . Note that the integrand of E µ+,µ- ε is convex in σ(x)
and the ∇ϕ i (x) as well as continuous in σ(x) and the ϕ i (x), thus E µ+,µε is lower semi-continuous along the sequence. Indeed, consider a subsequence along which each term L ε [ϕ k i ] converges and along which the ϕ k i converge pointwise almost everywhere (so that also γ k ε (x) = min i=1,...,N ϕ k i (x) 2 + α 2 i ε 2 /β i converges for almost all x ∈ Ω). By Mazur's lemma, a sequence of convex combinations m k j=k λ k j σ j of the σ k converges strongly (and up to another subsequence again pointwise) so that by Fatou's lemma we have

inf E µ+,µ- ε = lim k→∞ E µ+,µ- ε [σ k , ϕ k 1 , . . . , ϕ k N ] = lim k→∞ Ω ω ε α 0 , γ k ε (x) ε , |σ k (x)| dx + N i=1 β i lim k→∞ L ε [ϕ k i ] ≥ lim k→∞ m k j=k λ k j Ω ω ε α 0 , γ j ε (x) ε , |σ j (x)| dx + N i=1 β i L ε [ϕ i ] ≥ Ω lim inf k→∞ m k j=k λ k j ω ε α 0 , γ j ε (x) ε , |σ j (x)| dx + N i=1 β i L ε [ϕ i ] ≥ Ω lim inf k→∞ m k j=k λ k j ω ε α 0 , inf i=k,...,m k γ i ε (x) ε , |σ j (x)| dx + N i=1 β i L ε [ϕ i ] ≥ Ω lim inf k→∞ ω ε   α 0 , inf i=k,...,m k γ i ε (x) ε , m k j=k λ k j |σ j (x)|   dx + N i=1 β i L ε [ϕ i ] = E µ+,µ- ε [σ, ϕ 1 , . . . , ϕ N ] ,
where we exploited the weak lower semi-continuity of L ε , the monotonicity of ω ε α 0 , γε(x) ε , |σ(x)| in its second argument, its convexity in its last argument, and its continuity in its latter two arguments.

Remark 1 (Regularization of σ). Note that the phase field cost functional E µ+,µε is L 2 (Ω; R 2 )-coercive in σ, which is essential to have sequentially weak compactness of subsets of X µ+,µε with finite cost (and as a consequence existence of minimizers). For α 0 < ∞ this is ensured by the regularization term ε p |σ| 2 (which has no other purpose). Without it, the functional would only feature weak- * coercivity for σ in M(Ω; R 2 ), however, the integral Ω ω ε α 0 , γε(x) ε , |σ(x)| dx with γ ε Lebesgue-measurable would in general not be well-defined for σ ∈ M(Ω; R 2 ).

Remark 2 (Motivation of ω ε α 0 , γε(x) ε , |σ(x)| via relaxation). Keeping the phase fields ϕ 1 , . . . , ϕ N fixed and ignoring the regularizing term ε p |σ| 2 , the in-

tegrand ω ε (α 0 , γε ε , |σ|) is the convexification in σ of min α 0 |σ|, (ϕ 2 1 + α 2 1 ε 2 /β 1 )|σ| 2 2ε , . . . , (ϕ 2 N + α 2 N ε 2 /β N )|σ| 2 2ε = min α 0 |σ|, γ ε ε |σ| 2 2 ,
which shows the intuition of the phase field functional much clearer. Indeed, the minimum over N + 1 terms parallels the minimum in the definition of τ , and the i th term for i = 0, . . . , N describes (part of) the transportation cost α i m + β i . However, since the above is not convex with respect to σ, a functional with this integrand would not be weakly lower semi-continuous in σ and consequently possess no minimizers in general. Taking the lower semi-continuous envelope corresponds to replacing the above by ω ε (α 0 , γε ε , |σ|) (note that this only ensures existence of minimizers, but will not change the Γ-limit of the phase field functional).

Statement of Γ-convergence and equi-coercivity

Let us extend both E µ+,µ-and E

µ+,µ- ε to M(Ω; R 2 ) × L 1 (Ω) N via E µ+,µ-[σ, ϕ 1 , . . . , ϕ N ] = E µ+,µ-[σ] if σ ∈ X µ+,µ-and ϕ 1 = . . . = ϕ N = 1 almost everywhere, ∞ else, E µ+,µ- ε [σ, ϕ 1 , . . . , ϕ N ] = E µ+,µ- ε [σ, ϕ 1 , . . . , ϕ N ] if (σ, ϕ 1 , . . . , ϕ N ) ∈ X µ+,µ- ε , ∞ else.
We have the following Γ-convergence result.

Theorem 2.1 (Convergence of phase field cost functional). For admissible µ + , µ -∈ P(Ω) we have

Γ -lim ε→0 E µ+,µ- ε = E µ+,µ-,
where the Γ-limit is with respect to weak- * convergence in M(Ω; R 2 ) and strong convergence in

L 1 (Ω) N .
The proof of this result is provided in the next section. Together with the following equicoercivity statement, whose proof is also deferred to the next section, we have that minimizers of the phase field cost functional E µ+,µε approximate minimizers of the original cost functional E µ+,µ-. Theorem 2.2 (Equicoercivity). For ε → 0 let (σ ε , ϕ ε 1 , . . . , ϕ ε N ) be a sequence with uniformly bounded phase field cost functional

E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ] < C < ∞.
Then, along a subsequence, σ ε * σ in M(Ω; R 2 ) for some σ ∈ M(Ω; R 2 ) and ϕ ε i → 1 in L 1 (Ω), i = 1, . . . , N . As a consequence, if µ + , µ -∈ P(Ω) are admissible and such that there exists F ∈ X µ+,µ-with E µ+,µ-[F] < ∞, then any sequence of minimizers of E µ+,µε contains a subsequence converging to a minimizer of E µ+,µ-as ε → 0.

Remark 3 (Phase field boundary conditions). Recall that we imposed boundary conditions ϕ i = 1 on ∂Ω. Without those, the recovery sequence from the following section could easily be adapted such that all full phase field profiles near the boundary will be replaced by half, one-sided phase field profiles. It is straightforward to show that the resulting limit functional would become

S F ∩Ω min{α 0 m F , α 1 m F + β 1 , . . . , α N m F + β N } + S F ∩∂Ω min{α 0 m F , α 1 m F + β 1 /2, . . . , α N m F + β N /2} + α 0 |F ⊥ |(Ω) ,
where fluxes along the boundary are cheaper and thus preferred.

Remark 4 (Divergence measure vector fields and flat chains). Any divergence measure vector field can be identified with a flat 1-chain or a locally normal 1-current (see for instance [START_REF] Šilhavý | Divergence measure vectorfields: their structure and the divergence theorem[END_REF]Sec. 5] or [START_REF] Brancolini | General transport problems with branched minimizers as functionals of 1-currents with prescribed boundary[END_REF]Rem. 2.29(3)]; comprehensive references for flat chains and currents are [START_REF] Whitney | Geometric integration theory[END_REF][START_REF] Federer | Geometric measure theory[END_REF]). Furthermore, for a sequence σ j , j ∈ N, of divergence measure vector fields with uniformly bounded div σ j M , weak- * convergence is equivalent to convergence of the corresponding flat 1chains with respect to the flat norm [START_REF] Brancolini | General transport problems with branched minimizers as functionals of 1-currents with prescribed boundary[END_REF]Rem. 2.29(4)]. Analogously, scalar Radon measures of finite total variation and bounded support can be identified with flat 0-chains or locally normal 0-currents [29, Thm. 2.2], and for a bounded sequence of compactly supported scalar measures, weak- * convergence is equivalent to convergence with respect to the flat norm of the corresponding flat 0-chains.

From the above it follows that in Theorems 2.1 and 2.2 we may replace weak- * convergence by convergence with respect to the flat norm. Indeed, for both results it suffices to consider sequences (σ ε , ϕ ε 1 , . . . , ϕ ε N ) with uniformly bounded cost E µ+,µε

. For those we have uniformly bounded σ ε M (by Theorem 2.2) as well as uniformly bounded div σ ε M = µ ε + -µ ε -M so that weak- * and flat norm convergence are equivalent.

The Γ-limit of the phase field functional

In this section we prove the Γ-convergence result. As is canonical, we begin with the lim inf-inequality, after which we prove the lim sup-inequality as well as equicoercivity.

The Γ -lim inf inequality for the dimension-reduced problem

Here we consider the energy reduced to codimension-1 slices of the domain Ω. In our particular case of a two-dimensional domain, each slice is just onedimensional, which will simplify notation a little (the procedure would be the same for higher codimensions, though). The reduced functional depends on the (scalar) normal flux ϑ through the slice as well as the scalar phase fields ϕ 1 , . . . , ϕ N restricted to the slice.

Definition 4 (Reduced functionals). Let I ⊂ R be an open interval.

1. The decomposition of a measure ϑ ∈ M(I) into its atoms and the remain-

der is denoted ϑ = m ϑ H 0 S ϑ + ϑ ⊥ ,
where S ϑ ⊂ I is the set of atoms of ϑ, m ϑ : S ϑ → R is H 0 S ϑ -measurable, and ϑ ⊥ contains no atoms.

2. We define the reduced cost functional G[•; I] :

M(I) → [0, ∞), G[ϑ; I] = x∈S ϑ ∩I τ (|m ϑ (x)|) + τ (0)|ϑ ⊥ |(I)
for α 0 < ∞ (above, τ (0) = α 0 denotes the right derivative in 0) and otherwise

G[ϑ; I] = x∈S ϑ ∩I τ (|m ϑ (x)|) if ϑ ⊥ = 0, ∞ else. Its extension to M(I) × L 1 (I) N is G[•; I] : M(I) × L 1 (I) N → [0, ∞), G[ϑ, ϕ 1 , . . . , ϕ N ; I] = G[ϑ; I] if ϕ 1 = . . . = ϕ N = 1 almost everywhere, ∞ else.
3. For any (ϑ, ϕ 1 , . . . , ϕ N ) ∈ L 2 (I) × W 1,2 (I) N we define the reduced phase field functional on I as

G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] = I ω ε α 0 , γε(x) ε , |ϑ(x)| dx + N i=1 β i L ε [ϕ i ; I] , L ε [ϕ; I] = 1 2 I ε|ϕ (x)| 2 + (ϕ(x) -1) 2 ε dx , with ω ε and γ ε from Definition 3. Likewise we define G ε [•; I] : M(I) × L 1 (I) N → [0, ∞), G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] = G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] if (ϑ, ϕ 1 , . . . , ϕ N ) ∈ L 2 (I) × W 1,2 (I) N , ∞ else.
For notational convenience, we next introduce the sets

K ε i on which the pointwise minimum inside G ε (or also E µ+,µ- ε ) is realized by the i th element. Definition 5 (Cost domains). For given (ϑ, ϕ 1 , . . . , ϕ N ) ∈ L 2 (I) × W 1,2 (I) N we set K ε 0 = K ε 0 (ϑ, ϕ 1 , . . . , ϕ N ; I) = x ∈ I |ϑ(x)| > α 0 ε γ ε , K ε i = K ε i (ϑ, ϕ 1 , . . . , ϕ N ; I) = x ∈ I \ i-1 j=0 K ε j ϕ i (x) 2 + α 2 i ε 2 βi = γ ε (x) , i = 1, . . . , N .
The sets are analogously defined for (σ,

ϕ 1 , . . . , ϕ N ) ∈ L 2 (Ω; R 2 ) × W 1,2 (Ω) N ,
where we use the same notation (which case is referred to will be clear from the context).

We now show the following lower bound on the energy, from which the Γ -lim inf inequality for the dimension-reduced situation will automatically follow.

Proposition 3.1 (Lower bound on reduced phase field functional). Let

I = (a, b) ⊂ R and 0 ≤ δ ≤ η ≤ 1. Let I η ⊂ {x ∈ I | ϕ 1 (x), . . . , ϕ N (x) ≥ η},
and denote the collection of connected components of I \I η by C η . Furthermore define the subcollection of connected components

C ≥ η = {C ∈ C η | inf C ϕ 1 , . . . , inf C ϕ N ≥ δ} and C ≥ = C∈C ≥ η C. Finally assume ϕ i (a), ϕ i (b) ≥ η for all i = 1, . . . , N . 1. If α 0 < ∞ we have G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] ≥ (η -δ) 2 Iη∪C ≥ α 0 |ϑ| dx + (η -δ) 2 C∈Cη\C ≥ η max β 1 , τ C |ϑ| dx -α 2 0 H 1 (I) ε δ 2 . 2. If α 0 = ∞ we have G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] ≥ δ 2 2εH 1 (I) Iη∪C ≥ |ϑ| dx 2 + (η -δ) 2 C∈Cη\C ≥ η max β 1 , τ C |ϑ| dx -α 2 1 H 1 (I) ε δ 2 .
Proof. 1. (α 0 < ∞) We first show that without loss of generality we may assume

|ϑ| ≥ α 0 ε 2γ ε 1 1 -(η -δ) 2 on K ε 0 . ( 2
)
The motivation is that there may be regions in which a phase field ϕ i is (still) small, but in which we actually have to pay α 0 |ϑ|. Thus, in those regions we would like ω ε (α 0 , γε(x) ε , |ϑ(x)|) to approximate α 0 |ϑ(x)| sufficiently well, and the above condition on ϑ ensures

ω ε α 0 , γ ε (x) ε , |ϑ(x)| = α 0 |ϑ(x)| - α 2 0 ε 2γ ε (x) ≥ (η -δ) 2 α 0 |ϑ(x)| for x ∈ K ε 0 .
(3) We achieve (2) by modifying ϑ while keeping the cost as well as Iη |ϑ| dx and C |ϑ| dx for all C ∈ C η the same so that the overall estimate of the proposition is not affected. The modification mimics the relaxation from Remark 2: the modified ϑ oscillates between small and very large values. Indeed, for fixed

C ∈ C η ∪ {I η } and x ∈ C set θ(x) =      max{ α0ε 2γε(x) 1 1-(η-δ) 2 , ϑ(x)} if x ∈ K ε 0 ∩ (-∞, t C ] , α0ε γε(x) if x ∈ K ε 0 ∩ (t C , ∞) , ϑ(x) else, where t C is chosen such that C | θ| dx = C |ϑ| dx (this is possible, since for t C = ∞
G ε [ θ, ϕ 1 , . . . , ϕ N ; I] -G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] = K ε 0 ω ε α 0 , γ ε (x) ε , | θ(x)| -ω ε α 0 , γ ε (x) ε , |ϑ(x)| dx = K ε 0 α 0 | θ(x)| - α 2 0 ε 2γ ε (x) -α 0 |ϑ(x)| - α 2 0 ε 2γ ε (x) dx = α 0 K ε 0 | θ(x)| dx -α 0 K ε 0 |ϑ(x)| dx = 0 .
Note that the modification θ has a different set K ε 0 ( θ, ϕ 1 , . . . , ϕ N ; I) than the original ϑ. Indeed, by definition of θ we have

K ε 0 ( θ, ϕ 1 , . . . , ϕ N ; I) ⊂ K ε 0 (ϑ, ϕ 1 , . . . , ϕ N ; I) and |ϑ(x)| ≥ α0ε 2γε(x) 1 1-(η-δ) 2 on K ε 0 ( θ, ϕ 1 , . . . , ϕ N ; I), as desired. Let us now abbreviate m 0 = (Iη∪C ≥ )\K ε 0 |ϑ| dx. Using the definition of ω ε as well as γ ε ≥ δ 2 on I η ∪ C ≥ we compute G ε [ϑ, ϕ 1 , . . . , ϕ N ; I η ∪C ≥ ] ≥ (Iη∪C ≥ )∩K ε 0 ω ε α 0 , γ ε ε , |ϑ| dx+ δ 2 2ε (Iη∪C ≥ )\K ε 0 |ϑ| 2 dx ≥ (η -δ) 2 Iη∪C ≥ α 0 |ϑ| dx -α 0 m 0 + δ 2 2εH 1 ((I η ∪ C ≥ ) \ K ε 0 ) m 2 0 ,
where we have employed (3) and Jensen's inequality. Upon minimizing in m 0 , which yields the optimal value α 0 εH 1 ((

I η ∪ C ≥ ) \ K ε 0 )/δ 2 for m 0 , we thus obtain G ε [ϑ, ϕ 1 , . . . , ϕ N ; I η ∪ C ≥ ] ≥ (η -δ) 2 Iη∪C ≥ α 0 |ϑ| dx - α 2 0 εH 1 (I η ∪ C ≥ ) 2δ 2 .
Next consider for each C ∈ C η \ C ≥ η and i = 1, . . . , N the subsets

C ≥ i = C ∩ K ε i ∩ {ϕ i ≥ δ} , C < i = C ∩ K ε i ∩ {ϕ i < δ} ,
and abbreviate m A = A |ϑ| dx for any A ⊂ I. Using Young's inequality, for i, j = 1, . . . , N we have

C ≥ i ϕ 2 i + α 2 i ε 2 /β i 2ε |ϑ| 2 dx ≥ C ≥ i δ 2 + α 2 i ε 2 /β i 2ε |ϑ| 2 dx ≥ C ≥ i α j |ϑ| - α 2 j ε/2 δ 2 + α 2 i ε 2 /β i dx ≥ α j m C ≥ i - α 2 j ε 2δ 2 H 1 (C ≥ i ) .
Similarly, using Jensen's inequality we have

C < i ϕ 2 i + α 2 i ε 2 /β i 2ε |ϑ| 2 + β i 2ε (ϕ i -1) 2 dx ≥ C < i α 2 i ε 2β i |ϑ| 2 + β i 2ε (δ -1) 2 dx ≥ α 2 i ε 2β i 1 H 1 (C < i ) C < i |ϑ| dx 2 + β i 2ε (1 -δ) 2 H 1 (C < i ) ≥ α i (1 -δ)m C < i ,
where in the last step we optimized for H 1 (C < i ). Finally, if inf C ϕ i ≤ δ we have (using Young's inequality)

C\C < i β i 2 ε|ϕ i | 2 + (ϕ i -1) 2 ε dx ≥ β i C\C < i |ϕ i | |1 -ϕ i | dx ≥ β i d c |ϕ i | |1 -ϕ i | dx + f e |ϕ i | |1 -ϕ i | dx ≥ β i δ η ϕ i -1 dϕ i + η δ 1 -ϕ i dϕ i = β i ((1 -δ) 2 -(1 -η) 2 ) ≥ β i (η -δ) 2 ,
where (c, d) and (e, f ) denote the first and the last connected component of

C \ C < i . Next, for C ∈ C η \ C ≥ η define j(C) = max{j ∈ {1, . . . , N } | inf C ϕ j < δ}.
Summarizing the previous estimates we obtain

G ε [ϑ, ϕ 1 , . . . , ϕ N ; I \ (I η ∪ C ≥ )] ≥ C∈Cη\C ≥ η C∩K ε 0 ω ε α 0 , γ ε ε , |ϑ| dx + N i=1 C ≥ i ϕ 2 i + α 2 i ε 2 /β i 2ε |ϑ| 2 dx + C < i ϕ 2 i + α 2 i ε 2 /β i 2ε |ϑ| 2 + β i 2ε (ϕ i -1) 2 dx + C\C < i β i 2 ε|ϕ i | 2 + (ϕ i -1) 2 ε dx ≥ C∈Cη\C ≥ η   (η -δ) 2 α 0 m C∩K ε 0 + N i=1 α j(C) m C ≥ i - α 2 j(C) ε 2δ 2 H 1 (C ≥ i ) + α i (1 -δ)m C < i + N i=1 inf C ϕi≤δ β i (η -δ) 2    ≥ C∈Cη\C ≥ η (η -δ) 2 (α j(C) m C + β j(C) ) - α 2 j(C) ε 2δ 2 H 1 (C) ≥ (η -δ) 2 C∈Cη\C ≥ η max β 1 , τ C |ϑ| dx - α 2 1 ε 2δ 2 H 1 (I) .
Finally, we obtain the desired estimate,

G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] = G ε [ϑ, ϕ 1 , . . . , ϕ N ; I η ∪C ≥ ]+G ε [ϑ, ϕ 1 , . . . , ϕ N ; I\(I η ∪C ≥ )] ≥ (η-δ) 2 Iη∪C ≥ α 0 |ϑ| dx- α 2 0 εH 1 (I) 2δ 2 +(η-δ) 2 C∈Cη\C ≥ η max β 1 , τ C |ϑ| dx - α 2 1 ε 2δ 2 H 1 (I) .
2. (α 0 = ∞) In this case the set K ε 0 is empty, and the cost functional reduces to

G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] = I γ ε (x)|ϑ(x)| 2 2ε + N i=1 β i 2 ε|ϕ i (x)| 2 + (ϕ i (x) -1) 2 ε dx.
With Jensen's inequality we thus compute

G ε [ϑ, ϕ 1 , . . . , ϕ N ; I η ∪ C ≥ ] ≥ Iη∪C ≥ γ ε 2ε |ϑ| 2 dx ≥ δ 2 2ε Iη∪C ≥ |ϑ| 2 dx ≥ δ 2 2ε Iη∪C ≥ |ϑ| dx 2 H 1 (I η ∪ C ≥ ) ≥ δ 2 2ε Iη∪C ≥ |ϑ| dx 2 H 1 (I) .
Furthermore, the same calculation as in the case α 0 < ∞ yields

G ε [ θ, ϕ 1 , . . . , ϕ N ; I \ (I η ∪ C ≥ )] ≥ (η -δ) 2 C∈Cη\C ≥ η max β 1 , τ C |ϑ| dx - α 2 1 ε 2δ 2 H 1 (I)
so that we obtain the desired estimate

G ε [ϑ, ϕ 1 , . . . , ϕ N ; I] = G ε [ϑ, ϕ 1 , . . . , ϕ N ; I η ∪C ≥ ]+G ε [ϑ, ϕ 1 , . . . , ϕ N ; I\(I η ∪C ≥ )] ≥ δ 2 2ε Iη∪C ≥ |ϑ| dx 2 H 1 (I) +(η-δ) 2 C∈Cη\C ≥ η max β 1 , τ C |ϑ| dx - α 2 1 ε δ 2 H 1 (I) .
Corollary 3.1 (Γ -lim inf inequality for reduced functionals). Let J ⊂ R be open and bounded, ϑ ∈ M(J), and ϕ 1 , . . . , ϕ N ∈ L 1 (Ω). Then

Γ -lim inf ε→0 G ε [ϑ, ϕ 1 , . . . , ϕ N ; J] ≥ G[ϑ, ϕ 1 , . . . , ϕ N ; J]
with respect to weak- * convergence in M(J) and strong convergence in L 1 (J) N .

Proof. Let (ϑ ε , ϕ ε 1 , . . . , ϕ ε N ) be an arbitrary sequence converging to (ϑ, ϕ 1 , . . . , ϕ N ) in the considered topology as ε → 0, and assume without loss of generality that the limit inferior of G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; J] is actually a limit and is finite (else there is nothing to show). Further we may assume

G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; J] = G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; J] along the sequence. It suffices to show the lim inf-inequality for each connected component Ĩ = (a, b) of J separately. Due to G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; Ĩ] ≥ βi 2ε ϕ ε i -1 2 L 2 for i = 1, . . . , N we must have ϕ ε i → 1 in L 2 ( Ĩ)
and thus also in L 1 ( Ĩ) so that ϕ i = 1. Even more, after passing to another subsequence, by Egorov's theorem all ϕ ε i converge uniformly to 1 outside a set of arbitrarily small measure. In particular, for any ξ > 0 we can find an open interval (a + ξ, b -ξ) ⊂ I ⊂ Ĩ such that ϕ ε i → 1 uniformly on ∂I, and for any η < 1 there is an open set

I η ⊂ I with H 1 (I \ I η ) ≤ 1 -η such that ϕ ε i ≥ η on I η ∪ ∂I for all i = 1, .
. . , N and all ε small enough.

We now choose δ = ε 1/3 and η = 1 -ε and denote by C η (ε) and C ≥ η (ε) the collections of connected components of I \ I η from Proposition 3.1 (which now depend on ε). Further we abbreviate

C < η (ε) = C η (ε) \ C ≥ η (ε)
. The bound of Proposition 3.1 implies that the number of elements in C < η (ε) is bounded uniformly in ε and η. Passing to another subsequence we may assume C < η (ε) to contain K sets C 1 (ε), . . . , C K (ε) whose midpoints converge to x 1 , . . . , x K ∈ I, respectively. Thus for an arbitrary ζ > 0 we have that for all ε small enough each

C ∈ C < η (ε) lies inside the closed ζ-neighbourhood B ζ ({x 1 , . . . , x K }) of {x 1 , . . . , x K }. Now for α 0 < ∞ we obtain from Proposition 3.1 lim inf ε→0 G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; I] ≥ lim inf ε→0 (η -δ) 2 Iη∪C ≥ (ε) α 0 |ϑ ε | dx + (η -δ) 2 K i=1 max β 1 , τ Ci(ε) |ϑ ε | dx ≥ lim inf ε→0 (η -δ) 2 I\B ζ ({x1,...,x K }) α 0 |ϑ ε | dx + (η -δ) 2 K i=1 τ B ζ ({xi}) |ϑ ε | dx ≥ α 0 |ϑ|(I \ B ζ ({x 1 , . . . , x K })) + K i=1 τ (|ϑ|(B ζ ({x i }))) ,
where in the second inequality we used

G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; I] ≥ α 0 |ϑ|(I \ {x 1 , . . . , x K }) + K i=1 τ (|m ϑ (x i )|) ≥ α 0 |ϑ ⊥ |(I) + S ϑ ∩I τ (|m ϑ |) dH 0 = G[ϑ; I] .
If on the other hand α 0 = ∞ we obtain from Proposition 3.1

lim inf ε→0 G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; I] ≥ lim inf ε→0 δ 2 2εH 1 (I) Iη∪C ≥ (ε) |ϑ ε | dx 2 ≥ lim inf ε→0 δ 2 2εH 1 (I) I\B ζ ({x1,...,x K }) |ϑ ε | dx 2 which implies |ϑ|(I \ {x 1 , . . . , x K }) = 0 and thus |ϑ ⊥ |(I) = 0 as well as S ϑ ∩ I ⊂ {x 1 , . . . , x K }. Next note that for all i ∈ {1, . . . , K} with |m ϑ (x i )| > 0 we have Ci(ε) |ϑ ε | dx > 0 for all ε small enough. Indeed, |m ϑ (x i )| ≤ lim inf ε→0 B ζ ({xi}) |ϑ ε | dx = lim inf ε→0 Ci(ε) |ϑ ε | dx+ B ζ ({xi})\Ci(ε) |ϑ ε | dx ,
lim inf ε→0 G ε [ϑ ε , ϕ ε 1 , . . . , ϕ ε N ; I] ≥ lim inf ε→0 δ 2 2εH 1 (I) Iη∪C ≥ (ε) |ϑ ε | dx 2 + (η -δ) 2 K i=1 |m ϑ (xi)|>0 τ Ci(ε) |ϑ ε | dx ≥ lim inf ε→0 δ 2 2εH 1 (I) K i=1 |m ϑ (xi)|>0 B ζ ({xi})\Ci(ε) |ϑ ε | dx 2 + K i=1 |m ϑ (xi)|>0 τ B ζ ({xi}) |ϑ ε | dx -α 1 B ζ ({xi})\Ci(ε) |ϑ ε | dx ≥ lim inf ε→0 K i=1 |m ϑ (xi)|>0 τ B ζ ({xi}) |ϑ ε | dx - α 2 1 εH 1 (I) 2δ 2 ≥ K i=1 |m ϑ (xi)|>0 τ (|ϑ|(B ζ ({x i }))) = S ϑ τ (|m ϑ |) dH 0 = G[ϑ; I] ,
where in the second inequality we used τ (m 1 + m 2 ) ≤ τ (m 1 ) + α 1 m 2 for any m 1 > 0, m 2 ≥ 0 and in the third we optimized in

B ζ ({xi})\Ci(ε) |ϑ ε | dx.
The proof is concluded by letting ξ → 0 and noting lim inf ξ→0 G[ϑ; I] ≥ G[ϑ; Ĩ].

Slicing of vector-valued measures

We derive now some technical construction for divergence measure vector fields which are needed to reduce the Γ -lim inf inequality to the lower-dimensional setting of the previous section. In particular, we will introduce slices of a divergence measure vector field, which in the language of geometric measure theory correspond to slices of currents. We will slice in the direction of a unitary vector ξ ∈ S n-1 with orthogonal hyperplanes of the form

H ξ,t = π -1 ξ (t) for the projection π ξ : R n → R, π ξ (x) = x • ξ .
The orthogonal projection onto H ξ,t is denoted

π H ξ,t (x) = (I -ξ ⊗ ξ)x + tξ .
The slicing will essentially be performed via disintegration. Let σ be a compactly supported divergence measure vector field. By the Disintegration Theorem [1, Thm. 2.28], for all ξ ∈ S n-1 and almost all t ∈ R there exists a unique measure

ν ξ,t ∈ M(H ξ,t ) such that ν ξ,t M = 1 and σ • ξ = ν ξ,t ⊗ π ξ # |σ • ξ|(t) .
We decompose π ξ # |σ • ξ| into its absolutely continuous and singular part according to

π ξ # |σ • ξ| = σ ξ (t) dt + σ ⊥ ξ
for dt the Lebesgue measure on R.

Lemma 3.1. For any ξ ∈ S n-1 and any compactly supported divergence measure vector field σ we have σ ⊥ ξ = 0, that is, the measure π ξ # |σ • ξ| = σ ξ (t) dt is absolutely continuous with respect to the Lebesgue measure on R. Moreover, for almost all t ∈ R and any compactly supported θ ∈ C ∞ (R n ) we have

σ ξ (t) H ξ,t θ dν ξ,t = {ξ•x<t} ∇θ • dσ + {ξ•x<t} θ d div σ . (4)
Proof. Abbreviate H = ξ ⊥ = H ξ,0 with corresponding orthogonal projection π H , let φ ∈ C ∞ (H) and ψ ∈ C ∞ (R) be compactly supported, and define

I(φ, ψ) = R n φ(π H (x))ψ(π ξ (x)) d(σ • ξ)(x) .
Introducing Ψ(t) = t -∞ ψ(s) ds we obtain via the chain and product rule

(φ • π H )(ψ • π ξ )ξ = (φ • π H )∇[Ψ • π ξ ] = ∇[(φ • π H )(Ψ • π ξ )] -∇[φ • π H ](Ψ • π ξ )
so that (denoting by χ A the characteristic function of a set A)

I(φ, ψ) = R n ∇[(φ • π H )(Ψ • π ξ )] • dσ - R n ∇[φ • π H ](Ψ • π ξ ) • dσ = - R n (φ • π H )(Ψ • π ξ ) d div σ - R n ∇[φ • π H ](Ψ • π ξ ) • dσ = - R n φ(π H (x)) π ξ (x) -∞ ψ(s) ds d div σ(x) - R n ∇[φ • π H ](x) π ξ (x) -∞ ψ(s) ds • dσ = - R n φ(π H (x)) R χ {ξ•x≥s} ψ(s) ds d div σ(x) - R n ∇[φ • π H ](x) R χ {ξ•x≥s} ψ(s) ds • dσ .
(Note that we could just as well have used χ {ξ•x>s} instead of χ {ξ•x≥s} , which would ultimately lead to integration domains {ξ • x ≤ t} in (4); for almost all t this will be the same.) Applying the Fubini-Tonelli Theorem we obtain

I(φ, ψ) = - R ψ(s) {ξ•x≥s} ∇[φ • π H ](x) • dσ(x) + {ξ•x≥s} φ(π H (x)) d div σ(x) ds = R ψ(s) {ξ•x<s} ∇[φ • π H ](x) • dσ(x) + {ξ•x<s} φ(π H (x)) d div σ(x) ds ,
where in the second step we just added 0

= R n ∇[φ•π H ]•dσ + R n φ•π H d div σ
in the square brackets. On the other hand, using the disintegration of σ • ξ we also have

I(φ, ψ) = R n φ(π H (x))ψ(π ξ (x)) d(σ • ξ)(x) = R ψ(s) H ξ,s φ(π H (y)) dν ξ,s (y) (σ ξ (s) ds + dσ ⊥ ξ (s)) .
Comparing both expressions for I(φ, ψ) we can identify

H ξ,s φ(π H (y)) dν ξ,s (y) (σ ξ (s) ds + dσ ⊥ ξ (s)) = {ξ•x<s} ∇[φ • π H ](x) • dσ(x) + {ξ•x<s} φ(π H (x)) d div σ(x) ds .
Since the right-hand side has no singular component with respect to the Lebesgue measure, we deduce H ξ,s φ(π H (y)) dν ξ,s (y) σ ⊥ ξ (s) = 0. Now note that any compactly supported function in C 0 (R n ) or C 1 (R n ) can be arbitrarily well approximated (in the respective norm) by finite linear combinations of tensor products (φ

• π H )(ψ • π ξ ) with φ ∈ C ∞ (H) and ψ ∈ C ∞ (R) with compact
support. Thus, the above implies R H ξ,s θ(y) dν ξ,s (y) dσ ⊥ ξ (s) = 0 for any compactly supported θ ∈ C 0 (R n ) so that ν ξ,s ⊗ σ ⊥ ξ (s) = 0 and thus σ ⊥ ξ (s) = 0 .

Summarizing, we have σ • ξ = σ ξ (s)ν ξ,s ⊗ ds and

H ξ,s φ(π H (x))σ ξ (s) dν ξ,s (x) ds = {ξ•x<s} ∇[φ • π H ](x) • dσ(x) + {ξ•x<s} φ(π H (x))d div σ(x)
for all compactly supported φ ∈ C ∞ (H). Note that the right-hand side is leftcontinuous in s so that the left-hand side is as well. Consequently, σ ξ (s)ν ξ,s is left-continuous in s with respect to weak-* convergence. Now let χ ∈ C ∞ (R) with χ = 1 on (-∞, 0], χ = 0 on [1, ∞), and 0 ≤ χ ≤ 1, and define for ρ > 0

χ ρ (x) = χ π ξ (x)-t ρ , σ ρ = χ ρ σ, µ ρ = χ ρ div σ, σ ρ ξ,s = 1 ρ χ π ξ (•)-t ρ σ ξ (s)ν ξ,s .
In the distributional sense we have

div σ ρ = µ ρ + σ ρ ξ,s ⊗ ds so that for any compactly supported θ ∈ C ∞ (R n ) we have R n ∇θ • dσ ρ + R n θ dµ ρ = - R H ξ,s
θ dσ ρ ξ,s ds .

Letting ρ → 0 and using the left-continuity of σ ξ (s)ν ξ,s in s we arrive at (4).

We now define the slice of a divergence measure vector field as the measure obtained via disintegration with respect to the one-dimensional Lebesgue measure.

Definition 6 (Sliced sets, functions, and measures). Let ξ ∈ S n-1 and t ∈ R.

1. For A ⊂ R n we define the sliced set A ξ,t = A ∩ H ξ,t .

For

f : A → R we define the sliced function f ξ,t : A ξ,t → R, f ξ,t = f | A ξ,t . For f : A → R n we define f ξ,t : A ξ,t → R n , f ξ,t = ξ • f | A ξ,t .
3. We define the sliced measure of a compactly supported divergence measure vector field σ as σ ξ,t = σ ξ (t) ν ξ,t .

By Lemma 3.1 it holds σ • ξ = σ ξ,t ⊗ dt.
Remark 5 (Properties of sliced functions and measures). 1. By Fubini's theorem it follows that for any function f of Sobolev-type W m,p the corresponding sliced function f ξ,t is well-defined and also of Sobolev-type W m,p for almost all ξ ∈ S n-1 and t ∈ R. For the same reason, strong convergence f j → j→∞ f in W m,p implies strong convergence (f j ) ξ,t → f ξ,t in W m,p on the sliced domain.

2. The definitions of sliced functions and measures are consistent in the following sense. If we identify a Lebesgue function f with the measure χ = f L for L the Lebesgue measure, then the same identification holds between f ξ,t and χ ξ,t for almost all ξ ∈ S n-1 and t ∈ R.

3. Let σ be a divergence measure vector field, then the properties [1, Thm. 2.28] of the disintegration σ

• ξ = ν ξ,t ⊗ π ξ # |σ • ξ|(t) = ν ξ,t ⊗ σ ξ (t) dt = σ ξ,t ⊗ dt
immediately imply the following. The map t → σ ξ,t M is integrable and satisfies

R σ ξ,t M dt = R σ ξ (t) dt = σ • ξ M .
Furthermore, for any measurable function f : R n → R, absolutely integrable with respect to |σ • ξ|, it holds

R n f (x) dσ•ξ = R H ξ,t f (x) dν ξ,t (x) dπ ξ # |σ•ξ|(t) = R H ξ,t f (x) dσ ξ,t (x) dt .
We briefly relate our definition of sliced measures to other notions of slices from the literature.

Remark 6 (Notions of slices).

1. Let Lip(A) denote the set of bounded Lipschitz functions on A ⊂ R n . An alternative definition of the slice of a divergence measure vector field σ was introduced by Šilhavý [START_REF] Šilhavý | Divergence measure vectorfields: their structure and the divergence theorem[END_REF] as the linear operator

σ ξ,t : Lip(H ξ,t ) → R, σ ξ,t (ϕ| H ξ,t ) = lim δ 0 1 δ {x∈R n | t-δ<x•ξ<t} ϕξ • dσ
(5) for all ϕ ∈ Lip(R n ) (the right-hand side is well-defined and only depends on ϕ| H ξ,t [START_REF] Šilhavý | Divergence measure vectorfields: their structure and the divergence theorem[END_REF]Thm. 3.5 & Thm. 3.6]). This σ ξ,t equals the so-called normal trace of σ on H ξ,t (see [START_REF] Šilhavý | Divergence measure vectorfields: their structure and the divergence theorem[END_REF] for its definition and properties). In general it is not a measure but continuous on Lip(H ξ,t ) in the sense

σ ξ,t (ϕ) ≤ ( σ M + div σ M ) ϕ W 1,∞ for all ϕ ∈ Lip(H ξ,t ) .
2. Interpreting a divergence measure vector field as a 1-current or a flat 1chain, Šilhavý's definition of σ ξ,t is identical to the classical slice of σ on H ξ,t as for instance defined in [START_REF] White | Rectifiability of flat chains[END_REF] or [16, 4.3.1] (note that Šilhavý's definition corresponds to [28, (3.8)], whose analogue for currents is [16, 4.3.2(5)]).

3. Our notion of a sliced measure from Definition 6 is equivalent to both above-mentioned notions. Indeed, (4) implies

σ ξ,t = (div σ) {x • ξ < t} -div(σ {x • ξ < t}) ,
which shows that the sliced measure represents the normal flux through the hyperplane H ξ,t = {x • ξ = t}. This, however, is the same characterization as given in [28, (3.6)] and [16, 4.2.1] for both above notions of slices.

We conclude the section with several properties needed for the Γ -lim inf inequality. The following result makes use of the Kantorovich-Rubinstein norm (see for instance [21, eq. ( 2) & ( 5)]; in geometric measure theory it is known as the flat norm) on M(R n ), defined by

µ KR = inf{ µ 1 M + µ 2 M | µ 1 ∈ M(R n ), µ 2 ∈ M(R n ; R n ), µ = µ 1 + div µ 2 } = sup Ω f dµ f Lipschitz with constant 1, |f | ≤ 1 .
For measures of uniformly bounded support and uniformly bounded mass it is known to metrize weak- * convergence (see for instance [8, Rem. 2.29(3)-( 4)]). We will furthermore make use of the following fact. Let T s : x → x -sξ be the translation by s in direction -ξ. It is straightforward to check that for any divergence measure vector field

µ ∈ M(R n ; R n ) we have div(π H ξ,t # (µ -µ • ξ ξ)) = π H ξ,t # ( div µ) .
As a consequence, for any µ ∈ M(H ξ,t ) and ν ∈ M(H ξ,t+s ) we have

µ -ν KR ≥ µ -T s# ν KR .
Indeed, let δ > 0 arbitrary and

µ 1 ∈ M(R n ), µ 2 ∈ M(R n ; R n ) with µ -ν = µ 1 + div µ 2 such that µ -ν KR ≥ µ 1 M + µ 2 M -δ, then μ1 = π H ξ,t # µ 1 and μ2 = π H ξ,t # (µ 2 -µ 2 • ξ ξ) satisfy µ -T s# ν = μ1 + div μ2 and thus µ -T s# ν KR ≤ μ1 M + μ2 M ≤ µ 1 M + µ 2 M ≤ µ -ν KR + δ .
Theorem 3.1 (Weak convergence of sliced measures). Let σ j * σ as j → ∞ for a sequence {σ j } of compactly supported divergence measure vector fields with uniformly bounded div σ j M . Then for almost all ξ ∈ S n-1 and t ∈ R we have

σ j ξ,t * σ ξ,t .
Proof. It suffices to show σ j ξ,t * σ ξ,t for a subsequence.

Consider the measures ν j = |σ j | + | div σ j |. Since ν j M is uniformly bounded, a subsequence converges weakly- * to some compactly supported nonnegative ν ∈ M(R n ) (the subsequence is still indexed by j). For I ⊂ R introduce the notation H ξ,I = t∈I H ξ,t . Then for almost all t ∈ R, ν(H ξ,[t-s,t+s] ) → 0 as well as (|σ| + | div σ|)(H ξ,[t-s,t+s] ) → 0 as s 0. For such a t we show convergence of σ j ξ,t -σ ξ,t to zero in the Kantorovich-Rubinstein norm which implies weak- * convergence. To this end fix some arbitrary δ > 0. Given

ζ > 0 let ρ ζ = ρ(•/ζ)/ζ for a nonnegative smoothing kernel ρ ∈ C ∞ (R) with support in [-1, 1]
and R ρ dt = 1. For any compactly supported divergence measure vector field λ we now define the convolved slice λ ξ,ζ,t by

H ξ,t g dλ ξ,ζ,t = R ρ ζ (-s) H ξ,t g dT s# λ ξ,t+s ds = R ρ ζ (-s) H ξ,t g • T s dλ ξ,t+s ds ∀g ∈ C(H ξ,t ) ,
where T s : x → x -sξ is the translation by s in direction -ξ. By Remark 5 [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] we have σ j ξ,ζ,t * σ ξ,ζ,t . Furthermore, there exist ζ > 0 and J ∈ N such that σ ξ,t -σ ξ,ζ,t KR ≤ δ 3 and σ j ξ,t -σ j ξ,ζ,t KR ≤ δ 3 for all j ≥ J. Indeed, for a compactly supported divergence measure vector field λ we have

λ ξ,t -λ ξ,ζ,t KR ≤ R ρ ζ (-s) λ ξ,t -T s# λ ξ,t+s KR ds ≤ R ρ ζ (-s) λ ξ,t -λ ξ,t+s KR ds = R ρ ζ (-s) div(λ H ξ,[t,t+s) ) -(div λ) H ξ,[t,t+s) KR ds ≤ R ρ ζ (-s) |λ|(H ξ,[t,t+s) ) + | div λ|(H ξ,[t,t+s) ) ds ≤ |λ|(H ξ,[t-ζ,t+ζ] ) + | div λ|(H ξ,[t-ζ,t+ζ] ) ,
where in the equality we employed Remark 6(3). Thus, we can simply pick

ζ such that |σ|(H ξ,[t-ζ,t+ζ] ) + | div σ|(H ξ,[t-ζ,t+ζ] ) ≤ δ 3 and ν(H ξ,[t-ζ,t+ζ] ) ≤ δ 6 , while we choose J such that (ν j -ν)(H ξ,[t-ζ,t+ζ] ) ≤ δ
6 for all j > J. Now let J ≥ J such that σ j ξ,ζ,t -σ ξ,ζ,t KR ≤ δ 3 for all j ≥ J, then we obtain

σ j ξ,t -σ ξ,t KR ≤ σ j ξ,t -σ j ξ,ζ,t KR + σ j ξ,ζ,t -σ ξ,ζ,t KR + σ ξ,ζ,t -σ ξ,t KR ≤ δ
for all j > J. The arbitrariness of δ concludes the proof.

Remark 7 (Flat convergence of sliced currents). The convergence from Theorem 3.1 is consistent with the following property of slices of 1-currents: If σ j , j ∈ N, is a sequence of 1-currents of finite mass with σ j → σ in the flat norm, then (potentially after choosing a subsequence) σ j ξ,t → σ ξ,t in the flat norm for almost every ξ ∈ S n-1 , t ∈ R (see [13, step 2 in proof of Prop. 2.5] or [START_REF] White | Rectifiability of flat chains[END_REF]Sec. 3]).

Remark 8 (Characterization of sliced measures).

1. Let the compactly supported divergence measure vector field σ be countably 1-rectifiable, that is, σ = θmH 1 S for a countably 1-rectifiable set S ⊂ R n and H 1 Smeasurable functions m : S → [0, ∞) and θ : S → S n-1 , tangent to S H 1 -almost everywhere. Then the coarea formula for rectifiable sets [START_REF] Federer | Geometric measure theory[END_REF]Thm. 3.2.22]

implies |θ • ξ|H 1 S = H 0 S ξ,t ⊗ H 1 (t) so that R n f dσ • ξ = S f mθ • ξ dH 1 = R S ξ,t f m sgn(ξ • θ) dH 0 dt
for any Borel function f . Hence, for almost all t, σ ξ,t = sgn(ξ • θ) mH 0 S ξ,t .

The choice f = τ (m) m sgn(ξ • θ) yields S τ (m)|θ • ξ| dH 1 = R S ξ,t
τ (m) dH 0 dt .

2.

Let the compactly supported divergence measure vector field σ be H 1diffuse, that is, it is singular with respect to the one-dimensional Hausdorff measure on any countably 1-rectifiable set. Then for almost all ξ ∈ S n-1 and t ∈ R, σ ξ,t is H 0 -diffuse, that is, it does not contain any atoms. Indeed, let σ ξ,t have an atom at x ∈ H ξ,t , then

x ∈ Θ(σ) = x ∈ R n lim inf ρ 0 |σ|(B ρ (x))/ρ > 0 ,
where B ρ (x) denotes the open ball of radius ρ centred at x. This can be deduced as follows. Let φ ∈ C ∞ (R) be smooth and even with support in (-1, 1) and φ(0) = sgn(σ ξ,t ({x})). Further abbreviate K = max x∈R |φ (x)| > 0 and φ ρ = φ(| • -x|/ρ) for any ρ > 0. Equation ( 4) now implies

H ξ,t φ ρ dσ ξ,t = {ξ•x<t} φ ρ d div σ + {ξ•x<t} ∇φ ρ • dσ ≤ {ξ•x<t} φ ρ d div σ + K |σ|(B ρ (x)) ρ .
Taking on both sides the limit inferior as ρ → 0 we obtain |σ ξ,t |({x}) ≤ K lim inf ρ 0 |σ|(B ρ (x))/ρ, as desired.

As a result, for a given ξ the set of t such that σ ξ,t is not H 0 -diffuse is a subset of π ξ (Θ). Thus it remains to show that for almost all ξ ∈ S n-1 the set π ξ (Θ) is a Lebesgue-nullset. Writing

Θ = p∈N Θ p for Θ p = x ∈ R n lim inf ρ 0 σ(B ρ (x))/ρ ≥ 1 p ,
it actually suffices to show that π ξ (Θ) is a Lebesgue-nullset for any p ∈ N. Now by the properties of the 1-dimensional density of a measure [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]Thm. 256],

H 1 (Θ p ) ≤ p 2 |σ|(Θ p )
so that Θ p can be decomposed into a countably 1-rectifiable and a purely 1-unrectifiable set [1, p. 83],

Θ p = Θ r p ∪ Θ u p (Θ u p purely 1-unrectifiable means H 1 (Θ u p ∩ f (R)) = 0 for any Lipschitz f : R → R n ).
By the H 1 -diffusivity assumption on σ we have (abbreviating the Lebesgue measure by L)

L(π ξ (Θ r p )) ≤ H 1 (Θ r p ) ≤ p 2 |σ|(Θ r p ) = 0 ,
and by a result due to Besicovitch [1, Thm. 2.65] we have

L(π ξ (Θ u p )) = 0
for almost all ξ ∈ S n-1 . Thus, for almost all ξ ∈ S n-1 we have L(π ξ (Θ p )) = 0, as desired.

Remark 9 (Characterization of divergence measure vector fields). By a result due to Smirnov [START_REF] Smirnov | Decomposition of solenoidal vector charges into elementary solenoids, and the structure of normal one-dimensional flows[END_REF], any divergence measure vector field σ can be decomposed into simple oriented curves σ γ = γ # γds [0, 1] with γ : [0, 1] → R n a Lipschitz curve and ds the Lebesgue measure, that is,

σ = J σ γ dµ σ (γ)
with J the set of Lipschitz curves and µ σ a nonnegative Borel measure. The results of this section can alternatively be derived by resorting to this characterization, since the slice of a simple oriented curve σ γ can be explicitly calculated.

The Γ -lim inf inequality

We now prove the desired lim inf-inequality, which as already anticipated will be obtained by slicing.

Proposition 3.2 (Γ -lim inf of phase field functional). Let µ + , µ -∈ P(Ω).

We have

Γ -lim inf ε→0 E µ+,µ- ε = E µ+,µ-
with respect to weak- * convergence in M(Ω; R 2 ) and strong convergence in

L 1 (Ω) N .
Proof. Let (σ ε , ϕ ε 1 , . . . , ϕ ε N ) converge to (σ, ϕ 1 , . . . , ϕ N ) in the considered topology. We first extend σ ε and σ to R 2 \ Ω by zero and ϕ ε i and ϕ i to R 2 \ Ω by 1 for i = 1, . . . , N . The phase field cost functional and the cost functional are extended to R 2 in the obvious way (their values do not change). Without loss of generality (potentially after extracting a subsequence) we may assume lim ε→0 E µ+,µε [σ ε , ϕ ε 1 , . . . , ϕ ε N ] to exist and to be finite (else there is nothing to show). As a consequence we have div σ ε = µ ε + -µ ε -as well as div σ = µ + -µ - and ϕ 1 ≡ . . . ≡ ϕ N ≡ 1 (since the phase field cost functional is bounded below by

N i=1 βi 2ε ϕ ε i -1 2 L 2 ). Now let A ⊂ R 2 open
and bounded; the restriction of the phase field cost functional to a domain A will be denoted E µ+,µε [•; A]. Choosing some ξ ∈ S 1 , by Fubini's decomposition theorem we have

E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ; A] = ∞ -∞ A ξ,t ω ε α 0 , γ ε (x) ε , |σ ε (x)| + N i=1 β i 2 ε|∇ϕ i | 2 + (ϕ i -1) 2 ε ξ,t dx dt ≥ ∞ -∞ A ξ,t ω ε α 0 , (γ ε ) ξ,t ε , |σ ε ξ,t | + N i=1 β i 2 ε|(ϕ ε i ) ξ,t | 2 + ((ϕ ε i ) ξ,t -1) 2 ε dx dt = ∞ -∞ G ε [σ ε ξ,t , (ϕ ε 1 ) ξ,t , . . . , (ϕ ε N ) ξ,t ; A ξ,t ] dt ,
where G ε is the dimension-reduced phase field energy from Definition 4 and for simplicity we identified the domain A ξ,t of the sliced functions with an open subset of the real line. Fatou's lemma thus implies

lim inf ε→0 E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ; A] ≥ ∞ -∞ lim inf ε→0 G ε [σ ε ξ,t , (ϕ ε 1 ) ξ,t , . . . , (ϕ ε N ) ξ,t ; A ξ,t ] dt .
By assumption, the left-hand side is finite so that the right-hand side integrand is finite for almost all t ∈ R as well. Pick any such t and pass to a subsequence such that lim inf turns into lim. Indeed σ ε ξ,t * σ ξ,t for every ξ and almost all t, as σ ε * σ and Theorem 3.1. Thus, Corollary 3.1 on the reduced dimension problem implies

lim inf ε→0 G ε [σ ε ξ,t , (ϕ ε 1 ) ξ,t , . . . , (ϕ ε N ) ξ,t ; A ξ,t ] ≥ G[σ ξ,t , (ϕ 1 ) ξ,t , . . . , (ϕ N ) ξ,t ; A ξ,t ]
for almost all t ∈ R so that

lim inf ε→0 E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ; A] ≥ ∞ -∞ G[σ ξ,t , (ϕ 1 ) ξ,t , . . . , (ϕ N ) ξ,t ; A ξ,t ] dt .
For notational convenience let us now define the auxiliary function κ, defined for open subsets A ⊂ R 2 , as

κ(A) = lim inf ε→0 E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ; A] .
Furthermore, introduce the nonnegative Borel measure

λ(A) = α 0 |σ ⊥ |(A) + Sσ∩A τ (m σ ) dH 1
as well as the |σ|-measureable Borel functions

ψ j : R 2 → R, ψ j = σ |σ| • ξ j
for some sequence ξ j , j ∈ N, dense in S 1 . Since σ is a divergence measure vector field, we have

κ(A) ≥ ∞ -∞ G[σ ξ j ,t , (ϕ 1 ) ξ j ,t , . . . , (ϕ N ) ξ j ,t ; A ξ j ,t ] dt = ∞ -∞ α 0 |(σ ξ j ,t ) ⊥ |(A ξ j ,t ) + Sσ ξ j ,t ∩A ξ j ,t τ (|m σ ξ j ,t |) dH 0 dt = α 0 |σ ⊥ • ξ j |(A) + Sσ∩A τ (|m σ |)|θ σ • ξ j | dH 1 ≥ A ψ j dλ
for all j ∈ N where we used Remark 8 in the last equality. By [START_REF] Braides | Approximation of free-discontinuity problems[END_REF]Prop. 1.16] the above inequality implies

κ(A) ≥ A sup j ψ j dλ
for any open A ⊂ R 2 . In particular, choosing A as the 1-neighbourhood of Ω we obtain

lim inf ε→0 E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ; Ω] = κ(A) ≥ A sup j ψ j dλ = α 0 |σ ⊥ |(A) + Sσ∩A τ (m σ ) dH 1 = E µ+,µ-[σ] ,
the desired result.

Equicoercivity

Proof of Theorem 2.2. Due to

C > E µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ] ≥ βi 2ε ϕ ε i -1 2
L 2 for all i = 1, . . . , N , we have ϕ ε i → 1 in L 2 (Ω) and thus also in L 1 (Ω). Furthermore, we will show further below that σ ε L 1 = σ ε M is uniformly bounded, which by the Banach-Alaoglu theorem implies existence of a weakly-* converging subsequence (still denoted σ ε ) with limit σ ∈ M(Ω; R 2 ). It is now a standard property of Γ-convergence that, due to the above equicoercivity, any sequence of minimizers of E µ+,µε contains a subsequence converging to a minimizer of E µ+,µ-.

To finish the proof we show uniform boundedness of

σ ε in M(Ω; R 2 ). Indeed, using ω ε (α 0 , γε(x) ε , |σ(x)|) ≥ α0 2 |σ(x)| for x ∈ K ε 0 (remember that K ε 0 = ∅ for α 0 = ∞) we obtain σ ε M = N i=0 K ε i |σ ε | dx ≤ 2C α 0 + N i=1 K ε i |σ ε | dx
(the first term is interpreted as zero for α 0 = ∞). Furthermore, by Hölder's inequality we have

K ε i |σ ε | dx 2 ≤ K ε i (ϕ ε i ) 2 + α 2 i ε 2 /β i 2ε |σ ε | 2 dx K ε i 2ε (ϕ ε i ) 2 + α 2 i ε 2 /β i dx ≤ C K ε i 2ε (ϕ ε i ) 2 + α 2 i ε 2 /β i dx .
Choosing now some arbitrary λ ∈ (0, 1) we can estimate

K ε i 2ε (ϕ ε i ) 2 + α 2 i ε 2 /β i dx = K ε i ∩{ϕ ε i <λ} 2ε (ϕ ε i ) 2 + α 2 i ε 2 /β i dx+ K ε i ∩{ϕ ε i ≥λ} 2ε (ϕ ε i ) 2 + α 2 i ε 2 /β i dx ≤ 4 α 2 i (1 -λ) 2 K ε i ∩{ϕ ε i <λ} β i (1 -ϕ ε i ) 2 2ε dx+ 2ε λ 2 H 2 (Ω) ≤ 4C α 2 i (1 -λ) 2 + 2ε λ 2 H 2 (Ω) . Summarizing, σ ε M ≤ C 2α0 + N i=1 4C 2 α 2 i (1-λ) 2 + 2εC λ 2 H 2 (Ω).
3.5 The Γ -lim sup inequality Proposition 3.3 (Γ -lim sup of phase field functional). Let µ + , µ -∈ P(Ω) be an admissible source and sink. We have

Γ -lim sup ε→0 E µ+,µ- ε = E µ+,µ-
with respect to weak- * convergence in M(Ω; R 2 ) and strong convergence in

L 1 (Ω) N .
Proof. Consider a mass flux σ between the measures µ + and µ -. We will construct a recovery sequence (σ ε , ϕ ε 1 , . . . , ϕ ε N ) such that σ ε * σ and ϕ ε 1 → 1, . . . , ϕ ε N → 1 in the desired topology as ε → 0 as well as lim sup ε→0 E

µ+,µ- ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ] ≤ E µ+,µ-[σ, 1, . . . , 1]
. Without loss of generality we may restrict our attention to fluxes for which

E µ+,µ-[σ] ≤ C < ∞
since otherwise there is nothing to prove. By [8, Def. 2.2 & Prop. 2.32] there exists a sequence

σ j = Mj k=1 m k,j θ k,j H 1 Σ k,j
of polyhedral divergence measure vector fields in Ω such that

σ j * σ, div σ j = µ j + -µ j -, µ j ± * µ ± , E µ j + ,µ j -[σ j ] → E µ+,µ-[σ].
If µ + and µ -are finite linear combinations of Dirac masses (which we have assumed in the case α 0 = ∞), we may even choose µ j ± = µ ± . We will construct the recovery sequence based on those polyhedral divergence measure vector fields. In the following we will use the notation

E ε [σ, ϕ 1 , . . . , ϕ N ] = Ω ω ε α 0 , γ ε (x) ε , |σ(x)| dx + N i=1 β i L ε [ϕ i ]
for the phase field cost functional without divergence constraints.

Step 1. Initial construction for a single polyhedral segment In this step we approximate a single line element m k,j θ k,j H 1 Σ k,j of σ j by a phase field version. To this end we fix j and 0 ≤ k ≤ M j and drop these indices from now on in the notation for the sake of legibility. Without loss of generality we may assume Σ = [0, L] × {0}, m > 0, and θ = e 1 the standard Euclidean basis vector. Set ῑ = argmin{α i m + β i | i = 0, . . . , N } to identify the phase field that will be active on Σ (ῑ = 0 means that no phase field is active). We first specify (a preliminary version of) the vector field σ ε . To this end let d Σ denote the distance function associated with Σ and define the width

a ε ῑ = α ῑmε 2β ῑ if ῑ > 0 and a ε ῑ = α 0 mε otherwise (6)
over which the vector field will be diffused. We now set

σ ε = m 2a ε ῑ χ {dΣ≤a ε ῑ } e 1 ,
where χ A shall denote the characteristic function of a set A. Note that this vector field encodes a total mass flux of m that is evenly spread over the a ε ῑenlargement of Σ. The corresponding active phase field will be zero in that region. Indeed, consider the auxiliary Cauchy problem whose solution φ ε (t) = 1-exp -t ε represents the well-known optimal Ambrosio-Tortorelli phase field profile. Then we set ϕ ε i = 0 for all i = ῑ and, if ῑ = 0,

φ = 1 ε (1 -φ), φ (
ϕ ε ῑ (x) = φ ε (max{0, d Σ (x) -a ε ῑ }) = 0, if d Σ (x) < a ε ῑ , 1 -exp a ε ῑ -dΣ(x) ε , if d Σ (x) ≥ a ε ῑ .
Let us now evaluate the corresponding phase field cost. In the case ῑ = 0 (which can only occur for α 0 < ∞) we obtain

E ε [σ ε (x), ϕ ε 1 , . . . , ϕ ε N ] = Ω ω ε α 0 , γ ε ε , |σ ε | dx ≤ Ω α 0 |σ ε | + ε p |σ ε | 2 dx = {dΣ≤a ε 0 } 1 + ε p-1 2α 2 0 2ε dx = (α 0 m L+π α 2 0 m 2 ε)(1+ ε p-1 2α 2 0 ) = α 0 m L+C(m, L)ε q ,
where we abbreviated q = min{1, p-1} > 0 and C(m, L) > 0 denotes a constant depending on m and L. In the case ῑ = 0 we have |σ ε | = β ῑ/(α ῑε) as well as

γ ε = α 2 ῑ ε 2 /β ῑ on the support of |σ ε | so that Ω ω ε α 0 , γ ε ε , |σ ε | dx = α 2 ῑ ε 2 β ῑ β ῑ α ῑε 2 1 2ε + ε p-2 β ῑ α ῑ 2 (2 a ε ῑ L + π (a ε ῑ ) 2 ) = α ῑ 2 mL + C(m, L)ε q .
Furthermore we have L ε (ϕ ε i ) = 0 for i = ῑ and, employing the coarea formula,

β ῑL ε (ϕ ε ῑ ) = β ῑ 2 Ω ε|∇ϕ ε ῑ | 2 + (1 -ϕ ε ῑ ) 2 ε dx = β ῑ 2ε 2La ε ῑ + π(a ε ῑ ) 2 + β ῑ 2 ∞ a ε ῑ {dΣ=t} ε|φ ε (t -a ε ῑ )| 2 + (φ ε (t -a ε ῑ ) -1) 2 ε dH 1 (x) dt = β ῑ 2ε 2La ε ῑ + π(a ε ῑ ) 2 + β ῑ ∞ a ε ῑ 1 ε exp 2(a ε ῑ -t) ε (2L + 2πt) dt = α ῑ 2 m + β ῑ L + C(m, L)ε .
Summarizing,

E ε [σ ε , ϕ ε 1 , . . . , ϕ ε N ] ≤ τ (m) L + C(m, L)ε q .
Step 2. Adapting sources and sinks of all polyhedral segments The vector field constructions σ ε k,j from the previous step for each polyhedral segment Σ k,j are not compatible with the divergence constraint associated with the measure σ j , that is, div

  Mj k=1 σ ε k,j   = ρ ε * µ j + -µ j -.
We remedy this by adapting the source and sink of each σ ε k,j . Set

r(j) = max k=1,...,Mj |m k,j | • α1 β1 if α 0 = ∞, max α1
β1 , α 0 else, then all vector fields σ ε k,j have support in a band around Σ k,j of width no larger than r(j)ε. Without loss of generality we assume r(j) ≥ 1 (else we just increase it). Again we concentrate on a single segment with fixed j and k and drop these indices in the following (we will also write r instead of r(j)). Denote by s + and s -the starting and ending point of the segment Σ with respect to the orientation induced by θ. Consider the elliptic boundary value problems

∆u ± (x) = ±m (ρ * δ 0 )(x) on B r (0), ∇u ± (x) • ν(x) = εσ ε (εx + s ± ) • ν(x) on ∂B r (0),
where δ y denotes a Dirac mass centered at y, B r (y) denotes the open ball of radius r around y, and ν denotes the outer unit normal to ∂B r (0). Note that the boundary value problems and their solutions u + and u -are independent of ε due to the definition of σ ε . Setting

σε (x) = ∇u ± ((x -s ± )/ε)/ε if x ∈ B εr (s ± ), σ ε (x) else,
(where we assume ε small enough such that B εr (s + ) and B εr (s -) do not intersect) it is straightforward to check

div(σ ε ) = m ρ ε * (δ s --δ s -) .
Furthermore, to have at least one phase field zero on the new additional support B εr (s + ) ∪ B εr (s -) of the vector field we set and (ϕ ε i ) j = min{(ϕ ε i ) j,1 , . . . , (ϕ ε i ) j,Mj } for i = 1, . . . , N . Obviously, we have, as desired,

ϕ ε 1 (x) = min ϕ ε 1 (x), P (|x -s -|), P (|x -s + |) with P (t) = 0 if t ≤ r ε,
div σj = div(ρ ε * σ j ) = ρ ε * µ j + -µ j -.
Let us now estimate the costs. Let us assume that ε is small enough so that all balls B εr(j) (s ± k,j ) are disjoint as are the supports supp σε k,j \ (B εr(j) (s + k,j ) ∪ B εr(j) (s - k,j )) for all k. An upper bound can then be achieved via

E ε [σ ε j , (ϕ ε 1 ) j , . . . , (ϕ ε N ) j ] ≤ Mj k=1 E ε [σ ε j,k , (ϕ ε 1 ) k,j , . . . , (ϕ ε N ) k,j ] + B εr(j) (s + j,k )∪B εr(j) (s - j,k ) ω ε α 0 , (γ ε ) j ε , |σ ε j | dx + β 1 2 (B 2εr(j) (s + j,k )∪B 2εr(j) (s - j,k ))∩{(ϕ ε 1 )j <(ϕ ε 1 )j } ε|∇(ϕ ε 1 ) j | 2 + 1 ε ((ϕ ε 1 ) j -1) 2 dx .
The last summand can be bounded above by

β 1 B 2εr(j) (0) ε|P (|x|)| 2 + 1 ε (P (|x|) -1) 2 dx ≤ Cε
for some constant C > 0. For the second summand, note that (γ ε ) j ≤ α 2 1 ε 2 /β 1 on B εr(j) (s ± k,j ) due to (ϕ ε 1 ) j = 0 there; furthermore,

ω ε (α 0 , (γ ε ) j /ε, |σ ε j |) ≤ (γ ε ) j 2ε + ε p |σ ε j | 2 ≤ ( α 2 1 ε 2β 1 + ε p )|σ ε j | 2 .
Thus, if we set S ± = {l ∈ {1, . . . , M j } | s ± j,l = s} for fixed s = s + k,j or s = s - k,j

we have

B εr(j) (s) ω ε α 0 , (γ ε ) j ε , |σ ε j | dx ≤ Cε B εr(j) (s) |σ ε j | 2 dx = Cε B εr(j) (s) l∈S + ∇u + j,l ( x-s ε ) ε + l∈S - ∇u - j,l ( x-s ε ) ε 2 dx ≤ CεM j   l∈S + B εr(j) (s) ∇u + j,l ( x-s ε ) ε 2 dx + l∈S -B εr(j) (s) ∇u - j,l ( x-s ε ) ε 2 dx   = CεM j l∈S + B r(j) (0) ∇u + j,l (x) 2 dx + l∈S -B r(j) (s) ∇u - j,l (x) 2 dx = C(s, σ j )ε
for some constant C(σ j ) > 0 depending on the polyhedral divergence measure vector field σ j and the considered point s. In summary, we thus have

E ε [σ ε j , (ϕ ε 1 ) j , . . . , (ϕ ε N ) j ] ≤ Mj k=1 τ (m j,k )H 1 (Σ j,k )+C(m j,k , L j,k )ε q +C(s + j,k , σ j )ε+C(s - j,k , σ j )ε+Cε ≤ E µ j + ,µ j
-[σ j ]+C(σ j )ε q j(ε). Indeed, for a monotonic sequence ε 1 , ε 2 , . . . approaching zero we set j(ε 1 ) = 1 and

j(ε i+1 ) = j(ε i ) if C(σ j(εi)+1 )ε q i > 1 j(εi)+1 , j(ε i ) + 1 else. Then j(ε i ) → ∞ and C(σ j(εi) )ε q i → 0 as i → ∞ so that E µ+,µ- εi [σ εi , ϕ εi 1 , . . . , ϕ εi N ] = E εi [σ εi j(εi) , (ϕ εi 1 ) j(εi) , . . . , (ϕ εi N ) j(εi) ] ≤ E µ j(ε i ) + ,µ j(ε i ) - [σ j(εi) ] + Cκ j(εi) + C(σ j(εi) )ε i q ≤ E µ j(ε i ) + ,µ j(ε i ) - [σ j(εi) ] + Cκ j(εi) + 1 j(ε i ) → E µ+,µ-[σ] .

Numerical experiments

Here we describe the numerical discretization with finite elements and the subsequent optimization procedure used in our experiments.

Discretization

The proposed phase field approximation allows a simple numerical discretization with piecewise constant and piecewise linear finite elements. Let T h be a triangulation of the space Ω of grid size h such that Ω = T ∈T h T . Denoting by P m the space of polynomials of degree m, we define the finite element spaces

X 0 h = {v h ∈ L ∞ (Ω) | v h | T ∈ P 0 ∀ T ∈ T h } , X 1 h = {v h ∈ C 0 (Ω) | v h | T ∈ P m ∀ T ∈ T h } . Using as before the notation E ε [σ, ϕ 1 , . . . , ϕ N ] = Ω ω ε α 0 , γ ε (x) ε , |σ(x)| dx + N i=1 β i L ε [ϕ i ] ,
the discretized phase field problem now reads min

(σ,ϕ1,...,ϕ N )∈X 0 h ×(X 1 h ) N ϕ1| ∂Ω =...=ϕ N | ∂Ω =1 E ε [σ, ϕ 1 , . . . , ϕ N ] such that Ω -σ • ∇λ dx = Ω ρ ε * (µ + -µ -)v h dx ∀λ ∈ X 1 h ,
where all integrals are evaluated using midpoint quadrature and the divergence constraint is enforced in weak form. In our numerical experiments we use Ω = (0, 1) 2 with a regular quadrilateral grid whose squares are all divided into two triangles.

Optimization

Here we describe the numerical optimization method used to find a minimizer of the objective functional. We first elaborate the simplest case in which the transport cost τ only features a single affine segment, that is, α 0 = ∞ and N = 1. Afterwards we consider the setting with N > 1 and subsequently also with α 0 < ∞, which requires more care due to the higher complexity of the energy landscape.

Single phase field and no diffuse mass flux (N = 1, α 0 = ∞) In this case the energy reads

E ε [σ, ϕ] = Ω γ ε (x) ε |σ(x)| 2 2 + β 1 2 ε|∇ϕ(x)| 2 + (ϕ(x) -1) 2 ε dx with γ ε (x) = ϕ(x) 2 + α 2 1 ε 2 /β 1 .
The employed optimization method is similar to the one presented in [START_REF] Chambolle | Variational approximation of size-mass energies for k-dimensional currents[END_REF] and updates the variables σ and ϕ alternatingly.

Let us abbreviate f ε = ρ ε * (µ + -µ -). For minimization with respect to σ, we use the dual variable λ ∈ X 1 h to enforce the divergence constraint and write min

σ∈X 0 h Ω σ•∇λ+λfε dx=0 ∀λ∈X 1 h Ω γ ε ε |σ| 2 2 dx = min σ∈X 0 h max λ∈X 1 h Ω γ ε ε |σ| 2 2 -σ • ∇λ -λf ε dx = max λ∈X 1 h min σ∈X 0 h Ω γ ε ε |σ| 2 2 -σ • ∇λ -λf ε dx ,
where the last step follows by standard convex duality. The minimization in σ can be performed explicitly, yielding σ = ε∇λ γε . Inserting this solution leads to a maximization problem in λ,

max λ Ω - ε|∇λ| 2 2γ ε -λf ε dx .
The corresponding optimality conditions,

Ω ε∇λ • ∇µ γ ε dx = - Ω µf ε dx ∀µ ∈ X 1 h , (7)
represent a linear system of equations that can readily be solved numerically for λ so that subsequently σ can be computed (note that Ω f ε dx = 0 so that the above equation has a solution). Fixing σ, the optimality condition for ϕ reads

Ω |σ| 2 ϕψ ε + β 1 ε∇ϕ • ∇ψ + β 1 ε (ϕ -1)ψ dx = 0 ∀ψ ∈ X 1 h with ψ| ∂Ω = 0 , (8)
which can again be solved under the constraint ϕ 1 | ∂Ω = 1 by a linear system solver.

In addition to the alternating minimization a stepwise decrease of the phase field parameter ε is performed, starting from a large value ε start for which the energy landscape shows fewer local minima. Algorithm 1 summarizes the procedure.

Multiple phase fields and no diffuse mass flux (N > 1, α 0 = ∞) Again we aim for an alternating minimization scheme. Fixing ϕ 1 , . . . , ϕ N , the optimization for σ is the same as before, since only γ ε changes. However, the optimization in the phase fields ϕ 1 , . . . , ϕ N is strongly nonconvex (due to the minimum in γ ε ) and thus requires a rather good initialization and care in the alternating scheme.

To avoid minimization for phase field ϕ i with the min-function inside γ ε we perform a heuristic operator splitting: in each iteration of the optimization we

Algorithm 1 Minimization for N = 1, α 0 = ∞ function SPFS(ε start , ε end , N iter , α 1 , β 1 , µ + , µ -, ρ ε end ) set f ε = (µ + -µ -) * ρ ε end , σ 0 = 0 for j = 1, . . . , N iter do set ε j = ε start -(j -1) εstart-ε end
Niter-1 set ϕ j to the solution of (8) for given fixed σ = σ j-1 set γ j ε = ϕ j 2 + α 2 1 ε 2 j /β 1 set λ j to the solution of ( 7) for given fixed γ ε = γ j ε set σ j = εj ∇λ j 2γ j ε end for end function return σ Niter , ϕ Niter , λ Niter first identify at each location which term inside

γ ε = min(ϕ 2 1 +α 2 1 ε 2 /β 1 , . . . , ϕ 2 N + α 2 N ε 2 /β N )
is the minimizer, which is equivalent to specifying the regions

R ε i = {x ∈ Ω | γ ε (x) = ϕ i (x) 2 + α 2 i ε 2 /β i } , i = 1, . . . , N. (9)
Afterwards we then optimize the energy E ε separately for each phase field ϕ i assuming the regions R ε i fixed, that is, we minimize

N i=1 Ω ϕ i (x) 2 + α 2 i ε 2 /β i ε |σ(x)| 2 2 χ R ε i (x) + β i 2 ε|∇ϕ i (x)| 2 + β i 2 (ϕ i (x) -1) 2 ε dx ,
where χ R ε i is the characteristic function of region R ε i . Similarly to the case N = 1 of a single phase field, this amounts to solving the linear system

Ω ϕ i ψ ε |σ| 2 χ R ε i + β i ε∇ϕ i • ∇ψ + β i ϕ i ψ ε dx ∀ψ ∈ X 1 h with ψ| ∂Ω = 0 (10)
for ϕ i ∈ X 1 h with ϕ i | ∂Ω = 1. Since in the above simple approach the regions R ε i and phase fields ϕ i can move, but cannot nucleate within a different region (indeed, imagine for instance R ε 1 = Ω, then ϕ 2 , . . . , ϕ N will be optimized to equal 1 everywhere so that in the next iteration again R ε 1 = Ω), a suitable initial guess is crucial. To provide such a guess for the initial regions R ε i , we proceed as follows. We first generate some flux network σ 0 consistent with the given source and sink. This can for instance be done using the previously described algorithm for just a single phase field: in our simulations we simply ignore ϕ 2 , . . . , ϕ N and pretend only the phase field ϕ 1 would exist (essentially this means we replace the cost function τ by m → α 1 m + β 1 for m > 0; of course, an alternative choice would be to take m → αm + β for some α, β > 0 that better approximate τ for a larger range of values m). We then identify the total mass flowing through each branch of σ 0 . To this end we convolve |σ 0 | with the characteristic function χ Br(0) of a disc of radius r. If r is sufficiently large compared to the width of the support of σ 0 , we obtain χ Br(0

) * |σ 0 | (x) = Br(x) |σ 0 |(y) dy ≈ 2rm(x) ,
where m(x) denotes the mass flux through the nearby branch of σ 0 . Now we can compute the regions

R ε i = {x ∈ Ω | i = argmin j=1,...,N {α j m(x) + β j }} (11)
and furthermore use σ 0 as initial guess of the vector field. Algorithm 2 summarizes the full procedure.

Algorithm 2 Minimization for N > 1, α 0 = ∞ function MPFS(ε start , ε end , N iter , α 1 , . . . , α N , β 1 , . . . , β N , µ + , µ -, ρ ε end ) set f ε = (µ + -µ -) * ρ ε end set σ 0 , •, • = SP F S(ε start , ε end , N iter , α 1 , β 1 , µ + , µ -, ρ ε end ) compute regions R ε 1 , .
. . , R ε N via (11) for j = 1, . . . , N iter do set ε j = ε start -(j -1) εstart-ε end Niter-1 set ϕ j i to the solution of (10) for given fixed σ = σ j-1 , i = 1, . . . , N update regions R ε 1 , . . . , R ε N via (9) set γ j ε = min i=1,...,N (ϕ j i ) 2 + α 2 i ε 2 /β i set λ j to the solution of [START_REF] Braides | Γ-convergence for beginners[END_REF] for given fixed γ ε = γ j ε set σ j = εj ∇λ j 2γ j ε end for end function return σ Niter , ϕ Niter 1 , . . . , ϕ Niter N Note that the estimate m(x) of the flowing mass is only valid in close proximity of the support of σ 0 so that the regions R ε i are only reliable near σ 0 . However, away from σ 0 all phase fields will be close to 1 anyway so that the regions R ε i do not play a role there. The effectiveness of the initialization can be further improved by an energy rescaling which we typically perform in our simulations: Recall that the optimal width (6) of the support of the vector field σ not only depends on the transported mass m(x), but also on which phase field is active at x. Thus, initializing with some vector field σ that was computed based on preliminary active phase fields may erroneously give slight preference to incorrect phase fields. This can be avoided by a small parameter change which assigns a different ε to each phase field. Indeed, setting ε i = β i ε/α i to be the phase field parameter associated with phase field ϕ i , equation [START_REF] Braides | Approximation of free-discontinuity problems[END_REF] shows that the support width of σ becomes mε and thus no longer depends on the phase field. Thus, in practice we usually work with the phase field cost

Ẽε [σ, ϕ 1 , . . . , ϕ N ] = Ω ω ε α 0 , γε (x) ε , |σ(x)| dx + N i=1 β i L εi [ϕ i ] (12)
with γε (x) = min i=1,...,N {ϕ i (x) 2 +α 2 i εε i /β i } = min i=1,...,N {ϕ i (x) 2 +α i ε 2 }. The Γ-convergence result can readily be adapted to this case.

Multiple phase fields and diffuse mass flux (N > 1, α 0 < ∞) The difference to the previous case is that now there may be nonnegligible mass flux σ in regions where no phase field ϕ 1 , . . . , ϕ N is active. Correspondingly, we adapt the previous alternating minimization scheme by introducing the set

R ε 0 = x ∈ Ω | |σ(x)| > α 0 γ ε (x)/ε
, which according to the form of ω ε in Definition 3 describes the region in which mass flux σ is penalized by α 0 |σ|. The regions in which the i th phase field is active are thus modified to

Rε i = R ε i \ R ε 0 . (13)
As before, we now separately minimize

N i=1 Ω ϕ i (x) 2 + α 2 i ε 2 /β i ε |σ(x)| 2 2 χ Rε i (x) + β i 2 ε|∇ϕ i (x)| 2 + β i 2 (ϕ i (x) -1) 2 ε dx , (14)
for each phase field ϕ i . The optimization for σ changes a little compared to the previous case since the problem is no longer quadratic and thus no longer reduces to solving a linear system. Instead we will perform a single step of Newton's method in each iteration. The optimization problem in σ reads min

where

ξ(|σ|) = 1 |σ| ∂ 3 ω ε α 0 , γ ε ε , |σ| = min γ ε ε , α 0 |σ| + 2ε p .
Letting σ and λ be the coefficients of σ and λ in some basis {b 0 i } i of X 0 h and {b 1 i } i of X 1 h , respectively, the optimality conditions can be expressed as

0 = R(σ, λ) = M [ξ(|σ|)] B B T 0 σ λ + 0 F ,
where the finite element matrices and vectors are defined as

M [ξ] ij = Ω ξb 0 i b 0 j dx , B ij = Ω b 0 i ∇b 1 j dx , F i = Ω b 1 i f ε dx .
In each iteration of the alternating minimization scheme we now take one Newton step for 0 = R(σ, λ). As before, the algorithm requires a suitable initial guess, which is determined in the same way as for the case without diffuse component. Algorithm 3 summarizes the alternating scheme.

Algorithm 3 Minimization for N > 1, α 0 < ∞ function MPFSD(ε start , ε end , N iter , α 1 , . . . , α N , β 1 , . . . , β N , µ + , µ -, ρ ε end ) set f ε = (µ + -µ -) * ρ ε end set σ 0 , •, λ 0 = SP F S(ε start , ε end , N iter , α 1 , β 1 , µ + , µ -, ρ ε end ) compute regions R ε 0 , Rε 1 , . . . , Rε N via (13) for j = 1, . . . , N iter do set ε j = ε start -(j -1) εstart-ε end Niter-1 set ϕ j i to the minimizer of ( 14) for given fixed σ = σ j-1 , i = 1, . . . , N update regions R ε 0 , Rε 1 , . . . , Rε N via ( 9) and ( 13) set γ j ε = min i=1,...,N (ϕ j i ) 2 + α 2 i ε 2 /β i set (σ j , λj ) = (σ j-1 , λj-1 ) -DR(σ j-1 , λj-1 ) -1 R(σ j-1 , λj-1 ) for γ ε = γ j-1 ε end for end function return σ Niter , ϕ Niter 1 , . . . , ϕ Niter N

Experimental results

The algorithms were implemented in MATLAB c ○ ; parameters reported in this section refer to the rescaled cost [START_REF] Chambolle | Variational approximation of size-mass energies for k-dimensional currents[END_REF]. We first present simulation results for a single phase field and no diffuse mass flux, N = 1 and α 0 = ∞. Figures 3 and4 show solutions for a source and a number of equal sinks arranged as a regular polygon. If α 1 is small as in Figure 3, the solution looks similar to the Steiner tree (which would correspond to α 1 = 0), while the solutions become much more asymmetric for larger α 1 as in Figure 4. More complex examples are displayed in Figure 5.

Figure 6 shows simulation results for the same source and sink configuration as in Figure 5, but this time with N = 3 different linear segments in τ and corresponding phase fields. It is clear that different phase fields become active on the different network branches according to the mass flux through each branch. This can be interpreted as having streets of three different qualities: the street ϕ 3 allows faster (cheaper) transport, but requires more maintenance than the others, while street ϕ 1 requires the least maintenance and only allows expensive transport.

The case α 0 < ∞ finally can be interpreted as the situation in which mass can also be transported off-road, that is, part of the transport may happen without a street network, thus having maintenance cost β 0 = 0, but at the price of large transport expenses α 0 per unit mass. Corresponding results for again the same source and sink configuration are shown in Figure 7. In contrast to the case α 0 = ∞ it is now also possible to have sources and sinks that are not concentrated in a finite number of points. A corresponding example is shown in Figure 8.

Figure 1 :

 1 Figure 1: Numerically optimized transport networks with τ (m) = 1 + 0.05m for transport from (a) one mass source to five sinks, (b) four mass sources to four mass sinks, (c) a central mass source to 16 sinks around it.

Proposition 2 . 1 (

 21 Existence of minimizers to the phase field functional). The phase field cost functional E µ+,µε has a minimizer (σ, ϕ 1 , . . . , ϕ N ) ∈ X µ+,µε . Proof. The functional is bounded below by 0 and has a nonempty domain. Indeed, choose φ1 ≡ . . . ≡ φN ≡ 1 and σ = ∇ψ, where ψ solves ∆ψ = µ ε + -µ ε in Ω with Neumann boundary conditions ∇ψ • ν ∂Ω = 0, ν ∂Ω being the unit outward normal to ∂Ω. (Since Ω µ ε + -µ ε -dx = 0, a solution ψ exists and lies in W 2,2 (Ω) by standard elliptic regularity.) Obviously, (σ, φ1 , . .

  we have | θ| ≥ |ϑ| and for t C = -∞ we have | θ| ≤ |ϑ| everywhere on C). The cost did not change by this modification since

  dx + τ B∪C |ϑ| dx for all measurable A, B, C ⊂ I (due to the subadditivity of τ ) and in the third inequality we used τ (m) ≤ α 0 m as well as the lower semi-continuity of the mass on an open set. Letting now ζ → 0 (so that by the σ-continuity of ϑ we have |ϑ|(I \ B ζ ({x 1 , . . . , x K })) → |ϑ|(I \ {x 1 , . . . , x K })) we obtain lim inf ε→0

  where B ζ ({xi})\Ci(ε) |ϑ ε | dx decreases to zero by Proposition 3.1. Therefore, Proposition 3.1 implies

1 Figure 2 :

 12 Figure 2: Left: Sketch of the optimal profile of |σ ε | and a phase field ϕ ε ῑ for some ῑ > 0 with m = 2, ε = 0.1, α ῑ = 1, β ῑ = 1. Right: Sketch of the numerical solution to the 1D problem with the same parameters.

t r ε else and ϕ ε 2 = ϕ ε 2 ,

 2 . . . , ϕ ε N = ϕ ε N .Reintroducing now the indices k and j,

  )σ • ψ -∇λ • ψ dx for all ψ ∈ X 0 h , 0 = Ω σ • ∇µ + µf ε dx for all µ ∈ X 1 h ,

Acknowledgements The authors thank Antonin Chambolle and Bénoit Merlet for various discussions and suggestions. B.W.'s research was supported by the Alfried Krupp Prize for Young University Teachers awarded by the Alfried Krupp von Bohlen und Halbach-Stiftung. The work was also supported by the Deutsche Forschungsgemeinschaft (DFG), Cells-in-Motion Cluster of Excellence (EXC1003 -CiM), University of Münster, Germany.

for some constant C(σ j ) depending on σ j .

Step 3. Correction of the global divergence Recall that the vector field σ ε to be constructed has to satisfy div σ ε = ρ ε * (µ + -µ -). In the case α 0 = ∞ the vector field σε j already has that property due to µ j ± = µ ± (thus we set σ ε j = σε j ). However, if α 0 < ∞ (so that admissible sources and sinks µ + and µ -do not have to be finite linear combinations of Dirac masses) we still need to adapt the vector field to achieve the correct divergence. To this end, let

) be the optimal Wasserstein-1 flux between µ j ± and µ ± , that is, λ ± j minimizes λ M among all vector-valued measures λ with div λ = µ j ± -µ ± . Setting

we thus obtain div σ ε j = ρ ε * (µ + -µ -), as desired. The additional cost can be estimated using the fact

as well as

, where µ j ± M is an upper bound for the total mass moved by λ j ± and the constant C > 0 depends on ρ. With those ingredients we obtain

since the Wasserstein-1 distance W 1 (•, •) metrizes weak- * convergence. Furthermore, in the previous steps we have already estimated

Summarizing,

for some constant C > 0 and C(σ j ) > 0 depending only on σ j .

Step 4. Extraction of a diagonal sequence We will set σ ε = σ j (ε) ε , ϕ ε 1 = (ϕ ε 1 ) j(ε) , . . . , ϕ ε N = (ϕ ε N ) j(ε) for a suitable choice