# Convergence to equilibrium for cross diffusion systems with nonlocal
interaction

Daniel Matthes  and  Christian Parsch  D. Matthes, C. Parsch – Zentrum
Mathematik/M8, Technische Universität München, Boltzmannstrasse 3, D-80538
Garching, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We study the existence and the rate of equilibration of weak solutions to a
two-component system of non-linear diffusion-aggregation equations, with small
cross diffusion effects. The aggregation term is assumed to be purely
attractive, and in the absence of cross diffusion, the flow is exponentially
contractive towards a compactly supported steady state. Our main result is
that for small cross diffusion, the system still converges, at a slightly
lower rate, to a deformed but still compactly supported steady state.

Our approach relies on the interpretation of the PDE system as a gradient flow
in a two-component Wasserstein metric. The energy consists of a uniformly
convex part responsible for self-diffusion and non-local aggregation, and a
totally non-convex part that generates cross diffusion; the latter is scaled
by a coupling parameter $\varepsilon>0$. The core idea of the proof is to
perform an $\varepsilon$-dependent modification of the convex/non-convex
splitting and establish a control on the non-convex terms by the convex ones.

This research was supported by the DFG Collaborative Research Center TRR 109,
“Discretization in Geometry and Dynamics”.

##  1\. Introduction

In this paper, we consider the following system of diffusion-aggregation
equations with (small) cross diffusion on ${\mathbb{R}}^{d}$:

$\begin{cases}\partial_{t}\rho_{1}=\mathrm{div}\mathopen{}\mathclose{{}\left(\rho_{1}\nabla[F_{1}^{\prime}(\rho_{1})+\varepsilon\partial_{r_{1}}h(\rho_{1},\rho_{2})+K*\rho_{2}]}\right),\\\
\partial_{t}\rho_{2}=\mathrm{div}\mathopen{}\mathclose{{}\left(\rho_{2}\nabla[F_{2}^{\prime}(\rho_{2})+\varepsilon\partial_{r_{2}}h(\rho_{1},\rho_{2})+K*\rho_{1}]}\right).\end{cases}$
(1)

Above, the pair $\bm{\rho}=(\rho_{1},\rho_{2})$ of time-dependent probability
densities is sought, $F_{1},\,F_{2}$ are given non-linearities for the
(self-)diffusions, $h$ induces cross diffusion of strength $\varepsilon\geq
0$, and $K:{\mathbb{R}}^{d}\to{\mathbb{R}}$ is an attractive interaction
kernel. We shall assume that $F_{1}$ and $F_{2}$ are convex and degenerate at
zero, $F_{1}^{\prime}(0)=F_{2}^{\prime}(0)=0$, that the derivatives of $h$ can
be controlled by those of $F_{1}$ and $F_{2}$ in a way precised further below,
and that $K$’s convexity is uniform, $\nabla^{2}K\geq\lambda\,\mathbf{1}$ with
some $\lambda>0$. The system is augmented with initial conditions,

$\displaystyle\rho_{1}(0,\cdot)=\rho_{1}^{0},\quad\rho_{2}(0,\cdot)=\rho_{2}^{0},$
(2)

which we assume to have finite energy, see below. The PDE system (1) falls
into a general class of macroscopic models for multi-species many-particle
systems in which long-range attraction forces compete with local diffusion —
the former tend to accumulate particles, the latter prevent overcrowding. Such
competition is fundamental in various models for interacting particles,
including chemotaxis [6, 14], granular media [7, 12], social dynamics [3, 5],
semi-conductors [15] and various further applications. Here, we shall not
focus on specific choices of the parameters $F_{1}$, $F_{2}$, $K$ and $h$
related to any of these models, but consider (1) under quite general
hypotheses.

###  1.1. Motivation: persistence of equilibration

Our goal is to show at the example of the system (1) the persistence of the
equilibration behaviour in contractive aggregation-(self-)diffusion flows
under the influence of cross diffusion. Specifically, we perform a rigorous
analysis of existence and long-time asymptotics of weak solutions for
sufficiently small $\varepsilon>0$, and prove that the solutions are global in
time and converge to the unique equilibrium at almost the same exponential
rate as in the unperturbed situation $\varepsilon=0$. This result is beyond
perturbation theory not only because of its global character, but also because
the perturbation is singular as we shall explain next.

The point of departure for our analysis is the gradient flow formulation of
(1): the respective energy functional is

$\mathcal{E}_{\varepsilon}(\bm{\rho}):=\int_{\mathbb{R}^{d}}[F_{1}(\rho_{1})+F_{2}(\rho_{2})+\varepsilon
h(\rho_{1},\rho_{2})+\rho_{1}K*\rho_{2}]\,\mathrm{d}x,$ (3)

the metric is a combined $L^{2}$-Wasserstein distance for $\rho_{1}$ and
$\rho_{2}$. In the absence of cross diffusion, i.e., for $\varepsilon=0$, that
functional is uniformly geodesically convex, and — thanks to the general
theory of metric gradient flows [2] — the generated gradient flow (1) is
uniformly exponentially contractive. In particular, there is a unique
stationary solution $\bar{\bm{\rho}}$ to (1), and that is approached at an
exponential rate $\lambda>0$, see e.g. [2, Chapter 11]. Now, the perturbation
induced by cross diffusion is singular in the sense that for any
$\varepsilon>0$, the energy $\mathcal{E}_{\varepsilon}$ is not even semi-
convex anymore, and consequently, the gradient flow (1) is not contractive.
That devastating effect of cross diffusion has been proven in a related
situation in [4, Prop. 2.4].

Despite the failure of the abstract machinery, we shall prove that for
sufficiently small $\varepsilon>0$, the equilibration behaviour persists in
the sense that there is a unique stationary solution
$\bar{\bm{\rho}}_{\varepsilon}$ to (1), which is a deformation of
$\bar{\bm{\rho}}$ above, and that arbitrary solutions converge to
$\bar{\bm{\rho}}_{\varepsilon}$ at an exponential rate
$\lambda_{\varepsilon}\geq\lambda-C_{0}\varepsilon$. Our result thus
contributes to a line of studies on the exponential equilibration in non-
convex metric gradient flows, see [1, 4, 16, 13, 17, 25, 22, 23]. Differently
from these previous works, the confining force leading to equilibration in the
unperturbed system (1) is not induced by external potentials, but is self-
induced via the non-local interaction of the species. This is very natural
from a modeling point of view, but induces additional difficulties both on the
technical and on the conceptual level. Generally speaking, the equilibrating
force of non-local interactions is more fragile than that of fixed external
potentials, and requires a more refined analysis. Most prominently, while the
integral $\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x$ is a
uniformly geodesically convex functional of the pair
$\bm{\rho}=(\rho_{1},\rho_{2})$ in the considered metric [24], this functional
is _not_ convex in the sense of affine interpolation, which leads to non-
convexity of the entire energy $\mathcal{E}_{\varepsilon}$. Thus we are, for
instance, not able to exclude the emergence of additional critical points in
the perturbed energy $\mathcal{E}_{\varepsilon}$ a priori by standard method
from the calculus of variations, but need to obtain the uniqueness of the
stationary solution $\bar{\bm{\rho}}_{\varepsilon}$ a posteriori from the
global equilibration estimates.

There has been a great recent interest in understanding the long-time
behaviour of multi-species cross-diffusion systems, and particularly those
with non-local effects. The focus has been on proving qualitative properties
of the steady states, in particular segregation phenomena, see [8, 9, 10, 11]
and references therein for analytical and numerical studies. It appears,
however, that there are no rigorous results on the global exponential
stability of these stationary solutions, not even in the non-segregated case
that we study here.

###  1.2. Main result

Below, we summarize the main results of this paper in Theorem 1.1. For a
concise statement, introduce the space ${\mathcal{P}_{2}}({\mathbb{R}^{d}})$
of probability measures on ${\mathbb{R}^{d}}$ with finite second moment, and
the following subspace
${\mathfrak{X}}_{0}\subset{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ of pairs of
measures whose combined center of mass is at the origin:

$\displaystyle{\mathfrak{X}}_{0}:=\mathopen{}\mathclose{{}\left\\{\bm{\rho}=(\rho_{1},\rho_{2})\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}:\int_{\mathbb{R}^{d}}x\,\mathopen{}\mathclose{{}\left(\rho_{1}+\rho_{2}}\right)\,\mathrm{d}x=0}\right\\}.$
(4)

We prove:

######  Theorem 1.1.

Assuming that $F_{1},F_{2},h$ and $K$ satisfy Hypotheses 2.1 and 2.2 further
below, there exists some $\bar{\varepsilon}>0$ such that for every
$\varepsilon\in(0,\bar{\varepsilon}]$, the following is true.

  * •

The energy functional $\mathcal{E}_{\varepsilon}$ possesses a unique
$\varepsilon$-dependent global minimizer
$\bar{\bm{\rho}}=(\bar{\rho}_{1},\bar{\rho}_{2})$ in ${\mathfrak{X}}_{0}$. The
densities $\bar{\rho}_{1}$ and $\bar{\rho}_{2}$ are continuous in
${\mathbb{R}^{d}}$, radially symmetric about the origin, and compactly
supported. They satisfy, with suitable constants $C_{1},C_{2}>0$, the Euler-
Lagrange system

$\begin{split}F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})&=\mathopen{}\mathclose{{}\left(C_{1}-K*\bar{\rho}_{2}}\right)_{+}\\\
F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\partial_{r_{2}}h(\bar{\rho}_{1},\bar{\rho}_{2})&=\mathopen{}\mathclose{{}\left(C_{2}-K*\bar{\rho}_{1}}\right)_{+}\,.\end{split}$
(5)

  * •

Given initial data $\bm{\rho}^{0}\in{\mathfrak{X}}_{0}$ of finite energy, i.e.
$\mathcal{E}_{\varepsilon}(\bm{\rho}^{0})<+\infty$, there exists a weakly
continuous curve $\bm{\rho}:[0,+\infty)\to{\mathfrak{X}}_{0}$ with
$\bm{\rho}(0)=\bm{\rho}^{0}$ that satisfies the system (1) in the sense that
for every compactly supported test function $\zeta\in
C^{\infty}_{c}({\mathbb{R}}_{>0}\times{\mathbb{R}^{d}})$,

$\begin{split}\int_{0}^{\infty}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\rho_{1}\,\partial_{t}\zeta-\rho_{1}\,\nabla\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1})+\varepsilon\partial_{r_{1}}h(\rho_{1},\rho_{2})+K\ast\rho_{2}}\right]\cdot\nabla\zeta}\right]\,\mathrm{d}x\mathrm{d}t&=0,\\\
\int_{0}^{\infty}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\rho_{2}\,\partial_{t}\zeta-\rho_{2}\,\nabla\mathopen{}\mathclose{{}\left[F_{2}^{\prime}(\rho_{2})+\varepsilon\partial_{r_{2}}h(\rho_{1},\rho_{2})+K\ast\rho_{1}}\right]\cdot\nabla\zeta}\right]\,\mathrm{d}x\mathrm{d}t&=0.\end{split}$

  * •

There is a constant $C_{0}>0$, independent of
$\varepsilon\in(0,\bar{\varepsilon}]$, such that the energy decays
exponentially fast along this solution $\bm{\rho}$ with rate $\lambda-
C_{0}\varepsilon$:

$\displaystyle\mathcal{E}_{\varepsilon}(\bm{\rho}(t))-\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})\leq
C_{1}\,\mathopen{}\mathclose{{}\left(\mathcal{E}_{\varepsilon}(\bm{\rho}^{0})-\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})}\right)\,e^{-2(\lambda-
C_{0}\varepsilon)t}$ (6)

at every $t\geq 0$ for some constant $C_{1}$ independent of $\varepsilon$ and
$\bm{\rho}^{0}$. Further, $\bm{\rho}$ converges to the unique global minimizer
$\bar{\bm{\rho}}$ exponentially in $L^{1}$ with the same rate.

###  1.3. Strategy of proof

We give cartoon of our approach to proving the persistence of long-time
asymptotics in ODE language. Assume that two potentials
$E_{0},\,H:{\mathbb{R}}\to{\mathbb{R}}$ are given, where $E_{0}\in
C^{2}({\mathbb{R}})$ is uniformly convex of modulus $\lambda>0$, i.e.,
$E_{0}^{\prime\prime}\geq\lambda$, and $H\in C^{1}({\mathbb{R}})$ is bounded
below but _not_ convex. For $\varepsilon\geq 0$, define the energy
$E_{\varepsilon}:=E_{0}+\varepsilon H$, which is in general not convex for any
$\varepsilon>0$. We wish to study the long-time asymptotics of the associated
gradient flow,

$\displaystyle\dot{x}=-E_{\varepsilon}^{\prime}(x)=-E_{0}^{\prime}(x)-\varepsilon
H^{\prime}(x).$ (7)

Thanks to the lower bound on $H$, there exists a global minimizer
$\bar{x}_{\varepsilon}$ of $E_{\varepsilon}$ for every $\varepsilon\geq 0$; a
priori, $\bar{x}_{\varepsilon}$ need not be unique. Now introduce the two
auxiliary functions

$\displaystyle
L_{\varepsilon}(x)=E_{0}(x)+\varepsilon(x-\bar{x}_{\varepsilon})H^{\prime}(\bar{x}_{\varepsilon}),\quad
N_{\varepsilon}(x)=H(x)-(x-\bar{x}_{\varepsilon})H^{\prime}(\bar{x}_{\varepsilon}),$

and note that $\bar{x}_{\varepsilon}$ is the global minimum of the
$\lambda$-convex function $L_{\varepsilon}$, and is a critical point of the
non-convex function $N_{\varepsilon}$. Finally, we add our core hypothesis on
the relation between $H$ and $E_{0}$: for $0<\varepsilon<1$, the following
inequality holds for all $x\in{\mathbb{R}}$:

$\displaystyle|N_{\varepsilon}^{\prime}(x)|\leq|L_{\varepsilon}^{\prime}(x)|.$
(8)

With that at hand, $L_{\varepsilon}$ becomes a Lyapunov functional: for any
solution $x$ to (7),

$\displaystyle-\frac{\,\mathrm{d}}{\mathrm{d}t}\big{[}L_{\varepsilon}(x)-L_{\varepsilon}(x_{\varepsilon})\big{]}=-L_{\varepsilon}^{\prime}(x)\dot{x}$
$\displaystyle=L_{\varepsilon}^{\prime}(x)^{2}+\varepsilon
L_{\varepsilon}^{\prime}(x)N_{\varepsilon}^{\prime}(x)$
$\displaystyle\geq(1-\varepsilon)L_{\varepsilon}^{\prime}(x)^{2}\geq
2\lambda(1-\varepsilon)\big{[}L_{\varepsilon}(x)-L_{\varepsilon}(x_{\varepsilon})\big{]},$

the last inequality being a consequence of $L_{\varepsilon}$’s
$\lambda$-convexity. It follows that $L_{\varepsilon}(x)$ tends to its minimal
value at the (diminished) exponential rate $2\lambda(1-\varepsilon)$. It then
further follows easily that $E_{\varepsilon}(x)$ decays at the same
exponential rate, and also that the solution $x$ itself goes to
$\bar{x}_{\varepsilon}$ at that rate. And in particular, the chosen minimizer
$\bar{x}_{\varepsilon}$ turns out to be the unique one.

For (1), the roles of $E_{\varepsilon}$ and of $H$ are played by
$\mathcal{E}_{\varepsilon}$, which is $\lambda$-uniformly geodesically convex
for $\varepsilon=0$, and by the integral of $h(\rho_{1},\rho_{2})$, that has
no geodesic convexity, respectively. The auxiliary functionals
$L_{\varepsilon}$ and $N_{\varepsilon}$ correspond to, respectively,

$\displaystyle\mathcal{L}_{\varepsilon}(\bm{\rho})$
$\displaystyle:=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[F_{1}(\rho_{1})+F_{2}(\rho_{2})+\rho_{1}K\ast\rho_{2}+\varepsilon(\rho_{1}V_{1}+\rho_{2}V_{2})}\right]\,\mathrm{d}x$
(9) $\displaystyle\mathcal{N}_{\varepsilon}(\bm{\rho})$
$\displaystyle:=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[h(\rho_{1},\rho_{2})-(\rho_{1}V_{1}+\rho_{2}V_{2})}\right]\,\mathrm{d}x,$
(10)

where the auxiliary potentials $V_{1}$ and $V_{2}$ are chosen to make a global
minimizer $\bar{\bm{\rho}}$ of $\mathcal{E}_{\varepsilon}$, which we obtain
via the direct methods, the unique global minimizer of
$\mathcal{L}_{\varepsilon}$, i.e.,

$\displaystyle
V_{1}:=\partial_{r_{1}}h\big{(}\bar{\rho}_{1},\bar{\rho}_{2}),\quad
V_{2}:=\partial_{r_{2}}h\big{(}\bar{\rho}_{1},\bar{\rho}_{2}).$ (11)

In contrast to $L_{\varepsilon}$ in the ODE cartoon above, the modulus of
convexity of $\mathcal{L}_{\varepsilon}$ is not precisely $\lambda$, but
$\lambda-\varepsilon K_{0}$, where $K_{0}$ is an upper bound on the second
derivatives of $V_{1}$, $V_{2}$, see Corollary 5.2. Finally, the core relation
(8) is a functional inequality, see Proposition 5.7.

###  1.4. Outline

We start in Section 2 by specifying the assumptions on the nonlinearities
$F_{j}$, the coupling term $h$ and the interaction potential $K$. A generic
example of functions which satisfy all assumptions is given in the Appendix.
In Section 3, we prove existence of a global minimizer $\bar{\bm{\rho}}$ of
the energy $\mathcal{E}_{\varepsilon}$ over ${\mathfrak{X}}_{0}$ by means of
the direct methods from the calculus of variations, derive the Euler-Lagrange
system, and obtain additional regularity properties. In Section 4, we analyze
the time-discrete minimizing movement scheme for the energy
$\mathcal{E}_{\varepsilon}$, which is then used to construct a weak solution
to (1)–(2). Finally, in Section 5, we analyze the long-time behaviour of this
weak solution and prove exponential convergence to the unique equilibrium for
sufficiently small $\varepsilon>0$.

##  2\. Preliminaries and hypotheses

###  2.1. Notations and definitions

Recall that ${\mathcal{P}_{2}}({\mathbb{R}^{d}})$ is the space of probability
measures with finite second moment on ${\mathbb{R}^{d}}$; we consider
${\mathcal{P}_{2}}({\mathbb{R}^{d}})$ with the topology induced by narrow
convergence. By abuse of notation, we do not distinguish between an absolutely
continuous element $\rho\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$ and its
(Lebesgue-)density $\rho\in L^{1}({\mathbb{R}^{d}})$. The $L^{2}$-Wasserstein
distance ${\mathbf{W}}$ is a metric on ${\mathcal{P}_{2}}({\mathbb{R}^{d}})$,
and makes ${\mathcal{P}_{2}}({\mathbb{R}^{d}})$ a complete metric space.
Convergence in ${\mathbf{W}}$ is equivalent to weak convergence of measures
and convergence of the second moment. For an introduction to Wasserstein
metrics and the related gradient flows, we refer to [20, 21]. From the
Wasserstein metric on ${\mathcal{P}_{2}}({\mathbb{R}^{d}})$, we directly
obtain a metric $\mathrm{d}$ on the space
${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ of pairs
$\bm{\rho}=(\rho_{1},\rho_{2})$ of measures,

$\displaystyle\mathrm{d}\mathopen{}\mathclose{{}\left(\bm{\rho},\bm{\eta}}\right)^{2}:={\mathbf{W}}(\rho_{1},\eta_{1})^{2}+{\mathbf{W}}(\rho_{2},\eta_{2})^{2}.$

The energy functional
$\mathcal{E}_{\varepsilon}:{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}\to{\mathbb{R}}\cup\\{+\infty\\}$
formally introduced in (3) is defined on pairs
$\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ with the interpretation
that $\mathcal{E}_{\varepsilon}(\bm{\rho})=+\infty$ unless $\rho_{1}$ and
$\rho_{2}$ are both absolutely continuous, and the integral expression is in
$L^{1}({\mathbb{R}^{d}})$; we shall see below that this definition makes
$\mathcal{E}_{\varepsilon}$ weakly lower semi-continuous, at least for
sufficiently small $\varepsilon>0$. For brevity, introduce
$\mathcal{F}_{\varepsilon}:{\mathbb{R}}_{\geq 0}^{2}\to{\mathbb{R}}_{\geq
0}^{2}$ by
$\mathcal{F}_{\varepsilon}({\mathbf{r}}):=F_{1}(r_{1})+F_{2}(r_{2})+\varepsilon
h({\mathbf{r}})$, which allows to rewrite $\mathcal{E}_{\varepsilon}$ as

$\displaystyle\mathcal{E}_{\varepsilon}(\bm{\rho})=\int_{\mathbb{R}^{d}}[\mathcal{F}_{\varepsilon}(\bm{\rho})+\rho_{1}K\ast\rho_{2}]\,\mathrm{d}x.$

The functional
$\bm{\rho}\mapsto\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x$ will
be called the interaction energy.

Next, recall the definition (4) of the set
${\mathfrak{X}}_{0}\subset{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ of pairs
$\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ with combined center of
mass at zero,

$\displaystyle{\mathfrak{m}}_{1}[\rho_{1}]+{\mathfrak{m}}_{1}[\rho_{2}]=0\quad\text{where}\quad{\mathfrak{m}}_{1}[\rho_{j}]:=\int_{\mathbb{R}^{d}}x\,\rho_{j}\,\mathrm{d}x.$

Since $\mathcal{E}_{\varepsilon}$ is invariant under translations, i.e.,

$\displaystyle\mathcal{E}_{\varepsilon}(\bm{\rho})=\mathcal{E}_{\varepsilon}(\sigma_{v}\bm{\rho})\quad\text{with}\
\sigma_{v}\bm{\rho}=\mathopen{}\mathclose{{}\left(\rho_{1}(\cdot+v),\rho_{2}(\cdot+v)}\right)\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$

for any translation $v\in{\mathbb{R}}^{d}$ and any
$\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$, and since (1) conserves the
combined center of mass, there is no loss of generality to restrict attention
to initial data (2) and solutions to (1) belonging to ${\mathfrak{X}}_{0}$.

Finally, we introduce the notation $j^{\prime}=3-j$ for $j\in\\{1,2\\}$, so
that $\rho_{1^{\prime}}=\rho_{2}$ and $\rho_{2^{\prime}}=\rho_{1}$.

###  2.2. Summary of hypotheses

Throughout this paper, we consider $F_{1},F_{2},h$ and $K$ fullfilling the
following assumptions:

######  Hypothesis 2.1.

Assumptions on the interaction potential
$K:{\mathbb{R}^{d}}\to{\mathbb{R}}_{\geq 0}$:

  * •

$K$ is radially symmetric about the origin and it holds $K(0)=0$.

  * •

$K$ is $\lambda$-convex with some modulus $\lambda>0$.

  * •

It holds $K\in C^{2}({\mathbb{R}^{d}})$, and there exists a constant
$C_{K}<+\infty$ with $\|\nabla^{2}K(z)\|\leq C_{K}$ for every
$z\in{\mathbb{R}^{d}}$.

Assumptions on the nonlinearities $F_{1},F_{2}:{\mathbb{R}}_{\geq
0}\to{\mathbb{R}}_{\geq 0}$ and coupling term $h:{\mathbb{R}}_{\geq
0}^{2}\to{\mathbb{R}}$:

  * •

$F_{1}$ and $F_{2}$ are $C^{2}$ and strictly convex with
$F_{j}^{\prime\prime}(r)>0$ for every $r>0$ and $j=1,2$.

  * •

It holds $F_{1}(0)=F_{2}(0)=F_{1}^{\prime}(0)=F_{2}^{\prime}(0)=0$.

  * •

Their derivatives $F_{1}^{\prime}$ and $F_{2}^{\prime}$ grow at least linearly
for large inputs, but not faster than $F_{1}$ and $F_{2}$, respectively. More
precisely, it holds for $j=1,2$:

$\liminf_{r\to+\infty}F_{j}^{\prime\prime}(r)>0,\quad\limsup_{r\to+\infty}\frac{F_{j}^{\prime}(r)}{F_{j}(r)}<+\infty$
(12)

  * •

For small inputs, their second derivatives behave similarly to powers of $r$.
Explicitly, for $j=1,2$ there exist constants $m_{j},M_{j}>0$, exponents
$\beta_{j}\geq 0$ and an $r_{0}>0$ such that for all $0\leq r\leq r_{0}$, it
holds

$m_{j}r^{\beta_{j}}\leq F_{j}^{\prime\prime}(r)\leq M_{j}r^{\beta_{j}}.$ (13)

  * •

They satisfy the McCann-condition for displacement convexity: For all $r\geq
0$,

$F_{j}(r)-rF_{j}^{\prime}(r)+r^{2}F_{j}^{\prime\prime}(r)\geq 0.$ (14)

  * •

The coupling function $h$ is $C^{2}$ and satisfies
$h(r_{1},0)=h(0,r_{2})=\partial_{r_{j}}h(r_{1},0)=\partial_{r_{j}}h(0,r_{2})=0$
for $j=1,2$ and every $r_{1},r_{2}\geq 0$.

It directly follows from the convexity of $F_{j}$ that the functions
$F_{j}^{\prime}:{\mathbb{R}}_{\geq 0}\to{\mathbb{R}}_{\geq 0}$ are strictly
increasing and thus invertible, since they vanish at $0$ and grow at least
linearly. We can thus introduce the following notations: For $j=1,2$ and
$\bm{u}=(u_{1},u_{2})\in{\mathbb{R}}_{\geq 0}^{2}$, we define

$\theta_{j}(\bm{u}):=\partial_{r_{j}}h\mathopen{}\mathclose{{}\left((F_{1}^{\prime})^{-1}(u_{1}),(F_{2}^{\prime})^{-1}(u_{2})}\right).$

In other words, $\theta_{j}:{\mathbb{R}}_{\geq 0}^{2}\to{\mathbb{R}}$ is such
that
$\theta_{j}(F_{1}^{\prime}(r_{1}),F_{2}^{\prime}(r_{2}))=\partial_{r_{j}}h({\mathbf{r}})$
for all ${\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$. Further, we define for
$i,j=1,2$ the functions $\theta_{j,i}:{\mathbb{R}}_{\geq
0}^{2}\to{\mathbb{R}}$ by

$\theta_{j,i}(\bm{u}):=\partial_{u_{i}}\theta_{j}(\bm{u})=\frac{\partial_{r_{i}}\partial_{r_{j}}h({\mathbf{r}})}{F_{i}^{\prime\prime}(r_{i})},$

where
${\mathbf{r}}=(r_{1},r_{2})=((F_{1}^{\prime})^{-1}(u_{1}),(F_{2}^{\prime})^{-1}(u_{2}))$.
In addition to Hypothesis 2.1, we impose the following conditions on the
$\theta_{j,i}$:

######  Hypothesis 2.2.

We assume that the $F_{j}$ and $h$ are such that

  * •

The $\theta_{j,i}$ are locally Lipschitz in ${\mathbb{R}}_{\geq 0}^{2}$.

  * •

There are constants $\kappa_{j,i}>0$ such that

$\mathopen{}\mathclose{{}\left|\theta_{j,i}(\bm{u})}\right|\leq\kappa_{j,i}\,\min\mathopen{}\mathclose{{}\left\\{1,u_{1},u_{2},\sqrt{\frac{(F_{i}^{\prime})^{-1}(u_{i})}{(F_{j}^{\prime})^{-1}(u_{j})}}}\right\\}$
(15)

for all $\bm{u}\in{\mathbb{R}}_{\geq 0}^{2}$, which can also be stated in
terms of $F_{1},F_{2}$ and $h$ as

$\mathopen{}\mathclose{{}\left|\partial_{r_{i}}\partial_{r_{j}}h({\mathbf{r}})}\right|\leq\kappa_{j,i}\,F_{i}^{\prime\prime}(r_{i})\,\min\mathopen{}\mathclose{{}\left\\{1,F_{1}^{\prime}(r_{1}),F_{2}^{\prime}(r_{2}),\sqrt{\frac{r_{i}}{r_{j}}}}\right\\}$
(16)

for all ${\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$.

######  Example 2.3.

A simple example for an interaction potential $K$ that satisfies Hypothesis
2.1 is given by the quadratic potential $K=\frac{\lambda}{2}|\cdot|^{2}$ with
some $\lambda>0$. An example for $F_{1},F_{2}$ and $h$ that fulfill Hypotheses
2.1 and 2.2 is given by Proposition A.1 in the appendix. Note that the
assumptions on $K$ are independent from those on $F_{1},F_{2}$ and $h$.

###  2.3. Basic implications of the hypotheses

We assume in all of the following that $F_{1},F_{2},h$ and $K$ are such that
Hypotheses 2.1 and 2.2 are satisfied. We start by investigating the
assumptions more closely and proving some general consequences of them, which
will be useful in our analysis. We first prove some more properties of $F_{j}$
and $h$.

######  Proposition 2.4.

The following hold:

  * •

There is $\varepsilon_{0}>0$ such that the function

${\mathbf{r}}\mapsto\mathcal{F}_{2\varepsilon}({\mathbf{r}})=F_{1}(r_{1})+F_{2}(r_{2})+2\varepsilon
h({\mathbf{r}})$ (17)

is convex on ${\mathbb{R}}_{\geq 0}^{2}$ for every
$0<\varepsilon\leq\varepsilon_{0}$, strictly if $\varepsilon<\varepsilon_{0}$.

  * •

For all ${\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$ and
$\varepsilon\leq\varepsilon_{0}$, it holds with the constants $\kappa_{j,i}$
from (15):

$\frac{1}{2}F_{1}(r_{1})+\frac{1}{2}F_{2}(r_{2})\leq\mathcal{F}_{\varepsilon}({\mathbf{r}})\leq\mathopen{}\mathclose{{}\left(1+\frac{\varepsilon}{2}\kappa_{1,1}}\right)F_{1}(r_{1})+\mathopen{}\mathclose{{}\left(1+\frac{\varepsilon}{2}\kappa_{2,2}}\right)F_{2}(r_{2}).$
(18)

  * •

With the constants $r_{0}$ and $m_{j},M_{j},\beta_{j}$ from (13), it holds for
$0\leq r\leq r_{0}$:

$\displaystyle\frac{m_{j}}{\beta_{j}+1}r^{\beta_{j}+1}$ $\displaystyle\leq
F_{j}^{\prime}(r)\leq\frac{M_{j}}{\beta_{j}+1}r^{\beta_{j}+1},$ (19)
$\displaystyle\frac{m_{j}}{(\beta_{j}+1)(\beta_{j}+2)}r^{\beta_{j}+2}$
$\displaystyle\leq
F_{j}(r)\leq\frac{M_{j}}{(\beta_{j}+1)(\beta_{j}+2)}r^{\beta_{j}+2}.$ (20)

  * •

For $j=1,2$, there are constants $\alpha_{j}$ such that for all $r\geq 0$, it
holds

$F_{j}^{\prime}(r)\leq\alpha_{j}\,\mathopen{}\mathclose{{}\left(\min\\{1,r\\}+F_{j}(r)}\right).$
(21)

  * •

For any bound $H>0$, there is a constant $\beta_{H}>0$ such that for all
$\bar{\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$ with
$\bar{r}_{1},\bar{r}_{2}\leq H$, all ${\mathbf{r}}\in{\mathbb{R}}_{\geq
0}^{2}$ and any $i,j\in\\{1,2\\}$, the estimate

$\max\\{r_{1},r_{2}\\}\,\mathopen{}\mathclose{{}\left|\theta_{j,i}(\bm{u})-\theta_{j,i}(\bar{\bm{u}})}\right|^{2}\leq\beta_{H}\,\mathopen{}\mathclose{{}\left(\mathrm{d}_{F_{1}}(r_{1}|\bar{r}_{1})+\mathrm{d}_{F_{2}}(r_{2}|\bar{r}_{2})}\right)$
(22)

holds, where
$\bm{u}=\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(r_{1}),F_{2}^{\prime}(r_{2})}\right)$,
$\bar{\bm{u}}=\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\bar{r}_{1}),F_{2}^{\prime}(\bar{r}_{2})}\right)$
and $\mathrm{d}_{F}(\cdot|\cdot)$ is the Bregman divergence of a convex
function $F$, given by

$\mathrm{d}_{F}(r|\bar{r}):=F(r)-F(\bar{r})-F^{\prime}(\bar{r})(r-\bar{r}).$

###### Proof.

We prove convexity of (17) by observing that for every
${\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$, it holds

$\displaystyle\nabla_{\mathbf{r}}^{2}\mathcal{F}_{2\varepsilon}({\mathbf{r}})=\begin{pmatrix}F_{1}^{\prime\prime}(r_{1})+2\varepsilon\partial_{r_{1}}\partial_{r_{1}}h({\mathbf{r}})&2\varepsilon\partial_{r_{2}}\partial_{r_{1}}h({\mathbf{r}})\\\
2\varepsilon\partial_{r_{1}}\partial_{r_{2}}h({\mathbf{r}})&F_{2}^{\prime\prime}(r_{2})+2\varepsilon\partial_{r_{2}}\partial_{r_{2}}h({\mathbf{r}})\end{pmatrix}.$

We show that for some $\varepsilon_{0}>0$, this matrix is positive semi-
definite at every ${\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$. Indeed, it holds
for any $\bm{v}=(v_{1},v_{2})\in{\mathbb{R}}^{2}$ by assumption (16) and
$2v_{1}v_{2}\geq-v_{1}^{2}-v_{2}^{2}$:

$\displaystyle\bm{v}^{T}\nabla_{\mathbf{r}}^{2}\mathcal{F}_{2\varepsilon}({\mathbf{r}})\bm{v}$
$\displaystyle=\mathopen{}\mathclose{{}\left(F_{1}^{\prime\prime}(r_{1})+2\varepsilon\partial_{r_{1}}\partial_{r_{1}}h({\mathbf{r}})}\right)v_{1}^{2}+\mathopen{}\mathclose{{}\left(F_{2}^{\prime\prime}(r_{2})+\varepsilon\partial_{r_{2}}\partial_{r_{2}}h({\mathbf{r}})}\right)v_{2}^{2}+4\varepsilon\partial_{r_{1}}\partial_{r_{2}}h({\mathbf{r}})v_{1}v_{2}$
$\displaystyle\geq\mathopen{}\mathclose{{}\left(1-2\varepsilon(\kappa_{1,1}+\kappa_{2,1})}\right)\,F_{1}^{\prime\prime}(r_{1})v_{1}^{2}+\mathopen{}\mathclose{{}\left(1-2\varepsilon(\kappa_{1,2}+\kappa_{2,2})}\right)\,F_{2}^{\prime\prime}(r_{2})v_{2}^{2}\geq
0$

for every sufficiently small $\varepsilon>0$ by non-negativity of
$F_{j}^{\prime\prime}({\mathbf{r}})$. Hence there exists $\varepsilon_{0}>0$
such that $F_{2\varepsilon_{0}}$ is convex on ${\mathbb{R}}_{\geq 0}^{2}$. For
every $\varepsilon<\varepsilon_{0}$, it follows that

$\displaystyle\mathcal{F}_{2\varepsilon}({\mathbf{r}})=\frac{\varepsilon}{\varepsilon_{0}}F_{2\varepsilon_{0}}({\mathbf{r}})+\mathopen{}\mathclose{{}\left(1-\frac{\varepsilon}{\varepsilon_{0}}}\right)\mathopen{}\mathclose{{}\left(F_{1}(r_{1})+F_{2}(r_{2})}\right)$

is strictly convex on ${\mathbb{R}}_{\geq 0}^{2}$ by strict convexity of
$F_{1}$ and $F_{2}$, proving the first claim. To show the first inequality in
(18), notice that $\mathcal{F}_{2\varepsilon}$ is non-negative in
${\mathbb{R}}_{\geq 0}^{2}$, since it is convex and its derivative vanishes at
0. Hence

$\displaystyle
2\mathcal{F}_{\varepsilon}({\mathbf{r}})=\mathcal{F}_{2\varepsilon}({\mathbf{r}})+F_{1}(r_{1})+F_{2}(r_{2})\geq
F_{1}(r_{1})+F_{2}(r_{2}).$

The second inequality in (18) will follow from the fact that
$|\partial^{2}_{r_{i}r_{i}}h({\mathbf{r}})|\leq\kappa_{i,i}F_{i}^{\prime\prime}(r_{i})$
for all ${\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$ and $i=1,2$ by assumption
(16). Integrated twice from $0$ to $r_{i}$ and using the assumption that $h$
and its first derivatives vanish on $\partial{\mathbb{R}}_{\geq 0}^{2}$, this
yields $|h({\mathbf{r}})|\leq\kappa_{i,i}F_{i}(r_{i})$ for all ${\mathbf{r}}$,
implying (18) by definition of $\mathcal{F}_{\varepsilon}$.

Inequalities (19) and (20) follow directly by integrating (13) from 0 to $r$
and using that $F_{j}(0)=F_{j}^{\prime}(0)=0$. The proof of (21) can be done
for $r\leq r_{0}$ and $r>r_{0}$ separately: For $r\leq r_{0}$, the inequality
follows from (19) since $\beta_{j}\geq 0$ and thus $r^{\beta_{j}+1}$ decays at
least linearly for small $r$. The case $r>r_{0}$ follows from the growth
assumption (12) and the fact that the right-hand side of (20) is bounded away
from 0 for $r>r_{0}$.

It remains to prove (22). Fix $H>0$ and let
${\mathbf{r}},\bar{\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2}$ with
$\bar{r}_{1},\bar{r}_{2}\leq H$. Define

$\displaystyle m_{H}:=\min_{s\geq H,\,j=1,2}F_{j}^{\prime\prime}(s),\quad
M_{H}:=\max_{s\in[0,3H],\,j=1,2}F_{j}^{\prime\prime}(s).$

It holds $0<m_{H}\leq M_{H}<+\infty$, by $F_{j}^{\prime\prime}(r)>0$ for $r>0$
and assumption (12). By symmetry of the assumptions and the claim, we can
assume w.l.o.g. $r_{1}\geq r_{2}$. We split the proof into two cases:
$r_{1}\geq 3H$ and $r_{1}\leq 3H$.

Consider first $r_{1}\geq 3H$. We have $\bar{r}_{1}\leq H$ and
$r_{1}\geq\frac{r_{1}+H}{2}$, thus it holds

$\displaystyle\mathrm{d}_{F_{1}}(r_{1}|\bar{r}_{1})=\int_{\bar{r}_{1}}^{r_{1}}F_{1}^{\prime\prime}(s)(r_{1}-s)\,\mathrm{d}s\geq\int_{H}^{\frac{r_{1}+H}{2}}F_{1}^{\prime\prime}(s)(r_{1}-s)\,\mathrm{d}s\geq\frac{r_{1}-H}{2}\,m_{H}H$

where we have used that $F_{1}^{\prime\prime}(s)\geq m_{H}$ and $r_{1}-s\geq
r_{1}-\frac{r_{1}+H}{2}\geq H$ for every
$s\in\mathopen{}\mathclose{{}\left[H,\frac{r_{1}+H}{2}}\right]$ by definition
of $m_{H}$ and $r_{1}\geq 3H$, respectively. From (15) and
$H\leq\frac{r_{1}-H}{2}$, it follows

$\displaystyle
r_{1}\mathopen{}\mathclose{{}\left|\theta_{j,i}(\bm{u})-\theta_{j,i}(\bar{\bm{u}})}\right|^{2}\leq
4\kappa_{j,i}^{2}r_{1}=4\kappa_{j,i}^{2}(r_{1}-H+H)\leq
6\kappa_{j,i}^{2}(r_{1}-H)\leq\frac{12\kappa_{j,i}^{2}}{m_{H}H}\,\mathrm{d}_{F_{1}}(r_{1}|\bar{r}_{1}).$

In the second case, we have $r_{1}\leq 3H$, and thus also $r_{2}\leq 3H$ since
we assumed $r_{1}\geq r_{2}$. By the monotonicity of $F_{1}^{\prime}$ and
$F_{2}^{\prime}$, this implies that
$\bm{u},\bar{\bm{u}}\in[0,F_{1}^{\prime}(3H)]\times[0,F_{2}^{\prime}(3H)]$,
which is a compact subset of ${\mathbb{R}}_{\geq 0}^{2}$. Hence by Hypothesis
2.2, there is a constant $L_{H}$ such that $\theta_{j,i}$ is $L_{H}$-Lipschitz
on this subset. Thus

$\displaystyle
r_{1}\mathopen{}\mathclose{{}\left|\theta_{j,i}(\bm{u})-\theta_{j,i}(\bar{\bm{u}})}\right|^{2}\leq
3HL_{H}^{2}\mathopen{}\mathclose{{}\left((u_{1}-\bar{u}_{1})^{2}+(u_{2}-\bar{u}_{2})^{2}}\right).$

We estimate $(u_{j}-\bar{u}_{j})^{2}$ by means of
$\mathrm{d}_{F_{j}}(r_{j}|\bar{r}_{j})$: We have

$\displaystyle(u_{j}-\bar{u}_{j})^{2}=\mathopen{}\mathclose{{}\left(F_{j}^{\prime}(r_{j})-F_{j}^{\prime}(\bar{r}_{j})}\right)^{2}=\int_{\bar{r}_{j}}^{r_{j}}2\mathopen{}\mathclose{{}\left(F_{j}^{\prime}(s)-F_{j}^{\prime}(\bar{r}_{j})}\right)F_{j}^{\prime\prime}(s)\,\mathrm{d}s\leq
2M_{H}\,\mathrm{d}_{F_{j}}(r_{j}|\bar{r}_{j}),$

since the term $F_{j}^{\prime}(s)-F_{j}^{\prime}(\bar{r}_{j})$ is negative if
and only if $r_{j}<\bar{r}_{j}$, in which case the endpoints $\bar{r}_{j}$ and
$r_{j}$ of the integral are swapped. This finishes the proof of (22). ∎

We additionally prove some important general properties of the energy
functional $\mathcal{E}_{\varepsilon}$ and its components. In all of the
following, we assume $\varepsilon\leq\varepsilon_{0}$ with $\varepsilon_{0}$
from Proposition 2.4.

######  Lemma 2.5.

Let $\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ be any pair and let
$\mathcal{S}\subset{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ be any family of
pairs.

  * •

The following estimates hold:

$\displaystyle\int_{\mathbb{R}^{d}}[F_{1}(\rho_{1})+F_{2}(\rho_{2})]\,\mathrm{d}x\leq
2\mathcal{E}_{\varepsilon}(\bm{\rho}),\quad\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x\leq\mathcal{E}_{\varepsilon}(\bm{\rho}).$
(23)

In particular, the energy $\mathcal{E}_{\varepsilon}(\bm{\rho})$ is uniformly
bounded for $\bm{\rho}\in\mathcal{S}$ if and only if both $F_{j}(\rho_{j})$
and $\rho_{1}K\ast\rho_{2}$ are uniformly bounded in $L^{1}({\mathbb{R}^{d}})$
for $\bm{\rho}\in\mathcal{S}$.

  * •

With the constants $\alpha_{j}$ from (21), it holds:

$\int_{\mathbb{R}^{d}}F_{j}^{\prime}(\rho_{j})\,\mathrm{d}x\leq\alpha_{j}\,(1+2\mathcal{E}_{\varepsilon}(\bm{\rho})).$
(24)

  * •

For $\bm{\rho}\in{\mathfrak{X}}_{0}$, the interaction energy
$\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x$ satisfies the bounds

$\frac{\lambda}{2}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right)\leq\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x\leq
C_{K}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right).$
(25)

In particular, for $\mathcal{S}\subset{\mathfrak{X}}_{0}$, $\bm{\rho}$’s
interaction energy is uniformly bounded for $\bm{\rho}\in\mathcal{S}$ if and
only if $\rho_{1}$ and $\rho_{2}$ have uniformly bounded second moments for
$\bm{\rho}\in\mathcal{S}$.

###### Proof.

The interaction energy is clearly non-negative by our assumptions on $K$, thus
the first estimate in (23) follows directly from (18). The second one holds
since $\mathcal{F}_{\varepsilon}$ is a convex function whose derivative
vanishes at 0, and thus is non-negative in ${\mathbb{R}}_{\geq 0}^{2}$. By
(21), it holds
$F_{j}^{\prime}(\rho_{j})\leq\alpha_{j}\,(\rho_{j}+F_{j}(\rho_{j}))$. From
this we obtain (24) by integrating and using (23) and the fact that $\rho_{j}$
has unit mass.

In order to prove (25), note that by $\lambda$-convexity of $K$ and the bound
on its second derivative, it holds $\frac{\lambda}{2}|z|^{2}\leq
K(z)\leq\frac{C_{K}}{2}|z|^{2}$. Since for $\bm{\rho}\in{\mathfrak{X}}_{0}$,
we have ${\mathfrak{m}}_{1}[\rho_{2}]=-{\mathfrak{m}}_{1}[\rho_{1}]$, we thus
obtain the lower bound

$\displaystyle\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x$
$\displaystyle=\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\rho_{1}(x)\rho_{2}(y)K(x-y)\,\mathrm{d}x\mathrm{d}y\geq\frac{\lambda}{2}\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\rho_{1}(x)\rho_{2}(y)|x-y|^{2}\,\mathrm{d}x\mathrm{d}y$
$\displaystyle=\frac{\lambda}{2}\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\rho_{1}(x)\rho_{2}(y)\mathopen{}\mathclose{{}\left[|x|^{2}+|y|^{2}-2x\cdot
y}\right]\,\mathrm{d}x\mathrm{d}y$
$\displaystyle=\frac{\lambda}{2}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right)-\lambda{\mathfrak{m}}_{1}[\rho_{1}]\cdot{\mathfrak{m}}_{1}[\rho_{2}]=\frac{\lambda}{2}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right)+\lambda|{\mathfrak{m}}_{1}[\rho_{1}]|^{2}$
$\displaystyle\geq\frac{\lambda}{2}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right).$

For the upper bound, the same argument yields

$\displaystyle\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x\leq\frac{C_{K}}{2}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right)+C_{K}|{\mathfrak{m}}_{1}[\rho_{1}]|^{2}\leq
C_{K}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}]}\right),$

where the last step follows from
$|{\mathfrak{m}}_{1}[\rho_{1}]|^{2}=\frac{1}{2}(|{\mathfrak{m}}_{1}[\rho_{1}]|^{2}+|{\mathfrak{m}}_{1}[\rho_{2}]|^{2})\leq\frac{1}{2}({\mathfrak{m}}_{2}[\rho_{1}]+{\mathfrak{m}}_{2}[\rho_{2}])$.
∎

In the next section, we shall use these properties to proof existence of a
global minimizer $\bar{\bm{\rho}}\in{\mathfrak{X}}_{0}$ of
$\mathcal{E}_{\varepsilon}$ via the direct method from the calculus of
variations.

##  3\. Existence and properties of the minimizer for the energy functional

###  3.1. Existence of minimizers

We start by proving the following weak coercivity result for the energy
$\mathcal{E}_{\varepsilon}$. The proof mainly uses the estimates from Lemma
2.5.

######  Lemma 3.1.

The sublevel sets of the energy functional $\mathcal{E}_{\varepsilon}$ in
${\mathfrak{X}}_{0}$ are pre-compact in the weak $L^{1}$-topology, i.e. any
sequence of pairs
$(\bm{\rho}^{n})=(\rho_{1}^{n},\rho_{2}^{n})\subset{\mathfrak{X}}_{0}$ with
$\mathcal{E}_{\varepsilon}(\bm{\rho}^{n})$ uniformly bounded contains a
subsequence converging weakly in $L^{1}({\mathbb{R}^{d}})$ to some
$\bm{\rho}\in{\mathfrak{X}}_{0}$. The same holds true for sequences
$(\bm{\rho}^{n})\subset{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ with
$\mathcal{E}_{\varepsilon}(\bm{\rho}^{n})$ uniformly bounded, such that their
combined center of mass
${\mathfrak{m}}_{1}[\rho_{1}^{n}]+{\mathfrak{m}}_{1}[\rho_{2}^{n}]$ converges
to $0$.

###### Proof.

Take any sequence $(\bm{\rho}^{n})\subset{\mathcal{P}_{2}}({\mathbb{R}^{d}})$
such that
${\mathfrak{m}}_{1}[\rho_{1}^{n}]+{\mathfrak{m}}_{1}[\rho_{2}^{n}]\to 0$ and
$\mathcal{E}_{\varepsilon}(\bm{\rho}^{n})$ is uniformly bounded. We consider
the shifted pair $\tilde{\bm{\rho}}^{n}:=\bm{\rho}^{n}(\cdot-v_{n})$ with the
shift
$v_{n}:=\frac{1}{2}({\mathfrak{m}}_{1}[\rho_{1}^{n}]+{\mathfrak{m}}_{1}[\rho_{2}^{n}])$.
We have $\tilde{\bm{\rho}}^{n}\in{\mathfrak{X}}_{0}$, and from the translation
invariance of $\mathcal{E}_{\varepsilon}$, it follows that
$\mathcal{E}_{\varepsilon}(\tilde{\bm{\rho}}^{n})=\mathcal{E}_{\varepsilon}(\bm{\rho}^{n})$
is uniformly bounded. It thus follows from (23) and (25) that the second
moments
${\mathfrak{m}}_{2}[\tilde{\rho}_{1}^{n}],{\mathfrak{m}}_{2}[\tilde{\rho}_{2}^{n}]$
are uniformly bounded. Together with $v_{n}\to 0$, this implies uniform
boundedness of the second moments of $\rho_{j}^{n}$:

$\displaystyle{\mathfrak{m}}_{2}[\rho_{j}^{n}]=\int_{\mathbb{R}^{d}}|x|^{2}\tilde{\rho}_{j}^{n}(x+v_{n})\,\mathrm{d}x=\int_{\mathbb{R}^{d}}|y-v_{n}|^{2}\tilde{\rho}_{j}^{n}(y)\,\mathrm{d}y\leq
2{\mathfrak{m}}_{2}[\tilde{\rho}_{j}^{n}]+2|v_{n}|^{2}.$

In particular, the sequence $(\rho_{j}^{n})$ is tight as probability measures
for $j=1,2$, and thus converges narrowly to a
$\rho_{j}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$ along some subsequence, again
denoted $(\rho_{j}^{n})$. Since narrow convergence and uniform boundedness of
the second moments imply
${\mathfrak{m}}_{1}[\rho_{j}^{n}]\to{\mathfrak{m}}_{1}[\rho_{j}]$, we have
$\bm{\rho}\in{\mathfrak{X}}_{0}$.

To prove $\rho_{j}^{n}\rightharpoonup\rho_{j}$ weakly in
$L^{1}({\mathbb{R}^{d}})$, observe that by (23), the integrals of
$F_{j}(\rho_{j}^{n})$ are uniformly bounded. Since $F_{j}(r)$ grows super-
linearly for $r\to+\infty$ by the at least linear growth of
$F_{j}^{\prime}(r)$, the de la Vallée-Poussin compactness criterion, see e.g.
[19, Theorem 1.28], yields compactness of the sequence $\rho_{j}^{n}$ with
respect to local weak $L^{1}$-convergence. Together with the narrow
convergence $\rho_{j}^{n}\rightharpoonup\rho_{j}$, this implies
$\rho_{j}^{n}\rightharpoonup\rho_{j}$ weakly in
$L^{1}_{loc}({\mathbb{R}^{d}})$. Together with tightness, this yields
$\rho_{j}^{n}\rightharpoonup\rho_{j}$ weakly in $L^{1}({\mathbb{R}^{d}})$,
proving the claim. ∎

In order to apply the direct method in the calculus of variations, we need to
prove lower semi-continuity of the functional $\mathcal{E}_{\varepsilon}$ in
the same topology.

######  Lemma 3.2.

The energy functional $\mathcal{E}_{\varepsilon}$ is lower semi-continuous in
${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ with respect to weak
$L^{1}$-convergence.

###### Proof.

Since $\mathcal{F}_{\varepsilon}$ is non-negative by (18), Fatou’s lemma
yields lower semi-continuity of the functional
$\bm{\rho}\mapsto\int_{\mathbb{R}^{d}}\mathcal{F}_{\varepsilon}(\bm{\rho})\,\mathrm{d}x$
with respect to strong $L^{1}$-convergence. Since $\mathcal{F}_{\varepsilon}$
is a convex function, this implies weak lower semi-continuity. It remains to
show that the interaction energy
$\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x$ is lower semi-
continuous with respect to weak $L^{1}$-convergence. Let
$\bm{\rho}^{n}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ with
$\rho_{j}^{n}\rightharpoonup\rho_{j}$ weakly in $L^{1}({\mathbb{R}^{d}})$, and
for any radius $R>0$, define

$\displaystyle K_{R}\in L^{\infty}({\mathbb{R}^{d}}),\
K_{R}(z):=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}K(z)&|z|\leq R\\\
0&|z|>R\end{array}}\right..$

Also note that $\rho_{j}^{n}\rightharpoonup\rho_{j}$ weakly in
$L^{1}({\mathbb{R}^{d}})$ for $j=1,2$ implies
$\rho_{1}^{n}\otimes\rho_{2}^{n}\rightharpoonup\rho_{1}\otimes\rho_{2}$ weakly
in $L^{1}({\mathbb{R}^{d}}\times{\mathbb{R}^{d}})$. Since $0\leq K_{R}\leq K$
in ${\mathbb{R}^{d}}$, this yields

$\displaystyle\liminf_{n\to\infty}\int_{\mathbb{R}^{d}}\rho_{1}^{n}K\ast\rho_{2}^{n}\,\mathrm{d}x\geq\liminf_{n\to\infty}\int_{\mathbb{R}^{d}}\rho_{1}^{n}K_{R}\ast\rho_{2}^{n}\,\mathrm{d}x+\liminf_{n\to\infty}\int_{\mathbb{R}^{d}}\rho_{1}^{n}(K-K_{R})\ast\rho_{2}^{n}\,\mathrm{d}x$
$\displaystyle\geq\liminf_{n\to\infty}\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\rho_{1}^{n}\otimes\rho_{2}^{n}(x,y)\,K_{R}(x-y)\,\mathrm{d}x\mathrm{d}y=\int_{\mathbb{R}^{d}}\rho_{1}K_{R}\ast\rho_{2}\,\mathrm{d}x.$

Since $R$ was arbitrary, we can take the limit $R\to\infty$. By the monotone
convergence theorem, the integral on the right-hand side converges to
$\int_{\mathbb{R}^{d}}\rho_{1}K\ast\rho_{2}\,\mathrm{d}x$, proving weak lower
semi-continuity of the interaction energy. ∎

Together, these two lemmas yield existence of a global minimizer of the energy
$\mathcal{E}_{\varepsilon}$ in ${\mathfrak{X}}_{0}$.

######  Corollary 3.3.

There exists a pair
$\bar{\bm{\rho}}=(\bar{\rho}_{1},\bar{\rho}_{2})\in{\mathfrak{X}}_{0}$,
depending on $\varepsilon\in(0,\varepsilon_{0}]$, that minimizes the energy
functional $\mathcal{E}_{\varepsilon}$ over
${\mathcal{P}_{2}}({\mathbb{R}^{d}})$.

###### Proof.

The existence of $\bar{\bm{\rho}}\in{\mathfrak{X}}_{0}$ that minimizes
$\mathcal{E}_{\varepsilon}$ over ${\mathfrak{X}}_{0}$ follows from the
previous two lemmas via the direct method from the calculus of variations.
Since $\mathcal{E}_{\varepsilon}$ is translation invariant, the minimality of
$\bar{\bm{\rho}}$ over ${\mathfrak{X}}_{0}$ implies that $\bar{\bm{\rho}}$
also minimizes $\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})$ over the whole
space ${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$. ∎

######  Remark 3.4.

Note that we do not obtain uniqueness of the minimizer $\bar{\bm{\rho}}$
directly, since the energy functional $\mathcal{E}_{\varepsilon}$ is not
convex on ${\mathfrak{X}}_{0}$, neither along geodesics, because of the
coupling term, nor in the flat sense, because of the interaction term. We
shall prove uniqueness of the minimizer $\bar{\bm{\rho}}$ for every
sufficiently small $\varepsilon>0$ in section 5, see Corollary 5.11.

###  3.2. Euler-Lagrange system

We now use a variational argument to prove that the minimizer
$\bar{\bm{\rho}}$ obtained from Corollary 3.3 satisfies the Euler-Lagrange
system (5).

######  Lemma 3.5.

Let $\psi_{1},\psi_{2}\in L^{1}({\mathbb{R}^{d}})\cap
L^{\infty}({\mathbb{R}^{d}})$ be compactly supported, and assume that
$\int_{\mathbb{R}^{d}}\psi_{1}=\int_{\mathbb{R}^{d}}\psi_{2}=0$ and
$\bar{\rho}_{j}+\psi_{j}\geq 0$ in ${\mathbb{R}^{d}}$ for $j=1,2$.
Additionally, assume that there exist constants $H_{j}>0$ such that
$\psi_{j}=0$ on the set $\\{\bar{\rho}_{j}>H_{j}\\}$. Then it holds

$\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K\ast\bar{\rho}_{2}}\right)\psi_{1}\,\mathrm{d}x\,+\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\partial_{r_{2}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K\ast\bar{\rho}_{1}}\right)\psi_{2}\,\mathrm{d}x\geq
0.$

###### Proof.

Let $\psi_{1}$ and $\psi_{2}$ as above. Then for any $s\in[0,1]$, the
perturbation
$\bm{\rho}^{s}:=\bar{\bm{\rho}}+s\mathopen{}\mathclose{{}\left(\psi_{1},\psi_{2}}\right)$
fulfills $\rho^{s}_{j}\geq 0$ and $\int\rho^{s}_{j}=1$ for $j=1,2$, and thus
defines a pair of probability distributions. Since $\bar{\bm{\rho}}$ minimizes
$\mathcal{E}_{\varepsilon}$ over ${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$, it
holds
$\mathcal{E}_{\varepsilon}(\bm{\rho}^{s})\geq\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})$
for all $s\in[0,1]$, hence

$\displaystyle\mathopen{}\mathclose{{}\left.\frac{\,\mathrm{d}}{\,\mathrm{d}s}}\right|_{s=0^{+}}\mathcal{E}_{\varepsilon}(\bm{\rho}^{s})=\mathopen{}\mathclose{{}\left.\frac{\,\mathrm{d}}{\,\mathrm{d}s}}\right|_{s=0^{+}}\mathcal{E}_{\varepsilon}(\bar{\rho}_{1}+s\psi_{1},\bar{\rho}_{2}+s\psi_{2})\geq
0\,.$

Using the properties of $\psi_{1}$, $\psi_{2}$, in particular the fact that
they are supported on a set where $\bar{\rho}_{j}$ is bounded, the dominated
convergence theorem applies and provides that

$\displaystyle 0$
$\displaystyle\leq\mathopen{}\mathclose{{}\left.\frac{\,\mathrm{d}}{\,\mathrm{d}s}}\right|_{s=0^{+}}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\mathcal{F}_{\varepsilon}(\bar{\rho}_{1}+s\psi_{1},\bar{\rho}_{2}+s\psi_{2})+(\bar{\rho}_{1}+s\psi_{1})K\ast(\bar{\rho}_{2}+s\psi_{2})}\right]\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(\partial_{r_{1}}\mathcal{F}_{\varepsilon}(\bar{\rho}_{1},\bar{\rho}_{2})+K\ast\bar{\rho}_{2}}\right)\psi_{1}\,\mathrm{d}x+\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(\partial_{r_{2}}\mathcal{F}_{\varepsilon}(\bar{\rho}_{1},\bar{\rho}_{2})+K\ast\bar{\rho}_{1}}\right)\psi_{2}\,\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K\ast\bar{\rho}_{2}}\right)\psi_{1}\,\mathrm{d}x$
$\displaystyle\qquad+\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\partial_{r_{2}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K\ast\bar{\rho}_{1}}\right)\psi_{2}\,\mathrm{d}x,$

which proves the asserted integral inequality. ∎

######  Lemma 3.6.

There are positive constants $C_{1},C_{2}\in{\mathbb{R}}$ such that
$\bar{\bm{\rho}}$ satisfies

$\begin{split}F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K*\bar{\rho}_{2}&\geq
C_{1}\\\
F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\partial_{r_{2}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K*\bar{\rho}_{1}&\geq
C_{2}\end{split}$ (26)

with equality a.e. on $\\{\bar{\rho}_{1}>0\\}$ and a.e. on
$\\{\bar{\rho}_{2}>0\\}$ respectively.

###### Proof.

For brevity, we introduce the notation

$EL_{1}:=F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})+K*\bar{\rho}_{2},$

and similarly $EL_{2}$. Since both inequalities in (26) are analogous to each
other, we only need to prove the first one. Suppose there does not exist a
constant $C_{1}$ such that $EL_{1}\equiv C_{1}$ a.e on
$\\{\bar{\rho}_{1}>0\\}$ and $EL_{1}\geq C_{1}$ a.e. on ${\mathbb{R}^{d}}$.
Then there exist constants $\alpha<\beta$, a height $H>0$ and a radius $R>0$
such that the bounded measurable sets

$M_{\alpha}:=\\{EL_{1}\leq\alpha\\}\cap\\{\bar{\rho}_{1}\leq H\\}\cap
B_{R},\qquad M_{\beta}:=\\{EL_{1}\geq\beta\\}\cap\\{0<\bar{\rho}_{1}\leq
H\\}\cap B_{R}$

both have strictly positive measure. Now define

$\psi_{1}:=-\bar{\rho}_{1}\mathds{1}_{M_{\beta}}+\frac{1}{|M_{\alpha}|}\mathopen{}\mathclose{{}\left(\int_{M_{\beta}}\bar{\rho}_{1}\,\mathrm{d}x}\right)\mathds{1}_{M_{\alpha}}.$

The properties $EL_{1}\leq\alpha$ on $M_{\alpha}$ and $EL_{1}\geq\beta$ on
$M_{\beta}$ yield

$\displaystyle\int_{\mathbb{R}^{d}}EL_{1}\,\psi_{1}\,\mathrm{d}x$
$\displaystyle=-\int_{M_{\beta}}EL_{1}\,\bar{\rho}_{1}\,\mathrm{d}x+\frac{1}{|M_{\alpha}|}\mathopen{}\mathclose{{}\left(\int_{M_{\beta}}\bar{\rho}_{1}\,\mathrm{d}x}\right)\int_{M_{\alpha}}EL_{1}\,\mathrm{d}x$
$\displaystyle\leq(\alpha-\beta)\int_{M_{\beta}}\bar{\rho}_{1}\,\mathrm{d}x<0,$

where the last inequality follows from $\bar{\rho}_{1}>0$ on $M_{\beta}$. This
is however a contradiction to Lemma 3.5 with $\psi_{2}\equiv 0$, since it is
easily checked that $\psi_{1}\in L^{1}({\mathbb{R}^{d}})$ satisfies all
properties needed to apply the lemma.

It remains to show that $C_{1}$ and $C_{2}$ are strictly positive. We observe
that the strict convexity of $\mathcal{F}_{\varepsilon}$ implies for $j=1,2$:

$F_{j}^{\prime}(r_{j})+\varepsilon\partial_{r_{j}}h({\mathbf{r}})=\partial_{r_{j}}\mathcal{F}_{\varepsilon}({\mathbf{r}})\geq
0\quad\text{for all }{\mathbf{r}}\in{\mathbb{R}}_{\geq 0}^{2},$ (27)

with equality if and only if $r_{j}=0$. The interaction terms
$K\ast\bar{\rho}_{j}$ are strictly positive on ${\mathbb{R}^{d}}$ by
positivity of $K$ in ${\mathbb{R}^{d}}\setminus\\{0\\}$ and the fact that
$\bar{\rho}_{j}$ is non-singular. Hence $EL_{j}(x)>0$ for all
$x\in{\mathbb{R}^{d}}$ and $j=1,2$, and thus $C_{1},C_{2}>0$. This concludes
the proof. ∎

######  Lemma 3.7.

For a.e. $x\in{\mathbb{R}^{d}}$, it holds

$\displaystyle K\ast\bar{\rho}_{2}(x)>C_{1}\quad$
$\displaystyle\Leftrightarrow\quad\bar{\rho}_{1}(x)=0$ $\displaystyle
K\ast\bar{\rho}_{2}(x)<C_{1}\quad$
$\displaystyle\Leftrightarrow\quad\bar{\rho}_{1}(x)>0$

and the same statements hold true for $\bar{\rho}_{2}(x)\geqq 0$ and
$K\ast\bar{\rho}_{1}(x)\lessgtr C_{2}$.

###### Proof.

If $K\ast\bar{\rho}_{2}(x)>C_{1}$, then it follows from (27) that (26) does
not hold with equality, hence up to a null-set, $\bar{\rho}_{1}(x)=0$. On the
other hand, if $K\ast\bar{\rho}_{2}(x)<C_{1}$, then (26) implies
$F_{1}^{\prime}(\bar{\rho}_{1}(x))+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1}(x),\bar{\rho}_{2}(x))>0$,
and thus $\bar{\rho}_{1}(x)>0$.

It follows immediately that $\bar{\rho}_{1}(x)=0$ implies
$K\ast\bar{\rho}_{2}(x)\geq C_{1}$ and that $\bar{\rho}_{1}(x)>0$ implies
$K\ast\bar{\rho}_{2}(x)\leq C_{1}$ for a.e. $x$. Since $K$, and thus
$K\ast\bar{\rho}_{2}$, is $\lambda$-convex on ${\mathbb{R}^{d}}$, the set
$\\{x:K\ast\bar{\rho}_{2}(x)=C_{1}\\}$ has measure zero, which concludes the
proof. ∎

From the previous two results, we obtain the Euler-Lagrange system and a
characterization of the supports of $\bar{\rho}_{1}$ and $\bar{\rho}_{2}$:

######  Theorem 3.8.

Any global minimizer $\bar{\bm{\rho}}\in{\mathfrak{X}}_{0}$ of the energy
functional $\mathcal{E}_{\varepsilon}$ satisfies the Euler-Lagrange system

$\begin{split}F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})&=\mathopen{}\mathclose{{}\left(C_{1}-K*\bar{\rho}_{2}}\right)_{+}\\\
F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\partial_{r_{2}}h(\bar{\rho}_{1},\bar{\rho}_{2})&=\mathopen{}\mathclose{{}\left(C_{2}-K*\bar{\rho}_{1}}\right)_{+}\end{split}$
(28)

in ${\mathbb{R}^{d}}$. Moreover, the supports of $\bar{\rho}_{1}$ and
$\bar{\rho}_{2}$ are given by the compact convex sets

$\begin{split}\mathrm{supp}\,\bar{\rho}_{1}&=\\{K\ast\bar{\rho}_{2}\leq
C_{1}\\}\\\ \mathrm{supp}\,\bar{\rho}_{2}&=\\{K\ast\bar{\rho}_{1}\leq
C_{2}\\}.\end{split}$ (29)

###### Proof.

The form (28) of the Euler-Lagrange system follows immediately from Lemmas 3.6
and 3.7, and the characterization (29) of the supports follows from Lemma 3.7
and the observation
$\overline{\\{K\ast\bar{\rho}_{2}<C_{1}\\}}=\\{K\ast\bar{\rho}_{2}\leq
C_{1}\\}$, which holds by $\lambda$-convexity of $K\ast\bar{\rho}_{2}$. ∎

######  Remark 3.9.

In the case that $K=\frac{\lambda}{2}|\cdot|^{2}$ is the quadratic interaction
potential, the Euler-Lagrange system can be written as

$\begin{split}F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\partial_{r_{1}}h(\bar{\rho}_{1},\bar{\rho}_{2})&=(\tilde{C}_{1}-\frac{\lambda}{2}|x|^{2})_{+}\\\
F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\partial_{r_{2}}h(\bar{\rho}_{1},\bar{\rho}_{2})&=(\tilde{C}_{2}-\frac{\lambda}{2}|x|^{2})_{+}.\end{split}$
(30)

with the constants $\tilde{C}_{1},\tilde{C}_{2}>0$ chosen uniquely such that
$\int_{\mathbb{R}^{d}}\bar{\rho}_{j}\,\mathrm{d}x=1$ for $j=1,2$. It holds
$\tilde{C}_{j}=C_{j}-\frac{\lambda}{2}{\mathfrak{m}}_{2}[\bar{\rho}_{j^{\prime}}]$
with the constants $C_{j}$ from (28); recall that $j^{\prime}:=3-j$. The form
(30) is considerably simpler to solve than (28), since the right-hand side
does not explicitly depend on $\bar{\bm{\rho}}$. It can be proven by showing
that in the case $K=\frac{\lambda}{2}|\cdot|^{2}$, the energy functional
$\mathcal{E}_{\varepsilon}$ is strictly convex in the flat sense on
${\mathfrak{X}}_{0}$, implying that the minimizer
$\bar{\bm{\rho}}\in{\mathfrak{X}}_{0}$ is unique. This uniqueness yields
radial symmetry of $\bar{\bm{\rho}}$, see Lemma 5.12, hence
${\mathfrak{m}}_{1}[\rho_{j}]=0$ for both components. A direct calculation
then yields
$K\ast\bar{\rho}_{j}=\frac{\lambda}{2}\mathopen{}\mathclose{{}\left[|\cdot|^{2}+{\mathfrak{m}}_{2}[\bar{\rho}_{j}]}\right]$,
thus (30) follows from (28).

###  3.3. Properties of the minimizers

In the next results, we prove some additional properties of the minimizers
$\bar{\bm{\rho}}$.

######  Lemma 3.10.

For any $\bar{\varepsilon}\in(0,\varepsilon_{0}]$, the minimal energy
$\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})$, the constants $C_{1},C_{2}$ from
the Euler-Lagrange system (28) and the diameters of
$\mathrm{supp}\,\bar{\rho}_{1}$ and $\mathrm{supp}\,\bar{\rho}_{2}$ are
$\varepsilon$-uniformly bounded for $\varepsilon\in(0,\bar{\varepsilon}]$. If
$\bar{\varepsilon}<1/\max\\{\kappa_{1,1},\kappa_{2,2}\\}$, there exists an
$\varepsilon$-independent constant $H_{0}>0$ such that $\bar{\rho}_{j}\leq
H_{0}$ in ${\mathbb{R}^{d}}$ for $j=1,2$ and any
$\varepsilon\in(0,\bar{\varepsilon}]$.

###### Proof.

Since for any fixed $\bm{\rho}\in{\mathfrak{X}}_{0}$, the energy
$\mathcal{E}_{\varepsilon}(\bm{\rho})$ is $\varepsilon$-uniformly bounded for
small $\varepsilon$, the minimal value
$\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})$ is uniformly bounded as well. To
prove the other claims, we write the Euler-Lagrange system using the functions
$\theta_{j}$:

$\begin{split}F_{1}^{\prime}(\bar{\rho}_{1})+\varepsilon\theta_{1}(F_{1}^{\prime}(\bar{\rho}_{1}),F_{2}^{\prime}(\bar{\rho}_{2}))=(C_{1}-K\ast\bar{\rho}_{2})_{+}\\\
F_{2}^{\prime}(\bar{\rho}_{2})+\varepsilon\theta_{2}(F_{1}^{\prime}(\bar{\rho}_{1}),F_{2}^{\prime}(\bar{\rho}_{2}))=(C_{2}-K\ast\bar{\rho}_{1})_{+}\end{split}$
(31)

It follows from assumption (15) that
$|\theta_{j}(\bm{u})|\leq\kappa_{j,i}u_{i}$ for all
$\bm{u}\in{\mathbb{R}}_{\geq 0}^{2}$ and any $i,j=1,2$. It thus follows from
(31) that

$\displaystyle(1-\varepsilon\kappa_{1,1})\,F_{1}^{\prime}(\bar{\rho}_{1})\leq(C_{1}-K\ast\bar{\rho}_{2})_{+}\leq(1+\varepsilon\kappa_{1,1})\,F_{1}^{\prime}(\bar{\rho}_{1}),$
(32)

and analogously for $F_{2}^{\prime}(\bar{\rho}_{2})$. For the integral of the
right-hand side, we have by (21)

$\displaystyle(1+\varepsilon\kappa_{1,1})\int_{\mathbb{R}^{d}}F_{1}^{\prime}(\bar{\rho}_{1})\,\mathrm{d}x\leq(1+\varepsilon\kappa_{1,1})\,\alpha_{1}\,\mathopen{}\mathclose{{}\left(1+2\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})}\right),$

which is $\varepsilon$-uniformly bounded for
$\varepsilon\in(0,\bar{\varepsilon}]$. Hence by (32), the integral of
$(C_{1}-K\ast\bar{\rho}_{2})_{+}$ is $\varepsilon$-uniformly bounded as well.
With $\bar{x}:=\operatorname*{argmin}_{x}K\ast\bar{\rho}_{2}$, we have, using
the bound $C_{K}$ on the Hessian of $K\ast\bar{\rho}_{2}$ and inequality (23),

$\displaystyle K\ast\bar{\rho}_{2}(x)\leq
K\ast\bar{\rho}_{2}(\bar{x})+\frac{C_{K}}{2}|x-\bar{x}|^{2}\leq\int_{\mathbb{R}^{d}}\bar{\rho}_{1}K\ast\bar{\rho}_{2}\,\mathrm{d}x+\frac{C_{K}}{2}|x-\bar{x}|^{2}\leq\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})+\frac{C_{K}}{2}|x-\bar{x}|^{2}$

for every $x\in{\mathbb{R}^{d}}$. Hence we can estimate the integral of
$(C_{1}-K\ast\bar{\rho}_{2})_{+}$ from below by

$\displaystyle\int_{\mathbb{R}^{d}}(C_{1}-K\ast\bar{\rho}_{2})_{+}\,\mathrm{d}x\geq\int_{\mathbb{R}^{d}}(C_{1}-\mathcal{E}_{\varepsilon}(\bar{\bm{\rho}})-C_{K}|x-\bar{x}|^{2})_{+}\,\mathrm{d}x.$

Since the left-hand side of this estimate is $\varepsilon$-uniformly bounded,
the same is true for the right-hand side, which implies that $C_{1}$ is
$\varepsilon$-uniformly bounded for $\varepsilon\in(0,\bar{\varepsilon}]$.
Since it holds $\frac{\lambda}{2}|x-\bar{x}|^{2}\leq K\ast\bar{\rho}_{2}\leq
C_{1}$ in $\mathrm{supp}\,\bar{\rho}_{1}$, this directly implies that the
diameter of $\mathrm{supp}\,\bar{\rho}_{1}$ is $\varepsilon$-uniformly bounded
for $\varepsilon\in(0,\bar{\varepsilon}]$. If
$\bar{\varepsilon}<1/\max\\{\kappa_{1,1},\kappa_{2,2}\\}$, implying that
$1-\varepsilon\kappa_{1,1}$ is bounded away from 0 uniformly for
$\varepsilon\in(0,\bar{\varepsilon}]$, the uniform boundedness of $C_{1}$ and
non-negativity of $K\ast\bar{\rho}_{2}$ together with (32) imply uniform
boundedness of $F_{1}^{\prime}(\bar{\rho}_{1})$, and thus also for
$\bar{\rho}_{1}$ by the growth of $F_{1}^{\prime}$. The proof for the other
component is analogous. ∎

In the next lemma, we prove an $\varepsilon$-uniform regularity result for the
minimizer $\bar{\bm{\rho}}$.

######  Lemma 3.11.

Assuming that $\bar{\varepsilon}>0$ is sufficiently small, the functions
$F_{j}^{\prime}(\bar{\rho}_{j})$ are Lipschitz in ${\mathbb{R}^{d}}$, and the
gradients $\nabla F_{j}^{\prime}(\bar{\rho}_{j})$ are Lipschitz in the common
support $\\{\bar{\rho}_{1}>0\\}\cap\\{\bar{\rho}_{2}>0\\}$, both with
Lipschitz constants independent of $\varepsilon\in(0,\bar{\varepsilon}]$.

###### Proof.

We consider the function
$\Gamma_{\varepsilon}:{\mathbb{R}}^{2}\to{\mathbb{R}}^{2}$, defined by

$\Gamma_{\varepsilon}(\bm{u}):=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}\bm{u}+\varepsilon\begin{pmatrix}\theta_{1}(\bm{u})\\\
\theta_{2}(\bm{u})\end{pmatrix}=\nabla\mathcal{F}_{\varepsilon}\mathopen{}\mathclose{{}\left({\mathbf{r}}}\right)&\text{
for }\bm{u}\in{\mathbb{R}}_{\geq 0}^{2}\\\ \bm{u}&\text{
else}\end{array}}\right.,$

where
${\mathbf{r}}=((F_{1}^{\prime})^{-1}(u_{1}),(F_{2}^{\prime})^{-1}(u_{2}))$ for
$\bm{u}\in{\mathbb{R}}_{\geq 0}^{2}$. It follows from the form (31) of the
Euler-Lagrange system that

$\Gamma_{\varepsilon}\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\bar{\rho}_{1}),F_{2}^{\prime}(\bar{\rho}_{2})}\right)=\begin{pmatrix}(C_{1}-K\ast\bar{\rho}_{2})_{+}\\\
(C_{2}-K\ast\bar{\rho}_{1})_{+}\end{pmatrix}.$ (33)

By strict convexity of $\mathcal{F}_{\varepsilon}$, the map
$\Gamma_{\varepsilon}$ is injective. Since both $\theta_{j}$ vanish on
$\partial{\mathbb{R}}_{\geq 0}^{2}$, we have
$\Gamma_{\varepsilon}=\mathrm{id}$ there, and $\Gamma_{\varepsilon}$ is
continuous in ${\mathbb{R}}^{2}$. To prove that $\Gamma_{\varepsilon}$ maps
${\mathbb{R}}_{>0}^{2}$ to itself surjectively, take any
$\bm{v}=(v_{1},v_{2})\in{\mathbb{R}}_{>0}^{2}$ and consider the minimizer
$\bar{\mathbf{r}}$ of the function
${\mathbf{r}}\mapsto\mathcal{F}_{\varepsilon}({\mathbf{r}})-\bm{v}\cdot{\mathbf{r}}$
over ${\mathbb{R}}_{\geq 0}^{2}$. The minimizer exists by the superlinear
growth of $\mathcal{F}_{\varepsilon}$, and it does not lie on the boundary
$\partial{\mathbb{R}}_{\geq 0}^{2}$, as for any ${\mathbf{r}}$ with $r_{i}=0$,
it holds
$\partial_{r_{i}}(\mathcal{F}_{\varepsilon}({\mathbf{r}})-\bm{v}\cdot{\mathbf{r}})=-v_{i}<0$.
We thus have $\bar{\mathbf{r}}\in{\mathbb{R}}_{>0}^{2}$, and thus by the
necessary first-order minimality condition, it holds
$\nabla\mathcal{F}_{\varepsilon}(\bar{\mathbf{r}})=\bm{v}$, implying that
$\Gamma_{\varepsilon}(F_{1}^{\prime}(\bar{r}_{1}),F_{2}^{\prime}(\bar{r}_{2}))=\bm{v}$.
This shows that $\Gamma_{\varepsilon}$ maps ${\mathbb{R}}^{2}$,
${\mathbb{R}}_{\geq 0}^{2}$ and ${\mathbb{R}}_{>0}^{2}$ onto themselves
bijectively, and we can rewrite (33) as

$\begin{pmatrix}F_{1}^{\prime}(\bar{\rho}_{1})\\\
F_{2}^{\prime}(\bar{\rho}_{2})\end{pmatrix}=\Gamma_{\varepsilon}^{-1}\begin{pmatrix}(C_{1}-K\ast\bar{\rho}_{2})_{+}\\\
(C_{2}-K\ast\bar{\rho}_{1})_{+}\end{pmatrix}.$ (34)

We claim that for $\bar{\varepsilon}$ sufficiently small, the inverse mapping
$\Gamma_{\varepsilon}^{-1}$ is of class $C^{1}$, and its Jacobian
${\mathrm{D}}\Gamma_{\varepsilon}^{-1}$ is globally bounded and locally
Lipschitz in ${\mathbb{R}}^{2}$, uniformly with respect to
$\varepsilon\in(0,\bar{\varepsilon}]$. We prove the claim by the inverse
function theorem. In ${\mathbb{R}}_{>0}^{2}$, we have with (15)

$\displaystyle\det{\mathrm{D}}\Gamma_{\varepsilon}(\bm{u})=\det\mathopen{}\mathclose{{}\left[\mathrm{Id}+\varepsilon\begin{pmatrix}\theta_{1,1}(\bm{u})&\theta_{1,2}(\bm{u})\\\
\theta_{2,1}(\bm{u})&\theta_{2,2}(\bm{u})\end{pmatrix}}\right]\geq
1-\varepsilon\mathopen{}\mathclose{{}\left(\kappa_{1,1}+\kappa_{2,2}}\right)-\varepsilon^{2}\mathopen{}\mathclose{{}\left(\kappa_{1,1}\kappa_{2,2}+\kappa_{1,2}\kappa_{2,1}}\right),$

which is $\varepsilon$-uniformly bounded away from 0 for
$\varepsilon\in(0,\bar{\varepsilon}]$, if $\bar{\varepsilon}>0$ is
sufficiently small. In ${\mathbb{R}}^{2}\setminus{\mathbb{R}}_{>0}^{2}$, we
have ${\mathrm{D}}\Gamma_{\varepsilon}=\mathrm{Id}$. Since the $\theta_{j,i}$
vanish on $\partial{\mathbb{R}}_{\geq 0}^{2}$ by hypothesis (15), and since
they are locally Lipschitz, ${\mathrm{D}}\Gamma_{\varepsilon}$ is locally
Lipschitz on ${\mathbb{R}}^{2}$, uniformly for
$\varepsilon\in(0,\bar{\varepsilon}]$. Hence the inverse function theorem is
applicable and yields $C^{1}$-regularity of the inverse
$\Gamma_{\varepsilon}^{-1}$ on ${\mathbb{R}}^{2}$, with Jacobian
${\mathrm{D}}\Gamma_{\varepsilon}^{-1}=\mathopen{}\mathclose{{}\left({\mathrm{D}}\Gamma_{\varepsilon}\circ\Gamma_{\varepsilon}^{-1}}\right)^{-1}$.
Since $\det{\mathrm{D}}\Gamma_{\varepsilon}$ is $\varepsilon$-uniformly
bounded away from 0 on ${\mathbb{R}}^{2}$, this implies that
$D\Gamma_{\varepsilon}^{-1}$ is $\varepsilon$-uniformly bounded in
${\mathbb{R}}^{2}$, implying that $\Gamma_{\varepsilon}^{-1}$ is globally
Lipschitz on ${\mathbb{R}}^{2}$ with an $\varepsilon$-independent Lipschitz
constant for all $\varepsilon\in(0,\bar{\varepsilon}]$. Since the inverse
mapping $A\mapsto A^{-1}$ is Lipschitz on any set on $2\times 2$-matrices $A$
with $\det A$ bounded away from 0, this also implies that
${\mathrm{D}}\Gamma_{\varepsilon}^{-1}$ is locally Lipschitz on
${\mathbb{R}}^{2}$, uniformly in $\varepsilon$, since it is a composition of
locally Lipschitz functions. This finishes the proof of the claim.

By Lemma 3.10, the diameter of $\mathrm{supp}\,\bar{\rho}_{j}$ is
$\varepsilon$-uniformly bounded. Together with the bound $|\nabla K(z)|\leq
C_{K}|z|$, this implies that $(C_{1}-K\ast\bar{\rho}_{2})_{+}$ and
$(C_{2}-K\ast\bar{\rho}_{1})_{+}$ are globally Lipschitz in
${\mathbb{R}}^{d}$, uniformly in $\varepsilon$. It thus follows from (34) and
the $\varepsilon$-uniform Lipschitz continuity of $\Gamma_{\varepsilon}^{-1}$
that the $F_{j}^{\prime}(\bar{\rho}_{j})$ are $\varepsilon$-uniformly
Lipschitz in ${\mathbb{R}}^{d}$. On the set
$\\{\bar{\rho}_{1}>0\\}\cap\\{\bar{\rho}_{2}>0\\}$, it holds
$(C_{1}-K\ast\bar{\rho}_{2})_{+}=C_{1}-K\ast\bar{\rho}_{2}$ and
$(C_{2}-K\ast\bar{\rho}_{1})_{+}=C_{2}-K\ast\bar{\rho}_{1}$ by Lemma 3.7,
which are of class $C^{2}$ by the regularity of $K$. Thus, their gradients are
$\varepsilon$-uniformly Lipschitz in
$\\{\bar{\rho}_{1}>0\\}\cap\\{\bar{\rho}_{2}>0\\}$, and (34) together with the
local Lipschitz-continuity of ${\mathrm{D}}\Gamma_{\varepsilon}^{-1}$ implies
$\varepsilon$-uniform Lipschitz-continuity of $\nabla
F_{j}^{\prime}(\bar{\rho}_{j})$ in the common support of $\bar{\rho}_{1}$ and
$\bar{\rho}_{2}$. Note that local Lipschitz-continuity of
${\mathrm{D}}\Gamma_{\varepsilon}^{-1}$ is sufficient, since the constants
$C_{1},C_{2}$ are $\varepsilon$-uniformly bounded by Lemma 3.10 and thus
${\mathrm{D}}\Gamma_{\varepsilon}^{-1}$ is $\varepsilon$-uniformly Lipschitz
on the compact set $[0,C_{1}]\times[0,C_{2}]\subset{\mathbb{R}}^{2}$. This
finishes the proof. ∎

An important consequence of this lemma is the following semiconvexity result,
which we shall use in section 5 to prove exponential convergence to the
equilibrium for system (1).

######  Corollary 3.12.

For the same $\bar{\varepsilon}>0$ as in the previous lemma, there exists an
$\varepsilon$-independent constant $K_{0}>0$ such that
$V_{j}:=\partial_{r_{j}}h(\bar{\rho}_{1},\bar{\rho}_{2})$ is
$(-K_{0})$-semiconvex in ${\mathbb{R}^{d}}$ for $j=1,2$ and every
$\varepsilon\in(0,\bar{\varepsilon}]$.

###### Proof.

We show that for some $\varepsilon$-independent constant $K_{0}$, the gradient
$\nabla V_{j}$ is $K_{0}$-Lipschitz in ${\mathbb{R}^{d}}$ for $j=1,2$ and
$\varepsilon\in(0,\bar{\varepsilon}]$. With
$\bar{\bm{u}}:=(F_{1}^{\prime}(\bar{\rho}_{1}),F_{2}^{\prime}(\bar{\rho}_{2}))$,
it holds

$\nabla V_{j}=\nabla\theta_{j}(\bar{\bm{u}})=\theta_{j,1}(\bar{\bm{u}})\nabla
F_{1}^{\prime}(\bar{\rho}_{1})+\theta_{j,2}(\bar{\bm{u}})\nabla
F_{2}^{\prime}(\bar{\rho}_{2}),$ (35)

which is nonzero only in $\\{\bar{\rho}_{1}>0\\}\cap\\{\bar{\rho}_{2}>0\\}$,
since the $\theta_{j,i}$ vanish on $\partial{\mathbb{R}}_{\geq 0}^{2}$ by
assumption (15). Inside the common support of $\bar{\rho}_{1}$ and
$\bar{\rho}_{2}$, the $\nabla F_{j}^{\prime}(\bar{\rho}_{j})$ are
$\varepsilon$-uniformly Lipschitz by Lemma 3.11. Since by Lemma 3.10, the
$\bar{\rho}_{j}$ are $\varepsilon$-uniformly bounded by a constant $H_{0}$,
the function $\bar{\bm{u}}$ only attains values in the compact set
$[0,F_{1}^{\prime}(H_{0})]\times[0,F_{2}^{\prime}(H_{0})]$ for all
$\varepsilon\in(0,\bar{\varepsilon}]$. By Hypothesis 2.2, the $\theta_{j,i}$
are Lipschitz on this set, and $\bar{\bm{u}}$ is globally Lipschitz by Lemma
3.11, uniformly in $\varepsilon$, hence $\theta_{j,i}(\bar{\bm{u}})$ is
$\varepsilon$-uniformly Lipschitz in ${\mathbb{R}^{d}}$. This proves that the
expression (35) is Lipschitz in
$\\{\bar{\rho}_{1}>0\\}\cap\\{\bar{\rho}_{2}>0\\}$. Additionally, it is
continuous across the boundary of the common support. To see this, take any
$x^{*}$ on the boundary, and any sequence $x_{n}\to x^{*}$ with
$x_{n}\in\\{\bar{\rho}_{1}>0\\}\cap\\{\bar{\rho}_{2}>0\\}$. Since $x^{*}$ lies
on the boundary and both $F_{j}^{\prime}(\bar{\rho}_{j})$ are continuous in
${\mathbb{R}^{d}}$, at least one of the two sequences
$\bar{u}_{k}(x_{n})=F_{k}^{\prime}(\bar{\rho}_{k}(x_{n}))$ for $k=1,2$
converges to 0, implying by (15) that $\theta_{j,i}(\bar{\bm{u}}(x_{n}))\to 0$
for both $i$. Since the $\nabla F_{j}^{\prime}$-terms are bounded along the
sequence $(x_{n})$, this implies that the expression (35) evaluated at $x_{n}$
converges to 0, yielding continuity across the boundary. This proves that the
$\nabla V_{j}$ are $\varepsilon$-uniformly Lipschitz in ${\mathbb{R}^{d}}$
with some Lipschitz constant $K_{0}$, yielding $(-K_{0})$-semiconvexity on
$V_{j}$ for every $\varepsilon\in(0,\bar{\varepsilon}]$. ∎

##  4\. Minimizing movement scheme for the energy functional

Fixing any $\tau>0$, we introduce the functional
$\mathcal{E}_{\varepsilon,\tau}$, defined for
$\bm{\rho},\hat{\bm{\rho}}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ as

$\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}|\hat{\bm{\rho}}):=\frac{1}{2\tau}\mathrm{d}(\bm{\rho},\hat{\bm{\rho}})^{2}+\mathcal{E}_{\varepsilon}(\bm{\rho}).$

Our goal is to get a time-discrete variational approximation for a weak
solution of system (1) by iteratively minimizing the functional
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$ over
${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$, with $\hat{\bm{\rho}}$ being the
minimizing pair obtained in the previous iteration.

###  4.1. Existence of minimizers for the modified functional

We start by proving existence of global minimizers of the functional
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$.

######  Lemma 4.1.

For every pair $\hat{\bm{\rho}}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$,
the functional
${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}\ni\bm{\rho}\mapsto\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}|\hat{\bm{\rho}})$
possesses a global minimizer $\bm{\rho}^{+}$. Every minimizer $\bm{\rho}^{+}$
satisfies the inequality

$\frac{1}{2\tau}\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}})^{2}\leq\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}).$
(36)

If $\hat{\bm{\rho}}\in{\mathfrak{X}}_{0}$, then also
$\bm{\rho}^{+}\in{\mathfrak{X}}_{0}$ for every global minimizer
$\bm{\rho}^{+}$.

###### Proof.

We first claim that for any
$\hat{\bm{\rho}}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$, the minimization
problem of $\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$ over
${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ is equivalent to the restricted
minimization problem over the set of all
$\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ which have the same
combined center of mass as $\hat{\bm{\rho}}$, i.e.
${\mathfrak{m}}_{1}[\rho_{1}]+{\mathfrak{m}}_{1}[\rho_{2}]={\mathfrak{m}}_{1}[\hat{\rho}_{1}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}]$.

To prove the claim, fix any
$\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ with
$\mathcal{E}_{\varepsilon}(\bm{\rho})<+\infty$, so in particular, the
$\rho_{j}$ for $j=1,2$ are absolutely continuous with respect to the Lebesgue
measure. Denote by $\Phi_{j}:{\mathbb{R}^{d}}\to{\mathbb{R}}$ convex functions
such that $\nabla\Phi_{j}:{\mathbb{R}^{d}}\to{\mathbb{R}^{d}}$ are optimal
transport maps from $\rho_{j}$ to $\hat{\rho}_{j}$, which exist by Brenier’s
theorem. Now consider the pair of translated densities
$\sigma_{v}\bm{\rho}=\mathopen{}\mathclose{{}\left(\rho_{1}(\cdot+v),\rho_{2}(\cdot+v)}\right)$,
for any direction $v\in{\mathbb{R}^{d}}$. Clearly, the maps
$\nabla\Phi_{j}(\cdot+v)$ push $\sigma_{v}\rho_{j}$ forward to
$\hat{\rho}_{j}$, and since they are the gradients of the convex functions
$\Phi_{j}(\cdot+v)$, they are optimal in the transport problem from
$\sigma_{v}\rho_{j}$ to $\hat{\rho}_{j}$. This implies

$\displaystyle{\mathbf{W}}(\sigma_{v}\rho_{j},\hat{\rho}_{j})^{2}$
$\displaystyle=\int_{\mathbb{R}^{d}}|\nabla\Phi_{j}(x+v)-x|^{2}\,\rho_{j}(x+v)\,\mathrm{d}x=\int_{\mathbb{R}^{d}}|\nabla\Phi_{j}(y)-y+v|^{2}\,\rho_{j}(y)\,\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}^{d}}|\nabla\Phi_{j}(y)-y|^{2}\,\rho_{j}(y)\,\mathrm{d}y+|v|^{2}+2\,\mathopen{}\mathclose{{}\left(\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(\nabla\Phi_{j}(y)-y}\right)\,\rho_{j}(y)\,\mathrm{d}y}\right)\cdot
v$
$\displaystyle={\mathbf{W}}(\rho_{j},\hat{\rho}_{j})^{2}+|v|^{2}+2\,\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{1}[\hat{\rho}_{j}]-{\mathfrak{m}}_{1}[\rho_{j}]}\right)\cdot
v\,.$

Combining this fact with the translation invariance of the energy
$\mathcal{E}_{\varepsilon}$, we obtain

$\displaystyle\mathcal{E}_{\varepsilon,\tau}(\sigma_{v}\bm{\rho}|\hat{\bm{\rho}})=\frac{{\mathbf{W}}(\sigma_{v}\rho_{1},\hat{\rho}_{1})^{2}+{\mathbf{W}}(\sigma_{v}\rho_{2},\hat{\rho}_{2})^{2}}{2\tau}+\mathcal{E}_{\varepsilon}(\sigma_{v}\bm{\rho})$
$\displaystyle=\frac{{\mathbf{W}}(\rho_{1},\hat{\rho}_{1})^{2}+{\mathbf{W}}(\rho_{2},\hat{\rho}_{2})^{2}+2|v|^{2}+2\,({\mathfrak{m}}_{1}[\hat{\rho}_{1}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}]-{\mathfrak{m}}_{1}[\rho_{1}]-{\mathfrak{m}}_{1}[\rho_{2}])\cdot
v}{2\tau}+\mathcal{E}_{\varepsilon}(\bm{\rho})$
$\displaystyle=\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}|\hat{\bm{\rho}})+\frac{|v|^{2}+\,({\mathfrak{m}}_{1}[\hat{\rho}_{1}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}]-{\mathfrak{m}}_{1}[\rho_{1}]-{\mathfrak{m}}_{1}[\rho_{2}])\cdot
v}{\tau}\,.$

This expression is quadratic in $v\in{\mathbb{R}^{d}}$, and is uniquely
minimized by $\bar{v}$ given by

$\displaystyle\bar{v}:=-\frac{1}{2}\mathopen{}\mathclose{{}\left({\mathfrak{m}}_{1}[\hat{\rho}_{1}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}]-{\mathfrak{m}}_{1}[\rho_{1}]-{\mathfrak{m}}_{1}[\rho_{2}]}\right).$

For the combined center of mass of this optimally shifted pair
$\sigma_{\bar{v}}\bm{\rho}$, we obtain

$\displaystyle{\mathfrak{m}}_{1}[\sigma_{\bar{v}}\rho_{1}]+{\mathfrak{m}}_{1}[\sigma_{\bar{v}}\rho_{2}]={\mathfrak{m}}_{1}[\rho_{1}]+{\mathfrak{m}}_{1}[\rho_{2}]-2\bar{v}={\mathfrak{m}}_{1}[\hat{\rho}_{1}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}].$

This proves the claimed equivalence of the minimization problems over the
whole ${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ and the restricted set of
pairs which preserve $\hat{\bm{\rho}}$’s combined center of mass.

To show existence of a global minimizer, we can thus w.l.o.g. assume
$\hat{\bm{\rho}}\in{\mathfrak{X}}_{0}$ and prove existence of a minimizer
$\bm{\rho}^{+}$ over ${\mathfrak{X}}_{0}$. Similarly as in the previous
section, we shall use the direct method from the calculus of variations. Since
the distance term is non-negative, pre-compactness of the sublevel sets of
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$ follows directly from
Lemma 3.1. Weak lower semicontinuity follows from Lemma 3.2 and lower semi-
continuity of the Wasserstein distance with respect to weak
$L^{1}$-convergence, proving the existence of a minimizer
$\bm{\rho}^{+}\in{\mathfrak{X}}_{0}$.

From the minimality of $\bm{\rho}^{+}$, we directly obtain the inequality

$\displaystyle\frac{1}{2\tau}\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}}^{+})^{2}+\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})=\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}^{+}|\hat{\bm{\rho}})\leq\mathcal{E}_{\varepsilon,\tau}(\hat{\bm{\rho}}|\hat{\bm{\rho}})=\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}),$

which implies (36), thus finishing the proof. ∎

As a first regularity result, we prove a global $L^{2}$-estimate for $\nabla
F_{j}^{\prime}(\rho_{j}^{+})$, restricted to a set where $\rho_{j}^{+}$ is
smaller than some fixed height $\alpha>0$.

######  Lemma 4.2.

Let $\hat{\bm{\rho}}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ be any fixed
pair, and let $\alpha>0$ be a constant. Denote by
$\mathcal{H}_{c}(\bm{\rho}):=\mathcal{H}(\rho_{1})+\mathcal{H}(\rho_{2}):=\int_{\mathbb{R}^{d}}[\rho_{1}\log\rho_{1}+\rho_{2}\log\rho_{2}]\,\mathrm{d}x$
the combined (negative) entropy functional for
$\rho\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$, and let
$[\rho]_{\alpha}:=\min\\{\rho,\alpha\\}$ denote a density $\rho$ cut off at
height $\alpha$. Then any minimizer $\bm{\rho}^{+}$ of the functional
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$ satisfies

$\frac{1}{2A}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}\mathopen{}\mathclose{{}\left([\rho_{1}^{+}]_{\alpha}}\right)}\right|^{2}+\mathopen{}\mathclose{{}\left|\nabla
F_{2}^{\prime}\mathopen{}\mathclose{{}\left([\rho_{2}^{+}]_{\alpha}}\right)}\right|^{2}}\right]\,\mathrm{d}x\leq\frac{\mathcal{H}_{c}(\hat{\bm{\rho}})-\mathcal{H}_{c}(\bm{\rho}^{+})}{\tau}+2dC_{K}\,,$
(37)

where $A:=\max_{j=1,2}\max_{r\leq\alpha}F_{j}^{\prime\prime}(r)<+\infty$.

###### Proof.

Let $\bm{\rho}^{+}$ be any minimizer of
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$ over
${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$. We prove (37) by analyzing the
behaviour of the functional
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$ along the heat flow
with initial datum $\bm{\rho}^{+}$: Denote by $\mathcal{K}^{s}$ the heat
kernel in ${\mathbb{R}^{d}}$, and let $\rho_{j}^{s}$ be the solution of the
heat equation $\partial_{s}\rho_{j}^{s}=\Delta\rho_{j}^{s}$ for
$(s,x)\in[0,\infty)\times{\mathbb{R}^{d}}$ with initial condition
$\rho_{j}^{0}=\rho_{j}^{+}$ for $j=1,2$. Explicitly, we have

$\displaystyle\mathcal{K}^{s}(x)=\frac{1}{(4\pi
s)^{d/2}}e^{-\frac{|x|^{2}}{4s}},\qquad\rho_{j}^{s}=\mathcal{K}^{s}\ast\rho_{j}^{+}.$

Clearly,
$\bm{\rho}^{s}=(\rho_{1}^{s},\rho_{2}^{s})\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$
is again a feasible pair of probability distributions, hence by minimality of
$\bm{\rho}^{+}$, it holds

$0\leq\liminf_{s\to
0^{+}}\frac{\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}^{s}|\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}^{+}|\hat{\bm{\rho}})}{s}\leq\limsup_{s\to
0^{+}}\frac{\mathrm{d}(\bm{\rho}^{s},\hat{\bm{\rho}})^{2}-\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}})^{2}}{2\tau
s}+\liminf_{s\to
0^{+}}\frac{\mathcal{E}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}{s}.$
(38)

It is well known that the heat flow is the metric gradient flow in the
Wasserstein space of the (negative) entropy functional
$\mathcal{H}(\rho)=\int_{\mathbb{R}^{d}}\rho\log\rho\,\mathrm{d}x$, which is
convex along generalized geodesics. By [2, Theorem 4.0.4], it thus satisfies
the evolutional variational inequality

$\frac{1}{2}\limsup_{s\to
0^{+}}\,\frac{{\mathbf{W}}(\rho_{j}^{s},\hat{\rho}_{j})^{2}-{\mathbf{W}}(\rho_{j}^{+},\hat{\rho}_{j})^{2}}{s}\leq\mathcal{H}(\hat{\rho}_{j})-\mathcal{H}(\rho_{j}^{+}).$

By combining the two components $j=1,2$ and dividing by $\tau$, this yields

$\limsup_{s\to
0^{+}}\,\frac{\mathrm{d}(\bm{\rho}^{s},\hat{\bm{\rho}})^{2}-\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}})^{2}}{2\tau
s}\leq\frac{\mathcal{H}_{c}(\hat{\bm{\rho}})-\mathcal{H}_{c}(\bm{\rho}^{+})}{\tau}.$
(39)

We now analyze the change in the energy term. We split it in the following
way:

$\displaystyle\mathcal{E}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})=\int_{\mathbb{R}^{d}}[\mathcal{F}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{F}_{\varepsilon}(\bm{\rho}^{+})]\,\mathrm{d}x+\int_{\mathbb{R}^{d}}[\rho_{1}^{s}K\ast\rho_{2}^{s}-\rho_{1}^{+}K\ast\rho_{2}^{+}]\,\mathrm{d}x.$

For the interaction term, we use the representation
$\rho_{j}^{s}=\mathcal{K}^{s}\ast\rho_{j}^{+}$ together with the symmetry and
semigroup property of the heat kernel $\mathcal{K}^{s}$ and the properties of
convolutions to obtain

$\displaystyle\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\rho_{1}^{s}K\ast\rho_{2}^{s}-\rho_{1}^{+}K\ast\rho_{2}^{+}}\right]\,\mathrm{d}x=\int_{\mathbb{R}^{d}}[(\mathcal{K}^{s}\ast\rho_{1}^{+})K\ast(\mathcal{K}^{s}\ast\rho_{2}^{+})-\rho_{1}^{+}K\ast\rho_{2}^{+}]\,\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\rho_{1}^{+}K\ast\mathcal{K}^{s}\ast\mathcal{K}^{s}\ast\rho_{2}^{+}-\rho_{1}^{+}K\ast\rho_{2}^{+}}\right]\,\mathrm{d}x=\int_{\mathbb{R}^{d}}\rho_{1}^{+}[K\ast\mathcal{K}^{2s}-K]\ast\rho_{2}^{+}\,\mathrm{d}x.$

Hence by dominated convergence and the equation
$\partial_{s}[K\ast\mathcal{K}^{s}]=\Delta[K\ast\mathcal{K}^{s}]$, it holds

$\begin{split}\lim_{s\to
0}\frac{1}{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\rho_{1}^{s}K\ast\rho_{2}^{s}-\rho_{1}^{+}K\ast\rho_{2}^{+}}\right]\,\mathrm{d}x&=\int_{\mathbb{R}^{d}}\rho_{1}^{+}\mathopen{}\mathclose{{}\left[\lim_{s\to
0}\frac{K\ast\mathcal{K}^{2s}-K}{s}}\right]\ast\rho_{2}^{+}\,\mathrm{d}x\\\
&=2\int_{\mathbb{R}^{d}}\rho_{1}^{+}\Delta K\ast\rho_{2}^{+}\,\mathrm{d}x\leq
2dC_{K}.\end{split}$ (40)

The dominated convergence is justified, since from the assumption that $K$ is
convex with $\|\nabla^{2}K\|\leq C_{K}$ and thus $0\leq\Delta K\leq dC_{K}$
pointwise in ${\mathbb{R}^{d}}$, it follows
$0\leq\Delta[K\ast\mathcal{K}^{s}]\leq dC_{K}$ pointwise in ${\mathbb{R}^{d}}$
for every $s\geq 0$, which implies that the difference quotients are globally
bounded by a constant.

For the $\mathcal{F}_{\varepsilon}$-term, we use that
$\mathcal{F}_{2\varepsilon}$ is convex by Proposition 2.4, and thus
$s\mapsto\int_{\mathbb{R}^{d}}\mathcal{F}_{2\varepsilon}(\bm{\rho}^{s})\,\mathrm{d}x$
is non-increasing along the heat flow $\bm{\rho}^{s}$: By Jensen and Fubini,
we have

$\begin{split}\int_{\mathbb{R}^{d}}\mathcal{F}_{2\varepsilon}(\bm{\rho}^{s})\,\mathrm{d}x&=\int_{\mathbb{R}^{d}}\mathcal{F}_{2\varepsilon}(\mathcal{K}^{s}\ast\bm{\rho}^{+})\,\mathrm{d}x=\int_{\mathbb{R}^{d}}\mathcal{F}_{2\varepsilon}\mathopen{}\mathclose{{}\left(\int_{\mathbb{R}^{d}}\bm{\rho}^{+}(y)\mathcal{K}^{s}(x-y)\,\mathrm{d}y}\right)\,\mathrm{d}x\\\
&\leq\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\mathcal{F}_{2\varepsilon}(\bm{\rho}^{+}(y))\mathcal{K}^{s}(x-y)\,\mathrm{d}x\mathrm{d}y=\int_{\mathbb{R}^{d}}\mathcal{F}_{2\varepsilon}(\bm{\rho}^{+})\,\mathrm{d}y.\end{split}$
(41)

Together with the identity
$\mathcal{F}_{\varepsilon}({\mathbf{r}})=\frac{1}{2}\mathcal{F}_{2\varepsilon}({\mathbf{r}})+\frac{1}{2}F_{1}(r_{1})+\frac{1}{2}F_{2}(r_{2})$,
this implies

$\displaystyle\int_{\mathbb{R}^{d}}[\mathcal{F}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{F}_{\varepsilon}(\bm{\rho}^{+})]\,\mathrm{d}x\leq\frac{1}{2}\int_{\mathbb{R}^{d}}[F_{1}(\rho_{1}^{s})-F_{1}(\rho_{1}^{+})]\,\mathrm{d}x+\frac{1}{2}\int_{\mathbb{R}^{d}}[F_{2}(\rho_{2}^{s})-F_{2}(\rho_{2}^{+})]\,\mathrm{d}x.$

We estimate the first integral (the second one is done in the same way) by
replacing $\rho_{1}^{+}$ with $\rho_{1}^{s_{0}}$ for some arbitrarily small
$s_{0}>0$ and then using a limiting argument. We use the fact that the
$\rho_{j}^{\sigma}$ solve the heat equation and Fubini to obtain for any
$s_{0}>0$:

$\frac{1}{2}\int_{\mathbb{R}^{d}}[F_{1}(\rho_{1}^{s})-F_{1}(\rho_{1}^{s_{0}})]\,\mathrm{d}x=\frac{1}{2}\int_{\mathbb{R}^{d}}\int_{s_{0}}^{s}F_{1}^{\prime}(\rho_{1}^{\sigma})\Delta\rho_{1}^{\sigma}\,\mathrm{d}\sigma\mathrm{d}x=\frac{1}{2}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}F_{1}^{\prime}(\rho_{1}^{\sigma})\Delta\rho_{1}^{\sigma}\,\mathrm{d}x\mathrm{d}\sigma.$
(42)

To justify the Fubini, it suffices to show that
$F_{1}^{\prime}(\rho_{1}^{\sigma})$ is uniformly bounded in
$L^{\infty}({\mathbb{R}^{d}})$, and that $\Delta\rho_{1}^{\sigma}$ is
uniformly bounded in $L^{1}({\mathbb{R}^{d}})$ for $\sigma\in[s_{0},s]$. The
first claim follows from global uniform boundedness of
$\rho_{1}^{\sigma}=\mathcal{K}^{\sigma}\ast\rho_{1}^{+}$, which holds since
$\mathcal{K}^{\sigma}(x)\leq\mathcal{K}^{s_{0}}(0)<+\infty$ for all
$x\in{\mathbb{R}^{d}}$ and all $\sigma\in[s_{0},s]$. The second claim follows
in a similar way from the fact that $\Delta\mathcal{K}^{\sigma}$ is uniformly
bounded in $L^{1}({\mathbb{R}^{d}})$ for $\sigma\in[s_{0},s]$.

By integrating (42) by parts, using that
$\nabla[\rho_{1}^{\sigma}]_{\alpha}=\nabla\rho_{1}^{\sigma}$ in
$\\{\rho_{1}^{\sigma}\leq\alpha\\}$ and $0$ in
$\\{\rho_{1}^{\sigma}>\alpha\\}$ for every $\sigma\geq s_{0}$ by regularity,
and that $0\leq F_{1}^{\prime\prime}([\rho_{1}^{\sigma}]_{\alpha})\leq A$ in
${\mathbb{R}^{d}}$, we obtain

$\begin{split}&\frac{1}{2}\int_{\mathbb{R}^{d}}[F_{1}(\rho_{1}^{s})-F_{1}(\rho_{1}^{s_{0}})]\,\mathrm{d}x=\frac{1}{2}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}F_{1}^{\prime}(\rho_{1}^{\sigma})\Delta\rho_{1}^{\sigma}\,\mathrm{d}x\mathrm{d}\sigma\\\
&=-\frac{1}{2}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}\nabla[F_{1}^{\prime}(\rho_{1}^{\sigma})]\cdot\nabla\rho_{1}^{\sigma}\,\mathrm{d}x\mathrm{d}\sigma=-\frac{1}{2}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}F_{1}^{\prime\prime}(\rho_{1}^{\sigma})\mathopen{}\mathclose{{}\left|\nabla\rho_{1}^{\sigma}}\right|^{2}\,\mathrm{d}x\mathrm{d}\sigma\\\
&\leq-\frac{1}{2A}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}F_{1}^{\prime\prime}\mathopen{}\mathclose{{}\left([\rho_{1}^{\sigma}]_{\alpha}}\right)^{2}\mathopen{}\mathclose{{}\left|\nabla[\rho_{1}^{\sigma}]_{\alpha}}\right|^{2}\,\mathrm{d}x\mathrm{d}\sigma=-\frac{1}{2A}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}([\rho_{1}^{\sigma}]_{\alpha})}\right|^{2}\,\mathrm{d}x\mathrm{d}\sigma.\end{split}$
(43)

We shall now pass to the limit $s_{0}\to 0$ in this estimate. On the right-
hand side, the integrand is non-negative, hence by monotone convergence, it
holds

$-\frac{1}{2A}\int_{s_{0}}^{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}([\rho_{1}^{\sigma}]_{\alpha})}\right|^{2}\,\mathrm{d}x\mathrm{d}\sigma\to-\frac{1}{2A}\int_{0}^{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}([\rho_{1}^{\sigma}]_{\alpha})}\right|^{2}\,\mathrm{d}x\mathrm{d}\sigma.$

On the left-hand side, we observe that by convexity of $F_{1}$, the map
$s\mapsto\int_{\mathbb{R}^{d}}F_{1}(\rho_{1}^{s})\,\mathrm{d}x$ is non-
increasing, which follows from Jensen’s inequality as in (41). In particular,

$\limsup_{s_{0}\to
0}\int_{\mathbb{R}^{d}}F_{1}(\rho_{1}^{s_{0}})\,\mathrm{d}x\leq\int_{\mathbb{R}^{d}}F_{1}(\rho_{1}^{+})\,\mathrm{d}x.$

Passing to the limit $s_{0}\to 0$ in (43) thus yields the estimate

$\displaystyle\frac{1}{2}\int_{\mathbb{R}^{d}}[F_{1}(\rho_{1}^{s})-F_{1}(\rho_{1}^{*})]\,\mathrm{d}x\leq-\frac{1}{2A}\int_{0}^{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}([\rho_{1}^{\sigma}]_{\alpha})}\right|^{2}\,\mathrm{d}x\mathrm{d}\sigma.$

By adding the term for the other component and dividing by $-s<0$, this yields

$\frac{1}{2As}\int_{0}^{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}([\rho_{1}^{\sigma}]_{\alpha})}\right|^{2}+\mathopen{}\mathclose{{}\left|\nabla
F_{2}^{\prime}([\rho_{2}^{\sigma}]_{\alpha})}\right|^{2}}\right]\,\mathrm{d}x\mathrm{d}\sigma\leq-\frac{1}{s}\int_{\mathbb{R}^{d}}[\mathcal{F}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{F}_{\varepsilon}(\bm{\rho}^{+})]\,\mathrm{d}x$
(44)

By inserting (39) and (40) into (38), we can control the right hand side of
(44) from above as $s\to 0$:

$\begin{split}\limsup_{s\to
0^{+}}\mathopen{}\mathclose{{}\left[-\frac{1}{s}\int_{\mathbb{R}^{d}}[\mathcal{F}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{F}_{\varepsilon}(\bm{\rho}^{+})]\,\mathrm{d}x}\right]&=-\liminf_{s\to
0^{+}}\frac{1}{s}\int_{\mathbb{R}^{d}}[\mathcal{F}_{\varepsilon}(\bm{\rho}^{s})-\mathcal{F}_{\varepsilon}(\bm{\rho}^{+})]\,\mathrm{d}x\\\
&\leq\frac{\mathcal{H}_{c}(\hat{\bm{\rho}})-\mathcal{H}_{c}(\bm{\rho}^{+})}{\tau}+2dC_{K}\end{split}$
(45)

By (44), this implies that $\nabla F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})$ is
bounded in $L^{2}({\mathbb{R}^{d}};{\mathbb{R}^{d}})$, uniformly along a
subsequence $s\to 0$. Hence by Alaoglu’s theorem, there exist $v_{j}\in
L^{2}({\mathbb{R}^{d}};{\mathbb{R}^{d}})$ such that $\nabla
F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})\rightharpoonup v_{j}$ weakly in
$L^{2}$ globally up to subsequence. Additionally, since
$F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})$ is uniformly bounded in $L^{\infty}$
by $F_{j}^{\prime}(\alpha)$, we get that
$F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})$ is locally bounded in $H^{1}$,
uniformly in $s$. Applying again Alaoglu’s theorem gives $u_{j}\in
H^{1}_{loc}({\mathbb{R}^{d}})$ such that up to subsequence,
$F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})\rightharpoonup u_{j}$ weakly in
$H^{1}_{loc}$ as $s\to 0$. In particular, $\nabla
F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})\rightharpoonup\nabla u_{j}$ weakly in
$L^{2}_{loc}$, thus we can identify $v_{j}=\nabla u_{j}$. Also, we have
$F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})\to u_{j}$ strongly in $L^{2}$
locally, which implies pointwise a.e. convergence to $u_{j}$ along a further
subsequence. It thus holds $u_{j}=F_{j}^{\prime}([\rho_{j}^{+}]_{\alpha})$, as
$\rho_{j}^{s}\to\rho_{j}^{+}$ a.e. pointwise. This identification shows that
the whole sequence $\nabla F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})$ converges
weakly in $L^{2}({\mathbb{R}^{d}})$. All together, we have shown $\nabla
F_{j}^{\prime}([\rho_{j}^{s}]_{\alpha})\rightharpoonup\nabla
F_{j}^{\prime}([\rho_{j}^{+}]_{\alpha})$ weakly in $L^{2}({\mathbb{R}^{d}})$
as $s\to 0$.

By lower semicontinuity of the $L^{2}$-norm with respect to weak convergence,
this implies

$\displaystyle\frac{1}{2A}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}\mathopen{}\mathclose{{}\left([\rho_{1}^{+}]_{\alpha}}\right)}\right|^{2}+\mathopen{}\mathclose{{}\left|\nabla
F_{2}^{\prime}\mathopen{}\mathclose{{}\left([\rho_{2}^{+}]_{\alpha}}\right)}\right|^{2}}\right]\,\mathrm{d}x\leq\liminf_{s\to
0^{+}}\frac{1}{2A}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}\mathopen{}\mathclose{{}\left([\rho_{1}^{s}]_{\alpha}}\right)}\right|^{2}+\mathopen{}\mathclose{{}\left|\nabla
F_{2}^{\prime}\mathopen{}\mathclose{{}\left([\rho_{2}^{s}]_{\alpha}}\right)}\right|^{2}}\right]\,\mathrm{d}x$
$\displaystyle\leq\limsup_{s\to
0^{+}}\frac{1}{2As}\int_{0}^{s}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\nabla
F_{1}^{\prime}([\rho_{1}^{\sigma}]_{\alpha})}\right|^{2}+\mathopen{}\mathclose{{}\left|\nabla
F_{2}^{\prime}([\rho_{2}^{\sigma}]_{\alpha})}\right|^{2}}\right]\,\mathrm{d}x\mathrm{d}\sigma\leq\frac{\mathcal{H}_{c}(\hat{\bm{\rho}})-\mathcal{H}_{c}(\bm{\rho}^{+})}{\tau}+2dC_{K},$

where the last estimate follows from (44) and (45). This finishes the proof of
(37). ∎

###  4.2. Euler-Lagrange system for the minimizing movement step

As a next step, we derive an Euler-Lagrange system that minimizers
$\bm{\rho}^{+}$ satisfy. For simplicity, we prove it under additional
assumptions on $\hat{\bm{\rho}}$, and later use approximation arguments to
generalize our results.

######  Lemma 4.3.

Assume in addition that both components of $\hat{\bm{\rho}}$ are absolutely
continuous with respect to the Lebesgue measure, and that there exist positive
radii $R_{j}>0$ such that $\hat{\rho}_{j}>0$ a.e. in $B_{R_{j}}$ and
$\hat{\rho}_{j}=0$ a.e. outside of $B_{R_{j}}$. Fix any minimizer
$\bm{\rho}^{+}$ of $\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$
over ${\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ obtained from Lemma 4.1. Then
there exist constants $C_{1},C_{2}$ and optimal pairs $(\varphi_{1},\psi_{1})$
and $(\varphi_{2},\psi_{2})$ of $c$-conjugate Kantorovich potentials for the
(dual) transport problems of $\rho_{1}^{+}$ to $\hat{\rho}_{1}$ and of
$\rho_{2}^{+}$ to $\hat{\rho}_{2}$, respectively, such that $\bm{\rho}^{+}$
satisfies the Euler-Lagrange system

$\begin{split}F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})&=(C_{1}-\frac{\varphi_{1}}{\tau}-K\ast\rho_{2}^{+})_{+}\\\
F_{2}^{\prime}(\rho_{2}^{+})+\varepsilon\partial_{r_{2}}h(\rho_{1}^{+},\rho_{2}^{+})&=(C_{2}-\frac{\varphi_{2}}{\tau}-K\ast\rho_{1}^{+})_{+}\,.\end{split}$
(46)

The optimal maps $T_{1},T_{2}:{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$,
transporting $\rho_{1}^{+}$ to $\hat{\rho}_{1}$ and $\rho_{2}^{+}$ to
$\hat{\rho}_{2}$, respectively, satisfy

$\begin{split}\frac{T_{1}(x)-x}{\tau}&=\nabla\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}}\right]\qquad\mathcal{L}^{d}\text{-a.e.
in }\\{\rho_{1}^{+}>0\\}\\\
\frac{T_{2}(x)-x}{\tau}&=\nabla\mathopen{}\mathclose{{}\left[F_{2}^{\prime}(\rho_{2}^{+})+\varepsilon\partial_{r_{2}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{1}^{+}}\right]\qquad\mathcal{L}^{d}\text{-a.e.
in }\\{\rho_{2}^{+}>0\\}.\end{split}$ (47)

###### Proof.

Take $\hat{\bm{\rho}}\in{\mathfrak{X}}_{0}$ and $\bm{\rho}^{+}$ as above. We
only need to show the first equation in (46) and (47), since the other one is
analogous in both cases.

Fix any $\eta_{0}\in L^{1}({\mathbb{R}^{d}})$ which is compactly supported,
bounded, satisfies $\int_{\mathbb{R}^{d}}\eta_{0}\,\mathrm{d}x=0$ and
$\rho_{1}^{s}:=\rho_{1}^{+}+s\eta_{0}\geq 0$ for every sufficiently small
$s>0$, and there exists a constant $H>0$ such that $\eta_{0}$ vanishes almost
everywhere on $\\{\rho_{1}^{+}>H\\}$, similarly as in the proof of Lemma 3.5.
Then for all sufficiently small $s>0$, it holds
$\rho_{1}^{s}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$, and thus by minimality
of $\rho_{1}^{+}$ for the functional
$\mathcal{E}_{\varepsilon,\tau}((\cdot,\rho_{2}^{+})|\hat{\bm{\rho}})$ over
${\mathcal{P}_{2}}({\mathbb{R}^{d}})$, we have with
$\bm{\rho}^{s}:=(\rho_{1}^{s},\rho_{2}^{+})$:

$\displaystyle 0$ $\displaystyle\leq\liminf_{s\to
0^{+}}\,\frac{\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}^{s}|\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}^{+}|\hat{\bm{\rho}})}{s}=\frac{1}{2\tau}\liminf_{s\to
0^{+}}\,\frac{{\mathbf{W}}(\rho_{1}^{s},\hat{\rho}_{1})^{2}-{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}}{s}+\mathopen{}\mathclose{{}\left.\frac{\,\mathrm{d}}{\,\mathrm{d}s}}\right|_{s=0^{+}}\mathcal{E}_{\varepsilon}(\bm{\rho}^{s})$
$\displaystyle=\frac{1}{2\tau}\liminf_{s\to
0^{+}}\,\frac{{\mathbf{W}}(\rho_{1}^{s},\hat{\rho}_{1})^{2}-{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}}{s}+\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}}\right]\,\eta_{0}\,\mathrm{d}x.$
(48)

The last equality follows as in the proof of Lemma 3.5 from the dominated
convergence theorem. Our goal is to construct a pair
$(\varphi_{1}^{0},\psi_{1}^{0})$ of optimal Kantorovich potentials for the
transport problem from $\rho_{1}^{+}$ to $\hat{\rho}_{1}$, such that

$\frac{1}{2}\liminf_{s\to
0^{+}}\,\frac{{\mathbf{W}}(\rho_{1}^{s},\hat{\rho}_{1})^{2}-{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}}{s}\leq\int_{\mathbb{R}^{d}}\varphi_{1}^{0}\,\eta_{0}\,\mathrm{d}x.$
(49)

Denote by $(\varphi_{1}^{s},\psi_{1}^{s})$ pairs of optimal Kantorovich
potentials for the transport from $\rho_{1}^{s}$ to $\hat{\rho}_{1}$. We can
fix some point $x_{0}\in\mathrm{supp}\,\rho_{1}^{+}$ and without loss of
generality assume $\varphi_{1}^{s}(x_{0})=0$ for every $s$. Since
$\\{\rho_{1}^{+}>0\\}\subset\\{\rho_{1}^{s}>0\\}$ for every $s\in[0,\bar{s}]$
with some $\bar{s}>0$, the pair $(\varphi_{1}^{s},\psi_{1}^{s})$ is admissible
for the (dual) transport problem from $\rho_{1}^{+}$ to $\hat{\rho}_{1}$.
Hence for every $s\in[0,\bar{s}]$

$\displaystyle\frac{1}{2}\mathopen{}\mathclose{{}\left[{\mathbf{W}}(\rho_{1}^{s},\hat{\rho}_{1})^{2}-{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}}\right]\leq\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\varphi_{1}^{s}\,\rho_{1}^{s}+\psi_{1}^{s}\,\hat{\rho}_{1}}\right]\,\mathrm{d}x-\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\varphi_{1}^{s}\,\rho_{1}^{+}+\psi_{1}^{s}\,\hat{\rho}_{1}}\right]\,\mathrm{d}x=s\int_{\mathbb{R}^{d}}\varphi_{1}^{s}\,\eta_{0}\,\mathrm{d}x.$

The proof of (49) thus reduces to the construction of an optimal potential
$\varphi_{1}^{0}$ for the transport problem from $\rho_{1}^{+}$ to
$\hat{\rho}_{1}$ such that

$\liminf_{s\to
0^{+}}\int_{\mathbb{R}^{d}}\varphi_{1}^{s}\,\eta_{0}\,\mathrm{d}x\leq\int_{\mathbb{R}^{d}}\varphi_{1}^{0}\,\eta_{0}\,\mathrm{d}x.$
(50)

Denote $\Phi_{1}^{s}:=\frac{1}{2}|\cdot|^{2}-\varphi_{1}^{s}$. Since
$(\varphi_{1}^{s},\psi_{1}^{s})$ is a pair of optimal $c$-conjugate
potentials, $\Phi_{1}^{s}$ is a convex function and satisfies
$|\nabla\Phi_{1}^{s}(x)|=|T^{s}(x)|$ for almost every
$x\in\\{\rho_{1}^{s}>0\\}$, where $T^{s}\\#\rho_{1}^{s}=\hat{\rho}_{1}$ is the
optimal transport map from $\rho_{1}^{s}$ to $\hat{\rho}_{1}$. Since by
assumption, $\hat{\rho}_{1}$ vanishes a.e. outside of $B_{R_{1}}$, this
implies $|\nabla\Phi_{1}^{s}|\leq R_{1}$ a.e. on $\\{\rho_{1}^{s}>0\\}$. Since
$\Phi_{1}^{s}$ is convex on ${\mathbb{R}^{d}}$, this implies existence of a
convex function $\tilde{\Phi}_{1}^{s}:{\mathbb{R}^{d}}\to{\mathbb{R}}$ with
$\tilde{\Phi}_{1}^{s}=\Phi_{1}^{s}$ in $\\{\rho_{1}^{s}>0\\}$ and
$|\nabla\tilde{\Phi}_{1}^{s}|\leq R_{1}$ in ${\mathbb{R}^{d}}$. Hence the
family of potentials $\\{\tilde{\Phi}_{1}^{s}\\}_{s\in[0,\bar{s}]}$ is equi-
Lipschitz in ${\mathbb{R}^{d}}$, and because
$\tilde{\Phi}_{1}^{s}(x_{0})=\Phi_{1}^{s}(x_{0})=0$ for all $s$, it is bounded
pointwise. The Arzelà-Ascoli Theorem thus yields existence of a locally
uniformly converging subsequence
$(\tilde{\Phi}_{1}^{s_{k}})_{k}\to\Phi_{1}^{0}$ in ${\mathbb{R}^{d}}$.

Clearly, $\varphi_{1}^{0}:=\frac{1}{2}|\cdot|^{2}-\Phi_{1}^{0}$ satisfies
(50). It remains to show that $\varphi_{1}^{0}$, together with its
$c$-conjugate $\psi_{1}^{0}$, forms an optimal pair for the transport from
$\rho_{1}^{+}$ to $\hat{\rho}_{1}$. Since $\Phi_{1}^{0}$ is convex, it is
sufficient to show $\nabla\Phi_{1}^{0}\\#\rho_{1}^{+}=\hat{\rho}_{1}$. This
however follows from the fact that the $\tilde{\Phi}_{1}^{s}$ are convex and
$\tilde{\Phi}_{1}^{s_{k}}\to\Phi_{1}^{0}$ locally uniformly, which implies
$\nabla\tilde{\Phi}_{1}^{s_{k}}\to\nabla\Phi_{1}^{0}$ almost everywhere: For
any $\zeta\in C^{0}_{c}({\mathbb{R}^{d}})$, it holds by dominated convergence

$\displaystyle\int_{\mathbb{R}^{d}}\zeta\,\nabla\Phi_{1}^{0}\\#\rho_{1}^{+}\,\mathrm{d}x=\int_{\mathbb{R}^{d}}\zeta(\nabla\Phi_{1}^{0}(x))\,\rho_{1}^{+}(x)\,\mathrm{d}x=\lim_{k}\int_{\mathbb{R}^{d}}\zeta(\nabla\tilde{\Phi}_{1}^{s_{k}}(x))\,\rho_{1}^{+}(x)\,\mathrm{d}x=\int_{\mathbb{R}^{d}}\zeta\,\hat{\rho}_{1}\,\mathrm{d}x.$

All together, we have proven the inequality (49) for the optimal potential
$\varphi_{1}^{0}$.

We can insert this result into (48) to obtain

$0\leq\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}+\frac{\varphi_{1}^{0}}{\tau}}\right]\,\eta_{0}\,\mathrm{d}x.$

In order to obtain the Euler-Lagrange system (46), we have to show that for
some optimal Kantorovich potential $\varphi_{1}$, this inequality holds true
for all $\eta$ which satisfy our assumptions at once, and not just for the
specific $\eta_{0}$ used above to construct the potentials $\varphi_{1}^{0}$
and $\psi_{1}^{0}$. We define, using the notations
$\Phi_{1}^{0}=\frac{1}{2}|\cdot|^{2}-\varphi_{1}^{0}$ and
$\Psi_{1}^{0}=\frac{1}{2}|\cdot|^{2}-\psi_{1}^{0}$, the functions

$\displaystyle\Psi_{1}:=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}\Psi_{1}^{0}&\text{
in }B_{R_{1}}\\\ +\infty&\text{
else}\end{array}}\right.,\quad\Phi_{1}:=(\hat{\Psi}_{1}^{0})^{*},$

where $(\cdot)^{*}$ denotes the convex conjugate. Since $\hat{\rho}_{1}$ is
supported in $B_{R_{1}}$, it follows that $\Psi_{1}$ is convex with
$\nabla\Psi_{1}\\#\hat{\rho}_{1}=\nabla\Psi_{1}^{0}\\#\hat{\rho}_{1}=\rho_{1}^{+}$,
implying that the pair of $c$-conjugate potentials
$(\varphi_{1},\psi_{1}):=(\frac{1}{2}|\cdot|^{2}-\Phi_{1},\frac{1}{2}|\cdot|^{2}-\Psi_{1})$
is optimal. We claim that for all $\eta$ that have the properties of
$\eta_{0}$ specified above, it holds

$0\leq\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}+\frac{\varphi_{1}}{\tau}}\right]\,\eta\,\mathrm{d}x.$
(51)

To prove this, take any such $\eta$. By the same construction as for
$\eta_{0}$, there exists a pair of optimal potentials
$(\tilde{\varphi}_{1},\tilde{\psi}_{1})$, such that (51) holds with
$\tilde{\varphi}_{1}$ instead of $\varphi_{1}$. Now, observe that
$\tilde{\Psi}_{1}:=\frac{1}{2}|\cdot|^{2}-\tilde{\psi}_{1}$ is a convex
function whose gradient agrees with the gradient of $\Psi_{1}^{0}$ a.e. on the
connected set $\\{\hat{\rho}_{1}>0\\}=B_{R_{1}}$, since by absolute continuity
of $\hat{\rho}_{1}$ with respect to the Lebesgue measure, the optimal
transport map from $\hat{\rho}_{1}$ to $\rho_{1}^{+}$ is unique and given by
the gradient of any one of the convex potentials $\Psi_{1}^{0}$ and
$\tilde{\Psi}_{1}$. Hence there exists a global constant $C$ such that
$\Psi_{1}^{0}=\tilde{\Psi}_{1}+C$ in $B_{R_{1}}$. By definition of $\Psi_{1}$,
this yields $\Psi_{1}\geq\tilde{\Psi}_{1}+C$ in ${\mathbb{R}^{d}}$ with
equality in $B_{R_{1}}$. By the properties of the convex conjugate, this
implies

$\displaystyle\frac{1}{2}|\cdot|^{2}-\varphi_{1}=\Phi_{1}=\Psi_{1}^{*}\leq\tilde{\Psi}_{1}^{*}-C=\frac{1}{2}|\cdot|^{2}-\tilde{\varphi}_{1}-C,$

or equivalently, $\varphi_{1}\geq\tilde{\varphi}_{1}+C$ in ${\mathbb{R}^{d}}$,
with equality almost everywhere on the image of the map $\nabla\Psi_{1}$ over
$B_{R_{1}}$, which is $\\{\rho_{1}^{+}>0\\}$ since $\nabla\Psi_{1}$ pushes
$\hat{\rho}_{1}$ to $\rho_{1}^{+}$. Using that $\eta\geq 0$ in
$\\{\rho_{1}^{+}=0\\}$ and $\int_{\mathbb{R}^{d}}\eta=0$, this gives

$\displaystyle\int_{\mathbb{R}^{d}}\varphi_{1}\eta\,\mathrm{d}x=\int_{\\{\rho_{1}^{+}>0\\}}\mathopen{}\mathclose{{}\left(\tilde{\varphi}_{1}+C}\right)\,\eta\,\mathrm{d}x+\int_{\\{\rho_{1}^{+}=0\\}}\varphi_{1}\eta\,\mathrm{d}x\geq\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(\tilde{\varphi}_{1}+C}\right)\,\eta\,\mathrm{d}x=\int_{\mathbb{R}^{d}}\tilde{\varphi}_{1}\eta\,\mathrm{d}x,$

proving (51), since we have already shown that it holds with
$\tilde{\varphi}_{1}$ instead of $\varphi_{1}$.

By the same argument as in the proof of Lemma 3.6, (51) yields the existence
of some constant $C_{1}$, such that everywhere in ${\mathbb{R}^{d}}$, it holds

$F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}+\frac{\varphi_{1}}{\tau}\geq
C_{1},$

with equality almost everywhere on $\\{\rho_{1}^{+}>0\\}$. The form (46) then
follows in the same way as Lemma 3.7 from the properties of
$\mathcal{F}_{\varepsilon}$. Hence the optimal potential $\varphi_{1}$
satisfies the first equation in (46). The equation for the other component is
derived analogously.

Since the left-hand side of (46) is strictly positive on
$\\{\rho_{j}^{+}>0\\}$, so is the right-hand side. The Kantorovich potentials
are differentiable almost everywhere in $\\{\rho_{j}^{+}>0\\}$ with
$-\nabla\varphi_{j}=T_{j}(x)-x$, hence the gradient of the right-hand side
exists in $\\{\rho_{j}^{+}>0\\}$, and the equations (47) follow from (46) by
taking the gradient on both sides. ∎

As a consequence of Lemma 4.3, we prove that under the same assumptions, the
transition from $\hat{\bm{\rho}}$ to $\bm{\rho}^{+}$ satisfies a time-discrete
version of the weak formulation of the evolution equations (1).

######  Corollary 4.4.

Let $\hat{\bm{\rho}}$ satisfy the assumptions of Lemma 4.3. Then for all
$\zeta\in C^{\infty}_{c}({\mathbb{R}^{d}})$, there hold the estimates

$\begin{split}\mathopen{}\mathclose{{}\left|\int_{\mathbb{R}^{d}}\frac{\rho_{1}^{+}-\hat{\rho}_{1}}{\tau}\,\zeta\,\mathrm{d}x+\int_{\mathbb{R}^{d}}\rho_{1}^{+}\,\nabla\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}}\right]\cdot\nabla\zeta\,\mathrm{d}x}\right|&\leq\frac{{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}}{2\tau}\,\|\zeta\|_{C^{2}}\\\
\mathopen{}\mathclose{{}\left|\int_{\mathbb{R}^{d}}\frac{\rho_{2}^{+}-\hat{\rho}_{2}}{\tau}\,\zeta\,\mathrm{d}x+\int_{\mathbb{R}^{d}}\rho_{2}^{+}\,\nabla\mathopen{}\mathclose{{}\left[F_{2}^{\prime}(\rho_{2}^{+})+\varepsilon\partial_{r_{2}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{1}^{+}}\right]\cdot\nabla\zeta\,\mathrm{d}x}\right|&\leq\frac{{\mathbf{W}}(\rho_{2}^{+},\hat{\rho}_{2})^{2}}{2\tau}\,\|\zeta\|_{C^{2}}.\end{split}$
(52)

Denoting by $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ the respective left-hand
sides of these estimates, it holds

$\mathcal{R}_{1}+\mathcal{R}_{2}\leq\mathopen{}\mathclose{{}\left(\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}\right)\,\|\zeta\|_{C^{2}}.$
(53)

###### Proof.

Let $T_{1}\\#\rho_{1}^{+}=\hat{\rho}_{1}$ be the optimal transport map, and
let $\zeta\in C^{\infty}_{c}({\mathbb{R}^{d}})$ be arbitrary. We obtain for
the difference quotient on the left-hand side

$\displaystyle\int_{\mathbb{R}^{d}}\frac{\rho_{1}^{+}-\hat{\rho}_{1}}{\tau}\,\zeta\,\mathrm{d}x$
$\displaystyle=\frac{1}{\tau}\mathopen{}\mathclose{{}\left[\int_{\mathbb{R}^{d}}\rho_{1}^{+}\,\zeta\,\mathrm{d}x-\int_{\mathbb{R}^{d}}T_{1}\\#\rho_{1}^{+}\,\zeta\,\mathrm{d}x}\right]=\int_{\mathbb{R}^{d}}\rho_{1}^{+}(x)\,\frac{\zeta(x)-\zeta\mathopen{}\mathclose{{}\left(T_{1}(x)}\right)}{\tau}\,\mathrm{d}x$
$\displaystyle=-\int_{\mathbb{R}^{d}}\rho_{1}^{+}(x)\,\mathopen{}\mathclose{{}\left[\nabla\zeta(x)\cdot\frac{T_{1}(x)-x}{\tau}+\frac{(T_{1}(x)-x)\cdot\nabla^{2}\zeta(\xi_{x})(T_{1}(x)-x)}{2\tau}}\right]\,\mathrm{d}x$

for some $\xi_{x}\in[x,T_{1}(x)]$. Applying (47) and
${\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}=\int_{\mathbb{R}^{d}}\rho_{1}^{+}(x)\,\mathopen{}\mathclose{{}\left|T_{1}(x)-x}\right|^{2}\,\mathrm{d}x$
thus yields

$\displaystyle\mathopen{}\mathclose{{}\left|\int_{\mathbb{R}^{d}}\frac{\rho_{1}^{+}-\hat{\rho}_{1}}{\tau}\,\zeta\,\mathrm{d}x+\int_{\mathbb{R}^{d}}\rho_{1}^{+}\,\nabla\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}}\right]\cdot\nabla\zeta\,\mathrm{d}x}\right|$
$\displaystyle=\mathopen{}\mathclose{{}\left|\int_{\mathbb{R}^{d}}\frac{\rho_{1}^{+}-\hat{\rho}_{1}}{\tau}\,\zeta\,\mathrm{d}x+\int_{\mathbb{R}^{d}}\rho_{1}^{+}(x)\,\nabla\zeta(x)\cdot\frac{T_{1}(x)-x}{\tau}\,\mathrm{d}x}\right|$
$\displaystyle\leq\int_{\mathbb{R}^{d}}\rho_{1}^{+}(x)\,\mathopen{}\mathclose{{}\left|\frac{(T_{1}(x)-x)\cdot\nabla^{2}\zeta(\xi_{x})(T_{1}(x)-x)}{2\tau}}\right|\,\mathrm{d}x\leq\frac{{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}}{2\tau}\,\|\zeta\|_{C^{2}}.$

This proves the first inequality in (52), the second one follows analogously.
By adding the two inequalities and applying (36), we obtain (53). ∎

###  4.3. Regularity estimates for the minimizers

In order to prove convergence of the minimizing movement scheme to a weak
solution of (1), we need to prove the estimates from Corollary 4.4 for more
general $\hat{\bm{\rho}}$, which will be done using an approximating argument.
For the convergence proof, we require further regularity results for the
minimizers $\bm{\rho}^{+}$, in particular, a global $H^{1}$-bound for
$F_{j}^{\prime}(\rho_{j}^{+})$. This will follow from the next results,
together with the bound on $\nabla F_{j}^{\prime}(\rho_{j}^{+})$ from Lemma
4.2.

######  Lemma 4.5.

Let $\hat{\bm{\rho}}$ satisfy the assumptions of Lemma 4.3. Then it holds

$\sum_{j=1,2}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla[F_{j}^{\prime}(\rho_{j}^{+})+\varepsilon\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})]}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x\leq
4\,\frac{\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}{\tau}+\frac{8C_{K}^{2}}{\lambda}\,\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})\,.$
(54)

###### Proof.

Since $K\ast\rho_{j}^{+}$ is differentiable everywhere, the Euler-Lagrange
equations in the form (47) imply that the gradients on the left-hand side are
well-defined a.e. in $\\{\rho_{j}^{+}>0\\}$. For $j=1$ (the other term is
analogous), it holds

$\displaystyle\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})]}\right|^{2}\,\rho_{1}^{+}\,\mathrm{d}x$
$\displaystyle\leq
2\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}]}\right|^{2}\,\rho_{1}^{+}\,\mathrm{d}x+2\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla[K\ast\rho_{2}^{+}]}\right|^{2}\,\rho_{1}^{+}\,\mathrm{d}x$
$\displaystyle\leq\frac{2}{\tau^{2}}\int_{\mathbb{R}^{d}}|T_{1}(x)-x|^{2}\,\rho_{1}^{+}\,\mathrm{d}x+\frac{4C_{K}^{2}}{\lambda}\int_{\mathbb{R}^{d}}\rho_{1}^{+}K\ast\rho_{2}^{+}\,\mathrm{d}x\leq\frac{2}{\tau^{2}}\,{\mathbf{W}}(\rho_{1}^{+},\hat{\rho}_{1})^{2}+\frac{4C_{K}^{2}}{\lambda}\,\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}),$

where the second to last inequality uses (47) for the first integral, and the
fact that $K$, and thus $K\ast\rho_{2}^{+}$, is non-negative and
$\lambda$-convex with its second derivative bounded by $C_{K}$ for the second
integral. Adding the analogous term for $j=2$ yields

$\displaystyle\sum_{j=1,2}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla[F_{j}^{\prime}(\rho_{j}^{+})+\varepsilon\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})]}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x\leq\frac{2}{\tau^{2}}\,\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}})^{2}+\frac{8C_{K}^{2}}{\lambda}\,\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}),$

which together with (36) yields the claimed estimate (54). ∎

In order to obtain a global $H^{1}$-estimate for the functions
$F_{j}^{\prime}(\rho_{j}^{+})$, we first need to prove existence of their
gradients. This follows directly from the following lemma, which is a
consequence of the Euler-Lagrange system (46).

######  Lemma 4.6.

Under the assumptions of Lemma 4.3, the $F_{j}^{\prime}(\rho_{j}^{+})$ are
locally Lipschitz in ${\mathbb{R}^{d}}$.

###### Proof.

This follows from (46) in the same way as in Lemma 3.11, since the right-hand
side of (46) is locally Lipschitz in ${\mathbb{R}^{d}}$. Note that the
Kantorovich potentials $\varphi_{j}$ are of the form
$\frac{1}{2}|\cdot|^{2}-\Phi_{j}$, where $\Phi_{j}$ are convex functions which
are finite in ${\mathbb{R}^{d}}$, since the left-hand side of (46) is finite
everywhere. ∎

With this regularity result, we can prove the following estimate for $\nabla
F_{j}^{\prime}(\rho_{j}^{+})$:

######  Lemma 4.7.

Let $\hat{\bm{\rho}}$ satisfy the assumptions of Lemma 4.3, and assume that
$\varepsilon$ is such that $3\kappa^{2}\varepsilon^{2}<1$, where
$\kappa^{2}:=\max\\{\kappa_{1,1},\kappa_{1,2}\\}^{2}+\max\\{\kappa_{2,1},\kappa_{2,2}\\}^{2}$
with the constants $\kappa_{j,i}$ from Hypothesis 2.2. Then there exists a
constant $B$, independent of $\hat{\bm{\rho}}$ and $\bm{\rho}^{+}$, such that

$\sum_{j=1,2}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x\leq
B\,\mathopen{}\mathclose{{}\left(\frac{\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}{\tau}+\frac{2C_{K}^{2}}{\lambda}\,\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}\right).$

###### Proof.

The gradient of $F_{j}^{\prime}(\rho_{j}^{+})$ exists almost everywhere by
local Lipschitz-continuity, and the same is true for
$F_{j}^{\prime}(\rho_{j}^{+})+\varepsilon\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})$,
see the proof of Lemma 4.5. Thus, the gradient of
$\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})$ is also well-defined in
${\mathbb{R}^{d}}$. The claim will now follow from Lemma 4.5 and Hypothesis
2.2: With
$\bm{u}^{+}:=(F_{1}^{\prime}(\rho_{1}^{+}),F_{2}^{\prime}(\rho_{2}^{+}))$, we
have

$\displaystyle\nabla\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})=\nabla\theta_{j}(\bm{u}^{+})=\theta_{j,1}(\bm{u}^{+})\nabla
F_{1}^{\prime}(\rho_{1}^{+})+\theta_{j,2}(\bm{u}^{+})\nabla
F_{2}^{\prime}(\rho_{2}^{+}).$

By splitting the integral, we thus obtain for $j=1,2$ the estimate

$\displaystyle\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla\mathopen{}\mathclose{{}\left[F_{j}^{\prime}(\rho_{j}^{+})+\varepsilon\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})}\right]-\varepsilon\sum_{i=1,2}\theta_{j,i}(\bm{u}^{+})\nabla
F_{i}^{\prime}(\rho_{i}^{+})}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x$
$\displaystyle\leq
3\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla[F_{j}^{\prime}(\rho_{j}^{+})+\varepsilon\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})]}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x+3\varepsilon^{2}\sum_{i=1,2}\int_{\mathbb{R}^{d}}\theta^{2}_{j,i}(\bm{u}^{+})\mathopen{}\mathclose{{}\left|\nabla
F_{i}^{\prime}(\rho_{i}^{+})}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x.$

Summing over $j=1,2$ and applying (54) for the first term and (15) for the
other two yields

$\displaystyle\sum_{j=1,2}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x$
$\displaystyle\leq
12\,\frac{\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}{\tau}+\frac{24C_{K}^{2}}{\lambda}\,\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})$
$\displaystyle+3\kappa^{2}\varepsilon^{2}\sum_{i=1,2}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{i}^{\prime}(\rho_{i}^{+})}\right|^{2}\,\rho_{i}^{+}\,\mathrm{d}x.$

Subtracting the last term and using $3\kappa^{2}\varepsilon^{2}<1$ yields the
claim with $B=\frac{12}{1-3\kappa^{2}\varepsilon^{2}}$. ∎

What makes Lemma 4.7 important is that it controls $|\nabla
F_{j}^{\prime}(\rho_{j}^{+})|^{2}$ on the set where $\rho_{j}^{+}$ is large.
Together with Lemma 4.2, which controls the same gradient when $\rho_{j}^{+}$
is small, this implies an explicit global $L^{2}$-bound for $\nabla
F_{j}(\rho_{j}^{+})$:

######  Corollary 4.8.

Let $\hat{\bm{\rho}}$ satisfy the assumptions of Lemma 4.3. Let $\alpha>0$ be
arbitrary, and denote by $A$ and $B$ the constants from Lemmas 4.2 and 4.7,
respectively. Then it holds

$\sum_{j=1,2}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\mathrm{d}x\leq
C\,\mathopen{}\mathclose{{}\left(\frac{\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}{\tau}+\frac{2C_{K}^{2}}{\lambda}\,\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})+\frac{\mathcal{H}_{c}(\hat{\bm{\rho}})-\mathcal{H}_{c}(\bm{\rho}^{+})}{\tau}+2dC_{K}}\right),$
(55)

with the constant $C:=\max\\{\frac{B}{\alpha},2A\\}$.

###### Proof.

For $j=1,2$, we split the integral into one part where $\rho_{j}^{+}$ is large
and one where it is small. By local Lipschitz-continuity, it holds $\nabla
F_{j}^{\prime}(\rho_{j}^{+})=\nabla F_{j}^{\prime}([\rho_{j}^{+}]_{\alpha})$
on $\\{\rho_{j}^{+}<\alpha\\}$, hence

$\displaystyle\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\mathrm{d}x$
$\displaystyle=\int_{\\{\rho_{j}^{+}\geq\alpha\\}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\mathrm{d}x+\int_{\\{\rho_{j}^{+}<\alpha\\}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\mathrm{d}x$
$\displaystyle\leq\frac{1}{\alpha}\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}(\rho_{j}^{+})}\right|^{2}\,\rho_{j}^{+}\,\mathrm{d}x+\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left|\nabla
F_{j}^{\prime}([\rho_{j}^{+}]_{\alpha})}\right|^{2}\,\mathrm{d}x.$

The claim now follows by summing over $j=1,2$ and applying Lemmas 4.7 and 4.2.
∎

###  4.4. Generalization of the time-discrete error estimates

With the help of the estimate from Corollary 4.8, we can now prove the time-
discrete error estimate (53) for much more general $\hat{\bm{\rho}}$.
Concretely, we prove the following theorem:

######  Theorem 4.9.

Let $\hat{\bm{\rho}}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ be such that
$\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})<+\infty$. Then there exists a
minimizer $\bm{\rho}^{+}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$ of the
functional $\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$, such that
for all $\zeta\in C^{\infty}_{c}({\mathbb{R}^{d}})$:

$\mathcal{R}_{1}+\mathcal{R}_{2}\leq\mathopen{}\mathclose{{}\left(\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})}\right)\,\|\zeta\|_{C^{2}}\,,$

where $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are the integral errors defined
as in Corollary 4.4.

###### Proof.

In Corollary 4.4, the claim has already been shown under the more restrictive
assumptions that both components of $\hat{\bm{\rho}}$ are strictly positive
almost everywhere inside balls $B_{R_{j}}$ and vanish almost everywhere else.
We prove the more general case by approximating $\hat{\bm{\rho}}$ by a
sequence of pairs that satisfy these stronger assumptions.

Let $\hat{\bm{\rho}}$ be as above. Without loss of generality, we can assume
$\hat{\bm{\rho}}\in{\mathfrak{X}}_{0}$. We start by showing that the
assumption $\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})<+\infty$ implies that
$\mathcal{H}_{c}(\hat{\bm{\rho}})<+\infty$ as well. To see this, observe that
by comparing the asymptotic growths of the derivatives, it holds $r\log r\leq
cF_{j}(r)$ for every $r\geq 0$ with some constant $c$. It thus holds for every
$\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$

$\displaystyle\mathcal{H}_{c}(\bm{\rho})=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[\rho_{1}\log\rho_{1}+\rho_{2}\log\rho_{2}}\right]\,\mathrm{d}x\leq
c\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left[F_{1}(\rho_{1})+F_{2}(\rho_{2})}\right]\,\mathrm{d}x\leq
2c\,\mathcal{E}_{\varepsilon}(\bm{\rho}),$

proving that $\mathcal{H}_{c}(\hat{\bm{\rho}})$ is finite.

For $j=1,2$, we now take a sequence of positive radii
$0<R_{j,1}<R_{j,2}<\dots$ such that $R_{j,n}\nearrow$
$R_{j,\infty}:=\sup\,\\{R>0:\mathrm{supp}\,\hat{\rho}_{j}\nsubseteq
B_{R}\\}\in(0,+\infty]$. Further, we fix constants $H_{j},C_{j}>0$ such that
for every $n$, it holds

$|B_{R_{j,n}}^{\prime}|\geq C_{j}\,,\quad\text{where
}B_{R_{j,n}}^{\prime}:=\\{\hat{\rho}_{j}\leq H_{j}\\}\cap B_{R_{j,n}}.$

Define for every $n$ the probability densities

$\hat{\rho}_{j,n}:=\mathds{1}_{B_{R_{j,n}}}\,\hat{\rho}_{j}+\frac{m_{j,n}}{|B_{R_{j,n}}^{\prime}|}\mathds{1}_{B_{R_{j,n}}^{\prime}}\,,\quad\text{with
}m_{j,n}:=1-\int_{B_{R_{j,n}}}\hat{\rho}_{j}\,\mathrm{d}x\in(0,1].$

For the densities $\hat{\rho}_{j,n}$, we have
$\\{\hat{\rho}_{j,n}>0\\}=B_{R_{j,n}}$. Hence, the pairs
$\hat{\bm{\rho}}_{n}:=(\hat{\rho}_{1,n},\hat{\rho}_{2,n})$ satisfy the
assumptions of Corollary 4.4, and thus for every minimizer $\bm{\rho}_{n}^{+}$
of $\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\rho}_{n})$, the estimate (53).
The same estimate for $\hat{\bm{\rho}}$, which is the claim, thus follows if
we prove for $n\to\infty$ the convergences

$\int_{\mathbb{R}^{d}}\frac{\rho_{1,n}^{+}-\hat{\rho}_{1,n}}{\tau}\,\zeta\,\mathrm{d}x\rightarrow\int_{\mathbb{R}^{d}}\frac{\rho_{1}^{+}-\hat{\rho}_{1}}{\tau}\,\zeta\,\mathrm{d}x,$
(56)
$\begin{split}&\int_{\mathbb{R}^{d}}\rho_{1,n}^{+}\,\nabla\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1,n}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})+K\ast\rho_{2,n}^{+}}\right]\cdot\nabla\zeta\,\mathrm{d}x\\\
&\rightarrow\int_{\mathbb{R}^{d}}\rho_{1}^{+}\,\nabla\mathopen{}\mathclose{{}\left[F_{1}^{\prime}(\rho_{1}^{+})+\varepsilon\partial_{r_{1}}h(\rho_{1}^{+},\rho_{2}^{+})+K\ast\rho_{2}^{+}}\right]\cdot\nabla\zeta\,\mathrm{d}x\,,\end{split}$
(57)

$\limsup_{n\to\infty}\,\mathopen{}\mathclose{{}\left(\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}_{n})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}_{n})}\right)\leq\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})-\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})$
(58)

for some minimizer $\bm{\rho}^{+}$ of
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$, up to subsequence.
The convergence of the other component’s terms will follow analogously.

We start by observing that by construction,
$\hat{\rho}_{j,n}\to\hat{\rho}_{j}$ in $L^{1}({\mathbb{R}^{d}})$, and with
respect to the ${\mathbf{W}}$-metric. Additionally, it holds that
$\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}_{n})\to\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})$
and $\mathcal{H}_{c}(\hat{\bm{\rho}}_{n})\to\mathcal{H}_{c}(\hat{\bm{\rho}})$,
which can be shown from the construction by dominated convergence. In
particular, since $\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})$ and
$\mathcal{H}_{c}(\hat{\bm{\rho}})$ are finite,
$(\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}_{n}))_{n}$ and
$(\mathcal{H}_{c}(\hat{\bm{\rho}}_{n}))_{n}$ are bounded sequences. Since it
holds
$\mathcal{E}_{\varepsilon}(\bm{\rho}_{n}^{+})\leq\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}_{n})$,
the sequence $(\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}_{n}))_{n}$ is thus
bounded as well. By Lemma 4.1, it holds
${\mathfrak{m}}_{1}[\rho_{1}^{+}]+{\mathfrak{m}}_{1}[\rho_{2}^{+}]={\mathfrak{m}}_{1}[\hat{\rho}_{1}^{n}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}^{n}]\to{\mathfrak{m}}_{1}[\hat{\rho}_{1}]+{\mathfrak{m}}_{1}[\hat{\rho}_{2}]=0$.
Lemma 3.1 thus yields
$\bm{\rho}_{n}^{+}\rightharpoonup\bm{\rho}^{+}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})^{2}$
weakly in $L^{1}$ up to subsequence. We check that $\bm{\rho}^{+}$ is a global
minimizer of the functional
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$. Indeed, it holds

$\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})\leq\liminf_{n\to\infty}\mathcal{E}_{\varepsilon}(\bm{\rho}_{n}^{+})$
(59)

by weak $L^{1}$-convergence of the $\bm{\rho}_{n}^{+}$ and Lemma 3.2. By lower
semicontinuity of the distance with respect to weak $L^{1}$-convergence, the
triangle inequality and $\mathrm{d}(\hat{\bm{\rho}}_{n},\hat{\bm{\rho}})\to
0$, we also get

$\displaystyle\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}})^{2}\leq\liminf_{n\to\infty}\mathrm{d}(\bm{\rho}_{n}^{+},\hat{\bm{\rho}})^{2}\leq\liminf_{n\to\infty}\mathrm{d}(\bm{\rho}_{n}^{+},\hat{\bm{\rho}}_{n})^{2}.$

Together with the minimality of $\bm{\rho}_{n}^{+}$ for
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}}_{n})$, this implies for
every $\bm{\rho}\in{\mathcal{P}_{2}}({\mathbb{R}^{d}})$:

$\displaystyle\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}^{+}|\hat{\bm{\rho}})=\frac{1}{2\tau}\mathrm{d}(\bm{\rho}^{+},\hat{\bm{\rho}})^{2}+\mathcal{E}_{\varepsilon}(\bm{\rho}^{+})\leq\liminf_{n\to\infty}\mathopen{}\mathclose{{}\left(\frac{1}{2\tau}\mathrm{d}(\bm{\rho}_{n}^{+},\hat{\bm{\rho}}_{n})^{2}+\mathcal{E}_{\varepsilon}(\bm{\rho}_{n}^{+})}\right)$
$\displaystyle\leq\liminf_{n\to\infty}\mathopen{}\mathclose{{}\left(\frac{1}{2\tau}\mathrm{d}(\bm{\rho},\hat{\bm{\rho}}_{n})^{2}+\mathcal{E}_{\varepsilon}(\bm{\rho})}\right)=\frac{1}{2\tau}\mathrm{d}(\bm{\rho},\hat{\bm{\rho}})^{2}+\mathcal{E}_{\varepsilon}(\bm{\rho})=\mathcal{E}_{\varepsilon,\tau}(\bm{\rho}|\hat{\bm{\rho}})\,,$

proving that $\bm{\rho}^{+}$ is a minimizer of
$\mathcal{E}_{\varepsilon,\tau}(\cdot|\hat{\bm{\rho}})$.

The convergence (56) is clear by narrow convergence of
$\hat{\rho}_{1,n}\to\hat{\rho}_{1}$ and $\rho_{1,n}^{+}\to\rho_{1}^{+}$. Since
(59) and
$\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}_{n})\to\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}})$
imply the inequality (58), it remains to prove (57). We show this convergence
by proving that in $L^{2}(S)$ with
$S:=\mathrm{supp}\,\zeta\subset\subset{\mathbb{R}^{d}}$, it holds for $j=1,2$:

$\displaystyle\begin{split}\nabla
F_{j}^{\prime}(\rho_{j,n}^{+})\rightharpoonup\nabla
F_{j}^{\prime}(\rho_{j}^{+}),\quad\nabla\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})&\rightharpoonup\nabla\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+}),\quad\nabla
K\ast\rho_{j,n}^{+}\rightharpoonup\nabla K\ast\rho_{j}^{+}\text{ weakly,}\\\
\rho_{j,n}^{+}&\to\rho_{j}^{+}\text{ strongly.}\end{split}$

We start by showing that $(F_{j}^{\prime}(\rho_{j,n}^{+}))_{n}$ is a bounded
sequence in $H^{1}(S)$. In order to apply Corollary 4.8, note that the
sequences $(\mathcal{E}_{\varepsilon}(\hat{\bm{\rho}}_{n}))_{n}$,
$(\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}_{n}))_{n}$ and
$(\mathcal{H}_{c}(\hat{\bm{\rho}}_{n}))_{n}$ are bounded from above, and that
$(\mathcal{E}_{\varepsilon}(\bm{\rho}^{+}_{n}))_{n}$ and
$(\mathcal{H}_{c}(\bm{\rho}^{+}_{n}))_{n}$ are bounded from below, which
follows from non-negativity of $\mathcal{E}_{\varepsilon}$ and the general
inequality

$\displaystyle\int_{\mathbb{R}^{d}}\rho\log\rho\,\mathrm{d}x\geq
1-{\mathfrak{m}}_{2}[\rho]-\pi^{d/2}\qquad\text{ for all
}\rho\in{\mathcal{P}_{2}}({\mathbb{R}^{d}}),$ (60)

and from the boundedness of the second moments of $\rho_{j,n}^{+}$. Hence, the
right hand sides of (55) are bounded from above by a constant for all $n$,
which implies uniform boundedness of $\nabla F_{1}^{\prime}(\rho^{+}_{1,n})$
in $L^{2}({\mathbb{R}^{d}})$. In order to also get uniform boundedness of
$F_{j}^{\prime}(\rho_{j,n}^{+})$ in $L^{2}(S)$, note that
$F_{j}^{\prime}(\rho_{j,n}^{+})$ is uniformly bounded in
$L^{1}({\mathbb{R}^{d}})$ by inequality (24), thus uniform boundedness in
$L^{2}(S)$ follows from the Poincaré inequality in the smooth and bounded set
$S$. Hence, the sequence $(F_{j}^{\prime}(\rho_{j,n}^{+}))_{n}$ is bounded in
$H^{1}(S)$. Alaoglu’s and Rellich’s theorems thus imply
$F_{j}^{\prime}(\rho_{j,n}^{+})\to u_{j}^{+}$ strongly in $L^{2}(S)$, weakly
in $H^{1}(S)$, for some $u_{j}^{+}\in H^{1}(S)$. The identification
$u_{j}^{+}=F_{j}^{\prime}(\rho_{j}^{+})$ now follows from the narrow
convergence $\rho_{j,n}^{+}\to\rho_{j}^{+}$, since
$\rho_{j,n}^{+}\to(F_{j}^{\prime})^{-1}(u_{j}^{+})$ pointwise almost
everywhere along a subsequence. This finishes the proof of the weak
$L^{2}(S)$-convergence $\nabla
F_{j}^{\prime}(\rho_{j,n}^{+})\rightharpoonup\nabla
F_{j}^{\prime}(\rho_{j}^{+})$.

In order to obtain weak $L^{2}(S)$-convergence
$\nabla\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})\rightharpoonup\nabla\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})$,
we write

$\displaystyle\nabla\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})=\nabla\theta_{j}\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\rho_{1,n}^{+}),F_{2}^{\prime}(\rho_{2,n}^{+})}\right)=\sum_{i=1,2}\theta_{j,i}\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\rho_{1,n}^{+}),F_{2}^{\prime}(\rho_{2,n}^{+})}\right)\nabla
F_{i}^{\prime}(\rho_{i,n}^{+}).$

Denoting
$\bm{u}_{n}^{+}:=(u_{1,n}^{+},u_{2,n}^{+}):=\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\rho_{1,n}^{+}),F_{2}^{\prime}(\rho_{2,n}^{+})}\right)$,
it thus follows from assumption (15):

$\displaystyle\|\nabla\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})\|_{L^{2}(S)}\leq\sum_{i=1,2}\mathopen{}\mathclose{{}\left(\int_{S}\theta^{2}_{j,i}\mathopen{}\mathclose{{}\left(\bm{u}_{n}^{+}}\right)|\nabla
u_{i,n}^{+}|^{2}\,\mathrm{d}x}\right)^{\frac{1}{2}}\leq\sum_{i=1,2}\kappa_{j,i}\,\|u_{i,n}^{+}\|_{L^{2}(S)},$

which is uniformly bounded for all $n$, since
$u_{i,n}^{+}=F_{i}^{\prime}(\rho_{i,n}^{+})$ is uniformly bounded in
$L^{2}(S)$. Hence up to subsequence,
$\nabla\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})\rightharpoonup v$ for
some $v\in L^{2}(S;{\mathbb{R}^{d}})$. In order to identify the limit
$v=\nabla\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})$, notice that for every
$\psi\in C^{\infty}_{c}(S;{\mathbb{R}^{d}})$, it holds

$\displaystyle\int_{S}\nabla\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})\cdot\psi\,\mathrm{d}x$
$\displaystyle=-\int_{S}\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})\
\mathrm{div}\,\psi\,\mathrm{d}x$
$\displaystyle\to-\int_{S}\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})\
\mathrm{div}\,\psi\,\mathrm{d}x=\int_{S}\nabla\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})\cdot\psi\,\mathrm{d}x,$

where we have used that
$\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})\to\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})$
strongly in $L^{2}(S)$, which follows from
$\partial_{r_{j}}h(\rho_{1,n}^{+},\rho_{2,n}^{+})=\theta_{j}\mathopen{}\mathclose{{}\left(F_{1}^{\prime}(\rho_{1,n}^{+}),F_{2}^{\prime}(\rho_{2,n}^{+})}\right)$
and the $L^{2}(S)$-convergence $F_{i}^{\prime}(\rho_{i,n}^{+})\to
F_{i}^{\prime}(\rho_{i}^{+})$, since $\theta_{j}$ is globally
$(\kappa_{j,1}+\kappa_{j,2})$-Lipschitz by (15). This implies
$v=\nabla\partial_{r_{j}}h(\rho_{1}^{+},\rho_{2}^{+})$.

To prove $\nabla K\ast\rho_{j,n}^{+}\rightharpoonup\nabla K\ast\rho_{j}^{+}$
weakly in $L^{2}(S)$, fix any $\varphi\in
L^{2}({\mathbb{R}^{d}};{\mathbb{R}^{d}})$ with $\varphi\equiv 0$ in
${\mathbb{R}^{d}}\setminus S$. It now holds

$\displaystyle\int_{S}\varphi\cdot\nabla
K\ast\rho_{j,n}^{+}\,\mathrm{d}x=-\int_{\mathbb{R}^{d}}(\nabla
K\ast\varphi)\,\rho_{j,n}^{+}\,\mathrm{d}x\to-\int_{\mathbb{R}^{d}}(\nabla