Aggregation-Diffusion Equations with saturation
Workshop on ‘New Perspectives in Nonlocal and Nonlinear PDEs’, Anacapri
Collaborators
\[ \newcommand{\diff}{d} \newcommand{\dr}{\diff r} \newcommnand{\R}{{\mathbb R}} \newcommand{\Rd}{{\mathbb R^d}} \newcommand{\diver}{\nabla \cdot} \newcommand{\calF}{\mathcal{F}} \newcommand{\calU}{\mathcal{U}} \newcommand{\defeq}{\overset{\mathrm{def}}{=}} \newcommand{\ee}{\varepsilon} \newcommand{\mob}{\mathrm{m}} \newcommand{\mobee}{\mathrm{m}_\varepsilon} \newcommand{\Wass}{\mathfrak{W}} \newcommand{\inversedU}{T_{0,\alpha} \circ (U')^{-1}} \newcommand{\inversedUee}{(U_{\ee}')^{-1}} \newcommand{\mobup}{\mob_{\uparrow}} \newcommand{\mobdown}{\mob_{\downarrow}} \newcommand{\mobupwind}{\mob_{\textrm{w}}} \newcommand{\defeq}{=} \newcommand{\Rho}{P} \newcommand{\bm}[1]{\mathbf{#1}} \newcommand{\supp}{\operatorname{supp}} \DeclareMathOperator*{\argmin}{arg\,min} \]
Aggregation-Diffusion
\[\begin{equation*} \label{eq:ADE} \tag{ADE$_\mob$} \frac{\partial \rho}{\partial t} = \diver \Bigg( \mob(\rho) \nabla \Big( U'(\rho) + V + \int K(\cdot, y) \rho(y) dy \Big) \Bigg). \end{equation*}\] The mobility \(\mob\) is usually linear, i.e., \(\mob(\rho) = \rho\).
For the introductory part we will follow [Gómez-Castro, 2024].
Continuity equations
Let \(\rho\) be a density and \(\omega \subset \Rd\) any control volume, if \(\mathbf F\) is the out-going flux \[\begin{equation} \label{eq:conservation control volume} \frac{\diff }{\diff t} \int_\omega \rho \diff x =- \int_{\partial \omega} \mathbf F \cdot \mathbf n \diff S= -\int_{ \omega} \diver \mathbf F \diff x \end{equation}\] Hence we arrive at the \[\begin{equation} \label{eq:conservation law} \frac{\partial \rho}{\partial t} = -\diver \mathbf F. \end{equation}\]
Diffusion
For the transport of heat, we use Fourier’s law to model the flux \(\mathbf F = -D \nabla \rho\) and yields the heat equation \[\begin{equation*} \tag{HE} \label{eq:HE} \frac{\partial \rho}{\partial t} = D \Delta \rho . \end{equation*}\]
Using Darcy’s law we get Porous Medium Equation [Vázquez, 2006] \[\begin{equation*} \tag{PME} \label{eq:PME} \frac{\partial \rho}{\partial t} = \Delta \rho^m. \end{equation*}\]
Particle systems
We can also understand transport from the particle perspective. Consider a known velocity field \(\mathbf v(t,x)\). Consider \(N\) particles with positions \(X_i (t)\) of equal masses \(1/N\) moving in the velocity field \[\begin{equation*} % \tag{AE$_N$} \frac{\diff X_i}{\diff t} = \mathbf v(t, X_i(t)) , \qquad i = 1, \cdots, N. \end{equation*}\] Define the empirical distribution \[\begin{equation} \label{eq:empirical distribution} \mu_t^N = \sum_{j=1}^N \frac 1 N \delta_{X_j(t)} , \end{equation}\] where \(\delta_x\) is the Dirac delta at a point \(x\). This is a weak solution to the continuity equation \[\begin{equation*} \partial_t \mu^N + \diver (\mu^N \mathbf v ) = 0. \end{equation*}\]
Aggregation and Confinement
Consider the particles interacting non-locally by \[\begin{equation*} % \tag{AE$_N$} \frac{\diff X_i}{\diff t} = - \sum_{\substack{ j=1 \\ j \ne i}}^N \frac 1 N \nabla W (X_i - X_j) - \nabla V (X_i) , \qquad i = 1, \cdots, N. \end{equation*}\]
Notice that \[ \sum_{\substack{ j=1 \\ j \ne i}}^N \frac 1 N \nabla W (x - X_j(t)) = \int_\Rd \nabla W (x-y) \diff \mu^N_t (y). \]
The empirical measure \(\mu^N\) is a distributional solution of the Aggregation-Confinement Equation \[\begin{equation} \partial_t \mu = \diver (\mu \nabla ( W * \mu + V ) ). \end{equation}\]
We can replace \(W * \rho\) by \(\int K(\cdot, y) \rho(y) dy\).
No-flux conditions
In a bounded domain \(\Omega\), we preserve mass if we set \(\mathbf F \cdot \mathbf n = 0\) on \(\partial \Omega\), i.e., \[\begin{equation*} \tag{BC} \label{eq:BC Omega} \rho \nabla ( U'(\rho) + V + W*\rho ) \cdot \mathbf n = 0 \qquad \text{ for } t > 0 \text{ and } x \in \partial \Omega \end{equation*}\] (and set \(\rho = 0\) in \(\Rd \setminus \Omega\) for the convolution).
Non-linear mobility
There are several motivations to discuss the problem with non-linear mobility \[\begin{equation*} \tag{ADE$_{\mob}$} \label{eq:ADE mob} \begin{aligned} &\frac{\partial \rho}{\partial t} = \diver \Big( \mob(\rho) \nabla ( U'(\rho) + V + \mathcal K \rho ) \Big), \\ &\quad \text{where } \mathcal K \rho (x) = \int_\Omega K(x,y) \rho(y) \diff y. \end{aligned} \end{equation*}\]
The power-type case \(\mob(\rho) =\rho^{m-1}\) was made popular by Caffarelli and Vázquez
We will discuss the case of saturation \(\mob(0) = \mob(\alpha) = 0\).
Free energy, gradient-flow structure, and asymptotics
Free energy: a Lyapunov functional
\(\eqref{eq:ADE mob}\) is related to the free energy \[\begin{equation*} \tag{FE} \label{eq:free energy star} \mathcal F(\rho) =\int_\Omega U(\rho) + \int_\Omega V \rho + \frac 1 2 \iint_{\Omega \times \Omega } K(x,y) \rho(x) \rho(y) \diff x \diff y . \end{equation*}\] For convenience, let us define the first variation \[\begin{equation} \label{eq:first variation} \frac{\delta \mathcal F}{\delta \rho} [\rho] \defeq U'(\rho) + V + \int_\Omega K(\cdot,y)\rho(y) \diff y. \end{equation}\]
We write \(\eqref{eq:ADE}\) as \[\begin{equation*} \partial_t \rho = \diver\left( \mob(\rho) \nabla \frac{\delta \mathcal F}{\delta \rho} [\rho]\right) \end{equation*}\]
Then, \[\begin{equation*} \frac{d}{dt} \mathcal F[\rho_t] = \int \frac{\delta \mathcal F}{\delta \rho} [\rho] \partial_t \rho_t = - \int \mob(\rho_t) \left|\frac{\delta \mathcal F}{\delta \rho} [\rho] \right|^2. \end{equation*}\]
Gradient-flow structure
When \(\mob(\rho) = \rho\) there a nice gradient-flow theory:
Otto calculus (see, e.g., [Ambrosio, Brué & Semola, 2021]) allows us to show that \(\eqref{eq:ADE}\) in \(\Rd\) is formally the \(2\)-Wasserstein gradient flow in the sense that \[ \frac{\partial \rho}{\partial t} = - \nabla_{\Wass_2} \mathcal F[\rho_t] \]
There is a good notion of convexity (McCann’s displacement convexity). Under certain condition we have \[ \Wass_2(\mu_t, \widehat \mu_t) \le e^{-\lambda t} \Wass_2(\mu_0, \widehat \mu_0). \]
Wasserstein-type space with general mobility
For general \(\mob(\rho)\), there is a generalisation
[Dolbeault, Nazaret & Savaré, 2009], [Carrillo, Lisini, Savaré & Slepčev, 2010]
through a generalised Benamou-Brenier using \(\partial_t \rho + \diver( \mob(\mu) v ) = 0\).
This is only a distance if \(\mob\) is concave.
It is (almost) never displacement convex in our settings.
Formal steady states
\[\begin{equation} \tag{ADE$_\mob$} \frac{\partial \rho}{\partial t} = \diver \Bigg(\mob(\rho) \nabla \Big( \underbrace{ U'(\rho) + V + \int K(\cdot, y) \rho(y) dy }_{\frac{\delta \mathcal F}{\delta \rho}} \Big) \Bigg). \end{equation}\]
Formally, to get a steady state it suffices that \[ \frac{\delta \mathcal F}{\delta \rho} = C_\omega \text{ for each }\omega\text{ connected component of }\supp \rho \] .
Euler-Lagrange condition for \(\mob(\rho) = \rho\)
Recall \(\frac{\delta \calF}{\delta\rho} [\rho] = U'(\rho) + V + \int K(\cdot, y) \rho(y) dy\).
We have that \[\begin{equation*} \widehat \rho \in \argmin_{\mathcal A_M} \mathcal F \implies \exists C \text{ such that } \begin{dcases} \frac{\delta \calF}{\delta\rho} [\widehat \rho] = C & \text{a.e. in } \supp \rho, \\ \frac{\delta \calF}{\delta\rho} [\widehat \rho] \ge C & \text{a.e. where } \rho = 0. \\ \end{dcases} \end{equation*}\]
The free-energy minimiser is the steady state with same constant in all connected components.
\[ \widehat \rho(x) = \max\Bigg\{0, (U')^{-1} \big(C-V(x) - W*\widehat \rho\big)\Bigg\}. \]
Notice that when \(W = 0\) we the solution is parametrised by \(C\).
Many steady states, only one free-energy minimiser (\(K = 0\))
You don’t always converge to the global minimizer
Formal steady states in the saturation case
We consider \(\mob(0) = \mob(\alpha) = 0\).
Numerical results obtained in [Bailo, Carrillo & Hu, 2023a] suggested the formation of free boundaries at levels \(0, \alpha\).
Furthermore, their numerical results suggest that the steady states are the doubly-truncated Barenblatts \[ \widehat \rho = \min \Bigg\{ \alpha, \max\Big\{ 0, (U')^{-1}(C - V - W*\widehat \rho) \Big\} \Bigg \} \]
In [Carrillo, Fernández-Jiménez & Gómez-Castro] we discuss this for \(W = 0\).
Analysis for drift-diffusion with saturation
We discuss \(\eqref{eq:ADE mob}\) when \(K = 0\) and \(\mob(0) = \mob(\alpha) = 0\). \[\begin{equation} \tag{DDE} \label{eq:the problem Omega} \frac{\partial \rho}{\partial t} = \diver \left( \mob (\rho) \nabla \left( U' (\rho) + V \right) \right) \end{equation}\] in a bounded domain \(\Omega\) with no-flux condition.
We present the results in [Carrillo, Fernández-Jiménez & Gómez-Castro].
Assumptions
We assume that the diffusion is continuous, in the sense that \[\begin{equation} \label{hyp:Phi' in L1} \mob U'' \in L^1(0,\alpha), \end{equation}\] and we define \[ \Phi(s) \coloneqq \int_{s_0}^s \mob(\tau) U''(\tau) \diff \tau. \] Furthermore, we also assume that \(\Phi\) is strictly increasing at \(0\) and \(\alpha\), i.e., \[\begin{equation} \label{eq:Assumption Existence} \Phi(0) < \Phi (s) < \Phi (\alpha) \quad \text{for all } s \in (0, \alpha). \end{equation}\] Lastly, we impose a technical regularity condition which will be suitable for compactness estimates \[\begin{align} \label{eq:hypothesis Phi/Phi' ee = 0} \sup_{s \in [0,\alpha] }\left|\frac{(\Phi (\alpha)-\Phi (s))\Phi(s)}{\Phi'(s)} \right| + \left|\frac{(\Phi (s)-\Phi (0))\Phi(s)}{\Phi'(s)} \right| & < \infty. \end{align}\]
To study the steady states we sometimes assume some of the following strict convexity to different degrees \[\begin{align} \tag{SC$_U$} \label{eq:U locally strictly convex} & U''(s) >0 , \qquad \text{for a.e. } s \in (0,\alpha) . \\ \tag{USC$_U$} \label{eq:U uniformly strictly convex} \inf_{s \in (0,\alpha)} & U''(s) >0 . \end{align}\]
For certain statements on numerical schemes, we will assume \[\begin{equation} \label{eq: U is C1[0,alpha]} U \in C^1([0,\alpha]). \end{equation}\]
Asymptotics diagram
We introduce regularised problems and we prove the following diagram 
Domains
Let \[\begin{align*} \mathcal{A} & \coloneqq \left\{ \rho \in L^1 (\Omega) : 0 \leq \rho \leq \alpha \right\}, \\ \mathcal{A}_+ & \coloneqq \left\{ \rho \in L^1 (\Omega) : \exists \delta > 0 \text{ s.t. } \delta \leq \rho \leq \alpha - \delta \right\}, \\ \mathcal A_M & \coloneqq \left\{ \rho \in \mathcal A : \|\rho\|_{L^1(\Omega)} = M \right\} . \end{align*}\]
We say that \(S_t : \mathcal A \to \mathcal A\) is a free-energy dissipating semigroup of solutions if
For \(\rho_0 \in \mathcal A\), \(\rho_t = S_{t}\rho_0\) is a weak solution
\(S_t\) is a \(C_0\)-semigroup in \(L^1\), i.e., for \(t, h > 0\) we have \[ S_{t+h} = S_t S_h, \qquad \lim_{t \to 0^+}\| S_t \rho_0 - \rho_0 \|_{L^1 (\Omega)} = 0 \text{ for all } \rho_0 \in \mathcal A. \]
\(S_t : \mathcal A \to \mathcal A\) is an \(L^1\)-contraction, i.e., for any \(\rho_0, \eta_0 \in \mathcal A\) we have that \(\|S_t \rho_0 - S_t \eta_0 \|_{L^1(\Omega)} \le \|\rho_0 - \eta_0\|_{L^1(\Omega)}.\)
Free-energy dissipation and \(C^{\frac 1 2}_{loc}([0,\infty), W^{-1,1}(\Omega))\) continuity: If \(\rho_0 \in \mathcal A_+\) then calling \(\rho_t = S_t \rho_0\) we have \[\begin{align} \nonumber \text{For all } 0 < t_1 < t_2 \text{ we get } \qquad & \\ \label{eq:free energy dissipation} \int_{t_1}^{t_2} \int_\Omega \mob(\rho_\sigma) |\nabla(U'(\rho_\sigma) + V)|^2 & \le \mathcal F[\rho_{t_1}] - \mathcal F[\rho_{t_2}], \\ \label{eq:W-11 continuity} \| \rho_{t_2} - \rho_{t_1} \|_{W^{-1,1}(\Omega)} & \leq \| \mob \|_{L^\infty(0,\alpha)}^{\frac{1}{2}} | \Omega |^{\frac{1}{2}} \left( \mathcal{F} [\rho_{t_1}] - \mathcal{F} [\rho_{t_2}] \right)^{\frac{1}{2}} |t_2-t_1|^{\frac{1}{2}}. \end{align}\] In particular, \(t \mapsto \mathcal F[S_t \rho_0]\) is non-increasing.
Approximating problem
\[\begin{equation} \tag{DDE$_\ee$} \label{eq:the problem regularised} \frac{\partial \rho}{\partial t} = \Delta \Phi_{\varepsilon} (\rho ) + \diver \left( \mobee (\rho ) \nabla V \right) . \end{equation}\] with \(0 < c_\ee \le \Phi_\ee' \le C_\ee\).
If \(\rho_0 \in \mathcal{A}_+ \cap C^2 (\overline{\Omega})\), then problem \(\eqref{eq:the problem regularised}\) has a unique classical solution.
These classical solutions can be uniquely extended to \(S_t^{(\ee)}\), a free-energy dissipating semi-group.
If \(\rho_0 \in \mathcal A \setminus \{0,\alpha\}\) then \(0 < S_t^{(\ee)} \rho < \alpha\) in \(\Omega\) for \(t > 0\).
\(S_t^{(\ee)} : \mathcal A_+ \to \mathcal A_+\).
Exist \(\ee_k \to 0\) and \(S\) such that
\(S^{(\ee_k)} \rho_0 \to S \rho_0\) in \(C_{loc} ( [0,\infty); L^1( \Omega ) )\) for all \(\rho_0 \in \mathcal A\).
Euler-Lagrange condition
If \(\widehat \rho\) is a local minimiser of \(\mathcal F\) on \(\mathcal A_M\) with the \(L^1\) topology, then there exists \(C \in \mathbb R\) such that \[\begin{equation} \label{eq:EL for ee = 0} \begin{aligned} U'(\widehat \rho(x)) + V(x) \ge C, & \qquad \text{ for a.e. } x \text{ such that } 0 \le \widehat \rho(x) <\alpha. \\ U'(\widehat \rho(x)) + V(x) \le C, & \qquad \text{ for a.e. } x \text{ such that } 0 < \widehat \rho(x) \le\alpha . \end{aligned} \end{equation}\] Furthermore, if \(U'\) is invertible \[\begin{equation} \tag{EL$_P$} \label{eq:Euler-Lagrange for P} \widehat \rho (x)= \inversedU (C- V(x)) \quad \text{a.e.~in } \Omega. \end{equation}\] Lastly, if we assume \(\eqref{eq:U locally strictly convex}\) and \(M \in (0,\alpha|\Omega|)\), there exists a unique \(C\) such that \(\eqref{eq:Euler-Lagrange for P}\) has mass \(M\).
Time-asymptotics via semi-groups
We say that a semigroup \(S\) for a problem \(\eqref{eq:the problem Omega}\) has a time limit operator \(S_\infty : \mathcal A \to \mathcal A\) if:
For any \(\rho_0 \in \mathcal A\) there exists a limit in time \[ S_t \rho_0 \to S_\infty \rho_0 \qquad \text{ strongly in } L^1(\Omega) \text{ as } t \to \infty. \]
\(S_\infty\) is stationary for the semigroup, i.e., \(S_t S_\infty = S_\infty\).
For any \(\rho_0 \in \mathcal A\), \(S_\infty \rho_0\) is a constant-in-time weak solution to \(\eqref{eq:the problem Omega}\).
Time-asymptotics via semi-groups
For \(\ee >0\) the semigroup for \(\eqref{eq:the problem regularised}\) has a time-limit operator, which we denote \(S_\infty^{(\ee)}\).
For \(\rho_0 \in \mathcal A_M\) we have \(S_\infty^{(\ee)} \rho_0 = \widehat\rho^{(\ee)}\) \[\begin{equation} \tag{A$_\varepsilon$} \widehat\rho^{(\ee)} (x) \coloneqq \inversedUee \left( C_{\varepsilon} - V(x) \right), \quad \text{in } \Omega , \end{equation}\] where \(C_\ee\) is uniquely determined by the mass condition \(\int_\Omega \inversedUee( C_{\varepsilon} - V) = M.\)
Any free energy-dissipating semi-group \(S\) for \(\eqref{eq:the problem Omega}\) has a time-limit operator \(S_\infty\).
Both \(S^{(\ee)}_\infty\) and \(S_\infty\) are \(L^1\)-contractions.
As for linear mobility, the \(\omega\)-limit \(\{S \rho_0 : \rho_0 \in \mathcal A_M\}\) may not a single-ton, and may even be a continuous manifold of steady states.
Free-energy minimiser of \(\eqref{eq:the problem Omega}\)
Assume \(\eqref{eq:U locally strictly convex}\) and that \(M \in (0,\alpha|\Omega|)\). Then we can define \[\begin{equation} \tag{A} \widehat\rho^{(0)} (x) \coloneqq \inversedU (C_0 - V(x) ), \quad \text{in } \Omega, \end{equation}\] where \(C_0\) is uniquely determined by the mass condition \(\int_{\Omega} \inversedU (C_0 - V ) = M.\) We have that:
\(S_t \widehat\rho^{(0)} = \widehat\rho^{(0)}\) for any \(t > 0\).
\(\widehat\rho^{(0)}\) is the unique \(L^1\)-local minimiser of the free energy \(\mathcal F\) over \(\mathcal A_M\).
\(\widehat\rho^{(0)}\) is the unique global minimiser over \(\mathcal A_M\).
If we also assume \(\eqref{eq:U uniformly strictly convex}\), then \(\widehat{\rho}^{(\ee)} \rightarrow \widehat{\rho}^{(0)} \quad \text{in } L^1(\Omega ) \text{ as } \ee \rightarrow 0\).
A convergent numerical scheme
An up-winded numerical scheme
Let \(\mobupwind(a,b) = \mobup(a) \mobdown(b)\) and let
\[\begin{equation} \label{eq:scheme} \begin{aligned} \frac{\Rho_{\bm i}^{n+1} - \Rho_{\bm i}^n}{\tau} &= - \sum_{k=1}^d \frac{F_{\bm i + \frac 1 2 \bm e_k}^{n+1} - F_{\bm i - \frac 1 2 \bm e_k}^{n+1} }{h}, \\ F_{\bm i + \frac 1 2 \bm e_k}^{n+1} & \defeq \mobupwind(\Rho_{\bm i}^{n+1} , \Rho_{\bm i + e_k}^{n+1} ) (v_{\bm i + \frac 1 2 \bm e_k}^{n+1})_+ + \mobupwind(\Rho_{\bm i + e_k}^{n+1}, \Rho_{\bm i}^{n+1}) (v_{\bm i + \frac 1 2 \bm e_k}^{n+1})_- \\ v_{\bm i + \frac 1 2 \bm e_k}^{n+1} & \defeq -\frac{\xi_{\bm i + \bm e_k}^{n+1}-\xi_{\bm i}^{n+1}}{h} \\ \xi_{\bm i}^{n+1} & \defeq U'(\Rho_{\bm i}^{n+1}) + V_{\bm i} + h^d \sum_{\bm j \in \bm I} K_{\bm i, \bm j} \Rho_{\bm j}^{n + \frac 1 2}, \\ \bm P^{n+\frac 1 2} &\defeq \frac{\bm P^{n+1} + \bm P^n}{2}. \end{aligned} \end{equation}\] where we use the convention \(a_+ = \max\{a,0\}\) and \(a_- = \min\{a,0\}\) so \(a = a_+ + a_-\).
We set no-flux conditions on the boundary of numerical domain.
Previous results
Linear mobility:
- Properties of the scheme [Bailo, Carrillo & Hu, 2020],
- Convergence by Fréchet-Kolmgorov [Bailo, Carrillo, Murakawa & Schmidtchen, 2020]
The case \(\mob(\rho) = \rho \mobdown (\rho)\)
- Properties of the scheme [Bailo, Carrillo & Hu, 2023b]
\(K = 0\)
[Carrillo, Fernández-Jiménez & Gómez-Castro]
- Existence of numerical solutions: topological degree argument
- The scheme is monotone
- Comparison principle
- \(L^1\)-contraction
- Decay of a discrete free energy
- Convergence when the solution is smooth
- Asymptotic behaviour of the numerical scheme \(n \to \infty\)
\(K \ne 0\)
[To appear]
- Convergence of solutions by Aubin-Simons theorem [Gallouët & Latché, 2012] with discrete versions of
- \(u \in L^2 (0,T; H^1)\)
- \(\partial_t u \in L^2(0,T; W^{-1.1})\)
Extension and open problems
Upcoming
Porous-Medium with nonlocal pressure \(v = -(-\Delta)^{-s} \rho\) (i.e., \(K = |x|^{-d+2s}\)), \(U = 0, V = 0\)
Work in preparation with F. del Teso and E. Jakobsen w
- \(u \in L^p(0,T; W^{s,p})\) and \(\partial_t u \in L^2(0,T; W^{-1.1})\).
Open problems
Complete characterisation of steady states with free boundary
\(C^\alpha\) regularity in general settings [Alamri, 2025]
Numerics on non-square meshes (e.g., Voronoi)
Questions?
