\documentclass[JEP,XML,SOM,Unicode, NoFloatCountersInSection,published]{cedram}
\datereceived{2023-03-16}
\dateaccepted{2024-09-25}
\dateepreuves{2024-10-07}

\usepackage[scr=rsfs,cal=euler]{mathalfa}
\let\line\xspace
\hyphenation{hypo-thesis hypo-theses repla-ced datum data}
\multlinegap0pt
\newenvironment{enumeratei}
{\bgroup\def\theenumi{\roman{enumi}}\def\theenumii{\arabic{enumii}}\begin{enumerate}}
{\end{enumerate}\egroup}
\newcommand\mto{\mathchoice{\longmapsto}{\mapsto}{\mapsto}{\mapsto}}

\makeatletter
\def \rightharpoonupfill{$\m@th \mathord-\mkern-6mu \cleaders\hbox{$\mkern-2mu \mathord- \mkern-2mu$}\hfill\mkern-6mu \mathord \rightharpoonup$}
\makeatother
\newdimen\lengtharrow
\newbox\exponentbox

\def\Rightharpoonupts#1%
{\setbox\exponentbox=\hbox{$#1$}
\lengtharrow=\wd\exponentbox
\ifdim \lengtharrow=0pt
		\mathrel{\smash{\mathop{\hbox to 6mm{\rightharpoonupfill}}}}
\else \ifdim\lengtharrow<6mm
\lengtharrow=6mm
\else \advance\lengtharrow by 1mm
\fi
		\overset{\raisebox{3pt}{$\scriptstyle#1$}}{\mathrel{\smash{\hbox to\lengtharrow{\rightharpoonupfill}}}}
	\fi}
	
\def\Rightharpoonupds#1%
{\setbox\exponentbox=\hbox{$#1$}
\lengtharrow=\wd\exponentbox
\ifdim \lengtharrow=0pt
		\mathrel{\smash{\mathop{\hbox to 6mm{\rightharpoonupfill}}}}
\else \ifdim\lengtharrow<6mm
\lengtharrow=6mm
\else \advance\lengtharrow by 1mm
\fi
		\overset{\raisebox{5pt}{$\scriptstyle#1$}}{\mathrel{\smash{\hbox to\lengtharrow{\rightharpoonupfill}}}}
	\fi}
	
\def\Rightharpoonup#1{\mathchoice{\Rightharpoonupds{#1}}{\Rightharpoonupts{#1}}{\Rightharpoonupts{#1}}{\Rightharpoonupts{#1}}}

\def\longrightharpoonup{\Rightharpoonup{}}

\newcommand{\Psfrac}[2]{(\sfrac{#1}{#2})}
\newcommand{\psfrac}[2]{\sfrac{(#1)}{#2}}
\newcommand{\spfrac}[2]{\sfrac{#1}{(#2)}}
\newcommand{\pspfrac}[2]{\sfrac{(#1)}{(#2)}}

\newcommand{\ini}{\mathrm{ini}}

\DeclareMathOperator{\Leb}{Leb}
\DeclareMathOperator{\Lip}{Lip}

\graphicspath{{lagoutiere-et-al_figures}}
\usepackage{pgfplots}
\pgfplotsset{compat=1.9}
\usepgfplotslibrary{external}
\tikzexternalize[prefix=tikz/]

\DeclareMathOperator{\dv}{div}

\def\NN{\mathbb{N}}
\def\OO{\mathbb{O}}
\def\RR{\mathbb{R}}
\def\SS{\mathbb{S}}
\def\ZZ{\mathbb{Z}}
\def\E{\mathbb{E}}
\def\PP{\mathbb{P}}
\def\WW{\mathbb{W}}
\def\intquad {\int\mbox{\hspace{-2mm}}\int\mbox{\hspace{-2mm}}\int\mbox{\hspace{-2mm}}\int}
\newcommand{\limi}[2]{\underset{#1 \to #2}{\longrightarrow}}
\newcommand{\wslimi}[2]{{\underset{#1 \to #2}{\Rightharpoonup{\ast}}}}

\let\ds\displaystyle

\def\g#1{{\bf #1}}
\def\eps{\varepsilon}
\def\bar#1{{\overline #1}}
\def\pa{\partial}
\def\lip{\mathrm{Lip}}
\def\dt{\Delta t}
\def\dx{\Delta x}
\def\dS{\pa_t S + v\pa_x S}
\def\ddS{\pa_t S + v\cdot\nabla_x S}
\def\calA{{\mathcal A}}
\def\calB{{\mathcal B}}
\def\calC{{\mathcal C}}
\def\calD{{\mathcal D}}
\def\calF{{\mathcal F}}
\def\calI{{\mathcal I}}
\def\calM{{\mathcal M}}
\def\calP{{\mathcal P}}
\def\calT{{\mathcal T}}
\def\calU{{\mathcal U}}
\def\calV{{\mathcal V}}
\def\calW{{\mathcal W}}

\def\smes{\mathcal{S}_{\mathcal{M}}}
\def\Dpetit{{\mbox{\tiny $\Delta$}}}
\def\achapo{\widehat{a}}
\def\bchapo{\widehat{b}}
\def\Xchapo{\widehat{X}}
\def\nabWchapo{\widehat{\nabla W}}
\def\widetilderho{\widetilde{\rho}}
\def\widetildeX{\widetilde{X}}
\def\widetildea{\widetilde{a}}
\def\Ytilde{\widetilde{Y}}

\theoremstyle{plain}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{prop}[theorem]{Proposition}
\theoremstyle{definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{defi}[theorem]{Definition}

\datepublished{2024-10-08}
\begin{document}
\frontmatter

\title{Vanishing viscosity limit for aggregation-diffusion equations}

\author[\initial{F.} \lastname{Lagoutière}]{\firstname{Frédéric} \lastname{Lagoutière}}
\address{Universite Claude Bernard Lyon 1, CNRS, École Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208\\
69622 Villeurbanne, France}
\email{lagoutiere@math.univ-lyon1.fr}
\urladdr{http://math.univ-lyon1.fr/homes-www/lagoutiere/}

\author[\initial{F.} \lastname{Santambrogio}]{\firstname{Filippo} \lastname{Santambrogio}}
\address{Universite Claude Bernard Lyon 1, CNRS, École Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208\\
69622 Villeurbanne, France}
\email{santambrogio@math.univ-lyon1.fr}
\urladdr{https://math.univ-lyon1.fr/~santambrogio/}

\author[\initial{S.} \lastname{Tran Tien}]{\firstname{Sébastien}\nobreakauthor \lastname{Tran Tien}}
\address{Universite Claude Bernard Lyon 1, CNRS, École Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208\\
69622 Villeurbanne, France}
\email{trantien@math.univ-lyon1.fr}

\thanks{FS acknowledges the support of the Lagrange Mathematics and Computation Research Center project on Optimal Transportation and of the European Union via the ERC AdG 101054420 EYAWKAJKOS project.}

\begin{abstract}
This article is devoted to the convergence analysis of the diffusive approximation of the measure-valued solutions to the so-called aggregation equation, which is now widely used to model collective motion of a population directed by an interaction potential. We~prove, over the whole space in any dimension, a uniform-in-time convergence in Wasserstein distance in all finite-time intervals, in the general framework of Lipschitz continuous potentials, and provide an $O(\sqrt{\varepsilon})$ rate, where $\varepsilon$ is the diffusion parameter, when the potential is $\lambda$-convex. We~give an extension to some repulsive potentials and prove sharp convergence rates of the steady states towards the Dirac mass, under some uniform attractiveness assumptions.
\end{abstract}

\subjclass{35K20, 35B40, 35A35, 35Q92, 49Q22}

\keywords{Aggregation-diffusion equations, asymptotic analysis, numerical analysis, Wasserstein distance}

\altkeywords{Équations d'agrégation-diffusion, analyse asymptotique, analyse numérique, distance de Wasserstein}

\alttitle{Limite de viscosité évanescente pour des équations d'agrégation-diffusion}

\begin{altabstract}
Cet article est consacré à l'analyse de convergence de l'approximation diffusive des solutions à valeur mesure de l'équation d'agrégation, largement utilisée pour modéliser le mouvement collectif d'une population dirigée par un potentiel d'interaction. Nous prouvons, dans l'espace entier en n'importe quelle dimension d'espace, dans le cadre général des potentiels lipschitziens, la convergence uniforme en temps sur tout intervalle borné, en distance de Wasserstein, et avec un ordre de convergence $1/2$ lorsque le potentiel est $\lambda$-convexe. Nous étendons ces résultats à certains potentiels répulsifs et prouvons des taux de convergence optimaux des états stationnaires vers la masse de Dirac, sous certaines hypothèses d'attractivité uniforme.
\end{altabstract}

\maketitle
\vspace*{-\baselineskip}
{\let\line\\ \tableofcontents}
\mainmatter

\section{Introduction}

This paper addresses the vanishing viscosity limit $\varepsilon \to 0$ for the following aggregation-diffusion problem on the whole space $\RR^d$, in any dimension $d$:
\begin{subequations}\label{eq:agg_diff}
\begin{align}
& \pa_t \rho^\eps + \nabla \cdot (a[\rho^\eps] \rho^\eps) = \eps \Delta \rho^\eps, \label{eq:rho_diff} \\
& a[\rho^\eps] = - \nabla W * \rho^\eps, \label{eq:arho_diff} \\
& \rho^\eps(0, \cdot) = \rho_0^\eps, \label{eq:rho0_diff}
\end{align}
\end{subequations}
where $\eps > 0$, $W : \RR^d \to \RR$ is a given interaction potential and the sequence of initial data $(\rho_0^\eps)_{\eps>0}$ belongs to $\calP_2(\RR^d)$ the set of probability measures with finite second order moment, and converges as $\eps$ goes to 0 towards a given $\rho^{\ini} \in \calP_2(\RR^d)$.

Equation \eqref{eq:rho_diff}--\eqref{eq:arho_diff} is often used in population dynamics to describe the collective motion of a population subject to Brownian diffusion and interacting through the interaction potential $W$. The term $\nabla W * \rho^\eps(x)$ models the combined contribution of the interaction of a particle located at point $x$ with particles at all other points. These equations appear in several applications arising from physics and biology to model, for instance, swarming, chemotaxis, crowd motion, bird flocks, or fish schools, see, e.g., \cite{morale05, burger07,topaz06,topaz04,dolak05,JVChemotaxis}. The potential $W$ depends on the model we consider. For example, the celebrated parabolic-elliptic Patlak-Keller-Segel model \cite{keller71, keller71traveling} for chemotaxis with an adequate set of parameters corresponds to the aggregation-diffusion equation in dimension $d=2$ for the logarithmic potential $W(x) = \frac{1}{2\pi} \ln(|x|)$.

In this work, we~assume that the interaction potential $W$ satisfies the following properties:
\begin{enumerate}\renewcommand{\theenumi}{\textup{A}\arabic{enumi}}\setcounter{enumi}{-1}
\item\label{A0}
For all $x \in \RR^d, \, W(x)=W(-x)$ and $W(0)=0$,
\item\label{A1}
$W\in \calC^1(\RR^d\setminus\{0\})$,
\item\label{A2}
$W$ is $a_\infty$-Lipschitz continuous, for some constant $a_\infty \geq 0$ (nevertheless this assumption is not done in Section \ref{section:stationary}).
\end{enumerate}
In addition, some of our results only hold under one of the following supplementary assumptions:
\begin{enumerate}\renewcommand{\theenumi}{\textup{A}\arabic{enumi}}\setcounter{enumi}{2}
\item\label{A3}
$W$ is $\lambda$-convex for some $\lambda \leq 0$, that is, $x \mto W(x) - \Psfrac{\lambda}{2} |x|^2$ is convex,
\end{enumerate}
\begin{enumerate}\renewcommand{\theenumi}{\textup{A}\arabic{enumi}\textup{-}$p$}\setcounter{enumi}{3}
\item\label{A4p}
There exists a constant $C>0$ such that, for all $x \in \RR^d, \,\nabla W(x) \cdot x \geq C |x|^p$, where $p \geq 1$.
\end{enumerate}
Potentials satisfying assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}
but not differentiable at the origin are often referred to as pointy \cite{cjlv15,dlv20,lv17}.
These hypothesis exclude the Patlak-Keller-Segel system from the analysis above, system in which the singularities are much more complicate to understand (and can appear also for
$\varepsilon > 0$).
\begin{remark}
Note that assumption \eqref{A2} is incompatible with assumption \eqref{A4p} whenever $p>1$, thus when the latter is done, it is done {\em instead of} \eqref{A2}. This is the reason why we only consider $\lambda \leq 0$ in \eqref{A3}, since \eqref{A3} with $\lambda>0$ implies (A4-2) (incompatible with \eqref{A2}). Still, when studying well-posedness of inviscid aggregation equations, the case $\lambda >0$ can be tackled considering compactly supported data for, in that case, the support decreases in time (see \cite[Th.\,2.1]{dlv20} and \cite[Rem.\,2.14]{figalli11}): as a consequence only the local behavior of $W$ matters. When $\eps >0$, it is not clear however that we can reproduce this argument.
\end{remark}
When the potential is pointy, finite time blowup of weak solutions occurs \cite{BertozziBrandman, bertozzi09} for the inviscid problem:
\begin{subequations}\label{eq:agg}
\begin{align}
& \pa_t \rho + \nabla \cdot (a[\rho] \rho) = 0, \label{eq:rho} \\
& a[\rho] = - \nabla W * \rho, \label{eq:arho} \\
& \rho(0, \cdot) = \rho^{\ini}, \label{eq:rho0}
\end{align}
\end{subequations}
After blowup time, the solutions being possibly singular measures, the product $a[\rho] \rho$ is no longer well-defined. For $\lambda$-convex potentials, the continuation of weak solutions valued in $\calP_2(\RR^d)$ has therefore been studied through three different approaches: gradient flow solutions in the Wasserstein space \cite{figalli11}, duality solutions \textit{à la} Bouchut-James \cite{GF_dual_jv, JVChemotaxis} in one dimension of space and Filippov solutions \cite{cjlv15, lv17}. These notions of solutions turn out to be equivalent to that of solutions in the sense of distributions provided the velocity field $a[\rho]$ is replaced by:
\begin{equation}\label{eq:achapo}
\achapo[\rho](x) = -\int_{y \neq x} \nabla W(x-y) \rho(dy) = -\nabWchapo * \rho (x),
\end{equation}
where $\nabWchapo$ is defined as
\begin{equation*}%\label{eq:nabWchapo}
\nabWchapo(x) = 
\begin{cases}
\nabla W(x) &\text{if } x \in \RR^d \setminus \{0\}, \\
0 &\text{if } x = 0.
\end{cases}
\end{equation*}

Our objective in this paper is to study the convergence of the viscous solutions $(\rho^\eps)_{\eps>0}$ towards such a weak measure solution to \eqref{eq:agg}. When $W$ is $\lambda$-convex, these asymptotics had previously been mentioned in \cite{carrilloReview18}, where the authors explain how to use the techniques for the $\Gamma$-convergence of gradient flows developed by Serfaty in \cite{serfaty10}. Our method basically relies on the same arguments which actually do not require the $\lambda$-convexity of the potential but only its Lipschitz continuity -- along with the standard assumptions \eqref{A0}--\eqref{A1} (for other works with potentials that are Lipschitz continuous but not $\lambda$-convex, called repulsive, see for example \cite{raoul} and reference therein). Starting from the so-called Energy Dissipation Equality (EDE) for the viscous problem \eqref{eq:agg_diff}, we~prove lower bounds of lower semicontinuity-type on each term of the EDE. This amounts to verifying the assumptions of \cite[Th.\,2]{serfaty10}; if, in addition, the initial data is well-prepared, then we meet all the hypotheses of this theorem. However, we~deliberately pass to the limit \textit{by hand}, so as not to invoke abstract gradient flow arguments. Therefore, our proof is self-contained for the reader with minimal background regarding optimal transport. In particular, in our Theorem~\ref{thm:cvLipschitz} we recover, at the limit $\eps \to 0$, the right definition of the velocity field for \eqref{eq:agg} as defined in \eqref{eq:achapo}.

We generalize this result in Corollary \ref{coroll:cvLambdaCvx} to arbitrary $\calP_2(\RR^d)$ initial data converging in Wasserstein distance towards the initial datum $\rho^{\ini}$ of the inviscid problem, when~$W$~is, in addition, $\lambda$-convex. This is done by smoothing out the initial data and estimating the distance to the modified solutions at time $t$, which is possible since the interaction energy is $\lambda$-geodesically convex. We~then provide a convergence rate based on the differentiation formula of the Wasserstein distance between two absolutely continuous curves on the Wasserstein space. Note that, for the Newtonian potential, the vanishing viscosity limit had been established in \cite{cozzi16} in dimension $d \geq 2$ and up to the time of existence of weak solutions in $L^1 \cap L^\infty$ but, to the best of our knowledge, without convergence rates.

This article is structured as follows.\enlargethispage{.5\baselineskip}%

\begin{itemize}
\item
In Section \ref{section:preliminaries} we recall some useful results and definitions regarding optimal transport and functionals defined over the Wasserstein spaces.

\item
The main results concerning the convergence as $\eps$ tends to 0 for the evolutive equation are contained in Sections \ref{section:lambdaConvex} and \ref{section:repulsive}.

\begin{itemize}
\item In Section \ref{section:lambdaConvex}, in the framework of Lipschitz potentials, we~begin with the general convergence result of the diffusive solutions $(\rho^\eps)_{\eps>0}$ towards a solution $\rho$ to the inviscid problem \eqref{eq:agg} for well-prepared initial data. We~then relax some of our assumptions on the initial data and focus on $\lambda$-convex potentials, for which we prove that convergence still holds for arbitrary initial data $(\rho_0^\eps)_{\eps >0}$ converging towards $\rho^{\ini}$. \item We then prove that convergence occurs at rate $O(\sqrt{\eps})$ in Wasserstein distance. This is done with two different methods. The first relies on differentiating in time the distance $\WW_2$ between two smooth solutions and exploits the previously proved unquantified convergence. The second is based on the convergence estimates of an upwind-type scheme for the inviscid problem due to the first author with Delarue and Vauchelet \cite{dlv17, dlv20}.

\item Later, in Section \ref{section:repulsive}, we~show that convergence (without convergence rate) still holds, up to an extraction, for repulsive potentials that behave like $W(x) = -|x|$. The idea is to estimate, as in the $\lambda$-convex case, the distance between solutions associated with smoothed out initial data and solutions associated with a fixed initial datum~$\rho^{\ini}$. This is done by differentiating the Wasserstein distance between solutions and proving appropriate estimates on the aggregation velocity field using an additional integrability assumption on $\nabla^2 W$.
\end{itemize}
\item
Section \ref{section:stationary} is devoted to the study of the stationary problem and, in particular, we~provide higher convergence rates for the viscous steady states towards the unique steady state of the aggregation equation, that is, up to translations, the Dirac mass, when the interaction potential satisfies the key assumption (A4-1) but is not necessarily Lipschitz continuous. Under assumption \eqref{A4p} for an arbitrary $p\geq 1$, estimates are also obtained and proved to be sharp for $p=2$.

\item
We eventually illustrate our convergence results in Section \ref{section:numerics} and observe all the proved convergence rates.
\end{itemize}

\subsubsection*{Acknowledgements} The authors are indebted to Benoît Fabrèges for crucial help with the numerical code.

\section{Preliminaries}\label{section:preliminaries}

\subsection{Notations}

We denote by $\calC(\RR^d)$ the space of continuous functions from $\RR^d$ to $\RR$, and by $\calC_0(\RR^d)$ (\resp $\calC_b(\RR^d)$, $\calC_c(\RR^d)$) the subspace of continuous functions vanishing at $\infty$ (\resp of bounded continuous functions, of continuous and compactly supported functions). We~also denote by $\calM_b(\RR^d)$ the space of Borel signed measures with finite total variation, equipped with the weak topology $\sigma(\mathcal{M}_b(\RR^d),\calC_0(\RR^d))$. For a sequence $(\rho_n)_{n \in \NN} \in \calM_b(\RR^d)^\NN$ and $\rho \in \calM_b(\RR^d)$, we~denote the weak convergence of $(\rho_n)_{n \in \NN}$ towards $\rho$ by $\rho_n \wslimi{n}{\infty} \rho$.

For $\rho\in \mathcal{M}_{b}(\RR^d)$ and $r \in [0,+\infty)$, we~also denote by $M_r(\rho)$ the $r$-th moment of~$\rho$, given by $M_r(\rho) = \int_{\RR^d} |x|^r \rho(dx)$, where $|\cdot|$ is the Euclidean norm. For $\rho \in \calM_b(\RR^d)$ and~$Z$ a measurable map, we~denote by $Z_\#\rho$ the pushforward measure of $\rho$ by $Z$, which satisfies, for any $\varphi \in \calC_b(\RR^d)$,\vspace*{-3pt}
\[
\int \varphi(x)\, Z_\#\rho(dx) = \int \varphi(Z(x))\,\rho(dx).
\]
Note that, in the above equality as in the whole article, whenever the integration domain is not specified, the integrals are considered over the whole space (which is~$\RR^d$ here). If $\mu \in \calM_b(\RR^d)$ is a positive measure, we~also note $\rho \ll \mu$ whenever~$\rho$ is absolutely continuous with respect to $\mu$.

We call $\calP(\RR^d)$ the subset of $\calM_b(\RR^d)$ of probability measures and we denote, for $p \in [1, +\infty)$, $\calP_p(\RR^d) := \left\{\rho \in \calP(\RR^d),\ M_p(\rho) < +\infty\right\}$. For $\mu, \nu \in {\mathcal P}_{p}(\RR^d)$, we~define the Wasserstein distance of order $p$ between $\mu$ and $\nu$ by (see \cite{ambrosio08, otam15, villani_topics}):\vspace*{-3pt}
\begin{equation}\label{defWp}
W_p(\mu,\nu) := \inf_{\gamma\in \Gamma(\mu,\nu)} \biggl\{\iint |x-y|^p\,\gamma(dx,dy)\biggr\}^{1/p},
\end{equation}
where $\Gamma(\mu,\nu)$ is the set of measures on $\RR^d\times\RR^d$ with marginals $\mu$ and $\nu$, \ie\vspace*{-5pt}
\begin{multline*}
\Gamma(\mu,\nu) = \Bigl\{ \textstyle\gamma\in \calP_p(\RR^d\times\RR^d); \ \forall\, \xi\in \calC_0(\RR^d), \int \xi(x)\gamma(dx,dy) = \int \xi(x) \mu(dx), \\[-3pt]
\textstyle\int \xi(y)\gamma(dx,dy) = \int \xi(y) \nu(dy) \Bigr\}.
\end{multline*}
Any measure that realizes the minimum in the definition \eqref{defWp}
of $W_p$ is called an optimal plan, and the set of optimal plans is denoted by $\Gamma_0(\mu,\nu)$. The space ${\mathcal P}_{p}(\RR^d)$ equipped with the distance $W_p$ is called Wasserstein space of order $p$ and denoted $\WW_p(\RR^d)$.

We recall that the Wasserstein distance $W_p$ metrizes the weak convergence of measures in the sense that, for $(\rho_n)_{n\in\NN} \in \calP_p(\RR^d)^\NN$ and $\rho \in \calP_p(\RR^d)$, $W_p(\rho_n, \rho) \limi{n}{+\infty} 0$ if and only if $\rho_n \wslimi{n}{+\infty} \rho$ and $M_p(\rho_n) \limi{n}{+\infty} M_p(\rho)$ (see \cite[Th.\,7.12]{villani_topics}).

We shall also denote the conjugate exponent of $p$ by $p' \in [1, +\infty]$ defined by $\sfrac{1}{p} + \sfrac{1}{p'} = 1$, with the usual convention $1' = +\infty$ and $\infty' = 1$. For $\alpha \in \RR$, the positive and negative part of $\alpha$ are denoted by $\alpha^+ := \max(0,\alpha)$ and $\alpha^- := \max(0,-\alpha)$. With that convention, both $\alpha^+$ and $\alpha^-$ are always nonnegative.

Throughout this paper, we~will use the same notation $C$ to denote any positive constant.

\subsection{Curves and functionals over the Wasserstein space}

Let $p \in [1,+\infty)$ and $T > 0$. We~call curve on the metric space $\WW_p(\RR^d)$ any continuous function $\rho \in \calC([0,T], \WW_p(\RR^d))$. We~say that $\rho$ is an absolutely continuous curve if there exists $b \in L^1([0,T])$ such that $W_p(\rho_s, \rho_t) \leq \int_s^t b(\tau)d\tau$ for every $0\leq s<t\leq T$, and we denote $AC([0,T], \WW_p(\RR^d))$ the set of absolutely continuous curves on $\WW_p(\RR^d)$. We~also define for $t \in [0,T]$, the metric derivative of $\rho$ at time $t$ as:
\begin{equation*}
|\rho'_t| := \lim_{h \to 0} \frac{W_p(\rho_{t+h},\rho_t)}{h}.
\end{equation*}
If $\rho$ is a Lipschitz curve on $\WW_p(\RR^d)$, then the above limit exists for a.e. $t \in [0,T]$. Now, up to a reparametrization in time, any absolutely continuous curve can become Lipschitz continuous and therefore admits a metric derivative for almost every time.

The fundamental property of absolutely continuous curves in $\WW_p(\RR^d)$ is the link with a continuity equation:
\begin{theorem}[{\cite[Th.\,8.3.1]{ambrosio08}}]\label{thm:ACcurves}
Let $p \!\in\! (1,+\infty)$ and $T\!>\!0$.
Let $\rho \!\in\! AC([0,T], \WW_p(\RR^d))$. Then, for a.e. $t \in [0,T]$ there exists a vector field $v_t \in L^p(\rho_t, \RR^d)$ such that:
\begin{itemize}
\item the continuity equation $\pa_t \rho + \nabla \cdot (\rho v) = 0$ is satisfied in the sense of distributions
\item for a.e. $t \in [0,T]$, $\|v_t\|_{L^p(\rho_t)} \leq |\rho_t'|$.
\end{itemize}
Conversely, if we take a curve $\rho \in \calC([0,T], \WW_p(\RR^d))$ such that, for each $t \in [0,T]$, there exists a vector field $v_t \in L^p(\rho_t, \RR^d)$ with $\int_0^T \|v_t\|_{L^p(\rho_t)} dt < +\infty $ solving the continuity equation $\pa_t \rho + \nabla \cdot \rho v = 0$, then $\rho \in AC([0,T], \WW_p(\RR^d))$ and for a.e. $t \in [0,T]$, we~have $|\rho_t'| \leq \|v_t\|_{L^p(\rho_t)}$.

As a consequence, the velocity field $v$ introduced in the first part of the statement actually satisfies $\|v_t\|_{L^p(\rho_t)} = |\rho_t'|$ for a.e. $t \in [0,T]$.
\end{theorem}
We now recall the definition of the first variation of a functional defined over probability measures.
\begin{defi}
Let $F : \calP(\RR^d) \longrightarrow \RR \cup \{+\infty\}$. Assume that $\rho \in \calP(\RR^d)$ is such that:
\[
\forall \delta \in [0,1], \, \forall \mu \in \calP(\RR^d) \cap L_c^\infty(\RR^d), \quad F((1-\delta) \rho + \delta \mu) < +\infty,
\]
then we call first variation of $F$ at $\rho$, denoted $\Psfrac{\delta F}{\delta \rho}(\rho)$, any measurable function $g$ such that:
\[
\frac{d F(\rho + \delta \chi)}{d \delta}\Big|_{\delta = 0} = \int g d\chi,
\]
whenever $\chi = \mu - \rho$ for some $\mu \in \calP(\RR^d) \cap L_c^\infty(\RR^d)$, where $L_c^\infty(\RR^d)$ denotes the set of compactly supported functions in $L^\infty(\RR^d)$. If it exists, the first variation is defined up to an additive constant.
\end{defi}
We now introduce two functionals that are essential to our study, the interaction energy $\calW$ and the entropy $\calU$, defined on $\calP(\RR^d)$ by:
\begin{align*}
& \calW(\rho) = \frac 12 \iint W(x-y) \rho(dx) \rho(dy), \\
& \calU(\rho) = \begin{cases}
\int \rho \ln(\rho) & \text{if } \rho \ll \Leb, \\
+\infty & \text{ otherwise,}
\end{cases}
\end{align*}
where $\Leb$ is the Lebesgue measure on $\RR^d$. Note that, under assumption \eqref{A2}, the interaction energy $\calW(\rho)$ is finite whenever $\rho \in \calP_1(\RR^d)$. For $\eps \geq 0$, we~shall also define the energy functional as $F^\eps = \calW + \eps \calU$. One can easily show that $\Psfrac{\delta \calW}{\delta \rho}(\rho) = W * \rho$ and $\Psfrac{\delta \calU}{\delta \rho}(\rho) = \ln \rho + 1$.

A key point in our proofs will be the lower semicontinuity (l.s.c) of the above functionals so that minimization arguments apply.

\skpt
\begin{lemma}\label{lemma:W_lsc}
\begin{enumerate}
\item If $W$ is l.s.c.\ on $\RR^d$ and bounded from below, then the interaction energy $\calW$ is l.s.c.\ for the weak convergence.
\item If $W$ is Lipschitz continuous, then the interaction energy $\calW$ is Lipschitz continuous for the $W_1$ distance.
\end{enumerate}
\end{lemma}
\begin{proof}
The first claim is contained in \cite[Prop.\,7.2]{otam15}.

For the second claim, we~will prove
\[
|\calW(\rho)-\calW(\mu)|\leq \mathrm{Lip}(W)\,W_1(\rho,\mu).
\]
Indeed, we~can write $\calW(\rho)=\frac12 \int (W*\rho)d\rho$, so that we have
\[
\calW(\rho)-\calW(\mu)=\frac12 \int(W*\rho)d(\rho-\mu)+\frac12\int(W*(\mu-\rho))d\mu.
\]
We then use
\[
\biggl|\int(W*\rho)d(\rho-\mu)\biggr|\leq\mathrm{Lip}(W*\rho)\,W_1(\rho,\mu)
\]
together with $\mathrm{Lip}(W*\rho)\leq \mathrm{Lip}(W)$, and
\[
|(W*(\mu-\rho))(x)|=\biggl|\int W(x-y)(\rho-\mu)(dy)\biggr|\leq \mathrm{Lip}(W(x-\cdot))\,W_1(\rho,\mu)
\]
together with $\mathrm{Lip}(W(x-\cdot))=\mathrm{Lip}(W)$.
\end{proof}

The following lemma is proved in \cite[Prop.\,2.1]{moment}. Recall that $M_p(\rho)$ is the $p$-th moment of $\rho$.
\begin{lemma}\label{lemma:U_lsc}
There exists a constant $C$ only depending on $d$ such that the entropy functional $\calU$ satisfies $\calU(\rho)\geq -C(M_1(\rho)^{1/2}+1)$. Moreover, if
$(\rho_n)_n \in \calP(\RR^d)$ is a sequence weakly converging towards some $\rho \in \calP(\RR^d)$ such that $M_1(\rho_n)$ is bounded, then we have $\calU(\rho) \leq \liminf_{n \to +\infty} \calU(\rho_n)$.

In particular, this means that the entropy is l.s.c.\ for the $W_q$ distance for all $q \geq 1$.
\end{lemma}

In order to obtain convergence of the moments of a weakly converging sequence of probability measures, we~will often make use of the following lemma, which is a particular case of \cite[Prop.\,7.1.5]{ambrosio08}, since our assumption implies uniform integrability of the $p$-moment:
\begin{lemma}\label{lemma:boundedWpImpliesCompactWq}
Let $1 \leq p < +\infty$ and $(\rho_n)_{n \in \NN}$ be a sequence of probability measures in $\calP_p(\RR^d)$.
Assume that, for some constant $C>0$, we~have for all $n \in \NN$, $M_p(\rho_n) \leq C$.
Then, there exist a subsequence of $(\rho_n)_{n \in \NN}$ converging towards some $\rho \in \calP_p(\RR^d)$ in~$W_q$ distance for all $q \in [1,p)$.
\end{lemma}

We finally define one last functional that will be useful in our proofs. Let $p \in (1, +\infty)$. We~set $K_p = \bigl\lbrace (a,b) \in \RR \times \RR^d \mid a + \Psfrac{1}{p'} |b|^{p'} \leq 0 \bigr\rbrace$ and, for $(t,x) \in \RR_+ \times \RR^d$,
\[
f_p(t,x) = \begin{cases}
\dfrac{1}{p} \dfrac{|x|^p}{t^{p-1}} & \text{if } t>0, \\
0 & \text{if } t=0, \, x=0, \\
+\infty & \text{if } t=0, \, x \neq 0.
\end{cases}
\]
Then, for $X$ a measurable space and for $(\rho, E) \in \calM_b(X) \times \calM_b(X)^d$, we~define the $p$-Benamou-Brenier functional by:
\[
\calB_p(\rho, E) = \sup \biggl\lbrace \int a d\rho + \int b \cdot dE \,; \, (a,b) \in \calC_b(X, K_p) \biggr\rbrace.
\]
The Benamou-Brenier functional satisfies the following properties (see \cite[Prop.\,5.18]{otam15}):

\skpt
\begin{lemma}\label{lemma:BenamouBrenier}
\begin{enumeratei}
\item $\calB_p$ is convex and l.s.c.\ on $\calP(X) \times \calM(X)^d$ for the weak convergence,
\item If $\rho$ and $E$ are absolutely continuous with respect to a positive measure $\mu$, then $\calB_p(\rho, E) = \int f_p(\rho, E) d\mu$,
\item $\calB_p(\rho, E) < +\infty$ only if $E \ll \rho$,
\item In that case, if we denote by $v$ the density of $E$ with respect to $\rho$, that is $E = \rho v$, then $\calB_p(\rho, E) = \int \Psfrac{|v|^p}{p} d\rho$.
\end{enumeratei}
\end{lemma}
We also have the following symmetrization lemma, which we will repeatedly use for $V = \nabla W$:
\begin{lemma}\label{lemma:symmetrization}
Let $V$ be a bounded odd vector field on $\RR^d$, $\rho \in \calP(\RR^d)$ and $v$ a vector field on $\RR^d$ such that $v \cdot (V * \rho)$ is integrable with respect to $\rho$. Then, one has:
\[
\int v(x) \cdot (V * \rho)(x) \rho(dx) = \frac{1}{2} \iint V(x-y) \cdot (v(x)-v(y)) \rho(dx) \rho(dy).
\]
\end{lemma}
\begin{proof}
Using the fact that $V$ is odd, we~can write thanks to the change of variables $x \leftrightarrow y$:
\[
\iint V(x-y) \cdot v(x) \rho(dx) \rho(dy) = - \iint V(x-y) \cdot v(y) \rho(dx) \rho(dy).
\]
Therefore, taking the half sum of the two quantities above:
\begin{multline*}
\int v(x) \cdot (V * \rho)(x) \rho(dx) = \iint V(x-y) \cdot v(x) \rho(dx) \rho(dy) \\
= \frac 12 \left(\iint V(x-y) \cdot v(x) \rho(dx) \rho(dy)
- \iint V(x-y) \cdot v(y) \rho(dx) \rho(dy) \right) \\
= \frac{1}{2} \iint V(x-y) \cdot (v(x)-v(y)) \rho(dx) \rho(dy).\qedhere
\end{multline*}
\end{proof}
We finish with a computation of the derivative of $\calW$ along a curve satisfying a continuity equation:
\begin{lemma}\label{lemma:derEnergieInteraction}
Let $\rho$ be a curve on $\calP(\RR^d)$ that solves in the weak sense $\pa_t \rho + \nabla \cdot \rho v = 0$ with $v_t \in L^2(\rho_t)$ for a.e. $t \in [0,T]$ and $\int_0^T \|v_t\|_{L^2(\rho_t)}^2 dt < +\infty$. Then:
\begin{equation*}\label{eq:derEnergieInteraction}
\forall t \in [0,T], \quad \calW(\rho_t) - \calW(\rho_0) = \int_0^t \int (\nabWchapo * \rho_s) \cdot v_s d\rho_s.
\end{equation*}
\end{lemma}
\begin{proof}
Let $(W^\delta)_{\delta > 0}$ be an approximation of $W$ such that $W^\delta \in \calC^1(\RR^d)$, $W^\delta \limi{\delta}{0} W$ uniformly on $\RR^d$, $W^\delta$ is even, $\nabla W^\delta$ is bounded by $a_\infty$, and $\nabla W^\delta \limi{\delta}{0} \nabla W$ pointwise on $\RR^d \setminus \{0\}$.

We necessarily have $\nabla W^\delta(0) = 0$ for all $\delta>0$ and therefore $\nabla W^\delta \limi{\delta}{0} \nabWchapo$ pointwise on $\RR^d$. On the other hand, for $\delta > 0$, since $W^\delta \in \calC^1(\RR^d)$ and $W^\delta$ is even, we~have, for $t \in [0,T]$:
\begin{multline}\label{eq:Wdelta}
\frac 12 \iint W^\delta(x-y) \rho_t(dx) \rho_t(dy) - \frac 12 \iint W^\delta(x-y) \rho_0(dx) \rho_0(dy) \\
= \int_0^t \iint \nabla W^\delta(x-y) \cdot v_s(x) \rho_s(dx) \rho_s(dy) ds.
\end{multline}
Now, we~can bound the integrand on the right-hand side writing $|\nabla W^\delta(x-y) \cdot v_s(x)| \leq a_\infty |v_s|$. Noting that we have
\[
\int_0^t \iint |v_s(x)| \rho_s(dx)\rho_s(dy)ds = \int_0^t \|v_s\|_{L^1(\rho_s)} ds \leq \sqrt{T} \Big(\int_0^T \|v_s\|_{L^2(\rho_s)}^2 ds\Big)^{1/2} < +\infty,
\]
we can then use Lebesgue's dominated convergence theorem with respect to the measure $\rho_s(dx)\rho_s(dy)ds$ to get that the right-hand side in equation \eqref{eq:Wdelta} converges to
\[
\int_0^t \iint \nabWchapo(x-y) \cdot v_s(x) \rho_s(dx) \rho_s(dy) ds,
\]
which is equal to
\[
\int_0^t \int (\nabWchapo * \rho_s) \cdot v_s d\rho_s.
\]
The uniform convergence of $W^\delta$ towards $W$ ensures convergence of the left-hand side, which concludes the proof.
\end{proof}

\subsection{Preliminary results}

We recall the following result of existence of a characteristic flow and well-posedness of measure-valued solutions to \eqref{eq:agg}:
\begin{theorem}[{\cite[Th.\,2.12 \& 2.13]{figalli11}, \cite[Th.\,2.5 \& 2.9]{cjlv15}}]\label{thm:cjlv15}
Assume $W$ satisfies hypotheses \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3} and let $\rho^{\ini}$ be given in $\calP_2(\RR^d)$.
Then, there exists a unique solution $\rho \in \calC([0,+\infty), \WW _2(\RR^d))$
satisfying, in the sense of distributions, the aggregation problem \eqref{eq:agg} where $a[\rho]$ is replaced by $\achapo[\rho]$ as defined in \eqref{eq:achapo}.

This solution may be represented as the family of pushforward measures
$(\rho_t:=Z_{\rho}(t,\cdot){}_\# \rho^{\ini})_{t \geq 0}$
where $(Z_{\rho}(t,\cdot))_{t \geq 0}$ is the unique Filippov characteristic flow associated with the one-sided Lipschitz velocity
field $\achapo[\rho]$.
Besides, if $\rho$ and $\mu$
are the respective solutions to \eqref{eq:agg}
with $\rho^{\ini}$ and $\mu^{\ini}$ as initial conditions in ${\mathcal P}_{2}(\RR^d)$,
then, for all $t \geq 0$,
\begin{equation}\label{ineq:contractionLambdaCvx:cjlv5}
W_2(\rho_t,\mu_t) \leq e^{- \lambda t}
W_2(\rho^{\ini},\mu^{\ini}).
\end{equation}
\end{theorem}

In \cite{carrilloAsymptotic21}, Carrillo, Gómez-Castro, Yao and Zeng proved the following well-posedness and regularity Theorem for aggregation-diffusion equations with Lipschitz symmetric potentials. They prove existence and uniqueness through a fixed-point argument and regularity applying a bootstrap argument in adequate fractional Sobolev spaces. The solutions they define are mild solutions, which are stronger than our definition of solutions, which is in the sense of distributions. We~recall the definition of the heat kernel used in the mild formulation:
\[
G_t(x) = \frac{1}{(4\pi t)^{d/2}}\, e^{-\sfrac{|x|^2}{4t}}.
\]
\begin{theorem}[{\cite[Th.\,1.1, 2.1 \& 2.2]{carrilloAsymptotic21}}]\label{thm:agg_diff:WP}
Assume that $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}. Let $\eps > 0$ and $\rho_0^\eps \in \calP(\RR^d)$.
\begin{enumerate}
\item For all $T>0$, there exists a unique solution $\rho^\eps \in \calC([0,T], \calP(\RR^d))$ to the aggregation-diffusion problem \eqref{eq:agg_diff} in the sense that:
\[
\forall t \in [0,T], \qquad \rho_t^\eps = G_{\eps t} * \rho_0^\eps + \int_0^t \big(\nabla G_{\eps(t-s)}\big) * \big((\nabla W * \rho_s^\eps) \rho_s^\eps\big) ds.
\]
\item This solution is actually a classical solution that belongs, for all $T>0$, to $\calC((0,T], W^{k,p}(\RR^d))$ for all $k \in \NN$ and $p \in [1,+\infty]$ in the general case, and to $\calC([0,T], W^{s,p}(\RR^d))$ for all $s \geq 0$ and $p \in [1,+\infty]$ if we assume that $\rho_0^\eps \in W^{s,p}(\RR^d)$.
\end{enumerate}
\end{theorem}
\begin{remark}
In \cite{carrilloAsymptotic21}, the authors state the second item of the above Theorem under the assumption that $W \in \calW^{1,\infty}(\RR^d)$ and assuming that the initial datum belongs to $L^1_+(\RR^d)$ with total unit mass instead of $\calP(\RR^d)$. It seems to us that $W \in L^\infty$ is only required to obtain sharp decay of the energy functional and that the $L^1$ assumption on $\rho_0^\eps$ is only useful to simplify the notations.
\end{remark}
In the above theorem, we~actually have $\rho^\eps \in \calC([0,+\infty[, \WW_2(\RR^d))$. Indeed, as we will see in the proof of our Theorem \ref{thm:cvLipschitz} (see equation \eqref{eq:rhoEpsHolder}), $\frac 12$-Hölder continuity in time follows automatically from a uniform bound with respect to $t \in [0,T]$ on $M_2(\rho_t^\eps)$. This in turn comes from the following computations, where we use, first, integration by parts, and, second, the symmetrization Lemma \ref{lemma:symmetrization}:
\[
\begin{cases}
\begin{aligned}
\frac{d}{dt} M_2(\rho_t^\eps) &= \int |x|^2 \pa_t \rho_t^\eps = \int |x|^2 \nabla \cdot \big( (\nabla W * \rho_t^\eps) \rho_t^\eps \big) + \eps \int |x|^2 \Delta \rho_t^\eps \\
&= - 2 \int x \cdot (\nabla W * \rho_t^\eps) d\rho_t^\eps + 2\eps d,
\end{aligned}
\\
\begin{aligned}
- 2 \int x \cdot (\nabla W * \rho_t^\eps) d\rho_t^\eps &= - \iint \nabla W(x-y) \cdot (x-y) \rho_t^\eps(dx) \rho_t^\eps(dy) \leq 2 a_\infty M_1(\rho_t^\eps) \\
&\leq 2 a_\infty \sqrt{M_2(\rho_t^\eps)}.
\end{aligned}
\end{cases}
\]
We thus get $\frac{d}{dt} M_2(\rho_t^\eps) \leq 2 a_\infty \sqrt{M_2(\rho_t^\eps)} + 2\eps d$ which implies, using a nonlinear Grönwall lemma, that $M_2(\rho_t^\eps)$ is bounded over a finite horizon.

We finish by mentioning the special case of the dimension $d = 1$, with potentials of the form $W(x) = a|x|$ for $a\in\RR \setminus \{0\}$, for which the vanishing viscosity limit can be obtained using the correspondence with Burgers' equation. Indeed, let us set, for $\eps \geq 0$, $u^\eps(t,x) = a\big(1 - 2 f^\eps(t)\big)$, where $f^\eps(t)$ is the cumulative distribution function of $\rho^\eps_t$. One can show that $\rho^\eps$ solves \eqref{eq:rho_diff} if and only if $u^\eps$ solves the viscous Burgers equation:
\begin{equation}\label{eq:burgersVisc}
\pa_t u^\eps + \pa_x \frac{(u^\eps)^2}{2} = \eps \pa_{xx} u^\eps,
\end{equation}
and, similarly, $\rho$ solves the aggregation equation \eqref{eq:rho} with the correct velocity field $\achapo[\rho]$ if and only if $u$ solves Burgers' equation (see \cite{bonaschi13, relax21, jv15}). Using the fact that, in dimension $d=1$, we~have the representation $W_1(\rho^\eps_t, \rho_t) = \|f^\eps(t) - f(t)\|_{L^1(\RR)}$ and combining with Kuznetsov's estimate hereafter for the viscous Burgers equation (see~\cite{kuznetsov76}):
\[
\|u^\eps(t) - u(t)\|_{L^1(\RR)} \leq C \mathrm{TV}(u_0) \sqrt{\eps t},
\]
where $\mathrm{TV}$ denotes the total variation and $C$ is a positive constant, we~deduce the following proposition:
\begin{prop}\label{prop:cvViaBurgers}
Assume $d=1$ and $W(x) = a|x|$ for some constant $a\in \RR \setminus \{0\}$. Let $\rho^{\ini} \in \calP_2(\RR)$, set $\rho_0^\eps = \rho^{\ini}$ for all $\eps>0$ and let $(\rho^\eps)_{\eps > 0}$ be the sequence of weak solutions to \eqref{eq:agg_diff}.

Then, for all $T>0$, $(\rho^\eps)_{\eps > 0}$ converges in $W_1$ distance and uniformly on $[0,T]$, towards a solution $\rho \in \calC([0,T], \WW _2(\RR))$ to \eqref{eq:agg} with the velocity field $a[\rho]$ being replaced by $\achapo[\rho]$ as defined in \eqref{eq:achapo}. More precisely, we~have:
\[
\forall t \in [0,T], \quad W_1(\rho_t^\eps, \rho_t) \leq C \sqrt{\eps t},
\]
where the constant $C>0$ depends on $a_\infty$ only.
\end{prop}
In the case of one initial Dirac mass $\rho^{\ini} = \delta_0$, one can even obtain convergence of $\rho^\eps$ towards $\rho$ at order 1 with respect to $\eps$ using simple scaling arguments. The initial data to the Burgers problem is $u^{\ini} = 1 - 2H_0(x)$, and the solution to the inviscid Burgers problem is stationary, given by $u(t) = u^{\ini}$. One can also show that there exists a stationary solution to equation \eqref{eq:burgersVisc} of the form $v^\eps(t,x) = V(\sfrac{x}{\eps})$, with $V(-\infty) = 1$, $V(+\infty) = -1$ and $V'(\pm\infty = 0)$. We~then have using a contraction property of the viscous Burgers equation and the stationarity of $v^\eps$ and $u$:
\begin{align*}
\|u^\eps(t) - u(t)\|_{L^1} &\leq \underbrace{\|u^\eps(t) - v^\eps(t)\|_{L^1}}_{\leq \|u^{\ini} - v^\eps(0)\|_{L^1}} + \underbrace{\|v^\eps(t) - u(t)\|_{L^1}}_{=\|v^\eps(0) - u^{\ini}\|_{L^1}}\\
&\leq 2 \int \big| u^{\ini}(\sfrac{x}{\eps}) - V(\sfrac{x}{\eps}) \big| dx
\leq 2 \eps \int \big| u^{\ini} - V \big|,
\end{align*}
which gives $W_1(\rho^\eps_t, \rho_t) \leq C\eps$ with $C>0$ independent of time. This result can be extended to the case of a finite sum of Dirac masses as initial datum, using the arguments of Teng and Zhang \cite{tengzhang97} to compare shocks with traveling waves. We~also refer to \cite{tang2001, tangteng97, teng98} for generalizations of this result.

\section{\texorpdfstring{$O(\eps^{1/2})$ convergence rate when the $W$ is $\lambda$-convex}{titlesection}}\label{section:lambdaConvex}

In this section, we~assume that $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}.

\subsection{Method 1: computing $\frac{d}{dt} W_2^2(\rho_t^\eps, \rho_t)$}\label{section:lambdaConvex:1}

So as to make integration by parts rigorous, we~actually compute $\frac{d}{dt} W_2^2(\rho_t^\eps, \rho_t^\delta)$ for $\eps,\delta >0$ so that $\rho^\eps$ and $\rho^\delta$ are regular (see Theorem \ref{thm:agg_diff:WP}), and then we let $\delta \to 0$. We~therefore need to know that $\rho_t^\delta$ converges in the sense of measures towards $\rho_t$.

\subsubsection{Convergence in $\calC([0,T], \WW_1(\RR^d))$ without convergence rate}

Let $T>0$ and let $\rho^\eps \in \calC([0,T], \WW_2(\RR^d))$ be the solution to the aggregation-diffusion problem \eqref{eq:agg_diff} on $[0,T] \times \RR^d$, as given by Theorem \ref{thm:agg_diff:WP}. Let us denote $v^\eps = - \nabla W * \rho^\eps - \eps \sfrac{\nabla \rho^\eps}{\rho^\eps}$ so that the continuity equation $\pa_t \rho^\eps + \nabla \cdot \rho^\eps v^\eps = 0$ is satisfied in the sense of distributions. We~formally have, by definition of the first variation and then by integration by parts:
\begin{equation}\label{eq:derFeps}
\frac{d}{dt} F^\eps(\rho_t^\eps) = \int \frac{\delta F^\eps}{\delta \rho}(\rho_t^\eps) \pa_t \rho_t^\eps = \int \nabla \frac{\delta F^\eps}{\delta \rho}(\rho_t^\eps) \cdot v_t^\eps d\rho_t^\eps = - \int \Big| \nabla \frac{\delta F^\eps}{\delta \rho}(\rho_t^\eps) \Big|^2 d\rho_t^\eps,
\end{equation}
where, in the last equality, we~used the identity $\Psfrac{\delta F^\eps}{\delta \rho}(\rho) = W * \rho + \eps (\ln \rho +1)$ to deduce that $v_t^\eps$ is nothing else than $-\nabla \Psfrac{\delta F^\eps}{\delta \rho}(\rho_t^\eps)$. Proving rigorously \eqref{eq:derFeps} can be made using an easy adaptation of Lemma \ref{lemma:derEnergieInteraction}. Integrating \eqref{eq:derFeps} over time then yields:
\[
\forall t \in [0,T], \quad F^\eps(\rho^\eps_0) = F^\eps(\rho^\eps_t) + \int_0^t \int \Big| \nabla \frac{\delta F^\eps}{\delta \rho}(\rho_s^\eps) \Big|^2 d\rho_s^\eps ds.
\]
Let us only use this equality as an inequality as it will turn out sufficient for passing to the limit, and let us write $| \nabla \Psfrac{\delta F^\eps}{\delta \rho}(\rho_s^\eps) |^2$ as the half-sum
\[
\frac 12 \Bigl(|v^\eps_s|^2 + \Big| \nabla \frac{\delta F^\eps}{\delta \rho}(\rho_s^\eps) \Big|^2 \Bigr)
\]
so as to recover a link between the velocity $v$ and the functional $F$ at the limit $\eps \to 0$. This way, we~recover the so-called energy dissipation equality (EDE, that we use as an inequality in our paper):
\begin{multline}\label{eq:EDE_eps}
\forall t \in [0,T], \quad F^\eps(\rho^\eps_0) \geq F^\eps(\rho^\eps_t) + \frac 12 \int_0^t \int |v^\eps_s|^2 d\rho^\eps_s ds\\[-5pt]
+ \frac 12 \int_0^t \int \Big| \nabla \frac{\delta F^\eps}{\delta \rho}(\rho_s^\eps) \Big|^2 d\rho_s^\eps ds.
\end{multline}
Showing a sort of lower semicontinuity, when $\eps \to 0$, of each term in \eqref{eq:EDE_eps}, we~will prove that up to successive extractions, $(\rho^\eps)_{\eps>0}$ converges towards a measure $\rho$ that satisfies a continuity equation and an EDE. Combining both, we~will prove that $\rho$ solves the aggregation problem \eqref{eq:agg}. In case the solution to such a Cauchy problem is unique, the whole sequence $(\rho^\eps)_{\eps>0}$ converges towards $\rho$. This method does not require the $\lambda$-convexity but only the Lipschitz continuity of the potential $W$.
\begin{theorem}\label{thm:cvLipschitz}
Assume $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}. Let $\rho^{\ini} \in \calP_2(\RR^d)$, and let $(\rho^\eps)_{\eps > 0}$ be a sequence of weak solutions to \eqref{eq:agg_diff}.

Assume that the sequence of initial data $(\rho_0^\eps)_{\eps > 0}$ satisfies the following assumptions:%
\begin{subequations}\label{assumptionInitialData}
\begin{align}
& \limsup_{\eps \to 0} F^\eps(\rho_0^\eps) \leq F(\rho^{\ini}), \label{assumptionInitialData:1} \\
& \lim_{\eps \to 0} W_2(\rho_0^\eps,\rho^{\ini}) = 0. \label{assumptionInitialData:4}
\end{align}
\end{subequations}
Then, for all $T>0$, $(\rho^\eps)_{\eps > 0}$ converges up to a subsequence, in $W_1$ distance and uniformly on $[0,T]$, towards a solution $\rho \in \calC([0,T], \WW _2(\RR^d))$ to \eqref{eq:agg} with the velocity field $a[\rho]$ being replaced by $\achapo[\rho]$ as defined in \eqref{eq:achapo}:
\[
\sup_{t \in [0,T]} W_1(\rho_t^\eps, \rho_t) \limi{\eps}{0} 0.
\]
If the solution to \eqref{eq:agg} is unique, then the whole sequence $(\rho^\eps)_{\eps > 0}$ converges towards~$\rho$.
\end{theorem}
\begin{remark}
Note that assumptions \eqref{assumptionInitialData} are automatically satisfied if the entropy $\calU(\rho_0^\eps)$ is uniformly bounded with respect to $\eps>0$. In case we take $\rho_0^\eps = \rho^{\ini}$, this corresponds to $\rho^{\ini}$ having finite entropy.
\end{remark}
The following lemma shows that it is possible to construct such a sequence of initial data:
\begin{lemma}\label{lemma:initialData}
Recall that $\rho^{\ini}$ is given in $\calP_2(\RR^d)$. For all $p \geq 1$ such that $\rho^{\ini} \in \calP_p(\RR^d)$ and for all $\alpha \in (-1,0)$, there exists a sequence $(\mu_0^\eps)_{\eps > 0}$ in $\calP_p(\RR^d)$ satisfying:
\begin{subequations}\label{smoothedInitialData}
\begin{align}
& \liminf_{\eps \to 0} F^\eps(\mu_0^\eps) \leq F(\rho^{\ini}), \label{smoothedInitialData:1} \\
& W_p(\mu_0^\eps,\rho^{\ini}) \leq Ce^{-\eps^\alpha},\label{smoothedInitialData:4}
\end{align}
\end{subequations}
where the constant $C>0$ depends on $p$ but not on $\eps$.
\end{lemma}
\begin{proof}
Let $\alpha \in (-1,0)$ and let $p \geq 1$ such that $\rho^{\ini} \in \calP_p(\RR^d)$. Let $(r_\eps)_{\eps > 0}$ be a sequence of positive real numbers to be specified later in the proof. Let $\eta \in L^1(\RR^d)$ be a nonnegative function supported on $B(0,1)$, with unit total mass, such that $\eta \ln \eta$ and $|x|^p \eta(x)$ are integrable on $\RR^d$. We~then set $\eta^\eps(x) = r_\eps^{-d} \eta(\sfrac{x}{r_\eps})$ and $\mu_0^\eps = \eta^\eps * \rho^{\ini}$. Because of the compact support of $\eta$ we have
\[
M_p(\eta^\eps * \rho^{\ini}) \leq C(M_p(\rho^{\ini})+M_p(\eta^\eps))\leq C,
\]
so that, in particular, $\mu_0^\eps \in \calP_p(\RR^d)$ for all $\eps>0$.

Firstly, let us choose $r_\eps$ so that $\eps \calU(\eta^\eps)$ goes to 0 as $\eps \to 0$. Since $\eta^\eps \ll \Leb$, we~have $\calU(\eta^\eps) = \int \eta^\eps \ln \eta^\eps$. Therefore, using the change of variables $x = r_\eps y$, one has:
\[
\begin{aligned}
\calU(\eta^\eps) &= r_\eps^{-d} \int \eta(\sfrac{x}{r_\eps}) \ln \bigl( r_\eps^{-d} \eta(\sfrac{x}{r_\eps}) \bigr) dx\\
& = \int \eta(y) \ln \big(r_\eps^{-d} \eta(y)\big) dy 
= \int \eta(y) \ln \eta(y) dy - d \ln r_\eps.
\end{aligned}
\]
Based on the above computation, we~choose $r_\eps = e^{-h_\eps / \eps}$ for some positive sequence $(h_\eps)_{\eps>0}$ such that $\lim_{\eps \to 0} h_\eps = 0$. More precisely, we~set $h_\eps = \eps^{\alpha + 1}$, that is $r_\eps = e^{-\eps^\alpha}$.

Now, using the convexity and the invariance under translation of $\calU$, we~have \hbox{$\calU(\eta^\eps * \rho^{\ini}) \leq \calU(\eta^\eps)$}, and therefore $F^\eps(\mu_0^\eps) \leq \calW(\mu_0^\eps) + \eps \calU(\eta^\eps)$. Since $\calW$ is continuous on $\WW_1(\RR^d)$ thanks to Lemma \ref{lemma:W_lsc}, we~just need the convergence $\mu_0^\eps\to\rho^{\ini}$ in $\WW_1(\RR^d)$ in order to have $\calW(\mu_0^\eps) \to \calW(\rho^{\ini})$ and hence
\[
\lim_{\eps \to 0} \calW(\mu_0^\eps) + \eps \calU(\eta^\eps) = \calW(\rho^{\ini}) = F(\rho^{\ini}).
\]
Then, \eqref{smoothedInitialData:1} will immediately follow.

We now use
\[
W_p^p(\mu_0^\eps,\rho^{\ini}) = W_p^p(\eta^\eps * \rho^{\ini}, \delta_0 * \rho^{\ini}) \leq W_p^p(\eta^\eps, \delta_0) = M_p(\eta^\eps)\to 0,
\]
where the last limit is justified by $M_p(\eta^\eps) = r_\eps^p M_p(\eta) = C e^{-p\eps^\alpha}$. This proves
\eqref{smoothedInitialData:4} since $\alpha < 0$, and in turn \eqref{smoothedInitialData:1}.
\end{proof}
Relaxing assumption \eqref{assumptionInitialData:1} can only be done under additional assumptions on the potential. In the case $W$ satisfies assumption \eqref{A3}, replacing the original initial data~$\rho^\eps_0$ by a smoothed out initial data $\mu_0^\eps$ that verifies assumptions \eqref{assumptionInitialData} and using the $\lambda$\nobreakdash-convexity of the potential to estimate the distance between $\rho^\eps$ and the new sequence of viscous solutions $\mu^\eps$, we~obtain as a byproduct of Theorem \ref{thm:cvLipschitz} the following corollary:
\begin{corollary}\label{coroll:cvLambdaCvx}
Assume $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}. Let $\rho^{\ini} \in \calP_2(\RR^d)$, and let $(\rho^\eps)_{\eps > 0}$ be the sequence of weak solutions to \eqref{eq:agg_diff}. Assume that the sequence of initial data $(\rho_0^\eps)_{\eps>0}$ converges in $\WW_2(\RR^d)$ to $\rho^{\ini}$ as $\eps \to 0$.

Then, for all $T>0$, the whole sequence $(\rho^\eps)_{\eps > 0}$ converges in $W_1$ distance, uniformly on $[0,T]$, towards the unique solution $\rho \in \calC([0,T], \WW _2(\RR^d))$ of \eqref{eq:agg} with the velocity field $a[\rho]$ being replaced with $\achapo[\rho]$ as defined in \eqref{eq:achapo}: $\sup_{t\in[0, T]} W_1(\rho_t^\eps,\rho_t) \to 0$.
\end{corollary}
\begin{proof}[Proof of Theorem \ref{thm:cvLipschitz}]
First of all, let us extract from $(\rho_\eps)_{\eps>0}$ a converging subsequence. For $\eps > 0$, recall that the continuity equation $\pa_t \rho^\eps + \nabla \cdot \rho^\eps v^\eps = 0$ is satisfied. Moreover, let us rewrite equation \eqref{eq:EDE_eps} using the identity
\[
\nabla \dfrac{\delta F^\eps}{\delta \rho}(\rho) = \nabla W * \rho + \eps \dfrac{\nabla \rho}{\rho}
\]
and split it into three terms:
\begin{equation}\label{eq:EDE_eps_split}
\begin{aligned}
\forall t \in [0,T], \quad F^\eps(\rho_0^\eps) &\geq F^\eps(\rho^\eps_t) + \frac 12 \int_0^t \int |v^\eps_s|^2 d\rho^\eps_s ds\\
&\hspace*{3cm}+ \frac 12 \int_0^t \int \Big| \nabla W * \rho^\eps_s + \eps \frac{\nabla \rho^\eps_s}{\rho^\eps_s} \Big|^2 d\rho^\eps_s ds \\
&= D_1^\eps + D_2^\eps + D_3^\eps,
\end{aligned}
\end{equation}
where $D_1^\eps$, $D_2^\eps$, $D_3^\eps$ are the 3 terms in the above right-hand side, in the same order.
Note that, if $M_2(\rho_t^\eps)$ is uniformly bounded, then $D_1^\eps$ is uniformly bounded from below thanks to the estimate in Lemma \ref{lemma:U_lsc}. In that case, using the fact that $D_3^\eps$ is nonnegative and the fact that $F^\eps(\rho_0^\eps)$ is bounded from above thanks to assumption \eqref{assumptionInitialData:1} on the initial data, we~can deduce that $\int_0^T \int |v^\eps_s|^2 d\rho^\eps_s ds \leq C$ for some constant $C>0$ independent of~$\eps$ and $t$. In particular, for all $t \in [0,T]$, $v_t^\eps \in L^2(\rho_t^\eps)$ and $\int_0^T \int |v^\eps_s|^2 d\rho^\eps_s ds < +\infty$. Using Theorem \ref{thm:ACcurves}, we~obtain that $\rho^\eps \in AC([0,T], \WW_2(\RR^d))$ and that its metric derivative exists and is bounded by the $L^2$ norm of $v^\eps_s$: $|(\rho^\eps)'_s| \leq \|v^\eps_s\|_{L^2(\rho^\eps_s)}$ for all $s \in [0,T]$. We~deduce the following bound, that is uniform with respect to $\eps$, by integration over time:
\[
\int_0^T |(\rho^\eps)'_s|^2 ds \leq C.
\]
Then, using a Cauchy-Schwarz inequality, we~get:
\begin{equation}\label{eq:rhoEpsHolder}
\begin{aligned}
\forall 0 \leq s \leq t \leq T, \quad W_2(\rho_t^\eps, \rho_s^\eps) &\leq \int_s^t |(\rho^\eps)'_\tau| d\tau \\
&\leq \biggl(\int_s^t |(\rho^\eps)'_\tau|^2 d\tau \biggr)^{1/2} \sqrt{t-s} \leq \sqrt{C(t-s)},
\end{aligned}
\end{equation}
which gives equicontinuity of $(\rho^\eps)_{\eps > 0}$ in $W_2$ distance (and therefore in $W_1$ distance). If we still assume that $M_2(\rho_t^\eps)$ is uniformly bounded, then the set $\{\rho_t^\eps, \, \eps > 0\}$ is relatively compact in $\WW_1(\RR^d)$ in virtue of Lemma \ref{lemma:boundedWpImpliesCompactWq}.
We can therefore apply Ascoli-Arzelà theorem in the space $\calC([0,T], \WW _1(\RR^d))$ to extract from $(\rho^\eps)_{\eps > 0}$ a subsequence converging in $\WW_1(\RR^d)$, uniformly in $t \in [0,T]$, towards some $\rho \in \calC([0,T], \WW _1(\RR^d))$. We~still denote this subsequence $(\rho^\eps)_{\eps>0}$. Moreover, the l.s.c.\ of the $W_2$ distance along with the weak convergence $\rho_t^\eps \wslimi{\eps}{0} \rho_t$ for all $t \in [0,T]$ allows to pass to the $\liminf_{\eps \to 0}$ in~\eqref{eq:rhoEpsHolder} to show that $\rho \in \calC([0,T], \WW _2(\RR^d))$. The limit $\rho$ is actually $1/2$-Hölder in time and satisfies the same estimate as $\rho^\eps$:
\[
\forall 0 \leq s \leq t \leq T, \quad W_2(\rho_t, \rho_s) \leq \sqrt{C(t-s)}.
\]
Let us come back to the boundedness of $M_2(\rho_t^\eps)$. This bound can actually be obtained from inequality \eqref{eq:EDE_eps_split}. Indeed, from \eqref{eq:EDE_eps_split} and assumption \eqref{assumptionInitialData:1}, we~get, since $D_3^\eps \geq 0$:
\begin{equation}\label{eq:interm}
F^\eps(\rho_t^\eps) + \frac 12 \int_0^t \int |v^\eps_s|^2 d\rho^\eps_s ds \leq C.
\end{equation}
Let us show that the second term controls $M_2(\rho_t^\eps)$ if $t\in [0,T]$. Differentiating $M_2(\rho_t^\eps)$ in time and integrating by parts, we~have:
\begin{align*}
\frac{d}{dt} M_2(\rho_t^\eps) = 2 \int x \cdot v_t^\eps(x) \rho_t^\eps(dx) \leq 2 M_2(\rho_t^\eps)^{1/2} \biggl(\int |v^\eps_t|^2 d\rho^\eps_t \biggr)^{1/2},
\end{align*}
using Cauchy-Schwarz inequality. Applying a Grönwall lemma, this implies, for all $t\in [0,T]$,
\begin{multline*}
M_2(\rho_t^\eps)^{1/2} \leq M_2(\rho_0^\eps)^{1/2} + \int_0^t \biggl(\int |v^\eps_s|^2 d\rho^\eps_s \biggr)^{1/2} ds \\
\leq M_2(\rho_0^\eps)^{1/2} + \sqrt{T} \biggl(\int_0^t \int |v^\eps_s|^2 d\rho^\eps_s ds \biggr)^{1/2},
\end{multline*}
where we used Jensen's inequality with respect to the measure $\sfrac{ds}{t}$. Finally, we~get:
\[
\int_0^t \int |v^\eps_s|^2 d\rho^\eps_s \geq \frac{1}{T}\big(M_2(\rho_t^\eps) - M_2(\rho_0^\eps)\big).
\]
Plugging this inequality into \eqref{eq:interm} and using the estimate in Lemma \ref{lemma:U_lsc} one obtains:
\[
- a_\infty M_2(\rho_t^\eps)^{1/2} - \eps (M_2(\rho_t^\eps)^{1/4} +C) + \frac{1}{2T}\big(M_2(\rho_t^\eps) - M_2(\rho_0^\eps)\big) \leq C,
\]
which provides a uniform bound on $M_2(\rho_t^\eps)$ for $t \in [0, T]$.

The point is now, for every $t \in [0,T]$, to show l.s.c.\ of each term $D_i^\eps$, $i = 1,2,3$, with respect to the $W_1$ convergence of $(\rho_t^\eps)_{\eps>0}$ towards $\rho_t$ that we just proved.

\subsubsection*{Dealing with $D_1^\eps = F^\eps(\rho^\eps_t)$}
Using Lemma \ref{lemma:W_lsc}, we~see that the $W_1$-convergence of $(\rho_t^\eps)_{\eps>0}$ towards $\rho_t$ ensures that $\lim_{\eps \to 0} \calW(\rho^\eps_t) = \calW(\rho_t)$. Besides, thanks to Lemma~\ref{lemma:U_lsc}, we~have for the entropy $\liminf_{\eps \to 0} \calU(\rho^\eps_t) \geq \calU(\rho_t)$, and we deduce in turn $\liminf_{\eps \to 0} \eps \calU(\rho^\eps_t) \geq 0$. Therefore:
\[
\liminf_{\eps \to 0} F^\eps(\rho^\eps_t) \geq F(\rho_t).
\]

\subsubsection*{Dealing with $D_2^\eps = \frac 12 \int_0^t \int |v^\eps_s|^2 d\rho^\eps_s ds$}

For $\eps > 0$, letting $E^\eps = \rho^\eps v^\eps$, a Cauchy-Schwarz inequality shows that the total variation of $E^\eps$ is uniformly bounded with respect to $\eps>0$:
\[
|E^\eps|([0,t] \times \RR^d) = \int_0^t \int |v^\eps_s| d\rho^\eps_s ds \leq \sqrt{t} \left(\int_0^t \int |v^\eps_s|^2 d\rho^\eps_s ds\right)^{1/2} \leq \sqrt{CT},
\]
Thus, up to another extraction, we~can assume that
\[
E^\eps \wslimi{\eps}{0} E
\]
for some $E \in \calM_b([0,t] \times \RR^d)^d$. Now, since $\rho^\eps$ and $E^\eps$ are absolutely continuous with respect to the Lebesgue measure on $[0,t] \times \RR^d$ as long as $\eps>0$, Lemma \ref{lemma:BenamouBrenier} ensures that $D_2^\eps$ rewrites as follows:
\[
D_2^\eps = \int_0^t \int f_2(\rho^\eps_s, E^\eps_s) dx ds = \calB_2(\rho^\eps, E^\eps).
\]
Then the lower semicontinuity of $\calB_2$ on $\calM_b([0,t] \times \RR^d) \times \calM_b([0,t] \times \RR^d)^d$ yields:
\[
\liminf_{\eps \to 0} D_2^\eps \geq \calB_2(\rho, E),
\]
which, in turn, implies that $\calB_2(\rho, E)$ is finite and therefore gives the existence of a vector-valued density $v$ verifying $E = \rho v$. Using Lemma \ref{lemma:BenamouBrenier} (iv), the above inequality rewrites:
\[
\liminf_{\eps \to 0} D_2^\eps \geq \frac 12 \int_0^t \int |v_s|^2 d\rho_s ds.
\]
In addition, this transformation also allows to pass to the limit in the continuity equation $\pa_t \rho^\eps + \nabla \cdot E^\eps = 0$, which is now linear. Indeed, letting $\eps \to 0$ in the weak formulation, one easily gets $\pa_t \rho + \nabla \cdot (\rho v) = 0$. This shows that the limit density $\rho$ is still a solution to a continuity equation, and the link between the velocity field $v$ and the functional $F$ will be made thorough when passing to the limit $\eps \to 0$ in the EDE~\eqref{eq:EDE_eps}.

\subsubsection*{Dealing with $D_3^\eps = \frac 12 \int_0^t \int \big| \nabla W * \rho^\eps_s + \eps \sfrac{\nabla \rho^\eps_s}{\rho^\eps_s} \big|^2 d\rho^\eps_s ds$}

As it is standard when dealing with terms belonging to $L^2(\rho^\eps_s)$, we~set
\[
G^\eps = (\nabla W * \rho^\eps) \rho^\eps + \eps \frac{\nabla \rho^\eps}{\rho^\eps} \rho^\eps,
\]
so that $D_3^\eps = \calB_2(\rho^\eps, G^\eps)$.

We deduce from \eqref{eq:EDE_eps_split} that $D_3^\eps$ is uniformly bounded with respect to $\eps$, which implies that~$G^\eps$ is uniformly bounded in $\calM_b([0,t] \times \RR^d)^d$. Therefore, up to another extraction, we~can assume that $G^\eps \wslimi{\eps}{0} G$ for some $G \in \calM_b([0,t] \times \RR^d)^d$. Since $W$ is Lipschitz, we~have
\[
\int_0^t \int \big|\nabla W * \rho^\eps_s\big| d\rho^\eps_s ds \leq a_\infty t,
\]
thus $(\nabla W * \rho^\eps) \rho^\eps$ is uniformly bounded too in $\calM_b([0,t] \times \RR^d)^d$.

As a consequence, the difference $\eps \Psfrac{\nabla \rho^\eps}{\rho^\eps} \rho^\eps$ is also uniformly bounded in \hbox{$\calM_b([0,t] \times \RR^d)^d$}. Now, its limit when $\eps \to 0$ is 0 in the sense of distributions. Indeed, for $\xi \in \calC_c^\infty(\RR^d)$,
\[
\left\langle\eps \nabla \rho^\eps, \xi\right\rangle = - \eps \int_0^t \int \nabla \xi d\rho^\eps,
\]
which can be bounded, for instance, by $\eps t \|\nabla \xi\|_{L^\infty}$ and therefore goes to 0 as $\eps \to 0$. We~deduce that $\eps \Psfrac{\nabla \rho^\eps}{\rho^\eps} \rho^\eps$ actually converges in the sense of measures towards 0, hence the limit, in the sense of measures, of $G^\eps$ is that of $(\nabla W * \rho^\eps) \rho^\eps$.

\subsubsection*{Limit in the sense of measures of $(\nabla W * \rho^\eps) \rho^\eps$}

Let $\xi \in \calC_0([0,t] \times \RR^d)$ and let us consider the duality bracket $\langle(\nabla W * \rho^\eps) \rho^\eps, \xi\rangle$ as $\eps$ goes to 0. That quantity equals, using Lemma \ref{lemma:symmetrization} applied to the even vector field $\nabla W$:
\begin{multline}\label{eq:passlimNLterm:1}
\int_0^t \iint \nabla W(x - y) \cdot \xi(s, x) \rho^\eps_s(dx) \rho^\eps_s(dy) ds \\
= \frac 12 \int_0^t \iint \nabla W(x - y) \cdot \left( \xi(s, x) - \xi(s,y) \right) \rho^\eps_s(dx) \rho^\eps_s(dy) ds.
\end{multline}
Now, since $W$ is Lipschitz, $\nabla W$ is bounded, therefore the map
\[
(s,x,y) \mto \nabla W(x - y) \cdot \big( \xi(s, x) - \xi(s,y) \big)
\] is continuous and the weak convergence $\rho^\eps \otimes \rho^\eps \wslimi{\eps}{0} \rho \otimes \rho$ (which is equivalent to narrow convergence since we deal with probability measures) allows to pass to the limit $\eps \to 0$ in the above quantity to obtain:\vspace*{-5pt}
\begin{multline}\label{eq:passlimNLterm:2}
\lim_{\eps \to 0} \int_0^t \iint \nabla W(x - y) \cdot \xi(s, x) \rho^\eps_s(dx) \rho^\eps_s(dy) ds \\
= \frac 12 \int_0^t \iint \nabla W(x - y) \cdot \left( \xi(s, x) - \xi(s,y) \right) \rho_s(dx) \rho_s(dy) ds.
\end{multline}
Note that, until now, the value of $\nabla W(0)$ does not matter. Actually, all the integrals when $\eps > 0$ hold with respect to to the Lebesgue measure and therefore the diagonal $\{x = y\}$ can be avoided. We~therefore only need $\nabla W(z) = - \nabla W(-z)$ for nonzero $z$ to apply Lemma \ref{lemma:symmetrization}, and this do not impose any value to $\nabla W(0)$.

Now, to come back to some duality bracket tested against $\xi$, one needs to unsymmetrize the resulting expression by writing:
\begin{equation}\label{eq:passlimNLterm:3}
\begin{aligned}
\frac 12 \int_0^t \iint \nabla W(x - y) &{}\cdot \left( \xi(s, x) - \xi(s,y) \right) \rho_s(dx) \rho_s(dy) ds \\
& = \frac 12 \int_0^t \iint \nabWchapo(x - y) \cdot \xi(s, x) \rho_s(dx) \rho_s(dy) ds \\
& \hspace*{1cm} - \frac 12 \int_0^t \iint \nabWchapo(x - y) \cdot \xi(s, y) \rho_s(dx) \rho_s(dy) ds \\
& = \frac 12 \int_0^t \iint \nabWchapo(x - y) \cdot \xi(s, x) \rho_s(dx) \rho_s(dy) ds \\
& \hspace*{1cm}+ \frac 12 \int_0^t \iint \nabWchapo(x - y) \cdot \xi(s, x) \rho_s(dx) \rho_s(dy) ds \\
& = \int_0^t \iint \nabWchapo(x - y) \cdot \xi(s, x) \rho_s(dx) \rho_s(dy) ds,
\end{aligned}
\end{equation}
where we used the fact that $\nabWchapo(z) = - \nabWchapo(-z)$ for all $z \in \RR^d$, which now imposes $\nabWchapo(0) = 0$.

\begin{remark}
These computations could hold against a test function $\xi$ that is only Lipschitz on $[0,t] \times \RR^d$ provided $\nabla W(z) \leq C / |z|^{1-\beta}$ for some $\beta > 0$. Indeed, the map $(s,x,y) \mto \nabla W(x - y) \cdot \left( \xi(s, x) - \xi(s,y) \right)$ would be continuous on the diagonal and hence everywhere on $[0,t] \times (\RR^d)^2$. This could provide a way to deal with the non Lipschitz potentials $W(x) = |x|^\beta$ for $0<\beta<1$, that are more singular than the Lipschitz potentials but are still less singular than the logarithmic potential. However, extra difficulties arise for the limit analysis when $W$ is not Lipschitz.
\end{remark}

We finally get that $G = (\nabWchapo * \rho) \rho$ and therefore
\[
\calB_2(\rho, G)\ = \frac 12 \int_0^t \int |\nabWchapo * \rho_s|^2 d\rho_s ds.
\]
Using the l.s.c.\ of $\calB_2$ we finally get:
\[
\liminf_{\eps \to 0} D_3^\eps \geq \int_0^t \int |\nabWchapo * \rho_s|^2 d\rho_s ds.
\]
\subsubsection*{Passing to the $\liminf_{\eps \to 0}$ to recover a limit EDE}

We can now pass to the $\liminf_{\eps \to 0}$ in \eqref{eq:EDE_eps} using the assumption \eqref{assumptionInitialData:1} for the left-hand side to get the following EDE (which, once again is written as an inequality):
\begin{equation}\label{eq:EDE}
F(\rho^{\ini}) \geq F(\rho_t) + \frac 12 \int_0^t \int |v_s|^2 d\rho_s ds + \frac 12 \int_0^t \int \big| \nabWchapo * \rho_s \big|^2 d\rho_s ds.
\end{equation}
Recall that $\rho$ still solves the continuity equation $\pa_t \rho + \nabla \cdot \rho v = 0$ in the sense of distributions. Identifying the velocity $v$ is made through Lemma \ref{lemma:derEnergieInteraction} which gives:
\[
\forall t \in [0,T], \quad F(\rho_t) - F(\rho_0) = \int_0^t \int (\nabWchapo * \rho_s) \cdot v_s d\rho_s.
\]
Since $(\rho_0^\eps)_{\eps >0}$ converges to both $\rho_0$ and $\rho^{\ini}$ in $\WW_1(\RR^d)$, we~have $\rho_0 = \rho^{\ini}$. Plugging the above identity into \eqref{eq:EDE} then yields:
\[
\frac 12 \int_0^t \int \Big|v_s + \nabWchapo * \rho_s \Big|^2 d\rho_s ds \leq 0,
\]
so that $v = - \nabWchapo * \rho = \achapo[\rho]$ almost everywhere. We~deduce that $\rho$ solves the aggregation equation \eqref{eq:agg} in the sense of distributions with the correct velocity field $\achapo[\rho]$, which concludes the proof. Incidentally, the identity $v = - \nabWchapo * \rho$ confirms that the limit EDE \eqref{eq:EDE} is actually an equality.
\end{proof}

\begin{proof}[Proof of Corollary \ref{coroll:cvLambdaCvx}]
We now come back to the case of arbitrary initial data $\rho_0^\eps$ \ie we do not assume anymore that assumptions \eqref{assumptionInitialData} hold. However, we~still assume that $W_2(\rho_0^\eps, \rho^{\ini}) \limi{\eps}{0} 0$ and in addition, we~now assume $W$ to be $\lambda$-convex.

Let $(\mu^\eps_0)_{\eps > 0}$ be a sequence of smoothed out initial data for which $W_2(\mu_0^\eps,\rho^{\ini}) \limi{\eps}{0} 0$ and the assumptions \eqref{assumptionInitialData} hold on $(\mu^\eps_0)_{\eps > 0}$. We~denote by $\mu^\eps$ a solution to \eqref{eq:agg_diff} for the modified initial data $\mu^\eps_0$. Applying Theorem \ref{thm:cvLipschitz}, we~know that $\mu^\eps$ converges in $\calC([0,T], \WW_1(\RR^d))$ towards $\rho$ solution to \eqref{eq:agg} as $\eps \to 0$, up to a subsequence. But since~$W$ satisfies the assumptions of Theorem \ref{thm:cjlv15}, such a solution is unique and we deduce that the whole sequence $(\mu^\eps)_{\eps>0}$ converges towards $\rho$.

It remains to show that $W_2(\rho^\eps_t, \mu^\eps_t)$ goes to 0 as $\eps \to 0$ by estimating this quantity thanks to the $\lambda$-convexity of $W$, which is encapsulated in the following lemma.
\begin{lemma}\label{lemma:contractionLambdaCvx}
Assume $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}. Let $\rho, \mu \in \calP_2(\RR^d)$ and denote $(\varphi, \psi)$ a pair of Kantorovitch potentials from $\rho$ to $\mu$ for the quadratic cost $c(x,y) = \frac 12 |x-y|^2$. In addition, we~assume that $\rho$ or $\mu$ is an absolutely continuous measure. Then,
\begin{equation*}%\label{ineq:contractionLambdaCvx1}
\int \nabla \varphi \cdot a[\rho] d\rho + \int \nabla \psi \cdot a[\mu] d\mu \leq -\lambda W_2^2(\rho, \mu).
\end{equation*}
\end{lemma}

\skpt
\begin{remark}
\begin{enumerate}
\item In particular, we~recover the last estimate in Theorem \ref{thm:cjlv15}: if
\[
\rho, \mu \in AC_\mathrm{loc}([0,+\infty), \WW _2(\RR^d))
\]
are solution to \eqref{eq:agg} with initial data $\rho^{\ini}, \mu^{\ini} \in \calP_2(\RR^d)$ and if $\rho_t$ or $\mu_t$ is an absolutely continuous measure, the following inequality holds:
\begin{equation}\label{ineq:contractionLambdaCvx2}
\frac{d}{dt} W_2^2(\rho_t, \mu_t) \leq -2\lambda W_2^2(\rho_t, \mu_t).
\end{equation}
Indeed, this is a direct consequence of Lemma \ref{lemma:contractionLambdaCvx} and of the following computation (see \cite[Th.\,5.25]{otam15} or \cite[Th.\,8.4.7]{ambrosio08})
\begin{equation}\label{eq:derW2}
\frac{d}{dt} \frac 12 W_2^2(\rho_t, \mu_t) = \int \nabla \varphi_t \cdot v_t d\rho_t + \int \nabla \psi_t \cdot w_t d\mu_t,
\end{equation}
whenever $\rho, \mu$ satisfy the continuity equations $\pa_t \rho + \nabla \cdot \rho v = 0$, $\pa_t \mu + \nabla \cdot \mu w = 0$. Inequality \eqref{ineq:contractionLambdaCvx2} then yields the aforementioned estimate using a Grönwall lemma:
\begin{equation}
W_2(\rho_t, \mu_t) \leq e^{-\lambda t} W_2(\rho^{\ini}, \mu^{\ini}).
\end{equation}
Relaxing the assumptions that either $\rho_t$ or $\mu_t$ is an absolutely continuous measure can be done replacing $\rho_t$ by $\rho^\eps_t$ for instance, and passing to the limit $\eps \to 0$ in the resulting estimate, thanks to Corollary \ref{coroll:cvLambdaCvx}.
\item Another way of proving Lemma \ref{lemma:contractionLambdaCvx} can be found in \cite[Lem.\,4.12]{santambrogio_bourbaki}.
\end{enumerate}
\end{remark}
\begin{proof}
Assume $\rho$ is an absolutely continuous measure. Then, there exists an optimal map from $\rho$ to $\mu$ for the cost $c(x,y) = \frac 12 |x-y|^2$, which we denote $T$. Since $\nabla \psi \circ T = - \nabla \varphi$, using $\mu = T_\# \rho$ yields:
\begin{align*}
\int \nabla \varphi \cdot a[\rho] d\rho & + \int \nabla \psi \cdot a[\mu] d\mu = \int \nabla \varphi \cdot (a[\rho] - a[\mu] \circ T) d\rho \\
& = - \iint \nabla \varphi(x) \cdot \nabla W(x-y) \rho(dy) \rho(dx) \\
& \hspace*{1cm} + \iint \nabla \varphi(x) \cdot \nabla W(T(x)-y)\mu(dy)\rho(dx) \\
& = - \iint \nabla \varphi(x) \cdot \Big(\nabla W(x-y) - \nabla W\big(T(x)-T(y)\big)\Big) \rho(dy) \rho(dx),
\end{align*}
where we used once more $\mu = T_\# \rho$. Symmetrizing the above integral as in Lemma~\ref{lemma:symmetrization}, since $\nabla W$ is odd, and using $\nabla \varphi = \textup{id} - T$, we~get:
\begin{multline*}
\int \nabla \varphi \cdot a[\rho] d\rho + \int \nabla \psi \cdot a[\mu] d\mu \\
= - \frac 12 \iint \big(\nabla \varphi(x) - \nabla \varphi(y)\big) \cdot \Big(\nabla W(x-y) - \nabla W\big(T(x)-T(y)\big)\Big) \rho(dy) \rho(dx) \\
= - \frac 12 \iint \Big(x - y - \big(T(x) - T(y)\big)\Big) \cdot \Big(\nabla W(x-y) - \nabla W\big(T(x)-T(y)\big)\Big) \rho(dy) \rho(dx) \\
\leq - \frac{\lambda}{2} \iint |x - T(x) - (y - T(y))|^2 \rho(dy) \rho(dx),
\end{multline*}
where we used the $\lambda$-convexity of $W$.
We then expand the square to obtain:
\begin{multline*}
\iint |x - T(x) - (y - T(y))|^2 \rho(dy) \rho(dx) \\
= 2 \int |x - T(x)|^2 \rho(dx) - 2 \left(\iint \big(x - T(x)\big) \rho(dx)\right)^2 \leq 2 W_2^2(\rho, \mu),
\end{multline*}
which concludes the proof, as we assumed in \eqref{A3} that $\lambda \leq 0$.
\end{proof}
We now come back to the proof of Corollary \ref{coroll:cvLambdaCvx}. Denoting $(\varphi_t^\eps, \psi_t^\eps)$ a pair of Kantorovitch potentials from $\rho_t^\eps$ to $\mu_t^\eps$, and using Lemma \ref{lemma:contractionLambdaCvx} along with equation \eqref{eq:derW2}, we~get:
\[
\frac{d}{dt} \frac 12 W_2^2(\rho_t^\eps, \mu_t^\eps) \leq -\lambda W_2^2(\rho_t^\eps, \mu_t^\eps) - \eps \int (\nabla \varphi_t^\eps \cdot \nabla \rho_t^\eps + \nabla \psi_t^\eps \cdot \nabla \mu_t^\eps).
\]
The last term above being nonpositive (see \cite[Exer.\,66]{otam15} for instance), we~obtain, using a Grönwall lemma, that $W_2(\rho_t^\eps, \mu_t^\eps) \leq e^{-\lambda t} W_2(\rho_0^\eps, \mu_0^\eps)$. We~then write, for $t \in [0,T]$,
\[
W_1(\rho_t^\eps, \rho_t) \leq W_1(\rho_t^\eps, \mu_t^\eps) + W_1(\mu_t^\eps, \rho_t) \leq e^{-\lambda T} W_2(\rho_0^\eps, \mu_0^\eps) + \sup_{s \in [0,T]} W_1(\mu_s^\eps, \rho_s),
\]
where we used the fact that $W_1 \leq W_2$. Since both sequences $(\rho_0^\eps)_{\eps>0}$ and $(\mu_0^\eps)_{\eps>0}$ converge in $\WW_2(\RR^d)$ to the same limit, $W_2(\rho_0^\eps, \mu_0^\eps)$ goes to 0 as $\eps \to 0$. Moreover, $(\mu^\eps)_{\eps>0}$ converges to $\rho$ in $W_1$ distance uniformly in $[0,T]$. These two facts along with the above inequality show that $(\rho^\eps)_{\eps>0}$ also converges to $\rho$ in $\calC([0,T], \WW_1(\RR^d))$, which ends the proof of the corollary.
\end{proof}

\subsubsection{Convergence rate under the $\lambda$-convexity assumption}

We are now in position to prove the following theorem:
\begin{theorem}\label{thm:cvOrder1/2}
Assume $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}. Let $\rho^{\ini} \in \calP_2(\RR^d)$, and let $(\rho^\eps)_{\eps > 0}$ be the sequence of weak solutions to \eqref{eq:agg_diff}. Here, we~assume that $(\rho_0^\eps)_{\eps >0}$ is an arbitrary sequence in $\calP_2(\RR^d)$.

Denoting $\rho \in \calC([0,+\infty), \WW _2(\RR^d))$ the unique solution of \eqref{eq:agg} with $a[\rho]$ being replaced by $\achapo[\rho]$ as defined in \eqref{eq:achapo}, we~have the following estimate:
\begin{equation}\label{ineq:cvgRateLambdaCvx}
\forall t>0, \qquad
W_2(\rho_t^\eps, \rho_t) \leq e^{-\lambda t} W_2(\rho_0^\eps, \rho^{\ini}) + \sqrt{\frac{1 - e^{-2\lambda t}}{\lambda}} \sqrt{d\eps}.
\end{equation}
\end{theorem}
Please note that in the above estimate $\lambda \leq 0$. If $\lambda < 0$, $1 - e^{-2\lambda t}$ and $\lambda$ are negative numbers so the ratio is positive and for $\lambda=0$ the expression should be extended by continuity.
\begin{remark}\label{rmk:burgers:1}
In dimension $d=1$ with the Newtonian potential $W(x) = |x|$, the correspondence with Burgers' equation stated in Proposition \ref{prop:cvViaBurgers}, gives convergence at rate $\sqrt{\eps t}$ in $W_1$ distance. Due to $W$ being $0$-convex, our estimate leads to the same estimate but in $W_2$ distance, since taking $\lambda = 0$ in \eqref{ineq:cvgRateLambdaCvx} gives $W_2(\rho_t^\eps, \rho_t) \leq \sqrt{2d\eps t}$ for any $t>0$.

If assumption \eqref{A4p} is satisfied for some $p\geq 1$ instead of assumption \eqref{A3} and if $\rho_0^\eps = \delta_0$ for all $\eps > 0$, one can also obtain the exact same estimate using a direct computation. Indeed, in that case, $\rho_t = \delta_0$ for all $t\geq 0$ and we have, using integration by parts and Lemma \ref{lemma:symmetrization}:
\begin{align*}
\frac{d}{dt} W_2^2(\rho_t^\eps, \delta_0) & = \frac{d}{dt} \int |x|^2 \rho_t^\eps(dx) \\
& = - \iint \nabla W(x-y) \cdot (x-y) \rho_t^\eps(dx) \rho_t^\eps(dy) + 2\eps \int \rho_t^\eps(dx) \\
& \leq -C\iint |x-y|^p \rho_t^\eps(dx) \rho_t^\eps(dy) + 2\eps d, \quad \text{using assumption \eqref{A4p}}, \\
& \leq 2\eps d.
\end{align*}
Hence $W_2(\rho_t^\eps, \delta_0) \leq \sqrt{2d\eps t}$ for all $t \geq 0$.
\end{remark}
\begin{proof}
Take a sequence of initial data $(\mu_0^\eps)_{\eps>0}$ converging in $\WW_2(\RR^d)$ to $\rho^{\ini}$ as $\eps \to 0$ and denote $(\mu_\eps)_{\eps >0}$ the sequence of solutions to \eqref{eq:agg_diff} with such initial data. Let $\eps >0$. For all $\delta >0$, using Lemma \ref{lemma:contractionLambdaCvx} along with equation \eqref{eq:derW2}, we~have, denoting $(\varphi_t, \psi_t)$ a pair of Kantorovitch potentials for the quadratic cost from~$\rho^\eps_t$ to~$\rho^\delta_t$ and integrating by parts:
\begin{align*}
\frac{d}{dt} \frac 12 W_2^2(\rho_t^\eps, \mu_t^\delta) &\leq -\lambda W_2^2(\rho_t^\eps, \mu_t^\delta) - \eps \int \nabla \varphi_t \cdot \nabla \rho_t^\eps - \delta \int \nabla \psi_t \cdot \nabla \mu_t^\delta \\
&\leq -\lambda W_2^2(\rho_t^\eps, \mu_t^\delta) + \eps \int \Delta \varphi_t \, \rho_t^\eps + \delta \int \Delta \psi_t \, \mu_t^\delta.
\end{align*}
The map $x \mto \varphi_t(x) - \sfrac{|x|^2}{2}$ being concave, $\nabla^2 \varphi_t \leq I_d$, hence $\Delta \varphi_t \leq d$ and the same holds for $\psi_t$. Therefore:
\[
\frac{d}{dt} W_2^2(\rho_t^\eps, \mu_t^\delta) \leq -2\lambda W_2^2(\rho_t^\eps, \mu_t^\delta) + 2(\eps + \delta)d,
\]
which gives the result after using a Grönwall lemma and passing to the limit $\delta \to 0$ thanks to Corollary \ref{coroll:cvLambdaCvx}.
\end{proof}

\subsection{Method 2: using a numerical scheme}\label{section:lambdaConvex:2}

We now turn to a different proof of the previous result. This alternate proof will also allow to illustrate the results and the behavior of solutions with numerical results. Our main idea is to let, for a fixed $\eps>0$, $\rho_{\dx}^\eps$ be a suitable numerical approximation of the viscous solution $\rho^\eps$ to the problem~\eqref{eq:agg_diff} with fixed initial data $\rho_0^\eps = \rho^{\ini}$, and then use the formalism of \cite{dlv20} to estimate the distance from this discretized solution to the solution $\rho$ to the aggregation problem \eqref{eq:agg} in terms of $\eps$:
\[
\forall n \in \NN, \quad W_2(\rho_{\dx}^{\eps,n}, \rho_{t^n}) \leq C(t^n) \sqrt{\dx + \eps},
\]
under suitable stability conditions for the numerical scheme, and where $\dt>0$ is the time step, $t^n = n\dt$ and $\dx>0$ denotes the maximal space step. Proving the convergence of the scheme with fixed $\eps$ beforehand using compactness arguments and a Lax-Wendroff-type theorem, then letting $\dx \to 0$, we~shall deduce:
\[
\forall t>0, \quad W_2(\rho_t^\eps, \rho_t) \leq C(t) \sqrt{\eps},
\]
where we shall specify the constant $C(t)$. Note that our method also allows to deal with the case of arbitrary $\calP_2(\RR^d)$ initial data $\rho_0^\eps$ as in Theorem \ref{thm:cvOrder1/2}, but we choose to present it with initial data not depending on $\eps$ for the sake of clarity.

Let us be more specific. We~consider a Cartesian mesh of $\RR^d$ where the space step in the $i$th direction is denoted by $\dx_i>0$. The nodes of the mesh are denoted by $x_J = (J_1 \dx_1, \ldots, J_d \dx_d)$ for any $J=(J_1, \ldots, J_d)\in \ZZ^d$, and the cell centered on $x_J$ is denoted by $C_J:=[(J_1-\frac 12)\dx_1,(J_1+\frac 12)\dx_1]\times \ldots \times [(J_d-\frac 12)\dx_d,(J_d+\frac 12)\dx_d]$. We~also denote by $e_i$ the $i$th vector of the canonical basis of $\RR^d$. We~initialize our discretization with:
\begin{equation}\label{dis:rho0}
\rho_{J}^0:= \int_{C_J} \rho^{\ini}(dx)\geq 0, \qquad J \in \ZZ^d,
\end{equation}
and we consider an upwind type discretization for the aggregative part \cite{dlv17, lv17, dlv20} and an explicit discretization for the diffusive part. It writes, for $n \in \NN$,
\begin{align}
\rho_{J}^{n+1} = \rho_{J}^n &- \sum_{i=1}^{d} \frac{\dt}{\dx_i}
\Big(({a_i}^n_{J})^+ \rho_{J}^n - ({a_i}^n_{J+e_i})^- \rho_{J+e_i}^n
-({a_i}^n_{J-e_i})^+ \rho_{J-e_i}^n + ({a_i}^n_{J})^- \rho_{J}^n \Big)\notag \\[-5pt]\label{dis:explicit}
\\[-8pt]
& + \eps \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^n - 2\rho_J^n + \rho_{J-e_i}^n \right),\notag
\end{align}
where the discrete velocity is defined by:
\begin{equation}
\label{def:aij}
{a_i}^n_{J} := -\sum_{K\in \ZZ^d} \rho_{K}^n \,D_iW_J^K,
\quad \mbox{ where } \quad
D_iW_J^K := \widehat{\pa_{x_i} W}(x_J-x_K).
\end{equation}
Note that, for the sake of simplicity, we~drop, in this section, the superscripts $\eps$ when it comes to the discrete unknowns $(\rho_J^n)_{J\in \ZZ^d, n \in \mathbb N}$ but these unknowns always solve numerical schemes for the aggregation equation with viscosity $\eps>0$.

Since $W$ is even, we~also have $D_iW_{J}^{K} = -D_iW^{J}_{K}$ for all $J, K \in \ZZ^d$ and $i =1,\ldots,d$. Using a symmetrization argument as in the continuous setting, we~deduce the discrete equivalent of Lemma \ref{lemma:symmetrization}:
\begin{lemma}\label{lemma:symmetrizationDis}
Denote, for $J, K \in \ZZ^d$, $DW_J^K = (D_1W_j^K, \ldots, D_dW_J^K)$ and whenever $(v_J)_{J \in \ZZ^d}$ is a discrete vector field on the mesh, $v_J = (v_{1J}, \ldots, v_{dJ}) \in \RR^d$. For any~$(v_J)_{J \in \ZZ^d}$, we~have:
\[
\forall i=1, \ldots, d, \qquad \sum_{J\in \ZZ^d} v_{iJ} \, {a_i}^n_{J} \, \rho_{J}^n = \frac 12 \sum_{J\in \ZZ^d} \sum_{K\in \ZZ^d} D_iW_J^K \, (v_{iJ} - v_{iK}) \, \rho_{J}^n \, \rho_{K}^n,
\]
and therefore:
\[
\sum_{J\in \ZZ^d} v_{J} \cdot {a}^n_{J} \, \rho_{J}^n = \frac 12 \sum_{J\in \ZZ^d} \sum_{K\in \ZZ^d} DW_J^K \cdot (v_{J} - v_{K}) \, \rho_{J}^n \, \rho_{K}^n.
\]
\end{lemma}

\begin{proof}
Using the definition of the macroscopic velocity and the fact that $D_iW_{J}^{K} = -D_iW^{J}_{K}$, we~have:
\begin{align*}
\sum_{J\in \ZZ^d} v_{iJ} \, {a_i}^n_{J} \, \rho_{J}^n = - \sum_{J\in \ZZ^d} \sum_{K\in \ZZ^d} D_iW_J^K \, v_{iJ} \, \rho_{J}^n \, \rho_{K}^n & = \sum_{J\in \ZZ^d} \sum_{K\in \ZZ^d} D_iW_K^J \, v_{iJ} \, \rho_{J}^n \, \rho_{K}^n \\
& = \sum_{J\in \ZZ^d} \sum_{K\in \ZZ^d} D_iW_J^K \, v_{iK} \, \rho_{J}^n \, \rho_{K}^n,
\end{align*}
thanks to exchanging $K$ and $J$ in the latter sum. Taking the half sum of the first sum and the latter, we~obtain:
\[
\sum_{J\in \ZZ^d} v_{iJ} \, {a_i}^n_{J} \, \rho_{J}^n = \frac 12 \sum_{J\in \ZZ^d} \sum_{K\in \ZZ^d} D_iW_J^K \, (v_{iJ} - v_{iK}) \, \rho_{J}^n \, \rho_{K}^n.
\]
Summing over $i=1, \ldots, d$ concludes the proof.
\end{proof}
It is also natural to consider, instead of the explicit discretization of the Laplacian, an implicit discretization:
\begin{align}\label{dis:implicit}
\rho_{J}^{n+1} = \rho_{J}^n &- \sum_{i=1}^{d} \frac{\dt}{\dx_i}
\Big(({a_i}^n_{J})^+ \rho_{J}^n - ({a_i}^n_{J+e_i})^- \rho_{J+e_i}^n
-({a_i}^n_{J-e_i})^+ \rho_{J-e_i}^n + ({a_i}^n_{J})^- \rho_{J}^n \Big)\notag \\[-5pt]
\\[-8pt]
& + \eps \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^{n+1} - 2\rho_J^{n+1} + \rho_{J-e_i}^{n+1} \right),\notag
\end{align}
However, for the sake of simplicity, we~only provide the proof of our convergence estimate for the explicit scheme \eqref{dis:explicit}, although our method would also works for the implicit discretization \eqref{dis:implicit} but the computations are a bit more involved. Naturally, both schemes are asymptotic-preserving in the limit $\eps \to 0$ since they degenerate towards the upwind-type scheme of \cite{dlv20} when $\eps$ goes to 0.

One could also consider the $\theta$-scheme, for $\theta \in [0,1]$, defined by:
\begin{align}\label{dis:theta}
\rho_{J}^{n+1} = \rho_{J}^n &- \sum_{i=1}^{d} \frac{\dt}{\dx_i}
\Big(({a_i}^n_{J})^+ \rho_{J}^n - ({a_i}^n_{J+e_i})^- \rho_{J+e_i}^n
-({a_i}^n_{J-e_i})^+ \rho_{J-e_i}^n + ({a_i}^n_{J})^- \rho_{J}^n \Big)\notag \\
& + \eps (1-\theta) \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^n - 2\rho_J^n + \rho_{J-e_i}^n \right) \\
& + \eps \theta \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \bigl(\rho_{J+e_i}^{n+1} - 2\rho_J^{n+1} + \rho_{J-e_i}^{n+1} \bigr).\notag
\end{align}
The point of defining such a scheme comes from the fact that, for the heat equation $\pa_t \rho = \eps \Delta \rho$, under a parabolic CFL condition
\[
\eps \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \leq \frac{1}{2(1-2\theta)}
\]
if $\theta \in [0, 1/2)$ and unconditionally if $\theta \in [1/2, 1]$, the $\theta$-scheme is known to be convergent in $L^2$ norm at rate $O(\dt + \dx^2)$. Moreover, for $\theta = 1/2$, one obtains the so-called Crank-Nicolson scheme, which is convergent at rate $O(\dt^2 + \dx^2)$.
However, the convergence order of the $\theta$-scheme \eqref{dis:theta} for the aggregation-diffusion equation \eqref{eq:rho_diff} will anyway be limited by the order of the upwind scheme. Also, the positivity of the density can only be guaranteed if the more restrictive parabolic CFL condition
\[
a_\infty \sum_{i=1}^{d} \frac{\dt}{\dx_i} + 2\eps (1-\theta) \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \leq 1
\]
holds. Preserving a hyperbolic CFL condition thus imposes taking $\theta = 1$, which corresponds to the implicit scheme \eqref{dis:implicit}.

\begin{prop}\label{prop:cvOrder1/2:scheme}
Assume $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3} and let $\rho \in \calC([0,+\infty), \WW_2(\RR^d))$ be the unique measure solution to the aggregation equation~\eqref{eq:agg} with initial data $\rho^{\ini} \in \calP_2(\RR^d)$ as
given by Theorem \ref{thm:cjlv15}. Assume in addition that the following strict CFL condition holds:
\begin{equation}
\label{CFL}
\sum_{i=1}^d \Bigl(a_\infty \frac{\dt}{\dx_i} + 2\eps \frac{\dt}{\dx_i^2}\Bigr) < \frac 12.
\end{equation}
Denote also the reconstruction:
\[
\rho_{\dx}^{\eps,n} := \sum_{J\in \ZZ^d} \rho_J^n \delta_{x_J},
\quad n \in {\mathbb N},
\]
where $(\rho_J^n)_{J\in \ZZ^d, n \in \mathbb N}$ is defined through the explicit discretization
\eqref{dis:rho0}--\eqref{dis:explicit}--\eqref{def:aij}.
Then, there exists a constant $C>0$, depending only on $a_\infty$ and $d$, such that, for all $n\in \NN^*$,
\begin{equation}
\label{eq:cvOrder1/2:scheme:bound:1}
W_2(\rho_{t^n},\rho_{\dx}^{\eps,n} ) \leq C \sqrt{\frac{1 - e^{-4\lambda t^n}}{\lambda}}\, \sqrt{\dx + \eps} + e^{-2\lambda t^n} \dx.
\end{equation}
\end{prop}

\begin{remark}\mbox{}
In estimate \eqref{eq:cvOrder1/2:scheme:bound:1}, the $\sqrt{\dx + \eps}$ term corresponds to the error induced by the scheme \eqref{dis:explicit} and the $\dx$ term corresponds to the finite volume discretization of the initial data \eqref{dis:rho0}. As in \cite{dlv20}, one can also improve the prefactor in the exponentials to get the slightly better estimate:
\[
W_2(\rho_{t^n},\rho_{\dx}^{\eps,n} ) \leq C \sqrt{\frac{1 - e^{-2\lambda t^n}}{\lambda}}\, \sqrt{\dx + \eps} + e^{-\lambda t^n} \dx,
\]
which is similar to the estimates of the continuous setting, for instance \eqref{ineq:contractionLambdaCvx:cjlv5}, when~$\dt$ is small.
\end{remark}
In the above proposition as in the whole paper, we~do as if our discrete reconstructions $(\rho_{\dx}^\eps)_{\dx >0}$ depended only on $\dx$. Rigorously speaking, they also depend on $\dt$, but under the CFL condition \eqref{CFL} $\dt$ goes to 0 as $\dx$ goes to 0. Setting $\dt$ to be an adequate function of $\dx$, we~can therefore consider $(\rho_{\dx}^\eps)_{\dx >0}$ as sequence labeled by $\dx$ only.
\begin{theorem}\label{thm:cvOrder1/2:scheme}
Assume $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}. Let $\rho \in \calC([0,+\infty), \WW_2(\RR^d))$ be the unique measure solution to the aggregation equation \eqref{eq:agg} with initial data $\rho^{\ini} \in \calP_2(\RR^d)$ as given by Theorem \ref{thm:cjlv15} and let $(\rho^\eps)_{\eps > 0}$ be the sequence of weak solutions to \eqref{eq:agg_diff} with initial data $\rho_0^\eps = \rho^{\ini}$.

Then, there exists a constant $C>0$, depending only on $a_\infty$ and $d$, such that, for all $t>0$ the following estimate holds:
\begin{equation}\label{eq:cvOrder1/2:scheme:bound:2}
W_2(\rho_t^\eps,\rho_t) \leq C \sqrt{\frac{1 - e^{-4\lambda t}}{\lambda}}\, \sqrt{\eps}.
\end{equation}
\end{theorem}
\begin{remark}
The estimate above is slightly worse than the estimate \eqref{ineq:cvgRateLambdaCvx} that we obtain using gradient flow arguments. Although, as in the previous remark, the exponential factor can be improved to $e^{-2\lambda t}$ with a bit more technical computations, we~do not manage to obtain the same constant $C=\sqrt{d}$. Nevertheless the important fact is that the dependence with respect to $\eps$ is the same in both proofs. The advantage of the numerical proof is that it confirms the convergence of the numerical scheme and its asymptotic preserving property.
\end{remark}

\Subsubsection{Properties of the scheme}

\begin{lemma}\label{lem:CFL}
As in the continuous setting, our discretization \eqref{dis:explicit} preserves
\begin{enumerate}
\item total mass:
\begin{equation}\label{eq:conservation}
\forall n \in \NN, \quad \sum_{J \in \ZZ^d} \rho_{J}^n=1;
\end{equation}
\item positivity of the density and the bound on the velocity field:
\[
\forall (n,J) \in \NN \times \ZZ^d, \, \forall i=1,\ldots, d, \quad \rho_{J}^n \geq 0, \qquad |{a_i}_{J}^n| \leq a_\infty,
\]
under the CFL condition:
\begin{equation}\label{CFLlarge}
\sum_{i=1}^d \Bigl(a_\infty \frac{\dt}{\dx_i} + 2\eps \frac{\dt}{\dx_i^2}\Bigr) \leq 1;
\end{equation}
\item the center of mass:
\[
\forall n\in \NN^*, \rho_{\Delta x}^{\eps,n} \in \mathcal{P}_1(\RR^d) \quad\text{and }
\quad \sum_{J\in \ZZ^d} x_J \rho_{J}^n = \sum_{J\in \ZZ^d} x_J \rho_{J}^0.
\]
\end{enumerate}
\end{lemma}
\begin{proof}
The first item comes from summing equation \eqref{dis:explicit} over $J \in \ZZ^d$. Moreover, using the following rewriting of $\rho_{J}^{n+1}$ as a positive combination of $\rho_J$ and $\rho_{J\pm e_i}$, $i=1,\ldots,d$:\vspace*{-3pt}
\begin{multline*}
\rho_{J}^{n+1} = \rho_{J}^n \biggl[ 1 - \sum_{i=1}^d \Bigl(\frac{\dt}{\dx_i} |{a_i}^n_{J}| + \frac{2\eps\dt}{\dx_i^2}\Bigr) \biggr]
+ \sum_{i=1}^d \rho_{J+e_i}^n \Bigl(\frac{\dt}{\dx_i}({a_i}^n_{J+e_i})^- + \frac{\eps\dt}{\dx_i^2}\Bigr)\notag\\[-5pt] \\[-8pt]
+ \sum_{i=1}^d \rho_{J-e_i}^n \Bigl(\frac{\dt}{\dx_i}({a_i}^n_{J-e_i})^+ + \frac{\eps\dt}{\dx_i^2}\Bigr), \notag
%\label{schemarho}
\end{multline*}
it is classical to prove the second item by induction on $n \in \NN$, under the CFL condition~\eqref{CFLlarge} under which $\rho_J^{n+1}$ is a convex combination of the $\rho_K^n$, see \cite{leveque} for example.
\item Let us now focus on the third item. One has\vspace*{-5pt}
\begin{multline*}
\sum_{J\in \ZZ^d} |x_J|\rho_{J}^{n+1} \\[-3pt]
=
\sum_{J\in \ZZ^d} |x_J|\biggl[ \rho_{J}^n- \sum_{i=1}^{d} \frac{\dt}{\dx_i}
\Big(({a_i}^n_{J})^+ \rho_{J}^n - ({a_i}^n_{J+e_i})^- \rho_{J+e_i}^n
-({a_i}^n_{J-e_i})^+ \rho_{J-e_i}^n + ({a_i}^n_{J})^- \rho_{J}^n \Big) \\[-3pt]
\eps \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^n - 2\rho_J^n + \rho_{J-e_i}^n \right)\biggr],
\end{multline*}
thus\vspace*{-5pt}
\begin{multline*}
\sum_{J\in \ZZ^d} |x_J|\rho_{J}^{n+1} \leq \sum_{J\in \ZZ^d} |x_J|\rho_{J}^{n}
\Bigl( 1 + \sum_{i = 1}^d (a_\infty \frac{\Delta t}{\Delta x_i} + 2 \eps \frac{\Delta t}{\Delta x_i^2})\Bigr) \\[-3pt]
+ \sum_{i = 1}^d \sum_{J\in \ZZ^d} |x_{J+e_i}|\rho_{J + e_i}^{n} \Bigl(a_\infty \frac{\Delta t}{\Delta x_i} + \eps \frac{\Delta t}{\Delta x_i^2}\Bigr)
+ \sum_{i = 1}^d \sum_{J\in \ZZ^d} \Delta x_i\rho_{J + e_i}^{n} \Bigl(a_\infty \frac{\Delta t}{\Delta x_i} + \eps \frac{\Delta t}{\Delta x_i^2}\Bigr) \\[-3pt]
+ \sum_{i = 1}^d \sum_{J\in \ZZ^d} |x_{J-e_i}|\rho_{J - e_i}^{n} \Bigl(a_\infty \frac{\Delta t}{\Delta x_i} + \eps \frac{\Delta t}{\Delta x_i^2}\Bigr)
+ \sum_{i = 1}^d \sum_{J\in \ZZ^d} \Delta x_i\rho_{J - e_i}^{n} \Bigl(a_\infty \frac{\Delta t}{\Delta x_i} + \eps \frac{\Delta t}{\Delta x_i^2}\Bigr),
\end{multline*}
which shows by induction that $\rho_{\Delta x}^{\eps,n} \in \mathcal{P}_1(\RR^d)$ if $\rho_{\Delta x}^{\eps,0} \in \mathcal{P}_1(\RR^d)$.
Now more precisely, using the discretization \eqref{dis:explicit} together with a discrete integration by parts, we~have:\vspace*{-5pt}
\begin{multline*}
\sum_{J\in \ZZ^d} x_J\rho_{J}^{n+1} \\[-3pt]
= \sum_{J\in
\ZZ^d} x_J\rho_{J}^n 
- \sum_{i=1}^d\frac{\dt}{\dx_i} \sum_{J\in
\ZZ^d} \left(({a_i}^n_{J})^+ \, \rho_{J}^n \big(x_J-x_{J+e_i}\big)
-({a_i}^n_{J})^- \, \rho_{J}^n \big(x_{J-e_i}-x_J\big) \right) \\[-3pt]
+ \eps\sum_{i=1}^d\frac{\dt}{\dx_i^2} \sum_{J\in
\ZZ^d} \big(x_{J-e_i}-x_J\big)(\rho_{J+e_i}^n - \rho_{J}^n).
\end{multline*}
By definition of $x_J$, we~have $x_{J-e_i}-x_J = -\dx_i$. Hence, we~deduce:\vspace*{-5pt}
\begin{multline*}
\sum_{J\in \ZZ^d} x_J\rho_{J}^{n+1} = \sum_{J\in \ZZ^d}
x_J\rho_{J}^n + \dt \sum_{i=1}^d \sum_{J\in \ZZ^d} {a_i}^n_{J} \, \rho_{J}^n - \eps\sum_{i=1}^d\frac{\dt}{\dx_i} \sum_{J\in
\ZZ^d}(\rho_{J+e_i}^n - \rho_{J}^n) \\[-3pt]
= \sum_{J\in \ZZ^d}
x_J\rho_{J}^n + \dt \sum_{i=1}^d \sum_{J\in \ZZ^d} {a_i}^n_{J} \, \rho_{J}^n.
\end{multline*}
Applying the symmetrization Lemma \ref{lemma:symmetrizationDis} to the constant vector field given by $v_J = (1, \ldots, 1) \in \RR^d$ for all $J \in \ZZ^d$, we~have $\sum_{J\in \ZZ^d} {a_i}^n_{J} \, \rho_{J}^n = 0$ for all $i=1,\ldots,d$, hence the result.
\end{proof}
The following lemma ensures that $M_2(\rho_{\dx}^{\eps,n})$ remains bounded over finite time. It turns out necessary for the proof of convergence of the scheme by compactness, in order to extract a converging subsequence.
\begin{lemma}[Bound on the second moment]\label{bounddismom}
For all $n\in \NN^*$, the following estimate holds:\vspace*{-3pt}
\begin{equation*}
M_{2,\dx}^n := \sum_{J\in \ZZ^d} |x_J|^2\rho_{J}^n \leq e^{-4 \lambda t^n} \Big(M_{2,\dx}^0 + a_\infty t^n \sum_{i=1}^d \dx_i + 2d\eps t^n\Big).
\end{equation*}
\end{lemma}

\begin{proof}
Using \eqref{dis:explicit} and a discrete integration by parts,
one can write:\vspace*{-3pt}
\begin{multline*}
\sum_{J\in \ZZ^d} |x_J|^2 \rho_{J}^{n+1} =\sum_{J\in
\ZZ^d} |x_J|^2\rho_{J}^n\\[-3pt]
- \sum_{i=1}^d
\frac{\dt}{\dx_i} \sum_{J\in
\ZZ^d} \Bigl[ ({a_i}^n_{J})^+ \, \rho_{J}^n \big(|x_J|^2-|x_{J+e_i}|^2\big)
-({a_i}^n_{J})^- \, \rho_{J}^n \big(|x_{J-e_i}|^2-|x_{J}|^2\big)
\Bigr]
\\[-3pt]
+ \eps \sum_{i=1}^d
\frac{\dt}{\dx_i} \sum_{J\in
\ZZ^d}\big(|x_{J-e_i}|^2 - |x_J|^2\big)\big(\rho_{J+e_i} - \rho_J\big).
\end{multline*}
By definition of $x_J$,\vspace*{-3pt}
\[
|x_J|^2-|x_{J+e_i}|^2=-2J_i\, \dx_i^2 - \dx_i^2\quad\text{and}\quad
|x_{J-e_i}|^2-|x_J|^2=-2J_i\, \dx_i^2 + \dx_i^2.
\]
Therefore, we~get:\vspace*{-3pt}
\begin{multline*}
\sum_{J\in \ZZ^d} |x_J|^2 \rho_{J}^{n+1} = \sum_{J\in \ZZ^d}
|x_J|^2 \rho_{J}^{n} + 2\dt \sum_{i=1}^d\sum_{J\in \ZZ^d}
J_i \dx_i \, {a_i}^n_{J} \, \rho_{J}^n + \dt\sum_{i=1}^d \dx_i \sum_{J\in \ZZ^d} \rho_{J}^n |{a_i}^n_{J}| \\
+ \eps \dt\sum_{i=1}^d \sum_{J\in \ZZ^d} (-2J_i + 1) \dx_i(\rho_{J+e_i} - \rho_J).
\end{multline*}
As a consequence of Lemma \ref{lem:CFL}, we~have $|{a_i}^n_{J}|\leq a_\infty$.
Using, in addition, the mass conservation, we~deduce that the penultimate term is bounded by $a_\infty \dt \sum_{i=1}^d \dx_i$.
As for the last term, another integration by parts shows that the last term equals $2d\eps \dt$. Finally, Lemma \ref{lemma:symmetrizationDis} applied to the discrete vector field given by $v_J = x_J$ yields:\vspace*{-3pt}
\begin{align*}
2\dt \sum_{i=1}^d\sum_{J\in \ZZ^d}
J_i \dx_i \, {a_i}^n_{J} \, \rho_{J}^n = 2\dt \sum_{J\in \ZZ^d} x_J \cdot a_j^n \, \rho_{J}^n & = - \dt \sum_{J, K \in \ZZ^d}\hspace*{-2mm} DW_J^K \cdot (x_J - x_K) \rho_J^n \rho_K^n\\
& \leq - \lambda \dt \sum_{J, K \in \ZZ^d} |x_J - x_K|^2 \rho_J^n \rho_K^n \\
& \leq - 4\lambda \dt \sum_{J\in \ZZ^d} |x_J|^2 \rho_{J}^{n},
\end{align*}
where we used the $\lambda$-convexity of $W$ and the inequality $|x_J - x_K|^2 \leq 2(|x_J|^2 + |x_K|^2)$ along with the fact that $\lambda$ is nonpositive. We~obtain\vspace*{-3pt}
\[
\sum_{J\in \ZZ^d} |x_J|^2 \rho_{J}^{n+1} \leq
(1-4\lambda\dt)\sum_{J\in \ZZ^d} |x_J|^2 \rho_{J}^{n}
+ a_\infty\dt\sum_{i=1}^d\dx_i + 2d\eps \dt.
\]
We conclude the proof using a discrete version of Grönwall's lemma.
\end{proof}

\subsubsection{Proof of Proposition \ref{prop:cvOrder1/2:scheme}}

Before going through the proof of Proposition \ref{prop:cvOrder1/2:scheme}, let us introduce, for $J\in\ZZ^d$ and $y \in \RR^d$ the following coefficients:
\begin{multline}\label{def:alpha}
\alpha_J(y)\\
=
\begin{cases}
\ds 1-\sum_{i=1}^d \Bigl(\frac{|\langle y-{x_J},e_i\rangle|}{\dx_i} - \frac{2\eps\dt}{\dx_i^2}\Bigr)
& \textup{when} \ y\in C_J,
\vspace{5pt}
\\
\ds \frac{1}{\dx_i}\bigl(\langle y-x_{J-e_i},e_i\rangle\bigr)^+ + \frac{\eps\dt}{\dx_i^2}
&\textup{when} \ y\in C_{J-e_i}, \ \ \textup{for} \ i= 1,\dots,d,
\vspace{5pt}
\\
\ds \frac{1}{\dx_i}\bigl(\langle y-x_{J+e_i},e_i\rangle\bigr)^- + \frac{\eps\dt}{\dx_i^2}
&\textup{when} \ y\in C_{J+e_i}, \ \ \textup{for} \ i=1,\dots,d,
\vspace{5pt}
\\
0 &\textup{otherwise}.
\end{cases}
\end{multline}
It then holds that, for any $J,L \in {\mathbb Z}^d$,
\begin{multline*}%\label{eq:transition:probas}
\alpha_{J} \bigl( x_{L} + \dt a^n_{L}\bigr)\\ =
\begin{cases}
\ds 1-\sum_{i=1}^d \Bigl(\vert {a_{i}}^n_{J} \vert \frac{\dt}{\dx_i} - \frac{2\eps\dt}{\dx_i^2}\Bigr)
&\textup{ when } \ L=J,
\vspace{5pt}
\\
\ds \frac{\dt}{\dx_i}\bigl( {a_{i}}^n_{J-e_{i}} \bigr)^{+} + \frac{\eps\dt}{\dx_i^2}
&\textup{ when } \ L = J-e_{i}, \ \ \textup{for} \ i= 1,\dots,d,
\vspace{5pt}
\\
\ds \frac{\dt}{\dx_i}\bigl( {a_{i}}^n_{J+e_{i}} \bigr)^{-} + \frac{\eps\dt}{\dx_i^2}
&\textup{ when } \ L = J+e_{i}, \ \ \textup{for} \ i= 1,\dots,d,
\vspace{5pt}
\\
0 & \textup{ otherwise},
\end{cases}
\end{multline*}
so that we have the key identity:
\begin{equation*}%\label{scheme3}
\forall\, J \in \ZZ^d, \quad
\rho^{n+1}_J = \sum_{L\in\ZZ^d} \rho^n_L \alpha_{J}\bigl(x_L+\dt a^n_L\bigr),
\end{equation*}
\begin{lemma}\label{lemma:alpha}
For any $y \in \RR^d$, we~have
\[
\sum_{L\in \ZZ^d} \alpha_{L}(y) = 1 \quad \mbox{ and } \quad \sum_{L\in \ZZ^d} x_L \alpha_{L}(y) = y.
\]
\end{lemma}
\begin{proof}
Let $J\in \ZZ^d$ such that $y\in C_J$. To prove the first claim, we~just use the definition of $\alpha_{L}(y)$:
\begin{multline*}
\sum_{L\in \ZZ^d} \alpha_{L}(y) = \alpha_{J}(y) + \sum_{i=1}^d \bigl(\alpha_{J+e_i}(y)+\alpha_{J-e_i}(y)\bigr)
\\
= 1 - \sum_{i=1}^d \Bigl(\frac{|\langle y-{x_L},e_i\rangle|}{\dx_i} - \frac{2\eps\dt}{\dx_i^2}\Bigr)
+ \sum_{i=1}^d \frac{1}{\dx_i} \left(\bigl(\langle y-x_{J},e_i\rangle \bigr)^+ + \bigl(\langle y-x_{J},e_i\rangle \bigr)^-\right) \\
+ 2 \sum_{i=1}^d \frac{\eps\dt}{\dx_i^2} = 1.
\end{multline*}
As for the preservation of the barycenter, we~once again using the definition of the coefficients $\alpha_{L}(y)$:
\begin{align*}
\sum_{L\in \ZZ^d}& x_L \alpha_{L}(y) = x_J \alpha_J(y) + \sum_{i=1}^d x_{J+e_i} \alpha_{J+e_i}(y)+ \sum_{i=1}^d x_{J-e_i} \alpha_{J-e_i}(y)
\\
&\hspace*{-4mm}= x_J \biggl[1-\sum_{i=1}^d \Bigl(\frac{|\langle y-{x_J},e_i\rangle|}{\dx_i} - \frac{2\eps\dt}{\dx_i^2}\Bigr)\biggr] + \sum_{i=1}^d x_{J} \Bigl(\frac{1}{\dx_i}\bigl(\langle y-x_{J},e_i\rangle\bigr)^+ + \frac{\eps\dt}{\dx_i^2}\Bigr) \\
&\hspace*{4cm}+ \sum_{i=1}^d x_{J} \Bigl(\frac{1}{\dx_i}\bigl(\langle y-x_{J},e_i\rangle\bigr)^- + \frac{\eps\dt}{\dx_i^2} \Bigr) \\
&\hspace*{-4mm}= x_J \biggl[1+\sum_{i=1}^d \Bigl(-\frac{|\langle y-{x_J},e_i\rangle|}{\dx_i} + \frac{1}{\dx_i}\bigl(\langle y-x_{J},e_i\rangle\bigr)^+ + \frac{1}{\dx_i}\bigl(\langle y-x_{J},e_i\rangle\bigr)^- \Bigr)\biggr] \\
&\hspace*{4cm}+ \sum_{i=1}^d e_i \left(\bigl(\langle y-x_J, e_i\rangle \bigr)^+ - \bigl(\langle y-x_J, e_i\rangle \bigr)^- \right) \\
&\hspace*{-4mm}= x_J + \sum_{i=1}^d \langle y-x_J, e_i\rangle e_i = y.\qedhere
\end{align*}
\end{proof}

We now turn to the proof of Proposition \ref{prop:cvOrder1/2:scheme}.
For $n \in \NN^*$, we~denote $D_n:= W_2 \bigl(\rho_{t^n},\rho_{\dx}^{\eps,n} \bigr)$. The point is, roughly speaking, to obtain an estimate of the type $D_{n+1}^2 \leq D_n^2 + C\dt(\dt + \dx + \eps)$ and then use a discrete Grönwall lemma to obtain estimate \eqref{eq:cvOrder1/2:scheme:bound:1}.

Let $\gamma$ be an optimal transport plan between $\rho_{t^n}$ and $\rho_{\dx}^{\eps,n}$, so that
\[
D_n^2 =\iint |x-y|^2 \gamma(dx,dy).
\]
We~also let $a^n_{\dx}$ be any continuous reconstruction of the discrete velocity defined in \eqref{def:aij}, for instance piecewise affine, such that $a_{\dx}^n(x_J) = a_J^n$ for all $J \in \ZZ^d$.

To construct an adequate coupling $\gamma' \in \Gamma\big(\rho_{t^{n+1}}, \rho_{\dx}^{\eps,n+1}\big)$, recall that Theorem \ref{thm:cjlv15} gives $\rho_{t^{n+1}} = Z(t^{n+1}, t^n, \cdot) \# \rho_{t^n}$, where $Z$ is the Filippov characteristic flow associated to $\achapo[\rho]$ given by Theorem \ref{thm:cjlv15}, except that here the second variable of $Z$ denotes the time of the Cauchy data (which is the third variable) whereas in Theorem \ref{thm:cjlv15} it was omitted as it was $0$. If the discrete measure $\rho_{\dx}^{\eps,n+1}$ was a pushforward measure of $\rho_{\dx}^{\eps,n}$, we~would also define $\gamma'$ as a pushforward of $\gamma$, but it is not the case in general as we are dealing with atomic measures. Instead, if we denote by $\nu$ the kernel on $(\RR^d, \calB(\RR^d))$ given by:
\[
\forall (y, B) \in \RR^d \times \calB(\RR^d), \qquad \nu(y, B) = \sum_{J \in \ZZ^d} \alpha_J(y + \dt a^n_{\Delta_x}(y)) \delta_{x_J}(B),
\]
we have the kernel representation:
\[
\forall B \in \calB(\RR^d), \qquad \rho_{\dx}^{\eps,n+1}(B) = \int \nu(y, B) \rho_{\dx}^{\eps,n}(dy).
\]
The pushforward $\rho_{t^{n+1}} = Z(t^{n+1}, t^n, \cdot) \# \rho_{t^n}$ can also be seen as a kernel representation. Indeed, setting $\mu(x, A) = \delta_{Z(t^{n+1}, t^n, x)}(A)$ for $(x,A) \in \RR^d \times \calB(\RR^d)$, we~have:
\begin{align*}
\forall A \in \calB(\RR^d), \quad \rho_{t^{n+1}}(A) &= \int {\mathbf 1}_{A}(Z(t^{n+1}, t^n, x)) \rho_{t^n} (dx) \\
&= \int \delta_{Z(t^{n+1}, t^n, x)}(A) \rho_{t^n} (dx) = \int \mu(x, A) \rho_{t^n} (dx).
\end{align*}
We now define the product kernel ${\mathcal K}$ on $\big(\RR^d \times \RR^d \big) \times \calB\big(\RR^d \times \RR^d \big)$ by:
\[
{\mathcal K} \bigl( (x,y), A \times B \bigr) = \mu(x,A) \nu(y,B) = \delta_{Z(t^{n+1}, t^n, x)}(A)
\sum_{L \in {\ZZ}^d} \alpha_{L}\bigl(y+\dt a^n_{\dx}(y) \bigr)
\delta_{x_L}(B)
\]
and then set $\gamma'(A \times B) = \iint_{\RR^d \times \RR^d} {\mathcal K} \bigl( (x,y),A \times B \bigr)
\gamma(dx,dy)$. Equivalently, for any $\theta \in \calC_b(\RR^d \times \RR^d)$, we~have:
\begin{align*}
\iint \theta(x,y) \gamma'(dx,dy) & = \intquad \theta(x',y') \mu(x, dx') \nu(y, dy') \gamma(dx,dy) \\
& = \iint \biggl[\sum_{L\in \ZZ^d} \theta\bigl(Z(t^{n+1};t^n,x),x_L\bigr)
\alpha_L\bigl(y+\dt a^n_{\dx}(y)\bigr)\biggr]\,\gamma(dx,dy).
\end{align*}
One can show as in \cite{dlv20} that the marginals of $\gamma'$ are $\rho_{t^{n+1}}$ and $\rho_{\dx}^{\eps,n+1}$. In particular, we~have:
\[
D_{n+1}^2 \leq \iint |x-y|^2 \gamma'(dx,dy).
\]
Using the definition of $\gamma'$, we~get:
\begin{equation}\label{eq1}
D_{n+1}^2 \leq \iint \sum_{L\in \ZZ^d} \bigl|Z(t^{n+1};t^n,x)-x_L\bigr|^2 \alpha_L \bigl(y+\dt a^n_{\dx}(y)\bigr) \gamma(dx,dy).
\end{equation}
Writing
\[
Z(t^{n+1};t^n,x)-x_L = Z(t^{n+1};t^n,x) - (y+\dt a^n_{\dx}(y)) - \big(x_L - (y+\dt a^n_{\dx}(y))\big)
\]
and expanding the square, we~obtain:\vspace*{-3pt}
\begin{multline}\label{eqZz}
\sum_{L\in \ZZ^d}\bigl|Z(t^{n+1};t^n,x) - x_L\bigr|^2 \alpha_L\bigl(y+\dt a^n_{\dx}(y)\bigr)\\
= \bigl|Z(t^{n+1};t^n,x) - y - \dt a^n_{\dx}(y) \bigr|^2 + \sum_{L\in \ZZ^d} \bigl|x_L - y - \dt a^n_{\dx}(y)
\bigr|^2 \alpha_L\bigl(y+\dt a^n_{\dx}(y)\bigr)
\\
-2\Bigl( Z(t^{n+1};t^n,x) - y - \dt a^n_{\dx}(y) \Bigr) \cdot
\biggl( \sum_{L \in \ZZ^d}
\bigl(x_L-y-\dt a^n_{\dx}(y)\bigr) \alpha_L\bigl(y+\dt a^n_{\dx}(y) \biggr).
\end{multline}
Now, the last term in scalar product vanishes as we have, using Lemma \ref{lemma:alpha},\vspace*{-3pt}
\[
\sum_{L\in \ZZ^d} \bigl(x_L-y-\dt a^n_{\dx}(y)\bigr) \alpha_L\bigl(y+\dt a^n_{\dx}(y)\bigr) = y+\dt a^n_{\dx}(y) - \big(y+\dt a^n_{\dx}(y)\big) = 0.
\]
Plugging \eqref{eqZz} into \eqref{eq1}, we~therefore deduce, using the fact that $\rho_{\dx}^{\eps,n}$ is the second marginal of $\gamma$:\vspace*{-3pt}
\begin{multline}\label{eq2}
D_{n+1}^2 \leq
\iint \bigl|Z(t^{n+1};t^n,x)-y-\dt a^n_{\dx}(y)\bigr|^2 \gamma(dx,dy) \\[-3pt]
+ \int \sum_{L\in \ZZ^d} \bigl|x_L-y-\dt a^n_{\dx}(y) \bigr|^2 \alpha_L\bigl(y+\dt a^n_{\dx}(y)\bigr) \rho_{\dx}^{\eps,n}(dy),
\end{multline}
Let us deal with the last term in the above inequality.
We have\vspace*{-3pt}
\[
\rho_{\dx}^{\eps,n}(y) = \sum_{J\in\ZZ^d} \rho_J^n \delta_{x_J}(y),
\]
therefore\vspace*{-3pt}
\begin{align*}
&\sum_{L\in \ZZ^d} \int \bigl|x_L-y-\dt a^n_{\dx}(y)\bigr|^2 \alpha_L\bigl(y+\dt a^n_{\dx}(y)\bigr) \rho_{\dx}^{\eps,n}(dy) \\[-5pt]
&\hspace{5cm} = \sum_{J\in \ZZ^d} \sum_{L\in \ZZ^d} \bigl|x_L-x_J-\dt a^n_J\bigr|^2 \alpha_L\bigl(x_J+\dt a^n_J\bigr) \rho_J^n.
\end{align*}
Moreover, using the definition of $\alpha_L$ in \eqref{def:alpha}, we~compute:\vspace*{-5pt}
\begin{multline*}
\sum_{L\in \ZZ^d} \bigl|x_L-x_J-\dt a^n_J\bigr|^2 \alpha_L\bigl(x_J+\dt a^n_J\bigr)
= \dt^2 |a_J^n|^2 \biggl(1-\sum_{i=1}^d \frac{\dt}{\dx_i} |{a_i}_J^n| - \sum_{i=1}^d \frac{2\eps\dt}{\dx_i^2}\biggr) \\
+ \sum_{i=1}^d \Bigl[ \bigl|\dx_i e_i -\dt a^n_J\bigr|^2 \Big(\frac{\dt}{\dx_i} ({a_i}_J^n)^+ + \frac{\eps\dt}{\dx_i^2}\Big)+
\bigl|\dx_i e_i +\dt a^n_J\bigr|^2 \Big(\frac{\dt}{\dx_i} ({a_i}_J^n)^- + \frac{\eps\dt}{\dx_i^2}\Big) \Bigr] \\
\leq C\dt(\dt + \dx + \eps),
\end{multline*}
where we used, for the last inequality, the CFL condition \eqref{CFL} and the fact that the velocity $a_J^n$ is uniformly bounded. Multiplying by $\rho_J^n$, summing over $J \in \ZZ^d$, and injecting into \eqref{eq2} yields:
\begin{equation}\label{eqD1}
D_{n+1}^2 \leq \iint \bigl|Z(t^{n+1};t^n,x)-y-\dt a^n_{\dx}(y)\bigr|^2 \gamma(dx,dy) + C \dt (\dt+\dx+\eps).
\end{equation}

Dealing with the first term amounts to estimating the distance between the exact characteristics $Z(t^{n+1};t^n,x)$ and the forward Euler discretization $y+\dt a^n_{\dx}(y)$. To~this end, we~write, on the one hand, using the definition of the Filippov characteristics \cite{filippov88, poupaud_rascle97}:
\begin{align*}
Z(t^{n+1};t^n,x) & = x + \int_{t^n}^{t^{n+1}} \achapo_\rho\bigl(s,Z(s;t^n,x)\bigr) ds
\\
& = x - \int_{t^n}^{t^{n+1}} \int\nabWchapo\bigl(Z(s;t^n,x)-Z(s;t^n,\xi)\bigr) \rho_{t^n}(d\xi) ds.
\end{align*}
On the other hand, we~have, whenever $y$ is a node of the mesh,
\begin{equation*}
y + \dt a^n_{\dx}(y)
= y - \dt \int
\nabWchapo( y - \zeta) \rho^n_{\dx}( d \zeta \bigr).
\end{equation*}
Thus, still for $y$ a node of the mesh, we~have:
\begin{multline*}
\bigl|Z(t^{n+1};t^n,x)-y-\dt a^n_{\dx}(y)\bigr|^2 \leq |x-y|^2 \\
- 2\int_{t^n}^{t^{n+1}}\!\!\iint(x-y) \cdot \Big(\nabWchapo\big(Z(s;t^n,x)-Z(s;t^n,\xi)\big)-\nabWchapo(y-\zeta)\Big) \rho_{t^n}(d\xi) \rho_{\dx}^{\eps,n}(d\zeta) \\
+ C\dt^2.
\end{multline*}
Since $\gamma\in \Gamma(\rho_{t^n},\rho_{\dx}^{\eps,n})$ and since the above integral can be decoupled using the linearity of the scalar product, we~also have:
\begin{align*}
\iint(x-y) \cdot & \Big(\nabWchapo\big(Z(s;t^n,x)-Z(s;t^n,\xi)\big)-\nabWchapo(y-\zeta)\Big) \rho_{t^n}(d\xi) \rho_{\dx}^{\eps,n}(d\zeta) \\
& = \iint(x-y) \cdot \Big(\nabWchapo\big(Z(s;t^n,x)-Z(s;t^n,\xi)\big)-\nabWchapo(y-\zeta)\Big) \gamma(d\xi,d\zeta).
\end{align*}
Injecting into \eqref{eqD1}, we~get:
\begin{multline*}
D_{n+1}^2 \leq D_n^2 + C \dt(\dt+\dx+\eps) \\
- 2 \int_{t^n}^{t^{n+1}} \!\!\! \intquad \!\!\Big(x-y\Big) \cdot \Big(\nabWchapo\big(Z(s;t^n,x)-Z(s;t^n,\xi)\big)-\nabWchapo(y-\zeta)\Big) \gamma(d\xi,d\zeta) \gamma(dx,dy).
\end{multline*}
Decomposing $x-y=x-Z(s;t^n,x)+Z(s;t^n,x)-y$ and using the fact that \\
$|Z(s;t^n,x)-x|\leq a_\infty |s-t^n|$, we~get:
\begin{multline*}
D_{n+1}^2 \leq D_n^2 + C \dt(\dt+\dx+\eps)
- 2 \int_{t^n}^{t^{n+1}}\intquad \Big( Z(s;t^n,x)-y\Big) \cdot \\
\Big( \nabWchapo\bigl(Z(s;t^n,x) - Z(s;t^n,\xi)\bigr)-\nabWchapo(y-\zeta)\Big)
\gamma(d\xi,d\zeta) \gamma(dx,dy).
\end{multline*}
Using the fact that $W$ is even to symmetrize the last term as in Lemma \ref{lemma:symmetrization}, we~obtain:
\begin{multline*}
D_{n+1}^2 \leq \ D_n^2 + C \dt(\dt+\dx+\eps) \\
- \int_{t^n}^{t^{n+1}}\intquad\Big( Z(s;t^n,x)-Z(s;t^n,\xi)-y+\zeta \Big) \cdot \\
\Big(\nabWchapo\bigl(Z(s;t^n,x)-Z(s;t^n,\xi)\bigr)-\nabWchapo(y-\zeta)\Big) \, \gamma(d\xi,d\zeta) \gamma(dx,dy).
\end{multline*}
The $\lambda$-convexity of $W$ then yields:
\begin{multline*}
D_{n+1}^2 \leq D_n^2 + C \dt(\dt+\dx+\eps) \\
- \lambda \int_{t^n}^{t^{n+1}}\intquad \bigl|Z(s;t^n,x)-y-Z(s;t^n,\xi)+\zeta\bigr|^2
\, \gamma(d\xi,d\zeta) \gamma(dx,dy).
\end{multline*}
Expanding the last term gives:
\begin{multline}\label{interm2}
D_{n+1}^2 \leq \ D_n^2 + C \dt(\dt+\dx+\eps)
- 2\lambda \int_{t^n}^{t^{n+1}} \iint \bigl|Z(s;t^n,x)-y\bigr|^2
\,\gamma(dx,dy) \\
+ 2\lambda \int_{t^n}^{t^{n+1}} \biggl\vert \iint \bigl(Z(s;t^n,x)-y\bigr)
\,\gamma(dx,dy) \biggr\vert^2. 
\end{multline}
Now, since $\lambda \leq 0$, the last term above is nonpositive. It remains to estimate the penultimate term. Writing:
\[
\bigl|Z(s;t^n,x)-y\bigr| \leq \bigl|Z(s;t^n,x)-x\bigr| + \bigl|x-y\bigr| \leq a_\infty\bigl|s-t^n\bigr| + \bigl|x-y\bigr|,
\]
we deduce:
\[
\bigl|Z(s;t^n,x)-y\bigr|^2 \leq 2\big(a_\infty^2\bigl|s-t^n\bigr|^2 + \bigl|x-y\bigr|^2\big) \leq 2a_\infty^2 \dt^2 + 2\bigl|x-y\bigr|^2
\]
whenever $s \in [t^n, t^{n+1}]$. Integrating in space with respect to $\gamma(dx, dy)$ and integrating over $s \in [t^n, t^{n+1}]$, we~obtain:
\[
- 2\lambda \int_{t^n}^{t^{n+1}} \iint \bigl|Z(s;t^n,x)-y\bigr|^2
\,\gamma(dx,dy) \leq -4\lambda a_\infty^2 \dt^3 -4\lambda \dt D_n^2.
\]
Together with \eqref{interm2}, this yields:
\[
D_{n+1}^2 \leq (1-4\lambda\dt)D_n^2 + C \dt(\dt+\dx+\eps).
\]
Using a discrete Grönwall lemma, we~finally get:
\[
D_n^2 \leq e^{-4\lambda t^n} D_0^2 + C\, \frac{1 - e^{-4\lambda t^n} }{4\lambda}\, (\dt+\dx+\eps).
\]
Now, one can easily prove that $D_0^2 = W_2^2(\rho^{\ini}, \rho_{\dx}^0) \leq \frac{\sqrt{d}}{2} \dx^2$ (see \cite{dlv17}). This, along with the CFL condition \eqref{CFL}, which implies that $\dt \leq C\dx$, gives the desired result.

\subsubsection{Proof of Theorem \ref{thm:cvOrder1/2}}

We are now in position to prove Theorem \ref{thm:cvOrder1/2} using estimate \eqref{eq:cvOrder1/2:scheme:bound:1} and passing to the limit $\dx \to 0$. To do so, we~must verify that, for a given $\eps >0$, the approximate solutions given by the numerical scheme \eqref{dis:explicit}--\eqref{dis:rho0} converge, say uniformly in time (over a finite horizon) and weakly, in the sense of measures, in space, towards the solution $\rho^\eps$ to the aggregation-diffusion problem \eqref{eq:agg_diff} with initial datum $\rho^{\ini}$, as $\dx \to 0$. In all this section, $\eps$ is a fixed positive real number.

Let $T > 0$ and let $N \in \NN$ be such that $N\dt = T$ where $\Delta t$ satisfies the CFL condition. We~consider the following piecewise affine reconstruction in time, defined for $t \in [0,T]$ by
\begin{subequations}
\begin{align}
& \rho_{\dx, t}^{\eps} := \sum_{n=0}^N \Bigl(\frac{t^{n+1} - t}{\dt} \rho_{\dx}^{\eps,n} + \frac{t-t^n}{\dt}\rho_{\dx}^{\eps,n+1}\Bigr) {\mathbf 1}_{[t^n, t^{n+1}[}(t), \label{rhoaff}\\
& \rho_{\dx}^{\eps,n} := \sum_{J\in \ZZ^d} \rho_J^n \delta_{x_J}, \quad n = 0, \ldots, N,
\end{align}
\end{subequations}
where, once again, $(\rho_J^n)_{J\in\ZZ^d}^{n=0,\ldots,N}$ is given by the explicit discretization \eqref{dis:explicit}--\eqref{dis:rho0} (it~actually depends on $\eps$ but we drop this dependence for convenience). Lemmas \ref{lem:CFL} and \ref{bounddismom} show that, for all $n \in \{0,\ldots,N\}$, $\rho_{\dx}^{\eps,n} \in \calP_2(\RR^d)$, hence $(\rho_{\dx}^\eps)_{\dx > 0}$ is a collection, indexed by $\dx$, of curves in $\calC([0,T], \WW_1(\RR^d))$ (they are actually curves on $\WW_2(\RR^d)$ but compactness arguments require to work in a smaller space).

\subsubsection*{Outline of the proof}
We begin with assuming that $\rho^{\ini} \in L^2(\RR^d)$. Then, from $(\rho_{\dx}^\eps)_{\dx > 0}$, we~shall extract a subsequence, that we still denote $(\rho_{\dx}^\eps)_{\dx > 0}$, converging in the $\calC([0,T], \calM_b(\RR^d))$ topology towards a limit $\rho \in \calC([0,T], \WW_2(\RR^d))$. To do so, we~apply the Ascoli-Arzelà theorem: the relative compactness assumption follows quite directly from the uniform bound on $M_2(\rho_{\dx}^{\eps,n})$ that we proved in Lemma \ref{bounddismom}; the equicontinuity assumption, however, is more involved and requires discrete $H^1$ estimates (Lemma \ref{lemma:DiscreteH1Estimates}) in order to control the diffusive term. Then, using classical discrete integration by parts, we~show that $\rho$ solves the aggregation-diffusion initial value problem, the solution of which is unique, hence the whole sequence actually converges. Passing to the limit $\dx \to 0$ in estimate \eqref{eq:cvOrder1/2:scheme:bound:1} will give us the desired estimate \eqref{eq:cvOrder1/2:scheme:bound:2} for $L^2(\RR^d)$ initial datum, and it will only remain to use a regularization argument to conclude in the case of arbitrary $\calP_2(\RR^d)$ initial datum.
\begin{lemma}\label{lemma:DiscreteH1Estimates}
For all $m \in \{0, \ldots, N\}$, we~have:
\[
\dt\sum_{n=0}^{m-1} \sum_{J \in \ZZ^d} \sum_{i=1}^d \Big|\frac{\rho_{J+e_i}^n - \rho_J^n}{\dx_i}\Big|^2 \leq C(a_\infty, d, \eps,T) \sum_{J\in\ZZ^d} \frac{\big(\rho_J^0\big)^2}{2},
\]
with
\[
C(a_\infty, d, \eps,T) =
\frac{1}{2 \eps}\biggl(1 + \frac{8dT a_\infty^2}{\eps} \sum_{J\in\ZZ^d}
\exp\Bigl(\frac{4(1 + d)T a_\infty^2}{\eps}\Bigr)\biggr).
\]
\end{lemma}
\begin{proof}
The idea is to perform a discrete version of the following rationale. If $\rho^\eps$ solves \eqref{eq:agg_diff} with $L^2(\RR^d)$ initial data, we~have:
\[%\label{pfff}
\frac{d}{dt} \int \frac{(\rho_t^\eps)^2}{2} = - \int \nabla \rho_t^\eps \cdot (\nabla W * \rho_t^\eps) \rho_t^\eps - \eps \int |\nabla \rho_t^\eps|^2.
\]
First, using an adequate Young inequality on the first term along with the fact that $\nabla W$ is bounded allows to absorb the $|\nabla \rho_t^\eps|^2$ term into the last one, getting:
\[
\frac{d}{dt} \int \frac{(\rho_t^\eps)^2}{2} \leq
- \frac{\eps}{2} \int |\nabla \rho_t^\eps|^2 + \frac{a_\infty^2}{\eps} \int \frac{(\rho_t^\eps)^2}{2} \leq \frac{a_\infty^2}{\eps} \int \frac{(\rho_t^\eps)^2}{2}.
\]
A Grönwall lemma then ensures that $\int \sfrac{(\rho_t^\eps)^2}{2}$ remains bounded over finite time, where the bound depends on $\eps$, but $\eps$ is fixed. Second, plugging back this bound into the above estimate gives a bound on $\int_0^T |\nabla \rho_t^\eps|_{H^1(\RR^d)}^2 dt$. Let us reproduce these computations in the discrete setting.

\subsubsection*{Step 1: bound on $\sum_{J \in \ZZ^d} \sfrac{\big(\rho_J^n\big)^2}{2}$}

For the sake of compactness, let us note
\[
F_{J+ \sfrac{e_i}{2}}^n = ({a_i}^n_{J})^+ \rho_{J}^n - ({a_i}^n_{J+e_i})^- \rho_{J+e_i}^n.
\]
Using twice the definition of the explicit scheme \eqref{dis:explicit}, we~have:
\begin{multline*}
\sum_{J \in \ZZ^d} \frac{\big(\rho_J^{n+1}\big)^2 -(\rho_J^n)^2}{2} = \sum_{J \in \ZZ^d} \frac{\rho_J^{n+1}+ \rho_J^n}{2} \left(\rho_J^{n+1}-\rho_J^n\right) \\
= \sum_{J \in \ZZ^d}\! \frac{\rho_J^{n+1} + \rho_J^n}{2} \biggl(\!- \sum_{i=1}^{d} \frac{\dt}{\dx_i}
(F_{J+ \sfrac{e_i}{2}}^n - F_{J- \sfrac{e_i}{2}}^n)+ \eps \sum_{i=1}^{d} \frac{\dt}{\dx_i^2}(\rho_{J+e_i}^n - 2\rho_J^n + \rho_{J-e_i}^n)\biggr) \\
= - \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i}
\big(F_{J+ \sfrac{e_i}{2}}^n - F_{J- \sfrac{e_i}{2}}^n \big) \rho_J^n + \eps \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^n - 2\rho_J^n + \rho_{J-e_i}^n \right)\rho_J^n \\
+ \frac 12 \sum_{J \in \ZZ^d} \biggl(- \sum_{i=1}^{d} \frac{\dt}{\dx_i}
\big(F_{J+ \sfrac{e_i}{2}}^n - F_{J- \sfrac{e_i}{2}}^n \big) + \eps \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^n - 2\rho_J^n + \rho_{J-e_i}^n \right)\biggr)^2\\
:= S_1^n + S_2^n.
\end{multline*}
Performing discrete integrations by parts and using Young's inequality
\[
|ab| \leq \frac{a^2}{2\eps} + \frac{\eps b^2}{2}
\]
with $a=F_{J+ \sfrac{e_i}{2}}^n$ and $b=\psfrac{\rho_{J+e_i}^n - \rho_J^n}{\dx_i}$, we~can estimate $S_1^n$ as follows:
\begin{align*}
S_1^n & = \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i}
F_{J+ \sfrac{e_i}{2}}^n(\rho_{J+e_i}^n - \rho_J^n) - \eps \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} |\rho_{J+e_i}^n - \rho_J^n|^2 \\
& \leq \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \dt \Bigl(
\frac{\big(F_{J+ \sfrac{e_i}{2}}^n\big)^2}{2\eps} + \frac{\eps}{2}\,\Bigl|\frac{\rho_{J+e_i}^n - \rho_J^n}{\dx_i}\Bigr|^2 \Bigr) - \eps \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} |\rho_{J+e_i}^n - \rho_J^n|^2 \\
& \leq \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \dt
\frac{\big(F_{J+ \sfrac{e_i}{2}}^n\big)^2}{2\eps} - \frac{\eps}{2} \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} |\rho_{J+e_i}^n - \rho_J^n|^2.
\end{align*}
As for $S_2^n$, straightforward computations and the repeated use of $(a\pm b)^2 \leq 2a^2+2b^2$ lead to:
\[
S_2^n \leq \sum_{i=1}^{d} \frac{4 d \dt^2}{\dx_i^2} \sum_{J \in \ZZ^d} \big(F_{J+ \sfrac{e_i}{2}}^n\big)^2 + \sum_{i=1}^{d} 4 d \Bigl(\frac{\eps\dt}{\dx_i^2}\Bigr)^2 \sum_{J \in \ZZ^d} |\rho_{J+e_i}^n - \rho_J^n|^2.
\]
Using the fact that:
\[
\big(F_{J+ \sfrac{e_i}{2}}^n\big)^2 \leq \big(a_\infty \rho_J^n + a_\infty \rho_{J+e_i}^n\big)^2 \leq 2a_\infty^2 \Big(\big(\rho_J^n\big)^2 +(\rho_{J+e_i}^n)^2\Big),
\]
we deduce that $\sum_{J \in \ZZ^d} \big(F_{J+ \sfrac{e_i}{2}}^n\big)^2 \leq 4a_\infty^2 \sum_{J\in\ZZ^d}(\rho_J^n)^2$. Transferring in both estimates we found on $S_1^n$ and $S_2^n$, and summing both, we~obtain:
\begin{multline}\label{eq:disH1estimate:1}
\sum_{J \in \ZZ^d} \frac{(\rho_J^{n+1})^2 -(\rho_J^n)^2}{2}\leq \biggl(\frac{4d\dt a_\infty^2}{\eps} + \sum_{i=1}^{d} \frac{32 d a_\infty^2\dt^2}{\dx_i^2}\biggr) \sum_{J\in\ZZ^d} \frac{(\rho_J^n)^2}{2} \\
+ \sum_{i=1}^{d} \Bigl(-\frac{\eps\dt}{2\dx_i^2} + 4 d \Bigl(\frac{\eps\dt}{\dx_i^2}\Bigr)^2\Bigr) \sum_{J \in \ZZ^d} |\rho_{J+e_i}^n - \rho_J^n|^2.
\end{multline}

Under the Courant-Friedrichs-Lewy (CFL) condition
\begin{equation*}
\eps d \frac{\dt}{\dx_i^2} \leq \frac 18 \quad \text{for any } i,
\end{equation*}
the last term in the above estimate is nonpositive, thus we get
\begin{multline*}
\sum_{J \in \ZZ^d} \frac{(\rho_J^{n+1})^2 -(\rho_J^n)^2}{2} \leq \frac{4d\dt a_\infty^2}{\eps}\Bigl(1 + \sum_{i=1}^{d} \frac{8\eps\dt}{\dx_i^2}\Bigr) \sum_{J\in\ZZ^d} \frac{(\rho_J^n)^2}{2} \\
\leq \frac{4d\dt a_\infty^2}{\eps}\Bigl(1 + \frac{1}{d}\Bigr) \sum_{J\in\ZZ^d} \frac{(\rho_J^n)^2}{2} = \frac{4\dt a_\infty^2}{\eps}\left(1 + d\right) \sum_{J\in\ZZ^d} \frac{(\rho_J^n)^2}{2}.
\end{multline*}
Using a discrete Grönwall lemma, we~deduce the following bound on the discrete $L^2$ norm of $\rho_{\dx}^{\eps,n}$:
\[
\sum_{J \in \ZZ^d} \frac{(\rho_J^n)^2}{2} \leq \exp\Bigl(\frac{4(1 + d)t^n a_\infty^2}{\eps}\Bigr) \sum_{J \in \ZZ^d} \frac{(\rho_J^0)^2}{2}.
\]

\subsubsection*{Step 2: discrete $H^1$ bound}
Assume a stricter CFL condition: there exists $\delta$ such that
\begin{equation}\label{eq:CFL:16}
\eps d \frac{\dt}{\dx_i^2} \leq \delta < \frac 18 \quad \text{for any } i.
\end{equation}
Then, for any $i$,
\[
4d\Bigl(\frac{\eps \dt}{\dx_i^2}\Bigr)^2 - \frac{\eps \dt}{2\dx_i^2}
= \frac{\eps \dt}{2\dx_i^2}\Bigl(8d\frac{\eps \dt}{\dx_i^2} - 1\Bigr) \leq \frac{\delta}{d}(8\delta - 1) < 0.
\]
Thus, thanks to \ref{eq:disH1estimate:1},
\begin{multline*}
\sum_{i=1}^{d} \sum_{J \in \ZZ^d} |\rho_{J+e_i}^n - \rho_J^n|^2 \\
\leq
\frac{d}{\delta(1 - 8\delta)}\biggl(\Bigl(\frac{4d\dt a_\infty^2}{\eps} + \sum_{i=1}^{d} \frac{32 d a_\infty^2\dt^2}{\dx_i^2}\Bigr) \sum_{J\in\ZZ^d} \frac{(\rho_J^n)^2}{2}
- \sum_{J \in \ZZ^d} \frac{(\rho_J^{n+1})^2 -(\rho_J^n)^2}{2} \biggr) \\
\leq
\frac{d}{\delta(1 - 8\delta)}\biggl(\Bigl(\frac{4d\dt a_\infty^2}{\eps} + \sum_{i=1}^{d} \frac{4 a_\infty^2\dt}{\eps} \Bigr) \sum_{J\in\ZZ^d} \frac{(\rho_J^n)^2}{2}
- \sum_{J \in \ZZ^d} \frac{(\rho_J^{n+1})^2 -(\rho_J^n)^2}{2} \biggr),
\end{multline*}
which implies, thanks to the $L^2$ estimate,
\begin{multline*}
\sum_{i=1}^{d} \sum_{J \in \ZZ^d} |\rho_{J+e_i}^n - \rho_J^n|^2 \leq
\frac{d}{\delta(1 - 8\delta)} \biggl(\frac{8d\dt a_\infty^2}{\eps}
\exp\Bigl(\frac{4(1 + d)t^n a_\infty^2}{\eps}\Bigr) \sum_{J \in \ZZ^d} \frac{(\rho_J^0)^2}{2}
\\
- \sum_{J \in \ZZ^d} \frac{(\rho_J^{n+1})^2
-(\rho_J^n)^2}{2} \biggr).
\end{multline*}
Summing over $n=0,\ldots,m-1$ yields
\begin{multline*}
\sum_{n = 0}^{m-1} \sum_{i=1}^{d} \sum_{J \in \ZZ^d} |\rho_{J+e_i}^n - \rho_J^n|^2 \leq
\frac{d}{\delta(1 - 8\delta)} \biggl(\frac{8dT a_\infty^2}{\eps}
\exp\Bigl(\frac{4(1 + d)T a_\infty^2}{\eps}\Bigr) \sum_{J \in \ZZ^d} \frac{(\rho_J^0)^2}{2}
\\
 - \sum_{J \in \ZZ^d} \frac{(\rho_J^{m})^2}{2}
+ \sum_{J \in \ZZ^d} \frac{(\rho_J^0)^2}{2} \biggr).
\end{multline*}
Finally
\begin{multline*}
\sum_{n = 0}^{m-1} \sum_{i=1}^{d} \sum_{J \in \ZZ^d} |\rho_{J+e_i}^n - \rho_J^n|^2 \\
\leq
\frac{d}{\delta(1 - 8\delta)} \biggl(1 + \frac{8dT a_\infty^2}{\eps} \sum_{J\in\ZZ^d}
\exp\Bigl(\frac{4(1 + d)T a_\infty^2}{\eps} \Bigr)\biggr) \sum_{J \in \ZZ^d} \frac{(\rho_J^0)^2}{2}.
\end{multline*}
This is the desired result choosing $\delta = 1/16$.
\end{proof}
We now resume the proof of Theorem \ref{thm:cvOrder1/2:scheme}. From now on, we~always assume condition \eqref{eq:CFL:16} to hold.

\subsubsection*{Step 1: Ascoli-Arzelà theorem} Let us denote, for $K \subset \RR^d$ any compact set, $\Lip_K := \calC_c(K) \cap W^{1,\infty}(\RR^d)$ the space of Lipschitz continuous functions supported in $K$ and $\|\cdot\|_{\lip}$ the Lipschitz semi-norm. We~then introduce the pseudo-distance defined in duality with $\|{\cdot}\|_{\lip}$ by:
\[
\forall \mu, \nu \in \calP_1(\RR^d), \quad W_{1,K}(\mu, \nu) := \sup_{\substack{\phi \in \Lip_K\\ \|\phi\|_{\lip} \leq 1}} \int \phi d(\mu-\nu).
\]
For $0\leq s < t \leq T$, we~have, thanks to the Cauchy-Schwarz inequality:
\begin{equation}\label{ineq:HolderW1K}
W_{1,K}\big(\rho_{\dx,t}^\eps, \rho_{\dx,s}^\eps\big) = \int_s^t \big|(\rho_{\dx, \tau}^\eps)'\big| d\tau \leq \sqrt{t-s} \, \sqrt{\int_s^t \big|(\rho_{\dx,\tau}^\eps)'\big|^2 d\tau}.
\end{equation}
Here, the metric derivative is the one associated to the pseudo-distance $W_{1,K}$. Since we chose $\rho_{\dx}^\eps$ to be the piecewise affine reconstruction of the $\rho_{\dx}^{\eps,n}$ for $n=0, \ldots N$, we~have, for $\tau \in [t^n, t^{n+1}[$, $ |(\rho_{\dx,\tau}^\eps)'| = \frac{1}{\dt} W_{1,K}(\rho_{\dx}^{\eps,n}, \rho_{\dx}^{\eps,n+1})$. Indeed, $\rho_{\dx}^\eps$ is a constant-speed geodesic in $\WW_1(K)$ from $\rho_{\dx}^{\eps,n}$ to $\rho_{\dx}^{\eps,n+1}$ (recall \eqref{rhoaff} and the fact that linear interpolations are geodesic for the $\WW_1$ distance, which is a norm), hence its length on $[t^n, t^{n+1}[$ equals $\dt\big|(\rho_{\dx,\tau}^\eps)'\big|$ by definition and $W_{1,K}(\rho_{\dx}^{\eps,n}, \rho_{\dx}^{\eps,n+1})$ by the Benamou-Brenier formula. Therefore:
\begin{multline}\label{eq:pff}
\int_s^t \big|(\rho_{\dx,\tau}^\eps)'\big|^2 d\tau \leq \int_0^T \big|(\rho_{\dx,\tau}^\eps)'\big|^2 d\tau \\
= \sum_{k=0}^{N-1} \int_{t^n}^{t^{n+1}} \big|(\rho_{\dx,\tau}^\eps)'\big|^2 d\tau = \sum_{k=0}^{N-1} \frac{W_{1,K}^2(\rho_{\dx}^{\eps,n}, \rho_{\dx}^{\eps,n+1})}{\dt}.
\end{multline}
Now, let $\phi \in \Lip_K$ such that $\|\phi\|_{\lip} \leq 1$. We~have, denoting $\phi_J = \phi(x_J)$ and using the definition of the scheme \eqref{dis:explicit} along with discrete integrations by parts in space:
\begin{multline*}
\int \phi d\big(\rho_{\dx}^{\eps,n+1}-\rho_{\dx}^{\eps,n}\big) = \sum_{J\in\ZZ^d} \phi_J \left(\rho_J^{n+1} - \rho_J^{n}\right) \\
= \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i}
F_{J+ \sfrac{e_i}{2}}^{n}(\phi_{J+e_i} - \phi_J) - \eps \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \frac{\dt}{\dx_i^2} \left(\rho_{J+e_i}^{n} - \rho_J^{n}\right)(\phi_{J+e_i} - \phi_J) \\
\leq 2da_\infty \dt + \eps \dt \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \Bigl|\frac{\rho_{J+e_i}^{n} - \rho_J^{n}}{\dx_i}\Bigr|.
\end{multline*}
Taking the supremum over $\phi$ and using $(a + b)^2 \leq 2a^2 + 2b^2$, we~get:
\begin{align*}
W_{1,K}^2(\rho_{\dx}^{\eps,n}, \rho_{\dx}^{\eps,n+1}) & \leq 8d^2 a_\infty^2 \dt^2 + 2\eps^2 \dt^2 \bigg(\sum_{J \in \ZZ^d} \sum_{i=1}^{d} \Big|\frac{\rho_{J+e_i}^{n} - \rho_J^{n}}{\dx_i}\Big|\bigg)^2 \\
& \leq 8d^2 a_\infty^2 \dt^2 + 2 \eps^2 \dt^2 \frac{d \Leb(K)}{\prod_{i=1}^d \dx_i} \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \Big|\frac{\rho_{J+e_i}^{n} - \rho_J^{n}}{\dx_i}\Big|^2,
\end{align*}
where we used a discrete Cauchy-Schwarz inequality so as to use the discrete $H^1$ estimate we proved in Lemma \ref{lemma:DiscreteH1Estimates}: indeed, summing for $n=0,\ldots, N-1$ and plugging into \eqref{eq:pff}, we~obtain, using the aforementioned Lemma:
\begin{equation}\label{ineq:boundDisDerRho}
\begin{aligned}
\int_s^t \big|(\rho_{\dx,\tau}^\eps)'\big|^2 d\tau & \leq 8d^2 a_\infty^2 T + 2 d \eps^2 \frac{\Leb(K)}{\prod_{i=1}^d \dx_i} \dt \sum_{n=0}^{N-1} \sum_{J \in \ZZ^d} \sum_{i=1}^{d} \Big|\frac{\rho_{J+e_i}^{n} - \rho_J^{n}}{\dx_i}\Big|^2 \\
& \leq C(a_\infty,d, \eps, T, K) \Big(1 + \frac{1}{\prod_{i=1}^d \dx_i}\sum_{J\in\ZZ^d}(\rho_J^0)^2 \Big).
\end{aligned}
\end{equation}
Now, since we assumed that $\rho^{\ini} \in L^2(\RR^d)$, the term $\Psfrac{1}{\prod_{i=1}^d \dx_i}\sum_{J\in\ZZ^d}(\rho_J^0)^2$ is bounded with respect to $\dx$. Indeed, a Cauchy-Schwarz inequality along with our initialization of the scheme \eqref{dis:rho0} yield:
\begin{multline*}
\sum_{J\in\ZZ^d}(\rho_J^0)^2 = \sum_{J\in\ZZ^d} \bigg(\int_{C_J} \rho^{\ini}\bigg)^2 \\
\leq \sum_{J\in\ZZ^d} \Leb(C_J) \int_{C_J}(\rho^{\ini})^2 = \Bigl(\prod_{i=1}^d \dx_i\Bigr) \sum_{J\in\ZZ^d} \int_{C_J}(\rho^{\ini})^2 = \Bigl(\prod_{i=1}^d \dx_i\Bigr) \|\rho^{\ini}\|_{L^2}.
\end{multline*}
Transferring into \eqref{ineq:boundDisDerRho}, we~obtain a bound on $\int_s^t \big|(\rho_{\dx,\tau}^\eps)'\big|^2 d\tau$ that is uniform with respect to $s, t$ and $\dx$. Combining with \eqref{ineq:HolderW1K}, we~deduce that $(\rho_{\dx}^\eps)_{\dx > 0}$ is equi-$\frac 12$\nobreakdash-Hölder and in particular, equicontinuous in $\calC([0,T], (\Lip_K)')$. Lemma \ref{bounddismom} ensures, in~addition, that $M_2\big(\rho_{\dx,t}^\eps\big)$ is uniformly bounded with respect to $t \in [0,T]$ and $\dx>0$. Using Lemma \ref{lemma:boundedWpImpliesCompactWq}, we~deduce that $(\rho_{\dx,t}^\eps)_{\dx > 0}$ lies in a relatively compact set for all $t \in [0,T]$ and $\dx>0$. We~can therefore apply the Ascoli-Arzelà theorem along with a diagonal extraction to extract a subsequence, that we still denote $(\rho_{\dx}^\eps)_{\dx > 0}$, converging in
$\calC([0,T], \WW_1(\RR^d))$.

\subsubsection*{Step 2: $\rho^\eps$ solves \eqref{eq:agg_diff}}

Using discrete integrations by parts as in \cite{cjlv15, lv17}, we~can prove that $\rho_{\dx}^\eps$ satisfies the following approximate weak form of \eqref{eq:agg_diff}, for any $\phi \in \calC([0,T[ \times \RR^d)$:
\begin{multline}\label{eq:weakFormulationApprox}
\int_0^T \int \pa_t \phi(t,x) \rho_{\dx,t}^\eps(dx) dt\\
+ \int_0^t \int \achapo[\rho_{\dx,t}^\eps] \cdot \nabla \phi(t,x) \rho_{\dx,t}^\eps(dx) dt + \int \phi(0,x) \rho^{\ini}(dx) \\
= \eps \int_0^T \int \Delta \phi(t,x) \rho_{\dx,t}^\eps(dx) + O(\dx) + O(\dt).
\end{multline}
Passing to the limit $\dx \to 0$ in \eqref{eq:weakFormulationApprox} is straightforward for the linear terms since $\rho_{\dx,t}^\eps \wslimi{\dx}{0} \rho_t^\eps$ uniformly in time. For the nonlinear term, this convergence also ensures that $\rho_{\dx,t}^\eps \otimes \rho_{\dx,t}^\eps \wslimi{\dx}{0} \rho_t^\eps \otimes \rho_t^\eps$. Then, passing to the limit is done using a symmetrization argument as in equations \eqref{eq:passlimNLterm:1}--\eqref{eq:passlimNLterm:2}--\eqref{eq:passlimNLterm:3} using the fact that $W$ is Lipschitz and even.

We deduce that $\rho^\eps$ solves in the sense of distributions the aggregation-diffusion problem \eqref{eq:agg_diff} with initial datum $\rho_0^\eps = \rho^{\ini}$. Since such a solution is unique (see Theorem \ref{thm:agg_diff:WP}), we~deduce that actually the whole initial sequence $(\rho_{\dx}^\eps)_{\dx > 0}$ converges towards $\rho^\eps$.

\Subsubsection*{Step 3: passing to the limit in \eqref{eq:cvOrder1/2:scheme:bound:1} and relaxing the assumption $\rho^{\ini} \in L^2(\RR^d)$}

Now, let $t>0$ and let $n \in \{0,\ldots,N\}$ such that $t \in [t^n, t^{n+1}[$. Estimate \eqref{eq:cvOrder1/2:scheme:bound:1} gives:
\[
W_2(\rho_{t},\rho_{\dx,t}^{\eps}) \leq C \sqrt{\frac{1 - e^{-4\lambda t}}{\lambda}}\, \sqrt{\dx + \eps} + e^{-2\lambda t} \dx.
\]
Passing to the limit $\dx \to 0$ in the above estimate using the semicontinuity of $W_2$ then gives the desired estimate \eqref{eq:cvOrder1/2:scheme:bound:2}, hence proving Theorem \ref{thm:cvOrder1/2:scheme} in case of $L^2(\RR^d)$ initial datum. The general case can be obtained by approximation, using Assumption~\eqref{A3} which guarantees stability of the solutions for both $\eps = 0$ and $\eps > 0$. This ends the proof of Theorem \ref{thm:cvOrder1/2:scheme}.

\begin{remark}
As a byproduct of this proof, we~obtain uniform in time convergence in the $W_1$ distance in space of the numerical scheme \eqref{dis:explicit}--\eqref{dis:rho0} towards the $\calC([0,T], \WW_2(\RR^d))$ distributional solution to the aggregation-diffusion initial value problem, in case of $L^2(\RR^d)$ initial datum, and under a $1/6$-CFL condition. In fact, we~expect this convergence result to hold for arbitrary $\calP_2(\RR^d)$ initial datum and under the standard CFL condition:
\[
\sum_{i=1}^d \Bigl(a_\infty \frac{\dt}{\dx_i} + 2\eps \frac{\dt}{\dx_i^2}\Bigr) \leq \frac 16.
\]
\end{remark}

\section{Convergence for repulsive potentials such that
\texorpdfstring{$\Delta W \leq 0$\nobreakspace and\nobreakspace $\nabla^2 W \in L^{p_0}(\RR^d)$}{titlesection2}}\label{section:repulsive}

For any Lipschitz potential satisfying assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}, Theorem \ref{thm:cvLipschitz} guarantees the convergence of $\rho^\eps$ towards a solution $\rho$ to the aggregation equation up to a subsequence if the initial data satisfies the assumptions \eqref{assumptionInitialData}. Then, Corollary \ref{coroll:cvLambdaCvx} extended this result to arbitrary initial data by an approximation procedure, and using $\lambda$-convexity to estimate the distance between two solutions. The goal of this section is to proceed similarly in the case of repulsive potentials, typically $W(x) = -|x|$, where $\lambda$-convexity will be replaced by some integrability of the Hessian. More precisely, we~focus on initial data equal to $\rho^{\ini}$, for which we only assume finiteness of moments.

The outline of the proof is the same as that of Corollary \ref{coroll:cvLambdaCvx}. However, we~can no more use the $\lambda$-convexity of $W$ but, using the additional assumption $\nabla^2 W \in L^{p_0}(\RR^d)$ for a suitable $p_0$, we~still manage to estimate the distance between $\rho_t^\eps$ and a sequence of viscous solutions associated with smoothed out initial data. More precisely, we~obtain the following result:
\begin{theorem}\label{thm:cvRepulsive}
Let $W$ be an interaction potential satisfying assumptions \eqref{A0}--\eqref{A1}--\eqref{A2} along with the additional assumption:
\begin{enumerate}\renewcommand{\theenumi}{\textup{A}\arabic{enumi}}\setcounter{enumi}{4}
\item
\label{A5} $\Delta W \leq 0$ and $\nabla^2 W \in L^{p_0}(\RR^d)$ for some $p_0 > \max(\sfrac d2, 1)$,
\end{enumerate}
and let $\rho^{\ini}$ be an initial datum belonging to $ \in \calP_2(\RR^d)$.
Denote $(\rho^\eps)_{\eps > 0}$ the sequence of weak solutions to \eqref{eq:agg_diff} where the initial data is set to $\rho_0^\eps := \rho^{\ini}$ for all $\eps > 0$.

Then, for all $T>0$, the sequence $(\rho^\eps)_{\eps > 0}$ converges in $\calC([0,T], \WW_1(\RR^d))$, up to an extraction, towards a solution $\rho \in \calC([0,T], \WW _2(\RR^d))$ to equation \eqref{eq:agg} with the velocity field $a[\rho]$ being replaced by $\achapo[\rho]$ as defined in \eqref{eq:achapo}.

If, in addition, $\rho^{\ini} \in L^{p_0'}(\RR^d) \cap L^{\spfrac{p_0}{p_0-p}}(\RR^d)$, then there exists a unique solution in $\calC([0,T], \WW_2(\RR^d)) \cap L^\infty([0,T], L^{p_0'}(\RR^d) \cap L^{\spfrac{p_0}{p_0-p}}(\RR^d))$ to \eqref{eq:agg} and actually the whole sequence $(\rho^\eps)_{\eps > 0}$ converges.
\end{theorem}

\skpt
\begin{remark}
\begin{enumerate}
\item For $W(x) = -|x|$, this result cannot be applied in dimension $d=1$, since $\nabla^2 W = - \delta_0$ is not integrable. When $d>1$, we~have
\[
\nabla^2 W(x) = \frac{\Psfrac{x}{|x|} \otimes \Psfrac{x}{|x|} - I_d}{|x|} \sim \frac{1}{|x|},
\]
hence $\nabla^2 W \in L^{p_0}$ if and only if $p_0 < d$ (up to cutting off the potential at infinity) and therefore we can find $p_0 \in(\sfrac d2, d)$ so as to apply our result.
\item In dimension $d=1$, for $W(x) = - |x|$, Proposition \ref{prop:cvViaBurgers} shows that the whole sequence $(\rho^\eps)_{\eps > 0}$ converges in $\calC([0,T], \WW_1(\RR))$ towards a solution to the aggregation equation that can be obtained as the derivative of the entropy solution to a Burgers-type equation since entropy solutions and viscosity solutions coincide for scalar conservation laws.
\item As a byproduct of our result, one obtains existence of a solution in
\[
\calC([0,T], \WW _2(\RR^d))
\]
to the aggregation problem \eqref{eq:agg} for potentials satisfying \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A5}.
\end{enumerate}
\end{remark}
\begin{proof}
Let $T>0$. As in the proof of Corollary \ref{coroll:cvLambdaCvx}, for $\eps > 0$, we~introduce $\mu^\eps \in \calC([0,T], \WW_2(\RR^d))$ solution to \eqref{eq:agg_diff} with smoothed out initial data $\mu_0^\eps$, that we now assume to satisfy assumptions \eqref{smoothedInitialData} for some $\alpha \in (-1,0)$. In particular, $(\mu_0^\eps)_{\eps > 0}$ satisfies assumptions \eqref{assumptionInitialData} and Theorem \ref{thm:cvLipschitz} applies to $(\mu^\eps)_{\eps>0}$ and guarantees convergence of a subsequence, in $\calC([0,T], \WW_1(\RR^d))$, towards a solution to the aggregation equation \eqref{eq:agg}. As for Corollary \ref{coroll:cvLambdaCvx}, the key ingredient is now to prove that the distance $W_p(\rho_t^\eps, \mu_t^\eps)$ goes to 0 as $\eps \to 0$, for some $p > 1$ that will be specified later.

For the sake of clarity, let us drop the superscripts $\eps$ for the remaining of this section.

Denoting $(\varphi_t, \psi_t)$ a pair of Kantorovitch potentials from $\rho_t$ to $\mu_t$ for the cost $\frac{1}{p}|x-y|^p$, we~can formally write (see \cite[Th.\,5.24]{otam15} or \cite[Th.\,8.4.7]{ambrosio08})
\begin{multline*}
\frac 1p \frac{d}{dt} W_p^p(\rho_t, \mu_t) = \int \nabla \varphi_t \cdot a[\rho_t] d\rho_t
+ \int \nabla \psi_t \cdot a[\mu_t] d\mu_t \\
- \eps \int \Big(\nabla \varphi_t \cdot \nabla \rho_t + \nabla \psi_t \cdot \nabla \mu_t \Big) dx.
\end{multline*}
The last term above is nonnegative thanks to the so-called five (actually four) gradients inequality proved in \cite{caillet21} for the $W_p$ case with $p>1$. Actually, \cite{caillet21} proves the inequality in a compact setting and a full treatment of this last term would require a suitable approximation procedure. Yet, the inequality we need, \ie
\[
\frac 1p \frac{d}{dt} W_p^p(\rho_t, \mu_t) \leq \int \nabla \varphi_t \cdot a[\rho_t] d\rho_t + \int \nabla \psi_t \cdot a[\mu_t] d\mu_t
\]
can also be justified in many different ways, for instance by the stochastic interpretation of $\rho_t$ and $\mu_t$ as laws of the solutions of suitable SDE where the choice of a common Brownian motion would allow to get rid of the term coming from diffusion (see, for instance, \cite{BolGenGui}); since the diffusion effect of the Laplacian in the equation could also be handled using convolution with the heat kernel, another possible way to prove the same inequality would be to approximate the solutions by a splitting method, alternating convolutions (which decrease the $W_p$ distance) and transport (which lets the other term appear).

We thus get, using a triangle inequality along with the fact that
\[
\nabla \varphi_t(x) = |x - T_t(x)|^{p-1} (\widehat{x - T_t(x)}) = - \nabla \psi_t(x),
\]
where $T_t$ is the optimal transport map from $\rho_t$ to $\mu_t$ (which exists since $\rho_t \ll \Leb$ whenever $\eps>0$):
\begin{subequations}\label{eq:I12}
\begin{align}
	& \frac 1p \frac{d}{dt} W_p^p(\rho_t, \mu_t) \leq |I_1| + |I_2|, \\
& I_1 = \int |x - T_t(x)|^{p-1} (\widehat{x - T_t(x)}) \cdot (a[\rho_t](x) - a[\rho_t] \circ T_t(x)) \rho_t(dx), \label{eq:I1} \\
& I_2 = \int |x - T_t(x)|^{p-1} (\widehat{x - T_t(x)}) \cdot (a[\rho_t] \circ T_t(x) - a[\mu_t] \circ T_t(x)) \rho_t(dx). \label{eq:I2}
\end{align}
\end{subequations}
To estimate $I_1$, we~use the following bound on the Lipschitz constant of $a[\rho_t]$:
\[
\lip (a[\rho_t]) = \|\nabla^2 W * \rho_t\|_{L^\infty} \leq \|\nabla^2 W\|_{L^{p_0}} \|\rho_t\|_{L^{p_0'}}.
\]
We deduce:
\[
|I_1| \leq \lip (a[\rho_t]) \int |x - T_t(x)|^p \rho_t(dx) \leq \|\nabla^2 W\|_{L^{p_0}} \|\rho_t\|_{L^{p_0'}} W_p^p(\rho_t, \mu_t).
\]
To estimate $I_2$, we~first apply a Hölder inequality with respect to the measure $\rho_t(dx)$ and with the exponents $(p',p)$. We~get, since $p'(p-1)=p$:
\begin{equation}\label{ineq:I2}
|I_2| \leq \bigg( \int |x - T_t(x)|^p \rho_t(dx) \bigg)^{1/p'} \bigg( \int \big|a[\rho_t] \circ T_t(x) - a[\mu_t] \circ T_t(x) \big|^p \rho_t(dx)\bigg)^{1/p}.
\end{equation}
We recognize that the first factor equals $W_p^{p-1}(\rho_t, \mu_t)$ since $\sfrac{p}{p'} = p-1$.
Let us deal with the second one. We~consider $\nu_s := \left((1-s)\text{id} + sT_t\right)_\#\rho_t$ the constant-speed geodesic from $\rho_t$ to $\mu_t$. Note that this curve implicitly depends on $t$. We~also denote by $b_s \in L^p(\nu_s)$ the velocity field associated with $\nu \in AC([0,1], \WW_p(\RR^d))$, as given by Theorem~\ref{thm:ACcurves}. We~have as a consequence of the Benamou-Brenier formula $\pa_s \nu_s + \nabla \cdot(b_s\nu_s) = 0$ and $\|b_s\|_{L^p(\nu_s)} = |(\nu_s)'| = W_p(\rho_t, \mu_t)$ for a.e. $s \in [0,1]$. Therefore, for any $y \in \RR^d$, one has:
\begin{align*}
a[\rho_t](y) - a[\mu_t](y) & = - \int \nabla W(y-z) (\rho_t(z) - \mu_t(z)) dz \\
& = - \int_0^1 \int \nabla W(y-z) \pa_s \nu_s(z) dz ds \\
& = \int_0^1 \int \nabla W(y-z) \nabla \cdot (b_s(z) \nu_s(z)) dz ds \\
& = \int_0^1 \int \nabla^2 W(y-z) b_s(z) \nu_s(dz) ds,
\end{align*}
so that the inequality \eqref{ineq:I2} rewrites:
\[
|I_2| \leq W_p^{p-1}(\rho_t, \mu_t) \bigg( \int \bigg| \int_0^1 ds \int \nabla^2 W(T_t(x) - z) b_s(z) \nu_s(dz) ds \bigg|^p \rho_t(dx)\bigg)^{1/p}.
\]
Besides, using a Jensen inequality with respect to the measure $\nu_s(dz) ds$ for the convex function $|\cdot|^p$, we~have:
\begin{multline*}
\int \bigg| \int_0^1 ds \int \nabla^2 W(T_t(x) - z) b_s(z) \nu_s(dz) ds \bigg|^p \rho_t(dx) \leq \\
\int \int_0^1 \int |\nabla^2 W(T_t(x) - z)|^p |b_s(z)|^p \nu_s(dz) ds \rho_t(dx) \\
\leq \int_0^1 \int |b_s(z)|^p \int |\nabla^2 W(T_t(x) - z)|^p \rho_t(dx) \nu_s(dz) ds.
\end{multline*}
Now, since $\mu_t = {T_t}_\# \rho_t$, we~have $\int |\nabla^2 W(T_t(x) - z)|^p \rho_t(dx) = \int |\nabla^2 W(y - z)|^p \mu_t(y) dy$. Applying a Hölder inequality with respect to $dy$ and the exponents $(q,q')$, where we will specify $q$ right afterward, we~get:
\begin{multline*}
\int |\nabla^2 W(y - z)|^p d\mu(y) \leq \bigg( \int |\nabla^2 W(y - z)|^{pq'} dy \bigg)^{1/q'} \bigg( \int |\mu_t(y)|^{q} dy \bigg)^{1/q} \\
= \|\nabla^2 W\|_{L^{pq'}}^p \|\mu_t\|_{L^q}.
\end{multline*}
We therefore have to take $q$ such that $pq' = p_0$, so that $\|\nabla^2 W\|_{L^{pq'}}$ remains finite. This requires that we choose $p$ such that $p \leq p_0$, which imposes $p_0 > 1$ since we also needed $p>1$. We~also need to choose $p$ such that $\rho^{\ini}\in\calP_p$, which means $p\leq 2$. Using $\int_0^1 \int |b_s(z)|^p \nu_s(dz) ds = W_p^p(\rho_t, \mu_t)$, we~finally obtain:
\[
|I_2| \leq \|\nabla^2 W\|_{L^{p_0}} \|\mu_t\|_{L^q}^{1/p} W_p^p(\rho_t, \mu_t),\qquad\mbox{ for } q=\frac{p_0}{p_0-p},
\]
where the value of $q$ is computed so that we have $q'=\sfrac{p_0}{p}$.
We therefore have the following Grönwall inequality on $W_p^p(\rho_t, \mu_t)$:
\begin{equation}\label{ineq:GWp}
\frac 1p \frac{d}{dt} W_p^p(\rho_t, \mu_t) \leq \|\nabla^2 W\|_{L^{p_0}} \left( \|\rho_t\|_{L^{p_0'}} + \|\mu_t\|_{L^q}^{1/p} \right) W_p^p(\rho_t, \mu_t).
\end{equation}
Now, we~need a bound on $\|\rho_t\|_{L^r}$. The following lemma implies that, if the interaction potential $W$ satisfies $\Delta W \leq 0$, then the bound on $\rho_t$ is not worse than the one we would obtain if $\rho$ solved the sole heat equation and does not depend on the initial datum.
\begin{lemma}\label{lemma:LpBoundRho}
Let $p \in (1,+\infty)$, $\eps > 0$ and let $\rho$ solve the following Fokker-Planck equation on the whole space $\RR^d$:
\begin{equation}\label{eq:FokkerPlanck}
\pa_t \rho + \nabla \cdot (\rho \nabla V) = \eps \Delta \rho,
\end{equation}
where the potential $V$ might depend on $\rho$ and satisfies $\Delta V \geq 0$. Assume that $\rho_t$ is smooth for any $t>0$, and that is has unit total mass. Then one has:
\[
\|\rho_t\|_{L^p} \leq C(\eps t)^{-d/2p'},
\]
for a positive constant $C=C(p,d)$ depending on $p$ only and not on the initial datum~$\rho_0$.
\end{lemma}
\begin{proof}
In the following, $C(p)$ stands for any positive constant depending only on $p$. For $t>0$, testing equation \eqref{eq:FokkerPlanck} against $\rho_t^{p-1}$ and integrating by parts yields:
\[
\frac{d}{dt} \frac 1p \int \rho_t^p = - \frac{p-1}{p} \int \rho_t^p \Delta V - 4 \eps \frac{p-1}{p^2} \int |\nabla \rho_t^{p/2}|^2 \leq - 4 \eps \frac{p-1}{p^2} \int |\nabla \rho_t^{p/2}|^2,
\]
since $\Delta V \geq 0$. Using the following Gagliardo-Nirenberg-Sobolev inequality \cite{gagliardo59, nirenberg59}:
\[
\int \rho^{p+\sfrac 2d} \leq C(p) \int |\nabla \rho_t^{p/2}|^2,
\]
and interpolating the $L^p$ norm between the $L^1$ and $L^{p+\frac 2d}$ norms, we~deduce that $y_t := \int \rho_t^p$ verifies the following nonlinear Grönwall inequality:
\[
y' - \eps C(p) y^{1+\sfrac{2}{d(p-1)}} \leq 0.
\]
Integrating this inequality on $[s,t]$ for $0<s<t$, we~get:
\[
y_t^{-2/d(p-1)} \geq y_s^{-2/d(p-1)} + \eps C(p) \geq \eps C(p),
\]
and therefore $\|\rho_t\|_{L^p} = y_t^{1/p} \leq C(p) (\eps t)^{-d(p-1)/2} = (\eps t)^{-d/2p'}$. This is the bound one would obtain using a $L^p \times L^1$ convolution inequality if $\rho$ solved the sole heat equation on the whole space, that is, if we had $\rho_t = G_{\eps t} * \rho_0$ where $G_t$ denotes the heat kernel.
\end{proof}

Using Lemma \ref{lemma:LpBoundRho} with the potential $V = -W*\rho$ which has a positive Laplacian under the assumption $\Delta W \leq 0$, we~get $\|\rho_t\|_{L^{p_0'}} + \|\mu_t\|_{L^q}^{1/p} \leq C(d,p_0)(\eps t)^{-d/2p_0}$ which, in turn, yields the Grönwall inequality:
\[
\frac{d}{dt} W_p^p(\rho_t, \mu_t) \leq C(\eps t)^{-d/2p_0} W_p^p(\rho_t, \mu_t),
\]
where $C$ is a positive constant that depends on $p$, $p_0$ and $\|\nabla^2 W\|_{L^{p_0}}$ only. We~deduce:
\begin{equation*}%\label{eq:contractionRepulsif}
W_p^p(\rho_t, \mu_t) \leq W_p^p(\rho_0, \mu_0) e^{\int_0^t C(\eps \tau)^{-d/2p_0} d\tau},
\end{equation*}
provided $p_0 > \sfrac d2$ so that $\tau^{-d/2p_0}$ is integrable on $(0,t]$. Under this assumption, using Lemma \ref{lemma:initialData} along with the fact that $\rho_0 = \rho^{\ini}$, we~get, for some constant $C>0$ depending on $d, p, p_0$ and $\|\nabla^2 W\|_{L^{p_0}}$ only:
\[
\forall t \in [0,T], \quad W_p^p(\rho_t, \mu_t) \leq C e^{-C(\eps^{\alpha} - \eps^{-d/2p_0})} e^{Ct^{1-d/2p_0}} \leq C e^{-C(\eps^{\alpha} + \eps^{-d/2p_0})} e^{CT^{1-d/2p_0}},
\]
which goes to 0 uniformly in $t \in [0,T]$, as $\eps \to 0$, provided $\alpha < -d/2p_0$. Since $-d/2p_0 > -1$, it is possible make such a choice while guaranteeing $\alpha \in (-1,0)$. To~finish, we~conclude the proof as in that of Corollary \ref{coroll:cvLambdaCvx}.

Now, note that $\Delta W \leq 0$ ensures that any $L^p$ norm of solutions to \eqref{eq:agg} is nonincreasing in time. Therefore, when the initial datum belong to $L^{p_0'}(\RR^d) \cap L^{\spfrac{p_0}{p_0-p}}(\RR^d)$, estimate \eqref{ineq:GWp} still holds for $\eps=0$ between any two solutions to \eqref{eq:agg} and gives uniqueness of the solution among the class of
\[
\calC\bigl([0,T], \WW_2(\RR^d)\bigr) \cap L^\infty\bigl([0,T], L^{p_0'}(\RR^d) \cap L^{\spfrac{p_0}{p_0-p}}(\RR^d)\bigl)
\]
solutions.
\end{proof}

\section{Higher convergence rate for steady states under assumptions\protect\line \texorpdfstring{\protect\eqref{A0}--\protect\eqref{A1}--\protect\eqref{A4p}}{Ap}}\label{section:stationary}

In this section, we~compare stationary solutions to the aggregation-diffusion equation \eqref{eq:rho_diff} for a given $\eps>0$ with stationary solutions to the aggregation equation~\eqref{eq:agg}. We~discard, in this section, the assumptions of $\lambda$-convexity and Lipschitz continuity on $W$ but still assume that assumptions \eqref{A0} and \eqref{A1} hold. In addition, we~require the potential to satisfy assumption \eqref{A4p}, that is, to be at least as attractive as $|x|^p$, for some $p \in [1, \infty)$.
These stationary solutions are in many cases long-time limits of the corresponding evolving solutions, but we will not insist on these aspects that are usually studied by $\lambda$-convexity techniques, and we discarded such an assumption in this section.

Note that this assumption along with \eqref{A0} implies $W(x) \geq C \sfrac{|x|^p}{p}$ for all $x \in \RR^d$. If, in addition, $W$ satisfies assumption \eqref{A1} then $W$ is l.s.c.\ on $\RR^d$ and this implies that $\calW$ is l.s.c.\ for the weak convergence thanks to Lemma \ref{lemma:W_lsc}.

Also, without loss of generality, we~only consider measures with 0 center of mass, that is, measures $\rho \in \calP(\RR^d)$ verifying:
\[
\int x\rho (dx) = 0.
\]
We define steady states for the aggregation-diffusion equation in the spirit of \cite{kang19}:
\begin{defi}
Let $\eps \geq 0$. A steady state for the aggregation-diffusion equation \eqref{eq:rho_diff} is a probability measure $\rho \in \calP_1(\RR^d)$ such that:
\[
\text{if } \eps = 0, \quad \nabWchapo * \rho = 0, \quad \text{on } \mathrm{supp}(\rho),
\]
and, if $\eps > 0$:
\[
\begin{cases}
\nabla W * \rho + \eps \dfrac{\nabla \rho}{\rho} = 0 &\text{on } \RR^d,\\
\rho > 0&\text{on } \RR^d.
\end{cases}
\]
\end{defi}
One can prove that this definition is equivalent to that of stationary solutions, in~the sense of distributions, to equation \eqref{eq:agg_diff}. Besides, if $\eps > 0$, one can show that a distributional solution to the elliptic problem $- \nabla \cdot (\nabla W * \rho) \rho = \eps \Delta \rho$ is necessarily regular and positive on $\RR^d$ (see Theorem \ref{thm:agg_diff:WP}).

The following lemma justifies why we compare steady states for the aggregation equation to the Dirac mass.
\begin{lemma}
Under assumptions \eqref{A0}--\eqref{A1}--\eqref{A4p} for $p \geq 1$, the unique steady state for the aggregation equation \eqref{eq:rho} is, up to a translation, the Dirac mass $\delta_0$.
\end{lemma}

\begin{proof}
Let $\rho$ be a steady state for \eqref{eq:agg} and assume that $\rho$ is centered. Since $\nabWchapo * \rho = 0$ on the support of $\rho$, testing against $\rho x$ and using Lemma \ref{lemma:symmetrization} with the odd vector field $\nabWchapo$ yields:
\[
\iint \nabWchapo(x-y) \cdot (x-y) \rho(dx) \rho(dy) = 0.
\]
Under assumption \eqref{A4p}, we~therefore have $\iint |x-y|^p \rho(dx) \rho(dy) = 0$. In particular $\rho\otimes\rho$ is concentrated on the diagonal. Now, if $\rho$ is not a Dirac mass, then there exists disjoint Borel sets $A$ and $B$ with $\rho(A)>0$ and $\rho(B)>0$. Then we have, since $A\times B$ is disjoint from the diagonal
\[
0 =\rho \otimes \rho(A \times B) =\rho(A)\rho(B)> 0,
\]
and this contradiction concludes the proof.
\end{proof}

Note that the Dirac mass is actually the only minimizer of the interaction energy $\calW$ under these assumptions. Conversely, Proposition 7.20 in \cite{otam15} ensures that minimizers of the energy $F^\eps$ are actually steady states. This provides a way to prove existence of steady states for \eqref{eq:rho_diff} when $\eps >0$.

\Subsection{Existence of minimizers of $F^\eps$ for $\eps >0$}

\begin{prop}
Assume that $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A4p} for some $p \geq 1$ and let $\eps \geq 0$ be fixed. The functional $F^\eps = \calW + \eps \calU$ admits a minimizer over $\calP(\RR^d)$ that actually has finite $p$-th order moment.
\end{prop}

\begin{remark}
We were not able to prove uniqueness of the minimizer under such assumptions on $W$ but it is likely to hold. Moreover, numerical illustrations will show that, if we remove assumption \eqref{A4p}, multiple steady states can coexist even though $\eps>0$ (in case $\eps=0$, it is easy to build explicit counterexamples).
\end{remark}
To prove this proposition, we~will use that under assumptions \eqref{A0} and \eqref{A4p}, controlling $\calW(\rho)$ gives control on $\iint |x-y|^p \rho(dx) \rho(dy)$, and this latter quantity is equivalent to $M_p(\rho)$ whenever $\rho$ is centered, thanks to the following lemma:
\begin{lemma}\label{lemma:equivalence_Wp}
Let $p \in [1, \infty)$ and $\rho \in \calP_p(\RR^d)$. Assume that the center of mass of $\rho$ is $0$. Then:
\[
M_p(\rho) \leq \iint |x-y|^p \rho(dx) \rho(dy) \leq 2^{p-1} M_p(\rho).
\]
\end{lemma}

\begin{proof}
Let $u(x) = \int |x-y|^p \rho(dy)$. Since $p \geq 1$, $u$ is a convex function and therefore, using a Jensen inequality, we~get:
\[
M_p(\rho)=u(0) = u \Big(\int x\rho(dx)\Big) \leq \int u(x) \rho(dx).
\]
In other terms, $M_p(\rho) \leq \iint |x-y|^p \rho(dx) \rho(dy)$. The upper bound comes from the inequality $|x-y|^p \leq 2^{p-1} (|x|^p + |y|^p)$.
\end{proof}

\begin{proof}
Let $(\rho_n)_{n \in \NN}$ be a sequence of probability measures that minimize $F^\eps$. We~can assume that these measures are centered because $F^\eps$ is invariant under translation. Up to an extraction, we~can assume that $(\rho_n)_{n \in \NN}$ converges weakly towards some $\rho \in \calM_b(\RR^d)$. To ensure that $\rho \in \calP(\RR^d)$, we~need to prove tightness of $(\rho_n)_{n \in \NN}$. To do so, let us find a bound on $M_p(\rho_n)$.

Since $(\rho_n)_{n \in \NN}$ is a minimizing sequence, $F^\eps(\rho_n) = \calW(\rho_n) + \eps \calU(\rho_n)$ is bounded from above by some constant that we still denote $C>0$. Moreover, using assumption \eqref{A0} and \eqref{A4p} and Lemma \ref{lemma:equivalence_Wp}, since $\rho_n$ is centered, we~have:
\[
\calW(\rho_n) \geq \frac{C}{2p} \iint |x-y|^p \rho_n(dx) \rho_n(dy) \geq \frac{C}{2p}\, M_p(\rho_n).
\]
In order to get a lower bound involving $M_p(\rho_n)$ on the entropy term, recall that, using a Legendre transform, $y\ln y + e^{z-1} \geq yz$ for all $y \geq 0$ and $z \in \RR$. Setting, for $x \in \RR^d$, $y = \rho_n(x)$ and $z = -|x|^{\alpha p}$ for some exponent $\alpha > 0$ to be specified later, and integrating over $x \in \RR^d$, we~get:
\[
\int \rho_n \ln \rho_n \geq - \int (|x|^p)^\alpha \rho_n(dx) + \int e^{-|x|^{\alpha p} -1} dx.
\]
Choosing $\alpha \in (0,1)$ so that $x \mto |x|^\alpha$ is concave, and using a Jensen inequality, we~deduce $\calU(\rho_n) \geq - M_p(\rho_n)^\alpha + C(p,\alpha)$, where $C(p,\alpha)$ depends on $\alpha$ and $p$ only. Finally, we~obtain:
\[
\frac{C}{2p}\, M_p(\rho_n) - \eps M_p(\rho_n)^\alpha + \eps C(p,\alpha) \leq C,
\]
which implies, since $\alpha < 1$, that $M_p(\rho_n)$ is uniformly bounded with respect to $n$.

On the one hand, this implies that $(\rho_n)_{n \in \NN}$ is tight, hence $\rho \in \calP(\RR^d)$. Since $M_p$ is l.s.c.\ on $\calP(\RR^d)$ and $\rho_n \wslimi{n}{+\infty} \rho$, we~also get $\rho \in\calP_p(\RR^d)$. On the other hand, the uniform bound on $M_p(\rho_n)$ along with Lemma \ref{lemma:boundedWpImpliesCompactWq} ensures that $\rho_n$ is compact in $W_q$ and hence we obtain $M_q(\rho_n) \limi{n}{+\infty} M_q(\rho)$ for any $q \in (0,p)$. Lemma \ref{lemma:U_lsc} then gives $\calU(\rho) \leq \liminf_{n \to +\infty} \calU(\rho_n)$, and, since $\calW$ is l.s.c.\ for the weak convergence, we~get $F^\eps(\rho) \leq \liminf_{n \to +\infty} F^\eps(\rho_n)$. This proves that $\rho$ minimizes $F^\eps$ since $(\rho_n)_{n \in \NN}$ is a minimizing sequence.
\end{proof}

\Subsection{$O(\eps)$ convergence rate in $W_p$ for potentials such that $\nabla W(x) \cdot x \geq C |x|$}

In this section, we~focus on assumption (A4-1) under which the potential is ``really pointy'' and the aggregation compensates the diffusion so that convergence occurs at rate $O(\eps)$:
\begin{theorem}\label{thm:W1_leq_Ceps}
Assume that $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\textup{(A4-1)}. There exists a constant $C > 0$ depending on $d$, such that for any $\eps>0$ and $\rho^\eps$ steady state for \eqref{eq:rho_diff} which center of mass is $0$, the following estimate holds:
\begin{equation}\label{W1_leq_Ceps}
W_1(\rho^\eps, \delta_0) \leq C\eps.
\end{equation}
\end{theorem}
\begin{proof}[Proof of Theorem \ref{thm:W1_leq_Ceps}]
Let $\eps>0$ and let $\rho^\eps$ be a steady state for \eqref{eq:agg_diff}, that is:
\begin{equation}\label{eq:Feps_crit}
\nabla W * \rho^\eps + \eps \frac{\nabla \rho^\eps}{\rho^\eps} = 0.
\end{equation}
Testing the above equation against $\rho^\eps x$ we obtain:
\[
\int \rho^\eps x \cdot \nabla W * \rho^\eps dx + \eps \int x \cdot \nabla \rho^\eps dx = 0.
\]
Integrating by parts and using Lemma \ref{lemma:symmetrization} with the odd vector field $\nabla W$ yields:
\[
\frac{1}{2} \iint \nabla W(x-y) \cdot (x-y) \rho^\eps(dx) \rho^\eps(dy) = \eps d.
\]
The desired result then follows from assumption (A4-1) and Lemma \ref{lemma:equivalence_Wp} with $p=1$, since $W_1(\rho^\eps, \delta_0) = M_1(\rho^\eps)$.
\end{proof}
Note that, from equation \eqref{eq:Feps_crit}, one has $\rho^\eps =C(\eps)e^{-W * \rho^\eps / \eps}$. The value of the constant $C(\eps)$ can be computed by imposing a total mass $1$, so that we get $\rho^\eps=\sfrac{e^{-W * \rho^\eps / \eps}}{\int e^{-W * \rho^\eps / \eps}}$. Using this equality along with estimate \eqref{W1_leq_Ceps}, we~obtain a bound in $W_p$ distance for $p \in [1, \infty)$ provided $W$ is also Lipschitz continuous:
\begin{theorem}
Assume that $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\textup{(A4-1)}. For any $p \in [1, \infty)$ there exists a constant $C > 0$ such that for any $\eps>0$ and $\rho^\eps$ steady state for \eqref{eq:rho_diff} which center of mass is 0, the following estimate holds:
\begin{equation*}%\label{Wp_leq_Ceps}
W_p(\rho^\eps, \delta_0) \leq C\eps.
\end{equation*}
\end{theorem}
\begin{remark}
At least in dimension one, this result is optimal. Indeed, we~can take for $W$ the Newtonian potential $W(x)=|x|$, for which, using the correspondence with Burgers' equation, $\rho^\eps$ can be written as $\rho^\eps(x) = \Psfrac{1}{\eps} \rho(\sfrac{x}{\eps})$, where $\rho(x) = \psfrac{1-\tanh^2\left(\sfrac{x}{2}\right)}{4}$, and a scaling argument then gives $W_p^p(\rho^\eps, \delta_0) = \eps^p M_p(\rho)$.
\end{remark}
\begin{proof}
Since
\[
\rho^\eps = \frac{e^{-W * \rho^\eps / \eps}}{\int e^{-W * \rho^\eps / \eps}},
\]
we have:
\[
W_p^p(\rho^\eps, \delta_0) = \frac{\int |x|^p e^{-W * \rho^\eps(x) / \eps} dx}{\int e^{-W * \rho^\eps / \eps}},
\]
Now, since $W$ is Lipschitz continuous, one has
\[
|W*\rho^\eps - W*\delta_0| \leq a_\infty \sup_{\Lip(\varphi) \leq 1} \int \varphi d(\rho^\eps - \delta_0) = a_\infty W_1(\rho^\eps, \delta_0) \leq C \eps,
\]
because of Theorem \ref{thm:W1_leq_Ceps}. Thus, $-W*\rho^\eps \leq C\eps - W$ and therefore:
\[
\int |x|^p e^{-W * \rho^\eps(x) / \eps} dx \leq C \int |x|^p e^{-W(x) / \eps}dx
\leq C\eps^{p+d} \int |y|^p e^{-W(\eps y) / \eps} dy,
\]
using the change of variables $x = \eps y$. Recall that Assumption (A4-1) ensures $W(x) \geq C|x|$ for all $x \in \RR^d$. This allows us to bound $\int |y|^p e^{-W(\eps y) / \eps} dy$ uniformly with respect to $\eps$.

On the other hand, since $W$ is $a_\infty$-Lipschitz continuous, we~have
\[
W(x) \leq a_\infty |x| + W(0) = a_\infty |x|.
\]
Integrating with respect to $\rho^\eps(dx)$, we~deduce that $W * \rho^\eps(0) \leq a_\infty W_1(\rho^\eps, \delta_0)$. Besides, $W * \rho^\eps$ is also $a_\infty$-Lipschitz continuous. Hence,
\[
W * \rho^\eps(x) \leq W * \rho^\eps(0) + a_\infty |x| \leq a_\infty W_1(\rho^\eps, \delta_0) + a_\infty |x| \leq C\eps + a_\infty |x|,
\]
thanks again to estimate \eqref{W1_leq_Ceps}. After another rescaling, we~deduce:
\[
\int e^{-W * \rho^\eps / \eps} \geq C \eps^d,
\]
thus getting $W_p^p(\rho^\eps, \delta_0) \leq C \sfrac{\eps^{p+d}}{\eps^d} = C\eps^p$, which concludes the proof.
\end{proof}

\Subsection{$O(\eps^{1/p})$ convergence rate in $W_p$ for potentials such that $\nabla W(x) \cdot x \geq C |x|^p$}

Assume $W$ satisfies assumptions \eqref{A0}, \eqref{A1} and \eqref{A4p} for some $p \in [1, \infty)$. Under this assumption, a straightforward adaptation of the proof of Theorem \ref{thm:W1_leq_Ceps} provides an estimate on $W_p(\rho^\eps, \delta_0)$:
\begin{theorem}\label{TWp}
Assume that $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A4p} for some $p \in [1, \infty)$. There exists a constant $C > 0$ such that for any $\eps>0$ and $\rho^\eps$ steady state for \eqref{eq:rho_diff} which is centered, the following estimate holds:
\begin{equation}\label{Wp_leq_Ceps_bis}
W_p(\rho^\eps, \delta_0) \leq C\eps^{1/p}.
\end{equation}
\end{theorem}
\begin{remark}
It is possible to prove optimality of this rate for $p=2$. Let us consider the quadratic potential $W(x) = |x|^2$, that satisfies assumption (A4-$2$). Recall that $\rho^\eps = \sfrac{e^{-W * \rho^\eps / \eps}}{\int e^{-W * \rho^\eps / \eps}}$. Expanding $W(x-y) = |x-y|^2$ and using both facts that the total mass of $\rho$ is 1 and that $\rho^\eps$ is centered, one has:
\begin{align*}
e^{-W * \rho^\eps / \eps} & = \exp\left\lbrace -\frac{1}{\eps} \left(\int|x|^2 \rho^\eps(y) dy - 2x \cdot \int y \rho^\eps(y) dy + \int|y|^2\rho^\eps(y) dy \right) \right\rbrace \\
& = e^{-|x|^2/\eps} e^{-W_2^2(\rho^\eps, \delta_0)/\eps}.
\end{align*}
Hence, $\rho^\eps(x) = \sfrac{e^{-|x|^2 / \eps}}{\int e^{-|x|^2 / \eps} dx}$, which in turn yields:
\[
W_2^2(\rho^\eps, \delta_0) = \frac{\int |x|^2 e^{-|x|^2 / \eps} dx}{\int e^{-|x|^2 / \eps} dx}.
\]
A change of variables in both integrals then gives $W_2^2(\rho^\eps, \delta_0) = C\eps$. Note the estimate $W_p^p(\rho^\eps, \delta_0) = C\eps^{p/2}$ can be proved in the same way, but is not relevant in the context of Theorem \ref{TWp}.
\end{remark}

\section{Numerical illustrations}\label{section:numerics}

This sections aims to illustrate our convergence results both in the evolutive case and in the stationary case. The implementation of the schemes has been done in Python and the code is available at \href{https://github.com/strantien/aggregation}{github.com/strantien/aggregation}. Tests are conducted on $[-1,1]$, with $2J+1$ cells, and the velocity field is always discretized by \eqref{def:aij}. Wasserstein distances between two arbitrary probability measures are computed using the POT package (see \cite{flamary2021pot}).

\begin{figure}[htb]
\centering
\includegraphics[scale=0.42]{aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_1.0_p_2_T_0.5_CFL_0.9_bis.png}
\caption{Order $1/2$ convergence in $W_2$ distance of $\rho^\eps_T$ towards $\rho_T$ for $\rho^{\ini}(x) = 2 \sqrt{\sfrac{5}{\pi}}\; e^{-20 x^2}$, $W(x)=|x|^2$.}
\label{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_1.0_p_2_T_0.5_CFL_0.9}
\end{figure}

\subsection{Evolutive solutions}

We begin with the convergence rate in Wasserstein distance of the viscous solutions $\rho^\eps$ associated with a fixed initial datum $\rho^{\ini}$ (not depending on $\eps$). In this subsection $\rho_{\dx}^\eps$ is computed using the implicit discretization \eqref{dis:implicit}, for which the CFL condition is less restrictive than the parabolic CFL condition of the explicit scheme. We~also implemented no-flux boundary conditions so as to preserve total mass. This condition is a discretization of
\[
(\eps \partial_x \rho - \achapo [\rho] \rho)(t,-1) = (\eps \partial_x \rho - \achapo [\rho] \rho)(t,1) = 0.
\]
Namely, \eqref{dis:implicit} is used for $j = 2, \dots, 2J$ and the system is closed with (here we omit the index $i$ as the space dimension is $1$, and we recall that $\theta = 1$)
\[
\rho_{1}^{n+1} = \rho_{1}^n - \frac{\dt}{\dx}
\big(({a}^n_{1})^+ \rho_{1}^n - ({a}^n_{2})^- \rho_{2}^n \big) + \eps \frac{\dt}{\dx^2}(\rho_{2}^{n+1} - \rho_1^{n+1}),
\]
and a similar equation for $j = 2J + 1$.
Since in some cases we do not have explicitly a reference solution for the inviscid equation ($\eps = 0$), the convergence rate with respect to $\eps$ is estimated taking $\dx$ small enough so that $\rho_{\dx}^\eps$ approximates $\rho^\eps$, and computing $W_p(\rho_{\dx, T}^{\eps_{i+1}}, \rho_{\dx, T}^{\eps_{i}})$: this quantity is called ``error'' in the $y$ axis on Figures \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_1.0_p_2_T_0.5_CFL_0.9}, \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_1_T_0.5_CFL_0.9} and \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_3_T_0.5_CFL_0.9}. Actually in the case of Figure \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_1.0_p_2_T_0.5_CFL_0.9} the velocity field has the form $- \nabla W * \rho (x) = -x$, which would allow for the computation of the reference solution; yet in order to use the same tools based on the POT package (more suitable for atomic measures) we do not exploit this property.

\begin{figure}[htb]
\centering
\includegraphics[scale=0.42]{aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_1_T_0.5_CFL_0.9_bis.png}
\caption{Order 1 convergence in $W_1$ distance of $\rho^\eps_T$ towards $\rho_T$ for $\rho^{\ini}(x) = 2 \sqrt{\sfrac{5}{\pi}}\; e^{-20 x^2}$, $W(x)=|x|$.}
\label{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_1_T_0.5_CFL_0.9}
\end{figure}

In Theorems \ref{thm:cvOrder1/2} and \ref{thm:cvOrder1/2:scheme}, when $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A3}, we~proved convergence at rate $O(\eps^{1/2})$ in $W_2$ distance, which is what we recover when~$W$ is smooth, as shows Figure \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_1.0_p_2_T_0.5_CFL_0.9}. In practice, for this test case, we~observe $O(\eps^{1/2})$ convergence rate in $W_p$ distance for any $p\in [1,+\infty[$. However, in case $W$ has a Lipschitz discontinuity at the origin (Figure \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_1_T_0.5_CFL_0.9}) we observe convergence at order~$1$ in~$W_1$ distance. This is the superconvergence phenomenon investigated by Tang, Teng and Zhang \cite{tangteng97, tengzhang97} in the framework of scalar conservation laws. In terms of aggregation, the interpretation is that, when $W$ is singular, the concentration is strong enough to compensate part of the diffusion. In other $W_p$ distances, convergence seems to occur at order 1 when $\eps$ is not too small, and then degenerates quite clearly towards order $1/p$ for any $p \in [1,+\infty[$ (see Figure \ref{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_3_T_0.5_CFL_0.9} for $p=3$). Note that, in every case, the convergence order is robust with respect to the test case (be it for smooth or singular initial data, e.g. Dirac masses).
\begin{figure}[htb]
\centering
\includegraphics[scale=0.42]{aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_3_T_0.5_CFL_0.9_bis.png}
\caption{Order $1/3$ convergence in $W_3$ distance, for small $\eps$, of $\rho^\eps_T$ towards $\rho_T$ for $\rho^{\ini}(x) = 2 \sqrt{\sfrac{5}{\pi}}\; e^{-20 x^2}$, $W(x)=|x|$.}
\label{fig:aggdiff_convergence_evolutif_upwind_diffImplicite_eps_J_5000_alpha_0.0_p_3_T_0.5_CFL_0.9}
\end{figure}

\subsection{Steady states}

In order to simulate the steady states for $\eps>0$, recall that they are characterized, over the whole space, by the following equation:
\begin{equation}\label{eq:fixedPoint}
\rho^\eps = \frac{e^{-W * \rho^\eps / \eps}}{\int e^{-W * \rho^\eps / \eps}}.
\end{equation}
We therefore use a fixed-point method on the map sending $\rho$ into $\sfrac{e^{-W * \rho / \eps}}{\int e^{-W * \rho / \eps}}$ in order to solve Equation \eqref{eq:fixedPoint}. Fixed point algorithm is stopped as soon as the $W_p$ distance between two iterations is below some tolerance. Numerically, we~observe that this method turns two symmetric Gaussian bumps almost immediately (after the first iteration) into a centered Gaussian whenever $W$ is attractive and Lipschitz.

\begin{figure}[htb]
\centering
\includegraphics[scale=0.42]{aggdiff_convergence_stationnaire_upwind_diffImplicite_potentielMixte_eps_J_80000_p_1_tol_1e-12.png}
\caption{Order of convergence in $W_1$ distance of $\rho^\eps$ towards $\delta_0$, for the non-Lipschitz potential $W(x) = \sqrt{|x|}+|x|$. The initial density for the fixed point algorithm is the centered Gaussian $2 \sqrt{\sfrac{5}{\pi}}\; e^{-20 x^2}$.}
\label{fig:aggdiff_convergence_stationnaire_upwind_diffImplicite_potentielMixte_eps_J_80000_p_1_tol_1e-12}
\end{figure}
We first investigate the convergence rate towards the Dirac mass, for centered steady states. The error is estimated computing the integral $\int |x|^p\rho(dx) = W_p^p(\rho, \delta_0)$. When $W$ satisfies assumptions \eqref{A0}--\eqref{A1}--(A4-1), we~proved $O(\eps)$ convergence rate in~$W_1$ distance, which we do recover in Table \ref{table} for $W(x)=|x|$. We~also explore the case when $W$ verifies \eqref{A0}--\eqref{A1}--(A4-1) but is not Lipschitz continuous, which is the case of $W(x) = \sqrt{|x|} + |x|$. For this potential, we~obtain, in Figure \ref{fig:aggdiff_convergence_stationnaire_upwind_diffImplicite_potentielMixte_eps_J_80000_p_1_tol_1e-12} convergence at order 1.82264413 which is slightly less than 2, in $W_1$ distance. This can be linked to the fact that $W$ satisfies a sort of assumption (A4-$\sfrac 12$) when $|x| \leq 1$.
Under assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--(A4-3), we~observe convergence at rate $1/3$ in $W_3$ distance as we proved in \eqref{Wp_leq_Ceps_bis}, as shows Figure \ref{fig:aggdiff_convergence_stationnaire_upwind_diffImplicite_alpha_2.0_eps_J_1000_p_3_tol_1e-06}.

\begin{figure}[htb]
\centering
\includegraphics[scale=0.46]{aggdiff_convergence_stationnaire_upwind_diffImplicite_alpha_2.0_eps_J_1000_p_3_tol_1e-06.png}
\caption{$O(\eps^{1/3})$ convergence in $W_3$ distance of $\rho^\eps$ towards $\delta_0$, $ W(x) = |x|^3$. The initial density for the fixed point algorithm is the centered Gaussian $2 \sqrt{\sfrac{5}{\pi}}\; e^{-20 x^2}$.}
\label{fig:aggdiff_convergence_stationnaire_upwind_diffImplicite_alpha_2.0_eps_J_1000_p_3_tol_1e-06}
\end{figure}

More generally, under assumptions \eqref{A0}--\eqref{A1}--\eqref{A2}--\eqref{A4p}, convergence at rate $1/p$ seems to occur in any $W_q$ distance, $q\in [1,+\infty[$, which is what we proved in for $p=1$ or for $p=q$. To illustrate this latter case, we~compute the convergence order in $W_p$ distance for $W(x)=|x|^p$, which seems indeed to be $1/p$, see Table \ref{table} (when $p=1$, since the potential is pointy, one has to refine the mesh so as to observe proper convergence at order $1$).

\begin{table}[htb]
\centering
\caption{Convergence order $\simeq \sfrac 1p$ of $\rho^\eps$ towards $\delta_0$ for $W(x)=|x|^p$, $tol = 10^{-6}$, $\eps_i = 2^{-i}$, $i=4,\ldots,16$, initial density $2 \sqrt{\sfrac{5}{\pi}}\; e^{-20 x^2}$}
\label{table}
\begin{tabular}{c|c|c}
$p$ & Order & $J$ \\
\hline
1 & 1.00205259 & 50000 \\
2 & 0.49999997 & 2000 \\
3 & 0.33333333 & 2000 \\
4 & 0.25000000 & 2000 \\
5 & 0.20000000 & 2000
\end{tabular}
\end{table}

\backmatter
\bibliographystyle{jepplain+eid}
\bibliography{lagoutiere-et-al}
\end{document}