A finite-volume scheme for aggregation-diffusion with non-linear mobility

PDE Seminar @ U. Oxford

Author

Affiliation

David Gómez-Castro

Universidad Autónoma de Madrid - ICMAT

Aim of the talk

\[ \newcommand{\mob}{\mathsf m} \DeclareMathOperator{\sign}{sign} \]

Construction, key properties, and convergence of a finite-volume scheme for the aggregation-diffusion equation \[ \partial_t \rho = \operatorname{div} (\mob(\rho) \nabla (U'(\rho) + V + W*\rho)). \] with no-flux conditions

Surveys: [Gómez-Castro, 2024], [Bailo, Carrillo & Gómez-Castro, 2026]

Results: [Carrillo, Fernández-Jiménez & Gómez-Castro, 2025], [Gómez-Castro, 2025]

Construction of the scheme

Finite-volume schemes

Consider a conservation equation in \(d=1\) \[ \partial_t \rho = - \partial_x F \] for \((t,x) \in (0,T)\times (-\infty, \infty)\).

Then we have \(\displaystyle \partial_t \int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} \rho(t,x) dx = - \Big(F(t,x_{i+\frac 1 2}) - F(t,x_{i-\frac 1 2})\Big).\)

Finite-volume numerical schemes for approximated quantities \(h_{i} = x_{i+\frac 1 2}-x_{i-\frac 1 2}\) \[\begin{equation*} \rho_{i}(t) \approx \frac 1{h_{i}}\int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} \rho(t,x) dx, \qquad F_{i + \frac 1 2}(t) \approx F(t, x_{i+\frac 1 2}). \end{equation*}\]

The schemes are of the form \(\displaystyle\frac{\partial \rho_{i}}{\partial t} (t) = - \frac{F_{i+\frac 1 2}(t) - F_{i-\frac 1 2}(t)}{h_{i}}.\)

Mass conservation

We have that \[\begin{equation*} \partial_t \int_{x_{\frac 1 2}}^{x_{N+\frac 1 2}} \rho(t,x) dx = - \Big(F(t,x_{N+\frac 1 2}) - F(t,x_{\frac 1 2})\Big). \end{equation*}\]

In discrete form \[\begin{equation*} \frac{\partial}{\partial t} \left(\sum_{i=1}^{N} h_i \rho_{i} (t)\right) = - {F_{N+\frac 1 2}(t) - F_{\frac 1 2}(t)}. \end{equation*}\]

Up-winding. Transport with constant velocity

The transport equation for \(a \in \mathbb R\) we take flux \(F(t,x) = a \rho(t,x)\), i.e. \[ \partial_t \rho = -\partial_x(a \rho) . \]

The solutions are \[ \rho_t (x) = \rho_0(x - at). \] The solution is moving with the “wind”, and you have look “up-winding” for the information. The solution “looks” left if \(a > 0\) and right if \(a > 0\).

Up-winding. Discrete flux

Naive choices of discrete flux are \[ F_{i + \frac 1 2}(t) = a \rho_{i}(t) \qquad \text{or} \qquad F_{i + \frac 1 2}(t) = a \rho_{i + 1} (t) \] Lead to the schemes \[ {\partial_t \rho_{i}} = - a \frac{\rho_{i}(t) - \rho_{i-1}(t)}{h_i} \qquad {\partial_t \rho_{i}} = - a \frac{\rho_{i+1}(t) - \rho_{i}(t)}{h_i} \]

Since we want out solution bring information from the “correct” side \[ F_{i + \frac 1 2}(t) = \begin{dcases} a \rho_{i}(t) & \text{if } a > 0 \\  a \rho_{i + 1} (t) & \text{if } a< 0. \end{dcases} \]

Up-winding. Stability and positivity

Let us solve the problem in \(\mathbb R\) with \(F_{i+\frac 1 2}(t) = a \rho_{i}(t)\). Using standard \[ {\partial_t \rho_{i}} + \frac{a}{h_i} \rho_i(t) = \frac{a}{h_i} \rho_{i-1}(t) . \] then we have that \[ {\rho_{i}}(t) = e^{-t\frac{a}{h_i}} \rho_i(0) +\int_0^t e^{ -(t-s) \frac{a}{h_i} } \frac{a}{h_i} \rho_{i-1}(s) ds. \]

• If \(a > 0\):
  • \(\rho_i(0) \ge 0\) for all \(i \implies \rho_i(t) \ge 0\) for all \(i\) and \(t > 0\).
  • \(\rho_i(0) \le \overline \rho_i(0)\) for all \(i \implies \rho_i(t) \le \overline \rho_i(t)\) for all \(i\) and \(t > 0\)
• If \(a < 0\) we take \(\rho_i(0) = \delta_{i0} \implies \rho_{\pm 1}(t) < 0\) for \(t \sim 0^+\).

Up-winding with variable velocity

If \(F(t,x) = v(t,x) \rho(t,x)\), and we have some discretization of the velocity \(v_{i+\frac 1 2}(t)\), then we take the up-winding
\[ F_i(t) = \begin{dcases} \rho_i(t) v_{i+\frac 1 2} (t) & \text{if } v_{i + \frac 1 2} \ge 0, \\ \rho_{i+1}(t) v_{i+\frac 1 2} (t) & \text{if } v_{i + \frac 1 2} < 0 \end{dcases} \] It is common to write this concisely with the convention $a_+ = {a,0} $ and \(a_- = \min\{a,0\} \le 0\) as \[ F_{i+\frac 1 2} (t) = \rho_i(t) \left(v_{i+\frac 1 2} (t)\right)_+ + \rho_i(t) \left(v_{i+\frac 1 2} (t)\right)_-. \]

Higher spacial dimensions

Uniform grids

\(x_{\mathbf i} = (i_1h_1, i_2 h_2, \cdots, i_d h_d)\) with \(\mathbf i \in \mathbb Z^d\).

Then given a cube \(x_{\mathbf i} + Q_{\mathbf h}\) where \[ Q_{\mathbf h} = \prod_{k=1}^d \left(-\frac {h_k} 2, \frac {h_k} 2\right) \] there is a normal flux on each face, which we denote \(F_{\mathbf i \pm \frac 1 2 e_k}\).

Non-uniform meshes

The domain is decomposed into polyhedra \(K\), and the fluxes over faces are denoted by \(F_{K|L}\). This is a very general approach that is nicely explained in [Eymard, Gallouët & Herbin, 1997]

The scheme for linear mobility [Bailo, Carrillo & Hu, 2020]

\[ \tag{ADE} \partial_t \rho = \operatorname{div} (\rho \nabla (U'(\rho) + V + W*\rho)). \]

.

Continuous	Discrete
\(\partial_t \rho = - \operatorname{div} \mathbf F\)	\(\displaystyle \frac{\rho_{\mathbf i}^{n+1} - \rho_{\mathbf i}^n}{\Delta t} = -\sum_{k=1}^d \frac{F_{\mathbf i + \frac 1 2 \mathbf e_k}^{n+1}-F_{\mathbf i - \frac 1 2 \mathbf e_k}^{n+1}}{h_k}\)
\(\mathbf F = \rho \mathbf v\)	\(F_{\mathbf i + \frac 1 2 \mathbf e_k}^{n+1} = \rho_{\mathbf i}^{n+1} (v_{\mathbf i + \frac 1 2 \mathbf e_k})_+ + \rho_{\mathbf i + \mathbf e_k}^{n+1} (v_{\mathbf i + \frac 1 2 \mathbf e_k})_-\)
\(\displaystyle \mathbf v = -\nabla\xi\)	\(v_{\mathbf i + \frac 1 2 \mathbf e_k} = -\frac{\xi_{\mathbf i + \mathbf e_k}-\xi_{\mathbf i}}{h_k}\)
\(\xi = U'(\rho) + V + W*\rho\)	\(\xi_{\mathbf i}^{n+1} = U'(\rho_{\mathbf i}^{n+1}) + V(x_{\mathbf i}) + (W \star \rho)^{n+1}_{\mathbf i}\)

We we work on a bounded domain \(\Omega\) with no-flux conditions \(\mathbf F \cdot \mathbf n = 0\) on \(\partial \Omega\).

For simplicity the results on the slides are presented in \(d=1\).

\(W = 0\): Comparison principle

For \(\partial_t \rho = \operatorname{div} (\rho \nabla U'(\rho) + V)\) in \(\mathbb R^d\) or in a bounded domain with Dirichlet/no-flux conditions we have that \[ \int(\overline \rho_t - \underline \rho_t)_+ \le \int(\overline \rho_0 - \underline \rho_0)_+. \] This implies that ordered initial data imply ordered solutions, and hence uniqueness.

\(W\ne0\): no comparison principle

\(\mob(\rho) = \rho, U = 0, V = 0, W(x) = -\exp(-x^2/10)\)

Linear mobility and \(W = 0\): A monotone formulation

We rewrite the scheme \[ \DeclareMathOperator{\sign}{sign} \begin{aligned} H_i(\rho_{i-1}^{n+1}, \rho_i^{n+1}, \rho_{i+1}^{n+1}) &= \rho_i^n \end{aligned} \] with the nonlinear functions \[ \begin{aligned} H_i(\rho_{i-1}, \rho_i, \rho_{i+1}) &= \rho_i + \frac{F_{i+\frac 1 2}(\rho_{i}, \rho_{i+1}) - F_{i-\frac 1 2}(\rho_{i-1}, \rho_i)}{h} \end{aligned} \]

and the usual \[ \begin{aligned} F_{i+\frac 12}( a , b ) &= a (v_{i+\frac 12}(a,b))_+ + b (v_{i+\frac 12}(a,b))_- \\ v_{i+\frac 12}( a , b) & = -\frac{\xi_{i+1}(b) - \xi_i(a)}{h} \\ \xi_i (a) &= U'(a) + V(x_i). \end{aligned} \]

Linear mobility and \(W = 0\): Monotonicity of \(H\)

\[ \begin{aligned} H_i(\rho_{i-1}, \rho_i, \rho_{i+1}) &= \rho_i + \frac{F_{i+\frac 1 2}(\rho_{i}, \rho_{i+1}) - F_{i-\frac 1 2}(\rho_{i-1}, \rho_i)}{h} \\ F_{i+\frac 12}( a , b ) &= a (v_{i+\frac 12}(a,b))_+ + b (v_{i+\frac 12}(a,b))_- \\ v_{i+\frac 12}( a , b) & = -\frac{\xi_{i+1}(b) - \xi_i(a)}{h} \\ \xi_i (a) &= U'(a) + V(x_i). \end{aligned} \]

Notice that \[ \begin{gathered} \xi_i'(a) = U''(a) \ge 0 \qquad \frac{\partial v_{i+\frac 12}}{\partial a} = h^{-1} \xi_{i}'(a) \ge 0 \qquad \frac{\partial v_{i+\frac 12}}{\partial b} = -h^{-1} \xi_{i+1}'(b) \le 0 \end{gathered} \]

The up-winding for interior points is such that \[ \begin{aligned} \frac{\partial F_{i+\frac 12}}{\partial a} &= (v_{i+\frac 12})_+ + a \sign_+(v_{i+\frac 12}) \frac{\partial v_i}{\partial a} \ge 0 \\ \frac{\partial F_{i+\frac 12}}{\partial b} &= (v_{i+\frac 12})_- + a \sign_-(v_{i+\frac 12}) \frac{\partial v_{i+\frac 12}}{\partial b} \le 0 \end{aligned} \] and if \(i+\frac 1 2\) then \(F_{i+\frac 1 2} = 0\) so the sign of the directional derivatives is preserved.

Linear mobility and \(W = 0\): A monotone formulation

\[ \begin{aligned} H_i(\rho_{i-1}^{n+1}, \rho_i^{n+1}, \rho_{i+1}^{n+1}) &= \rho_i^n \end{aligned} \]

Eventually \[ \begin{aligned} \frac{\partial H_i}{\partial \rho_i} &= 1 + h^{-1} \frac{\partial F_{i+\frac 1 2}}{\partial a} - h^{-1} \frac{\partial F_{i-\frac 1 2}}{\partial b} \ge 1, \\ \frac{\partial H_i}{\partial \rho_{i-1}} &= - h^{-1} \frac{\partial F_{i-\frac 1 2}}{\partial a} \le 0 \\ \frac{\partial H_i}{\partial \rho_{i+1}} &= h^{-1} \frac{\partial F_{i+\frac 1 2}}{\partial b} \le 0 \end{aligned} \]

The up-winding gives this monotonocity

Linear mobility and \(W = 0\): \(L^1\)-comparison principle

Assume that

• \(\frac{\partial H_i}{\partial \rho_i} \ge 0\) and \(\frac{\partial H_i}{\partial \rho_j} \le 0\) for \(j \ne i\)
• The conservation of mass \(\sum_{i=1}^N H_i (\rho) = \sum_{i=1}^n \rho_i\)
• \(H(\overline \rho_i^{n+1}) \ge \overline \rho_i^{n+1}\) and \(H(\underline \rho^{n+1}) \le \underline \rho_i^{n}\).

Then, we have that
\[ \sum_i (\overline \rho_i^{n+1} - \underline \rho_i^{n+1})_+ \le \sum_i (\overline \rho_i^{n} - \underline \rho_i^{n})_+ \]

Proof.

Due to the conservation of mass we have \(\sum_i \frac{\partial H_i}{\partial s_j} = 1,\) so \(\frac{\partial H_j}{\partial s_j} \ge 0\).

Abusing slightly the notation, we write \(H_i (\rho) = H_i (\rho_i, (\rho_j)_{j\ne i})\).

With the monotonicity indicated \[\begin{equation} \begin{aligned} H_i \Big(\underline\rho_i, (\min\{ \underline\rho_j , \overline \rho_j \} )_{j\ne i}\Big) &\ge H_i (\underline\rho) \ge \underline f_i - \underline g_i, \\ H_i \Big(\overline \rho_i, (\min\{ \underline \rho_j , \overline \rho_j \} )_{j\ne i}\Big) &\ge H_i (\overline \rho) \ge \overline f_i - \overline g_i. \end{aligned} \end{equation}\] We have that \[\begin{equation*} H_i \Big((\min\{ \underline\rho_j , \overline \rho_j \} )_{j \in I}\Big) = \begin{dcases} H_i \Big(\underline\rho_i, (\min\{ \underline\rho_j , \overline \rho_j \} )_{j\ne i}\Big), & \text{if } \underline\rho_i \le \overline \rho_i , \\ H_i \Big(\overline \rho_i, (\min\{\underline \rho_j , \overline \rho_j \} )_{j\ne i}\Big), & \text{if } \overline\rho_i \le \underline\rho_i . \end{dcases} \end{equation*}\] In either case \(H_i (\min\{\underline\rho,\overline \rho\}) \ge \min\{\underline f_i, \overline f_i\} - \max \{\underline g_i, \overline g_i\}\).

Similarly, we get \(H_i (\max\{\underline \rho,\overline \rho\}) \le \max\{ \underline f_i, \overline f_i \} + \max \{\underline g_i, \overline g_i\}\).

Then, we write \[\begin{align*} \sum_{i \in I} &\left( \max\{ \underline\rho_i, \overline \rho_i \} - \min\{\underline \rho_i, \overline \rho_i \} \right) \\ & = \sum_{i \in I} \Big( H_i (\max\{\underline \rho,\overline \rho\}) - H_i (\min\{\underline \rho,\overline \rho\}) \Big) \\ & \le \sum_{i \in I} (\max\{\underline f_i, \overline f_i\} - \min\{\underline f_i, \overline f_i\} + 2 \max\{\underline g_i, \overline g_i\} ). \end{align*}\] Taking into account that \(|a-b| = \max\{a,b\} - \min \{a,b\}\) we recover \[ \sum_{i \in I} |\underline \rho_i - \overline \rho_i| \le \sum_{i \in I} |\underline f_i - \overline f_i| + 2\sum_{i \in I} \max\{\underline g_i, \overline g_i\} . \] Similarly to before, we also have due to mass conservation that \[ \sum_{i \in I} (\underline \rho_i - \overline \rho_i) = \sum_{i \in I} (H_i (\underline \rho) - H_i(\overline \rho_i)) \le \sum_{i \in I} (\underline f_i - \overline f_i) + \sum_{i \in I} (\underline g_i + \overline g_i) . \] Using that \(a^+ = \frac{|a| + a}{2}\) we deduce that \[ \sum_{i \in I} (\underline \rho_i - \overline \rho_i)^+ \]

Non-linear mobility and \(W = 0\): a monotone scheme by up-winding

We still have comparison principle, by taking \(P(\rho) = \mob(\rho) U''(\rho)\).

Consider the prototypical case \(\mob(\rho) = \rho^\alpha (1-\rho)^\beta\).

We have to pick an up-winding so that \[ \frac{\partial F}{\partial a} \ge 0, \qquad \frac{\partial F}{\partial b} \le 0. \]

In this case a choice of up-winding that works is \[ \begin{aligned} F_{i+\frac 1 2} (a,b) &= a^\alpha (1-b)^{\beta} (v_{i+\frac 1 2}(a,b))_+ \\ &\quad + b^{\alpha} (1-a)^{\alpha} (v_{i+\frac 1 2}(a,b))_- \end{aligned} \] where \(a = \rho_i\) and \(b = \rho_{i+1}\).

Up-winding general mobilities

\[ \DeclareMathOperator{\mobup}{\mob_\uparrow} \DeclareMathOperator{\mobdown}{\mob_\downarrow} \] We want to write \[ \mob(\rho) = \mobup(\rho) \mobdown (\rho) \] where \(\mobup\) is non-decreasing and \(\mobdown\) is non-increasing.

If we take \(\log\) we want \(\log \mob = \log \mobup + \log \mobdown\) where the first is non-decreasing and the second non-increasing.

We therefore take \[\begin{align*} \log \mobup(\rho) &= \frac{\log \mob (\rho_0)}{2} + \int_{\rho_0}^\rho \left( \frac{d}{d\rho} \log \mob \rho \right)_+, \\ \log \mobdown (\rho) &= \frac{\log \mob (\rho_0)}{2} + \int_{\rho_0}^\rho \left( \frac{d}{d\rho} \log \mob (\rho) \right)_- . \end{align*}\]

Scheme with general mobilities [Bailo, Carrillo & Hu, 2023]

\[ \tag{ADE} \begin{aligned} \partial_t \rho &= \operatorname{div} (\mob(\rho) \nabla (U'(\rho) + V + W*\rho)) \\ &\qquad \text{where } \mob(\rho) = \mobup(\rho) \mobdown(\rho). \end{aligned} \]

Continuous	Discrete
\(\partial_t \rho = - \operatorname{div} \mathbf F\)	\(\displaystyle \frac{\rho_{\mathbf i}^{n+1} - \rho_{\mathbf i}^n}{\Delta t} = -\sum_{k=1}^d \frac{F_{\mathbf i + \frac 1 2 \mathbf e_k}^{n+1}-F_{\mathbf i - \frac 1 2 \mathbf e_k}^{n+1}}{h_k}\)
\(\mathbf F = \mob(\rho) \mathbf v\)	\[\begin{aligned}F_{\mathbf i + \frac 1 2 \mathbf e_k}^{n+1} &= \mobup(\rho_{\mathbf i}^{n+1})\mobdown(\rho_{\mathbf i + \mathbf e_k}^{n+1}) (v_{\mathbf i + \frac 1 2 \mathbf e_k})_+ \\ &\qquad + \mobup(\rho_{\mathbf i + \mathbf e_k}^{n+1})\mobdown(\rho_{\mathbf i}^{n+1}) (v_{\mathbf i + \frac 1 2 \mathbf e_k})_- \end{aligned}\]
\(\mathbf v = -\nabla\xi\)	\(\displaystyle v_{\mathbf i + \frac 1 2 \mathbf e_k} = -\frac{\xi_{\mathbf i + \mathbf e_k}-\xi_{\mathbf i}}{h_k}\)
\(\xi = U'(\rho) + V + W*\rho\)	\(\xi_{\mathbf i}^{n+1} = U'(\rho_{\mathbf i}^{n+1}) + V(x_{\mathbf i}) + (W \star \rho)^{n+1}_{\mathbf i}\)

Convergence by compactness

A priori estimate: evolution of the \(\max/\min\) when \(W \ne 0\)

Let \(M(t) = \max \rho_t\). Formally \(M(t) = \rho_t(\xi(t))\) and
\[\begin{align*} \frac{d M(t)}{dt} &= \partial_t \rho(\xi(t)) \\ &= \mob'(\rho(\xi(t))) {\nabla \rho(\xi(t))} \cdot v + \mob(\rho(t)) \Delta v(\xi(t)) \\ &\le \mob(\rho_t(\xi(t))) \Delta ( V + (W*\rho))(\xi(t)). \end{align*}\] Since \(\mob(\alpha) = 0\), if \(M(0) \le \alpha\) then \(M(t) \le \alpha\) for all \(t>0\).

Similarly \(m(t) \ge 0\).

Justification: [Constantin & Escher, 2002].

A priori estimate: evolution of the \(\max/\min\) when \(W \ne 0\)

\(\newcommand{\Rho}{P} \newcommand{\bm}[1]{\mathbf{\mathrm{#1}}}\) Lemma. Let \(\bm i_0 \in \bm I_{\bm h} (\Omega)\) be such that \(\Rho^{n+1}_{\bm i_0} = \max_{\bm i } \Rho^{n+1}_{\bm i}\), then

• If \(\bm i_{0} + \bm e_k \in \bm I_{\bm h} (\Omega)\) then \[\begin{equation*} F_{\bm i_0 + \frac 1 2 \bm e_k}^{n+1} \ge \mob(\Rho^{n+1}_{\bm i_0} ) v_{\bm i_0 + \frac 1 2 e_k}^{n+1}. \end{equation*}\]
• If \(\bm i_{0} - \bm e_k \in \bm I_{\bm h} (\Omega)\) then \[\begin{equation*} F_{\bm i_0 - \frac 1 2 \bm e_k}^{n+1} \le \mob(\Rho^{n+1}_{\bm i_0} ) v_{\bm i_0- \frac 1 2 e_k}^{n+1}. \end{equation*}\]
• Therefore, if \(\bm i_{0} + \bm e_k , \bm i_{0} - \bm e_k \in \bm I_{\bm h} (\Omega)\) then \[\begin{equation} \label{eq:sup as limsup} \max_{\bm i } \Rho^{n+1}_{\bm i} \le \max_{\bm i } \Rho^{n}_{\bm i} - \tau \mob \left( \Rho_{\bm i_0}^{n+1} \right) \sum_{k=1}^d \frac{v^{n+1}_{\bm i_0 + \frac 1 2 \bm e_k} - v^{n+1}_{\bm i_0 - \frac 1 2 \bm e_k}}{h_k} \end{equation}\] The equivalent results holds for \(\min\).

\(\newcommand{\thesup}{{P_{\bm i_0}^{n+1}}} \newcommand{\mobupwind}{\mob_w}\)

Sketch of Proof

For convenience we write \(\mobupwind(a,b) = \mobup(a)\mobdown(b)\).

If \(\bm x_{\bm i} , x_{\bm i + \bm e_k} \in \Omega\) then, noticing that \(\mob(a) = \mobupwind(a,a)\) we can write \[\begin{align*} F_{\bm i + \frac 1 2 \bm e_k}^{n+1} & = \mob(\Rho_{\bm i}^{n+1}) v^{n+1}_{\bm i + \frac 1 2 \bm e_k} + \Big( \mobupwind (\Rho_{\bm i}^{n+1}, \Rho_{\bm i + \bm e_k}^{n+1})- \mobupwind(\Rho_{\bm i}^{n+1},\Rho_{\bm i}^{n+1}) \Big) (v^{n+1}_{\bm i + \frac 1 2 \bm e_k})_+ \\ & \quad + \Big( \mobupwind (\Rho_{\bm i + \bm e_k}^{n+1},\Rho_{\bm i}^{n+1}) - \mobupwind (\Rho_{\bm i}^{n+1}, \Rho_{\bm i}^{n+1}) \Big) (v^{n+1}_{\bm i + \frac 1 2 \bm e_k})_-. \end{align*}\]

Recalling that \(\mobupwind(a,b)\) is non-decreasing in \(a\) and non-increasing in \(b\), that \(a_+ \ge 0\) and \(a_- \le 0\) we have that \[\begin{align*} F_{\bm i + \frac 1 2 \bm e_k}^{n+1} & \ge \mob(\Rho_{\bm i}^{n+1}) v^{n+1}_{\bm i + \frac 1 2 \bm e_k} + \Big( \mobupwind (\Rho_{\bm i}^{n+1}, \thesup)- \mobupwind(\Rho_{\bm i_0}^{n+1},\Rho_{\bm i}^{n+1}) \Big) (v^{n+1}_{\bm i + \frac 1 2 \bm e_k})_+ \\ & \quad + \Big( \mobupwind (\thesup,\Rho_{\bm i}^{n+1}) - \mobupwind (\Rho_{\bm i}^{n+1}, \Rho_{\bm i_0}^{n+1}) \Big) (v^{n+1}_{\bm i + \frac 1 2 \bm e_k})_- \end{align*}\] Evaluating at \(i = i_0\) we prove the first claim. The second result can be proven analogously. Joining these two facts we get the last claim.

Discrete sub and super-solution

Lemma. Assume \(0 < \rho_0 < \alpha\), \(U \in C^2\) is convex, and \(V\) and replace \(W*\rho\) with \(\int_\Omega K(x,y) \rho(y) dy\) with \(K\) compactly supported. Let \[\begin{equation*} \lambda = \| D^2 V \|_{L^\infty} + \|D^2 K\|_{L^\infty} \| \rho_0 \|_{L^1}. \end{equation*}\]

Then, letting \(\underline \Rho^n = \inf_{\bm i} \Rho_{\bm i}^n\) we have \[\begin{equation*} \underline \Rho^{n+1} + 2\tau \lambda d \mob(\underline \Rho^{n+1}) \ge \underline \Rho^n. \end{equation*}\]

Corollary. If \(\mob'\) is Lipschitz and \(\Rho^n \ge 0\) we can uniformly bound \(\Rho^{n+1} \ge (1 + c\tau)^{-1} \min \Rho_i^n\). This leads to exponential lower bounds.

Remark. If we remove the condition compact support condition, the same estimate holds without rates.

Remark. When \(\mob(\alpha) = 0\) the equivalent result holds for \(\alpha - \max_i \Rho_i^n\).

A priori estimate: decay of the free energy

At the continuous, if \(W(x) = W(-x)\) simple integration by parts ensures that \[\begin{equation*} \frac{d}{dt} \underbrace{\int \left(U(\rho_t) + \rho_t V + \frac 1 2 \rho_t (W*\rho_t)\right)}_{\mathcal F[\rho_t]} = - \int \mob(\rho) |\underbrace{\nabla (U'(\rho) + V + W*\rho)}_{-v}|^2. \end{equation*}\]

Sketch of proof and assumptions.

We use the convexity inequality \(U(b) - U(a) \le U'(b) (b-a)\) and

• If \(W\) is positive-definite then \((W\star\Rho)_i^n = \sum_{j=1}^n W(x_i - x_j) \Rho_j^{n+1}\).
• If \(W\) is negative-definite then \((W\star\Rho)_i^n = \sum_{j=1}^n W(x_i - x_j) \Rho_j^{n}\).
• In a general setting \((W\star\Rho)_i^n = \sum_{j=1}^n W(x_i - x_j)\frac{\Rho_j^{n+1} + \Rho_{j}^n}{2}\).

A priori estimate: evolution of other energies

\[\begin{align*} \frac{d}{dt} \int \Lambda(\rho_t) = - \int \mob(\rho_t) \Lambda''(\rho_t) \nabla \rho_t \cdot (U'(\rho) + V + W*\rho) = - \int \nabla H'(\rho) \cdot v. \end{align*}\] if \(\Lambda''(s) \mob(s) = H''(s)\).

At the discrete level we can an adapted “chain rule” to prove \[\begin{equation} \label{eq:energy control with general H} h^d \sum_{\bm i \in \bm I} \Lambda(\Rho_{\bm i}^{n+1}) - \tau h^d \sum_{k=1}^d \sum_{\substack{\bm i \in \bm I \text{ s.t.} \\ \bm i + \bm e_k \in \bm I}} \frac{ H'(\Rho_{\bm i + \bm e_k}^{n+1}) - H'(\Rho_{\bm i}^{n+1})}{h_k}v^{n+1}_{\bm i + \frac 1 2 \bm e_{k}} \le h^d \sum_{\bm i \in \bm I} \Lambda(\Rho_{\bm i}^{n}). \end{equation}\]

Compactness theorem: Riesz-Fréchet-Kolmogorov

A sequence \(\rho_n\) in \(L^1(\mathbb R^d)\) is pre-compact if and only if

• It has controlled shifts \[ \int |\rho_n(x+y) - \rho_n(x)| dx \le \omega(|y|). \]

• It has controlled tails \[ \int_{\mathbb R^d \setminus B_R} |\rho_n| \le \omega(R^{-1}). \]

Compactness theorem: Aubin-Lions theorem

If there exist spaces such that \(X_0 \overset{\text{compact}}\subset X \overset{\text{cont.}}\subset Y\) and we have uniform bounds \[\begin{equation*} u_n \in L^p (0,T; X_0) , \qquad \partial_t u_n \in L^p (0,T; Y) \end{equation*}\] then \(u_n\) is pre-compact in \(L^p(0,T; X)\).

Typically \(X_0 = W^{1,p}\), \(X = L^p\), and \(Y = W^{-1,p}\)

Notice that we have that \[ \partial_t \rho_t = \operatorname{div} (\mob(\rho) \underbrace{\nabla (U'(\rho) + V + W*\rho)}_{-v}) \]

and we estimate \[\begin{align*} \int_t^s \int \left|\mob(\rho) v\right|^2 &\le \|\mob(\rho)\|_{L^\infty} \int_t^s \int \mob(\rho) \left| v \right|^2 \le \|\mob(\rho)\|_{L^\infty}(\mathcal F[\rho_t] - \mathcal F[\rho_s]) \end{align*}\]

Compactness theorem: Discrete Aubin-Simon theorem [Gallouët & Latché, 2012]

\(\newcommand{\seqi}{\mathsf k}\) Let \((B, \|\cdot \|_B)\) be a Banach space, \(B_\seqi\) a sequence of finite-dimensional subspaces of \(B\), and \(\| \cdot \|_{X_\seqi}, \| \cdot \|_{Y_\seqi}: B_\seqi \to \mathbb R\) be norms of \(B_\seqi\) such that:

• If \(u_\seqi \in B_\seqi\) is such that \(\|u_\seqi\|_{X_\seqi}\) are bounded, then a subsequence converges strongly in \(B\).
• If \(u_\seqi \in B_\seqi\) for each \(\seqi\), \(u_\seqi \to u\) strongly in \(B\), and \(\|u_\seqi\|_{Y_\seqi} \to 0\), then \(u = 0\). \end{enumerate}

Let \(H_\seqi \subset L^\infty(0,T; B_\seqi)\) be the space of functions which are constant in time for \(t \in (n \tau_\seqi, (n+1)\tau_\seqi)\). By notation, for \(u \in H_\seqi\) we define \(u^n = u((n+\frac 1 2)\tau_\seqi)\).
Fix \(T > 0\) and assume that we have a sequence \(u_\seqi \in H_\seqi\) such that, for some \(q \in [1,\infty)\) and all \(\seqi \in \mathbb N\) \[\begin{equation*} \tau_n \sum_{n=0}^{ \lfloor T / \tau_n \rfloor + 1 } \| u_\seqi^n \|_{X_\seqi}^q + \tau_n \sum_{n=0}^{ \lfloor T / \tau_n \rfloor } \left\| \frac{ u_\seqi^{n+1} - u_\seqi^n}{\tau_\seqi} \right\|^q_{Y_\seqi} \le C. \end{equation*}\] Then \(u_\seqi\) has a subsequence that converges in \(L^q(0,T; B)\).

Convergence results

[Bailo, Carrillo, Murakawa & Schmidtchen, 2020]:

• \(\mob(\rho) = \rho\)
• \(0 < c_1 \le \rho_t \le c_2\)
• Free boundary does not form for some \(T_0 > 0\)
• Convergence up to \(T_0\)

[Gómez-Castro, 2025] Setting 1 (no free-boundary):

• There exists constants \(C_i\) such that \(0 < C_1 \le \rho_0 \le C_2 < \alpha\).
• \(\mobup, \mobdown \in W^{1,\infty}(0,\alpha)\).
• \(U \in C^2((0,\alpha))\) with \(U'' > 0\) in \((0,\alpha)\).
• \(V \in W^{2,\infty} (\Omega)\), \(K \in W^{2,\infty} (\Omega \times \Omega)\), and \(K(x,y) = K(y,x)\) with \(\nabla V = 0\) and \(\nabla_x K = 0\) on \(\partial \Omega\).
• Free boundary will not form for any \(T > 0\).
• Convergence for any \(T > 0\).

More convergence results

[Gómez-Castro, 2025] Setting 2 (free-boundary)

• \(0 \le \rho_0 \le \alpha\).
• \(\mobup, \mobdown \in C ([0,\alpha])\).
• \(U \in C^2((0,\alpha)) \cap C^1([0,\alpha])\) with \(\inf_{s \in (0,\alpha)} U'' (s) > 0\), and \(\Lambda_U \in C([0,\alpha])\).
• \(V \in W^{1,\infty} (\Omega)\), \(K \in W^{1,\infty} (\Omega \times \Omega)\), and \(K(x,y) = K(y,x)\).
• Convergence for any \(T > 0\).

Simulations

\(d=1\) convergence towards a stady state when \(V(x) = x^2/2\)

\(d=2\)

Upcoming work

The Porous-Medium Equation with non-local pressure

\[ \partial_t \rho_t = \operatorname{div}(\rho^{m-1} \nabla \xi) , \qquad (-\Delta)^s \xi = \rho. \]

\(m = 2\): [Caffarelli & Vázquez, 2011]

\(m \ne 2\): [Stan, Del Teso & Vázquez, 2016],[Stan, Del Teso & Vázquez, 2019]

Radial solutions \(s = 1\): [Carrillo, Gómez-Castro & Vázquez, 2022a], [Carrillo, Gómez-Castro & Vázquez, 2022b]

Numerics \(d=1\): [del Teso & Jakobsen, 2023] studying \(v(t,x) = \int_{-\infty}^x \rho(t,y) dy\)

We can use the techniques above for the case \(d > 1\).

Bibliography

Bailo, R., Carrillo, J.A. & Gómez-Castro, D. (2026) Aggregation-Diffusion Equations for Collective Behaviour in the Sciences. In: S. Oishi, H. Okamoto, & K. Hayami (eds.). Recent Developments in Industrial and Applied Mathematics. Singapore, Springer Nature Singapore. pp. 177–200. doi:10.1007/978-981-95-1446-5_9.

Bailo, R., Carrillo, J.A. & Hu, J. (2023) Bound-preserving finite-volume schemes for systems of continuity equations with saturation. SIAM J. Num. Anal. 83 (2). doi:10.1137/22M1488703.

Bailo, R., Carrillo, J.A. & Hu, J. (2020) Fully discrete positivity-preserving and energy-dissipating schemes for aggregation-diffusion equations with a gradient-flow structure. Communications in Mathematical Sciences. 18 (5), 1259–1303. doi:10.4310/cms.2020.v18.n5.a5.

Bailo, R., Carrillo, J.A., Murakawa, H. & Schmidtchen, M. (2020) Convergence of a fully discrete and energy-dissipating finite-volume scheme for aggregation-diffusion equations. Mathematical Models and Methods in Applied Sciences. 30 (13), 2487–2522. doi:10.1142/s0218202520500487.

Caffarelli, L. & Vázquez, J.-L. (2011) Asymptotic behaviour of a porous medium equation with fractional diffusion. Discrete & Continuous Dynamical Systems - A. 29 (4), 1393–1404. doi:10.3934/dcds.2011.29.1393.

Carrillo, J.A., Fernández-Jiménez, A. & Gómez-Castro, D. (2025) Drift-diffusion equations with saturation. SIAM J. of Math. Anal. 57 (6), 6300–6360. doi:10.1137/24M1706694.

Carrillo, J.A., Gómez-Castro, D. & Vázquez, J.L. (2022a) A fast regularisation of a Newtonian vortex equation. Annales de l’Institut Henri Poincaré C, Analyse non linéaire. http://arxiv.org/abs/1912.00912.

Carrillo, J.A., Gómez-Castro, D. & Vázquez, J.L. (2022b) Vortex formation for a non-local interaction model with Newtonian repulsion and superlinear mobility. Advances in Nonlinear Analysis. 11 (1), 937–967. doi:10.1515/anona-2021-0231.

Constantin, A. & Escher, J. (2002) Global solutions for quasilinear parabolic problems. Journal of Evolution Equations. 2 (1), 97–111. doi:10.1007/s00028-002-8081-2.

Deimling, K. (1985) Nonlinear Functional Analysis. Berlin, Heidelberg, Springer Berlin Heidelberg. doi:10.1007/978-3-662-00547-7.

Eymard, R., Gallouët, T. & Herbin, R. (1997) Finite volume methods. In: Handbook of numerical analysis.

Gallouët, T. & Latché, J.-C. (2012) Compactness of discrete approximate solutions to parabolic PDEs - Application to a turbulence model. Communications on Pure and Applied Analysis. 11 (6), 2371–2391. doi:10.3934/cpaa.2012.11.2371.

Gómez-Castro, D. (2024) Beginner’s guide to aggregation-diffusion equations. SeMA Journal. doi:10.1007/s40324-024-00350-y.

Gómez-Castro, D. (2025) Convergence of a finite-volume scheme for aggregation-diffusion equations with saturation. 2025. https://arxiv.org/abs/2507.11132.

Stan, D., Del Teso, F. & Vázquez, J.L. (2019) Existence of Weak Solutions for a General Porous Medium Equation with Nonlocal Pressure. Archive for Rational Mechanics and Analysis. 233 (1), 451–496. doi:10.1007/s00205-019-01361-0.

Stan, D., Del Teso, F. & Vázquez, J.L. (2016) Finite and infinite speed of propagation for porous medium equations with nonlocal pressure. Journal of Differential Equations. 260 (2), 1154–1199. doi:10.1016/j.jde.2015.09.023.

del Teso, F. & Jakobsen, E.R. (2023) A convergent finite difference-quadrature scheme for the porous medium equation with nonlocal pressure. 9 March 2023. http://arxiv.org/abs/2303.05168 [Accessed: 29 May 2023].

Appendix: Existence by topological degree

Topological degree: definition

\(\DeclareMathOperator{\degree}{deg}\)

We recall the notion of topological degree in \(\mathbb R^{|I|}\), see e.g., [Deimling, 1985]. The topological degree is a function \(\degree :X \to \mathbb Z\) where \[ \begin{aligned} X = \{ (f, D, y) : & D \subset \mathbb R^d \text{ is open and bounded, } \\ & f: \overline D \to \mathbb R^d \text{ continuous } ,\\ & y \in \mathbb R^d \setminus f(\partial D) \}. \end{aligned} \]

The construction is by extension of this definition of the case \(f \in C^1 (\overline \Omega)\), \(y \not \in f(\partial \Omega)\) a regular value¹ of \(f\):

\[ \degree(f,\Omega,y) = \sum_{x \in f^{-1}(y)} \operatorname{sign} \det Df(x). \] by convention \(0\) if \(f^{-1}(y)\) is empty.

Topological degree: key properties

We use three key properties.

• \(\degree(\mathrm{id},D, y) = 1\).

• if \(h: [0,1] \times \overline D \to \mathbb R^d\) is continuous and \(y \notin h(\lambda, \partial D)\) then \[ \degree( h(\lambda, \cdot), D, y ) \text{ is constant in } \lambda. \]

• Theorem 3.1 in [Deimling, 1985]
  If \(\degree(f, D, y) \ne 0\), then there exists \(x \in D\) such that \(f(x) = y\).

Non-linear mobility and \(W=0\):
Existence of solutions by topological degree

Lemma. Assume that

• \(U \in C^1([0,\alpha])\),
• \(\mob(0) = \mob(\alpha) = 0\), and
• \(0 \le \rho^n \le \alpha\).

Then, there exists \(0 \le \rho^{n+1} \le \alpha\).

Sketch of proof.

We write an extended formulation
\[ \widetilde H_i(\lambda, \rho) = \rho_i + \lambda \Delta t \frac{\widetilde F_{i+\frac{1} 2 }(\rho) - \widetilde F_{i-\frac{1} 2 }(\rho)}{\Delta x} \qquad \text{ for } \lambda \in [0,1] \text{ and } \rho \in \mathbb R^{|I|}. \]

For some \(R > 0\) to be determined \(D_R =\left\{ \rho \in \mathbb R^{|I|} : \sum_{i\in I}|\rho_i| < R \right\}.\)

Consider \(\rho^{n} = \widetilde H(\lambda, \rho)\) with \(\rho \in \partial D_R\) for some \(\lambda \in [0,1]\).

By the comparison principle, since \(0 \le \rho_i^n \le \alpha\) then \(0\le \rho_i^n \le \alpha\).

Since the fluxes at zero at the boundary \(\sum_i |\rho_i| = \sum_i \rho_i = \sum_i \widetilde H_i(\lambda, \rho) = \sum_i \rho^n_i.\)

Hence, for \(R > \sum_{i \in I} \rho_i^n\) it is clear that \(\rho_i^n \notin \widetilde H(\lambda, \partial D_R)\) for any \(\lambda > 0\).

\(\degree ( \widetilde H(\lambda, \cdot) , D_R, \rho^n ) = \degree ( \widetilde H(0,\cdot), D_R, \rho^n ) = \degree(\mathrm{id}, D_R, \rho^n ) = 1 \ne 0. \square\)

Appendix: \(L^1\) comparison principle

\(L^1\) comparison priciple

When \(W = 0\) the problem is local, and we expect a comparison principle.

Notice that for any two solutions, let \(P'(s) = s U''(s)\) and use \(\varphi = \Phi'(P(\overline \rho) - P(\underline \rho))\) as a test function \[ \begin{aligned} & \int \partial_t (\overline \rho - \underline \rho) \Phi'(P(\overline \rho) - P(\underline \rho)) \\ & { =- \int \Phi''\Big(P(\overline \rho) - P(\underline \rho)\Big) \nabla \Big(P(\overline \rho) - P(\underline \rho)\Big) \cdot \Big( \overline \rho \nabla (U'(\overline \rho) + V ) - \underline \rho \nabla (U'(\underline \rho) + V )\Big) } \\ &{ = - \int \Phi''\Big(P(\overline \rho) - P(\underline \rho)\Big) \nabla \Big(P(\overline \rho) - P(\underline \rho)\Big) \cdot \nabla \Big(P(\overline \rho) - P(\underline \rho)\Big) } \\ &{ \qquad - \int \Phi''\Big(P(\overline \rho) - P(\underline \rho)\Big) \nabla \Big(P(\overline \rho) - P(\underline \rho)\Big) \cdot (\overline \rho - \underline \rho) \nabla V }\\ &{ \le - \int (\overline \rho - \underline \rho) \Phi''\Big(P(\overline \rho) - P(\underline \rho)\Big) \nabla \Big(P(\overline \rho) - P(\underline \rho)\Big)\cdot \nabla V. } \end{aligned} \]

We deduce that \[ \begin{aligned} & \int \partial_t (\overline \rho - \underline \rho) \Phi'(P(\overline \rho) - P(\underline \rho)) \\ & \le - \int (\overline \rho - \underline \rho) \Phi''\Big(P(\overline \rho) - P(\underline \rho)\Big) \nabla \Big(P(\overline \rho) - P(\underline \rho)\Big)\cdot \nabla V. \end{aligned} \]

If we take \(\Phi''(s) = \varepsilon^{-1} \chi_{[0,\varepsilon]}\) so that \(\Phi'(s) \to \sign_+(s)\), we have (under suitable bounds on \(\rho\)) that
\[ \left|(\overline \rho - \underline \rho) \Phi''\Big(P(\overline \rho) - P(\underline \rho)\Big) \right| \le \chi_{ \{ 0 \le \overline \rho - \underline \rho \le \varepsilon \} } . \]

As \(\varepsilon \to 0\) we recover \(\partial_t \int (\overline \rho_t - \underline \rho_t)_+ \le 0.\)

This is the well-known \(L^1\) comparison principle \(\int (\overline \rho_t - \underline \rho_t)_+ \le \int (\overline \rho_0 - \underline \rho_0)_+ .\)

It implies, among other things, uniqueness of solutions.

---

1. for all \(x \in f^{-1}(y)\) we have \(\det Df(x) \ne 0\))