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

HAL
open science

# Minimum Flow Steepest Descent Approach for Nonlinear PDE

Noureddine Igbida

► To cite this version:

Noureddine Igbida. Minimum Flow Steepest Descent Approach for Nonlinear PDE. 2024. hal-04436356

HAL Id: hal-04436356

https://hal.science/hal-04436356v1

Preprint submitted on 3 Feb 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Minimum Flow Steepest Descent Approach for
Nonlinear PDE

Noureddine Igbida *

February 3, 2024

Abstract

This paper presents a minimum flow approach applicable to a wide range of doubly nonlinear diffusion problems. We introduce a minimum flow steepest descent algorithm that seeks an optimal traffic flow by minimizing an internal energy function, while incorporating a specific minimum flow constraint for the transition work. This flexible framework surpasses traditional methods by generalizing (among other things) the established $H^{-1}$-theory for linear diffusion to the nonlinear setting. It offers distinct advantages in handling diverse applications, leveraging both the intrinsic internal energy and the inherent traffic flow concepts. The approach demonstrably tackles various applications, including fluid flow in porous media and many other scenarios like the Stefan/Hele-Shaw problem. Notably, it can handle diverse differential operators, encompassing even nonlinear ones like the $p$-Laplacian and Leray-Lions operators. Furthermore, it allows for the simultaneous handling convection/transport phenomena and complex boundary conditions within both the energy and transition work functions, contrasting favorably with other steepest descent algorithm like the JKO scheme in Wasserstein spaces. The paper delves into details of this comparison and offers additional theoretical insights involving $p'$-curve and metric gradient flow in Sobolev dual spaces.

Contents

<table><tr><td>1</td><td>Introduction and preliminaries</td><td>2</td></tr><tr><td>1.1</td><td>Introduction</td><td>2</td></tr><tr><td>1.2</td><td>Preliminaries</td><td>4</td></tr><tr><td>2</td><td>Plan of the paper</td><td>6</td></tr><tr><td>3</td><td>Overview : Wasserstein distance vs $H^{-1}$-norm</td><td>6</td></tr><tr><td>3.1</td><td>Proximity through $W_2$-Wasserstein distance</td><td>6</td></tr><tr><td>3.2</td><td>Proximity through $H^{-1}$-norm</td><td>7</td></tr><tr><td>3.3</td><td>Minimum flow concept and related steepest descent</td><td>9</td></tr><tr><td>4</td><td>Notations and assumptions</td><td>11</td></tr><tr><td>4.1</td><td>Reminder</td><td>11</td></tr><tr><td>4.2</td><td>Assumptions</td><td>13</td></tr></table>

<footnote>*Institut de recherche XLIM, UMR-CNRS 7252, Faculté des Sciences et Techniques, Université de Limoges, France.
Email: noureddine.igbida@unilim.fr</footnote>
<table>
  <tbody>
    <tr>
      <td><b>5</b></td>
      <td><b>Proximal energy minimization</b></td>
      <td><b>14</b></td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Main results</td>
      <td>14</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Proofs of Theorem 1 and Theorem 2</td>
      <td>18</td>
    </tr>
    <tr>
      <td><b>6</b></td>
      <td><b>Steepest descent algorithm</b></td>
      <td><b>24</b></td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Main results</td>
      <td>26</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>The case where V ≡ 0</td>
      <td>26</td>
    </tr>
    <tr>
      <td>6.1.2</td>
      <td>Diffusion-Transport problem with linear diffusion</td>
      <td>26</td>
    </tr>
    <tr>
      <td>6.1.3</td>
      <td>Diffusion-Transport problem with nonlinear diffusion</td>
      <td>27</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Proofs of Theorem 3, Theorem 4 and Theorem 5</td>
      <td>27</td>
    </tr>
    <tr>
      <td><b>7</b></td>
      <td><b>Comments, remarks and extensions</b></td>
      <td><b>32</b></td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>W<sub>Γ<sub>D</sub></sub><sup>-1,p'</sup>(Ω)-Gradient flow and p'-curves of maximal slope</td>
      <td>32</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Prediction-correction algorithm</td>
      <td>35</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Uniqueness</td>
      <td>36</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Linear growth cost F : granular diffusion</td>
      <td>37</td>
    </tr>
    <tr>
      <td>7.5</td>
      <td>Linear growth cost F*: total variation flow</td>
      <td>37</td>
    </tr>
    <tr>
      <td>7.6</td>
      <td>Numerical computation</td>
      <td>37</td>
    </tr>
    <tr>
      <td><b>A</b></td>
      <td><b>Appendix : technical tools</b></td>
      <td><b>38</b></td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>On W<sub>Γ<sub>D</sub></sub><sup>-1,p'</sup>(Ω) characterization</td>
      <td>38</td>
    </tr>
    <tr>
      <td>A.2</td>
      <td>PDE chain rule formula</td>
      <td>40</td>
    </tr>
    <tr>
      <td>A.3</td>
      <td>Week Aubin's result</td>
      <td>41</td>
    </tr>
  </tbody>
</table>

# 1 Introduction and preliminaries

## 1.1 Introduction

Dynamic systems generally aim to minimize their internal energy over time. By examining the dynamics at a discrete level, we can develop strategies that systematically guide the system towards the state of minimal energy using a step-by-step approach called steepest descent. One efficient approach involves managing the process of proximal energy minimization. This entails a dynamic procedure that continuously updates the system's state, aiming to decrease its internal energy while simultaneously minimizing the effort needed to transition from the current state to the subsequent one.

This approach has proven to be a robust tool for both describing continuous dynamics and establishing new dynamic systems that utilize their internal energy and the effort required to transition between states. Additionally, this approach has proven to be a valuable tool for devising dynamic strategies that aim to achieve energy-minimized states in various applications, including control theory, machine learning, and physics. However, most existing steepest descent algorithms often assess transition effort based on topological properties like gradient flow within specific spaces (e.g., Hilbert, Banach, metric Wasserstein), which may limit their flexibility. Our approach goes through the so-called minimum flow problem and offers a way to overcome this limitation, enabling steepest descent to minimize internal energy in a broader class of nonlinear PDEs like

$$ (1.1) \qquad \partial_t \rho - \nabla \cdot (\phi) = f, \quad \phi = \phi(x, \nabla \eta), \eta = \eta(x, \rho) \quad \text{for } (t, x) \in Q := (0, T) \times \Omega, $$
in an open and bounded regular domain $\Omega$, subject to non-homogeneous boundary conditions of Dirichlet, Neumann, or mixed type. Here $T > 0$ is an horizon time and $f$ is source term. The applications $A \in \mathbb{R}^N \to \phi(x, A) \in \mathbb{R}^N$ and $r \in \mathbb{R} \to \eta(x, r) \in \mathbb{R}$ are given maps, possibly multivalued, and assumed to be maximal monotone. One sees that, when the mobility function $\phi$ is linear, specifically $\phi(x, \nabla\eta) = \nabla\eta$, the resulting diffusion process is characterized by linear diffusion. For nonlinear mobility functions like $\phi(x, \nabla\eta) = |\nabla\eta|^{p-2}\nabla\eta$, the diffusion process exhibits $p$-Laplacian characteristics.

In addition to its theoretical value in proving existence and providing numerical approximations, our approach offers a meaningful interpretation that connects the solution to the system's internal energy and the various parameters, such as the source term and boundary conditions. It also establishes a direct link between the Euler-Implicit scheme, commonly used in the Crandall-Liggett $L^1$-theory, and a steepest descent algorithm that aims to minimize a proximal internal energy. This connection is possible despite the fact that the primary operator driving the dynamic is not necessarily a sub-differential operator in any usual Banach space ($L^p$, Sobolev or its dual space). In certain cases, our approach enables to restore the sub-differential operator process through the use of its metric version in dual Sobolev space. For homogeneous problems (homogeneous $\phi$ and boundary conditions), the approach can be simply linked to classical Hilbert gradient flow theory (cf. [14, 8]) in $H^{-1}$ for linear $\phi$; i.e., $\phi(x, \nabla\eta) = \nabla\eta$. Furthermore, for general $\phi$ of $p$-Laplacian type, it can be connected to metric gradient flow and $p'$-curve of maximal slope theory in the metric space $W^{-1,p'}$

We show, as well, that the approach provides a valuable tool for understanding and solving systems with dynamics that can be decomposed into two processes: a forcing process connected for instance to a given transport vector field, and minimizing internal energy process. More precisely; we show how one can use the approach to handle nonlinear diffusion-transport equation of the form

$$ (1.2) \qquad \partial_t \rho - \nabla \cdot (\phi + \rho V) = f, \quad \phi = \phi(x, \nabla\eta), \quad \eta = \eta(x, \rho) \quad \text{in } Q, $$

where $V = V(t, x)$ is a given transport fields. Roughly speaking, the dynamics of the system here can be decomposed into two processes: the transport process describes how the system moves under the action of the forcing vector field $V$, while the minimization energy process describes how the system tends to reach a stable equilibrium. Degenerate nonlinear parabolic problems like (1.2) received a lot of attention in the past couple of decades. For a non-exhaustive list of classical works on this subject we refer to [2, 15, 53, 49] and the references therein. The approach we present provides a way to understand how these two processes interact to produce the system's dynamics. In particular, it shows that the solution of the system is the state that minimizes the system's internal energy subject to the constraints imposed by the transport process and the boundary conditions. This interpretation has several implications. First, it provides a way to understand how the system's dynamics are affected by the different parameters. For example, the source term can be interpreted as a force that acts on the system, and the boundary conditions can be interpreted as constraints that limit the system's motion. Second, the interpretation can be used to design new algorithms for solving the system's dynamics from theoretical and numerical point of view.

Let us recall that the approach using the JKO algorithm in Wasserstein spaces has already been successfully applied to partially address this type of question (cf. [53, 4, 52, 41]). However, this theory's limitations to positive data, mass conservation, specific internal energies, and its technical complexity did not allow to broaden the scope of applications. The approach we present here draws inspiration from the JKO algorithm and the minimal flow problem (cf. [9]). It allows to cover a wide range of applications that go beyond $H^{-1}$-gradient flow framework while keeping the Wasserstein spirit for a large field of applications.
## 1.2 Preliminaries

More formally, let $X$ be a state space and $\mathcal{E}(\rho)$ be the internal energy of a given system at state $\rho \in X$. At each time $t > 0$, the dynamic system tends simultaneously to decrease its internal energy and to minimize the work required to move from state $\rho(t)$ to state $\rho(t+h)$. This can be achieved by solving the following optimization problem (proximal energy) at each time step:

$$ \min_{\rho \in X} \left\{ \mathcal{E}(\rho) + W_h(\rho(t), \rho) \right\}, $$

where $W_h(\rho(t), \rho)$ is a measure of the work required to move from state $\rho(t)$ to state $\rho$, throughout small time step $h > 0$. The solution to this optimization problem $\rho(t+h)$ is the next state of the system. To the point, one can consider the $\tau$-time steps defined by $t_0 = 0 < t_1 < \dots < t_n < T$, and the sequence of piecewise constant curves

$$ \rho_{\tau} = \sum_{i=1}^{n} \rho_{i} \chi_{[t_{i-1}, t_{i}]}, $$

where

$$ (1.3) \qquad \rho_i = \underset{\rho \in X}{\operatorname{argmin}} \left\{ \mathcal{E}(\rho) + W_{\tau}(\rho_{i-1}, \rho) \right\}, \quad \text{for } i = 1, \dots, n, \text{ with } \rho_0 = \rho(0). $$

Then, we expect the limit as $\tau \to 0$ of $\rho_{\tau}$ converges to the solution of the evolution problem (1.1). Working within the framework of Hilbert space $(X, \|\cdot\|)$, the general theory suggest to connect $W(\rho(t), \rho)$ to $\|\rho(t) - \rho\|^2/2\tau$ to obtain a representation of the solution through the gradient flow equation:

$$ (1.4) \qquad \begin{cases} \rho_t(t) + \partial \mathcal{E}(\rho(t)) \ni 0 & \text{in } (0, T) \\ x(0) & \text{is given,} \end{cases} $$

where $\partial \mathcal{E}(\rho)$ denote the sub-differential (possibly multi-valued map) of the functional $\rho \in X \to \mathcal{E}(\rho) \in (-\infty, \infty]$. The expression (1.3) can be related to both to the resolvent associated with the operator $\partial \mathcal{E}$ and Euler-Implicit discretisation of (1.4). The evolution problem in turn is a gradient flow in a Hilbert space which describes the evolution of a continuous curve $t \to \rho(t) \in X$ that follows the direction of steepest descent of the functional $\mathcal{E}$ (cf. [14]). One can consider a more general situation where $(X, d)$ is metric space and $W_h$ is build on the distance $d$. In this case, the expression (1.3) reveals the minimizing movement scheme à la De Giorgi (see [4]). The approach resulted in a broad theory of gradient flows in a metric space where the derivatives $\rho_t$ and $\partial_{\rho} \mathcal{E}(\rho)$ need to be interpreted in an appropriate way in $(X, d)$ (cf. the book [4]).

For instance, the porous medium-like equation

$$ (1.5) \qquad \partial_t \rho - \Delta \rho^m = 0 $$

in a bounded domain arises in various applications, especially in biology, to represent the evolution of species as they strive to minimize their internal energy $\mathcal{E}(\rho)$, which depends inherently on the density $\rho$. Since the works [51] and [39], it has been established that the PDE (1.5) subject to homogeneous Neuman boundary condition can be derived by employing $\mathcal{E}(\rho) = \frac{1}{m-1} \int \rho^m dx$ in the context of
W₂-Wasserstein distance by setting $W_τ(ρ, \tilde{ρ}) = \frac{1}{2τ}W_2(ρ, \tilde{ρ})^2$ (see Section 3.1 for the definition of W₂). It is well known that it can be also derived by utilizing the steepest descent algorithm for the internal energy $E(ρ) = \frac{1}{m+1} \int ρ^{m+1} dx$ in the context of $H^{-1}$-norm. This is achieved by setting $W_τ(ρ, \tilde{ρ}) = \frac{1}{2τ}||ρ - \tilde{ρ}||_{H^{-1}}^2$ (cf. [8]). Yet, one must carefully define the $H^{-1}$-norm to adequately address boundary conditions (see Section 3.2). These concepts can be expanded to accommodate some monotone nonlinearity $η(r)$ instead of $r^m$, both for W₂-Wasserstein distance (cf. [41] and the refs therein) and also $H^{-1}$-norm (cf. [21, 22, 40]) frameworks. Owing to the significant advancements in optimal mass transportation theory, the work of [1] permits handling doubly nonlinear PDE

$$ \partial_t ρ - \nabla \cdot φ(\nabla η(ρ)) = 0 $$

in the context of steepest descent algorithm and gradient flow with Wassertsein distance $W_c$, where $c$ is closely linked to the application $φ$ assumed in these works to be regular and satisfying the classical assumption of ellipticity. In particular, this enables providing an intriguing interpretation via steepest descent algorithm in Wasserstein space for porous medium-like equations governed by p-Laplacian nonlinear diffusion operator instead of the Laplacian. However, as far as we are aware, the extension of the concept to more general operators as in (1.1) remains unclear. Typical example may be given by space-dependent $φ$ and/or more general boundary conditions like Dirichlet, Neumann, or mixed homogeneous/non-homogeneous boundary conditions. The crux of the issue lies in the absence of a gradient flow framework to handle these types of extensions. Neither the $H^{-1}$ approach nor other possible extension to Sobolev dual approach appear to be capable of addressing these issues, even for the p-Laplacian operator. Our objective is to develop an alternative steepest descent-type approach for an internal energy that generalizes the $H^{-1}$-approach, in the spirit of Wasserstein distance. The interest of this approach stems from its physical interpretation of the dynamics studied and its interrelation with the system's internal energy and boundary conditions, even if the evolution PDE may not be related directly to a gradient flow dynamic.

Coming back to the class of PDE (1.1), remember that $ρ ≥ 0$ represents in many situations an unknown density of a fluid or a population, but $ρ$ could also be a sign-changing solution. In this case, the equation models the evolution of two types of materials/species whose densities are represented by the positive and negative parts of $ρ$, respectively. The equation highlights a second-order term to include linear or nonlinear diffusion phenomena through the term $∇ · φ(x, ∇η)$. However, as we can see, the diffusion does not act directly on $ρ$, but on the parameter $η$, which depends on $ρ$. In fluid mechanics, $η$ can be related to pressure, but in certain situations, it can simply appear as a corrective potential to adjust or reorient the transport phenomena of the model like in crowed motion (cf. [44], [45], [46], [27]). The relationship between $ρ$ and $η$, as well as the relationship between $φ$ and $∇η$, depends strongly on the physical phenomenon being modeled. This relationship can be linear, nonlinear, or multivalued, but it must always be monotone to ensure convexity of the energy and well-posedness. Typical example for the applications (possibly multi-valued) $(x, A) → φ(x, A) ∈ R^N$ and $(x, r) → η(x, r) ∈ R^N$ includes the following problems:

* Parabolic equation governed by a Leray-Lions type operator: $φ(x, A)$ is a Caratheodory applications satisfying classical ellipticity assumptions.

* Parabolic p-Laplacian equation: $φ(x, A) = k(x)|A|^{p-2}A$, where $0 < k$ is a given continuous function.
* Porous-Medium equation : $\eta(x, \rho) = \rho^m$

* Stefan problem : $\eta(x, r) = (r - a_1)^+ - (a_2 + r)^-$, where $a_1$ and $a_2$ are two given non-negative functions.

* Hele-Shaw problem $\eta(x, r) = a(x) \text{Sign}^{-1}(r)$, where Sign denotes the usual sign graph, and a is a given positive function.

Our aim in this paper is to explore this class of PDE through a “steepest descent” of an intrinsic internal energy functional we associate with the system against a minimum flow work. To the best of our knowledge, the proposed approach is a new and innovative perspective that is comparable to, but not necessarily identical to, the optimal mass transport approach (as we will see). Additionally, one of the key features of the proposed approach is its ability to handle changing sign solutions and a variety of boundary conditions, including mixed non-homogeneous and Robin boundary conditions. It is important to note that Wasserstein distances capture the movement of non-negative masses while maintaining mass conservation, making them suitable for describing the movement of particles within a given domain. However, the translation of Dirichlet boundary conditions into Wasserstein gradient flows remains a challenge (cf. [30]).

## 2 Plan of the paper

This paper focuses on an approach for solving evolution problems of the type (1.1) and (1.2) using the steepest descent method through minimum flow process. In Section 3, we provide a brief overview of this algorithm in the Wasserstein and $H^{-1}$-frameworks, followed by a formal description of our proposed approach. In Section 4, we introduce assumptions and review fundamental concepts related to Sobolev functions, their traces, and normal traces of vector-valued Lebesgue functions with Lebesgue-integrable divergence. Section 5 delves into a comprehensive analysis of the proximal energy minimization problem, which serves as the foundation for our steepest descent algorithm. Section 6 demonstrates how to employ the proximal energy minimization problem and steepest descent algorithm to investigate evolution problems of the type (1.1) and (1.2). Section 7 offers general remarks, perspectives, and open questions related to the proposed approach. At last, Appendix section presents some technical lemmas which we use in the paper.

## 3 Overview : Wasserstein distance vs $H^{-1}$-norm

The aim of this section is to remind the reader formally the main idea of steepest descent algorithms for nonlinear PDE in the context of $W_2$-Wasserstein distance and $H^{-1}$-Hilbert space.

### 3.1 Proximity through $W_2$-Wasserstein distance

Since the work of [39], the approach of steepest descent algorithms is well known and well used in the study of linear and nonlinear diffusion-convection phenomena. Here we focus the description on PDE of the type (1.5).

Remember, that $W_2$-Wasserstein distance between two Radon measures $\mu_1$ and $\mu_2$ with equal masses is defined as

$$W_2(\mu, \nu) = \left( \inf_{\gamma \in \mathcal{P}(\Omega \times \Omega)} \left\{ \int \|x - y\|^2 d\gamma(x, y) : \pi_x\#\gamma = \mu, \pi_y\#\gamma = \nu \right\} \right)^{1/2}$$
The $W_2$-Wasserstein space is the metric space of probability measures endowed with the metric $W_2$. The $W_2$-Wasserstein distance between probability measures $\mu \in \mathcal{P}(\Omega)$ and $\nu \in \mathcal{P}(\Omega)$ is given also by

$$ (3.1) \qquad W_2(\mu, \nu)^2 = \inf \left\{ \int_0^1 \int |V(t,x)|^2 d\rho(t,x) : (\rho, V) \in \mathcal{A}(\mu, \nu) \right\}, $$

where

$$ \begin{aligned} \mathcal{A}(\mu, \nu) &:= \{ (\rho, V) : \rho \in AC([0, 1]; \mathcal{P}_2(\Omega)), V \in L_{\rho}^2([0, 1] \times \Omega), \rho(0) = \mu, \rho(1) = \nu, \\ & \text{and } \partial_t \rho + \nabla \cdot (V \rho) = 0 \text{ in } [0, 1] \times \bar{\Omega} \}. \end{aligned} $$

Let us consider an internal energy functional $\mathcal{E}: \mathcal{P}_2(\Omega) \to [0, \infty]$ given by

$$ (3.2) \qquad \mathcal{E}(\rho) = \begin{cases} \int_{\Omega} f(\rho) \, dx & \text{if } f(\rho) \in L^1(\Omega) \\ +\infty & \text{otherwise} \end{cases} $$

where $f$ is a given convex function assumed here to be regular enough. The steepest descent algorithm associated with $\mathcal{E}$ through Wasserstein distance $W_2$, usually known as the JKO scheme since the seminal work [39], is given by the proximal algorithm

$$ \rho_{i+1} = \operatorname{argmin}_{\rho} \left\{ \frac{1}{2\tau} W_2(\rho, \rho)^2 + \mathcal{E}(\rho) \right\} \quad \text{for , } i = 0, 1, \dots, n $$

Then, taking an horizon time $T > 0$, $\tau$-time discretization $t_0 = 0 < t_1 < \dots < t_n < T$ and letting $\tau \to 0$ in the approximation

$$ \rho_{\tau} = \sum_{i} \rho_{i} \chi_{[t_{i}, t_{i+1}]} $$

one covers the solution of the evolution problem

$$ \left\{ \begin{array}{ll} \partial_t \rho - \nabla \cdot (\rho \nabla f'(\rho)) = 0 & \text{in } Q := (0, T) \times \Omega \\ \\ \rho \nabla f'(\rho) \cdot \nu = 0 & \text{on } \Sigma := (0, T) \times \partial \Omega \\ \\ \rho(0) = \rho_0. & \end{array} \right. $$

For classical results of this approach and related issues we refer the reader to the books [4, 52] for a thoroughgoing results as well as a corresponding literature

## 3.2 Proximity through $H^{-1}$-norm

To keep the same context as the Wasserstein case, which conserves the total mass, we present here the core concept of an appropriate framework for the $H^{-1}$ approach (for more in-depth discussions, we refer to [40, 22]). The conservation of mass involves considering the quotient vector space

$$ H := H^1(\Omega)/\mathbb{R}. $$
That is the quotient space associated with the binary relation defined in $H^1(\Omega)$, by

$$u_1 \mathcal{R} u_2 \iff u_1 - u_2 = c \text{ a.e in } \Omega, \text{ with } c \in \mathbb{R}.$$

To avoid any confusion, we will identify the vectors $\bar{u} \in H$, with their corresponding elements $u \in \bar{u}$. The space $H$ is Hilbert space when equipped with the inner product

$$a(u, v) = \int_{\Omega} \nabla u \cdot \nabla v \, dx, \quad \text{for any } u, v \in H,$$

and its associate norm

$$\|u\|_{H^{-1}} = \|\nabla u\|_{L^2(\Omega)}, \quad \text{for any } u \in H.$$

Then, let us consider $H^{-1}$, the dual space (topological) of $H$. Let us denote by $J$ the duality map of $H$ and use the notation $\langle \mu, z \rangle$ for the value of $\mu \in H^{-1}$ at $u \in H$. We denote by $\| \cdot \|_{H^{-1}}$ the dual norm in $H^{-1}$ induced by $\| \cdot \|_{H}$ in $H$. The duality map $J$ may be defined by $J: H^{-1} \to H$ and $J(\mu) = u$, where $u$ is a solution of the Neuman boundary problem governed by the Laplace operator

$$\begin{cases} -\Delta u = \mu & \text{in } \Omega \\ \partial_{\eta} u = 0 & \text{on } \partial\Omega. \end{cases}$$

Moreover, the duality bracket $\langle \cdot, \cdot \rangle$ is given by

$$\langle \mu, u \rangle = \int_{\Omega} \nabla J(\mu) \cdot \nabla u \, dx, \quad \text{for any } \mu \in H^{-1} \text{ and } u \in H.$$

In particular this implies that

$$\|\mu\|_{H^{-1}} = \|\nabla J(\mu)\|_{L^2(\Omega)}.$$

The steepest descent algorithm associated with a given internal energy $\mathcal{E}$ through $H^{-1}$-norms appears through the usual Euler-Implicit time discretization associated with the evolution problem

$$(3.3) \qquad \frac{du}{dt} + \partial \mathcal{E}(u) \ni 0 \quad \text{in } (0,T),$$

where $\partial \mathcal{E}$ is given by the usual definition of sub-differential in Hilbert space which is $H^{-1}$ here. Indeed, taking a $\tau$-time discretization $t_0 = 0 < t_1 < \dots < t_n < T$ the Euler-Implicit time discretization aims to solve the stationary problem

$$u_{i+1} + \tau \partial \mathcal{E}(u_{i+1}) \nu u_i, \quad \text{for each } i = 0, \dots, n,$$

where $u_0$ is the initial data associated with the evolution problem (3.3). In particular, $u_{i+1}$ may be given by

$$u_{i+1} = \operatorname*{argmin}_{\rho} \left\{ \frac{1}{2\tau} \|u_i - u\|_{H^{-1}}^2 + \mathcal{E}(u) \right\} \quad \text{for } i = 0, 1, \dots, n$$

Then, letting $\tau \to 0$ in the approximation

$$u_{\tau} = \sum_{i} u_{i} \chi_{[t_{i}, t_{i+1}]}$$
one covers the solution of the evolution problem (3.3). For the internal energy (3.2), the associate PDE read

$$ egin{cases} 
\partial_t u - \Delta f'(u) = 0 & \text{in } Q \\
\nabla f'(u) \cdot \nu = 0 & \text{on } \Sigma \\
u(0) = u_0
\end{cases} $$

**Remark 1** In the non-conservativ case, the usual $H^{-1}$-approach consists to take simply the dual space of $H_0^1(\Omega)$ equipped with its usual norm. In this case, the associate nonlinear PDE is subject to homogeneous Dirichlet boundary condition ([8], one can see also remark 7.1 in Section 5).

### 3.3 Minimum flow concept and related steepest descent

The divergence differential operator can be used to model the movement of particles between regions with different concentrations. The equation

$$ (3.4) \qquad -\nabla \cdot \phi = \mu - \nu \quad \text{in } \Omega $$

characterizes the traffic flow $\phi$ between two distributions represented by a source $\mu$ and a target $\nu$ within an open bounded domain $\Omega$. The minimum flow problem, originally formulated by Beckman (cf. [9]), involves finding an optimal traffic flow by minimizing a specified cost function on $\phi$. Inspired by the concept presented in [9], we propose a generalization of the minimum flow problem that can accommodate the study of relevant practical scenarios for (1.2). For each given traffic flow, there is a corresponding traffic cost involving two outlays: $F(x, \phi(x))$ in $\Omega$ and $G(x, -\phi(x) \cdot \nu(x))$ on $\Gamma_D$ a specific portion of the boundary. Here, $F: \Omega \times \mathbb{R}^N \to \mathbb{R}$ and $G: \partial\Omega \times \mathbb{R} \to \mathbb{R}$ are given functions. In this context, the optimal traffic cost is defined as the minimum total cost of the traffic, which corresponds to the minimum flow of the traffic cost

$$ (3.5) \qquad \int_{\Omega} F(\cdot, \phi) \, dx + \int_{\Gamma_D} G(\cdot, -\phi \cdot \nu) \, d\sigma. $$

The goal is to find the traffic scheme that results in the lowest overall cost, taking into account the cost functions $F$ and $G$ applied to the flow quantities $\phi$ and $\phi \cdot \nu$ at each location in the domain $\Omega$ and on the boundary $\Gamma_D$, respectively.

In other words, we define an intrinsic quantity, denoted for the moment by $I(\mu, \nu)$, for each distribution of mass $\mu$ and $\nu$. This quantity is defined as the minimum of the value (3.5) among all vector valued flux $\phi$ satisfying the balance equation (3.4). It characterizes the movement of mass from $\mu$ to $\nu$. The concept of $I(\mu, \nu)$ is reminiscent of optimal mass transportation, where the goal is to find the most efficient way to transport mass from one distribution to another while minimizing the total traffic cost. Similarly, in our context, $I(\mu, \nu)$ can be interpreted as a transfer fee that captures the essential characteristics of the mass transfer from $\mu$ to $\nu$, providing insight into the optimal traffic process. Observe that equation (3.4) can alternatively be viewed as a transport equation governing the transfer of mass $\mu$ into $\nu$. This can be seen by introducing the following variables:

$$ V = -\phi/\rho, $$

where

$$ \rho(t) = t\nu + (1-t)\mu, \quad \text{for any } t \in [0, 1]. $$
By adopting this formulation, (3.4) shares formally the same objective as a transport equation,

$$ \begin{cases} \partial_t \rho + \nabla \cdot (\rho V) = 0 & \text{in } (0,1) \times \Omega \\ \rho(0) = \mu, \quad \rho(1) = \nu, \end{cases} $$

which is to describe the movement of a quantity, represented by $\rho$, across a given domain through the vector field $V$. Yet, this is just formal since the velocity $V$ may be here infinite in general. See in this case that the total cost (3.5) may be written as

$$ (3.6) \qquad \int_0^1 \int_{\Omega} F(x, \rho(t, x) V(t, x)) dxdt + \int_0^1 \int_{\Gamma_D} G(x, -\rho(t, x)V(t, x) \cdot \nu) d\sigma dt. $$

While the two formulas (3.6) and (3.1) are not the same, it is noteworthy that the expressions (3.6) bears resemblance to the Wasserstein distance expression (3.1). This similarity suggests that the approach we are developing here shares some of the characteristics of optimal transport, which is already more or less established for the $H^{-1}$-approach (cf. [53]). This observation motivates the use of the quantity $I$ to measure transition work in a steepest descent algorithm for minimizing internal energy. Another compelling reason to incorporate $I$ into the proximal energy minimization algorithm is that the transfer phenomenon we are addressing arises directly and instantaneously from diffusion phenomena with different boundary conditions. Indeed, formal duality arguments can demonstrate that the Euler-Lagrange equation associated with $I(\mu, \nu)$ is closely related to elliptic problem :

$$ (3.7) \qquad \left\{ \begin{array}{ll} -\nabla \cdot \dot{\phi} = \mu - \nu & \text{in } \Omega \\ \phi \in \partial_{\xi} F^*(x, \nabla u) & \text{on } \Gamma_D \\ \mu/\Gamma - \dot{\phi} \cdot \nu \in \partial G^*(x, u) & \text{on } \Gamma_D \\ \phi \cdot \nu = 0 & \text{on } \partial\Omega \setminus \Gamma_D, \end{array} \right. $$

where $F^*$ and $G^*$ denote the usual Legendre transforms of $F$ and $G$ respectively (see Section 4.2 for the definitions).

To simplify the presentation, we focus in this paper on the case where $G(x, -\phi \cdot \nu)$ is proportional to $-\phi \cdot \nu$. More precisely, we assume that

$$ G(x, r) = g(x) r, \quad \text{for any } r \in \mathbb{R} \text{ and a.e. } x \in \Gamma_D, $$

for a given $g$ assumed to be the trace of a Sobolev function. In some sense, $g$ represents some given charge when the mass $\rho$ moves onto the boundary $\Gamma_D$. As we will see, the boundary condition on $\Gamma_D$ simply turns into non-homogeneous Dirichlet boundary condition in this case. The right assumptions on $\Omega$, $\Gamma_D$, $F$, $g$ rigorous definition of $I$, and its connection with (3.7) typed PDE will be given in details in the next section.

This being said, we can consider a dynamic system in $\Omega$ with an internal energy $\int_{\Omega} \beta(., \rho) dx$, associated with densities $\rho$. The proximal optimization algorithm aims to construct a sequence $\rho^i$ initialized with a given initial density $\rho_0$. For $i = 0, 1, ...n$, the algorithm establishes $\rho^i$, by minimizing the problem:

$$ (3.8) \qquad \operatorname{argmin}_{\rho} \left\{ \int_{\Omega} \beta(., \rho) \, dx + I^{\tau}(\rho^{i-1}, \rho) \right\}, $$
where $\tau > 0$ is a given time step and, $I^\tau$ corresponds to the $\tau$-scaled transfer fees of $\rho^{i-1}$ into $\rho$ given through $I(\rho^{i-1}, \rho)$. It is constructed up to the minimum of the amount

$$ \int_{\Omega} F(., \phi) dx - \int_{\Gamma_D} g \phi \cdot \nu d\sigma, $$

among the vector valued fluxes given by the balance PDE

$$ (3.9) \qquad -\nabla \cdot \phi = \rho^{i-1} - \rho \quad \text{in } \Omega. $$

Here, on the Dirichlet boundary $\Gamma_D$, we'll keep the normal trace of $\phi$ free to vary. And, we fix it on the remaining part $\partial\Omega \setminus \Gamma_D$. This leads to Neumann boundary condition (not necessarily homogeneous) for the PDE (3.9) (see more details in Section 5).

**Remark 2** We observe that the term $I(\rho^{i-1}, \rho)$ depends only on the quantity $\rho^{i-1} - \rho$. Without any ambiguity, we will consider $I(\rho^{i-1}, \rho)$ definitely as $I(\rho^{i-1} - \rho)$. Furthermore, we will utilize the balance equation (3.4) under general considerations, including non-homogeneous boundary conditions and dual Sobolev source terms. That is, we will consider equation (3.4) in its most general form to accommodate a wider range of practical applications, including the case of non-homogeneous boundary conditions and the case of a forcing term of the type $\nabla \cdot \chi$ (see the exact definition of $I$ in (5.2) and Remark 5).

# 4 Notations and assumptions

Assume that $\Omega$ is bounded regular domain with Lipschitz boundary $\partial\Omega$ satisfying

$$ \partial\Omega = \Gamma_D \cup \Gamma_N, $$

where $\Gamma_D$ and $\Gamma_N$ are assumed to be disjoint regular connected components. We assume moreover that

$$ \mathcal{H}^{N-1}(\Gamma_D) > 0. $$

## 4.1 Reminder

Let $1 < p < \infty$ and $1 < p' < \infty$ be such that

$$ \frac{1}{p} + \frac{1}{p'} = 1. $$

We focus in this paper on the case where $1 < p < \infty$.

* $L^p(\Omega)$ and $L^1(\Omega)$ denotes the usual Lebesgue spaces endowed with their natural norms. Unless otherwise stated, and without abuse of notation, we will use $\int f$ to denote the Lebesgue integral over $\Omega$, i.e.

$$ \int_{\Omega} h := \int_{\Omega} h(x) dx, \quad \text{for any } h \in L^1(\Omega). $$

* $W^{1,p}(\Omega)$ denotes the usual Sobolev space endowed with its natural norm

$$ \|z\|_{W^{1,p}(\Omega)} = \left( \int_{\Omega} |z|^p + \int_{\Omega} |\nabla z|^p \right)^{1/p}. $$
* $W_{\Gamma_D}^{1,p}(\Omega)$ is the closure, in $W^{1,p}(\Omega)$, of $C^1(\bar{\Omega})$ which are null on $\Gamma_D$. In particular, for any $z \in W_{\Gamma_D}^{1,p}(\Omega)$, $\gamma(z) = 0$ on $\Gamma_D$, where $\gamma$ denotes the usual trace application defined from $W^{1,p}(\Omega) \to L^p(\partial\Omega)$.

* We denote by $W_{\Gamma_D}^{-1,p'}(\Omega)$, the dual space of $W_{\Gamma_D}^{1,p}(\Omega)$. Thanks to Lemma A.7, we know that, for any $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$, there exists a couple $(f_0, \bar{f}) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$, such that

$$ \langle f, \xi \rangle_{W_{\Gamma_D}^{-1,p'}(\Omega), W_{\Gamma_D}^{1,p}(\Omega)} = \int_{\Omega} f_0 \xi - \int_{\Omega} \bar{f} \cdot \nabla \xi, \quad \text{for any } \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Without abusing, we will use the couple $(f_0, \bar{f})$ to identify $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$ and denote the duality bracket by $\langle \cdot, \cdot \rangle_{\Omega}$; i.e.

$$ \langle f, \xi \rangle_{\Omega} := \langle f, \xi \rangle_{W_{\Gamma_D}^{-1,p'}(\Omega), W_{\Gamma_D}^{1,p}(\Omega)}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

* For any $\tilde{\Gamma} \subset \partial\Omega$, we denote by $W^{1-1/p,p}(\tilde{\Gamma}) := \gamma_{/\tilde{\Gamma}}(W^{1,p}(\Omega))$, where $\gamma_{/\tilde{\Gamma}}$ is the trace application restricted to $\tilde{\Gamma}$. We need to define moreover

$$ W_{00}^{1-1/p,p}(\tilde{\Gamma}) := \{ \kappa \in L^{p'}(\tilde{\Gamma}) : \exists \tilde{\kappa} \in W^{1,p}(\Omega), \gamma(\tilde{\kappa}) = \kappa \text{ on } \tilde{\Gamma} \text{ and } \tilde{\kappa} = 0 \text{ on } \partial\Omega \setminus \tilde{\Gamma} \}. $$

Thanks to the assumptions on $\partial\Omega$, taking $\tilde{\Gamma} = \Gamma_D$ or $\tilde{\Gamma}_N$, the space $W_{00}^{1-1/p,p}(\tilde{\Gamma})$ coincides with the set of functions belonging to $W^{1-1/p,p}(\tilde{\Gamma})$ and vanishing on the boundary $\partial\Omega \setminus \tilde{\Gamma}$. Remember that this not automatically true for any $\tilde{\Gamma} \subset \partial\Omega$.

* For any $\tilde{\Gamma} \subset \partial\Omega$, we denote by $W^{-1/p',p'}(\tilde{\Gamma}) := (W_{00}^{1-1/p,p}(\tilde{\Gamma}))'$, the topological dual space of $W_{00}^{1-1/p,p}(\tilde{\Gamma})$. We denote the duality bracket simply by the formal expression $\langle \pi, \kappa \rangle_{\tilde{\Gamma}}$; i.e.

$$ \langle \pi, \kappa \rangle_{\tilde{\Gamma}} := \langle \pi, \kappa \rangle_{W^{-1/p',p'}(\tilde{\Gamma}), W_{00}^{1-1/p,p}(\tilde{\Gamma})}, \quad \forall \kappa \in W_{00}^{1-1/p,p}(\tilde{\Gamma}) \text{ and } \pi \in W^{-1/p',p'}(\tilde{\Gamma}). $$

* We can define simultaneously the space

$$ H^{p'}(div, \Omega) := \{ v \in L^{p'}(\Omega)^N : \nabla \cdot v \in L^{p'}(\Omega) \}, $$

where $\nabla \cdot v$ is taken in $D'(\Omega)$. The space $H^{p'}(div, \Omega)$ endowed with the norm

$$ \|v\|_{H^{p'}(div, \Omega)} := \|v\|_{L^{p'}(\Omega)} + \|\nabla \cdot v\|_{L^{p'}(\Omega)} $$

is a Banach space.

* Let $v \in H^{p'}(div, \Omega)$. For any $\tilde{\Gamma} \subset \partial\Omega$, the normal trace $v \cdot \nu$ is well defined on $\tilde{\Gamma}$ by duality. More precisely, $v \cdot \nu \in W^{-1/p',p'}(\tilde{\Gamma})$, and

$$ (4.1) \qquad \langle v \cdot \nu, \kappa \rangle_{\tilde{\Gamma}} = \int_{\Omega} v \cdot \nabla \tilde{\kappa} \, dx + \int_{\Omega} \tilde{\kappa} \nabla \cdot v \, dx, $$

for any $\kappa \in W^{1-1/p,p}(\tilde{\Gamma})$ and $\tilde{\kappa} \in W^{1,p}(\Omega)$ such that $\tilde{\kappa} = \kappa$ on $\tilde{\Gamma}$, and $\tilde{\kappa} = 0$ on $\partial\Omega \setminus \tilde{\Gamma}$.
We conclude this reminder section with the following result, which is relatively well-known in the field of elliptic PDEs.

**Lemma 1** For any $\mu \in L^{p'}(\Omega)$, $\chi \in L^{p'}(\Omega)^N$ and $\eta \in W^{-1/p', p'}(\Gamma_N)$, there exists $\phi \in H^{p'}(\text{div}, \Omega)$ solution of the following PDE :

$$ (4.2) \qquad \begin{cases} -\nabla \cdot \phi = \mu + \nabla \cdot \chi & \text{in } \Omega \\ (\phi + \chi) \cdot \nu = \pi & \text{on } \Gamma_N, \end{cases} $$

in the sense that

$$ \int_{\Omega} (\phi + \chi) \cdot \nabla \xi \, dx = \int_{\Omega} \mu \xi \, dx + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \text{for any } \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

**Proof** : It is enough to work with the optimization problem

$$ (4.3) \qquad \min_{z \in W_{\Gamma_D}^{1,p}(\Omega)} \left\{ \frac{1}{p} \int_{\Omega} |\nabla z|^p \, dx - \int_{\Omega} \chi \cdot \nabla z \, dx + \int_{\Omega} \mu z \, dx + \langle \pi, z \rangle_{\Gamma_N} \right\}. $$

Using standard argument of calculus of variation, one can prove that (4.3) has a solution that we denote by $u_p$. Moreover, we can prove that $u_p$ satisfies

$$ \int_{\Omega} |\nabla u_p|^{p-2} \nabla u_p \cdot \nabla \xi \, dx = \int_{\Omega} \mu \xi \, dx - \int_{\Omega} \chi \cdot \nabla \xi \, dx + \int_{\Gamma_N} [\pi, \xi], \quad \text{for any } \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Taking $\phi = |\nabla u_p|^{p-2} \nabla u_p$, the proof of the lemma is complete. $\square$

**Remark 3** In terms of $W_{\Gamma_D}^{-1,p'}(\Omega)$, the equation (4.2) is equivalent to write $-\nabla \cdot \phi = f$ in $W_{\Gamma_D}^{-1,p'}(\Omega)$, where $f$ is given by

$$ \langle f, \xi \rangle_{\Omega} = \int_{\Omega} \mu \xi \, dx - \int_{\Omega} \chi \cdot \nabla \xi \, dx + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

## 4.2 Assumptions

We consider the applications $F: \Omega \times \mathbb{R}^N \to \mathbb{R}$ and $\beta: \Omega \times \mathbb{R} \to \mathbb{R}$ to be two Caratheodory applications; i.e. the applications $(x, \xi) \in \Omega \times \mathbb{R}^N \to F(x, \xi) \in \mathbb{R}$ and $(x, r) \in \Omega \times \mathbb{R} \to \beta(x, r) \in \mathbb{R}$ are measurable with respect to $x \in \Omega$, and continuous with respect to $\xi \in \mathbb{R}^N$ and $r \in \mathbb{R}$, respectively. We assume moreover that

(H1) there exists two constants $C_1, C_2 > 0$, such that

$$ C_1 \|\xi\|^{p'} \le F(x, \xi) \le C_2 \|\xi\|^{p'}, \quad \text{for any } \xi \in \mathbb{R}^N, $$

(H2) the application $\xi \in \mathbb{R}^N \to F(x, \xi) \in \mathbb{R}$ is lower semi-continuous (l.s.c. for short) and convex, for a.e. $x \in \Omega$,

(H3) for a.e. $x \in \Omega$ and for a.e. $x \in \Omega$, the application

$$ r \in \mathbb{R} \to \beta(x, r) \text{ is l.s.c., convex and } \beta(x, 0) = 0, $$
(H4) there exists $C_3$, $M > 0$ such that

$$C_3 (|r| - M)^{+p'} \le \beta(x, r), \quad \text{for any } r \in \mathbb{R} \text{ and a.e. } x \in \Omega.$$

Under the assumptions (H1)-(H4), we can define the Legendre transforms of $F$ and $\beta$ in the second variable as usual by

$$F^*(x, \xi^*) = \max_{\xi \in \mathbb{R}^N} \{\xi^* \cdot \xi - F(x, \xi)\}, \quad \text{for a.e. } x \in \Omega \text{ and } \xi^* \in \mathbb{R}^N$$

and

$$\beta^*(x, r^*) = \max_{r \in \mathbb{R}} \{r^* \cdot r - \beta(x, r)\}, \quad \text{for a.e. } x \in \Omega \text{ and } r^* \in \mathbb{R}.$$

Let us remind the reader here that $F^*$ and $\beta^*$ are again l.s.c. and convex with respect to the second variable. Moreover,

$$F^*(x, \xi^*) = \xi^* \cdot \xi - F(x, \xi) \Leftrightarrow \xi^* \in \partial_{\xi} F(x, \xi)$$

and

$$\beta^*(x, r^*) \in r^* \cdot r - \beta(x, r) \Leftrightarrow r^* \in \partial_r \beta(x, r).$$

These are also equivalent to $\xi \in \partial_{\xi^*} F^*(x, \xi^*)$ and $r \in \partial_{r^*} \beta^*(x, r^*)$, respectively. Moreover, thanks to the assumptions (H1), it is not difficult to see that

(H'1) there exists two constants $C'_1, C'_2 > 0$, such that

$$C'_1 \|\xi^*\|^p \le F^*(x, \xi^*) \le C'_2 \|\xi^*\|^p, \quad \text{for any } \xi^* \in \mathbb{R}^N \text{ and a.e. } x \in \Omega.$$

For the main facts and classical results about conjugate convex functions, one can see for instance the book [26].

**Remark 4** The results of the paper remains to be true if we replace the assumption (H3) by the following one

$$C_1 \|\xi\|^{p'} - c_1(x) \le F(x, \xi) \le C_2 \|\xi\|^{p'} + c_2(x), \quad \text{for any } \xi \in \mathbb{R}^N, \text{ and a.e. } x \in \Omega,$$

where $c_1$ and $c_2$ are two $L^1$ functions on $\Omega$. Also, they remains to be true if we replace the assumption (H4) by

$$C_3 (|r| - M)^{+p'} - c_3(x) \le \beta(r), \quad \text{for any } r \in \mathbb{R} \text{ and a.e. } x \in \Omega,$$

where $c_3$ is a $L^1$ function on $\Omega$.

# 5 Proximal energy minimization

## 5.1 Main results

Let $1 < p < \infty$

$$ (5.1) \qquad \pi \in W^{-\frac{1}{p'}, p'}(\Gamma_N) \quad \text{and} \quad g \in W_{00}^{1-\frac{1}{p'}, p}(\Gamma_D) $$
be given and fixed throughout this paper. We denote by $\tilde{g}$ the $W^{1,p}(\Omega)$ function such that $\tilde{g} = 0$ on $\Gamma_N$ and $\tilde{g} = g$ on $\Gamma_D$. To define the transition work in (3.8) rigorously, for any $\mu \in L^{p'}(\Omega)$ and $\chi \in (L^{p'}(\Omega))^N$, we consider

$$ \mathcal{A}_{\pi}^{\chi}(\mu) = \{ \phi \in L^{p'}(\Omega)^N, -\nabla \cdot \phi = \mu + \nabla \cdot \chi \text{ in } \Omega \text{ and } (\phi + \chi) \cdot \nu = \pi \text{ on } \Gamma_N \}. $$

That is $\phi \in \mathcal{A}_{\pi}^{\chi}(\mu, h)$ if and only if $\phi \in L^{p'}(\Omega)^N$ and

$$ \int_{\Omega} \phi \cdot \nabla \xi \, dx = \int_{\Omega} \mu \xi \, dx - \int_{\Omega} \chi \cdot \nabla \xi \, dx + \langle \pi, \xi \rangle_{\Gamma_N} \quad \text{for any } \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Thanks to Lemma 1, under the assumptions (5.1), we have

$$ \mathcal{A}_{\pi}^{\chi}(\mu) \neq \emptyset, \quad \text{for any } \mu \in L^{p'}(\Omega) \text{ and } \chi \in L^{p'}(\Omega)^N. $$

Now, for any $\mu \in L^{p'}(\Omega)$ and $\chi \in L^{p'}(\Omega)^N$, fixed, we define transition work as follows

$$ (5.2) \qquad \mathcal{I}_{\pi,g}^{\chi}(\mu) := \inf_{\phi} \left\{ \int_{\Omega} F(x, \phi(x)) \, dx - \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D} : \phi \in \mathcal{A}_{\pi}^{\chi}(\mu) \right\}. $$

**Remark 5** As we said in Section 2.2, the quantity $\mathcal{I}_{\pi,g}^{\chi}(\mu)$ may be considered as the transition work associated with the distribution of mass $\mu$. However, revisiting the formulation (3.4) we note that in the definitions of $\mathcal{A}$ and $\mathcal{I}$, the source term of (3.4) is generalized to a more comprehensive form, $\mu + \nabla \cdot \chi$. Indeed, the term $\nabla \cdot \chi$ can be utilized to represent a significant external influence exerted on the system by an external force, $\chi$ (along with boundary conditions). In some sense $\mathcal{I}_{\pi,g}^{\chi}(\mu)$ may be interpreted as a transition work balancing the distribution $\mu$ under external force $\chi$. The main future of this interpretation is to encompass a broad spectrum of applications, including PDEs of the type (1.1) with a fixed forcing term incorporated into the source term, f. Moreover, it can handle also PDEs of the type (1.2), in which the forcing term is represented by a transport phenomenon clearly identified in the first-order term, $\nabla \cdot (\rho V)$.

Then, to solve the steepest descent algorithm associated with transition work $\mathcal{I}_{\pi,g}^{\chi}$ and the internal energy $\int_{\Omega} \beta(., \rho)$, we consider the proximal energy minimization problem

$$ (5.3) \qquad \mathcal{N}_{\pi,g}^{\chi}(\mu) := \inf_{\rho} \left\{ \int_{\Omega} \beta(., \rho) + \mathcal{I}_{\pi,g}^{\chi}(\mu - \rho), \rho \in L^{p'}(\Omega) \right\}; $$

i.e.

$$ \mathcal{N}_{\pi,g}^{\chi}(\mu) := \inf_{\rho, \phi} \left\{ \int_{\Omega} \beta(., \rho) + \int_{\Omega} F(., \phi) \, dx - \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D} : \phi \in \mathcal{A}_{\pi}^{\chi}(\mu - \rho), \rho \in L^{p'}(\Omega) \right\}. $$

See that all the terms in $\mathcal{N}_{\pi,g}^{\chi}(\mu)$ are well defined. Indeed, for any $\phi \in \mathcal{A}_{\pi}^{\chi}(\mu - \rho)$, we have $\phi + \chi \in H^{p'}(\text{div}, \Omega)$, so that the normal trace of $\phi + \chi$ is well defined on $\Gamma$. Remember that $\langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D}$ needs to be understood in the sense of

$$ \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D} = \int_{\Omega} (\phi + \chi) \cdot \nabla \tilde{g} \, dx + \int_{\Omega} \tilde{g} \nabla \cdot (\phi + \chi) \, dx. $$
**Remark 6** The optimization problem $\mathcal{N}_{\pi,g}^{\chi}(\mu)$ seeks to identify the proximal state, $\rho$ for $\mu$, relative to the total energy, which simultaneously minimizes the internal energy $\int_{\Omega} \beta(x, \rho) dx$ and the transition work $\mathcal{I}_{\pi,g}^{\chi}$ required to reach a proximal state of $\mu$ driven by the forcing term $\chi$. Notably, no assumption is made regarding the signs of either $\mu$ or $\rho$. Nevertheless, practical considerations may necessitate imposing sign constraints on $\rho$. This must be carefully considered in subsequent analysis and interpretation in terms of PDE (see Remark 7).

Our first main result of this section is the following characterization in terms of PDE of the solution to (5.2).

**Theorem 1** For any $\mu \in L^{p'}(\Omega)$ and $\chi \in L^{p'}(\Omega)^N$, we have:

$$ (5.4) \qquad \mathcal{I}_{\pi,g}^{\chi}(\mu) = \underbrace{\max_{\eta \in W^{1,p}(\Omega), \eta/\Gamma_D = g} \left\{ \int_{\Omega} \mu \eta - \int_{\Omega} \chi \cdot \nabla \eta - \int_{\Omega} F^*(., \nabla \eta) + \langle \pi, \eta \rangle_{\Gamma_N} \right\}}_{=: \mathcal{J}_{\pi,g}^{\chi}(\mu)} $$

Moreover,

1. $\mathcal{I}_{\pi,g}^{\chi}(\mu)$ has a solution $\phi \in L^{p'}(\Omega)^N$.

2. $\phi$ and $\eta$ are solutions of $\mathcal{I}_{\pi,g}^{\chi}(\mu)$ and $\mathcal{J}_{\pi,g}^{\chi}(\mu)$, respectively, if and only if they satisfies the following PDE:

$$ (5.5) \qquad \left\{ \begin{array}{ll} -\nabla \cdot \phi = \mu + \nabla \cdot \chi & \text{in } \Omega, \\ \phi = \partial_{\xi} F^*(., \nabla \eta) & \text{on } \Gamma_N, \\ (\phi + \chi) \cdot \nu = \pi & \text{on } \Gamma_D, \\ \eta = g & \text{on } \Omega \end{array} \right. $$

in the sense that $\eta \in W^{1,p}(\Omega)$, $\eta = g$ on $\Gamma_D$, $\phi \in L^{p'}(\Omega)^N$, $\phi \in \partial F^*(., \nabla \eta)$ a.e. in $\Omega$, and we have

$$ \int_{\Omega} \phi \cdot \nabla \xi \, dx = \int_{\Omega} \mu \xi \, dx - \int_{\Omega} \chi \cdot \nabla \xi \, dx + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

As to (5.3), we have

**Theorem 2** For any $\mu \in L^{p'}(\Omega)$ and $\chi \in L^{p'}(\Omega)^N$, we have:

$$ \mathcal{N}_{\pi,g}^{\chi}(\mu) = \underbrace{\max_{\eta \in W^{1,p}(\Omega), \eta/\Gamma_D = g} \left\{ \int_{\Omega} \mu \eta - \int_{\Omega} \chi \cdot \nabla \eta - \int_{\Omega} F^*(., \nabla \eta) - \int_{\Omega} \beta^*(., \eta) + \langle \pi, \eta \rangle_{\Gamma_N} \right\}}_{=: \mathcal{D}_{\pi,g}^{\chi}(\mu)} $$

Moreover,
1. $\mathcal{N}_{\pi,g}^{\chi}(\mu)$ has a solution $(\rho, \phi) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$.

2. $(\rho, \phi)$ and $\eta$ are solutions of $\mathcal{N}_{\pi,g}^{\chi}(\mu)$ and $D_{\pi,g}^{\chi}(\mu)$, respectively, if and only if they satisfies the following PDE :

$$ (5.6) \qquad \left\{ \begin{array}{ll} \rho - \nabla \cdot \phi = \mu + \nabla \cdot \chi, & \text{in } \Omega, \\ \rho \in \partial \beta^*(x, \eta), \quad \phi = \partial_{\xi} F^*(., \nabla \eta) & \\ (\phi + \chi) \cdot \nu = \pi & \text{on } \Gamma_N, \\ \eta = g & \text{on } \Gamma_D \end{array} \right. $$

in the sense that $\rho \in L^{p'}(\Omega)$, $\eta \in W^{1,p}(\Omega)$, $\phi \in L^{p'}(\Omega)^N$, $\eta = g$ on $\Gamma_D$, $\phi \in \partial F^*(., \nabla \eta)$, $\rho \in \partial \beta^*(x, \eta)$, a.e. in $\Omega$, and

$$ \int_{\Omega} \rho \xi + \int_{\Omega} (\phi + \chi) \cdot \nabla \xi = \int_{\Omega} \mu \xi + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

**Remark 7** 1. We do not treat uniqueness concerns of the solution for the problem (5.6). However, one sees that primal issues in this direction can be deduced by using the optimization problems $\mathcal{N}_{\pi,g}^{\chi}$, $\mathcal{J}_{\pi,g}^{\chi}$ and strict convexity of $\beta$, $F$, $\beta^*$ and/or $F^*$.

2. In certain practical situations, the requirement that the density $\rho$ in the proximal energy optimization problem $\mathcal{N}_{\pi,g}^{\chi}$ remains non-negative arises. The feasibility of this constraint hinges on the specific conditions imposed on the data $\mu$, $\chi$, $g$ and $\pi$. For instance, if $g \equiv 0$, $\mu + \nabla \cdot \chi \ge 0$ and $\pi \ge 0$, it is possible to prove that the solution $\rho$ to the problem $\mathcal{N}_{\pi,0}(\mu)$ is conceivable non-negative, by utilizing the associated PDE (5.6). In the absence of these conditions, one must explicitly incorporate the constraint $\rho \ge 0$ into the optimization problem $\mathcal{N}$, thereby necessitating the introduction of Lagrangian multipliers.

The main future of the PDE (5.5) deserves the evolution dynamic of the type

$$ \left\{ \begin{array}{ll} \partial_t \rho - \nabla \cdot (\phi - V(x, \rho, \eta)) = \mu, & \text{in } \Omega, \\ \rho \in \partial \beta^*(x, \eta), \quad \phi = \partial_{\xi} F^*(., \nabla \eta) & \\ (\phi - V) \cdot \nu = \pi & \text{on } \Gamma_N, \\ \eta = g & \text{on } \Gamma_D, \\ \rho(0) = \rho_0 & \end{array} \right. $$

In Section 5, we show how to use this analysis to solve this kind of problems in the case where

$$ V(x, \rho, \eta) = \rho V, $$

where $V$ is a given vector field.
## 5.2 Proofs of Theorem 1 and Theorem 2

The proofs of Theorem 1 and Theorem 2 are connected. To simplify the presentation of both proofs, we introduce the problems

$$ \tilde{\mathcal{N}}_{\alpha}(\mu) := \inf_{\rho, \phi} \left\{ \alpha \int_{\Omega} \beta(\cdot, \rho) + \int_{\Omega} F(\cdot, \phi) \, dx - \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D} : \phi \in \mathcal{A}_{\pi}^{\chi}(\mu - \alpha \rho), \rho \in L^{p'}(\Omega) \right\}. $$

and

$$ \tilde{D}_{\alpha} := \max_{\eta \in W^{1,p}(\Omega), \eta/\Gamma_D = g} \left\{ \int_{\Omega} \mu \, \eta - \int_{\Omega} \chi \cdot \nabla \eta - \int_{\Omega} F^*(\cdot, \nabla \eta) - \alpha \int_{\Omega} \beta^*(\cdot, \eta) + \langle \pi, \eta \rangle_{\Gamma_N} \right\}, $$

where $\alpha \in \{0, 1\}$. Then, it is clear that $\mathcal{I}_{\pi,g}^{\chi}(\mu) = \tilde{\mathcal{N}}_0(\mu)$, $\mathcal{N}_{\pi,g}^{\chi}(\mu) = \tilde{\mathcal{N}}_1(\mu)$, $\mathcal{J}_{\pi,g}^{\chi}(\mu) = \tilde{D}_0(\mu)$ and $D_{\pi,g}^{\chi}(\mu) = \tilde{D}_1(\mu)$. So, it is enough to prove that $\tilde{\mathcal{N}}_{\alpha}(\mu) = \tilde{D}_{\alpha}(\mu)$, for $\alpha \in \{0, 1\}$.

To begin with, we prove the following lemma.

**Lemma 2** For any $\mu \in L^{p'}(\Omega)$, $\chi \in L^{p'}(\Omega)^N$ and $\alpha \in \{0, 1\}$, the problems $\tilde{\mathcal{N}}_{\alpha}(\mu)$ and $\tilde{D}_{\alpha}(\mu)$ have solutions.

**Proof:** Let $(\rho_n^\alpha, \phi_n^\alpha)$ be a minimizing sequence for $D(\mu)$. Here, we take $\rho_n^\alpha$ equal to 0 in the case where $\alpha = 0$. Thanks to the assumptions on $F$ and $\alpha\beta$, one sees that $\rho_n$ and $\phi_n$ are bounded in $L^{p'}(\Omega)$ and $L^{p'}(\Omega)^N$, respectively. So, there exists a sub-sequence that we denote again by $(\rho_n^\alpha, \phi_n^\alpha) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$, and $(\rho^\alpha, \phi^\alpha) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$, such that

$$ \rho_n^\alpha \to \rho^\alpha, \quad \text{in } L^{p'}(\Omega)\text{-weak}^* $$

and

$$ \phi_n^\alpha \to \phi^\alpha, \quad \text{in } L^{p'}(\Omega)^N\text{-weak,} $$

where $\rho^0 \equiv 0$. In addition, combining (4.1) with the fact that $\nabla \cdot \phi_n^\alpha = \mu + \nabla \cdot \chi - \rho_n^\alpha$ in $\Omega$ and $(\phi_n^\alpha + \chi) \cdot \nu = \pi$ on $\Gamma_N$, we see that

$$ \langle (\phi_n^\alpha + \chi) \cdot \nu, g \rangle_{\Gamma_D} \to \langle (\phi^\alpha + \chi) \cdot \nu, g \rangle_{\Gamma_D}. $$

Clearly $(\rho^\alpha, \phi^\alpha)$ is an admissible test function for the optimization problem $\tilde{\mathcal{N}}_\alpha(\mu)$; i.e. $(\rho^\alpha, \phi^\alpha) \in \mathcal{A}_{\pi}^{\chi}(\mu - \alpha \rho^\alpha)$. Then, using the l.s.c. and convexity of $\beta$ and $F$, we deduce that $(\rho^\alpha, \phi^\alpha)$ is a solution of the optimization problem $\tilde{\mathcal{N}}_\alpha(\mu)$. The proof for $\tilde{D}_\alpha(\mu)$ follows more or less the same arguments. To avoid redundancy, we will omit the details of the proof. $\square$

Now, for any $\alpha \in \{0, 1\}$, let us consider the application $K_\alpha : W_{\Gamma_D}^{-1,p'}(\Omega) \to \mathbb{R}$ given by

$$ K_{\alpha}(f) := \int_{\Omega} f_0 \tilde{g} - \int_{\Omega} \bar{f} \cdot \nabla \tilde{g} + \inf_{\rho, \phi} \left\{ \alpha \int_{\Omega} \beta(\cdot, \rho) + \int_{\Omega} F(\cdot, \phi) \, dx - \langle (\phi + \chi + \bar{f}) \cdot \nu, g \rangle_{\Gamma_D} \right. \\ \left. : \phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0 - \alpha \rho), \rho \in L^{p'}(\Omega) \right\}, $$
where $f_0$ and $\bar{f}$ are given by the decomposition $f = f_0 + \nabla \cdot \bar{f}$ in $W_{\Gamma_D}^{-1,p'}(\Omega)$. Then, it is clear that

$$ \mathcal{N}_{\pi,g}^{\chi}(\mu) = K_1(0) \quad \text{and} \quad \mathcal{I}_{\pi,g}^{\chi}(\mu) = K_0(0). $$

Thank to (4.1), we see that

$$ \langle (\phi + \chi + \bar{f}) \cdot \nu, g \rangle_{\Gamma_D} = \int_{\Omega} (\phi + \chi + \bar{f}) \cdot \nabla \tilde{g} - \int_{\Omega} (f_0 + \mu - \alpha \rho) \tilde{g}, $$

for any $\rho \in L^{p'}(\Omega)$ and $\phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0 - \alpha \rho)$. This implies that $K_{\alpha}$ may be written as

$$ (5.7) \qquad K_{\alpha}(f) = \inf_{\rho, \phi} \left\{ \begin{array}{l} \alpha \int_{\Omega} \beta(\cdot, \rho) + \int_{\Omega} F(\cdot, \phi) \, dx - \int_{\Omega} (\phi + \chi) \cdot \nabla \tilde{g} \, dx + \int_{\Omega} (\mu - \alpha \rho) \tilde{g} \\ \qquad : \phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0 - \alpha \rho), \rho \in L^{p'}(\Omega) \end{array} \right\}. $$

The following duality result is the key to proving Theorem 1 and Theorem 2. After, we need simply to compute the Legendre transform of $K_{\alpha}$ to complete the proofs.

**Lemma 3** For any $\alpha \in \{0, 1\}$, we have

$$ (5.8) \qquad -K_{\alpha}(0) = -K_{\alpha}^{**}(0) = \min_{\eta \in W^{1,p}(\Omega)} K_{\alpha}^{*}(\eta), $$

**Proof**: To prove (5.8), it is enough to prove that $K_{\alpha}$ is convex and l.s.c., for any $\alpha \in \{0, 1\}$, and conclude by classical duality results (cf. [26]).

**Convexity**: For any $f, h \in W_{\Gamma_D}^{-1,p'}(\Omega)$, taking $(\rho_f, \phi_f)$ and $(\rho_h, \phi_h)$ the solutions corresponding to the optimization problems $K_{\alpha}(f)$ and $K_{\alpha}(h)$ respectively, one sees that

$$ \underbrace{t\phi_h + (1-t)\phi_f}_{\phi_t} \in \mathcal{A}_{\pi}^{\bar{h}_t}(\mu - \alpha \underbrace{(t\rho_h + (1-t)\rho_f)}_{\rho_t} + \underbrace{(th_0 + (1-t)f_0)}_{h_0t}), $$

where $\bar{h}_t = t\bar{h} + (1-t)\bar{f}$, for any $t \in [0, 1]$. Moreover, using the convexity of $\beta(x, \cdot)$ and $F(x, \cdot)$, we have

$$ \alpha \int \beta(\cdot, \rho_t) + \int_{\Omega} F(\cdot, \phi_t) - \langle (\phi_t + \chi + \bar{f}_t) \cdot \nu, g \rangle_{\Gamma_D} \le t \left( \alpha \int_{\Omega} \beta(\cdot, \rho_h) + \int_{\Omega} F(\cdot, \phi_h) - \langle (\phi_h + \chi + \bar{h}) \cdot \nu, g \rangle_{\Gamma_D} \right) \\ +(1-t) \left( \alpha \int_{\Omega} \beta(\cdot, \rho_f) + \int_{\Omega} F(\cdot, \phi_f) - \langle (\phi_f + \chi + \bar{f}) \cdot \nu, g \rangle_{\Gamma_D} \right), $$

which implies that

$$ \int_{\Omega} (th_0 + (1-t)f_0) \tilde{g} - \int_{\Omega} (t\bar{h} + (1-t)\bar{f}) \cdot \nabla \tilde{g} + \alpha \int_{\Omega} \beta(\cdot, \rho_t) + \int_{\Omega} F(\cdot, \phi_t) - \langle (\phi_t + \chi + \bar{f}_t) \cdot \nu, g \rangle_{\Gamma_D} \\ \le tK_{\alpha}(h) + (1-t)K_{\alpha}(f), $$

and

$$ K_{\alpha}(th + (1-t)f) \le tK_{\alpha}(h) + (1-t)K_{\alpha}(f). $$
**Lower semi-continuity:** Let us consider $f_n$ a sequence of $W_{\Gamma_D}^{-1,p'}(\Omega)$ which converges to $f$. That is a sequence of $L^{p'}(\Omega)$ and $L^{p'}(\Omega)^N$ functions $f_{0n}$ and $\bar{f}_n$ respectively, such that

$$ \int_{\Omega} f_{0n} \xi - \int_{\Omega} \bar{f}_n \cdot \nabla \xi \rightarrow \int_{\Omega} f_0 \xi - \int_{\Omega} \bar{f} \cdot \nabla \xi, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Let us prove that $K_\alpha(f) \le \liminf_{n \to \infty} K_\alpha(f_n)$. Then, let us consider $(\rho_n, \phi_n)$ be the solution corresponding to $K_\alpha(f_n)$; i.e.

$$ (5.9) \qquad \phi_n \in \mathcal{A}_{\pi}^{\chi+\bar{f}_n}(\mu + f_{0n} - \alpha \rho_n). $$

and

$$ K_{\alpha}(f_n) = \alpha \int_{\Omega} \beta(., \rho_n) + \int_{\Omega} F(., \phi_n) dx - \int_{\Omega} (\phi_n + \chi) \cdot \nabla \tilde{g} dx + \int_{\Omega} (\mu - \alpha \rho_n) \tilde{g}, $$

where we use (5.7). We can assume that $K_\alpha(f_n)$ is bounded. So, there exists $C < \infty$ such that

$$ \alpha \int_{\Omega} \beta(., \rho_n) + \int_{\Omega} F(., \phi_n) dx \le C + \int_{\Omega} (\phi_n + \chi) \cdot \nabla \tilde{g} dx - \int_{\Omega} (\mu - \alpha \rho_n) \tilde{g} $$

Using assumptions (H1) and (H4) with Young formula, we see that $\rho_n$ and $\phi_n$ are bounded in $L^{p'}(\Omega)$ and $L^{p'}(\Omega)^N$, respectively. So, there exists $(\rho, \phi) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$ and a sub-sequence that we denote again by $(\rho_n, \phi_n)$, such that

$$ \rho_n \to \rho, \quad \text{in } L^{p'}(\Omega)\text{-weak}^* $$

and

$$ \phi_n \to \phi, \quad \text{in } L^{p'}(\Omega)^N\text{-weak.} $$

Using (5.9), we have

$$ \int_{\Omega} (\phi_n - \chi) \cdot \nabla \xi = \int_{\Omega} (\mu - \rho_n) \xi + \langle f_n, \xi \rangle_{\Omega} + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}, $$

and, by letting $n \to \infty$, we get

$$ \int_{\Omega} (\phi - \chi) \cdot \nabla \xi = \int_{\Omega} (\mu - \rho) \xi + \langle f, \xi \rangle_{\Omega} + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}. $$

This implies that $\phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0 - \alpha \rho)$. Then, using the l.s.c. and convexity of $\beta$ and $F$, we have

$$ \begin{align*} K_{\alpha}(f) & \le \alpha \int_{\Omega} \beta(., \rho) + \int_{\Omega} F(., \phi) dx - \int_{\Omega} (\phi + \chi) \cdot \nabla \tilde{g} dx + \int_{\Omega} (\mu - \alpha \rho) \tilde{g} \\ & \le \liminf_{n \to \infty} \left\{ \alpha \int_{\Omega} \beta(., \rho_n) + \int_{\Omega} F(., \phi_n) dx - \int_{\Omega} (\phi_n + \chi) \cdot \nabla \tilde{g} dx + \int_{\Omega} (\mu - \alpha \rho_n) \tilde{g} \right\} \\ & = \liminf_{n \to \infty} K_{\alpha}(f_n). \end{align*} $$

Thus the result.

□
**Proof of Theorem 1** : As we said above, it is enough to compute $K_0^*$. By definition of $K_0^*$, for any $\eta \in W_{\Gamma_D}^{1,p}(\Omega)$, we have

$$
\begin{aligned}
K_0^*(\eta) &= \max \left\{ \langle f, \eta \rangle_{\Omega} - K_0(f) : f \in W_{\Gamma_D}^{-1,p'}(\Omega) \right\} \\
&= \max_{f_0, \bar{f}, \phi} \left\{ \underbrace{\int_{\Omega} f_0 \eta - \int_{\Omega} \bar{f} \cdot \nabla \eta - \int_{\Omega} F(\cdot, \phi) + \int_{\Omega} (\phi + \chi) \cdot \nabla \tilde{g} \, dx - \int_{\Omega} \mu \tilde{g}}_{H(f_0, \bar{f}, \phi)} \right. \\
&\qquad \left. : (f_0, \bar{f}) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N, \phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0) \right\}.
\end{aligned}
$$

Since $\phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0)$ and $\eta \in W_{\Gamma_D}^{1,p}(\Omega)$, we have

$$ \int_{\Omega} f_0 \eta - \int_{\Omega} \bar{f} \cdot \nabla \eta = \int_{\Omega} (\phi + \chi) \cdot \nabla \eta - \int_{\Omega} \mu \eta - \langle \pi, \eta \rangle_{\Gamma_N}. $$

This implies that

$$ H(f_0, \bar{f}, \phi) = \int_{\Omega} (\phi \cdot \nabla(\eta + \tilde{g}) - F(\cdot, \phi)) - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N}. $$

for any $(f_0, \bar{f}) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$ and $\phi \in \mathcal{A}_{\pi}^{\chi+\bar{f}}(\mu + f_0)$, so that

$$ K_0^*(\eta) = \max_{\phi \in L^{p'}(\Omega)^N} \left\{ \int_{\Omega} (\phi \cdot \nabla(\eta + \tilde{g}) - F(\cdot, \phi)) - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N} \right\}. $$

Under the assumptions (H1), we deduce

$$ K_0^*(\eta) = \int_{\Omega} F^*(\cdot, \nabla(\eta + \tilde{g})) - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N}. $$

Then by using (5.8), we conclude that

$$
\begin{aligned}
-K_0(0) &= \min_{\eta \in W_{\Gamma_D}^{1,p}(\Omega)} \int_{\Omega} F^*(\cdot, \nabla(\eta + \tilde{g})) - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N} \\
&= \min_{\eta \in W^{1,p}(\Omega), \eta/\Gamma_D = g} \int_{\Omega} F^*(\cdot, \nabla\eta) - \int_{\Omega} \mu\eta + \int_{\Omega} \chi \cdot \nabla\eta - \langle \pi, \eta \rangle_{\Gamma_N}.
\end{aligned}
$$

Thus the duality

$$ \mathcal{I}_{\pi,g}^{\chi}(\mu) = \mathcal{J}_{\pi,g}^{\chi}(\mu). $$

Thanks to Lemma 2, let us consider $\phi$ and $\eta$ the solutions of $\mathcal{I}_{\pi,g}^{\chi}(\mu)$ and $\mathcal{J}_{\pi,g}^{\chi}(\mu)$. Our aim now is to prove that the couple $(\eta, \phi)$ solve the PDE (5.5). The duality (5.4) implies

$$ \int_{\Omega} \mu \eta - \int_{\Omega} \chi \cdot \nabla \eta - \int_{\Omega} F^*(\cdot, \nabla \eta) + \langle \pi, \eta \rangle_{\Gamma_N} = \int_{\Omega} F(\cdot, \phi) \, dx - \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D}. $$

Since $\phi \in \mathcal{A}_{\pi}^{\chi}(\mu)$, we have

$$ \int_{\Omega} (\phi + \chi) \cdot \nabla \eta = \int_{\Omega} \mu \eta + \langle \pi, \eta \rangle_{\Gamma_N} + \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D}. $$
Combining both equations we deduce that

$$ \int_{\Omega} F(., \phi) dx + \int_{\Omega} F^{*}(., \nabla \eta) dx = \int_{\Omega} \phi \cdot \nabla \eta. $$

Then, by using the fact that $\phi \cdot \nabla \eta \le F^{*}(., \nabla \eta) + F(., \phi)$, we deduce that $F(., \phi) + F^{*}(., \nabla \eta) = \phi \cdot \nabla \eta$ and $\phi \in \partial F^{*}(., \nabla \eta)$ a.e. in $\Omega$. Thus $(\eta, \phi)$ solve the PDE (5.5). For the converse part, it is not difficult to see first that

$$ \mathcal{I}_{\pi,g}^{\chi}(\mu) \geq \mathcal{J}_{\pi,g}^{\chi}(\mu). $$

Then, working with $(\eta, \phi)$ being a solution of PDE (5.5), on proves easily the converse inequality and then concludes $\mathcal{I}_{\pi,g}^{\chi}(\mu) = \mathcal{J}_{\pi,g}^{\chi}(\mu)$. $\square$

**Proof of Theorem 2** : As in the proof of Theorem 1, we only need to compute $K_{1}^{*}$. By definition, for any $\eta \in W_{\Gamma_D}^{1,p}(\Omega)$, we have

$$ K_{0}^{*}(\eta) = \max \left\{ \langle f, \eta \rangle_{\Omega} - K_{0}(f) : f \in W_{\Gamma_D}^{-1,p'}(\Omega) \right\} \\ = \max_{f_0, \bar{f}, \phi} \left\{ \underbrace{\int_{\Omega} f_0 \eta - \int_{\Omega} \bar{f} \cdot \nabla \eta - \int_{\Omega} \beta(., \rho) - \int_{\Omega} F(., \phi) + \int_{\Omega} (\phi + \chi) \cdot \nabla \tilde{g} dx - \int_{\Omega} (\mu - \rho) \tilde{g}}_{H(f_0, \bar{f}, \rho, \phi)} : (f_0, \bar{f}) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N, \phi \in \mathcal{A}_{\pi}^{\chi}(\mu - \rho + f_0, \chi + \bar{f}), \rho \in L^{p'}(\Omega) \right\}. $$

Working as in the proof of Theorem 1, we get

$$ H(f_0, \bar{f}, \rho, \phi) = \int_{\Omega} (\phi \cdot \nabla(\eta + \tilde{g}) - F(., \phi)) + \int_{\Omega} (\rho(\eta + \tilde{g}) - \beta(., \rho)) \\ - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N}. $$

Then we deduce that

$$ K_{1}^{*}(\eta) = \max_{\rho, \phi} \left\{ \int_{\Omega} (\phi \cdot \nabla(\eta + \tilde{g}) - F(., \phi)) + \int_{\Omega} (\rho(\eta + \tilde{g}) - \beta(., \rho)) - \int_{\Omega} \mu(\eta + \tilde{g}) \right. \\ \left. + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N} : (\rho, \phi) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N \right\} \\ = \max_{\phi \in L^{p'}(\Omega)^N} \int_{\Omega} (\phi \cdot \nabla(\eta + \tilde{g}) - F(., \phi)) + \max_{\rho \in L^{p'}(\Omega)} \left\{ \int_{\Omega} (\rho(\eta + \tilde{g}) - \beta(., \rho)) \right. \\ \left. - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N} \right\}. $$

Using the assumptions (H1) and (H4), we deduce that

$$ K_{1}^{*}(\eta) = \int_{\Omega} F^{*}(., \nabla(\eta + \tilde{g})) + \int_{\Omega} \beta^{*}(., \eta + \tilde{g}) - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N}. $$
Combining this with (5.8), we obtain

$$ 
\begin{aligned}
-K_1(0) &= \min_{\eta \in W^{1,p}(\Omega)} \int_{\Omega} F^*(., \nabla(\eta + \tilde{g})) + \int_{\Omega} \beta^*(., \eta + \tilde{g}) - \int_{\Omega} \mu(\eta + \tilde{g}) + \int_{\Omega} \chi \cdot \nabla(\eta + \tilde{g}) - \langle \pi, \eta \rangle_{\Gamma_N} \\
&= \min_{\eta \in W^{1,p}(\Omega), \eta/\Gamma_D = g} \int_{\Omega} F^*(., \nabla\eta) + \int_{\Omega} \beta^*(., \eta) - \int_{\Omega} \mu\eta + \int_{\Omega} \chi \cdot \nabla\eta - \langle \pi, \eta \rangle_{\Gamma_N}.
\end{aligned}
$$

Thus the duality $N_{\pi,g}^{\chi}(\mu) = D_{\pi,g}^{\chi}(\mu)$. Now, thanks to Lemma 2, let us consider $(\rho, \phi)$ and $\eta$ be the solutions of $N_{\pi,g}^{\chi}(\mu) = D_{\pi,g}^{\chi}(\mu)$. We have

$$ 
\begin{aligned}
& \int_{\Omega} \mu \eta - \int_{\Omega} \chi \cdot \nabla \eta - \int_{\Omega} F^*(., \nabla \eta) - \int_{\Omega} \beta^*(., \eta) + \langle \pi, \eta \rangle_{\Gamma_N} \\
&= \int_{\Omega} \beta(., \rho) + \int_{\Omega} F(., \phi) \, dx - \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D}.
\end{aligned}
$$

Since $\phi \in A_{\pi}^{\chi}(\mu - \rho)$, $\eta \in W^{1,p}(\Omega)$ and $\eta = g$ on $\Gamma_D$, we also have

$$ \int_{\Omega} \rho \eta + \int_{\Omega} (\phi + \chi) \cdot \nabla \eta = \int_{\Omega} \mu \eta + \langle \pi, \eta \rangle_{\Gamma_N} + \langle (\phi + \chi) \cdot \nu, g \rangle_{\Gamma_D}. $$

Combining both equation we get

$$ \int_{\Omega} F(., \phi) + \int_{\Omega} F^*(., \nabla \eta) + \int_{\Omega} \beta^*(., \eta) + \int_{\Omega} \beta(., \rho) = \int_{\Omega} \phi \cdot \nabla \eta + \int_{\Omega} \eta \rho. $$

This implies that $F(., \phi) + F^*(., \nabla \eta) = \phi \cdot \nabla \eta$ and $\beta^*(., \eta) + \beta(., \rho) = \rho \eta$ a.e. in $\Omega$. Thus $(\rho, \eta, \phi)$ is a solution of the PDE (5.6). The proof of the converse part follows in the same way of the proof of Theorem 1. Indeed, first one sees directly that $D(\mu) \le N(\mu)$. The, by working with the solution of (5.6) one proves $D(\mu) = N(\mu)$. $\square$

To end up this section, we prove the following results which will be useful for the sequel

**Corollary 1** For any $\mu \in L^{p'}(\Omega)$, $\pi \in W^{-1/p',p'}(\Gamma_N)$ and $g \in W_{00}^{1-1/p,p}(\Gamma_D)$], the problem (5.6) has a solution $(\rho, \eta, \phi)$ in the sense of: $\rho \in L^{p'}(\Omega)$, $\eta \in W^{1,p}(\Omega)$, $\eta = g$ on $\Gamma_D$, $\phi \in L^{p'}(\Omega)^N$, $\phi \in \partial F^*(., \nabla \eta)$, $\rho \in \partial \beta^*(x, \eta)$ a.e. in $\Omega$, and

$$ \int_{\Omega} \rho \xi \, dx + \int_{\Omega} \phi \cdot \nabla \xi \, dx = \int_{\Omega} \mu \xi \, dx - \int_{\Omega} \chi \cdot \nabla \xi - \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Moreover, the triplet $(\rho, \eta, \phi)$ satisfies

$$ (5.10) \qquad \int_{\Omega} \rho (\eta - \tilde{g}) + \int_{\Omega} \phi \cdot \nabla (\eta - \tilde{g}) = \int_{\Omega} \mu (\eta - \tilde{g}) + \int_{\Omega} \chi \cdot \nabla (\eta - \tilde{g}) + \langle \pi, \eta \rangle_{\Gamma_N}. $$

**Proof:** The existence result is a direct consequence of Theorem 1. Let $(\rho, \eta, \phi)$ be a solution of the problem. See that $\eta - \tilde{g}$ is an admissible test function. Thus (5.10). $\square$
## 6 Steepest descent algorithm

Thanks to Section 4, the main application of the optimization problem $\mathcal{I}_{\pi,g}^{\chi}$ and $\mathcal{N}_{\pi,g}^{\chi}$ is to introduce a steepest descent algorithm to study the evolution problem of the type

$$ (6.1) \qquad \left\{ \begin{array}{ll} \partial_t \rho - \nabla \cdot (\phi + \rho V) = f, & \text{in } Q, \\ \\ \phi = \partial_\xi F^*(., \nabla \eta), \rho \in \partial \beta^*(x, \eta) & \\ \\ (\phi + \rho V) \cdot \nu = \pi & \text{on } \Sigma_N \\ \\ \eta = g & \text{on } \Sigma_D \\ \\ \rho(0) = \rho_0 & \text{in } \Omega, \end{array} \right. $$

where we assume that $\rho_0 \in L^{p'}(\Omega)$, with $V \in L^\infty(Q)$ and $f = f_0 + \nabla \cdot \bar{f} \in L^{p'}(0, T; W_{\Gamma_D}^{-1, p'}(\Omega))$; i.e. $(f_0, \bar{f}) \in L^{p'}(Q) \times L^{p'}(Q)^N$.

**Definition 1** A weak solution for the problem (6.1) is a triplet $(\rho, \eta, \phi)$ such that $\rho \in L^\infty(0, T; L^{p'}(\Omega)) \cap W^{1, p'}(0, T; W_{\Gamma_D}^{-1, p'}(\Omega))$, $\eta \in L^p(0, T; W^{1, p}(\Omega))$, $\phi \in L^{p'}(Q)$, $\eta = g$ on $\Sigma_D$, $\eta \in \partial\beta(x, \rho)$, $\nabla\eta \in \partial F(x, \phi)$, a.e. in $Q$, and we have

$$ \frac{d}{dt} \int_{\Omega} \rho \xi dx + \int_{\Omega} (\phi + \rho V) \cdot \nabla \xi dx = \langle f, \xi \rangle_{\Omega} + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \text{in } \mathcal{D}'([0, T]), $$

for any $\xi \in W_{\Gamma_D}^{1, p}(\Omega)$, with $\rho(0) = \rho_0$ in $\Omega$.

We proceed by discretizing time with discrete steps of size $\tau$, denoted as $t_0 = 0 < t_1 < \dots < t_n = T$, with $t_i - t_{i-1} = \tau$, for any $i = 1, \dots, n$. Then, we consider the limit as $\tau$ approaches zero in the sequence of piecewise constant curves

$$ (6.2) \qquad \rho^{\tau} = \sum_{i=1}^{n} \rho^{i} \chi_{[t_{i-1}, t_{i}]} \quad \text{with } \rho^{\tau}(0) = \rho_{0}, $$

where $\rho^i$ is given, for each $i = 1, \dots, n$, by

$$ \rho^i = \underset{\rho}{\operatorname{argmin}} \left\{ \int_{\Omega} \beta(., \rho) + \tau \int_{\Omega} F(., \phi/\tau) dx - \tau \langle (\phi - \rho^{i-1}V^i) \cdot \nu, g \rangle_{\Gamma_D} \right. $$

$$ (6.3) \qquad : \quad \phi \in \mathcal{A}_{\pi}^{\tau \rho^{i-1} V^i + \tau \bar{f}^i} (\tau f_0^i + \rho^{i-1} - \rho), \rho \in L^{p'}(\Omega) \right\}. $$

Here

$$ V^i := \frac{1}{\tau} \int_{t_{i-1}}^{t_i} V(t, \cdot) dt, \quad f_0^i(\cdot) := \frac{1}{\tau} \int_{t_{i-1}}^{t_i} f_0(t, \cdot) dt \quad \text{and} \quad \bar{f}^i(\cdot) := \frac{1}{\tau} \int_{t_{i-1}}^{t_i} \bar{f}(t, \cdot) dt \quad \text{a.e. in } \Omega. $$
In terms of $\mathcal{I}_{\pi,g}$, the optimization problem (6.3) may be written as follows

$$ \rho^i = \underset{\rho \in L^{p'}(\Omega)}{\operatorname{argmin}} \left\{ \int_{\Omega} \beta(\cdot, \rho) + \tau \mathcal{I}_{\pi,g} \left( \frac{\tau f_0^i + \rho^i - \rho}{\tau} \right) \right\}. $$

Thanks to Theorem 2, we have

**Corollary 2** The problem (6.3) has a solution $(\rho_i, \phi_i) \in L^{p'}(\Omega) \times L^{p'}(\Omega)^N$. Moreover, there exists $\eta_i \in W^{1,p}(\Omega)$ such that $\eta^i = g$ on $\Gamma_D$, $\rho^i \in \partial\beta^*(x, \eta^i)$, $\phi^i \in \partial F^*(\cdot, \nabla\eta^i)$ a.e. in $\Omega$, and the triplet $(\rho_i, \phi_i, \eta_i)$ satisfies in a weak sense the following system PDE

$$ (6.4) \qquad \left\{ \begin{array}{ll} \rho^i - \tau \nabla \cdot (\phi^i + \rho^{i-1}V^i + \bar{f}^i) = \tau f_0^i + \rho^{i-1}, & \quad \text{in } \Omega, \\ \\ \rho^i \in \partial\beta^*(x, \eta^i), \quad \phi^i \in \partial\xi F^*(\cdot, \nabla\eta^i) & \\ \\ (\phi^i + \rho^{i-1}V^i + \bar{f}^i) \cdot \nu = \pi & \quad \text{on } \Gamma_N, \\ \\ \eta^i = g, & \quad \text{on } \Gamma_D; \end{array} \right. $$

$$ \text{i.e.} \qquad \int_{\Omega} (\rho^i - \rho^{i-1}) \xi \,dx + \tau \int_{\Omega} (\phi^i + \rho^{i-1}V^i + \bar{f}^i) \cdot \nabla \xi \,dx = \tau \int_{\Omega} f_0^i \xi \,dx + \tau \int_{\Gamma_N} \pi \,\xi \,dx, $$ 
for any $\xi \in W_{\Gamma_D}^{1,p}(\Omega)$.

**Proof:** It is enough to apply Theorem 2 with $\tilde{F}_{\tau}(A) = \tau F(A/\tau)$, for any $A \in \mathbb{R}^N$. Then, it is not difficult to see that $\tilde{F}_{\tau}^*(A^*) = \tau F^*(\tau A^*)$, for any $A^* \in \mathbb{R}^N$. Moreover $A \in \partial \tilde{F}_{\tau}^*(A^*)$ if and only if $A^* \in \tau \partial F(A)$. $\square$

To let $\tau \to 0$, we consider the sequences

$$ \eta^{\tau} = \sum_{i=1}^{n} \eta^{i} \chi_{[t_{i-1}, t_{i}]}, \quad \phi^{\tau} = \sum_{i=1}^{n} \phi^{i} \chi_{[t_{i-1}, t_{i}]}, \quad f_{0}^{\tau} = \sum_{i=1}^{n} f^{i} \chi_{[t_{i-1}, t_{i}]}, \quad \bar{f}^{\tau} = \sum_{i=1}^{n} \bar{f}^{i} \chi_{[t_{i-1}, t_{i}]} $$

and

$$ \tilde{\rho}^{\tau}(t) = \frac{(t - t_{i-1})\rho^{i} - (t - t_{i})\rho^{i-1}}{\tau}, \quad \text{for any } t \in [t_{i-1}, t_{i}], \quad i = 1, \dots, n. $$

By construction, the applications $\tilde{\rho}^{\tau}$, $\rho^{\tau}$, $\eta^{\tau}$ and $\phi^{\tau}$ satisfy the following PDE

$$ \left\{ \begin{array}{ll} \partial_t \tilde{\rho}^{\tau} - \nabla \cdot (\phi^{\tau} + \rho^{\tau} \cdot (-\tau)V^{\tau} + \bar{f}^{\tau}) = f_0^{\tau} & \quad \text{in } Q, \\ \\ \rho^{\tau} \in \partial \beta^*(x, \eta^{\tau}), \quad \phi^{\tau} \in \partial \xi F^*(\cdot, \nabla \eta^{\tau}) & \\ \\ (\phi^{\tau} + \rho^{\tau} \cdot (-\tau)V^i + \bar{f}^i) \cdot \nu = \pi & \quad \text{on } \Sigma_N, \\ \\ \eta^{\tau} = g. & \quad \text{on } \Sigma_D, \end{array} \right. $$

where we extend $\rho^{\tau}$ by $\rho_0$ on $(-\infty, 0)$. Our goal now is to prove that, as $\tau \to 0$, the $\tau$-approximation solutions $\tilde{\rho}^{\tau}$, $\rho^{\tau}$, $\eta^{\tau}$ and $\phi^{\tau}$ converge to the solutions of the continuous problem (6.1). We present the main results following this process in the subsequent sections, where we differentiate between different scenarios based on the structural assumptions imposed on $F$ and $\beta$.
**Remark 8** 1. Note that the discrete approximation (6.4) is widely used in the context of nonlinear semigroup theory (at least in the case $V \equiv 0$). The time step approximation $\rho^\tau$ is typically referred to as the Crandall-Liggett approximation (cf. [20]). Its convergence can be attributed to the $L^1$-accretivity of the fundamental operator that drives the dynamic. This approach of "monotonicity" enables to come with strong compactness of $\rho^\tau$ and $\tilde{\rho}^\tau$. We do not use this approach in this paper since the $L^1$-accretive property is hard to obtain for general $V$ and $F$, even if we do believe that it should be a nice issue for the proof of existence (and uniqueness) of a solution. On can see for instance [34] and [35] for linear diffusion with Hele-Shaw and PME non-linearity respectively. Consequently, the approach presented in this paper can be viewed as a fresh perspective on the discrete approximation process typically employed in Crandall-Liggett theory for nonlinear evolution problems in $L^1(\Omega)$.

2. Observe that in the situation where $V \equiv 0$, the steepest descent algorithm through minimum flow we employ in (6.3) directly derives from Euler-Implicit discretization for the evolution problem (6.1). It is worth noting that this is not the case for steepest descent algorithm through Wasserstein metric, even though both methods can capture the same continuous dynamics (at least in the case of homogeneous Neumann boundary conditions and $F$ of $p$-Laplacian type). This observation may be linked to the fact that the linearized regime of Wasserstein distance is closely associated with dual Sobolev norm (cf. Section 7.6 of [53]), and the fact that our approach is intimately connected to dual Sobolev metric gradient flow in the homogeneous case (see Section 7.1).

# 6.1 Main results

## 6.1.1 The case where $V \equiv 0$

Assuming $V \equiv 0$, the evolution problem (1.4) corresponds to

$$ (6.5) \qquad \left\{ \begin{array}{ll} \partial_t \rho - \nabla \cdot (\phi + \bar{f}) = f_0 & \text{in } Q, \\ \rho \in \partial \beta^*(., \eta), \quad \phi \in \partial_\xi F^*(., \nabla \eta) & \\ (\phi + \bar{f}) \cdot \nu = \pi & \text{on } \Sigma_N \\ \eta = g & \text{on } \Sigma_D \\ \rho(0) = \rho_0 & \text{in } \Omega, \end{array} \right. $$

where $(f_0, \bar{f}) \in L^{p'}(Q) \times L^{p'}(Q)^N$.

**Theorem 3** Under the assumptions (H1)-(H4), for any $\rho_0 \in L^{p'}(\Omega)$ and $(f_0, f) \in L^{p'}(Q) \times L^{p'}(Q)^N$, the problem (6.5) has a weak solution $(\rho, \eta, \phi)$.

## 6.1.2 Diffusion-Transport problem with linear diffusion

For linear diffusion, we select quadratic $F$ as follows:

$$ F(x, A) = \frac{k(x)}{2} |A|^2, \quad \text{for a.e. } x \in \Omega \text{ and } A \in \mathbb{R}^N, $$
where $k \in L^\infty(\Omega)$, is such that $0 < \text{Infess}_\Omega(k)$. In this case, the evolution problem reads

$$ (6.6) \qquad \left\{ \begin{array}{ll} \partial_t \rho - \nabla \cdot (k \nabla \eta + \rho V + \bar{f}) = f_0, & \rho \in \partial \beta^*(x, \eta) \quad \text{in } Q, \\ (k \nabla \eta + \rho V + \bar{f}) \cdot \nu = \pi & \text{on } \Sigma_N \\ \eta = g & \text{on } \Sigma_D \\ \rho(0) = \rho_0 & \text{in } \Omega. \end{array} \right. $$

**Theorem 4** Assume (H1) – (H4) are fulfilled. For any $\rho_0 \in L^2(\Omega)$, $V \in L^\infty(Q)$ and $(f_0, f) \in L^2(Q) \times L^2(Q)^N$,

### 6.1.3 Diffusion-Transport problem with nonlinear diffusion

For general $F$, not necessary quadratic as in Theorem 4, we need to work with additional assumptions on $\beta$. For technical reason, we assume in this case the following condition on $\beta$:

(H5) for a.e. $x \in \Omega$, the application $r \in \mathbb{R} \to \beta^*(x, r)$ is differentiable.

**Theorem 5** Under the assumptions (H1) – (H5), for any $\rho_0 \in L^{p'}(\Omega)$, $V \in L^\infty(Q)$ and $(f_0, \bar{f}) \in L^{p'}(\Omega)^N \times L^{p'}(\Omega)$, the problem (6.6) has a weak solution $(\rho, \eta, \phi)$.

## 6.2 Proofs of Theorem 3, Theorem 4 and Theorem 5

Building upon the notation established at the beginning of this section, our primary goal is the (weak) compactness of the sequences $\rho^\tau, \tilde{\rho}^\tau, \eta^\tau$ and $\phi^\tau$. Next, we need to establish the connections between the limits of $\rho^\tau$ and $\eta^\tau$, and also between the limits of $\eta^\tau$ and $\phi^\tau$. The proofs follow from the subsequent lemmas.

**Lemma 4** For any $t \in [0, T)$, we have

$$ (6.7) \qquad \begin{aligned} & \frac{d}{dt} \int_{\Omega} (\beta(., \rho^\tau) - \rho^\tau \tilde{g}) + \int_{\Omega} \phi^\tau \cdot \nabla (\eta^\tau - \tilde{g}) \le \langle f^\tau, \eta^\tau - \tilde{g} \rangle_{\Omega} \\ & - \int_{\Omega} \rho^\tau (t - \tau) V^\tau \cdot \nabla (\eta^\tau - \tilde{g}) + \langle \pi, \eta^\tau \rangle_{\Gamma_N}, \quad \text{in } \mathcal{D}'([0, T]). \end{aligned} $$

**Proof:** Thanks to Corollary 1, we know that

$$ \int_{\Omega} \rho^i (\eta^i - \tilde{g}) + \tau \int_{\Omega} \phi^i \cdot \nabla (\eta^i - \tilde{g}) = \int_{\Omega} \rho^{i-1} (\eta^i - \tilde{g}) + \tau \langle f^i, \eta^i - \tilde{g} \rangle_{\Omega} \\ - \int_{\Omega} \rho^{i-1} V^i \cdot \nabla (\eta^i - \tilde{g}) + \langle \pi, \eta^i \rangle_{\Gamma_N}. $$

Since $(\rho^i - \rho^{i-1})\eta^i \ge \beta(., \rho^i) - \beta(., \rho^{i-1})$, this implies that

$$ \int_{\Omega} (\beta(., \rho^i) - \rho^i \tilde{g}) - \int_{\Omega} (\beta(., \rho^{i-1}) - \rho^{i-1} \tilde{g}) + \tau \int_{\Omega} \phi^i \cdot \nabla (\eta^i - \tilde{g}) \le \tau \langle f^i, \eta^i - \tilde{g} \rangle_{\Omega} \\ - \tau \int_{\Omega} \rho^{i-1} V^i \cdot \nabla (\eta^i - \tilde{g}) + \tau \langle \pi, \eta^i \rangle_{\Gamma_N}. $$
Thus (6.7).

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

**Lemma 5** We have

1. There exists a constant $C = C(p, N, \Omega, \|V\|_{\infty})$ such that :

$$ (6.8) \quad C \left\{ \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} + \int_{0}^{T} \int_{\Omega} |\nabla \eta^{\tau}|^{p} \right\} \le B \left( 1 + \frac{T \|V^{\tau}\|_{\infty}^{p'}}{C} \exp \left( T \|V^{\tau}\|_{\infty}^{p'} / C \right) \right), $$

where $B$ is given by

$$ B = \int_{\Omega} \beta(., \rho_0) + \int_{\Omega} |\rho_0 - M|^{p'} + \int_{\Omega} |\tilde{g}|^p + T \|V\|_{\infty}^{p'} \int_{\Omega} |\nabla \tilde{g}|^p + T \|\pi\|_{W^{-1/p', p'}(\Gamma_N)}^{p'} \\ + \tau \int_{\Omega} (|\rho_0| - M)^{+p} + \sup_{0<\tau<1} \int_{0}^{T} \int_{\Omega} (|f_0^{\tau}|^{p'} + |\bar{f}^{\tau}|^{p'}). $$

2. The sequences $\rho^\tau$ and $\tilde{\rho}^\tau$ are bounded in $L^\infty(0, T; L^p(\Omega))$, $\eta^\tau$ is bounded in $L^p(0, T; W^{1,p}(\Omega))$ and $\phi^\tau$ in bounded in $L^p'(Q)^N$.

3. The sequence $\partial_t \tilde{\rho}^\tau$ is bounded in $L^{p'}(0, T; W_{\Gamma_D}^{-1, p'}(\Omega))$.

**Proof** : Using Lemma 4, the assumption (H1) and (H4), and the fact that $\phi^\tau \cdot \nabla(\eta^\tau - \tilde{g}) \ge F^*(\nabla\eta^\tau) - F^*(\nabla\tilde{g})$, we get

$$ C_3 \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} + C_1 \int_{0}^{T} \int_{\Omega} |\nabla \eta^{\tau}|^{p} \le \int_{\Omega} \beta(., \rho_0) - \int_{\Omega} \rho_0 \tilde{g} + \int_{\Omega} \rho^{\tau}(t) \tilde{g} + T \int_{\Omega} F^*(., \nabla \tilde{g}) \\ + \int_{0}^{T} \langle f^{\tau}, \eta^{\tau} - \tilde{g} \rangle_{\Omega} - \int_{0}^{T} \int_{\Omega} \rho^{\tau}(t-\tau) V^{\tau} \cdot \nabla (\eta^{\tau} - \tilde{g}) + \int_{0}^{T} \langle \pi, \eta^{\tau} \rangle_{\Gamma_N}, $$

for any $t \in [0, T)$, and then

$$ C_3 \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} + C_1 \int_{0}^{T} \int_{\Omega} |\nabla \eta^{\tau}|^{p} \le \int_{\Omega} \beta(., \rho_0) + \int_{\Omega} (M - \rho_0) \tilde{g} + \int_{\Omega} (|\rho^{\tau}(t)| - M)^+ \tilde{g} \\ + T \int_{\Omega} F^*(., \nabla \tilde{g}) + \int_{0}^{T} \langle f^{\tau}, \eta^{\tau} - \tilde{g} \rangle_{\Gamma_N} - \int_{0}^{T} \int_{\Omega} \rho^{\tau}(t-\tau) V^{\tau} \cdot \nabla (\eta^{\tau} - \tilde{g}) + \int_{0}^{T} \langle \pi, \eta^{\tau} \rangle_{\Gamma_N}, $$

for any $t \in [0, T)$. Using Young and Poincaré inequalities, we can find a constant $C = C(p, N, \Omega) > 0$ such that
(6.9)

$$ C \left\{ \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} + \int_{0}^{T} \int_{\Omega} |\nabla \eta^{\tau}|^{p} \right\} \le \int_{\Omega} \beta(., \rho_0) + \int_{\Omega} |\rho_0 - M|^{p'} + \int_{\Omega} |\tilde{g}|^{p} + T \int_{\Omega} |\nabla \tilde{g}|^{p} \\ + T \|\pi\|_{W^{-1/p', p'}(\Gamma_N)}^{p'} + \int_{0}^{T} \int_{\Omega} (|f_0^{\tau}|^{p'} + \|\bar{f}\|^{p'}) + \int_{0}^{T} \int_{\Omega} \rho^{\tau}(t-\tau) V^{\tau} \cdot \nabla (\eta^{\tau} - \tilde{g}), \quad \text{for any } t \in [0, T), $$
where we use also the fact that $F^*(., \nabla \tilde{g}) \le C_2' \|\nabla \tilde{g}\|^p$. On the other hand, we see that the last term satisfies

$$ \int_0^T \int_{\Omega} \rho^{\tau}(t-\tau) V^{\tau} \cdot \nabla(\eta^{\tau} - \tilde{g}) \le \int_0^T \int_{\Omega} (|\rho^{\tau}(t-\tau)| - M)^{+} |V^{\tau}| (|\nabla \eta^{\tau}| + |\nabla \tilde{g}|) \\ + M \int_0^T \int_{\Omega} |V^{\tau}| (|\nabla \eta^{\tau}| + |\nabla \tilde{g}|). $$

Using again Young inequality and the fact that

$$ \int_0^T \int_{\Omega} (|\rho^{\tau}(t-\tau)| - M)^{+p} \le \int_0^T \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p} + \tau \int_{\Omega} (|\rho_0| - M)^{+p}, $$

for any $\epsilon > 0$ small enough, we can find a constant that we denote by $C_{\epsilon} = C_{\epsilon}(p, N, \Omega, \epsilon) > 0$, such that

$$ C_{\epsilon} \int_{\tau}^{T} \int_{\Omega} \rho^{\tau}(t-\tau) V^{\tau} \cdot \nabla(\eta^{\tau} - \tilde{g}) \le \|V^{\tau}\|_{\infty}^{p'} \left\{ \int_{0}^{T} \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} + \epsilon \int_{0}^{T} \int_{\Omega} |\nabla \eta^{\tau}|^{p} \right. \\ \left. + \epsilon T \int_{\Omega} |\nabla \tilde{g}|^{p} + \tau \int_{\Omega} (|\rho_0| - M)^{+p} \right\}. $$

Working with a fixed small $\epsilon$, and combining this with (6.9), we can find a constant $C = C(p, N, \Omega, \|V\|_{\infty}) > 0$, such that

$$ (6.10) \qquad C \left\{ \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} + \int_{0}^{T} \int_{\Omega} |\nabla \eta^{\tau}|^{p} \right\} \le B + \|V^{\tau}\|_{\infty}^{p'} \int_{0}^{T} \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'}. $$

This implies that

$$ C \int_{\Omega} (|\rho^{\tau}(t)| - M)^{+p'} \le B \exp \left( T \|V^{\tau}\|_{\infty}^{p'}/C \right) $$

and then, (6.8) follows by (6.10). As to $\partial_t \tilde{\rho}^{\tau}$, we recall that

$$ \partial_t \tilde{\rho}^{\tau} = f_0^{\tau} + \nabla \cdot (\rho^{\tau}(. - \tau)V^{\tau} + \phi^{\tau} + \bar{f}^{\tau}) $$

with $(\phi^{\tau}+\rho^{\tau}(. - \tau)V^{\tau}+\bar{f}^{\tau}) \cdot \nu = \pi$ on $\Sigma_N$. Using (6.8) and (H2), one sees easily that $\phi^{\tau}+\rho^{\tau}(. - \tau)V^{\tau}+\bar{f}^{\tau}$ is bounded in $L^{p'}(Q)$. This implies that $\partial_t \tilde{\rho}^{\tau}$ is bounded in $L^{p'}(0, T; W_{\Gamma_D}^{-1, p'}(\Omega))$. $\square$

**Lemma 6** There exists $\rho \in L^{\infty}(0, T; L^{p'}(\Omega))$, $\eta \in L^p(0, T; W_{\Gamma_D}^{1, p}(\Omega))$, $\phi \in L^{p'}(Q)^N$, and subsequences that we denote by again $\rho^{\tau}$, $\tilde{\rho}^{\tau}$ and $\eta^{\tau}$, such that

$$ (6.11) \qquad \rho^{\tau} \to \rho, \quad \text{in } L^{p'}(Q) - \text{weak,} $$

$$ \tilde{\rho}^{\tau} \to \rho, \quad \text{in } L^{p'}(Q) - \text{weak,} $$
$$ (6.12) \qquad \eta^{\tau} \to \eta, \quad \text{in } L^p(0, T; W_{\Gamma_D}^{1,p}(\Omega)) - \text{weak} $$

and

$$ \phi^{\tau} \to \phi, \quad \text{in } L^{p'}(Q)^N - \text{weak.} $$

Moreover, we have

1. $\eta \in \partial\beta(., \rho)$, a.e. in $Q$

2. Under the assumption (H5), we have

$$ (6.13) \qquad \rho^{\tau} \to \rho, \quad \text{in } L^{p'}(Q). $$

3. If

$$ (6.14) \qquad \phi \in \partial F^*(., \nabla \eta), \quad \text{a.e. in } Q, $$

then the triplet $(\rho, \eta, \phi)$ is a solution of the problem 6.5.

Thanks to Lemma 5, the weak compactness of the approximations $\tilde{\rho}^{\tau}$, $\rho^{\tau}$, $\eta^{\tau}$ and $\phi^{\tau}$ is clear. Now, in order to connect the limit of $\rho^{\tau}$, and $\eta^{\tau}$ through the non-linearity $\partial\beta$, we use the very weak Aubin's result due to Andreanov and al. (cf. Theorem A.1 in the appendix, the proof may be found in [6] or [48]).

**Proof of Lemma 6** : The weak compactness of $\rho^{\tau}$, $\tilde{\rho}^{\tau}$ and $\eta^{\tau}$ follows clearly by Lemma 5. Moreover, since

$$ \tilde{\rho}^k(t) - \rho^\tau(t) = (t - t_i) \partial_t \tilde{\rho}^k(t), \quad \text{for any } t \in [t_{i-1}, t_i], \quad i = 1, \dots, n. $$

and $\partial_t \tilde{\rho}^{\tau}$ is bounded in $L^{p'}(0, T; W_{\Gamma_D}^{-1,p'}(\Omega)^*)$, we deduce that $\rho^{\tau}$ and $\tilde{\rho}^{\tau}$ have the same $L^{p'}(Q)$-weak limits. Now, having in mind Theorem A.5, we combine (6.11), (6.12) with the fact that $\partial_t \tilde{\rho}^{\tau}$ is bounded in $L^{p'}(0, T; W_{\Gamma_D}^{-1,p'}(\Omega)^*)$ to deduce first that, as $\tau \to 0$, we have

$$ \int_0^T \int_{\Omega} \rho^{\tau} \eta^{\tau} \varphi \,dxdt \to \int_0^T \int_{\Omega} \rho \eta \varphi \,dxdt \quad \text{for any } \varphi \in \mathcal{D}(Q). $$

Then, using the fact that $\eta^{\tau} \in \partial\beta(\rho^{\tau})$, we deduce that

$$ \eta \in \partial\beta(., \rho), \quad \text{a.e. in } Q. $$

Thanks to (6.7) and (6.8), we see that $\int F^*(., \phi^{\tau})$ is bounded. Combining this with the assumption (H2), we deduce that $\phi^{\tau}$ bounded in $L^{p'}(Q)$, and (6.13) follows. Remember that, for any $t \in (0, T)$ and $\xi \in W_{\Gamma_D}^{1,p}(\Omega)$, we have

$$ \int_{\Omega} \partial_t \tilde{\rho}^{\tau}(t) \xi + \int_{\Omega} \phi^{\tau}(t) \cdot \nabla \xi = \langle f^{\tau}(t), \xi \rangle_{\Omega} - \int_{\Omega} (\rho^{\tau}(t - \tau) V^{\tau}(t) \cdot \nabla \xi + \langle \pi, \xi \rangle_{\Gamma_N}, $$
and then

$$ (6.15) \qquad \frac{d}{dt} \int_{\Omega} \tilde{\rho}^{\tau}(t) \xi + \int_{\Omega} (\phi^{\tau}(t) + \rho^{\tau}(t) V^{\tau}(t)) \cdot \nabla \xi = \langle f^{\tau}(t), \xi \rangle + \langle \pi, \xi \rangle_{\Gamma_N} \\ \qquad \qquad + \underbrace{\int_{\Omega} (\rho^{\tau}(t) - \rho^{\tau}(t-\tau)) V^{\tau}(t) \cdot \nabla \xi}_{J_{\tau}(t)}, \quad \text{in } \mathcal{D}'([0,T)). $$

Using the fact that $\rho^\tau$ is bounded in $L^{p'}(Q)$ and $V^\tau$ is relatively compact in $L^p(Q)$, we see that $J_\tau \to 0$ weakly. So, passing to the limit in (6.15), we deduce that the triplet $(\rho, \eta, \phi)$ satisfies $\rho \in L^\infty(0, T; L^{p'}(\Omega))$, $\partial_t \rho \in W^{1,p'}(0, T; W_{\Gamma_D}^{-1,p'}(\Omega))$, $\rho(0) = \rho_0$ in $\Omega$, $\eta \times L^p(0, T; W^{1,p}(\Omega))$, $\eta = g$ on $\Sigma_D$, $\eta \in \partial\beta(x, \rho)$, a.e. in $Q$, $\phi \times L^{p'}(Q)$ and we have

$$ \frac{d}{dt} \int_{\Omega} \rho \xi \,dx + \int_{\Omega} (\phi + \rho V) \cdot \nabla \xi \,dx = \langle f, \xi \rangle_{\Omega} + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \text{in } \mathcal{D}'([0, T)), \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

At last it is clear if $\phi = \partial F(\cdot, \nabla \eta)$, then the triplet is definitely a solution of the problem (6.1) in a weak sense. $\square$

**Proof of Theorem 4** : Clearly, the proof of Theorem 4 follows by Lemma 6. Indeed, since in this case $\phi^\tau = k(x) \nabla \eta^\tau$ and $\nabla \eta^\tau \to \nabla \eta$ in $L^2(Q)$-weak, then (6.14) is obviously fulfilled. $\square$

The main difficulty in proving Theorems 3 and 5 lies in the justification of $\phi \in \partial F(x, \nabla \eta)$. To address this, we employ a Minty-type argument and establish a sufficient condition that encompasses weak/strong convergences of $\rho^\tau$, $\nabla \eta^\tau$ and $V^\tau$.

**Lemma 7** Under the assumptions of Lemma 6, if

$$ (6.16) \qquad \int_0^T \int_{\Omega} \rho V \cdot \nabla \eta \, dtdx \le \limsup_{\tau \to 0} \int_0^T \int_{\Omega} \rho^\tau(t-\tau) V^\tau \cdot \nabla \eta^\tau \, dtdx, $$

then (6.14) is fulfilled.

**Proof** : Having in mind $\phi^\tau \in \partial F^*(\cdot, \nabla \eta^\tau)$ and the weak convergence of $\phi^\tau$ and $\nabla \eta^\tau$ in $L^{p'(Q)^N}$, we use standard monotonicity arguments of Minty's type to prove this lemma. Indeed, it is enough to prove that $\liminf_{\tau \to 0} \iint_Q \phi^\tau \cdot \nabla (\eta^k - \eta) \le 0$. This is equivalent to prove that

$$ (6.17) \qquad \liminf_{\tau \to 0} \iint_Q \phi^\tau \cdot \nabla \eta^k \, dtdx \le \iint_Q \phi \cdot \nabla \eta \, dtdx. $$

First, one sees that letting $\tau \to 0$, in (6.7), and using the l.s.c. of the convex function $\beta$ with respect to the $L^{p'}(Q)$-weak topology, we have

$$ \liminf_{\tau \to 0} \int_0^T \int_{\Omega} \phi^\tau \cdot \nabla \eta^\tau \le \int_0^T \int_{\Omega} \phi \cdot \nabla \tilde{g} - \int_{\Omega} \beta(\cdot, \rho(t)) + \int_{\Omega} \rho(t) \tilde{g} + \int_{\Omega} \beta(\cdot, \rho_0) - \int_{\Omega} \rho_0 \tilde{g} \\ + \int_0^T \langle f, \eta - \tilde{g} \rangle_{\Omega} + \int_0^T \langle \pi, \eta \rangle_{\Gamma_N} - \limsup_{\tau \to 0} \int_0^T \int_{\Omega} \rho^\tau(t-\tau) V^\tau \cdot \nabla (\eta^\tau - \tilde{g}). $$
Assuming (6.16), we obtain in one hand

$$ (6.18) \qquad \liminf_{\tau \to 0} \int_0^T \int_\Omega \phi^\tau \cdot \nabla \eta^\tau \le \int_0^T \int_\Omega \phi \cdot \nabla \tilde{g} - \int_\Omega \beta(\cdot, \rho(t)) + \int_\Omega \rho(t) \tilde{g} + \int_\Omega \beta(\cdot, \rho_0) - \int_\Omega \rho_0 \tilde{g} \\ \qquad \qquad \qquad + \int_0^T \langle f, \eta - \tilde{g} \rangle_\Omega + \int_0^T \langle \pi, \eta \rangle_{\Gamma_N} - \int_0^T \int_\Omega \rho V \cdot \nabla (\eta - \tilde{g}). $$

On the other hand, since $\partial_t \rho = \nabla \cdot (\phi + \rho V + \bar{f}) + f_0$ in $Q$ and $(\phi + \rho V + \bar{f}) \cdot \nu = \pi$ on $\Sigma_N$, by using Lemma A.7, we have

$$ \int_0^T \int_\Omega \phi \cdot \nabla \eta = \int_0^T \int_\Omega \phi \cdot \nabla \tilde{g} - \int_\Omega \beta(\cdot, \rho(t)) + \int_\Omega \rho(t) \tilde{g} + \int_\Omega \beta(\cdot, \rho_0) - \int_\Omega \rho_0 \tilde{g} \\ \qquad + \int_0^T \langle f, \eta - \tilde{g} \rangle_\Omega + \int_0^T \langle \pi, \eta \rangle_{\Gamma_N} - \int_0^T \int_\Omega \rho V \cdot \nabla (\eta - \tilde{g}). $$

Combining this with (6.18), we deduce (6.17). Thus the result of the lemma. □

**Proof of Theorem 3** : Since we are assuming here that $V \equiv 0$, (6.16) is obviously true and the proof follows directly from Lemma 7. □

**Proof of Theorem 5** : Using Lemma 6.13, we know that under the assumption (H5), $\rho^\tau$ converges in $L^{p'}(Q)$. This implies that (6.16) is fulfilled. Indeed, it is enough to rewrite

$$ \int_0^T \int_\Omega \rho^\tau(t-\tau) V^\tau \cdot \nabla \eta^\tau = \underbrace{\int_0^T \int_\Omega (\rho^\tau(t-\tau) - \rho^\tau(t))V^\tau \cdot \nabla \eta^\tau}_{J_\tau^1} + \underbrace{\int_0^T \int_\Omega \rho^\tau(t)V^\tau \cdot \nabla \eta^\tau}_{J_\tau^2(t)}. $$

Then using the fact that $V^\tau$ and $\nabla \eta^\tau$ are bounded in $L^\infty(Q)$ and $L^p(Q)$, respectively, we have $J_\tau^1 \to 0$, and

$$ \lim_{\tau \to 0} J_\tau = \int_0^T \int_\Omega \rho V \cdot \nabla \eta. $$

# 7 Comments, remarks and extensions

## 7.1 $W_{\Gamma_D}^{-1,p'}(\Omega)$-Gradient flow and $p'$-curves of maximal slope

Despite the fact that (6.4) coincides with Euler-Implicit discretization in $L^1(\Omega)$, the primary differential operator driving the dynamic no longer adheres to the definition of a conventional subgradient operator and cannot be linked to any $L^q$-minimizing energy process (despite quadratic $\beta$). However, we show here formally that, under specific conditions, this operator may be interpreted as a metric subgradient operator within the dual Sobolev space $W_{\Gamma_D}^{-1,p'}(\Omega)$ equipped with its natural distance.
Assume that

$$F(x, A) = \frac{1}{p'} |A|^{p'}, \quad \text{for any } A \in \mathbb{R}^N,$$

with $1 < p' < \infty$, and moreover we assume that $g \equiv 0$ and $\pi \equiv 0$. It is clear in this case that the application

$$f := f_0 + \nabla \cdot \bar{f} \in W_{\Gamma_D}^{-1, p'}(\Omega) \to \mathcal{I}_{0,0}^{\bar{f}}(f_0)^{1/p'},$$

coincides with the usual norm in $W_{\Gamma_D}^{-1, p'}(\Omega)$. Yet, one sees that in general the application $\mathcal{I}_{\pi,g}^{\chi}$ is not connected to a norm, unless one assumes some additional assumptions on $F$. Thanks to Theorem 1, for any $f \in W_{\Gamma_D}^{-1, p'}(\Omega)$, there exists $\phi_f \in L^{p'}(\Omega)^N$, given by

$$\phi_f = \operatorname{argmin} \left\{ \frac{1}{p'} \int_{\Omega} |\phi|^{p'} : \phi \in L^{p'}(\Omega)^N, -\nabla \cdot \phi = f_0 + \nabla \cdot \bar{f} \text{ in } \Omega \text{ and } (\phi + \bar{f}) \cdot \nu = 0 \text{ on } \Gamma_N \right\},$$

such that

$$\|f\|_{W_{\Gamma_D}^{-1, p'}(\Omega)}^{p'} = \frac{1}{p'} \int_{\Omega} |\phi_f|^{p'} dx,$$

and

$$\phi_f = |\nabla u|^{p-2} \nabla u \quad \text{a.e. in } \Omega,$$

where $u$ is the solution of the $p$-Laplacian equation with mixed boundary condition

$$\left\{
\begin{array}{ll}
-\Delta_p u = f_0 + \nabla \cdot \bar{f} & \text{in } \Omega \\
(\phi + \bar{f}) \cdot \nu = 0 & \text{on } \Gamma_N \\
u = 0 & \text{on } \Gamma_D.
\end{array}
\right.$$

In addition, one can prove that the duality bracket $\langle \cdot, \cdot \rangle_{\Omega}$ coincides with

$$\langle f, z \rangle_{\Omega} = \int_{\Omega} \phi_f \cdot \nabla z \, dx \quad \text{for any } f \in W_{\Gamma_D}^{-1, p'}(\Omega) \text{ and } z \in W_{\Gamma_D}^{1, p}(\Omega).$$

It is clear that if $p' = 2$, $W_{\Gamma_D}^{-1,2}(\Omega)$ is Hilbert space and the approach we develop in this paper may be connected to the earlier works for quadratic $F$, [21, 22, 40], in the contest of gradient flow in $H^{-1}$-Hilbert spaces (one can see also the books [14, 8] for the general theory). For the case where $p \neq 2$, the natural duality bracket $\langle \cdot, \cdot \rangle_{\Omega}$ do not allow to treat doubly nonlinear PDEs (6.5) in the contest of gradient flow. To introduce $W_{\Gamma_D}^{-1, p'}(\Omega)$-Gradient flow for (6.5), we consider the metric framework where $\mathcal{X} := W_{\Gamma_D}^{-1, p'}(\Omega)$ is equipped with its natural distance

$$d(f_1, f_2) = \|f_1 - f_2\|_{W_{\Gamma_D}^{-1, p'}(\Omega)}.$$

Then, under the assumptions of Section 4, we consider the application (internal energy) $\mathcal{E}: \mathcal{X} \to (-\infty, +\infty]$ given by

$$\mathcal{E}(\rho) = \begin{cases} \int_{\Omega} \beta(\cdot, \rho) \, dx & \text{if } \beta(\cdot, \rho) \in L^1(\Omega) \\ +\infty & \text{otherwise.} \end{cases}$$
Revisiting the steepest distance algorithm for the PDE (5.5) with $f \equiv 0$, we consider again the time $\tau$-discretization $t_0 = 0 < t_1 < \dots < t_n = T$, with $t_i - t_{i-1} = \tau > 0$, for any $i = 1, \dots, n$. Remember that the piecewise constant curves (6.2) is given by

$$ \rho^i = \underset{\rho \in L^{p'(\Omega)}}{\operatorname{argmin}} \left\{ \int_{\Omega} \beta(\cdot, \rho) + \tau \mathcal{I}_{0,0}^0 \left( \frac{\rho^{i-1} - \rho}{\tau} \right) \right\}. $$

Then, it is clear that

$$ \rho^i = \underset{\rho \in X}{\operatorname{argmin}} \left\{ \mathcal{E}(\rho) + \tau^{p'-1} d(\rho^i, \rho)^{p'} \right\}, $$

which corresponds to the $p'$-scheme for the approximation of $p'$-curves of maximal slope for $\mathcal{E}$ (cf. [4]). So, letting $\tau \to 0$ in the discrete curves

$$ \rho^{\tau} = \sum_{i=1}^{n} \rho^{i} \chi_{[t_{i-1}, t_{i}]} \quad \text{with } \rho^{\tau}(0) = \rho_{0}, $$

we cover the $p'$-curve of maximal slope for the energy $\mathcal{E}$. This implies in particular that the PDE (6.5) may be interpreted as a "metric" gradient flow of $\mathcal{E}$ in $(X, d)$; i.e.

$$ \frac{d\rho}{dt} + \partial_{\mathcal{X}} \mathcal{E}(\rho) = 0 \quad \text{in } (0, T), $$

in $(\mathcal{X}, d)$. The solution may be interpreted in turn as the $p'$-curve of maximal slope in $(\mathcal{X}, d)$ for the energy $\mathcal{E}$ (cf. [4]). That is, $\rho$ satisfies the $p'$-energy identity

$$ (7.1) \qquad \frac{d}{dt} \mathcal{E}(\rho) + \frac{1}{p'} |\rho'|^p + \frac{1}{p} |\partial \mathcal{E}|^p(\rho) = 0 \quad \text{in } \mathcal{D}'([0, T]), $$

where $|\partial \mathcal{E}|$ denotes the slopes of $\mathcal{E}$ and $|\rho'|$ the metric derivative of $\rho$. In particular, Theorem 3 gives an interpretation in terms of PDE of the $p'$-curve of maximal slope given by (7.1), and vice versa. Moreover, one can associate the stationary PDE

$$ \left\{ \begin{array}{ll} -\Delta_p \eta = f, & \rho \in \partial \beta^*(x, \eta) \quad \text{in } \Omega, \\ |\nabla \eta|^{p-2} \nabla \eta \cdot \nu = 0 & \text{on } \Gamma_N \\ \eta = 0 & \text{on } \Gamma_D \end{array} \right. $$

to the metric gradient flow of the internal energy $\int_{\Omega} \beta(\cdot, \rho)$ in $W_{\Gamma_D}^{-1, p'}(\Omega)$. It is worth noting that, to the best of our knowledge, this type of results is well-established for quadratic $F$ with $p = 2$ and homogeneous boundary, but has not been addressed for other cases. We believe that the approach developed in this paper could shed more light on the $p'$-curve of maximal slope solutions for such cases. We do not go into depth on this direction in order to avoid making the paper much more longer.
## 7.2 Prediction-correction algorithm

As we said in the introduction the dynamic in (6.1) results from the superposition of two processes: a transport process and a nonlinear energy minimization process. The transport process represents the external forces acting on the system, while the nonlinear energy minimization process represents the system's internal dynamics. Remember that the Wasserstein steepest descent framework enables to handle the transport process within the internal energy and keeps the transition work out of it. For instance, in the case where $V = \nabla \varphi$, the solution of diffusion-transport porous medium equation

$$ \partial_t \rho - \nabla \cdot (\rho^m - \rho V) = 0, $$

can be obtained through $W_2$-Wasserstein steepest descent framework, by working with $W_\tau(\rho, \tilde{\rho}) = \frac{1}{2\tau}W_2(\rho, \tilde{\rho})^2$, and employing the internal energy

$$ \mathcal{E}(\rho) = \frac{1}{m-1} \int_{\Omega} \rho^m dx + \int_{\Omega} \rho \varphi dx. $$

In our approach, the transport process is incorporated within the traffic cost entity $\mathcal{I}$. Indeed, we retain the same internal energy $\mathcal{E}(\rho) = \frac{1}{m+1} \int_{\Omega} \rho^{m+1} dx$ as in the case without a transport $V$, but the transition cost is modified to a $\chi$-forced one $\mathcal{I}^\chi$, where $\chi$ involves $\rho V$. Remember that $\mathcal{I}$ involves also boundary conditions $g$ and $\pi$ which can be considered also as boundary forced terms processes.

An alternative way to address the diffusion-transport dynamics in Equation (6.1) which may be used successfully is connected to the so-called prediction-correction algorithm (see some variants in [44, 27] in the contest of crowed motion and the references therein). The idea consists initially separates the dynamic into its primary processes: transport and diffusion. In the first step, the transport process is applied to the system through the transport equation $\partial_t \rho + \nabla \cdot (\rho V) = 0$. And in the second step, one can address the diffusion process. For this second step, we can use the proximal minimizing internal energy process as given in Section (5) within a transition cost $\mathcal{I}$ kick without the forcing term $V$.

More precisely, let us see how the algorithm turns on formally for the evolution dynamic (6.1). Again, we consider $T > 0$ a given time horizon, and for a given time step $\tau > 0$, we consider a uniform partition of $[0, T]$ given by $t_k = k\tau$, $k = 0, \dots, n-1$. Supposing that we know the density of the population $\rho_k$ at a given step $k$, starting by $\rho_0$. Then, we superimpose successively the following two steps:

**Prediction:** In the predictive time step $[t_k, t_{k+\frac{1}{2}}]$, where $t_{k+\frac{1}{2}} = k\tau + \tau/2$, the density of population trends to move into

$$ \rho^{k+\frac{1}{2}} = \rho(t_{k+\frac{1}{2}}), $$

where $\rho$ is the solution of the transport equation

$$ \partial_t \rho + \nabla \cdot (\rho V) = 0 \quad \text{in } [t_k, t_{k+\frac{1}{2}}), $$

with $\rho(t_k) = \rho^k$.

**Correction:** In general it is not expected that $\rho_{k+\frac{1}{2}}$ to be an allowable since the forcing term may increase the internal energy. Then, one can proceed minimizing instantaneously the energy
against the transition work $\mathcal{I}_{\pi,g}^{\chi}$, with $\chi \equiv \tau \bar{f}^k$ (keeping $\chi$ out of $V$). More precisely, one considers $\rho_{k+1}$ given by the following optimization problem

$$ \rho_{k+1} = \underset{\rho}{\operatorname{argmin}} \left\{ \int_{\Omega} \beta(., \rho) + \tau \int_{\Omega} F(., \phi/\tau) dx - \tau \langle \phi \cdot \nu, g \rangle_{\Gamma_D} \right. \\ \left. \qquad : \phi \in \mathcal{A}_{\pi}^{\tau \bar{f}^k} (\tau f_0^k + \rho^{k+1/2} - \rho], \rho \in L^{p'}(\Omega) \right\}. $$

Thanks to Theorem 2, one sees that $\rho_{k+1}$ is the solution of the PDE

$$ \left\{ \begin{array}{ll} \rho^{k+1} - \tau \nabla \cdot (\phi^{k+1} + \bar{f}^k) = \tau f_0^k + \rho^{k+1/2}, & \text{in } \Omega, \\ \\ \rho^{k+1} \in \partial \beta^*(x, \eta^{k+1}), \quad \phi^{k+1} \in \partial_{\xi} F^*(., \nabla \eta^{k+1}) & \\ \\ (\phi^{k+1} + \bar{f}^k) \cdot \nu = \pi & \text{on } \Gamma_N, \\ \\ \eta^{k+1} = g, & \text{on } \Gamma_D. \end{array} \right. $$

As in the algorithm of Section (6), the main future is to prove that the approximation

$$ \rho^{\tau} = \sum_{k=0}^{N-1} \rho^{k+1/2} \chi_{[t_k, t_{k+1/2}]} + \rho^{k+1} \chi_{[t_{k+1/2}, t_{k+1}]} \quad \text{with } \rho^{\tau}(0) = \rho_0, $$

converges to the solution of (6.1). This type of algorithm is used in [27] for the study of crowd motion. Other version based on $W_2$-Wasserstein steepest descent algorithm may be found in [44, 45, 46, 47].

Both the algorithm presented in Section 6 and the prediction-correction algorithm differ slightly in their approach. Consider, for instance, a system of particles in a fluid, where the particles are driven by the fluid's velocity and simultaneously attempt to minimize their energy by transitioning to a lower energy state. While the first algorithm focuses on minimizing the energy state by incorporating the force term into the transition cost, the second algorithm aims to first propel the system with the fluid's velocity and then guide it towards a lower energy state using a diffusive transition cost. In essence, the second algorithm employs a prediction-correction strategy, whereas the first algorithm can be considered a "transport forcing transition" type.

At last, it is worth mentioning that we employ a "transport forcing transition" type algorithm in this paper to address Equation (6.1) due to its ability to simultaneously recover compactness and effectively enforce Neumann boundary conditions consistent with those of (6.1) without the need for additional constraints on $V$ at the boundary. Otherwise, careful consideration must be given to the Neumann boundary condition when utilizing the prediction-correction algorithm described above.

## 7.3 Uniqueness

The uniqueness of the solution for the PDE (6.1) is a very interesting question that can quickly become very complicated depending on the assumptions on $\beta$, $F$ and $V$. In the case where $V \equiv 0$, one might expect the uniqueness either through convexity arguments for the optimization problem or by using
PDE techniques to prove $L^1$-contraction principle for the solutions. However, this is no longer true for more general $V$ and the question of uniqueness becomes challenging even for some concrete choice of $F$ or $\beta$ (see for instance Remark 2 in [34]) with respect to the assumptions on $V$ as well as the choice of boundary conditions. To the best of our knowledge, except particular cases involving linear type diffusion (cf. [34, 41, 23, 47, 15, 35] for linear diffusion and [38] for the $p$-Laplacian operator), the question of uniqueness of weak solutions with nonlinear diffusion and general non linearity is not addressed in the literature. We will not treat these questions in this paper.

## 7.4 Linear growth cost $F$: granular diffusion

The case where $F$ satisfies (H1) with $p' = 1$ and $p = \infty$, appears when one deals with granular diffusion phenomena like in sandpile or crowd motion (cf. [27]). In this case, one needs to use vector valued fields whose divergence are Radon measures (cf. [17]). One can see the papers [28, 29] and [36, 37] for the study of some particular cases connected to $I_{\pi,g}^{\chi}$. The steepest descent algorithm of the type $\mathcal{N}$ has been used in [27] with linear growth $F$ and Hele-Shaw type nonlinearity to model congestion phenomena in crowed motion. Formally, the evolution problem of the type (6.1) associated with linear growth $F$, may be given by

$$ 
\left\{ 
\begin{array}{ll} 
\partial_t \rho - \nabla \cdot (\phi + \rho V) = f, \\ 
\rho \in \partial \beta^*(x, \eta), \phi \in \partial \Pi_{B_F}(\nabla \eta) \\ 
(\phi + \rho V) \cdot \nu = \pi \\ 
\eta = g \\ 
\rho(0) = \rho_0 
\end{array} 
\right. 
\quad \text{in } Q, \\ 
\quad \text{on } \Sigma_N \\ 
\quad \text{on } \Sigma_D \\ 
\quad \text{in } \Omega, 
$$

where

$$ B_{F(x, \cdot)} := \{A \in \mathbb{R}^N : F^*(x, A) \le 1\}, \quad \text{for a.e. } x \in \Omega. $$

For instance, if $F(x, A) = k(x) |A|$, for any $x \in \Omega$ and $A \in \mathbb{R}^N$, $B_{F(x, \cdot)}$ coincides with the closed bull $\overline{B}(0, \mathbb{R}^N, k(x))$. Yet, to be well posed one needs to assume some necessary conditions on $g$ and $\pi$. Moreover, the involved flux $\phi$ is no more a Lebesgue integrable function (see [25] and the references therein for Lebesgue regularity issues of such flux $\phi$). It is a Radon measure and one needs to use more involved techniques to treat this case. We will back on this case on forthcoming works.

## 7.5 Linear growth cost $F^*$: total variation flow

A common scenario that we do not address in this paper corresponds to $F$ such that $F^*$ has linear growth function (cf. [32]). This case arises in instances like total variation flow (cf. [7]) and doubly nonlinear total variation flow (cf. [43]), where $p = 1$. In this case, the functional framework outlined in Section 5 is no longer applicable. Consequently, a thorough analysis of this scenario necessitates a careful examination in the appropriate functional spaces.

## 7.6 Numerical computation

Gradient descent algorithms are a powerful class of iterative optimization methods that leverage the gradient of a function to guide the search for its optimal solution. While they are typically
employed for nonlinear optimization problems, as we demonstrate in Section 6, they can also effectively approximate the solutions of PDEs. This algorithmic approach has gained prominence in the context of Wasserstein optimization, particularly through the JKO algorithm, a cornerstone of numerical approximation in this framework. For a comprehensive understanding of these methods and their numerical considerations in this contest, we refer the reader to the works [10, 11] and Chapter 6 of the book [52].

To our knowledge, the recent contribution [27] pioneers the use of steepest descent through minimum flow to numerically tackle evolution dynamics of the form (6.1). The authors highlight the algorithm's adaptability in handling both quadratic and linear growth cost F, incorporating Hele-Shaw type nonlinearities to model congestion phenomena in crowd motion. Moreover, the algorithm's flexibility enables the incorporation of diverse practical scenarios as well as mixed and non-homogeneous boundary conditions, paving the way for more realistic simulations.

So, basically the steepest descent through minimum flow algorithm introduced in this paper presents a novel approach for numerically addressing a broad spectrum of composite PDEs. By utilizing sophisticated numerical optimization techniques known for their effectiveness, even in the presence of low-regularity nonlinearities, this algorithm opens up new methods for PDE analysis and simulation. Specifically, for numerical approximation of the solution of (6.1), Euler-Implicit discretization proves to be a valuable tool. This leads to the numerical solution of the PDE (6.4). By connecting this PDE to the steepest descent formulation of the type $\mathcal{N}$ and exploiting the duality between these problems, we gain access to numerical convex optimization techniques, particularly primal-dual methods such as the alternating direction method of multipliers (ADMM), the augmented Lagrangian method (cf. [31]), and the Chambolle-Pock algorithm (cf. [16]). These methods offer a robust framework for handling the challenges of PDE approximation, particularly when dealing with complex nonlinearities and non-homogeneous boundary conditions.

# A Appendix : technical tools

## A.1 On $W_{\Gamma_D}^{-1,p'}(\Omega)$ characterization

The aim of the following Lemma is to give a practical characterization of the dual bracket $\langle \cdot, \cdot \rangle_{W_{\Gamma_D}^{-1,p'}(\Omega), W_{\Gamma_D}^{1,p'}(\Omega)}$, that we denote by $\langle \cdot, \cdot \rangle_{\Omega}$, in terms of PDE. To this aim, we prove a decomposition of elements of $W_{\Gamma_D}^{-1,p'}(\Omega)$ through a divergence formulation that includes a trace of these elements on $\Gamma_N$. More involved proofs for this type of results that remain valid for a more general class of $\Omega$ can be found in [42] for $p = 2$ and [5] for general $1 < p < \infty$.

**Lemma A.7** Under the assumptions of Section 4, $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$ if and only of one the following equivalent formulation is fulfilled

(A1) there exists $f_0 \in L^{p'}(\Omega)$, $\bar{f} \in L^p(\Omega)^N$, such that

$$ (1.1) \qquad \langle f, \xi \rangle_{\Omega} = \int_{\Omega} f_0 \xi \, dx - \int_{\Omega} \bar{f} \cdot \nabla \xi \, dx, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

(A2) there exists $f_0 \in L^{p'}(\Omega)$, $\tilde{f} \in L^p(\Omega)^N$ and $f_{\Gamma_N} \in W_{\Gamma_N}^{-1,p'}(\Omega)$, such that

$$ (1.2) \qquad \langle f, \xi \rangle_{\Omega} = \int_{\Omega} f_0 \xi \, dx - \int_{\Omega} \tilde{f} \cdot \nabla \xi \, dx + \langle f_{\Gamma_N}, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$
**Proof**: It is clear that for a given $f_0 \in L^{p'}(\Omega)$ and $\bar{f} \in L^{p'}(\Omega)^N$, if $f$ satisfies (1.1) then $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$. Similarly, for a given $f_0 \in L^{p'}(\Omega)$, $\tilde{f} \in L^{p'}(\Omega)^N$ and $f_{\Gamma_N} \in W^{-1/p',p'}(\Gamma_N)$, if $f$ satisfies (1.2) then $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$. To prove that (A1) is a necessary condition, for a given $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$, we consider the the minimization problem

$$ (1.3) \qquad \min_{z \in W_{\Gamma_D}^{1,p}(\Omega)} \left\{ \frac{1}{p} \int_{\Omega} |\nabla z|^p dx - \langle f, z \rangle_{\Omega} \right\} $$

Working as in Lemma 1, we can prove that this problem has a solution $z$ which satisfies

$$ \int_{\Omega} |\nabla z|^{p-2} \nabla z \cdot \nabla \xi = \langle f, \xi \rangle_{\Omega}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Then, we deduce that (A1) is fulfilled by taking simply $f_0 = f$ and $\bar{f} = 0$ whenever $f \in L^{p'}(\Omega)$, otherwise $f_0 = 0$ and $\bar{f} = |\nabla z|^{p-2}\nabla z$.

Now let us prove that (A2) is also a necessary condition. Since $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$, we know that there exists $f_0 \in L^{p'}(\Omega)$ and $\tilde{f} \in L^{p'}(\Omega)^N$ such that $f$ coincides with $f_0 + \nabla \cdot \tilde{f}$ in the dual Sobolev space $(W_0^{1,p}(\Omega))^* ; \text{i.e.}$

$$ \langle f, \xi \rangle_{\Omega} = \int_{\Omega} f_0 \xi \, dx - \int_{\Omega} \tilde{f} \cdot \nabla \xi \, dx, \quad \forall \xi \in W_{0}^{1,p}(\Omega). $$

Keeping in mind that $f$ satisfies (1.3), we have

$$ \int_{\Omega} \phi \cdot \nabla \xi \, dx = \langle f, z \rangle_{\Omega} = \int_{\Omega} f_0 \xi \, dx - \int_{\Omega} \tilde{f} \cdot \nabla \xi \, dx, \quad \forall \xi \in W_{0}^{1,p}(\Omega), $$

so that that $\nabla \cdot (\phi + \tilde{f}) = f_0$ in $D'(\Omega)$. In particular, this implies that $\phi + \tilde{f} \in H^{p'}(\text{div}, \Omega)$, $(\phi + \tilde{f}) \cdot \nu \in W^{-\frac{1}{p'},p'}(\Gamma_N)$ is well defined, and we have

$$ \int_{\Omega} (\phi + \tilde{f}) \cdot \nabla \xi \, dx = \int_{\Omega} f_0 \xi \, dx + \langle (\phi + \tilde{f}) \cdot \nu, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Thus (1.2) is fulfilled by taking $f_{\Gamma_N} := (\phi + \tilde{f}) \cdot \nu$. $\square$

**Remark A.8** 1. Lemma A.7 aims to give a practical characterization of the dual bracket $\langle \cdot, \cdot \rangle_{W_{\Gamma_D}^{-1,p'}(\Omega), W_{\Gamma_D}^{1,p}(\Omega)}$ in terms of PDE. This characterizations emerges a contribution of any element $f \in W_{\Gamma_D}^{-1,p'}(\Omega)$ on the boundary $\Gamma_N$. In particular, Lemma A.7 presents two equivalent manners to handle this contribution.

(a) In (A1), this contribution is linked to the normal trace of $\bar{f}$ in (1.1), if it ever exists in a some "good sense". Indeed, solving $-\nabla \cdot \phi = f \in W_{\Gamma_D}^{-1,p'}(\Omega)$ is equivalent to find $\phi$ such that the PDE

$$ \left\{ \begin{array}{ll} \nabla \cdot \phi = f_0 + \nabla \cdot \bar{f} & \text{in } \Omega \\ (\phi + \bar{f}) \cdot \nu = 0 & \text{on } \Gamma_D \end{array} \right. $$
is fulfilled in a weak sense; i.e.

$$ \int_{\Omega} (\phi + \bar{f}) \cdot \nabla \xi = \int_{\Omega} f_0 \xi, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

In some sense the boundary condition $(\phi + \bar{f}) \cdot \nu$ needs to be understood in $W^{-1/p',p'}(\Gamma_N)$. Yet, despite particular cases, the normal trace of $\bar{f}$ is not $W^{-1/p',p'}(\Gamma_N)$ in general.

(b) The characterization (A2) aims to exhibit explicitly a $W^{-1/p',p'}(\Gamma_N)$ contribution of $f$ on the boundary $\Gamma_N$. Yet, one sees that the formulation still involve the normal boundary trace of $\tilde{f}$ which is not $W^{-1/p',p'}(\Gamma_N)$ in general. See that, solving $-\nabla \cdot \phi = f \in W_{\Gamma_D}^{-1,p'}(\Omega)$ with this formulation is equivalent to find $\phi$ such that the PDE

$$ \left\{ \begin{array}{ll} \nabla \cdot \phi = f_0 + \nabla \cdot \tilde{f} & \text{in } \Omega \\ (\phi + \tilde{f}) \cdot \nu = \pi & \text{on } \Gamma_D \end{array} \right. $$

in a weak sense; i.e.

$$ \int_{\Omega} (\phi + \bar{f}) \cdot \nabla \xi = \int_{\Omega} f_0 \xi + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Again, the boundary condition $(\phi + \bar{f}) \cdot \nu$ needs to be understood in $W^{-1/p',p'}(\Gamma_N)$. In some sense, the formulation (A2) aims just to rewrite $\bar{f}$ by emerging explicitly an $W^{-1/p',p'}(\Gamma_N)$ term. This makes it easier to write $f$ in situations where $\bar{f}$ can be expressed clearly using a $\tilde{f}$, which has a null normal trace on the boundary $\Gamma_N$, and a $W^{-1/p',p'}(\Gamma_N)$ term. Yet, this is not fulfilled in general.

2. As we see in the proof, the term $f_0$ may be considered equals to 0. Nevertheless, we keep this term in the decomposition to make it easier to write $f$ in situations where $f$ is an $L^{p'}$ function.

## A.2 PDE chain rule formula

**Lemma A.7** Let $\rho$, $f \in L^p'(0, T; W_{\Gamma_D}^{-1,p'}(\Omega))$ and $\Psi \in L^p'(Q)$ be such that $\partial_t \rho = \nabla \cdot \Psi + f$ in $Q$ and $\Psi \cdot \nu = \pi$ on $\Sigma_N$ in the sense that

$$ \frac{d}{dt} \int_{\Omega} \rho \xi \, dx + \int_{\Omega} \Psi \cdot \nabla \xi \, dx = \langle f, \xi \rangle_{\Omega} + \langle \pi, \xi \rangle_{\Gamma_N}, \quad \text{in } \mathcal{D}'([0, T)), \quad \forall \xi \in W_{\Gamma_D}^{1,p}(\Omega). $$

Then, for any $\eta \in L^p(0, T; W^{1,p}(\Omega))$ such that $\eta = g$ on $\Gamma_D$ and $\eta \in \partial\beta(., \rho)$ a.e. in $Q$, we have

$$ (1.4) \quad \frac{d}{dt} \int_{\Omega} (\beta(x, \rho) - \tilde{g} \rho) + \int_{\Omega} \Psi \cdot \nabla (\eta - \tilde{g}) = \langle f, \eta - \tilde{g} \rangle_{\Omega} + \langle \pi, \eta \rangle_{\Gamma_N}, \quad \text{in } \mathcal{D}'([0, T)). $$

**Proof:** Let $\psi \in \mathcal{D}([0, T[ \times \mathbb{R}^N)$ be such that $\psi \ge 0$. For any $h > 0$, we consider

$$ \eta_h^n(t) := \frac{1}{h} \int_t^{t+h} \underbrace{(\eta(s) - \tilde{g})}_{=: \eta_g(s)} \psi(s) ds, $$
See that $\eta_h^n$ can be used as test function in (1.4) and therefore

$$ \underbrace{\int_0^T \int_{\Omega} \rho(t) \frac{\eta_g(t+h)\psi(t+h) - \eta_g(t)\psi(t)}{h}}_{=:I_h^n} + \int_0^T \int_{\Omega} \Psi \cdot \nabla \eta_h^n = \int_0^T \langle f, \eta_h^n \rangle_{\Omega} + \int_0^T \langle \pi, \eta \rangle_{\Gamma_N}. $$

We see that

$$ \begin{aligned} I_h^n &= \int_0^T \int_{\Omega} \rho(t) \frac{\eta(t+h)\psi(t+h) - \eta(t)\psi(t)}{h} - \int_0^T \int_{\Omega} g \rho(t) \frac{\psi(t+h) - \psi(t)}{h} \\ &= \int_0^T \int_{\Omega} \frac{\rho(t-h) - \rho(t)}{h} \eta_g(t)\psi(t) - \int_0^T \int_{\Omega} g \rho(t) \frac{\psi(t+h) - \psi(t)}{h} \\ &\le \int_0^T \int_{\Omega} \frac{\beta(\rho(t-h)) - \beta(\rho(t))}{h} \psi(t) - \int_0^T \int_{\Omega} g \rho(t) \frac{\psi(t+h) - \psi(t)}{h} \\ &\le \int_0^T \int_{\Omega} \frac{\psi(t+h) - \psi(t)}{h} \beta(\rho(t)) - \int_0^T \int_{\Omega} g \rho(t) \frac{\psi(t+h) - \psi(t)}{h}. \end{aligned} $$

So, taking limits as $\tau \to 0^+$ we get

$$ \int_0^T \int_{\Omega} \partial_t^n \psi (\beta(x, \rho(t)) - \tilde{g}) + \int_0^T \int_{\Omega} \Psi \cdot \nabla (\eta - \tilde{g}) \psi \ge \int_0^T \psi \langle f, \eta - \tilde{g} \rangle_{\Omega} + \int_0^T \psi \langle \pi, \eta \rangle_{\Gamma_N} $$

Taking now $\tilde{\eta}_{\tau}(t) = \frac{1}{\tau} \int_{t}^{t+\tau} \eta(s-\tau)\psi(s)ds$, and arguing as above we get the converse inequality. That is

$$ \int_0^T \int_{\Omega} \partial_t^n \psi (\beta(x, \rho(t)) - \tilde{g}) + \int_0^T \int_{\Omega} \Psi \cdot \nabla (\eta - \tilde{g}) \psi = \int_0^T \psi \langle f, \eta - \tilde{g} \rangle_{\Omega} + \int_0^T \langle \pi, \eta \rangle_{\Gamma_N}. $$

Thus (1.4). □

## A.3 Week Aubin's result

**Theorem A.5 (cf. [6] and [48])** Let $(u_n)_n$ and $(v_n)_n$ be two bounded sequences in $L^p(0, T; W^{1,p}(\Omega))$ and $L^p(Q)$, respectively. Assume moreover that $v_n$ satisfies the following weak estimates

$$ \int_0^T \int_{\Omega} v_n \partial_t^n \varphi \, dxdt \le C \|\nabla \varphi\|_{L^\infty(Q)}, \quad \text{for any } \varphi \in \mathcal{D}(Q). $$

If $u \in L^p(0, T; W^{1,p}(\Omega))$ and $v \in L^p(Q)$, are such that, up to extraction of a subsequence,

$$ u_n \to u \text{ in } L^p(0, T; W^{1,p}(\Omega)) - \text{weak, and } v_n \to v \text{ in } L^p(Q) - \text{weak} $$

then

$$ \int_0^T \int_{\Omega} u_n v_n \varphi \, dxdt \to \int_0^T \int_{\Omega} uv \varphi \, dxdt \quad \text{for any } \varphi \in \mathcal{D}(Q). $$

If moreover, $v_n \in \alpha(., u_n)$ a.e. in $Q$, where $\alpha(z, .)_{z \in Q}$ is a sequence of measurable maximal monotone graph with $0 \in \alpha(., 0)$, a.e. in $Q$, then the following assertions are fulfilled
1. $v \in \alpha(., u)$ a.e. in $Q$,

2. a.e. $(t,x) \in Q$, either $u_n(t,x) \to u(t,x)$ or $v_n(t,x) \to v(t,x)$,

3. if $\alpha$ (resp. $\alpha^{-1}$) is single valued then $v_n \to v$ (resp. $u_n \to u$) a.e. in $Q$.

For the proof of this theorem one can see the papers [6] and [48].

## Acknowledgments

We would like to thank Chérif Amrouche for many interesting discussions about Sobolev dual spaces.

## References

[1] M. AGUEH Existence of solutions to degenerate parabolic equations via the Monge-Kantorovich theory. *Adv. Differential Equations*, 10(3): 309-360 (2005).

[2] H. W. ALT and H. W. LUCKHAUS. Quasilinear elliptic-parabolic differential equations. *Math. Z.*, 183(1983), pp. 311–341.

[3] L. AMBROSIO Minimizing movements. *Rend. Accad. Naz. Sci. XL Mem. Mat. Appl.*, (5), 19(1995), pp. 191-246.

[4] L. AMBROSIO N. GIGLI and G. SAVARÉ. Gradient flows in metric spaces and in the space of probability measures. 2nd edn, Lecture notes in Mathematics ETH Zurich (Birkhauser, Basel, , 2008

[5] Ch. AMROUCHE Preprint, 2024.

[6] B. ANDREIANOV and C. CANCES and A. MOUSSA. A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs. J. Func. Anal., 2017, Vo 273(12), 3633-3670.

[7] F. ANDREU and V. CASELLES and J. MAZÓN. Parabolic Quasilinear Equations Minimizing Linear Growth Functionals. Progress in Mathematics, Birkhauser, Basel, vol. 223. (2004)

[8] V. BARBU. Nonlinear differential equations of monotone types in Banach spaces. Springer Monographs in Mathematics, Springer, New York, 2010.

[9] M. BECKMANN. A continuous model of transportation. Econometrica, 20 (1952), 643-660.

[10] J.D. BENAMOU, G. CARLIER and M. LABORD. An augmented Lagrangian approach to Wasserstein gradient flows and applications. ESAIM: Proceedings and Surveys, Volume 54, Issue 1: 1 - 17 (2016).

[11] J.D. BENAMOU, G. CARLIER, M. M. CUTURI L. NENNA and G. PEYRÉ Iterative Bregman projections for regularized transportation problems. SIAM Journal on Scientific Computing, 37(2), (2015), A1111-A1138.
[12] G. BOUCHITTÉ and G. BUTTAZZO. Characterization of optimal shapes and masses through Monge-Kantorovich equation. J. Eur. Math. Soc., 3 (2), 139-168, 2001.

[13] Y. BRENIER, Extended Monge-Kantorovich theory, in Optimal transportation and applications, (Martina Franca, 2001), vol. 1813 of Lecture Notes in Math., Springer, Berlin, 2003, pp. 91–121

[14] H. BRÉZIS, Opérateurs maximaux monotones et semigroups de contractions dans les espaces de Hilbert. (French), North-Holland Publishing Co., Amsterdam-London; American Elsevier Publishing Co., Inc., New York, 1973.

[15] J. CARILLO, Entropy solutions for nonlinear degenerate problems. Arch. Ration. Mech. Anal., 147(4):269-361, (1999).

[16] A. CHAMBOLLE and T. POCK, A first-order primal-dual algorithm for convex problems with applications to imaging, J. Math. Imaging Vis. 40(1):120-145, 2011.

[17] G. Q. CHEN, Some recent methods for partial differential equations of divergence form. Bulletin of the Brazilian Mathematical Society, 34(1):107-144, 2003.

[18] G. Q. CHEN and H. FRID, Divergence-measure fields and hyperbolic conservation laws. Archive for Rational Mechanics and Analysis, 147(2):89-118, 1999.

[19] G. Q. CHEN and H. FRID, Extended divergence-measure fields and the Euler equations for gas dynamics. Communications in mathematical physics, 236(2):251-280, 2003.

[20] M. G. CRANDALL and T.M. LIGGETT, Generation of semigroups of nonlinear transformations on general Banach spaces. Am. J. Math. 93, 265-298 (1971)

[21] A. DAMLAMIAN, Some results on the multiphase Stefan problem Comm. P.D.E, 2(1977),1017-1044.

[22] A. DAMLAMIAN and N. KENMOCHI, Evolution equations generated by subdifferentials in the dual space of $H^1(\Omega)$ Discrete Contin. Dynam. Systems, 5 (1999) 269-278.

[23] N. DAVID and M. SCHMIDTCHEN, On the incompressible limit for a tumour growth model incorporating convective effects Com. PDE. October 2023, https://doi.org/10.1002/cpa.22178.

[24] S. DI MARINO and A.R. MÉSZÁROS. Uniqueness issues for evolution equations with density constraints. Mathematical Mod. Meth. App. Sci., 26 (09), 1761-1783.

[25] S. DWEIK. Summability of the transport density in the import-export transport problem with Riemannian cost. Preprint 2024.

[26] I. EKELAND and R. TÉMAM. Convex Analysis and Variational Problems. Classics in Mathematics, Society for Industrial and Applied Mathematics, 1999.

[27] H. ENNAJI N. IGBIDA and G. JRADI. Prediction-Correction Pedestrian Flow by Means of Minimum Flow Problem. Math. Models Methods Appl. Sci., to appear in 2024.

[28] H. ENNAJI N. IGBIDA and V. Th. NGUYEN. Beckmann-Type Problem for degenerate Hamilton-Jacobi equations. Quart. Appl. Math, 80 (2022), 201-220.
[29] H. ENNAJI N. IGBIDA and V. Th. NGUYEN. Augmented Lagrangian method for Hamilton-Jacobi equations. Calculus of Variations and PDE, Vol. 60(23), 2021.

[30] A. FIGALI and N. GIGLI. A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions. Journal de Mathématiques Pures et Appliquées Vo. 94(2), 2010, 107-130.

[31] R. GLOWINSKI and P. LE TALLEC. Augmented Lagrangian and operator-splitting methods in nonlinear mechanics. volume 9. SIAM, 1989.

[32] W. GORNY and J. MAZON. A duality-based approach to gradient flows of linear growth functionals. ArXiv preprint, https://arxiv.org/pdf/2212.08725.pdf

[33] D. KWON and A.R. MÉSZÁROS. Degenerate nonlinear parabolic equations with discontinuous diffusion coefficients. J. Lond. Math. Soc., 104(2), (2021), pp. 688-746.

[34] N. IGBIDA. $L^1$-Theory for Hele-Shaw flow with linear drift. M3AS, Math. Mod. Meth. Appl. Sci., Vol. 33, No. 07, 1545-1576 (2023).

[35] N. IGBIDA. $L^1$-Theory for Incompressible limit of Porous Medium Equation with linear Drift. ArXiv Preprint, https://arxiv.org/pdf/2107.02606.pdf

[36] N. IGBIDA and V. Th NGUYEN. Optimal Partial Mass Transportation and Obstacle Monge-Kantorovich Equation. J. Differential Equations, 264(10), 6380-6417, 2018.

[37] N. IGBIDA and V. Th NGUYEN. Augmented Lagrangian method for optimal partial transportation. IMA Journal of Numerical Analysis, 38(1):156-183, 2018.

[38] N. IGBIDA and J. M. URBANO. A granular approach to crowd motion. Preprint 2024.

[39] R. JORDAN, D. KINDERLEHRER, F. OTTO, The variational formulation of the Fokker-Planck equation SIAM J. Math. Anal. 29 (1998) 1-17.

[40] N. KENMOCHI, Neumann problems for a class of nonlinear degenerate parabolic equations, Differential Integral Equations, 3 (1990) 253-273.

[41] D. KWON and A. R. MÉSZÁROS. Degenerate nonlinear parabolic equations with discontinuous diffusion coefficients. J. Lond. Math. Soc., 104(2), (2021), pp. 688-746.

[42] J. L. LIONS and E. MAGENES. Problèmes aux Limites Non-Homogènes et Applications Vol. 1, Dunod Paris, 1968

[43] J. MAZON, A. MOLINO and J. TOLEDO. Doubly nonlinear equations for the 1-Laplacian J. Evol. Equ., (2023) 23:67, 26 pp.

[44] B. MAURY, A. ROUDNEFF-CHUPIN, F. SANTAMBROGIO, A macroscopic crowd motion model of gradient flow type. Math. Models and Meth. in Appl. Sci., 20 (2010), No. 10, 1787-1821.

[45] B. MAURY, A. ROUDNEFF-CHUPIN, F. SANTAMBROGIO, Congestion-driven dendritic growth. Discrete Con-tin. Dyn. Syst., 34 (2014), no. 4, 1575-1604.
[46] B. MAURY, A. ROUDNEFF-CHUPIN, F. SANTAMBROGIO and J. VENEL. Handling congestion in crowd motion modeling, Netw. Heterog. Media, 6 (2011), No. 3, 485-519.

[47] A. R. MÉSZÁROS and F. SANTAMBROGIO. Advection-diffusion equations with density constraints Analysis and PDE, 9(3):615-644, (2016).

[48] A. MOUSSA. Some variants of the classical Aubin-Lions Lemma J. of Evolution Equations, Vo. 16, 65-93, (2016).

[49] L. A. PELETIER. The porous media equations Applications of Nonlinear Analysis in Physical Sciences (Edited by H.Amann, N.Bazley, and Kirchgassner), pp. 229 - 241(1981).

[50] F. OTTO. $L^1$ contraction and uniqueness for quasilinear elliptic-parabolic equations. *J. Diff. Eq.*, 131(1996), pp. 20-38.

[51] F. OTTO. The geometry of dissipative evolution equations: the porous medium equation. Comm. Partial Differential Equations, 26 (2001) 101-174.

[52] F. SANTAMBROGIO. Optimal Transport for Applied Mathematicians. Progress in Nonlinear Differential Equations and Their Applications, no 87, Birkhauser Basel (2015).

[53] J. L. VÁZQUEZ. The porous medium equation : mathematical theory. Oxford University Press, 2007.

[54] A. VISINTIN. Models of Phase Transitions. Progress in nonlinear Differential Equations and their Applications. Birkhauser, 1996.