{'fraction_non_alphanumeric': 0.05844091849807825, 'fraction_numerical': 0.024177917939620906, 'mean_word_length': 3.7550765385816933, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Density-constrained Chemotaxis and Hele-Shaw flow

We consider a model of congestion dynamics with chemotaxis, where the density of cells follows the chemical signal it generates, while observing an incompressibility constraint. We show that when the chemical diffuses slowly and attracts the cells strongly, then the dynamics of the congested cells is well approximated by a surface-tension driven free boundary problem. More precisely, we show that in this limit the density of cell converges to the characteristic function of a set whose evolution is described by a Hele-Shaw free boundary problem with surface tension. Our problem is set in a bounded domain, which leads to an interesting analysis on the limiting boundary conditions for the density function. Namely, we prove that the assumption of Robin boundary conditions for the chemical potential leads to a contact angle condition for the free interface.


1.1.
A model for chemotaxis with density constraint. The classical parabolic-elliptic Patlak-Keller-Segel model for chemotaxis reads: where ρ denotes the cell density and φ the concentration of some chemical. The nonnegative parameters µ and η are the cell and chemical diffusivity, χ is the cell sensitivity, and θ and σ describe the production and degradation of the chemical (see [10], [16], [7]).
In this model, the diffusion competes with the aggregating potential φ, leading to the well-known phenomena of concentration and finite time blow-up of the density (see e.g. [4], [8]). In order to investigate the behavior of the density ρ after saturation occurs we take into account the incompressibility of the cells by imposing a constraint ρ ≤ ρ M . We replace ρ with ρ/ρ M and φ with φ/(ρ M θ) so that ρ M = θ = 1 and denoteχ = χρ M θ. We are then led to the equation (see [18,11] for details): where the pressure p is a Lagrange multiplier for the contraint ρ ≤ 1, and satisfies Similar models have been used in particular in the study of congested crowd motion (see [18]). The conditions on p can also be expressed by writing p ∈ P (ρ) with In a companion paper [11], we proved the existence and uniqueness of a weak solution for (1.1) and investigated its relation to some free boundary problems. In this paper, we will investigate the singular limit of strong attraction (χ µ) and small chemical diffusion (η 1) and prove that the model is asymptotically close to a Hele-Shaw free boundary problem with surface tension (1.8). This establishes the first rigorous link between a general Chemotaxis system and Hele-Shaw flow with surface tension (to the best of our knowledge).
We also aim to analyze the behavior of the solutions of (1.1) near a fixed boundary, by setting the problem in a bounded domain Ω ⊂ R d . In particular we are interested in the effect on the dynamics of different absorption rates of the chemical at the boundary. For full generality, we will use Robin boundary conditions for φ with a fixed parameter for absorption rate. For the density, we impose Neumann boundary conditions which ensure the conservation of cell density.
Note that the scaling of the continuity equation can also be obtained by rescaling the time variable so that we observe the evolution of ρ at a time scalet ∼ ε −1 /χ, under the assumption that µ = O(εχ).
The assumption σ > 0 is important. When σ = 0, the functionφ = ε 2 φ is the usual Newtonian potential (up to the boundary condition on ∂Ω), which does not localize in the limit ε → 0: see the discussion below (1.8). By contrast, when σ > 0, we have φ ∼ 1 σ ρ when ε 1 and the effect of ε −1 ∇φ on the dynamic of the saturated regions is akin to that of surface tension.
In other words, in this fully saturated regime, the chemotaxis system (1.3)-(1.4) can be seen as a free boundary problem describing the motion of the region occupied by the cell, driven by the chemical concentration φ and the pressure variable p. Since we obtained (1.3)-(1.4) by imposing the constraint ρ ≤ 1, but without requiring ρ ∈ {0, 1}, it is not clear that we should actually have ρ(x, t) = χ Ωs(t) (x) in general. In [11], we proved that if ρ is a characteristic function at t = 0, then this remains true at positive times for (1.3)-(1.4) when µ = 0 . On the other hand when µ > 0 the density is never a characteristic function. Indeed in this case the saturated set interacts with the unsaturated part of the density by a Richards-type problem, as shown in [11].
Nevertheless, we will show in this paper that, in the limit ε → 0, the effect of the attractive potential is strong enough to ensure the convergence of ρ to a characteristic function χ Ωs(t) (x) for all µ ≥ 0. We will then show that the asymptotic dynamic of Ω s (t) is described by the Hele-Shaw free boundary problem with surface tension where κ denotes the mean curvature of the free boundary Σ(t) (taken to be positive when Ω s (t) is convex). Formally, we can get (1.8) from (1.7) by proving that the quantity ε −1 1 2σ − φ is an approximation of the mean-curvature of Ω s when ε 1.
Note that, as σ tends to zero, the weight on surface tension grows to infinity in (1.8). Thus heuristically we expect that the limit density support will re-adjust itself into a ball at time scale of order σ 3/2 (and instantly when σ = 0). This is consistent with the convergence to radial solutions of (1.1) when σ = µ = 0: see [3] and [7] for further discussions.
1.3. The presence of bounded domain. Our result for µ = 0 bears similarities with [9], where the emergence of surface tension and derivation of a Muskat problem is studied via a variational approximation. In that paper, the potential φ solves φ t − ∆φ = ρ (in R d ) and instead of the Keller-Segel system, the authors considers a discrete-time approximations constructed via a JKO scheme. A similar variational analysis is performed in [12] for the L 2 -based thresholding scheme.
Note that both [9] and [12] consider the setting of periodic torus or entire R d for the interaction energy, in which case φ can be written as a convolution with the heat kernel. Such a representation of φ, as well as the symmetry of the heat kernel in space variables, played an important role in the analysis of the aforementioned papers, in particular when deriving the weak limit equation. The fact that our problem is set in a bounded domain presents an interesting challenge to this analysis. In particular this necessitates a more PDE-oriented proof of Proposition 5.2, replacing corresponding proofs in [12] and [9]. Our result appears to be the first that links a Keller-Segel system with a Hele-Shaw flow with surface tension, regardless of the choice of the domain. This connection was also suggested in the very recent paper [5], where the incompressible limit of a generalized version of Keller-Segel system is investigated. The model is a variant of Cahn-Hilliard equation, which can be seen as a diffuse-interface approximation of our Hele-Shaw flow with surface tension.
Another novel feature of our analysis, also related to the bounded domain, is the characterization of the free boundary behavior near the fixed boundary ∂Ω. Of particular interest, in the context of the singular limit ε → 0, is how the Robin boundary conditions imposed on φ plays a role on the dynamic of Ω s (t). We will show that (1.8) must be supplemented by the contact angle condition (1.9) cos(θ) = γ := − min 1, 2α where θ is the angle formed by the free surface Σ(t) and the fixed boundary ∂Ω, measured from inside of the set Ω s (t), and the fixed boundary ∂Ω at the triple junction Σ(t) ∩ ∂Ω (see Figure 1). In particular, for Neuman condition α = 0 (zero absorption of chemicals), the contact must be orthogonal, while for Dirichlet condition β = 0 (and whenever the absortion rate α is bigger or equal than √ σβ), the contact must be tangential.
1.4. Notations and definitions. Throughout the paper Ω is a smooth bounded domain in R d .
In [11] we prove the existence and uniqueness of a weak solution in the sense of Definition 1.1, using the fact that it is a gradient flow with respect to the Wasserstein metric. Here the free energy is given by This energy structure of the equation will also play a key role in this paper. Because this energy does not behave well when ε 1, we will work instead with the functional Since´Ω ρ dx is preserved by the equation we are only adding a constant to the energy, but this constant is important when ε 1 (it was proved in [15] that F ε (ρ) is bounded uniformly in ε when ρ = χ E ∈ BV(Ω; {0, 1})). The following result was proved in [11]: 11]). For any ε > 0 and any initial condition ρ in satisfying 0 ≤ ρ in ≤ 1 a.e. in Ω, there exists a unique (ρ ε , p ε ) weak solution of (1.3)-(1.4) in the sense of Definition 1.1. Furthermore, ρ ε satisfies the energy inequality with v ε defined as in Definition 1.1.
The goal of this paper is to show that when ε 1, the solution of (1.3)-(1.4) given by Theorem 1.2 converges to the solution of the following Hele-Shaw problem with surface tension: together with the contact angle condition (1.9). Recall that n and ν respectively denote the outward normal of Ω and Ω s (t) at their boundary points.
The definition of a weak solution of (1.13)-(1.9) is parallel to the Definition 1.1: for some s > 0 and the followings hold: for any vector field ξ ∈ C ∞ c (Ω × (0, ∞); R d ) such that ξ · n = 0 on ∂Ω.
This definition, similar to the one given in [12] and [9], warrant several comments.
(3) Note that q has very low regularity in this definition. Since q ∼ κ along Σ(t), we cannot expect much more regularity on q without improving the regularity of the free boundary Σ(t). (4) As in Definition 1.1, (1.14) is simply the weak formulation for the continuity equation ∂ t ρ + div (ρv) = 0 with Neumann boundary conditions and initial data ρ in . Since ρ = χ Ωs(t) , it encodes the velocity law V = v · ν, the incompressibility condition div v = 0 in Ω s (t) and the Neumann condition v · ν = 0 on ∂Ω ∩ Ω s (t). (5) By taking test functions ξ supported in either {ρ = 0} or {ρ = 1}, we see that Equation (1.15) implies ∇q = 0 in {ρ(t) = 0} and v = −∇q in Ω s (t) = {ρ(t) = 1}. Subtracting a constant if needed, we can in particular assume that q = 0 in {ρ(t) = 0}. For general test functions ξ, and taking into account the right hand side of (1.15) we further get the surface tension condition q = κ 4σ 3/2 on Σ(t) and the contact angle condition (1.9). This can be seen by using the classical formula (for a smooth interface Σ): where ν is the normal vector to Σ, κ denotes the mean curvature of Σ and b is the conormal vector along Γ = ∂Σ. Indeed, formally at least, the right hand side of (1.15) is (using the fact that ξ · n = 0 on ∂Ω): where b and c are unit conormal vectors along Γ = ∂Σ ∩ ∂Ω: b is tangent to Σ while c is tangent to ∂Ω (see Figure 1).
Integration by parts in (1.15) thus reveals that the jump of q across Σ must be equal to 1 4σ 3/2 κ. Since q = 0 in {ρ(t) = 0}, we get q = 1 4σ 3/2 κ along Σ = ∂E ∩ Ω. Finally, the cancellation of the lower dimensional integral requires that the component of the vector b − γ c that is tangential to ∂Ω must vanish. In particular, we must have b − γ c · c = 0 and so c · b = γ, which gives the contact angle condition cos θ = γ.
Passing to the limit in the right hand side of this equation to derive (1.15) will be the main result of the second part of the paper (see Proposition 5.2) and is at the heart of the relation between the Hele-Shaw model with active potential (1.6) and the Hele-Shaw model with surface tension (1.13). We will see in particular that the contact angle γ depends on the boundary condition for the potential φ and is given by 1.5. Energy. The proof of the convergence will require some assumptions about the convergence of the energy. Before stating this assumption, we need to recall a few important facts about the singular part of the energy The properties of J ε when ε 1 were studied by two of the authors in [15]. The first important result is the following: Let Ω be a bounded open set with C 1,α boundary. Given a set E ⊂ Ω with finite perimeter P (E, Ω) < ∞, we have We point out that while the result is proved only for σ = 1 in [15], but it can be easily extended to σ = 1 by scaling. More precisely withφ = σφ, which is the equation studied in [15].
Above proposition identifies the limit of J ε (χ E ). However, this functional is not lower-semicontinuous when 2α α+ √ σβ > 1 and cannot be the Γ-limit of J ε . We can in fact prove: Theorem 1.5. Let Ω be a bounded open set with C 1,α boundary. The functional J ε Γ-converges, when ε → 0 to This theorem is proved in [15] when J ε is restricted to characteristic functions, so we show in Appendix C how the proof can be generalized to our more general framework. This extension requires a new formulation of the energy J ε , see (2.2).
This Γ-convergence result suggests that the solution of the gradient flow associated to the energy J ε (which corresponds to equation (1.3)) converges when ε → 0 to a solution of the gradient flow associated to J 0 which is formally a Hele-Shaw flow (1.13)-(1.9). This is indeed the result that we want to make precise in the present paper.
1.6. Main results. We are now able to state the main result: Theorem 1.6. Given an initial data ρ in = χ Ein ∈ BV (Ω; {0, 1}), µ ≥ 0 and a sequence ε n → 0, let (ρ εn , p εn ) be the unique solution of (1.3)-(1.4) given by Theorem 1.2. Then along a subsequence the density ρ εn (x, t) converges strongly in L ∞ ((0, T ); L 1 (Ω)) to and the modified pressure variable q εn (defined by (5.1)) converges to q weak- * in L 2 ((0, T ); (C s (Ω)) * ) for any s > 0. Furthermore, ρ satisfies the continuity equation (1.14) for some velocity function v(x, t) as well as the energy dissipation property Finally, if the following energy convergence assumption holds: then the limit (ρ, q) also satisfies the pressure equation (1.15) on (0, T ). Thus it follows that (ρ, p) is a weak solution of (1.13)-(1.9) in the sense of Definition 1.3, with initial condition ρ in and contact angle The result also holds if we consider a sequence of initial data ρ ε in = χ E ε in bounded in BV (Ω), converging strongly to ρ in = χ Ein in L 1 and satisfying lim J εn (ρ ε in ) → J 0 (ρ in ). The existence of such a sequence, for any finite perimeter set E in , is proved in [15] (Proposition 5.3) as part of the Γ convergence result.
We note that Theorem 1.6 is a conditional result, since it requires the energy convergence assumption (1.20). The analysis of [15] implies that we always have so (1.20) ensures that there is no loss of boundary between phases in the limit. This assumption is rather natural and is similar to the one required for instance in [14,12,9]. It is likely that by proceeding as in [2] one could also obtain a weaker notion of solutions using the theory of varifolds without this assumption, but we do not pursue this direction here. 1.7. Outline of the paper. We begin with deriving two alternative formulas for the energy J ε in the next section, which play a crucial role in our analysis. Section 3 collects the main a priori estimates for the ε-solutions. The proof of Theorem 1.6 is then split between sections 4 and 5. In Section 4, we prove Proposition 4.1 which gives the first part of the theorem, namely the strong convergence in L 1 of the density toward a characteristic function which satisfies the continuity equation (1.14) and the energy inequality (1.19). Section 5 completes the proof of Theorem 1.6 by deriving equation (1.15) under condition (1.20). The main step is Proposition 5.2 which shows that the convergence of the energy (1.20) implies the convergence of the first variation. In the last section, we briefly recall the construction of the JKO scheme used in [11] to prove the existence of weak solutions to (1.3)-(1.4) (Theorem 1.2) and we state a convergence result similar to Theorem 1.6 for a discrete-time approximation: such a result is of independent interest for numerical applications.

Alternate formulas for J ε
A crucial tool in our analysis will be a couple of alternate formula for the energy J ε . We recall that the total energy of the model, F ε , is given by where J ε is defined by (1.17) and plays key role in the analysis when ε 1. For ρ satisfying the constraint 0 ≤ ρ ≤ 1, we have for some kernel G ε . A similar energy functional is used in [12,9] with G ε is the heat kernel in R d . However, we will rely on some different formulations for J ε which make use of the particular equation solved by the function φ ε in our model: First, we write when β = 0. Alternatively, we can write the more symmetric formula (for any α, β ≥ 0 with α + β > 0): This formula played a key role in the proof of Proposition 1.4 and Theorem 1.5 in [15]. Thanks to the constraint 0 ≤ ρ ≤ 1, all the terms in this formula are non-negative (without the constraint, the first term will favor values of ρ larger than 1). Furthermore, in the regime ε 1, the first term will be bounded only for characteristic functions. This observation will be crucial when proving that the limiting density is a characteristic function (even though ρ ε may not be). We also note that the scaling of the following two terms is consistent with the scaling of the classical Modica-Mortola regularization of the perimeter functional.
Our analysis will also require a slight variation of this formula: we can write Here also we note that all the terms are non-negative when 0 ≤ ρ ≤ 1. This formula will in particular be crucial in proving the strong convergence of ρ ε and when deriving the pressure equation (1.15) (see Section 5).

A priori estimates
We now derive the a priori estimates that will be used to prove the convergence of ρ ε . We have: such that for all ε > 0 we have: Proof. We recall that where (since 0 ≤ ρ ≤ 1) −C ≤ ρ log ρ ≤ 0. The energy inequality (1.12) thus implies Using Proposition 1.4, we see that when ρ in = χ Ein ∈ BV (Ω; {0, 1}), we have J ε (ρ in ) ≤ C for some constant C independent on ε. We deduce The energy inequality also giveŝ and since ρ ε ≤ 1 (ii) follows immediately.

Strong Convergence of ρ ε and continuity equation.
The main result of this section is the following proposition which proves the first part of Theorem 1.6: Consider a sequence such that ε n → 0. The followings hold: (i) There exists a subsequence (still denoted ε n ) along which ρ εn (t) converges uniformly with respect to t in H −1 (Ω) to ρ(t) and E εn converges weakly in L 2 (Ω × (0, ∞)) to E. (ii) There exists v ∈ (L 2 (Ω × (0, ∞), dρ)) d such that E = ρv and the continuity equation (1.14) holds.
Note that (i) and (ii) are classical. The most important statement is thus (iii).
Proof. The a priori estimates of Lemma 3.1 implies (i). Furthermore, we can pass to the limit in (1.10) to getˆΩ . This is the continuity equation To complete the proof of (ii) and derive (1.14), we just need to show that E can be written in the form ρv. This is a classical argument which we recall here: First recall that the function +∞ otherwise is lower semi-continuous for the weak convergence of measure. Together with the uniform bound Θ(ρ εn , E εn ) =´T 0´Ω ρ εn |v εn | 2 ≤ C (see Lemma 3.1 (ii)), it implies that E is absolutely continuous with respect to ρ and that there exists v(t, ·) ∈ L 2 (dρ(t)) such that E = ρv. Inserting this in (4.1) yields (1.14).
The rest of the proof is devoted to (iii). The fact that we can get stronger convergence for the density is non trivial and is due to the fact that the energy J ε controls the BV norm of φ ε , which is close to ρ ε when ε 1. To see this, we introduce the function We then have F (σφ) = 2 min(σφ, 1 − σφ) and so 1 (as long as 0 ≤ ρ ≤ 1). This inequality, together with the formula (2.2) for J ε implies (note that in [15] a similar inequality was derived when ρ is a characteristic function. The computation above extends this important property of J ε to the case where 0 ≤ ρ ≤ 1 by using the formula (2.2)).
Next, we can write and we are going to show that the first two terms in the right hand side go to zero (uniformly in t): • Formula (2.1) and the energy bound (Lemma 3.1 (i)) imply • When ρ is a characteristic function, we have 2F (ρ) = ρ and so the second term in (4.4) vanishes.
It remains to show that ρ is a characteristic function. We note that given t > 0, we can extract a subsequence which converges a.e. in Ω and (4.6) implies that |{δ ≤ ρ(t) ≤ 1 − δ} = 0 for all δ > 0. We deduce that ρ(x, t) ∈ {0, 1} a.e. x ∈ Ω (for all t > 0). Finally, we can pass to the limit in (1.12) to get the energy inequality (1.19): The lower-semicontinuity of Θ allows us to pass to the limit in the dissipation and Proposition 1.4 gives (since F ε (ρ in ) = J ε (ρ in ) when ρ in is a characteristic function) The liminf property of Proposition C.1 and the strong convergence of ρ ε then yield: This completes the proof of Proposition 4.1.

5.
Convergence of the first variation 5.1. Convergence of the first variation and proof of Theorem 1.6. The only thing left to do to prove Theorem 1.6 is to pass to the limit in equation (1.11) (under the convergence assumption (1.20)) to derive (1.15). We recall (1.11) here: and we note that when ε 1, neither the term ε −1 ρ ε ∇φ ε nor the function p ε (or its gradient ∇p ε ) are bounded. As explained in the introduction, it is the modified pressure, defined by p ε + ε −1 1 2σ − φ ε ρ ε which is expected to converge to the pressure in the Hele-Shaw model with surface tension. We can also include the term µρ ε in this modified pressure and normalize it to have zero average. We thus set: After a straightforward computation, we can rewrite (1.11) as Passing to the limit in (5.3) is the objective of this section. The main challenge being the term which gives rise to the mean-curvature and the contact angle condition in the limit ε → 0. This term is related to the first variation of the energy J ε defined by (1.17). The key result of this section is Proposition 5.7 below, which states that the assumption on the convergence of the energy (1.20) implies the convergence of the first variation. Similar results have been proved for different energy functionals. In particular, a classical result of Reshetnyak [17] gives that if χ Eε converges to χ E strongly in L 1 , then the convergence of the perimeter implies the convergence of the first variation tô A similar result was proved by Luckhaus and Modica [13] for the Ginzburg-Landau functional and by Laux and Otto [12] with when G ε is the Gaussian kernel of variance ε 2 .
Crucially, both E 1 and E 2 are regularizations of the perimeter functional. In our case, we recall (see Theorem 1.5) that the energy functional J ε Γ-converges to the perimeter functional as well, together with a boundary term. Furthermore, it is easy to check that (5.4) is the first variation of J ε (ρ) for a perturbation defined by Remark 5.1. Note that (5.5) preserves the constraint ρ s ≤ 1 but not the condition´ρ s = 1 (unless div ξ = 0). Alternatively, we could consider the first variation of J ε (ρ) for a perturbation defined by which leads to the integral ε −1´Ω ρ ξ∇( 1 2σ − φ ε ) dx. The two integrals are the same when div ξ = 0, but this second integral does not, in general, converge when ε → 0. This is due to the fact that (5.6) does not preserve characteristic functions and the energy J ε blows up in that case for ε 1.
End of the proof of Theorem 1.6. We can now pass to the limit in (5.3) and thus complete the proof of Theorem 1.6. We recall that E εn converges weakly in L 2 and the convergence of the first term in the right hand side follows from Corollary 5.3, which we can use here since Proposition 4.1 gives the strong convergence of ρ εn in L 1 and we are assuming condition (1.20). We thus only need to explain why the pressure q ε converges.

Proof of Proposition As noted earlier, if we setφ
It is therefore enough to prove the result when σ = 1 and use the fact that to get the result when σ > 0.
The proof of Proposition 5.2 makes use of the following lemma, the proof of which is elementary and presented at the end of this section: Lemma 5.4. Given ρ(x) such that 0 ≤ ρ ≤ 1 and φ solution of (1.4) with σ = 1, we have the following formulas, for all ξ ∈ C 1 (Ω, R d ): (i) If α = 0 or β = 0 (Neumann or Dirichlet boundary condition for φ) then Idea of the proof of Proposition 5.2: While the proof of this Proposition 5.2 may appear long and technical, the idea is quite simple. We recall that when α = 0 (Neumann boundary condition), we have (see (2.2)): The boundedness of the energy implies that u ε and v ε are bounded in L 2 and in order to pass to the limit in the first two terms in (5.12), we need to show the convergence of´[u 2 ε + v 2 ε ]div ξ dx. Using the notations of Proposition 4.1, and in particular the function F such that F (φ) = 2 min(φ, 1− φ), we now write Since the limit ρ of ρ ε is a characteristic function and we know that φ ε → ρ, we have F (φ ε ) → F (ρ) = 1 2 ρ in L 1 (see the proof of Proposition 4.1), and so lim inf On the other hand, the assumption of convergence of the energy, (5.7), implies Together, these inequalities imply that there is equality in (5.14) when ε → 0, which means that To pass to the limit in the last term in (5.12), we note that the (asymptotic) equality in Young's inequality in (5.14) also implies that u ε − v ε → 0 in L 2 and so 2ε −1 |∇φ ε | 2 = 2v 2 ε ∼ u 2 ε + v 2 ε ∼ |∇F (φ ε )| which proves the convergence of 2ε −1 |∇φ ε | 2 . A simple (if somewhat technical) result (see Proposition A.2) then shows that the convergence of ε −1 |∇φ ε | 2 implies that of ε −1 ∇φ ε ⊗ ∇φ ε .
Additional care will be needed to take care of the boundary condition when α = 0, which is why we will first give the detailed proof for Neumann boundary conditions, then Dirichlet conditions (which requires extending φ ε to R d by 0) and finally general Robin boundary conditions (which combine the Neumann and Dirichlet case).
The difficult part of the proof is, of course, to establish the limit (5.9), and we will first give the proof in the simpler case of Neumann boundary conditions (α = 0). We have to pass to the limit in (5.12). As above, we use the function F such that F (φ) = 2 min(φ, 1 − φ) and note that Next, the strong convergence of ρ ε and (4.5) imply that F (φ ε ) converges to F (ρ) strongly in L 1 . Since ρ is a characteristic function and F (0) = 0, F (1) = 1/2, we deduce F (φ ε ) → F (ρ) = 1 2 ρ strongly in L 1 and so: On the other hand, the convergence assumption (5.7) implieŝ Inequality (5.15) thus implies: These facts allow us to pass to the limit in the first two terms of (5.12). Indeed, using first the definition of u ε and v ε , then the limit (5.16) and finally (5.17) (together with Proposition A.1), we can write: Furthermore, (5.16) and (5.18) yields: (Ω) which we use to pass to the limit in the term involving ∇φ ε ⊗ ∇φ ε . Indeed, we can write Using the fact that F (φ) ≥ 0, we can also write and using (5.17) and Proposition A.2 we deduce Using (5.12) together with (5.19) and (5.21) we get which gives (5.9) in the case of Neumann boundary condition α = 0 (and when σ = 1) For the case of Dirichlet conditions (β = 0), we proceed similarly, but we first extend φ ε to R d by setting it equal to 0 in R d \ Ω. Denotingφ ε this extension, we find (using inequality (5.15) and the Dirichlet boundary condition for φ ε ): The lower semicontinuity of the BV norm and assumption (5.7) give Indeed, when β = 0 (and σ = 1), we have We can thus proceed as before to show that These are the same convergences as (5.16)-(5.18), except for (5.17) which involves the extensionφ ε . We can now write: Similarly, we can show (as before) 2εˆΩ ∇φ ε ⊗ ∇φ ε : Dξ dx −→ 1 2ˆRd This gives the result when β = 0 (Dirichlet boundary conditions).
For the case of Robin conditions, we need to combine the two cases above. We first choose ϕ : Introducing the function G(t) = α β t 2 − F (t), we then write: The right hand side satisfies (with the same notations as above): so in order to proceed as before, we need to show that the lim inf of the left hand side is greater than 2J 0 (ρ). For this, we notice that the function thanks to our particular choice of ϕ. Going back to (5.22), we see that we just showed that and since this holds for any δ > 0, we get We can now conclude as in the previous cases: Using (5.23) and the assumption that lim ε→0 J ε (ρ ε ) = J 0 (ρ) to conclude that Furthermore, using (5.24), we also get lim sup We then write (using (5.25)) and using the same function ϕ as above, we decompose the integral in (5.28) as follows: Since the left hand side is independent of ϕ, we can take δ → 0 and use (5.28) to get It then only remains to show that which can be easily done by combining the argument above with how that term was handled in the case of Dirichlet boundary conditions.
In the case of general Robin boundary conditions (with in particular β = 0), we can write the last term as: and so which is (5.13).

JKO Scheme and convergence of the discrete time approximation
The main result of this paper can also be proved at the level of the discrete time approximation constructed in [11]. Such a result can be relevant to some numerical applications so we will state it here. First, we briefly recall the construction of the JKO scheme: As usual, P(Ω) denotes the set of probability measures on Ω and we denote (6.1) K = {ρ ∈ P(Ω), ρ(x) ≤ 1 a.e. in Ω} (in particular all ρ ∈ K are absolutely continuous with respect to the Lebesgue measure and we can identify the measure with its density). The set P(Ω) is equipped with the usual Wasserstein distance, defined by where Π(ρ 1 , ρ 2 ) denotes the set of all probability measures π ∈ P(Ω × Ω) with marginals ρ 1 and ρ 2 . The idea of the JKO scheme is to construct a time-discrete approximation of the solution by successive applications of a minimization problem: For a given initial data ρ in ∈ K, we fix a time step τ > 0 (destined to go to zero) and define the sequence ρ n by: The fact that this problem has a minimizer is proved in [11,Proposition 2.1]. Furthermore, if T n is the unique optimal transport map from ρ n to ρ n−1 (that is T n #ρ n = ρ n−1 and W 2 2 (ρ n , ρ n−1 ) = |x − T n (x)| 2 ρ n ), we define the velocity v n (x) = x − T n (x) τ and the pressure variable p n (x) such that ρ n v n = ε −1 ρ n ∇φ n − ∇p n , p n ∈ H 1 ρ n . with φ n solution of (1.4). The existence of p n is shown in [11,Proposition 2.6].