# Discretization and convergence of the ballistic Benamou-Brenier formulation of the porous medium and Burgers' equations

Jean-Marie Mirebeau, Erwan Stampfli

► To cite this version:

Jean-Marie Mirebeau, Erwan Stampfli. Discretization and convergence of the ballistic Benamou-Brenier formulation of the porous medium and Burgers' equations. 2025. hal-05005367v2

HAL Id: hal-05005367

https://hal.science/hal-05005367v2

Preprint submitted on 17 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discretization and convergence of the ballistic Benamou-Brenier formulation of the porous medium and Burgers' equations.

Jean-Marie Mirebeau,\* Erwan Stampfli†

April 6, 2025

## Abstract

We study the discretization, convergence, and numerical implementation of recent reformulations of the quadratic porous medium equation (multidimensional and anisotropic) and Burgers' equation (one-dimensional, with optional viscosity), as forward in time variants of the Benamou-Brenier formulation of optimal transport. This approach turns those evolution problems into global optimization problems in time and space, of which we introduce a discretization, one of whose originalities lies in the harmonic interpolation of the densities involved. We prove that the resulting schemes are unconditionally stable w.r.t. the space and time steps, and we establish a quadratic convergence rate for the dual PDE solution, under suitable assumptions. We also show that the schemes can be efficiently solved numerically using a proximal splitting method and a global space-time fast Fourier transform, and we illustrate our results with numerical experiments.

## 1 Introduction

A number of non-linear partial differential equations (PDEs) can be reformulated as global convex optimization problems in time and space [Bre20], revealing so-called *hidden convexity* properties. We follow in this paper an approach whose properties are particularly appealing for evolution problems featuring quadratic non-linearities [Vor22], already investigated numerically for Burgers' equation [Bre20, KA24] and for non-linear elasticity [SGA24]. In the continuous setting, its benefits include new notions of weak solutions, and related methods for proving their existence and/or uniqueness. In the discrete setting, it leads to a radically new method for solving evolution problems, of which this paper provides the first numerical convergence analysis, with schemes free of any Courant-Friedrichs-Levy (CFL) condition, second-order accurate in time and space, are solved using proximal optimization algorithms and global space-time fast-Fourier transforms. As discussed below, this approach shares similarities with the Benamou-Brenier formulation of optimal transport, first discretized in [PPO14] and whose numerical analysis is studied in e.g. [NT21] or [LZ24] in the matrix-valued setting. While [Vor22, Bre20] mostly focus on *systems* of PDEs, often related to fluid mechanics, we shall limit in this paper our attention to two *scalar* PDEs for simplicity, see (1). Before that, let us further motivate the approach and replace it within the literature.

<footnote>\*University Paris-Saclay, ENS Paris-Saclay, CNRS, Centre Borelli, F-91190 Gif-sur-Yvette, France.</footnote>

<footnote>†University Paris-Saclay, Laboratoires de Mathématiques d'Orsay, 91400 Orsay, France.</footnote>

<footnote n="0">⁰Mathematics Subject Classification (2020) : 65M06, 65M12, 49M29, 35A15, 35L65</footnote>
The initial value problem for many non-linear evolution PDEs may admit infinitely many weak solutions for most initial conditions [DiP79, LJ09, LS10]. For the Euler equations of incompressible fluids and, also, for the class of hyperbolic systems of conservation laws with a convex entropy, Brenier [Bre18] considered the following problem: given an initial condition $u_0$ and a finite time interval $[0, T]$, find among all weak solutions starting from $u_0$ one that minimizes the time integral of the kinetic energy (for Euler) or the entropy (for systems of conservation laws). This problem might have no solution, but it generates a dual variational problem which is automatically convex and usually admits a solution, that we call dual solution. Brenier obtained an explicit formula relating this dual solution to any tentative smooth solution of the initial value problem, at least when $T$ is not too large, see theorem 1.8. It was also pointed out that dual problems of this kind are strongly related to optimal transport problems (in their ‘Benamou-Brenier’ formulation [BB09, VC03]) and can be interpreted as generalized (with matrix-valued densities) Mean Field Games à la Lasry-Lions. The method of [Bre18] was further developed by Vorotnikov [Vor22] for a very large class of quadratic evolution PDEs and a similar duality approach has been used by Acharya and collaborators [Ach23, AS24, ASZ24]. Let us finally mention [BM22] where the dual technique is applied to a multi-stream formulation of the Euler-Poisson system and [Bre20] where various parabolic PDEs are considered, including the QPME, for which there is no restriction on the size of $T$ and which is the main subject of the present paper.

As announced, this paper is devoted to the discretization and convergence analysis of dual reformulations of two scalar evolution PDEs featuring quadratic non-linearities: the quadratic porous medium equation (QPME, multidimensional and anisotropic), and Burgers' equation with optional viscosity (one dimensional).

$$ \text{QPME}: \partial_t u = \frac{1}{2} \text{div}(\mathcal{D} \nabla u^2), \qquad \text{Burgers': } \partial_t u + \frac{1}{2} \partial_x u^2 = \nu \partial_{xx} u. \qquad (1) $$

They are posed for simplicity on the periodic unit box $\mathbb{T}^d$, with $d = 1$ for Burgers' equation, where $\mathbb{T} := \mathbb{R}/\mathbb{Z}$ is the periodic unit interval. See [KA24] for a closely related approach implementing non-periodic boundary conditions in Burgers' case. The PDE parameters are a smooth positive definite tensor field $\mathcal{D} \in C^\infty(\mathbb{T}^d, S_d^{++})$, following [Li18], and a non-negative diffusion coefficient $\nu \ge 0$ respectively. A smooth and positive initial condition $u_0 \in C^\infty(\mathbb{T}^d, ]0, \infty[)$ and a bounded time interval $[0, T]$ are also given. We make these strong assumptions for the sake of the numerical analysis, even though the PDEs of interest admit solutions under weaker assumptions [Vá07]. Our numerical experiments section 4 also show that the proposed schemes are robust and appear to work without these assumptions.

The two PDEs (1) share a common structure, with a quadratic non-linearity:

$$ \partial_t u + \frac{1}{2} L(u^2) + LAu = 0 \quad \text{on } [0, T], \qquad u(0) = u_0, \qquad (2) $$

where $L = -\text{div}(\mathcal{D}\nabla\cdot)$ and $A = 0$ for the QPME, and $L = \partial_x$ and $A = -\nu\partial_x$ for Burgers' equation. Following [Bre20, Vor22] we formally consider the (apparently trivial) problem (3) of minimizing the $L^2$ norm among all weak solutions of (2); variants featuring a time-dependent weight or a base state can also be considered, see [Vor25, ASZ24] and theorem 1.8. Formally, again, we exchange in (4) the sup and inf,

$$ \inf_{u(0)=u_0} \sup_{\phi(T)=0} \int_{[0,T]\times\mathbb{T}^d} \frac{1}{2} u^2 + (\partial_t u + \frac{1}{2} L(u^2) + LAu)\phi \qquad (3) $$

$$ \stackrel{?}{=} \sup_{\phi(T)=0} \inf_{u(0)=u_0} \int_{[0,T]\times\mathbb{T}^d} \frac{1}{2} u^2 (1 + L^*\phi) - (u - u_0)\partial_t \phi + u A^* L^* \phi. \qquad (4) $$
Here and below, given a bivariate mapping $u(t, x)$, we denote by $u(t) := u(t, \cdot)$ the univariate function of the space variable alone, where the time argument has been frozen to some given value $t$. Note that $(3) \ge (4)$ always, and the difference $(3) - (4)$ is referred to as the duality gap. We make at this point the critical assumption that $1 + L^*\phi > 0$ everywhere on $[0, T] \times \mathbb{T}^d$, which is closely related with the lack of duality gap, see theorem 1.8. Minimizing pointwise (4) w.r.t. the primal variable $u$, we obtain the non-trivial dual problem [Bre20, Vor22, KA24]:

$$ \inf_{\phi(T)=0} \int_{[0,T] \times \mathbb{T}^d} \frac{(\partial_t \phi - A^* L^* \phi)^2}{2(1 + L^* \phi)} - u_0 \partial_t \phi. \quad (5) $$

The primal and dual variables $u$ and $\phi$ obey the following optimality relation

$$ \partial_t \phi - A^* L^* \phi = (1 + L^* \phi)u, \qquad \phi(T) = 0, \qquad (6) $$

which can be regarded as a reversed in time linear evolution problem w.r.t. $\phi$. (For large $T$, this may contradict our assumption $1 + L^*\phi > 0$, see theorem 1.8.) Letting

$$ m := \partial_t \phi, \qquad \rho := 1 + L^* \phi, \qquad (7) $$

we obtain as announced a reformulation of the evolution PDE (2) as an optimization problem which is global in time and space, and is close in spirit to the Benamou-Brenier formulation of optimal transport:

$$ \inf_{m, \rho} \int_{[0,T] \times \mathbb{T}^d} \frac{(m - A^* \rho)^2}{2\rho} - u_0 m \quad \text{subject to } \partial_t \rho = L^* m \text{ and } \rho(T) = 1. \quad (8) $$

The solution of the original problem (1) can be recovered as $u = (m - A^*\rho)/\rho$ at all points where $\rho > 0$, in view of (6) and (7), in other words

$$ \text{QPME: } u = m/\rho, \qquad \text{Burgers': } u = (m - \nu \partial_x \rho)/\rho. $$

We refer to [Bre20] and [Vor25, Appendix B] for further discussion of the analogy with optimal transport, and the ‘ballistic’ terminology. See also [HT+19] for the connection of the QPME with martingale optimal transport.

As announced, this paper is devoted to a numerical convergence analysis of the energetic formulation (8) of the QPME and Burgers’ equation (1). For simplicity, we discretize the unknowns on Cartesian grids. For stability and accuracy, we use staggered finite differences in time and space, and for that purpose we define $\mathcal{T}_\tau, \mathcal{T}_\tau' \subset [0, T]$ and $\mathbb{T}_h, \mathbb{T}_h' \subset \mathbb{T}$ as

$$ \begin{aligned} \mathcal{T}_\tau &:= \{0, 2\tau, \dots, T\}, & \mathcal{T}_\tau' &:= \{\tau, 3\tau, \dots, T - \tau\}, \\ \mathbb{T}_h &:= \{0, 2h, \dots, 1 - 2h\}, & \mathbb{T}_h' &:= \{h, 3h, \dots, 1 - h\}. \end{aligned} $$

where $\tau > 0$ and $h > 0$ respectively denote the *half* timestep, and *half* gridscale. The grids $\mathcal{T}_\tau$ and $\mathbb{T}_h$ are *centered*, whereas $\mathcal{T}_\tau'$ and $\mathbb{T}_h'$ are *staggered*. The following assumption is implicit throughout the paper.

**Assumption 1.1.** The grid sizes $N_\tau := T/(2\tau)$ and $N_h := 1/(2h)$ are positive integers.

The staggered first order time finite difference, and the first order spatial finite difference in the direction $e \in \mathbb{Z}^d$ (which may or may not be staggered depending on the parity of the components of $e$), are defined as

$$ \partial_{\tau} u(t, x) := \frac{u(t + \tau, x) - u(t - \tau, x)}{2\tau}, \qquad \partial_{h}^{e} u(t, x) := \frac{u(t, x + he) - u(t, x - he)}{2h}. \quad (9) $$
The one-dimensional spatial finite difference operators are denoted $\partial_h := \partial_h^{(1)}$ and $\partial_{hh} := \partial_h \partial_h$.
We also use the averaged $\ell^1(\mathbb{T}_h^d)$ norm, and the perspective function $\mathcal{P}: \mathbb{R}^2 \to [0, \infty]$ [CM18]

$$ \|f\|_{\ell^1(\mathbb{T}_h^d)} := (2h)^d \sum_{x \in \mathbb{T}_h^d} |f(x)|, \qquad \mathcal{P}(m, \rho) := \begin{cases} \frac{m^2}{2\rho} & \text{if } \rho > 0, \\ 0 & \text{if } m = \rho = 0, \\ \infty & \text{else.} \end{cases} \qquad (10) $$

The ratios appearing in the energies (13) and (14) should be understood as instances of $\mathcal{P}$, which is convex and lower semi-continuous in view of the identity $\mathcal{P}(m, \rho) = \sup\{a\rho+bm \mid a+\frac{1}{2}b^2 \le 0\}$.

**Discretization of the QPME.** A preliminary step is the construction of a finite differences approximation of the anisotropic divergence-form Laplacian operator $-L := \text{div}(D\nabla\cdot)$. In the isotropic case where $\mathcal{D} = \text{Id}$ identically, the usual finite differences Laplacian $-L_h := \sum_{1\le i\le d} (\partial_h^{e_i})^2$ should be used, where $(e_1, \dots, e_d)$ denotes the canonical basis, in expanded form $-L_h u(x) := (2h)^{-2} \sum_{1\le i\le d} [u(x+2he_i) - 2u(x) + u(x-2he_i)]$. In the anisotropic case, we use an adaptive finite differences scheme with similar structural properties: for any $u: \mathbb{T}_h^d \to \mathbb{R}$

$$ -L_h u := \sum_{e \in E} \partial_h^e (\lambda^e \partial_h^e u), \qquad \text{where } \mathcal{D} = \sum_{e \in E} \lambda^e e e^\top. \qquad (11) $$

We denoted by $E \subset \mathbb{Z}^d$ a finite set referred to as the stencil, and by $\lambda^e \in C^\infty(\mathbb{T}^d, [0, \infty])$ some smooth non-negative coefficients obeying (11, right), see e.g. [FM14, K$^+$05, BBM23] for the construction and the convergence analysis of numerical schemes based on such positive definite matrix decompositions. Expanding (11, left) we obtain for all $x \in \mathbb{T}_h^d$ (note that $x+2he \in \mathbb{T}_h^d$ for any $e \in E$)

$$ -L_h u(x) = (2h)^{-2} \sum_{e \in E} \left[ \lambda^e(x+he)(u(x+2he)-u(x)) + \lambda^e(x-he)(u(x-2he)-u(x)) \right]. \quad (12) $$

**Lemma 1.2.** The set $\Phi_{\tau h}^P := \{\phi \in \mathbb{R}^{T_\tau \times \mathbb{T}_h^d} \mid \phi(T) = 0\}$ is in bijection with $\{(m, \rho) \in \mathbb{R}^{T_\tau \times \mathbb{T}_h^d} \times \mathbb{R}^{T_\tau \times \mathbb{T}_h^d} \mid \partial_\tau \rho = L_h m, \rho(T) = 1\}$, via $\phi \mapsto (\partial_\tau \phi, 1+L_h \phi)$.

We omit the proof, which is a direct consequence of the time-independence of the coefficients $\lambda^e$, for all $e \in E$, and of the commutation of the time and space finite differences (9). Using these notations, we discretize the dual energy (8, left) as follows

$$ \mathcal{E}_{\tau h}^{P}(m, \rho) := 2\tau(2h)^d \sum_{\substack{t \in T_{\tau} \\ x \in \mathbb{T}_h^d}} \left( \frac{1}{2} \sum_{\sigma=\pm} \frac{m(t, x)^2}{2\rho(t + \sigma\tau, x)} - m(t, x)u_0(x) \right). \qquad (13) $$

The second sum symbol is features two contributions, corresponding to the two possible signs $\sigma = +$ and $\sigma = -$. By convention we set $\mathcal{E}_{\tau h}^P(\phi) := \mathcal{E}_{\tau h}^P(\partial_\tau \phi, 1+L_h\phi)$, and we show in theorem 2.3 that this energy is a strictly convex function of $\phi \in \Phi_{\tau h}^P$. Let us emphasize that the energy (13) implicitly involves the harmonic mean of the density field $\rho$, at the times $t+\tau$ and $t-\tau$ associated with the two possible signs $\sigma = \pm$. This choice is deliberate, and we show in section B that using the arithmetic mean instead as in [PPO14, NT21, LZ24] leads to lesser stability properties, and to a degeneracy of the optimization problem, which becomes local in time rather than global (49).
**Assumption 1.3.** We assume that $u \in C^\infty([0, T] \times \mathbb{T}^d, ]0, \infty[)$ is a smooth positive solution of the QPME (1, left), associated with a smooth and symmetric positive definite diffusion tensor $D \in C^\infty(\mathbb{T}^d, S_d^{++})$. We denote by $\phi \in C^\infty([0, T] \times \mathbb{T}^d)$ the solution of the reversed in time parabolic PDE: $\partial_t\phi = (1 - \text{div}(D\nabla\phi))u$ with terminal boundary condition: $\phi(T) = 0$, and we assume that $1 - \text{div}(D\nabla\phi)$ is positively bounded below. The numerical scheme stencil $E \subset \mathbb{Z}^d \setminus \{0\}$ is finite, and the coefficients $\lambda^e \in C^\infty(\mathbb{T}^d, ]0, \infty[)$ are smooth for all $e \in E$.

**Theorem 1.4.** Under theorem 1.3. For all $0 < \tau \le \tau_*$ and $0 < h \le h_*$, there exists a unique minimizer $\tilde{\phi} \in \Phi_{\tau h}^P$ of the energy $\mathcal{E}_{\tau h}^P$. It satisfies

$$ \max_{t \in T_{\tau}} \|\tilde{\phi}(t) - \phi(t)\|_{\ell^1(\mathbb{T}_h^d)} \le C_* (\tau^2 + h^2). $$

The constants $\tau_*, h_* > 0$ and $C_*$ only depend on $(\lambda, \phi)$.

A uniformly converging approximation of the primal PDE solution $u = m/\rho = \partial_t\phi/(1 - \text{div}(D\nabla\phi))$ can also be recovered from $\tilde{\phi}$, using a regularization by convolution to estimate its derivatives, see section C. In the numerical experiments section 4, a second order convergence rate is empirically observed for the primal variable $u$, without the need for such a convolution.

**Discretization of Burgers' equation.** This equation is built around a first order, rather than second order, differentiation operator in space, and for this reason the problem density $\rho$ is discretized on an appropriately staggered grid.

**Lemma 1.5.** The set $\Phi_{\tau h}^B := \{\phi \in \mathbb{R}^{T_\tau \times T_h} | \phi(T) = 0\}$ is in bijection with $\{(m, \rho) \in \mathbb{R}^{T'_\tau \times T_h} \times \mathbb{R}^{T_\tau \times T'_h} | \partial_\tau \rho + \partial_h m = 0, \rho(T) = 1\}$, via $\phi \mapsto (\partial_\tau \phi, 1 - \partial_h \phi)$.

We again omit the proof, which directly follows from the commutation of the finite difference operators $\partial_\tau$ and $\partial_h$. Using these notations we discretize (8, left) as follows

$$ \mathcal{E}_{\tau h}^{B}(m, \rho) := 4\tau h \sum_{\substack{t \in T'_{\tau} \\ x \in T_h}} \left( \frac{1}{4} \sum_{\substack{\sigma_t = \pm \\ \sigma_x = \pm}} \frac{(m(t, x) - \nu \partial_h \rho(t + \sigma_t \tau, x))^2}{2\rho(t + \sigma_t \tau, x + \sigma_x h)} - m(t, x)u_0(x) \right). \quad (14) $$

By convention we set $\mathcal{E}_{\tau h}^B(\phi) := \mathcal{E}_{\tau h}^B(\partial_\tau \phi, 1 - \partial_h \phi)$, and we show in theorem A.2 that this energy is a strictly convex function of $\phi \in \Phi_{\tau h}^B$. The energy (14) implicitly involves the harmonic mean of four values of the density field $\rho$, at all combinations of the times $t \pm \sigma_t \tau$ and positions $x \pm \sigma_x h$, and this is again the key to the proof of stability and convergence. In the discretization of the Benamou-Brenier formulation of optimal transport, a closely related energy is considered, with $\nu = 0$ and $u_0 = 0$, and a Cauchy rather than ballistic problem in time. Somewhat curiously, the arithmetic mean is usually chosen in this context [PPO14, NT21, LZ24] for the interpolation of the density w.r.t. time. Let us mention that [NT21] observes improved stability and less oscillations in the solution when using harmonic (as opposed to arithmetic) interpolation in space, yet only arithmetic interpolation is considered in time. The viscous Burgers' equation is also related with Schrödinger's bridge problem and entropic optimal transport [Bre20, §5.2].

**Assumption 1.6.** We assume that $u \in C^\infty([0, T] \times \mathbb{T}, ]0, \infty[)$ is a smooth positive solution of Burgers' equation (1, right), with non-negative viscosity $\nu \ge 0$ (possibly zero). We denote by $\phi \in C^\infty([0, T] \times \mathbb{T})$ the solution of the reversed in time first order linear PDE: $\partial_t\phi + \nu\partial_{xx}\phi = (1 - \partial_x\phi)u$ with b.c.: $\phi(T) = 0$, and we assume that $1 - \partial_x\phi$ is positively bounded below.
The assumption of a smooth solution is not too restrictive in the case of the QPME, see theorem 1.3, which has a regularizing effect on positive solutions [Vá07], and likewise for Burgers' equation with positive viscosity. In contrast, it appears highly questionable for the inviscid Burgers' equation ($\nu = 0$), which develops a shock in finite time $T_{\text{shock}} > 0$, even if the initial condition is smooth. However, strikingly and unfortunately, it turns out that the BBB approach on the interval $[0, T]$ does not characterize the (viscosity) solution $u$ to Burgers' equation if $T > T_{\text{shock}}$, but rather a modified solution $u^T$ depending on the final time $T$. Let us emphasize that, nevertheless, the correct entropy solution is still obtained at the final time $t = T$, although it is not for $t < T$. In other words $u^T(T, x) = u(T, x)$, yet in general $u^T(t, x) \neq u(t, x)$ for $t \in [0, T[$, and in particular $u^T(0, x) \neq u_0(x)$, see [Bre20, Theorem 5.2.2] or [Bre18], as well as fig. 5 20 for a numerical illustration. (The critical assumption $1 - \partial_x \phi > 0$ is also violated if $T > T_{\text{shock}}$.) A numerical method for Burgers' equation proposed in [KA24] addresses this issue by concatenating the numerical solutions of the BBB formulation discretized over small subintervals of time, see also theorem 1.8. Proving convergence in this setting is however outside the scope of this paper, and we thus stick to theorem 1.6, which implicitly implies $T < T_{\text{shock}}$ in the inviscid case.

**Theorem 1.7.** Under theorem 1.6. For all $0 < \tau \le \tau_*$ and $0 < h \le h_*$, there exists a unique minimizer $\tilde{\phi} \in \Phi_{\tau h}^B$ of the energy $\mathcal{E}_{\tau h}^B$. It satisfies

$$ \max_{t \in T_{\tau}} \|\tilde{\phi}(t) - \phi(t)\|_{\ell^1(\mathbb{T}_h)} \le C_{*}(\tau^2 + h^2). $$

The constants $\tau_*, h_* > 0$ and $C_*$ only depend on $\phi$.

**Minimization algorithm and implementation.** We numerically minimize the energy (13) (resp. (14)), w.r.t. the momentum $m$ and density $\rho$ variables subject to the continuity equation discretized as in theorem 1.2 (resp. theorem 1.5). For that purpose, we use the classical Chambolle-Pock [CP11] proximal splitting primal-dual optimization algorithm, similarly to [PPO14], which requires a reformulation of the objective functional using carefully chosen additional variables. (A Newton solver can also be used for small problem sizes, following [KA24].)

For instance we introduce an anti-symmetric (resp. symmetric) extension of the optimization unknown $m$ (resp. $\rho$) to negative times, which has two benefits. First, it decouples the perspective functions appearing in (13), so that the corresponding proximal operator is computable independently for each point in time and space. Second, it allows using a global FFT in time and space, in order to compute the orthogonal projection onto the linear subspace associated with the discretized continuity equation, see theorems 1.2 and 1.5. In a similar spirit, additional unknowns $m_e = \partial_h^e m$, $e \in E$, containing the directional derivatives of the momentum in the stencil directions, are required for the anisotropic QPME.

The proposed numerical method, being embarrassingly parallel in space and time, is easily ported to GPUs. Experiments confirm second order accuracy, and illustrate the use of large time steps in contrast with e.g. the explicit scheme, which is subject to a parabolic type Courant-Friedrichs-Levy (CFL) stability condition, or the semi-implicit scheme. See sections 3 and 4.

**Contributions and comparisons with related work.** The ballistic Benamou-Brenier (BBB) formulation of non-linear PDEs is a subject of active research in the continuous setting [Bre20, Vor22, Vor25, Ach23, AS24, ASZ24]. The discretization of this approach has so far only been considered for one-dimensional problems and without proof of convergence under mesh refinement: for Burgers' equation [Bre20, KA24], and for non-linear elasticity [SGA24]. This paper
extends the numerical study to the QPME which is a multidimensional and possibly anisotropic problem, in addition to Burgers' equation, and establishes second order convergence rates when the solution is smooth, see theorems 1.4 and 1.7.

Our results can also be related to the discretization [PPO14] and convergence analysis [NT21, LZ24] of the Benamou-Brenier formulation of optimal transport [BB09]. In contrast with some of these works, we assume the existence of a smooth solution of the PDE, which is justified in view of the regularizing effect of the QPME on positive solutions (numerical experiments section 4 nevertheless show that the scheme also works on non-smooth compactly supported solutions), and likewise for the viscous Burgers' equation (the inviscid case being subject to an obstruction discussed below theorem 1.6). Our numerical scheme involves an harmonic (rather than arithmetic) local interpolation of the density w.r.t. time, see section B, in contrast with earlier works. The implementation of the proximal optimization method is modified accordingly, and also to address the ballistic time-boundary conditions, and the possibly anisotropic spatial finite difference operators.

The discretization of the anisotropy (12), via adaptive decompositions of the anisotropic diffusion tensor field, is in the vein of [FM14, BBM22, BBM23].

**Summary.** Our main convergence result theorem 1.4 for the QPME is established in section 2, the proximal optimization method for the energies (13) and (14) is described in section 3, and numerical experiments are presented in section 4. Section A is devoted to the convergence result theorem 1.7 for Burgers' equation, section B discusses a modification of the energy (13) involving the arithmetic interpolation of the density field, and section C deduces from theorems 1.4 and 1.7 and a mollification argument a uniform convergence result for the solution $u$ of the non-linear PDE of interest (1). Finally, an explicit solution to the BBB formulation of the QPME, corresponding to the Barenblatt profile, is presented in section D, and a well-posedness result for the BBB formulation of Burgers' viscous equation is given in section E.

**Remark 1.8** (Positivity of the density and final evolution time). The BBB formulation requires that the dual potential $\phi$ obeying the time-reversed evolution PDE (6) satisfies (at the very least) $1 + L^*\phi \ge 0$ almost everywhere in time and space, see (5). Depending on the PDE of interest, this condition is typically satisfied for small $T$ in view of the terminal condition $\Phi(T, \cdot) = 0$, but may fail for larger time intervals. See [Bre20, Vor22] for more discussion on this topic, which is directly related with the absence of duality gap in (4).

In the case of the QPME, a non-trivial argument involving Aronson-Bénilan estimates [Bre20, Theorem 5.1.1] fortunately shows that an arbitrarily large time $T$ is admissible. Likewise for Burgers' equation with positive viscosity, see [Bre20, Proposition 5.2.1] and section E. However, in the case of the inviscid Burgers' equation, the largest admissible time $T$ is typically finite and corresponds to the time of formation of the first shock, as discussed below theorem 1.6.

Two approaches aimed at extending the time $T$ of applicability of the BBB formulation, for a selection of PDEs, have been proposed. Vorotnikov [Vor25, Remark 4.2] minimizes $\int_{[0,T] \times \mathbb{T}^d} u(t,x)^2 \exp(-\gamma t) dxdt$ among all weak solutions of the PDE of interest, where $\gamma > 0$ is a sufficiently large discount factor, rather than the unweighted kinetic energy (3) corresponding to $\gamma = 0$. Acharya et al [ASZ24] minimize $\int_{[0,T] \times \mathbb{T}^d} (u - \bar{u})^2$, where $\bar{u}$ is a well chosen base state, which can be regarded as an initial guess for the solution of the primal problem; in numerical applications they also replace the full interval $[0, T]$ with a succession of overlapping sub-intervals [KA24].
## 2 Convergence analysis, quadratic porous medium equation

This section is devoted to the proof of our main result: the unconditional second order convergence rate theorem 1.4 for the proposed discretization of the quadratic porous medium equation. We begin by recalling an exact expansion satisfied by the perspective function [CM18].

**Proposition 2.1.** Let $m, \hat{m}, \rho, \hat{\rho} \in \mathbb{R}$ be such that $P(m, \rho) < \infty$. Then

$$P(m + \hat{m}, \rho + \hat{\rho}) = P(m, \rho) + (u\hat{m} - \frac{1}{2}u^2\hat{\rho}) + P(\hat{m} - \hat{\rho}u, \rho + \hat{\rho}), \quad \text{with } u := \begin{cases} m/\rho & \text{if } \rho > 0, \\ 0 & \text{if } m = \rho = 0. \end{cases}$$

We omit the proof, which follows from direct substitution. The perspective function is convex and lower-semi-continuous, as already observed below (10), and the vector $(u, -\frac{1}{2}u^2)$ is its gradient if $\rho > 0$, and a subgradient element if $\rho = m = 0$.

**Notation 2.2.** We fix in the following the half timestep $\tau > 0$ and half gridscale $h > 0$. For readability, and in order to avoid multiple sub- or super-scripts, we denote by $\mathcal{E} := \mathcal{E}_{\tau h}^P$ the energy of interest and by $\Phi := \Phi_{\tau h}^P$ its domain, see theorem 1.2. Likewise, we denote by $\mathcal{T} := \mathcal{T}_{\tau}$ and $\mathcal{T}' := \mathcal{T}'_{\tau}$ the centered and staggered grids in time, and by $\mathcal{X} := \mathbb{T}_h^d$ the centered grid in space. Finally, we denote $\| \cdot \|_1 := \| \cdot \|_{\ell^1(\mathcal{X})}$ the normalized $\ell^1$ norm, see (10).

Using theorem 2.1 we obtain

$$\mathcal{E}(\phi + \hat{\phi}) = \mathcal{E}(\phi) + \mathcal{L}(\hat{\phi}; \phi) + \mathcal{Q}(\hat{\phi}; \phi), \quad (15)$$

for any $\phi, \hat{\phi} \in \Phi$ such that $\mathcal{E}(\phi) < \infty$. The linear and second order (non-quadratic) terms are

$$\mathcal{L}(\hat{\phi}; \phi) := \tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}'} \sum_{\sigma=\pm} \left[ u^{\sigma}(t, x) \hat{m}(t, x) - \frac{1}{2} u^{\sigma}(t, x)^2 \hat{\rho}(t + \sigma \tau, x) \right] + (2h)^d \sum_{x \in \mathcal{X}} u_0(x) \hat{\phi}(0, x), \\ \mathcal{Q}(\hat{\phi}; \phi) := \tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}'} \sum_{\sigma=\pm} \mathcal{P}(\hat{m}(t, x) - \hat{\rho}(t + \sigma \tau, x) u^{\sigma}(t, x), \rho(t + \sigma \tau, x) + \hat{\rho}(t + \sigma \tau, x)),$$

$$\text{where } m := \partial_{\tau} \phi, \rho := 1 + L_h \phi, \hat{m} := \partial_{\tau} \hat{\phi}, \hat{\rho} := L_h \hat{\phi}, \text{ and } u^{\sigma}(t, x) := \frac{m(t, x)}{\rho(t + \sigma \tau, x)}. \quad (16)$$

Note that $\phi, \hat{\phi}, \rho, \hat{\rho}$ are defined on $\mathcal{T} \times \mathcal{X}$, whereas $m, \hat{m}, u^{\sigma}$ are defined on $\mathcal{T}' \times \mathcal{X}$ for any sign $\sigma = \pm$. The denominator $\rho$ of $u^{\sigma}$ is positive for sufficiently small $h > 0$, by theorems 1.3 and 2.8 below. We used the telescopic sum $2\tau \sum_{t \in \mathcal{T}'} m(t, x) = -\phi(0, x)$ to obtain $\mathcal{L}$. The proof of theorem 1.4 is based on a quantitative coercivity property of the energy functional $\mathcal{E}$, improving on the following qualitative preliminary result.

**Proposition 2.3.** The energy $\mathcal{E}$ is strictly convex on the domain $\Phi$.

*Proof.* This energy is convex since the perspective function (10) is convex. Consider $\phi, \hat{\phi} \in \Phi$ such that $s \in [0, 1] \mapsto \mathcal{E}(\phi + s\hat{\phi})$ is finitely valued and affine. We prove by induction on $n = 0, \dots, N_{\tau} := T/(2\tau)$ that $\hat{\phi}(T - 2\tau n, x) = 0$ for all $x \in \mathcal{X}$. This property holds by definition of the domain $\Phi$ for $n = 0$. Assuming $\hat{\phi}(T - 2\tau n, x) = 0$, for all $x \in \mathcal{X}$, and denoting $t := T - 2\tau n - \tau$ we obtain $\hat{m}(t, x) = \hat{\phi}(t - \tau, x)/(2\tau)$ and $\hat{\rho}(t + \tau, x) = 0$ with the notation (16). Since $\mathcal{Q}(\hat{\phi}; \phi)$ is a sum of convex terms w.r.t. $\hat{\phi}$, and $\mathcal{L}(\hat{\phi}; \phi)$ is linear w.r.t. $\hat{\phi}$, the following term

$$\mathcal{P}(s\hat{m}(t, x) - s\hat{\rho}(t + \tau, x)u^{+}(t, x), \rho(t + \tau, x) + s\hat{\rho}(t + \tau, x)) = s^2\mathcal{P}(\hat{\phi}(t - \tau, x)/(2\tau), \rho(t + \tau, x))$$

(corresponding to $\sigma = +$) must be an affine function of $s \in [0, 1]$. This implies $\hat{\phi}(t - \tau, x) = 0$, which proves the induction. Thus $\hat{\phi} = 0$, and therefore $\mathcal{E}$ is strictly convex as announced. $\square$
## 2.1 Upper estimate of the linear term

We estimate in this section the linear term $\mathcal{L}(\hat{\phi}; \phi)$ in terms of the $\ell^1(\mathcal{X})$ norm of the perturbation $\hat{\phi}$, see the estimate (24) which holds under theorem 1.3. For that purpose, we first rewrite $\mathcal{L}(\hat{\phi}; \phi)$ in a suitable form using two discrete summation by parts identities, in time and space.

**Lemma 2.4.** For any $f: \mathcal{T} \to \mathbb{R}$ and $g: \mathcal{T}' \to \mathbb{R}$, one has denoting $\mathcal{T} := \mathcal{T} \setminus \{0, T\}$

$$2\tau \sum_{t \in \mathcal{T}'} g(t) \partial_{\tau} f(t) = f(T)g(T - \tau) - f(0)g(\tau) - 2\tau \sum_{t \in \mathcal{T}} f(t) \partial_{\tau} g(t). \quad (17)$$

The proof of this basic and classical identity is omitted, and we turn to the spatial counterpart (19), which incidentally shows that $L_h$ is positive semi-definite. This property is used in [BBM23, Appendix C] to establish convergence rates for the discretization of an elliptic equation, but not here since we instead rely on degenerate ellipticity, see theorem 2.9.

**Lemma 2.5.** The operator $\partial_h^e$ is skew-adjoint, for all $e \in E$, in the sense that

$$\sum_{x \in \mathcal{X}} u(x) \partial_h^e w(x) = \frac{1}{2\tau} \sum_{x \in \mathcal{X}} u(x) w(x + he) - \frac{1}{2\tau} \sum_{x \in \mathcal{X}} u(x) w(x - he) = - \sum_{y \in \mathcal{X} + he} w(y) \partial_h^e u(y) \quad (18)$$

for all $u \in \mathbb{R}^{\mathcal{X}}$, $w \in \mathbb{R}^{\mathcal{X}+he}$. The operator $L_h$ is self-adjoint: for all $u, v \in \mathbb{R}^{\mathcal{X}}$

$$\sum_{x \in \mathcal{X}} u(x) L_h v(x) = \sum_{e \in E} \sum_{y \in \mathcal{X} + he} \lambda^e(y) \partial_h^e u(y) \partial_h^e v(y) = \sum_{x \in \mathcal{X}} v(x) L_h u(x). \quad (19)$$

*Proof.* Equation (18) is obtained using the change of variables $y = x + he$ (resp. $y = x - he$) in the intermediate sums terms, which maps $\mathcal{X} \to \mathcal{X} + he$ (resp. likewise, since $\mathcal{X} - he = \mathcal{X} + he$ for any $e \in E \subset \mathbb{Z}^d$). Equation (19) follows from the defining formula $-L_h u := \sum_{e \in E} \partial_h^e (\lambda^e \partial_h^e u)$. $\square$

As announced, we rewrite $\mathcal{L}(\hat{\phi}; \phi)$ using successively theorems 2.4 and 2.5, and recalling (16):

$$\begin{align}
-\mathcal{L}(\hat{\phi}; \phi) &= (2h)^d \sum_{x \in \mathcal{X}} \left( \frac{u^{-}(\tau, x) + u^{+}(\tau, x)}{2} \hat{\phi}(0, x) + \frac{\tau}{2} u^{-}(\tau, x)^2 \hat{\rho}(0, x) - u_0(x) \hat{\phi}(0, x) \right) \nonumber \\
&\quad + \tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T} \atop \sigma=\pm} \left( \partial_{\tau} u^{\sigma}(t, x) \hat{\phi}(t, x) + \frac{1}{2} u^{\sigma}(t - \sigma\tau, x)^2 \hat{\rho}(t, x) \right) \nonumber \\
&= (2h)^d \sum_{x \in \mathcal{X}} \left( \frac{u^{-}(\tau, x) + u^{+}(\tau, x)}{2} + \frac{\tau}{2} L_h [(u^{-})^2](\tau, x) - u_0(x) \right) \hat{\phi}(0, x) \tag{20} \\
&\quad + \tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T} \atop \sigma=\pm} \left( \partial_{\tau} u^{\sigma}(t, x) + \frac{1}{2} L_h [(u^{\sigma})^2](t - \sigma\tau, x) \right) \hat{\phi}(t, x) \tag{21}
\end{align}$$

In the rest of this subsection, we work under theorem 1.3, and thus assume that $\phi \in C^{\infty}([0, T] \times \mathbb{T}^d, \mathbb{R})$ is a smooth function. We distinguish the variables

$$m = \partial_{\tau} \phi, \quad \rho := 1 + L_h \phi, \quad u^{\sigma}(t, x) := \frac{m(t, x)}{\rho(t + \sigma \tau, x)}, \quad m := \partial_{\tau} \phi, \quad \rho := 1 - \operatorname{div}(\mathcal{D} \nabla \phi), \quad u := \frac{m}{\rho},$$

and note that $m$ and $u^\sigma$ are defined on the time-restricted domain $[\tau, T - \tau] \times \mathbb{T}^d$, whereas the other quantities are defined over $[0, T] \times \mathbb{T}^d$. By theorem 1.3, one has $\rho > 0$ uniformly, and $u$ obeys the QPME: $\partial_t u = \frac{1}{2} \operatorname{div}(\mathcal{D} \nabla u^2)$ and $u(0, \cdot) = u_0$.

The following two basic lemmas, on the integral representation of the finite difference operator, and the Taylor expansion of a smooth function which is symmetric w.r.t. some variables, should be clear by themselves and are thus stated without proof.
**Lemma 2.6.** Let $F \in C^k(\mathbb{T}^d \times [-1, 1]^n)$, where $k \ge 1$, and let $e \in \mathbb{R}^d$. Define $G \in C^{k-1}(\mathbb{T}^d \times [-1, 1]^{n+1})$ as follows: for all $x \in \mathbb{T}^d$ and all $h_1, \dots, h_n, h \in [-1, 1]$

$$G(x; h_1, \dots, h_n, h) := \frac{1}{2} \int_{-1}^{1} e^{\top} \nabla F(x + she; h_1, \dots, h_n) ds,$$

where the gradient is w.r.t. the spatial variable $x$. Then $G(x; h_1, \dots, h_n, h) = \partial_h^e F(x; h_1, \dots, h_n)$ for all $h \ne 0$. In addition $g(x) = e^{\top} \nabla f(x)$, where $f(x) := F(x; 0, \dots)$ and $g(x) := G(x; 0, \dots)$.

**Lemma 2.7.** Let $F \in C^2(\mathbb{T}^d \times [-1, 1]^2)$. Define $a = 1$ (resp. $a = 2$ if $F(x; \tau, h) = F(x; -\tau, h)$) for all $x \in \mathbb{T}^d$ and $\tau, h \in [-1, 1]$, and likewise $b = 1$ (resp. $b = 2$ if $F(x; \tau, h) = F(x; \tau, -h)$) for all $x, \tau, h$. Then $\max_{x \in \mathbb{T}^d} |F(x; \tau, h) - f(x)| = O(|\tau|^a + |h|^b)$, where $f(x) := F(x; 0, 0)$.

**Proposition 2.8.** Under theorem 1.3 and with the notations $m, \rho, u^\sigma, m, \rho, u$ above one has

$$m = m + \mathcal{O}(\tau^2), \quad \rho = \rho + \mathcal{O}(h^2), \quad u^\sigma = u + \mathcal{O}(\tau + h^2), \quad u^+ + u^- = 2u + \mathcal{O}(\tau^2 + h^2), \quad (22)$$

$$-L_h u^- = \operatorname{div}(\mathcal{D}\nabla u) + \mathcal{O}(\tau + h^2), \quad -L_h [(u^-)^2] = \operatorname{div}(\mathcal{D}\nabla u^2) + \mathcal{O}(\tau + h^2), \quad (23)$$

$$\partial_{\tau} \sum_{\sigma=\pm} u^{\sigma} = 2\partial_{t} u + \mathcal{O}(\tau^2 + h^2), \quad \sum_{\sigma=\pm} L_h [(u^{\sigma})^2](t - \sigma\tau, x) = -2 \operatorname{div}(\mathcal{D}\nabla u^2)(t, x) + \mathcal{O}(\tau^2 + h^2),$$

where the $\mathcal{O}$ notation is understood in the $L^\infty$ norm over the common domain of definition.

*Proof.* It is sufficient in this proof to assume that $\phi \in C^6$, and the coefficients $\lambda^e \in C^5$ for all $e \in E$, rather than $C^\infty$, and thus theorem 1.3 may be slightly weakened. We consider an arbitrary extension of $\phi$ into a periodic function of the time variable $t \in \mathbb{R}/(T+1)\mathbb{Z}$, in addition to the space variable $x \in \mathbb{T}^d$, with the same smoothness, so as to fit the assumptions of theorems 2.6 and 2.7. We remove here theorem 1.1 on the timestep and gridscale, and let $\tau$ and $h$ be possibly zero, positive or negative. The announced estimates can be reformulated as $F(t, x; \tau, h) = f(t, x) + \mathcal{O}(|\tau|^a + |h|^b)$, for some exponents $a, b \in \{1, 2\}$. Expressing successive finite differences as integral operators as in theorem 2.6, we find that $F(t, x; \tau, h) = \mathfrak{F}(t, x; \tau, \dots, \tau, h, \dots, h)$ where $\mathfrak{F}(t, x; \tau_1, \dots, \tau_m, h_1, \dots, h_n)$ is at least $C^2$, and that $f(t, x) = \mathfrak{F}(t, x; 0, \dots, 0)$ holds by design. Theorem 2.7 thus applies, and the announced estimates follow from the symmetry properties of $F$ w.r.t. the arguments $\tau$ and $h$.

The staggered time finite difference $\partial_\tau$ defined in (9) is indeed symmetric w.r.t. $\tau$, thus $\partial_\tau \phi = \partial_t \phi + \mathcal{O}(\tau^2)$ which is (22, i). The scheme $L_h$ expanded in (12) is symmetric w.r.t. $h$, hence is the second order accurate: recalling the decomposition (11, right) of $\mathcal{D}$ we obtain

$$\operatorname{div}(\mathcal{D}\nabla \phi) = \sum_{e \in E} \operatorname{div}(\lambda^e e e^\top \nabla \phi) = \sum_{e \in E} e^\top \nabla (\lambda^e e^\top \nabla \phi) = -L_h \phi + \mathcal{O}(h^2)$$

which is (22, ii), see also [BBM23, Proposition C.15] for a similar discussion. In particular, $\rho$ is like $\rho$ positively bounded below for sufficiently small $h$, and thus $u^\sigma$ is well defined. The time shift in the denominator of $u^\sigma$ breaks the symmetry w.r.t. $\tau$, but not w.r.t. $h$, hence the $\mathcal{O}(\tau + h^2)$ error in (22, iii) and likewise (23). The sum $u^+ + u^-$ restores the symmetry w.r.t. $\tau$, hence (22, iv) and the remaining estimates. $\square$
Inserting the expansions of theorem 2.8 in (21) yields as desired

$$
\begin{align}
-\mathcal{L}(\hat{\phi}; \phi) &= (2h)^d \sum_{x \in \mathcal{X}} \left( \mathbf{u}(\tau, x) - \frac{\tau}{2} \operatorname{div}(\mathcal{D}\nabla \mathbf{u}^2)(\tau, x) - u_0(x) + \mathcal{O}(\tau^2 + h^2) \right) \hat{\phi}(0, x) \nonumber \\
&\quad + 2\tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \bar{\mathcal{T}}} \left( \partial_t \mathbf{u}(t, x) - \frac{1}{2} \operatorname{div}(\mathcal{D}\nabla \mathbf{u}^2)(t, x) + \mathcal{O}(\tau^2 + h^2) \right) \hat{\phi}(t, x) \nonumber \\
&\le C_{\mathcal{L}}(\tau^2 + h^2) \left( \|\hat{\phi}(0)\|_1 + \tau \sum_{t \in \bar{\mathcal{T}}} \|\hat{\phi}(t)\|_1 \right), \tag{24}
\end{align}
$$

where the constant $C_{\mathcal{L}}$ only depends on $(\lambda, \phi)$.

## 2.2 Lower estimate of the second order term

We conclude in this subsection the proof of our convergence result theorem 1.4, using the following two steps. We first estimate (27) the solution of an implicit discretization scheme for a time-reversed PDE (26). We then deduce from this estimate a lower bound for the second order term $Q(\hat{\phi}; \phi)$ in (15), and obtain from this coercivity-like property that the minimizer of the discretized energy and of the continuous energy functional are close to each other (30).

The term $Q(\hat{\phi}; \phi)$ defined in (15) is a sum of perspective functions, whose first argument reads

$$r^{\sigma}(t) := \hat{m}(t) - \hat{\rho}(t + \sigma\tau)u^{\sigma}(t) = \frac{\hat{\phi}(t + \tau) - \hat{\phi}(t - \tau)}{2\tau} - u^{\sigma}(t)L_{h}\hat{\phi}(t + \sigma\tau), \quad (25)$$

for all times $t \in \mathcal{T}'$ and signs $\sigma = \pm$. Here and below, given a bivariate function $u(t, x)$, we denote by $u(t) := u(t, \cdot)$ the univariate function where the time argument has been frozen to some given value $t$. The term (25) can be regarded as the r.h.s. of a timestep for the discretization of a time-reversed linear PDE closely related to (6), see also (26) below, with the terminal boundary condition $\hat{\phi}(T) = 0$. This timestep is explicit (resp. implicit) if $\sigma = +$ (resp. $\sigma = -$).

**Estimate of $\hat{\phi}$ in terms of $r^{-}$.** We begin with a stability lemma, based on the discrete degenerate ellipticity of the adaptive finite difference scheme $L_h$ for the anisotropic Laplacian $-\operatorname{div}(\mathcal{D}\nabla \cdot)$, see [BBM23]. We say that a linear operator on $\mathbb{R}^{\mathcal{X}}$ is inverse positive, if it is invertible and the inverse matrix has non-negative entries in the canonical basis.

**Proposition 2.9.** For any $u \in [0, \infty]^{\mathcal{X}}$ the operator $\operatorname{Id} + uL_h$ on $\mathbb{R}^{\mathcal{X}}$ is inverse positive. If in addition $\|L_h u\|_{\infty} < 1$, then $\|(\operatorname{Id} + uL_h)^{-1} f\|_1 \le (1 - \|L_h u\|_{\infty})^{-1} \|f\|_1$, for all $f \in \mathbb{R}^{\mathcal{X}}$.

*Proof.* The matrix of the operator $L_h$ has non-positive off-diagonal entries, and its rows (or columns since it is symmetric) sum to zero, as is clear from the expanded form (12). Therefore, and since $u$ is non-negative, the matrix of $\operatorname{Id} + uL_h$ has non-positive off-diagonal entries, and its rows sum to one. This fits the conditions of an M-matrix, thus an inverse positive matrix, as announced.

Now consider $f \ge 0$, and let $g := (\operatorname{Id} + uL_h)^{-1} f$, which is thus non-negative. One has

$$\frac{\|f\|_1}{(2h)^d} = \sum_{x \in \mathcal{X}} f(x) = \sum_{x \in \mathcal{X}} [g(x) + u(x)L_h g(x)] = \sum_{x \in \mathcal{X}} g(x)(1 + L_h u(x)) \ge \frac{\|g\|_1}{(2h)^d} (1 - \|L_h u\|_{\infty}),$$

using the self-adjointness of $L_h$, see theorem 2.5. Observing that any $f \in \mathbb{R}^{\mathcal{X}}$ can be expressed as the difference $f = f_+ - f_-$ of its positive and negative part, with $\|f\|_1 = \|f_+\|_1 + \|f_-\|_1$, we conclude the proof by linearity of $(\operatorname{Id} + uL_h)^{-1}$. $\square$
In the following we choose the sign $\sigma = -$ corresponding to the *implicit* scheme in (25), for stability considerations in view of theorem 2.9, and thus for readability we denote $u := u^-$ and $r := r^-$. In addition to the explicit and implicit schemes, we discuss in section B the semi-implicit scheme, which arises when considering a slightly modified energy functional. Rearranging terms in (25) we obtain

$$ (1 + 2\tau u(t)L_h)\hat{\phi}(t - \tau) = \hat{\phi}(t + \tau) - 2\tau r(t), \qquad (26) $$

where the spatial variable $x \in \mathcal{X}$ is omitted for readability.

In the rest of this subsection, we work under theorem 1.3, in such way that $\phi \in C^\infty([0, T] \times \mathbb{T}^d, \mathbb{R})$ is a smooth function, and the QPME solution $u > 0$ is smooth and positive.

By theorem 2.8 one has $u = u + O(\tau + h^2)$. Therefore $u \ge 0$ on $[\tau, T - \tau] \times \mathbb{T}^d$ for any half timestep $0 < \tau \le \tau_*$ and half gridscale $0 < h \le h_*$, where $\tau_*, h_* > 0$ depend only on $(\lambda, \phi)$. Using (26) and theorem 2.9 we obtain

$$ \|\hat{\phi}(T - 2\tau N)\|_1 \le (1 - \tau K)^{-1} \left( \|\phi(T - 2\tau (N - 1))\|_1 + 2\tau \|r(T - (2N - 1)\tau)\|_1 \right) $$

for any $0 < N \le N_\tau := T/(2\tau)$, and provided $\tau K < 1$. Thus by an immediate induction

$$ \|\hat{\phi}(T - 2\tau N)\|_1 \le 2\tau \sum_{0 \le n < N} \|r(T - (2n + 1)\tau)\|_1 (1 - \tau K)^{-(N-n)}, \quad \text{where } K := 2 \max_{t \in \mathcal{T}'} \|L_h u(t)\|_{\infty}. $$

Note that $K$ is bounded independently of $\tau \in [0, \tau_*]$ and $h \in [0, h_*]$, since $-L_h u = \text{div}(\mathcal{D}\nabla u) + O(\tau + h^2)$ by theorem 2.8. We assume that $\tau_* K \le 1/2$, up to reducing the maximum timestep $\tau_* > 0$, which implies $(1 - \tau K)^{-1} \le \exp(2\tau K)$, hence $(1 - \tau K)^{-N_\tau} \le \exp(TK)$. Therefore $\forall t \in \mathcal{T}$

$$ \|\hat{\phi}(t)\|_1 \le \exp(\kappa T) R, \qquad \text{where } \kappa := 2\tau \sum_{s \in \mathcal{T}'} \|r(s)\|_1. \qquad (27) $$

**Comparison of $Q(\hat{\phi}; \phi)$ with $\mathcal{L}(\hat{\phi}; \phi)$, and proof of theorem 1.4.** In order to estimate $Q(\hat{\phi}; \phi)$, which is a sum of perspective functions (16), we begin with a convexity lemma.

**Lemma 2.10.** One has $\|\mathcal{P}(m, \rho)\|_1 \ge \mathcal{P}(\|m\|_1, \|\rho\|_1)$, for all $m, \rho \in \mathbb{R}^\mathcal{X}$.

*Proof.* Follows from Jensen's inequality applied to the perspective function, which is convex. $\square$

**Lemma 2.11.** Let $\phi \in \mathbb{R}^\mathcal{X}$ and $\rho := 1 + L_h \phi$. If $\rho \ge 0$ then $\|\rho\|_1 = 1$.

*Proof.* We have $\sum_{x \in \mathcal{X}} L_h \phi(x) = 0$, since $L_h$ is self-adjoint, see theorem 2.5, and since it vanishes on constants, see (12). Thus $\|\rho\|_1 = (2h)^d \sum_{x \in \mathcal{X}} (1 + L_h \phi(x)) = 1$, since $\rho \ge 0$ and $\#(\mathcal{X}) = \#(\mathbb{T}_h)^d = (2h)^{-d}$. $\square$

We obtain

$$ \frac{1}{T} R^2 \le 2\tau \sum_{t \in \mathcal{T}'} \|r(t)\|_1^2 = 2\tau \sum_{t \in \mathcal{T}'} \mathcal{P}(\|r(t)\|_1, \|\rho(t - \tau) + \hat{\rho}(t - \tau)\|_1) \le 2Q(\hat{\phi}; \phi), \quad (28) $$

using successively (i) Cauchy-Schwartz's inequality, (ii) the fact that $\rho(t - \tau) + \hat{\rho}(t - \tau) \ge 0$, otherwise $Q(\hat{\phi}; \phi) = \infty$ in which case the announced inequality is clear, and thus $\|\rho(t - \tau) + \hat{\rho}(t - \tau)\|_1 = \|2 + L_h(\phi + \hat{\phi})\|_1 = 2$ by theorem 2.11, and (iii) theorem 2.10 together with the
non-negativity of the perspective function (hence of the contributions to $Q(\phi; \hat{\phi})$ associated with the other sign choice $\sigma = +$).

On the other hand, we obtain from (24) and (27) that

$$ -\mathcal{L}(\hat{\phi}; \phi) \le C'_{\mathcal{L}}(\tau^2 + h^2)R, \quad \text{with } C'_{\mathcal{L}} := (1 + T) \exp(KT)C_{\mathcal{L}} \quad (29) $$

In the rest of this section, we assume that $0 \ge \mathcal{E}(\phi + \hat{\phi}) - \mathcal{E}(\phi) = \mathcal{L}(\hat{\phi}; \phi) + \mathcal{Q}(\hat{\phi}; \phi)$. The estimates (28) and (29) together yield an upper bound on $R$, since

$$ R^2/(2T) \le \mathcal{Q}(\hat{\phi}; \phi) \le -\mathcal{L}(\hat{\phi}; \phi) \le C'_{\mathcal{L}}(\tau^2 + h^2)R. $$

Therefore $\|\hat{\phi}(t)\|_1 \le C_*(\tau^2 + h^2)$ for all $t \in \mathcal{T}$, where $C_* := 2TC'_\mathcal{L} \exp(KT)$, by (27). Let us conclude the proof of theorem 1.4: we have by the above

$$ \begin{aligned} \inf\{\mathcal{E}(\tilde{\phi}) | \tilde{\phi} \in \Phi\} &= \inf\{\mathcal{E}(\phi + \hat{\phi}) | \hat{\phi} \in \Phi, \mathcal{E}(\phi + \hat{\phi}) \le \mathcal{E}(\phi)\} \\ &= \inf\{\mathcal{E}(\phi + \hat{\phi}) | \hat{\phi} \in \Phi, \forall t \in \mathcal{T}, \|\hat{\phi}(t)\|_1 \le C_*(\tau^2 + h^2)\}. \end{aligned} \quad (30) $$

The latter infimum is attained, since $\mathcal{E}$ is lower-semi-continuous (as a sum of perspective functions) and the optimization domain is compact. The corresponding minimizer is global, is unique since $\mathcal{E}$ is strictly convex by theorem 2.3, and it obeys the estimate announced in theorem 1.4.

### 3 Efficient implementation

We present in this section an implementation of the primal-dual proximal algorithm [CP11, Algorithm 2] for the minimization of the energies (13) and (14) associated with the BBB formulations of the QPME and Burgers's equation, subject to the relevant constraints. (For problem instances of small size, a damped Newton method can also be used, following [KA24]. See the discussion in section 4.) Our approach is closely related to the proximal implementation of the Benamou-Brenier formulation of optimal transport [PPO14], with mild differences related to: the change in boundary conditions w.r.t. the time variable, the additional linear term in the Burgers viscous equation, the use of the harmonic rather than the arithmetic mean of densities, and the anisotropic tensor field $\mathcal{D}$ in the QPME.

Our implementation relies for simplicity on a global Fourier transform in time and space, as opposed to space only in [PPO14]. For that purpose, we introduce the extended time grids

$$ \mathcal{T}_{\text{per}} := \{-T + 2\tau, \dots, T - 2\tau, T\}, \quad \mathcal{T}'_{\text{per}} := \{-T + \tau, T + 3\tau, \dots, T - \tau\}, $$

where the half timestep $\tau > 0$ is fixed throughout this section. The mapping $\rho$ defined on $\mathcal{T}_{\text{per}}$ is implicitly extended by $2T$-periodicity, in such way that $\partial_{\tau}\rho$ is well defined on $\mathcal{T}'_{\text{per}}$ in (32, iii), (34, iii) and (35, left) below. In this section only, we denote by $\epsilon(t) \in \{-1, 1\}$ the sign of $t \in \mathcal{T}'_{\text{per}}$. We fix likewise the half gridscale $h > 0$ and denote by $\mathcal{X} := \mathbb{T}_h^d$ the discrete domain (not to be confused with the characteristic function $\chi$ below). Focusing as a first step on the QPME, we consider the minimization of

$$ \mathcal{E}_{\text{per}}(m, \rho) := \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}'_{\text{per}}} \left[ \mathcal{P}(m(t, x), \rho(t + \tau, x)) - \epsilon(t)m(t, x)u_0(x) \right] + \sum_{x \in \mathcal{X}} \chi(\rho(T, x) - 1). \quad (31) $$

The characteristic function is defined as $\chi(s) = 0$ if $s = 0$, $\chi(s) = \infty$ otherwise, and implements the terminal boundary condition $\rho(T) = 1$. The objective function should be minimized among
all $m \in \mathbb{R}^{\mathcal{T}_{\text{per}}' \times \mathcal{X}}$ and $\rho \in \mathbb{R}^{\mathcal{T}_{\text{per}} \times \mathcal{X}}$ subject to the linear constraints of oddness of $m$, evenness of $\rho$, and continuity equation:

$$m(t) = -m(-t), \ t \in \mathcal{T}_{\text{per}}', \quad \rho(t) = \rho(-t), \ t \in \mathcal{T}_{\text{per}}, \quad \partial_{\tau} \rho(t) = L_h m(t), \ t \in \mathcal{T}_{\text{per}}'. \quad (32)$$

Let us emphasize that the time periodization, the parity constraints, and the time shift $\rho(t+\tau, x)$ of the density in $\mathcal{E}_{\text{per}}$, together, are equivalent to the sum over the two signs $\sigma = \pm$ in the original energy $\mathcal{E}_{\tau h}^{\text{P}}$, see (13). Thus minimizing $\mathcal{E}_{\text{per}}$ subject to (32) is equivalent to minimizing $\mathcal{E}_{\tau h}^{\text{P}}$ on the domain of theorem 1.2.

The chosen minimization algorithm is the iterative numerical method [CP11, Algorithm 2], which requires computing the proximal operator of the energy $\mathcal{E}_{\text{per}}$, and the orthogonal projection onto the linear space defined by (32). As described below, these two steps are made simple by the fully separable structure of (31), and the Fourier-compatible structure of (32). Recall that the proximal operator of a proper convex lower-semi-continuous $f: \mathbb{R}^n \to [-\infty, \infty]$ is defined for any $\lambda > 0$ as

$$\text{prox}_{\lambda f}(x) := \underset{y \in \mathbb{R}^n}{\operatorname{argmin}} \frac{1}{2} \|x - y\|^2 + \lambda f(y).$$

The proximal operator of the perspective function is obtained by solving a univariate third degree polynomial equation, see [CM18, Example 3.8].

**Proposition 3.1.** The proximal value $(m', \rho') = \text{prox}_{\lambda \mathcal{E}_{\text{per}}}(m, \rho)$ is expressed independently for each point and time in terms of the proximal operator of the perspective function:

$$(m'(t, x), \rho'(t + \tau, x)) = \text{prox}_{\lambda P}(m(t, x) + \lambda \epsilon(t) u_0(x), \rho(t + \tau, x))$$

for all $x \in \mathcal{X}$ and $t \in \mathcal{T}_{\text{per}} \setminus \{T - \tau\}$, except for the special case $m'(T - \tau) = (m(T - \tau) + \lambda u_0)/(1 + \lambda)$ and $\rho'(T) = \text{Id}$.

*Proof.* The energy $\mathcal{E}_{\text{per}}$ is a sum of contributions depending on disjoint sets of variables, namely the pairs $(m(t, x), \rho(t + \tau, x))$ for each $x \in \mathcal{X}$ and $t \in \mathcal{T}_{\text{per}}'$. Thanks to this separable structure, the corresponding proximal values can be computed independently. The additional linear term $u_0(x)m(t, x)$ is equivalent to a translation of the projected point. The case $t = T - \tau$ is specific due to the characteristic function enforcing the constraint $\rho(T, x) = 1$, in such way that the objective is quadratic w.r.t. $m(t, x)$. $\square$

We compute the orthogonal projection onto the linear space defined by the constraints (32) in two steps: we first project onto the subspace $E$ defined by the parity constraints (32, i and ii) using the explicit expressions $m_{\text{odd}}(t) = (m(t) - m(-t))/2$ and $\rho_{\text{even}}(t) := (\rho(t) + \rho(-t))/2$, and then project onto the subspace $F$ defined by the continuity equation (32, iii). Indeed if $(m, \rho)$ obey the continuity equation, then so do $(-m(-\cdot), \rho(-\cdot))$, hence also $(m_{\text{odd}}, \rho_{\text{even}}) \in E$ and $(m_{\text{even}}, \rho_{\text{odd}}) \in E^{\perp}$, where $m_{\text{even}} := m - m_{\text{odd}}$ and $\rho_{\text{odd}} := \rho - \rho_{\text{even}}$, in such way that the next lemma applies.

**Lemma 3.2.** Let $E, F$ be closed vector subspaces of a Hilbert space $H$, such that $F = (F \cap E) \oplus (F \cap E^{\perp})$. Then the orthogonal projections on $E, F$ and $E \cap F$, obey $P_{E \cap F} = P_F P_E$.

*Proof.* If $x \in E^{\perp}$, then $P_F P_E(x) = P_F(0) = 0$ and $P_{E \cap F}(x) = 0$ since $E \cap F \subset E$. Let us now consider $x \in E$, and $y := P_{E \cap F}(x)$. Consider $z \in F$, which by assumption can be expressed as $z = z_1 + z_2$ where $z_1 \in E \cap F$ and $z_2 \in E^{\perp} \cap F$. Then $\langle x - y, z_1 \rangle = 0$ by definition of the orthogonal projection on $E \cap F$, and $\langle x - y, z_2 \rangle = 0$ since $x, y \in E$ and $z_2 \in E^{\perp}$. Thus $\langle x - y, z \rangle = 0, \forall z \in F$, and therefore $y = P_F(x)$ as announced, which concludes the proof. $\square$
The projection onto the discretized continuity equation constraint (32, iii) is computed using a global space-time Fourier transform. More precisely, denote for all $\zeta \in \mathcal{T}_{\text{per}}$ and $\xi \in \mathcal{X}$

$$M(\zeta, \xi) := \sum_{\substack{t \in \mathcal{T}_{\text{per}} \\ x \in \mathcal{X}}} m(t + \tau, x) e^{-i\pi(\frac{\zeta t}{\tau} + \frac{(\xi, x)}{h})}, \quad P(\zeta, \xi) := \sum_{\substack{t \in \mathcal{T}_{\text{per}} \\ x \in \mathcal{X}}} \rho(t, x) e^{-i\pi(\frac{\zeta t}{\tau} + \frac{(\xi, x)}{h})}, \quad (33)$$

assuming for simplicity that $2T = 1$. A half time step $m(t + \tau, x)$ is introduced in (33, left) to account for the fact that the momentum $m$ is defined on a staggered grid in time; similarly, the discretized Burgers' equation below involves a density defined on a staggered grid in space, whose Fourier transform is expressed in terms of $\rho(t, x + h)$.

In the case of the standard isotropic Laplacian $-L_h u(x) := \Delta_h u(x) = (2h)^{-2} \sum_{1 \le i \le d} [u(x + 2he_i) - 2u(x) + u(x - 2he_i)]$, where $(e_1, \dots, e_d)$ denotes the canonical basis, the discretized continuity equation (32, iii) reads as follows in Fourier coordinates

$$P(\zeta, \xi) \frac{\exp(2\pi i \zeta) - 1}{2\tau} = M(\zeta, \xi) \sum_{1 \le i \le d} \frac{\sin^2(\pi \xi_i)}{h^2}$$

The projection can be computed independently for each pair of Fourier coordinates $(\zeta, \xi)$ using the following basic lemma in dimension 2. This elementary result, stated without proof, also addresses all projections required in the numerical implementation of the variants below.

**Lemma 3.3.** Let $x, \in \mathbb{C}^d$ and $v \in \mathbb{C}^d \setminus \{0\}$. The orthogonal projection of $x$ onto the line $\mathbb{C}v$ (resp. orthogonal space $(\mathbb{C}v)^\perp$) is obtained as $(v^*x)x/\|v\|^2$ (resp. $x - (v^*x)x/\|v\|^2$), with $v^* := \bar{v}^\top$.

On the other hand, if the coefficients $\lambda^e$ vary over the domain, then we introduce additional variables $m_e, n_e \in \mathbb{R}^{\mathcal{T}_{\text{per}}' \times (\mathcal{X}+he)}$ for each offset $e \in E$. The constraint (32, iii) is removed, and replaced with

$$m_e = \partial_h^e m, \quad n_e = \lambda^e m_e, \quad \partial_\tau \rho = \sum_{e \in E} \partial_h^e n_e, \quad (34)$$

as well as the oddness constraints $m_e(-t) = -m_e(t)$ and $n_e(-t) = -n_e(t)$, $t \in \mathcal{T}_{\text{per}}'$. The constraints (34, i and iii) are added to (32) and treated in space-time Fourier coordinates using theorem 3.3 in dimension $\#(E)+1$ (projecting on a line for (34, i), and on a hyperplane for (34, iii)), after a preliminary projection onto the oddness and evenness constraints. On the other hand, the constraint (34, ii), for all $e \in E$, is added to the objective functional $\mathcal{E}_{\text{per}}(m, \rho)$ in the form of the characteristic function $\chi(n_e - \lambda^e m_e)$, and treated using theorem 3.3 in dimension 2.

**Burgers' equation.** The energy (14) is reformulated as follows:

$$\sum_{x \in \mathcal{X}'} \sum_{\sigma = \pm} \sum_{t \in \mathcal{T}_{\text{per}}'} \left[ \mathcal{P}(\epsilon(t)m_{\sigma}(t, x) + a_{\sigma}(t+\tau, x), \rho(t+\tau, x+\sigma h)) - \epsilon(t)m_{\sigma}(t, x)u_0(x) \right] + \chi(\rho(T, x) - 1),$$

up to the multiplicative constant $\tau h$. We duplicated in this expression the momentum variable: $m_+ = m_-$, and we introduced an auxiliary variable which is also duplicated: $a_+ = a_-$, see (35, right) below. As a result, the energy has a separable structure: it is a sum of contributions depending on the disjoint 5-tuples of variables $(\rho(t + \tau, x + h), m_+(t, x), a_+(t + \tau, x), m_-(t, x + 2h), a_-(t + \tau, x + 2h))$ for all $x \in \mathcal{X}'$ and $t \in \mathcal{T}_{\text{per}}'$. The corresponding proximal values are thus computable independently, using the explicit expression of the proximal operator of the two-dimensional perspective function $(m, \rho) \in \mathbb{R}^2 \times \mathbb{R} \mapsto \mathcal{P}_2(|m|, \rho)$ [CM18, Example 3.8] (indeed
$\mathcal{P}(m_0, \rho) + \mathcal{P}(m_1, \rho) = \mathcal{P}_2((m_0, m_1), \rho)$ by construction), and the following lemma, which is stated without proof since it immediately follows from the co-variance of the proximal operator under the isometric change of coordinates $(m, \rho, a) \mapsto (\frac{m+a}{\sqrt{2}}, \rho, \frac{m-a}{\sqrt{2}})$.

**Lemma 3.4.** Let $f: \mathbb{R}^n \times \mathbb{R} \to [-\infty, \infty]$ be proper, convex, lower-semi-continuous, and positively 2-homogeneous w.r.t. the first variable, and let $F(m, \rho, a) := f(m + a, \rho)$. Then

$$ \text{prox}_{\tau F}(m, \rho, a) = \left( \frac{u+v}{\sqrt{2}}, \tilde{\rho}, \frac{u-v}{\sqrt{2}} \right), \quad \text{where } (u, \tilde{\rho}) := \text{prox}_{2\tau f}\left(\frac{m+a}{\sqrt{2}}, \rho\right), \quad \text{and } v := \frac{m-a}{\sqrt{2}}. $$

The optimization variables obey the discretized continuity equation, and defining constraint,

$$ \partial_{\tau} \rho + \partial_{h} m_{\sigma} = 0, \qquad a_{\sigma} = -\nu \partial_{h} \rho, \qquad (35) $$

for any sign $\sigma = \pm$ (recall that we use duplicated variables $m_+ = m_-$ and $a_+ = a_-$), as well as the parity constraints $m_\sigma(-t) = -m_\sigma(t)$, $t \in \mathcal{T}_{\text{per}}'$, and $\rho(-t) = \rho(t)$ and $a_\sigma(-t) = a_\sigma(t)$, $t \in \mathcal{T}_{\text{per}}$. For the implementation we successively project thanks to theorem 3.2 on the duplication constraints, then the parity constraints, and finally the PDE constraints (35) using space-time Fourier coordinates and theorem 3.3.

Finally, let us mention that the auxiliary variables $a_\sigma \in \mathbb{R}^{\mathcal{T}_{\text{per}}' \times \mathcal{X}}$, $\sigma = \pm$, and the rotation trick theorem 3.4, can be discarded for simplicity and efficiency in the non-viscous case $\nu = 0$.

# 4 Numerical results

We present in this subsection numerical results for the QPME, one- and two-dimensional, isotropic or anisotropic, and Burgers' viscous equation, obtained with the proposed numerical methods, see section 3. Let us immediately acknowledge that, unfortunately and consistently with [Bre20, Theorem 5.2.2] discussed below theorem 1.6, the proposed method appears unsuited to for Burgers' inviscid equation after shock formation, unless we are only interested in the final value at $t = T$ where the entropy solution is still recovered in the continuous setting, while it is not as $t < T$. An extensive literature has been devoted to the numerical solution of the QPME and Burgers' equation, see [BAL$^{+}$18, Li18] and references therein; we not compare our results with or even review the state of the art methods for these PDEs in this paper, which should instead be regarded as an early investigation of the radically new BBB numerical approach. The numerical implementation is available on demand and relies on the Python language and the Taichi library for just-in-time compilation of CPU and GPU code.

The computational cost and the numerical accuracy of the numerical solutions are dictated by three parameters: the half timestep $\tau > 0$, the half gridscale $h > 0$, as well as the number of iterations $N_{\text{prox}}$ of the primal dual proximal method. Indeed, proximal algorithms applied to constrained and/or non-smooth problems are subject to the rather slow $O(1/N_{\text{prox}})$ convergence rate, see the discussion of numerical cost below; in practice we often choose $N_{\text{prox}}$ between 2000 and 12000. For comparison purposes, a damped Newton method (whose step is adjusted using a line search) was also implemented, following [KA24]. It minimizes the same energy functional $\mathcal{E}_{\tau h}^P$ or $\mathcal{E}_{\tau h}^B$ but expressed in terms of the dual variable $\phi$ rather than $(m, \rho)$, and augmented with a logarithmic barrier penalty enforcing the positivity of the problem density $1+L_h\phi$ or $1-\partial_h\phi$, see theorems 1.2 and 1.5. When applicable, this second approach in practice minimizes the energy functional to machine precision in 6 to 270 iterations depending on the experiment; however, it does not scale to large problem instances since it relies on the inversion of ill-conditioned linear systems of with $N_\tau N_h^d$ unknowns, which quickly becomes intractable as $N_\tau$ and $N_h$ increase.
<figure><img src="image_2.png" /><figcaption>Figure 1: Top left : Barenblatt profile (36, left), with parameter $\gamma = 1$, on the time interval $[10^{-4}, 2 \times 10^{-3}]$. Bottom left : explicit solution to Burgers' equation (36, right), with $\nu = 10^{-2}$ and $\delta = \exp(\text{Re}) - 1$ with $\text{Re} = 5$ the Reynolds number. Top and bottom, center and right: numerically computed auxiliary variables $m$ and $\rho$, for the corresponding BBB formulation (8).</figcaption></figure>

We use periodic boundary conditions for simplicity, and consistently with the theoretical analysis. The extension to Neumann boundary conditions for the QPME is straightforward, by doubling the physical domain size and extending the solution by reflection (in practice, this amounts to replacing the FFT with a cosine transform in section 3). Dirichlet-type boundary conditions for Burgers' equation are discussed in [KA24].

For reference, let us recall the Barenblatt [Vá07] explicit solution $u^P$ to the QPME, which is non-smooth and compactly supported, and a classical explicit solution $u^B$ to Burgers viscous equation obtained via the Hopf-Cole transform [Whi11]:

$$u^P(t, x) := \frac{2}{t^\alpha} \max \left\{ 0, \gamma - \frac{\beta \|x\|^2}{4 t^{2\beta}} \right\}, \qquad u^B(t, x) := \sqrt{\frac{\nu}{\pi t}} \frac{\delta \exp\left(\frac{-x^2}{4\nu t}\right)}{1 + \frac{\delta}{2} \operatorname{erfc}\left(\frac{x}{\sqrt{4\nu t}}\right)}, \quad (36)$$

where $\alpha := d/(d+2)$, $\beta := 1/(d+2)$, and $\gamma, \delta > 0$ are arbitrary positive constants. We denoted by $\operatorname{erfc}(x) := \frac{2}{\sqrt{\pi}} \int_x^\infty \exp(-s^2) ds$ the complementary error function. In the process of running these experiments, we conjectured, and then proved, a closed form expression for the dual potential associated with the Barenblatt solution, in the BBB formulation of the QPME, see section D.

**Discussion of numerical cost.** The proposed BBB numerical method for the QPME has computational cost $\tilde{O}(N_{\text{prox}}N_{\tau}N_h^d)$ where the “tilde” ignores the logarithmic factors associated with the fast Fourier transform. Memory usage is $O(N_{\tau}N_h^d)$, and typical accuracy $O(N_{\text{prox}}^{-1} + N_{\tau}^{-2} + N_h^{-2})$. For comparison, the simple explicit scheme (49, with $\theta = 1$) has computational cost $O(N_{\tau}N_h^d)$, memory usage $O(N_h^d)$, and typical accuracy $O(N_{\tau}^{-1} + N_h^{-2})$, but is subject to the CFL condition $N_{\tau} \gtrsim N_h^2$. Let us acknowledge that the slow convergence of the primal-dual solver, and the higher memory usage, make the BBB approach the costlier one usually, except in specific edge cases such as the use of very large timesteps discussed in the second paragraph of section 4.1. Some optimizations of the BBB numerical implementation are proposed in [KA24]
<figure><img src="image_3.png" /><figcaption>Figure 2: Left: log-log plot of numerical error, between the Barenblatt profile with the parameters of fig. 1, and the numerical solution of the discretized QPME, using $4 \le N_\tau \le 200$ timesteps and $N_h = 5N_\tau$ discretization points; note the second order convergence rate in the $L^1$ norm. Center: likewise for Burgers viscous equation, using $4 \le N_\tau \le 200$ timesteps and $N_h = 5N_\tau$ discretization points; note the second order convergence in the $L^\infty$ norm. Right: numerical solution of the QPME obtained with $N_h = 128$ discretization points, and very few $N_\tau = 10$ timesteps, see the discussion of large timesteps in section 4.1; note that a Barenblatt-like profile is obtained at the final time.</figcaption></figure>

(for Burgers' equation), such as concatenating the numerical solutions obtained over small time sub-intervals, and using a well chosen initial guess and "base state", see theorem 1.8.

## 4.1 Anisotropic Quadratic Porous medium equation

**Accuracy validation.** In order to assess the second order accuracy of the proposed numerical method, we reproduce the known Barenblatt solution (36, left), with parameter $\gamma = 1$ and from the time $T_0 = 10^{-4}$ to $T = 10^{-3}$, whose spatial support is contained in $[-1/2, 1/2]$. We numerically solve a sequence of small problem instances, where the number $N_\tau$ of timesteps ranges from 4 to 200, the number of spatial discretization points is chosen as $N_h = 5N_\tau$; the energy functional of interest (13) is minimized to machine precision using a damped Newton method (60 to 270 iterations). In contrast with our convergence result theorem 1.4, the problem solution is *non-smooth, non-positive* since it has compact support, and we evaluate the numerical error on the *primal variable u* rather than the dual potential $\phi$. We nevertheless observe a second order convergence rate in the $L^1$ norm, which is consistent with the scheme order; convergence in the $L^\infty$ norm is slower, as expected in view of the non-smoothness of the solution, see fig. 2. The exact density $\rho(t, x)$ features two discontinuities w.r.t. $x \in [0, \infty[$ for each $t \in [T_0, T[$: it is constant on $[0, R(t)][$, then affine on $]R(t), R_T[$, then again constant on $]R_T, \in[$, with jumps at $R(t)$ and $R_T$, whose closed form expression is obtained in section D. Convergence for the variable $\rho$ is expectedly slow, and we see some spurious numerical oscillations near the first discontinuity, see fig. 1, but this does not affect the ratio $u = m/\rho$ whose support ends at $R(t)$.

**Large time steps.** The solution of the QPME, with compactly supported initial data $u_0 \in C_0^0(\mathbb{R})$, converges (when suitably rescaled) as $t \to \infty$ to the Barenblatt profile [Vá07]. We illustrate this phenomenon on fig. 2 (right), where the initial data array $u_0: \mathbb{T}_h \to \mathbb{R}$ vanishes except for three arbitrary positive values ($\cdots, 3/10, 1, 0, 0, 0, 7/10, \cdots$) close to the domain center. Setting the final time to $T = 1$, using only $N_\tau = 10$ timesteps and (arbitrarily) $N_h = 128$ discretization points, the final time solution resembles as expected the bell-shaped Barenblatt profile.

For comparison, the explicit scheme (49, $\theta = 1$) is subject to the CFL condition $\max(u)\tau \le h^2$
<figure><img src="image_4.png" /><figcaption>Figure 3: I: Image of the Cartesian grid by the chosen diffeomorphism $\phi$. II: Visualization of the anisotropy via the ellipses $\mathcal{E}(x) := \{v \in \mathbb{R}^2 \mid \|D\phi(x)v\| = r\}$, for some $r > 0$. III: Numerical solution of the anisotropic QPME, on a small grid with $N_{\tau} = 12$ and $N_{h} = 48$. IV: Difference between the numerical and the exact solution, which is obtained as the composition of the Barenblatt profile $u^P$ with the diffeomorphism $\phi$.</figcaption></figure>

(equivalently $N_{\tau} \ge 2 \max(u) N_h^2$), and thus diverges unless $N_{\tau} \gtrsim 33000$ timesteps are used in this example. The semi-implicit scheme (49, $\theta = 1/2$), which can be solved one timestep at a time but also corresponds to the choice of the arithmetic rather than harmonic mean of the density in the discretized energy functional, compare (13) with (45), diverges unless $N_{\tau} \ge 5000$ timesteps are used in this example (larger timesteps can be used with smooth initial conditions it appears). Finally, the implicit scheme (49, $\theta = 0$) is unconditionally stable, hence $N_{\tau} = 10$ timesteps can be used as well, but it is only first order accurate w.r.t. time.

**Two-dimensional anisotropic.** Let $\psi \in C^2(\mathbb{R}^d, \mathbb{R}^d)$ be a diffeomorphism whose Jacobian matrix obeys $|\det(D\psi)| = 1$ identically, i.e. $\psi$ preserves the Lebesgue measure, and let $u(t, x)$ obey the isotropic QPME. Then $(t, x) \mapsto u(t, \psi(x))$ obeys the anisotropic QPME associated with the diffusion tensor field $D(x) := (D\psi(x)^{\top}D\psi(x))^{-1}$, by the usual change of variable formulas for the gradient and divergence. In this experiment, we choose $u$ as the Barenblatt profile with $\gamma = 1$ on the time interval $[10^{-5}, 10^{-4}]$, and we define $\psi(x_0, y_0) = (x_2, y_2)$ via $x_1 = x_0 + \epsilon \sin(2\pi y_0 + 1)$, $y_1 = y_0 + \epsilon \sin(4\pi x_1 + 5)$, $x_2 = x_1 + \epsilon \sin(4\pi y_1 + 3)$, $y_2 = x_2 + \epsilon \sin(2\pi x_2 + 2)$, with $\epsilon = 0.035$. By construction, $\psi$ is a 1-periodic measure preserving diffeomorphism, illustrated on fig. 3, with maximum distortion ratio $\approx 3.5$ (ratio of the largest to the smallest singular value of $D\phi$). Selling's decomposition of this tensor field $D$ yields the collection $E := \{(1, 0), (0, 1), (1, 1), (1, -1)\}$ of offsets, and a piecewise smooth [BBM22, Proposition B.5] weight $\lambda^e$ for each $e \in E$ obeying (11, right), as required for the construction of the proposed scheme (12). A modified decomposition introduced in [BBM23, Theorem 1.8] yields smooth weights $\tilde{\lambda}^e$, fitting theorem 1.3, but uses slightly more offsets $e \in \tilde{E} = E \cup \{(2, 1), (2, -1), (1, 2), (-1, 2)\}$, which increases the memory usage and the numerical cost for no observable benefit in this experiment. A numerical solution obtained on a small grid $N_{\tau} = 12, N_h = 48$, is shown on fig. 3.

**Scalability.** The proposed BBB discretization of the QPME works in arbitrary space dimension $d$, but its numerical cost is of course subject to the curse of dimensionality, especially so since the numerical method manipulates space-time $(d+1)$-dimensional arrays. The primal-dual proximal optimization algorithm is embarrassingly parallel in time and space, assuming a suitable implementation of the FFT, and thanks to the Taichi library our numerical implementation can be executed on the CPU or GPU without modification. In order to investigate the scalability of the proposed method, we ran the above two-dimensional anisotropic QPME experiment
<figure><img src="image_5.png" /><figcaption>Figure 4: Top: Burgers' equation with viscosity $\nu = 10^{-3}$, initial condition $u_0 = \mathbb{1}_{0.15 \le x \le 0.65}$, evolution time $T = 0.2$, $N_\tau = 512$, $N_h = 1024$, $N_{\text{prox}} = 6000$. Bottom: Burgers' equation with viscosity $\nu = 10^{-3}$, initial condition $u_0 = \frac{1}{2}(\mathbb{1}_{0.15 \le x \le 0.65} + \mathbb{1}_{0.3 \le x \le 0.5})$, evolution time $T = 0.4$, $N_\tau = 512$, $N_h = 1024$, $N_{\text{prox}} = 6000$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Failed attempt at the numerical solution of Burgers inviscid equation ($\nu = 0$), initial condition $u_0 = \mathbb{1}_{0.15 \le x \le 0.65}$, evolution time $T = 0.2$, $N_\tau = 256$, $N_h = 512$, $N_{\text{prox}} = 40000$. The numerical solution is reasonably good at the final time $t = 0.2$, but is incorrect, or non-converging, at earlier times. In particular, the initial condition ($t = 0.0$, dark blue) is not reproduced.</figcaption></figure>

on a larger grid: with $N_\tau = 64$ and $N_h = 128$, computations take $7.9s$ per proximal iteration on the CPU, and $0.069s$ on the GPU, which is a $\approx 110 \times$ performance improvement. Similarly, the isotropic two-dimensional QPME with $N_\tau = 128$ and $N_h = 512$ takes $10.5s$ per proximal iteration on the CPU, and $0.085s$ on the GPU, which is a $\approx 120 \times$ performance improvement. Approximately 2000 proximal iterations are needed to achieve good accuracy. Laptop equipped with an Intel Core i7, 8 core, 2.3Ghz CPU, and an Nvidia 2060 Max-Q GPU. The chosen dimensions $N_\tau \times N_h \times N_h$ are the largest powers of two such that the computations fit in the GPU 8GB memory.

## 4.2 Burgers' equation

**Large viscosity.** In order to assess the accuracy of the proposed numerical method, we reproduce the known solution (36, right), with viscosity $\nu = 10^{-2}$ and parameter $\delta = \exp(\text{Re}) - 1$ where $\text{Re} = 5$ is regarded as a Reynolds number, from the time $T_0 = 0.1$ to $T = 1$, and on the interval $[-0.6, 0.9]$ outside of which the solution vanishes up to machine precision. We solve a
<figure><img src="image_7.png" /><figcaption>Figure 6: Left: exact solution of Burgers’ equation on the time interval $[0.1, 1]$ with $\nu = 10^{-3}$ and $\delta = \exp(\text{Re}) - 1$ where $\text{Re} = 50$ is the Reynolds number. Center: numerical solution obtained with $N_{\tau} = 64$ timesteps and $N_h = 200$ discretization points. Right: numerical solution obtained with only $N_{\tau} = 8$ timesteps, and $N_h = 200$, see the discussion on large timesteps in section 4.2.</figcaption></figure>

sequence of small problem instances, where the number of timesteps ranges from 4 to 400 and the number of spatial discretization points is chosen as $N_h = 5N_{\tau}$, and we minimize the energy functional of interest to machine precision using a Newton method, which only requires 7 iterations here independently of the problem size. In contrast with the assumptions of theorem 1.7, the problem solution is not really positively bounded below, since the quadratic exponential $\exp(-\frac{x^2}{4\nu t})$ in (36) falls below machine precision on the domain boundary, and error is measured on the primal variable $u$ rather than the dual potential $\phi$. Nevertheless, we observe a second order convergence rate in the $L^{\infty}$ norm, see fig. 2, which could be expected since the solution is smooth.

**Small viscosity.** We illustrate on fig. 4 the numerical solution obtained with a small viscosity $\nu = 10^{-3}$, and when the initial condition is the characteristic function of a compact interval, or the average of two such characteristic functions. This creates rarefaction waves on the left, and shocks on the right, which behave and merge as expected. From the theoretical standpoint, the density $\rho^T$ is positively bounded below, for any final evolution time $T$ in the BBB formulation, see theorem 1.8. Let us acknowledge however that if $T$ is chosen too large, then the values of $\rho$ fall below machine precision (32bit floats are used for these GPU experiments), which raises numerical convergence issues, not illustrated here, see the discussion at the end of section E.

Depending on one’s perspective, the last experiments may be regarded as negative results.

**Non-viscous.** The solution $u^T$ of the BBB formulation of Burgers’ equation on the interval $[0, T]$, differs in general from the viscosity solution $u$ if $T > T_{\text{shock}}$, except at the final time: $u^T(T, \cdot) = u(T, \cdot)$. More precisely, [Bre20, Theorem 5.5.2] characterizes $u^T$ as the unique shock free solution of Burgers’ equation on $[0, T]$ such that $u^T(T, \cdot) = u(T, \cdot)$. In our experiment (using the proximal solver on the GPU), the numerically computed approximation of $u^T$ does match the viscosity solution at the final time $T$, but typically fails to converge for smaller values $t \in [0, T[$, likely because the density $\rho$ falls below machine precision, see fig. 5.

A strategy around this issue is proposed in [KA24]. In that work, the BBB formulation is successively discretized on a collection of small time sub-intervals, with some overlapping, rather a single large interval. In addition, crucially, a well chosen base state is introduced in the energy functional, see theorem 1.8. It may either be generated numerically, or defined a-priori e.g. as the solution of Burgers’ equation with a small but positive viscosity parameter [KA24, §5.4].
**Large time steps.** We consider the exact solution (36, right), with smaller viscosity $\nu = 10^{-3}$ and higher Reynolds number $Re = 50$, on the time interval $[0.1, 1]$ and the domain $[-0.15, 0.55]$, see fig. 6. This solution is reasonably well reproduced using the proposed BBB formulation with $N_\tau = 64$ timesteps, $N_h = 200$ discretization points, here solved using the damped Newton method. In order to illustrate the unconditional stability of the scheme, we also solve the BBB formulation with excessively few $N_\tau = 8$ timesteps. The numerical solution accuracy is degraded, as could be expected, but interestingly this degradation includes the initial condition which is not well reproduced: the ‘right angle’ is turned into a symmetric ‘hat’, see fig. 6 (right). It is not clear to us whether this is related to the phenomenon of the previous paragraph.

# 5 Conclusion

In this paper, we presented the first numerical convergence analysis of the ballistic Benamou-Brenier (BBB) formulation of non-linear evolution PDEs, applied to the (possibly anisotropic) quadratic porous medium equation (QPME) and the (possibly viscous) Burgers’ equation. The obtained schemes are second order accurate w.r.t. the timestep and the gridscale, and are free of any Courant-Friedrichs-Levy condition relating these parameters. They take the form of an optimization problem, whose structure is reminiscent of a mean field game, and which can be solved using a proximal primal-dual method and space-time fast Fourier transforms.

Numerical experiments with the QPME show that the approach efficiently handles large timesteps, and strong anisotropies, while retaining second order accuracy and being easily parallelizable. Limitations come from memory usage, due to the use of space-time global solution arrays, and from the slow convergence of the proximal optimizer. The approach applies to Burgers’ equation with a (possibly small) positive viscosity, but suffers from delicate theoretical obstructions in the inviscid case [Bre18].

Numerous perspectives are open for future work. The PDE models may be generalized, by introducing a right hand side or non-periodic boundary conditions. Alternative discretizations may be considered, using e.g. triangulations instead of a Cartesian grid. The convergence analysis which is currently limited to smooth and uniformly positive solutions might be extended so as to accommodate e.g. the Barenblatt profile. The primal-dual optimization method, which currently limits the numerical efficiency of our approach, might be improved using e.g. the uniform convexity properties underlying our convergence analysis. The introduction of a weight in the energy functional, decaying exponentially with time following [Vor25], should be investigated and the numerical method suitably adapted. Finally, and most importantly, the numerical BBB approach will be investigated in the case of systems of PDEs, related to e.g. fluid mechanics, which are a fundamental motivation of this theoretical framework [Bre20, Vor22].

**Acknowledgement.** The authors would like to thank Yann Brenier and Thomas Gallouët for numerous constructive discussions.

# References

[Ach23] Amit Acharya. Variational principles for nonlinear PDE systems via duality. *Quarterly of Applied Mathematics*, 81(1):127–140, 2023.

[AS24] Amit Acharya and Ambar N Sengupta. Action principles for dissipative, non-holonomic Newtonian mechanics. *Proceedings of the Royal Society A*, 480(2293):20240113, 2024.
[ASZ24] Amit Acharya, Bianca Stroffolini, and Arghir Zarnescu. Variational dual solutions for incompressible fluids. *arXiv preprint arXiv:2409.04911*, 2024.

[BAL$^{+}$18] Mayur P Bonkile, Ashish Awasthi, C Lakshmi, Vijitha Mukundan, and VS Aswin. A systematic literature review of Burgers' equation with recent advances. *Pramana*, 90:1–21, 2018.

[BB09] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. *Numerische Mathematik*, 84(3):375–393, 2009.

[BBM22] J Frederic Bonnans, Guillaume Bonnet, and Jean-Marie Mirebeau. A linear finite-difference scheme for approximating Randers distances on Cartesian grids. *ESAIM: Control, Optimisation and Calculus of Variations*, 28:45, 00 2022.

[BBM23] Frédéric Bonnans, Guillaume Bonnet, and Jean-Marie Mirebeau. Monotone discretization of anisotropic differential operators using voronoi's first reduction. *Constructive Approximation*, pages 1–61, 2023.

[BM22] Yann Brenier and Iván Moyano. Relaxed solutions for incompressible inviscid flows: A variational and gravitational approximation to the initial value problem. *Philosophical Transactions of the Royal Society A*, 380(2219):20210078, 2022.

[Bre18] Yann Brenier. The initial value problem for the Euler equations of incompressible fluids viewed as a concave maximization problem. *Communications in Mathematical Physics*, 364:579–605, 2018.

[Bre20] Yann Brenier. Examples of hidden convexity in nonlinear PDEs. 2020.

[CM18] Patrick L Combettes and Christian L Müller. Perspective functions: Proximal calculus and applications in high-dimensional statistics. *Journal of Mathematical Analysis and Applications*, 457(2):1283–1306, 2018.

[CP11] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. *Journal of Mathematical Imaging and Vision*, 2011.

[DiP79] Ronald J DiPerna. Uniqueness of solutions to hyperbolic conservation laws. *Indiana University Mathematics Journal*, 28(1):137–188, 1979.

[FM14] Jérôme Fehrenbach and Jean-Marie Mirebeau. Sparse non-negative stencils for anisotropic diffusion. *Journal of Mathematical Imaging and Vision*, 49(1):123–147, 2014.

[HT$^{+}$19] M Huesmann, D Trevisan, et al. A Benamou-Brenier formulation of martingale optimal transport. *BERNOULLI*, 25(4 A):2729–2757, 2019.

[K$^{+}$05] Nicolai V Krylov et al. The rate of convergence of finite-difference approximations for bellman equations with lipschitz coefficients. *Applied Mathematics and Optimization*, 52, 2005.

[KA24] Uditnarayan Kouskiya and Amit Acharya. Inviscid Burgers as a degenerate elliptic problem. *Quarterly of Applied Mathematics*, 2024.
[Li18] Xianping Li. Anisotropic mesh adaptation for finite element solution of anisotropic porous medium equation. *Computers & Mathematics with Applications*, 75(6):2086–2099, 2018.

[LJ09] Camillo De Lellis and László Székelyhidi Jr. The Euler equations as a differential inclusion. *Annals of mathematics*, pages 1417–1436, 2009.

[LS10] Camillo De Lellis and László Székelyhidi. On admissibility criteria for weak solutions of the Euler equations. *Archive for rational mechanics and analysis*, 195:225–260, 2010.

[LZ24] Bowen Li and Jun Zou. On the convergence of discrete dynamic unbalanced transport models. *ESAIM: Mathematical Modelling and Numerical Analysis*, 58(3):957–992, 2024.

[NT21] Andrea Natale and Gabriele Todeschi. Computation of optimal transport with finite volumes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55(5):1847–1871, 2021.

[PPO14] Nicolas Papadakis, Gabriel Peyré, and Edouard Oudet. Optimal transport with proximal splitting. *SIAM Journal on Imaging Sciences*, 7(1):212–238, 2014.

[SGA24] Siddharth Singh, Janusz Ginster, and Amit Acharya. A hidden convexity of nonlinear elasticity. *Journal of Elasticity*, 156(3):975–1014, 2024.

[VC03] Villani and Cédric. Topics in Optimal Transportation. 58, 2003.

[Vor22] Dmitry Vorotnikov. Partial differential equations with quadratic nonlinearities viewed as matrix-valued optimal ballistic transport problems. *Archive for Rational Mechanics and Analysis*, 243(3):1653–1698, 2022.

[Vor25] Dmitry Vorotnikov. Hidden convexity and Dafermos’ principle for some dispersive equations. *arXiv*, 2025.

[Vá07] Juan Luis Vázquez. *The porous medium equation: mathematical theory*. Oxford university press, 2007.

[Whi11] Gerald Beresford Whitham. *Linear and nonlinear waves*. John Wiley & Sons, 2011.

## A Burgers’ equation

This section is devoted to the proof theorem 1.7, our convergence result for the discretization of Burgers’ equation. Incidentally, this shows that the proof template of section 2 for the QPME can accommodate variations, in this case (i) a non-linearity associated with a first order operator (rather than second order), (ii) an optional viscous linear term, and (iii) a staggered grid discretization in space, in addition to time.

**Notation A.1.** We fix in the following the half timestep $\tau > 0$ and the half gridscale $h > 0$. For readability and in order to avoid double sub- or super-scripts, we denote by $\mathcal{E} := \mathcal{E}_{\tau h}^B$ the energy of interest, and by $\Phi := \Phi_{\tau h}^B$ its domain, see theorem 1.5. Likewise we denote by $\mathcal{T} := \mathcal{T}_{\tau}$, $\mathcal{T}' := \mathcal{T}'_{\tau}$, $\mathcal{X} := \mathbb{T}_h$ and $\mathcal{X}' := \mathbb{T}'_h$ the centered and staggered grids, in time and space, see (9). Finally, we denote $\| \cdot \|_1 = \| \cdot \|_{l^1(\mathcal{X})}$, see (10), and $\tilde{\mathcal{T}} := \mathcal{T} \setminus \{0, T\} = \{2\tau, \dots, T - 2\tau\}$.
Using the exact expansion of the perspective function, presented in theorem 2.1, we express the energy of interest as $\mathcal{E}(\phi+\hat{\phi}) = \mathcal{E}(\phi)+\mathcal{L}(\hat{\phi}; \phi)+\mathcal{Q}(\hat{\phi}; \phi)$ for any $\phi, \hat{\phi} \in \Phi$ such that $\mathcal{E}(\phi) < \infty$. The linear and second order (non-quadratic) terms are

$$ \mathcal{L}(\hat{\phi}; \phi) := \tau h \sum_{t \in T'} \sum_{\substack{\sigma_x = \pm \\ \sigma_t = \pm}} \left[ u_{\sigma_x}^{\sigma_t}(t, x) (\hat{m}(t, x) - \nu \partial_h \hat{\rho}(t + \sigma_t \tau, x)) - \frac{1}{2} u_{\sigma_x}^{\sigma_t}(t, x)^2 \hat{\rho}(t + \sigma_t \tau, x + \sigma_x h) \right] + 2h \sum_{x \in \mathcal{X}} u_0(x) \hat{\phi}(0, x), $$

$$ \mathcal{Q}(\hat{\phi}; \phi) := \tau h \sum_{t \in T'} \sum_{\substack{\sigma_x = \pm \\ \sigma_t = \pm}} \mathcal{P}(\hat{m}(t, x) - \nu \partial_h \hat{\rho}(t + \sigma_t \tau, x) - \hat{\rho}(t + \sigma_t \tau, x + \sigma_x h) u_{\sigma_x}^{\sigma_t}(t, x), \quad (37) \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \rho(t + \sigma_t \tau, x + \sigma_x h) + \hat{\rho}(t + \sigma_t \tau, x + \sigma_x h)), $$

$$ \text{where } m = \partial_{\tau} \phi, \quad \rho = 1 - \partial_h \phi, \quad \hat{m} := \partial_{\tau} \hat{\phi}, \quad \hat{\rho} := -\partial_h \hat{\phi}, \quad u_{\sigma_x}^{\sigma_t}(t, x) := \frac{m(t, x) - \partial_h \rho(t + \sigma_t, x)}{\rho(t + \sigma_t \tau, x + \sigma_x h)}, $$

whose denominator is positive for sufficiently small $h > 0$ by theorems A.3 and 1.6 below, and where $\sigma_t, \sigma_x \in \{+, -\}$ are signs. Note that $\phi, \hat{\phi}, \partial_h \rho, \partial_h \hat{\rho} \in \mathbb{R}^{T \times \mathcal{X}}$ are centered in time and space, $m, \hat{m}, u_{\sigma_x}^{\sigma_t} \in \mathbb{R}^{T \times \mathcal{X}}$ are staggered in time but centered in space, and $\rho, \hat{\rho} \in \mathbb{R}^{T \times \mathcal{X}'}$ are centered in time but staggered in space, since we use staggered finite difference operators $\partial_{\tau}$ and $\partial_h$.

**Proposition A.2.** The energy $\mathcal{E}$ is strictly convex over the domain $\Phi$.

*Proof.* We proceed similarly to theorem 2.3 and note that $\mathcal{E}$ is convex since it is a sum of perspective functions (composed with affine mappings), which are convex. Consider $\phi, \hat{\phi} \in \Phi$ such that $s \in [0, 1] \mapsto \mathcal{E}(\phi + s\hat{\phi})$ is finitely valued and affine. We note that $\hat{\phi}(T, x) = 0$ for all $x \in \mathcal{X}$, by definition of the domain $\Phi$, and assume for induction that $\hat{\phi}(T - 2n\tau, x) = 0$ for all $x \in \mathcal{X}$, for some $0 \le n < N_{\tau} = T/(2\tau)$. Letting $t := T - (2n + 1)\tau$ we obtain $\hat{m}(t, x) = \hat{\phi}(t - \tau, x)/(2\tau)$ and $\hat{\rho}(t + \tau, x + h) = 0$ for all $x \in \mathcal{X}$. Consider the contribution to $\mathcal{Q}$ associated with the signs $\sigma_t = +$ and $\sigma_x = +$ (arbitrarily in the case of $\sigma_x$), which reads

$$ \mathcal{P}(s\hat{m}(t, x) - \nu s \partial_h \hat{\rho}(t + \tau, x) - s\hat{\rho}(t + \tau, x + h)u_{+}^{\tau}(t, x), \rho(t + \tau, x + h) + s\hat{\rho}(t + \tau, x + h)) \\ = \frac{s^2}{2\tau} \mathcal{P}(\hat{\phi}(t - \tau, x), \rho(t + \tau, x + h)). $$

By assumption, this is an affine function of $s \in [0, 1]$. Therefore $0 = \hat{\phi}(t - \tau, x) = \hat{\phi}(T - 2\tau(n + 1), x)$, hence $\hat{\phi} = 0$ by an immediate induction, thus the announced strict convexity. $\square$

## A.1 Upper estimate of the linear term

We estimate in this subsection the linear term $\mathcal{L}(\hat{\phi}; \phi)$ in terms of the $\ell^1(\mathcal{X})$ norm of the perturbation $\hat{\phi}$, see the estimate (40) which holds under theorem 1.6. As a first step, we recall that the finite difference operator $\partial_h$ is skew-adjoint in the sense that

$$ \sum_{x \in \mathcal{X}} u(x) \partial_h v(x + h) = \frac{1}{2h} \sum_{x \in \mathcal{X}} u(x) v(x + 2h) - \frac{1}{2h} \sum_{x \in \mathcal{X}} u(x) v(x) = - \sum_{x \in \mathcal{X}} v(x) \partial_h u(x - h), \quad (38) $$
for all $u, v \in \mathbb{R}^{\mathcal{X}}$, similarly to theorem 2.5. We rewrite $\mathcal{L}(\hat{\phi}; \phi)$ using (38) and the discrete summation by parts w.r.t. time of theorem 2.4: $-\mathcal{L}(\hat{\phi}; \phi) =$

$$
\begin{align}
& h \sum_{x \in \mathcal{X}} \sum_{\sigma_x = \pm} \left[ \frac{u_{\sigma_x}^{-}(\tau, x) + u_{\sigma_x}^{+}(\tau, x)}{2} \hat{\phi}(0, x) + \tau \nu u_{\sigma_x}^{-}(\tau, x) \partial_h \hat{\rho}(0, x) + \frac{\tau}{2} u_{\sigma_x}^{-}(\tau, x)^2 \hat{\rho}(0, x + \sigma_x h) - u_0(x) \hat{\phi}(0, x) \right] \nonumber \\
& \quad + \tau h \sum_{\substack{t \in \hat{\tau} \\ x \in \mathcal{X}}} \sum_{\substack{\sigma_x = \pm \\ \sigma_t = \pm}} \left[ \partial_{\tau} u_{\sigma_x}^{\sigma_t}(t, x) \hat{\phi}(t, x) + u_{\sigma_x}^{\sigma_t}(t - \sigma_t \tau, x) \nu \partial_h \hat{\rho}(t, x) + \frac{1}{2} u_{\sigma_x}^{\sigma_t}(t - \sigma_t \tau, x)^2 \hat{\rho}(t, x + \sigma_x h) \right] \nonumber \\
&= h \sum_{x \in \mathcal{X}} \sum_{\sigma_x = \pm} \left( \frac{u_{\sigma_x}^{-}(\tau, x) + u_{\sigma_x}^{+}(\tau, x)}{2} - \tau \nu \partial_{hh} u_{\sigma_x}^{-}(\tau, x) + \frac{\tau}{2} \partial_h [(u_{\sigma_x}^{-})^2](\tau, x - \sigma_x h) - u_0(x) \right) \hat{\phi}(0, x) \nonumber \\
& \quad + \tau h \sum_{\substack{t \in \hat{\tau} \\ x \in \mathcal{X}}} \sum_{\substack{\sigma_x = \pm \\ \sigma_t = \pm}} \left( \partial_{\tau} u_{\sigma_x}^{\sigma_t}(t, x) - \nu \partial_{hh} u_{\sigma_x}^{\sigma_t}(t - \sigma_t \tau, x) + \frac{1}{2} \partial_h [(u_{\sigma_x}^{\sigma_t})^2](t - \sigma_t \tau, x - \sigma_x h) \right) \hat{\phi}(t, x). \tag{39}
\end{align}
$$

In the rest of this section, we work under theorem 1.6, and thus assume that $\phi \in C^\infty([0, T] \times \mathbb{T}^d, \mathbb{R})$ is a smooth function. We distinguish the variable $m := \partial_\tau \phi$ from $m := \partial_t \phi$, and likewise

$$
\rho := 1 - \partial_h \phi, \quad \rho := 1 - \partial_x \phi, \quad u_{\sigma_x}^{\sigma_t}(t, x) := \frac{m(t, x) - \partial_h \rho(t + \sigma_t, x)}{\rho(t + \sigma_t \tau, x + \sigma_x h)}, \quad \mathbf{u} := \frac{\mathbf{m} - \partial_x \rho}{\rho},
$$

and note that $m$ and $u_{\sigma_x}^{\sigma_t}$ are defined over $[\tau, T - \tau] \times \mathbb{T}$ whereas the other quantities are defined over the full domain $[0, T] \times \mathbb{T}$. By theorem 1.6 one has $\rho > 0$ uniformly, and $\mathbf{u}$ obeys Burger's equation $\partial_t \mathbf{u} + \frac{1}{2} \partial_x \mathbf{u}^2 = \nu \partial_{xx} \mathbf{u}$ for some non-negative viscosity coefficient $\nu \ge 0$.

**Proposition A.3.** Under theorem 1.6 and with the notations $m, \rho, u_{\sigma_x}^{\sigma_t}, \mathbf{m}, \rho, \mathbf{u}$ above:

$$
\begin{align*}
& m = \mathbf{m} + \mathcal{O}(\tau^2), \quad \rho = \rho + \mathcal{O}(h^2), \quad u_{\sigma_x}^{\sigma_t} = \mathbf{u} + \mathcal{O}(\tau + h), \quad \sum_{\sigma_t, \sigma_x = \pm} u_{\sigma_x}^{\sigma_t} = 4\mathbf{u} + \mathcal{O}(\tau^2 + h^2), \\
& \partial_{\tau} u_{\sigma_x}^{\sigma_t} = \partial_t \mathbf{u} + \mathcal{O}(\tau + h), \quad \partial_{hh} u_{\sigma_x}^{\sigma_t} = \partial_{xx} \mathbf{u} + \mathcal{O}(\tau + h), \quad \partial_h (u_{\sigma_x}^{\sigma_t})^2 = \partial_x \mathbf{u}^2 + \mathcal{O}(\tau + h), \\
& \sum_{\substack{\sigma_x=\pm \\ \sigma_t=\pm}} \partial_{\tau} u_{\sigma_x}^{\sigma_t} = 4\partial_t \mathbf{u} + \mathcal{O}(\tau^2 + h^2), \quad \sum_{\substack{\sigma_x=\pm \\ \sigma_t=\pm}} \partial_{hh} u_{\sigma_x}^{\sigma_t} (t-\sigma_t\tau, x) = 4\partial_{xx} \mathbf{u}(t,x) + \mathcal{O}(\tau^2 + h^2),
\end{align*}
$$

and $\sum_{\substack{\sigma_x=\pm \\ \sigma_t=\pm}} \partial_h [(u_{\sigma_x}^{\sigma_t})^2](t - \sigma_t \tau, x - \sigma_x h) = 4\partial_x \mathbf{u}^2 + \mathcal{O}(\tau^2 + h^2)$. The $\mathcal{O}$ notation is understood in the $L^\infty$ norm over the common domain of definition. The signs $\sigma_t = \pm$ and $\sigma_x = \pm$ are arbitrary in the above expressions unless they are bound to a sum.

*Proof.* It is sufficient in this proof to assume that $\phi \in C^5$, rather than $\phi \in C^\infty$, hence theorem 1.6 may be slightly weakened. We follow the proof template of theorem 2.8, to which we refer for details. As before, the result follows from the consistency of the finite differences schemes $\partial_\tau$ and $\partial_h$, and the second order error terms are obtained when the expressions are invariant under the reflections $\tau \mapsto -\tau$ and $h \mapsto -h$. The time and space shifts in the definition of $u_{\sigma_x}^{\sigma_t}$ break these symmetries, but they are restored when considering averages of over the four possible sign combinations $\sigma_t = \pm$ and $\sigma_x = \pm$. Note also that $\rho > 0$ by assumption, hence $\rho > 0$ for sufficiently small $|h|$, and thus the ratio $u_{\sigma_x}^{\sigma_t}$ is well defined. $\square$
Inserting the expansions of theorem A.3 in (39) yields as desired

$$
\begin{align}
-\mathcal{L}(\hat{\phi}; \phi) &= 2h \sum_{x \in \mathcal{X}} \left( \mathbf{u}(\tau, x) - \nu \tau \partial_{xx} \mathbf{u}(\tau, x) + \frac{\tau}{2} \partial_x \mathbf{u}^2(\tau, x) - u_0(x) + \mathcal{O}(\tau^2 + h^2) \right) \hat{\phi}(0, x) \nonumber \\
&\quad + 4\tau h \sum_{t \in \mathcal{T}} \sum_{x \in \mathcal{X}} \left( \partial_t \mathbf{u}(t, x) - \nu \partial_{xx} \mathbf{u}(t, x) + \frac{1}{2} \partial_x \mathbf{u}^2(t, x) + \mathcal{O}(\tau^2 + h^2) \right) \hat{\phi}(t, x) \nonumber \\
&\le C_{\mathcal{L}}(\tau^2 + h^2) \left( \|\hat{\phi}(0)\|_1 + \tau \sum_{t \in \mathcal{T}} \|\hat{\phi}(t)\|_1 \right), \tag{40}
\end{align}
$$

where the constant $C_{\mathcal{L}}$ only depends on $\phi$.

## A.2 Lower estimate of the quadratic term

The second order term $\mathcal{Q}(\hat{\phi}; \phi)$ defined in (37) is a sum of perspective functions, whose first argument reads

$$r_{\sigma_x}^{\sigma_t}(t, x) := \hat{m}(t, x) - \nu \partial_h \hat{\rho}(t + \sigma \tau, x) - u_{\sigma_x}^{\sigma_t}(t, x) \hat{\rho}(t + \sigma_t \tau, x + \sigma_x h), \quad (41)$$

for all $t \in \mathcal{T}'$, $x \in \mathcal{X}$, and any signs $\sigma_t, \sigma_x = \pm$. Similarly to section 2.2, we regard $r_{\sigma_x}^{\sigma_t}$ as the r.h.s. of the discretization a time reversed linear PDE (43) below, closely related to (6); the key step of the proof, achieved in (44) below, is to estimate the unknown $\hat{\phi}$ in terms of these r.h.s. The time discretization is explicit if $\sigma_t = +$, and implicit if $\sigma_t = -$. Define

$$\partial_h^+ u(x) := \frac{u(x + 2h) - u(x)}{2h} = \partial_h u(x + h), \quad \partial_h^- u(x) := \frac{u(x) - u(x - 2h)}{2h} = \partial_h u(x - h), \quad (42)$$

the non-centered first order finite differences operator. We focus on the parameters $\sigma_t = \sigma_x = -$ (implicit scheme, upwind finite difference) which yield the most stable discretization, and rename for readability $u := u_-$ and $r := r_-$. Inserting the expressions $\hat{m} := \partial_\tau \hat{\phi}$ and $\hat{\rho} := -\partial_h \hat{\phi}$ in (41) and rearranging terms we obtain (omitting the space variable $x$ for readability)

$$\hat{\phi}(t - \tau) - 2\tau \nu \partial_{hh} \hat{\phi}(t - \tau) + 2\tau u(t) \partial_h^- \hat{\phi}(t - \tau) = \hat{\phi}(t + \tau) - 2\tau r(t), \quad \hat{\phi}(T) = 0. \quad (43)$$

**Proposition A.4.** For any $u \in [0, \infty]^{\mathcal{X}}$ and any $\nu \ge 0$, the operator $A := \text{Id} + u\partial_h^- - \nu\partial_{hh}$ is inverse positive. If in addition $\|\partial_h u\|_{\infty} < 1$, then $\|A^{-1}f\|_1 \le (1 - \|\partial_h u\|_{\infty})^{-1}\|f\|_1$, $\forall f \in \mathbb{R}^{\mathcal{X}}$.

*Proof.* We proceed similarly to the proof of theorem 2.9, and observe that the matrix of the finite differences operator $A := \text{Id} + u\partial_h^- - \nu\partial_{hh}$ has non-positive off-diagonal entries, and its rows sum to one. This fits the condition of an M-matrix, thus an inverse positive matrix, as announced. Now consider $f \ge 0$, and let $g := A^{-1}f$, which is thus non-negative. One has

$$\frac{\|f\|_1}{2h} = \sum_{x \in \mathcal{X}} f(x) = \sum_{x \in \mathcal{X}} [g(x) + u(x)\partial_h^- g(x) - \nu\partial_{hh}g(x)] = \sum_{x \in \mathcal{X}} g(x)(1 - \partial_h^+ u(x)) \ge \frac{\|g\|_1}{2h}(1 - \|\partial_h u\|_{\infty}),$$

using the adjointness (38) of $\partial_h^-$ and $-\partial_h^+$, the fact that $\partial_{hh}$ is self-adjoint and vanishes on constants, and the observation that $\|\partial_h^+ u\|_{\infty} = \|\partial_h u\|_{\infty}$ by (42). We conclude, similarly to theorem 2.9, by splitting an arbitrary $f \in \mathbb{R}^{\mathcal{X}}$ into a positive and a negative part and using the linearity of $A^{-1}$. $\square$
In the rest of this section, we work under theorem 1.6, in such way that $\phi \in C^\infty([0, T] \times \mathbb{T})$ is a smooth function, and the Burgers' equation solution $u$ is smooth and positive over the time interval of interest. Recalling that $u_{\sigma_t}^{\sigma_x} = u + O(\tau + h)$ by theorem A.3, we choose $\tau_* > 0$ and $h_* > 0$ such that $u_{\sigma_t}^{\sigma_x} \ge 0$ on $[\tau, T - \tau] \times \mathbb{T}$ for all $\tau \in [0, \tau_*]$ and $h \in [0, h_*]$, and thus in particular $u := u_- \ge 0$. Using (43) and theorem A.4 we obtain

$$ \|\hat{\phi}(T - 2\tau N)\|_1 \le 2\tau \sum_{0 \le n < N} \|r(T - (2n + 1)\tau)\|_1 (1 - \tau K)^{N-n}, \quad K := 2 \max_{t \in \mathcal{T}'} \|\partial_h u(t)\|_{\infty}, \quad (44) $$

for all $0 \le N \le N_\tau := T/(2\tau)$, by an immediate induction argument and provided $\tau K < 1$. By theorem A.3, one has $\partial_h u(t) = \partial_x u + O(\tau + h)$, hence $K$ is bounded independently of $\tau$ and $h$. From this point, the proof is identical to section 2.2, with the sole exception of theorem 2.11 which needs to be replaced with the following counterpart, whose proof is trivial hence omitted.

**Lemma A.5.** Let $\phi \in \mathbb{R}^\mathcal{X}$ and $\rho := 1 - \partial_h \phi$. If $\rho \ge 0$ then $\|\rho\|_1 = 1$.

## B Arithmetic mean based scheme

We consider in this appendix a variant $\mathcal{E}_{\tau h}^\theta$, where $\theta \in [0, 1]$, of the energy (13) associated with the discretization of the QPME. It is defined as follows, with the notations of theorem 1.2

$$ \mathcal{E}_{\tau h}^{\theta}(m, \rho) := 2\tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}'} \frac{m(t, x)^2}{2\mathcal{I}_{\tau}^{\theta}\rho(t, x)} - m(t, x)u_0(x), \quad (45) $$

$$ \text{where } \mathcal{I}_{\tau}^{\theta} \rho(t, x) := (1 - \theta)\rho(t - \tau, x) + \theta\rho(t + \tau, x). \quad (46) $$

The ratio in (45) should be understood as an instance of the perspective function (10), and $\mathcal{I}_{\tau}^{\theta}$ thus denotes a linear time interpolation operator. By convention, $\mathcal{E}_{\tau h}^{\theta}(\hat{\phi}) := \mathcal{E}_{\tau h}^{\theta}(\partial_{\tau} \phi, 1 + L_h \phi)$, for all $\phi \in \Phi_{\tau h}^P$. For readability, we denote $\mathcal{E}^{\theta} := \mathcal{E}_{\tau h}^{\theta}$ and $\mathcal{I}^{\theta} := \mathcal{I}_{\tau}^{\theta}$ in the following, and we also adopt theorem 2.2: $\Phi := \Phi_{\tau h}^P$, $\mathcal{T} := \mathcal{T}_{\tau}$, $\mathcal{T}' := \mathcal{T}'_{\tau}$, $\mathcal{X} := \mathbb{T}_h^d$, etc.

Using the exact expansion of the perspective function, see theorem 2.1, we obtain $\mathcal{E}^{\theta}(\hat{\phi} + \hat{\phi}) = \mathcal{E}^{\theta}(\phi) + \mathcal{L}^{\theta}(\hat{\phi}; \phi) + \mathcal{Q}^{\theta}(\hat{\phi}; \phi)$ for any $\phi, \hat{\phi} \in \Phi$ such that $\mathcal{E}^{\theta}(\phi) < \infty$, where

$$ \mathcal{L}^{\theta}(\hat{\phi}; \phi) := 2\tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}'} \left( u(t, x)\hat{m}(x) - \frac{1}{2}u(t, x)^2 \mathcal{I}^{\theta}\hat{\rho}(t, x) \right) + (2h)^d \sum_{x \in \mathcal{X}} u_0(x)\hat{\phi}(0, x), $$

$$ \mathcal{Q}^{\theta}(\hat{\phi}; \phi) := 2\tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}'} \mathcal{P}(\hat{m}(t, x) - u(t, x)\mathcal{I}^{\theta}\hat{\rho}(t, x), \mathcal{I}^{\theta}\rho(t, x) + \mathcal{I}^{\theta}\hat{\rho}(t, x)), $$

$$ \text{with } m := \partial_{\tau} \phi, \rho := 1 + L_h \phi, \hat{m} := \partial_{\tau} \hat{\phi}, \hat{\rho} := L_h \hat{\phi}, \text{ and } u(t, x) := m(t, x)/\mathcal{I}^{\theta} \rho(t, x). $$

The latter denominator $\mathcal{I}^{\theta}\rho(t, x)$ is positive for sufficiently small $h > 0$ under theorem 1.3. The arithmetic interpolation operator w.r.t. time satisfies

$$ \sum_{t \in \mathcal{T}'} \mathcal{I}_{\tau}^{\theta} f(t) g(t) = (1 - \theta)f(0)g(\tau) + \theta f(T)g(T - \tau) + \sum_{t \in \mathcal{T}} f(t) \mathcal{I}_{\tau}^{1-\theta} g(t), $$

for any centered $f: \mathcal{T} \to \mathbb{R}$ and staggered $g: \mathcal{T}' \to \mathbb{R}$, recalling that $\mathcal{T}' := \mathcal{T} \setminus \{0, T\}$. In
combination with the summation by parts formulas of theorems 2.4 and 2.5 we obtain

$$
\begin{align}
-\mathcal{L}^{\theta}(\hat{\phi}; \phi) &= (2h)^d \sum_{x \in \mathcal{X}} \left( u(\tau, x)\hat{\phi}(0, x) + \tau(1-\theta) u(\tau, x)^2 \hat{\rho}(0, x) - u_0(x)\hat{\phi}(0, x) \right) \nonumber \\
&\quad + 2\tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}} \left( \partial_{\tau}u(t, x)\hat{\phi}(t, x) + \frac{1}{2}\hat{\rho}(t, x) \mathcal{I}^{1-\theta}u(t, x)^2 \right) \nonumber \\
&= (2h)^d \sum_{x \in \mathcal{X}} \left( u(\tau, x) + \tau(1-\theta)L_h u^2(\tau, x) - u_0(x) \right) \hat{\phi}(0, x) \tag{47} \\
&\quad + 2\tau(2h)^d \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}} \left( \partial_{\tau}u(t, x) + \frac{1}{2}\mathcal{I}^{1-\theta}L_h u^2(t, x) \right) \hat{\phi}(t, x). \tag{48}
\end{align}
$$

The map $\phi \in \Phi$ is a critical point of $\mathcal{E}^\theta$ iff all the coefficients of $\hat{\phi}$ vanish in (47) and (48), i.e.

$$
u(\tau) + \tau(1 - \theta)L_h u^2(\tau) = u_0, \qquad \partial_{\tau}u(t) + \frac{1}{2}\mathcal{I}^{1-\theta}L_h u^2(t) = 0, \forall t \in \mathcal{T}. \quad (49)
$$

This amounts, it turns out, to an implicit ($\theta = 0$), explicit ($\theta = 1$), or second order midpoint ($\theta = 1/2$) discretization of the QPME (1). The numerical solution $u$ can thus be computed one timestep at a time, and there is no point in minimizing the energy (45) using a global space-time method. In contrast, the coefficient of $\hat{\phi}(t,x)$ is more complex in the proposed discretization (21), and the corresponding criticality condition can be reformulated as

$$
\partial_{\tau} \tilde{u}(t) + \frac{1}{4} L_h \left( \tilde{u}(t - \tau)^2 \frac{\mathcal{H} \rho(t)^2}{\rho(t - \tau)^2} \right) + \frac{1}{4} L_h \left( \tilde{u}(t + \tau)^2 \frac{\mathcal{H} \rho(t)^2}{\rho(t + \tau)^2} \right) = 0, \forall t \in \mathcal{T}, \quad (50)
$$

where $\mathcal{H}\rho(t) := 2(\rho(t - \tau)^{-1} + \rho(t + \tau)^{-1})^{-1}$ denotes the harmonic interpolation of densities w.r.t. time, and $\tilde{u} := m/\mathcal{H}\rho$. If $\rho$ varies slowly w.r.t. the time variable $t$, then (50) is close in spirit with (49, right). However, in contrast with (49), it does not seem possible to solve (50) successively one timestep at a time, especially in view of the terminal boundary condition $\rho(T) = 1$, and of the initial condition for $u(\tau)$ derived from (20). This problem formulation is reminiscent of mean field games, see the discussion in [Bre20, Chapter 5].

Turning our attention to the second order term $\mathcal{Q}^\theta$, we see that the first argument

$$
r(t) := \hat{m}(t) - u(t) \mathcal{I}^{\theta} \hat{\rho}(t) = \frac{\phi(t + \tau) - \phi(t - \tau)}{2\tau} - (1 - \theta) u(t) L_h \phi(t - \tau) - \theta u(t) L_h \phi(t + \tau),
$$

can be regarded as the r.h.s. of a reversed in time discretized linear PDE closely related with (6), recall the definition (46) of the time interpolation operator $\mathcal{I}^\theta$. Equivalently

$$
\phi(t - \tau) + 2\tau(1 - \theta)u(t)L_h\phi(t - \tau) = \phi(t + \tau) - \theta u(t)L_h\phi(t + \tau) - r(t), \quad \phi(T) = 0. \quad (51)
$$

The scheme is implicit if $\theta = 0$, explicit if $\theta = 1$, and second order semi-implicit if $\theta = 1/2$. The proof template of section 2.2 requires the scheme (51) to be stable in the $\ell^1(\mathcal{X})$ norm. Choosing $\theta = 0$, and recalling from theorem 2.9 that the implicit scheme (51) is indeed $\ell^1$ stable, we can obtain a convergence result for the implicit scheme (49) for the QPME (convergence of the dual potential $\phi$, similarly to theorem 1.4 but with a first order rate $\mathcal{O}(\tau + h^2)$). This approach unfortunately does not extend to the explicit scheme or to the midpoint scheme (51), which are unstable in the $\ell^1(\mathcal{X})$ norm unless a restrictive Courant-Friedrichs-Levy condition is enforced (even though the midpoint scheme is stable in the $\ell^2(\mathcal{X})$ norm without CFL).
C Uniform convergence of the primal variable, by mollification

Our convergence results theorems 1.4 and 1.7 state that the proposed BBB numerical method produces<fnref n="1" /> an approximation $\phi_{\tau,h}$ of the dual potential $\phi$, in the rather weak $L^\infty([0, T], L^1(\mathbb{T}^d))$ norm, with a second order $O(\tau^2+h^2)$ convergence rate w.r.t. the half timestep $\tau$ and half gridscale $h$. Using a standard regularization by convolution w.r.t. the spatial variable, we construct in this appendix an approximation of the primal variable $u$ of the PDE of interest, the QPME or Burgers' viscous equation, which converges in the $L^\infty([0, T] \times \mathbb{T}^d)$ norm, but with a slower $O(h^{\frac{2}{d+2}})$ convergence rate, when $\tau = h$. This convolution step is not implemented in our numerical experiments, since it does not appear to be necessary in practice and only degrades the results.

**Definition C.1 (Spatial convolution).** A mollifier is a non-negative, radially symmetric mapping $\rho \in C^2([0, T] \times \mathbb{T}^d)$ with compact support and unit integral. Spatial mollification is defined as

$$[\phi \star \rho_{\varepsilon}](t, x) := \int_{\mathbb{R}^d} \phi(t, x - y)\rho_{\varepsilon}(y)dy, \quad \text{where } \rho_{\varepsilon}(y) := \frac{1}{\varepsilon^d}\rho\left(\frac{y}{\varepsilon}\right), \text{ with } \varepsilon > 0.$$

Given $\phi \in C^2([0, T] \times \mathbb{T}^d)$, one obtains using the symmetry of $\rho$ and a Taylor expansion

$$|\phi \star \rho_{\varepsilon}(t, x) - \phi(t, x)| \le \int_{\mathbb{R}^d} |\phi(t, x - y) - \phi(t, x) + y^{\top}\nabla\phi(t, x)|\rho_{\varepsilon}(y)dy \le \varepsilon^2 C(\rho)\|\nabla^2\phi\|_{\infty},$$

where $\nabla$ and $\nabla^2$ denote the gradient and Hessian operators w.r.t. the space variable $x$, and $C(\rho) := \frac{1}{2} \int_{\mathbb{R}^d} \rho(x)\|x\|^2dx$.

**Proposition C.2.** Let $\phi \in C^4([0, T] \times \mathbb{T}^d)$, and let $\phi_{\tau,h} \in L^\infty([0, 1], L^1(\mathbb{T}^d))$ be such that $\sup_{t \in [0, T]} \|\phi(t, \cdot) - \phi_{\tau,h}(t, \cdot)\|_{L^1(\mathbb{T}^d)} = O(\tau^2 + h^2)$. Then with the notations of theorem C.1

$$\|\nabla^2\phi - \nabla_h^2(\phi_{\tau,h} \star \rho_{\varepsilon})\|_{L^\infty} = O\left(\varepsilon^2 + h^2 + \frac{\tau^2 + h^2}{\varepsilon^{2+d}}\right), \quad (52)$$

$$\|\nabla\phi - \nabla_h(\phi_{\tau,h} \star \rho_{\varepsilon})\|_{L^\infty} = O\left(\varepsilon^2 + h^2 + \frac{\tau^2 + h^2}{\varepsilon^{1+d}}\right), \quad (53)$$

$$\|\partial_t\phi - \partial_\tau(\phi_{\tau,h} \star \rho_{\varepsilon})\|_{L^\infty} = O\left(\varepsilon^2 + \tau^2 + \frac{\tau^2 + h^2}{\tau\varepsilon^d}\right), \quad (54)$$

for all $\tau, h, \varepsilon > 0$, where the $L^\infty$ norm is over the common domain of definition, and where $\nabla_h$ and $\nabla_h^2$ denote the standard centered finite difference schemes for the gradient and Laplacian.

*Proof.* Note that $\nabla_h^2(\phi_{\tau,h} \star \rho_{\varepsilon}) = (\nabla_h^2\phi_{\tau,h}) \star \rho_{\varepsilon} = \phi_{\tau,h} \star \nabla_h^2(\rho_{\varepsilon})$, and likewise $\partial_{\tau}(\phi_{\tau,h} \star \rho_{\varepsilon}) = (\partial_{\tau}\phi_{\tau,h}) \star \rho_{\varepsilon}$. Focusing first on (52), we estimate: $\|\nabla^2\phi - \nabla_h^2(\phi_{\tau,h} \star \rho_{\varepsilon})\|_{L^\infty}$

$$\begin{aligned} &\le \|\nabla^2\phi - (\nabla^2\phi) \star \rho_{\varepsilon}\|_{L^\infty} + \| (\nabla^2\phi - \nabla_h^2\phi) \star \rho_{\varepsilon}\|_{L^\infty} + \| (\phi - \phi_{\tau,h}) \star (\nabla_h^2\rho_{\varepsilon})\|_{L^\infty} \\ &\le C\varepsilon^2 \|\nabla^4\phi\|_{L^\infty} + \|\nabla^2\phi - \nabla_h^2\phi\|_{L^\infty} \|\rho_{\varepsilon}\|_{L^1(\mathbb{R}^d)} + \|\phi - \phi_{\tau,h}\|_{L^\infty([0,T], L^1(\mathbb{T}^d))} \|\nabla_h^2\rho_{\varepsilon}\|_{L^\infty(\mathbb{R}^d)}. \end{aligned}$$

We conclude the proof of (52) using the second order consistency of the centered finite differences Hessian $\nabla_h^2$, and the identities $\|\rho_{\varepsilon}\|_{L^1} = \|\rho\|_{L^1} = 1$ and $\|\nabla_h^2\rho_{\varepsilon}\|_{L^\infty} \le C\|\nabla^2\rho_{\varepsilon}\|_{L^\infty} =$

<footnote n="1">¹Precisely, the mapping $\phi_{\tau,h} : [0, T] \times \mathbb{T}^d \to \mathbb{R}$ can be obtained as the extension by multilinear interpolation of the minimizer $\phi_{\tau,h}^0 : \mathcal{T}_{\tau} \times \mathbb{T}^d \to \mathbb{R}$ of the discrete energy.</footnote>
$C\varepsilon^{-(2+d)}\|\nabla^2\rho\|_{\infty}$, for some absolute constant $C$. The first order case (53) is similar. Considering the time finite difference (54) we compute

$$
\begin{aligned}
\|\partial_t\phi - \partial_\tau(\phi_{\tau,h} \star \rho_\varepsilon)\|_{L^\infty} &\le \|\partial_t\phi - \partial_\tau\phi\|_{L^\infty} + \|\partial_\tau\phi - (\partial_\tau\phi) \star \rho_\varepsilon\|_{L^\infty} + \|[\partial_\tau(\phi - \phi_{\tau,h})] \star \rho_\varepsilon\|_{L^\infty} \\
&\le \|\partial_t\phi - \partial_\tau\phi\|_{L^\infty} + \varepsilon^2 C(\rho) \|\nabla^2\partial_\tau\phi\|_{L^\infty} + \frac{1}{\tau} \|\phi - \phi_{\tau,h}\|_{L^\infty([0,T], L^1(\mathbb{T}^d))} \|\rho_\varepsilon\|_{L^\infty(\mathbb{R}^d)},
\end{aligned}
$$

and (54) follows using the second order consistency of the staggered finite difference $\partial_\tau$ and the identity $\|\rho_\varepsilon\|_{L^\infty} = \varepsilon^{-d}\|\rho\|_{\infty}$. $\square$

Letting $\phi_h := \phi_{\tau,h} \star \rho_\varepsilon$ with $\tau := h$ and $\varepsilon := h^{\frac{1}{2+d}}$, we obtain in the QPME case

$$ \|u - u_h\|_{\infty} = \mathcal{O}(h^{\frac{2}{d+2}}), \quad \text{with } u := \frac{\partial_t \phi}{1 - \operatorname{div}(\mathcal{D} \nabla \phi)}, \quad \text{and } u_h := \frac{\partial_\tau \phi_h}{1 - \operatorname{Tr}(\mathcal{D} \nabla_h^2 \phi_h) - \operatorname{div}(\mathcal{D})^\top \nabla_h \phi_h}. $$

Likewise, in Burgers' case, with $\tau = h$ and $\varepsilon = h^{\frac{1}{3}}$ (since $d = 1$)

$$ \|u - u_h\|_{\infty} = \mathcal{O}(h^{\frac{2}{3}}), \quad \text{with } u = \frac{\partial_t \phi + \nu \partial_{xx} \phi}{1 - \partial_x \phi}, \quad u_h = \frac{\partial_\tau \phi_h + \nu \partial_{hh} \phi_h}{1 - \partial_h \phi_h}. $$

## D Closed form of the Barenblatt Potential

The Barenblatt profile is a closed form solution of the QPME in a homogeneous medium, which is self-similar in time, and compactly supported in space for each time, see (36). We describe in this section the corresponding dual variables $m, \rho, \phi$ in the BBB formulation of the QPME, on some time interval $[T_0, T]$. Note that $m, \rho, \phi$ depend on the final time $T$, through the terminal constraints $\rho(T) = 1$ and $\phi(T) = 0$, and are thus not self-similar in time. Interestingly, $x \in \mathbb{R}^d \mapsto \rho(t, x)$ has discontinuities across two interfaces: $|x| = R(t)$ and $|x| = R_T$, see below.

**Proposition D.1.** Consider the Barenblatt profile $u(t, x) := \frac{2}{t^\alpha} \max\{0, \gamma - \frac{\beta|x|^2}{4t^{2\beta}}\}$, where $t > 0$, $x \in \mathbb{R}^d$, $\alpha := d/(d+2)$, $\beta := 1/(d+2)$, and $\gamma > 0$ is an arbitrary parameter. Fix a compact time interval $[T_0, T] \ni t$, where $T > T_0 > 0$, and define

$$ R(t) := \nu t^\beta, \quad \text{with } \nu := 2\sqrt{(d+2)\gamma}. $$

By construction, $\text{supp}(u(t, \cdot)) = \bar{B}(0, R(t))$. Denote by $R_0 := R(T_0)$ and $R_T := R(T)$ the radius of the support of the Barenblatt solution at the initial and final times. For $T_0 < t < T$ we set

$$ \phi(t, x) := \begin{cases} \frac{2}{T^\alpha} \gamma (t-T) + \frac{|x|^2}{2d} \left(1 - \left(\frac{t}{T}\right)^\alpha\right) & \text{if } |x| < R(t), \\ -2\gamma T^{1-\alpha} + \frac{|x|^2}{2d} - \frac{1}{d(d+2)R_T^d} |x|^{d+2} & \text{if } R(t) < |x| < R_T, \\ 0 & \text{if } R_T < |x|. \end{cases} \quad (55) $$

Then $\phi$ is continuous on $[T_0, T] \times \mathbb{R}^d$, and is $C^\infty$ outside $\Gamma := \{(t, x) \in [T_0, T] \times \mathbb{R}^d \mid |x| = R(t) \text{ or } |x| = R_T\}$. The distributional derivative $m := \partial_t\phi$ is a continuous function, whereas the distributional derivative $\rho := 1 - \Delta\phi$ is a bounded and piecewise continuous function with jumps across $\Gamma$. One has the closed form expressions

$$ m(t, x) = \frac{2}{T^\alpha} \max \left\{ 0, \gamma - \frac{\beta |x|^2}{4t^{2\beta}} \right\}, \quad \rho(t, x) = \begin{cases} \left(\frac{t}{T}\right)^\alpha, & \text{if } |x| < R(t), \\ 2\left(\frac{|x|}{R_T}\right)^d, & \text{if } R(t) < |x| < R_T, \\ 1, & \text{otherwise.} \end{cases} \quad (56) $$
Moreover one has $\partial_t \rho = \Delta m$ in the sense of distributions, $\rho(T, \cdot) = 0$, and denoting $u_0 := u(T_0, \cdot)$ one has

$$ \int_{[T_0, T] \times \mathbb{R}^d} \frac{u^2}{2} = \int_{[T_0, T] \times \mathbb{R}^d} \left( -\frac{m^2}{2\rho} + m u_0 \right), \quad (57) $$

showing that there is no duality gap in (4), and therefore $(\rho, m)$ solves the ballistic Benamou-Brenier formulation of the QPME. Finally, we observe that $m = \rho u$, as expected from (6).

**Proof.** Since $\phi$ has a radially symmetric profile, we express it as $\phi(t, x) = \varphi(t, |x|)$, where $\varphi: [T_0, T] \times [0, \infty[ \to \mathbb{R}$ is piecewise polynomial. Our first objective is to show that $\varphi$ and $\partial_r\varphi$ are continuous across the interfaces $\{r = R(t)\}$ and $\{r = R_T\}$; once this is proved, it follows that the distributional derivatives $m(t, x) = \partial_t\varphi(t, |x|)$ and $\rho(t, x) = 1 - r^{-(d-1)}\partial_r(r^{d-1}\partial_r\varphi)(t, |x|)$ have no singular part (i.e. jump term), hence can be evaluated piecewise by direct differentiation, which yields (56). Denoting by $[\varphi(t, \cdot)]_r := \lim_{\epsilon \to 0^+} \varphi(t, r + \epsilon) - \varphi(t, r - \epsilon)$ the jump (if any) at $r$ we compute

$$ \begin{aligned} -[\varphi(t, \cdot)]_{R(t)} &= \frac{2}{T^\alpha} \gamma t - \frac{R(t)^2}{2d} \left(\frac{t}{T}\right)^\alpha + \frac{1}{d(d+2)R_T^d} R(t)^{d+2} = \frac{2\gamma t}{T^\alpha} - \frac{\nu^2 t}{2dT^\alpha} + \frac{\nu^2 t}{d(d+2)T^\alpha} = 0, \\ -[\partial_r\varphi(t, \cdot)]_{R(t)} &= -\frac{R(t)}{d} \left(\frac{t}{T}\right)^\alpha + \frac{R(t)^{d+1}}{dR_T^d} = \frac{R(t)}{d} \left(\left(\frac{R(t)}{R_T}\right)^d - \left(\frac{t}{T}\right)^\alpha\right) = 0, \\ [\varphi(t, \cdot)]_{R_T} &= 2\gamma T^{1-\alpha} - \frac{R_T^2}{2d} + \frac{R_T^{d+2}}{d(d+2)R_T^d} = T^{1-\alpha} \left(2\gamma - \frac{4\gamma(d+2)}{2d} + \frac{4\gamma(d+2)}{d(d+2)}\right) = 0, \end{aligned} $$

and finally $[\partial_r\varphi(t, \cdot)]_{R_T} = -R_T/d + R_T^{d+1}/(dR_T^d) = 0$. We simply substituted the definition of $R(t)$ and $\nu$, and used the relations $d\beta = \alpha$ and $1 - \alpha = 2\beta$ relating the exponents. The identities $\phi(T, \cdot) = 0$ and $m = \rho u$ are clear by construction, and $\partial_t\rho = -\partial_t\Delta\phi = -\Delta m$.

We now turn to the second part of the proof, and obtain for any $x \in \mathbb{R}^d$, using $m = \rho u = \partial_t\phi$

$$ \int_{T_0}^{T} \frac{m(t, x)^2}{\rho(t, x)} - m(t, x)u_0(x)dt = \int_{T_0}^{T} m(t, x)(u(t, x) - u_0(x))dt = - \int_{T_0}^{T} \phi(t, x)\partial_t u(t, x)dt. \quad (58) $$

The integration by parts is justified since $\partial_t\phi = m$ is continuous by (56), and $\partial_t u(\cdot, x) \in L_{\text{loc}}^\infty$ is locally bounded on $[T_0, T]$ in view of the explicit expression of the Barenblatt profile (36).

Likewise, for any fixed $t \in [0, T]$, we obtain omitting for readability the arguments $(t, x)$

$$ 2 \int_{\mathbb{R}^d} \phi \partial_t u = \int_{B(t)} \phi \Delta(u^2) = - \oint_{\partial B(t)} u^2 \langle \nabla \phi, n \rangle + \int_{B(t)} u^2 \Delta \phi = \int_{B(t)} u^2 (1 - \rho) = \int_{\mathbb{R}^d} \left( u^2 - \frac{m^2}{\rho} \right), $$

where the open unit ball $B(t) := B(0, R(t))$ is the interior of the support of $u(t, \cdot)$, and $n(x) := x/\|x\|$ denotes the outward unit vector. We used the QPME $2\partial_t u = \Delta(u^2)$, the fact that $u$ vanishes on $\partial B(0, R(t))$, the defining identity $1 - \rho = \Delta\phi$, and the equality $m = \rho u$. The integration by parts is justified since $u(t, \cdot)$ and $\phi(t, \cdot)$ are $C^\infty$ on the open ball $B(0, R(t))$, and $u$ and $\nabla\phi(t, x) = \partial_r\varphi(t, |x|)n(x)$ are continuous on the closed ball $\bar{B}(0, R(t))$, see above.

Integrating the identity (58) w.r.t. $x \in \mathbb{R}^d$, and the last equation w.r.t. $t \in [T_0, T]$, we obtain (57) by linear combination, which concludes the proof. $\square$

## E Validity of the BBB formulation of Burgers' viscous equation

We show that the BBB formulation of Burgers' viscous equation is well posed on any time interval $[0, T]$, by a straightforward adaptation of [Bre20, Proposition 5.2.1] which is devoted to
the closely related Hamilton-Jacobi equation $\partial_t u + \frac{1}{2}|\nabla u|^2 = \nu \Delta u$. In numerical applications one may nevertheless need to restrict the time $T$, for reasons discussed below.

**Proposition E.1.** Given a smooth initial condition $u_0 : \mathbb{T} \to \mathbb{R}$, and a positive viscosity coefficient $\nu > 0$, consider the solution $u : [0, T] \times \mathbb{T} \to \mathbb{R}$ to the viscous Burgers' equation $\partial_t u + \partial_x \frac{u^2}{2} = \nu \Delta u$. The BBB formulation (8) of this PDE admits a solution $(\rho, m)$, obeying

$$ \int_{[0,T]\times\mathbb{T}} \frac{u^2}{2} = \int_{[0,T]\times\mathbb{T}} \left( -\frac{(m - \nu \partial_x \rho)^2}{2\rho} + m u_0 \right), \qquad (59) $$

showing that there is no duality gap in (4). Moreover, $m = \rho u + \nu \partial_x \rho$ as expected from (6).

*Proof.* Let us define $\rho$ as the solution of the following backward in time Cauchy problem, which is positive by general properties of transport diffusion equations

$$ \partial_t \rho + \partial_x (\rho u) = -\nu \Delta \rho, \qquad \rho(T, \cdot) = 1. \qquad (60) $$

We then set $m := \rho u + \nu \partial_x \rho$ and $\phi(t, x) := -\int_t^T m(s, x)ds$. By construction, there holds

$$ \partial_t \rho + \partial_x m = 0, \qquad u = \frac{m - \nu \partial_x \rho}{\rho}, \qquad \partial_x \phi = 1 - \rho, \qquad \partial_t \phi = m. $$

By integration by parts, successively in time and space, we obtain

$$ \int_{[0,T]\times\mathbb{T}} m(u - u_0) = - \int_{[0,T]\times\mathbb{T}} \phi \partial_t u = \int_{[0,T]\times\mathbb{T}} \left[ \phi \partial_x \frac{u^2}{2} - \nu \phi \partial_{xx} u \right] = \int_{[0,T]\times\mathbb{T}} \left[ (\rho - 1) \frac{u^2}{2} + \nu u \partial_x \rho \right] $$

Using the defining expression of $m$, and the above identity, we conclude that

$$ \int_{[0,T]\times\mathbb{T}} \rho u^2 = \int_{[0,T]\times\mathbb{T}} (m - \nu \partial_x \rho) u = \int_{[0,T]\times\mathbb{T}} [m u_0 + (\rho - 1) \frac{u^2}{2}]. $$

which is equivalent to (59), given that $\rho u^2 = (m - \nu \partial_x \rho)^2 / \rho$. This shows that there is no duality gap in (4), and therefore that $(\rho, m)$ solves the ballistic Benamou-Brenier formulation. $\square$

The positivity of the solution $\rho$ of the transport diffusion equation (60) is the key ingredient of the proof of theorem E.1. We describe below a heuristic estimate of $\min(\rho)$, which is interesting for numerical considerations, since our BBB implementation is not expected to reliably work if this quantity falls below machine floating point precision. Denoting by $x(t) \in \mathbb{T}$ the minimizer of $\rho(t, \cdot)$, in such way that $\partial_x \rho(t, x(t)) = 0$ and $\partial_{xx} \rho(t, x(t)) \ge 0$, we obtain formally

$$ \frac{d}{dt}[\rho(t, x(t))] = \partial_t \rho(t, x(t)) = -\partial_x(\rho u)(t, x(t)) - \nu \Delta \rho(t, x(t)) \le -\rho(t, x(t))\partial_x u(t, x(t)), $$

using the envelope theorem for the first identity, and (60) for the second one. Therefore, $\rho(t, x(t)) \ge \exp(\int_t^T \partial_x u(s, x(s))ds)$ in view of the terminal boundary condition $\rho(T) = 1$. From a dimensional analysis, one typically expects that $\|\partial_x u\|_{\infty} \lesssim (\max(u) - \min(u))^2 / (8\nu)$, where the factor $1/8$ is introduced to match the steadily propagating wave example $u(t, x) = 2c/(1 + \exp(c(x - ct)/\nu))$ with speed $c > 0$. Note also that $\max(u) = \max(u_0)$, and $\min(u) = \min(u_0)$. Therefore, we heuristically expect $\min(\rho) \gtrsim \exp(-T/T_*)$ with $T_* := 8\nu / (\max(u_0) - \min(u_0))^2$.