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*}\]

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*}\]

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: 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}\]

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*}\]

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\).

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

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\).

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\).

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)