# Numerical approximation of a class of constrained Hamilton-Jacobi equations

Benoît Gaudeul, Hélène Hivert

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

## To cite this version:

Benoît Gaudeul, Hélène Hivert. Numerical approximation of a class of constrained Hamilton-Jacobi equations. 2024. hal-04510185

HAL Id: hal-04510185

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

Preprint submitted on 18 Mar 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
Numerical approximation of a class of constrained Hamilton-Jacobi equations

Benoît Gaudeul\*, Hélène Hivert†

Abstract

In this paper, we introduce a framework for the discretization of a class of constrained Hamilton-Jacobi equations, a system coupling a Hamilton-Jacobi equation with a Lagrange multiplier determined by the constraint. The equation is non-local, and the constraint has bounded variations. We show that, under a set of general hypothesis, the approximation obtained with a finite-differences monotonic scheme, converges towards the viscosity solution of the constrained Hamilton-Jacobi equation.

Constrained Hamilton-Jacobi equations often arise as the long time and small mutation asymptotics of population models in quantitative genetics. As an example, we detail the construction of a scheme for the limit of an integral Lotka-Volterra equation. We also construct and analyze an Asymptotic-Preserving (AP) scheme for the model outside of the asymptotics. We prove that it is stable along the transition towards the asymptotics.

The theoretical analysis of the schemes is illustrated and discussed with numerical simulations. The AP scheme is also used to conjecture the asymptotic behavior of the integral Lotka-Volterra equation, when the environment varies in time.

Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>Acknowledgements</td><td>5</td></tr><tr><td>2 Construction of the scheme and main results</td><td>5</td></tr><tr><td>3 An example: Lotka-Volterra integral equations</td><td>9</td></tr><tr><td>3.1 Examples of schemes</td><td>9</td></tr><tr><td>3.2 Asymptotic-preserving scheme</td><td>10</td></tr><tr><td>4 Convergence of scheme (S)</td><td>14</td></tr><tr><td>4.1 Proof of Prop. 2.3</td><td>14</td></tr><tr><td>4.2 Proof of Prop. 2.4</td><td>19</td></tr><tr><td>4.3 Relaxation of the CFL condition</td><td>26</td></tr><tr><td>5 Numerical tests and discussion</td><td>27</td></tr><tr><td>5.1 Properties of (S)</td><td>28</td></tr><tr><td>5.2 Asymptotic-preserving property of scheme (S<sub>ε</sub>)</td><td>29</td></tr><tr><td>5.3 Scheme (S<sub>ε</sub>) when the population vanishes</td><td>31</td></tr><tr><td>5.4 Generalization</td><td>33</td></tr><tr><td>A Asymptotic-preserving property of (S<sub>ε</sub>)</td><td>35</td></tr><tr><td>A.1 Proof of Prop. 3.3</td><td>35</td></tr><tr><td>A.2 Proof of Prop. 3.4</td><td>37</td></tr></table>

<footnote>\*Université Paris-Saclay, CNRS, Laboratoire de mathématiques d'Orsay, 91405, Orsay, France.
benoit.gaudeul@universite-paris-saclay.fr</footnote>

<footnote>†Univ Rennes, Inria, Géosciences Rennes - UMR 6118, F-35000 Rennes, France. helene.hivert@inria.fr</footnote>
# 1 Introduction

We are interested in the design and analysis of numerical schemes for a class of constrained Hamilton-Jacobi equations,

$$ \left\{ \begin{array}{ll} \partial_t u(t,x) + b(x, I(t)) \mathcal{H}(\nabla_x u(t,x)) + \mathcal{R}(t,x, I(t)) = 0, & x \in \mathbb{R}^d, t > 0 \\ \min_{x \in \mathbb{R}^d} u(t,x) = 0, & t \ge 0, \end{array} \right. \qquad (HJ) $$

supplemented with an initial condition $u^{in}$, and where the Hamiltonian $p \mapsto \mathcal{H}(p)$ is strictly convex. Equations such as (HJ) often arise as long time and small mutations asymptotics of models of population structured by a phenotypic trait [15, 36, 9, 13, 35, 30, 29]. This is a particular case of models of the theory of adaptive evolution [33, 21, 20, 32, 14]. They describe the size $n_\varepsilon(t, x)$ of a population, where $t$ denotes the time, $x$ the phenotypic trait, and $\varepsilon \in (0, 1]$ a scaling parameter. In what follows, we will consider as an example the following Lotka-Volterra integral equation

$$ \varepsilon \partial_t n_\varepsilon(t, x) - \frac{1}{\varepsilon^d} \int_{\mathbb{R}^d} \mathcal{K}\left(\frac{y-x}{\varepsilon}\right) b(y, J_\varepsilon(t)) n_\varepsilon(t, y) dy - \mathcal{R}(t, x, J_\varepsilon(t)) n_\varepsilon(t, x) = 0, \quad x \in \mathbb{R}^d, t > 0, \quad (1.1) $$

where $J_\varepsilon$ denotes a weighted size of the population. The weight $\psi$ being given, it is defined by

$$ J_{\varepsilon}(t) = \int_{\mathbb{R}^d} \psi(x) n_{\varepsilon}(t, x) dx, $$

for all $t \ge 0$. Equation (1.1) is supplemented with an initial condition $n_\varepsilon^{in} \in L^1(\mathbb{R}^d)$. The parameters in (1.1) and (HJ) are chosen according to the biological context. Indeed, with the above notations, the evolution of the population is driven by births, through the birth rate $b(x, J_\varepsilon)$, and deaths. The net growth rate is denoted here $\mathcal{R}(t, x, J_\varepsilon)$. Note that $b$ and $\mathcal{R}$ depend on the phenotypic trait $x$, meaning that some individuals may be advantaged because they are better adapted. They also depend on the population burden on the environment, through the size of the population $J_\varepsilon$. The phenotypic trait of the parent is transmitted to its offspring, possibly with mutations. In [5, 4] and in (1.1), it is modelled by an integral, where $\mathcal{K}$ denotes a probability kernel. It could also have been represented by a Laplacian [37, 4, 31].

Considering the limit $\varepsilon \to 0$ in models such as (1.1), stands for the study of the population, in an asymptotic regime of long time and small mutations. It is usually referred to, as the separation of ecological and evolutionary time scales. If the population does not extinct or grow uncontrolled, when $\varepsilon \to 0$, the population $n_\varepsilon(t, x)$ is expected to concentrate around a set of dominant traits, see for instance [37, 9]. The dynamics of this concentration is usually studied using the Hopf-Cole transform, a logarithmic transform of the unknown, as in [5, 9, 15, 31]. Namely, $v_\varepsilon = -\varepsilon \ln(n_\varepsilon)$ is introduced. With the example (1.1), it satisfies

$$ \left\{ \begin{array}{ll} \partial_t v_\varepsilon(t, x) + \int_{\mathbb{R}^d} \mathcal{K}(z) b(x + \varepsilon z, J_\varepsilon(t)) e^{-(v_\varepsilon(t,x+\varepsilon z)-v_\varepsilon(t,x))/\varepsilon} dz + \mathcal{R}(t, x, J_\varepsilon(t)) = 0 & (P_\varepsilon) \\ J_\varepsilon(t) = \int_{\mathbb{R}^d} \psi(x) e^{-v_\varepsilon(t,x)/\varepsilon} dx, & \end{array} \right. $$

with initial condition $v_\varepsilon^{in} = -\varepsilon \ln(n_\varepsilon^{in})$. Under suitable assumptions, see [5, 4], $v_\varepsilon$ is expected to converge, when $\varepsilon \to 0$, to the solution $u$ of a constrained Hamilton-Jacobi equation such as (HJ), with

$$ \forall p \in \mathbb{R}^d, \mathcal{H}(p) = \int_{\mathbb{R}^d} \mathcal{K}(z) e^{-z \cdot p} dz. \qquad (1.2) $$

The size of the population $J_\varepsilon$, also converges, its limit is denoted $I$ in (HJ). In the asymptotic regime, it is not defined as the population size anymore, but it behaves as a Lagrange multiplier regarding the constraint $\min u = 0$.

The main difficulty in the analysis of (HJ) comes from the regularity of its solution. Indeed, $u$ is expected to enjoy no more than Lipschitz regularity, as it can be expected for viscosity solutions of Hamilton-Jacobi equation [3, 11, 17], but $I$ may have jumps discontinuities [37, 36, 4]. Because of these discontinuities, the appropriate functional space for the well-posedness of (HJ) is $W_{\text{loc}}^{1,\infty}(\mathbb{R}^d) \times BV_{\text{loc}}$. The uniqueness of the pair $(u, I)$ solution of (HJ) has been addressed in some particular cases [37, 34, 26], and then in a more general setting close to the problem under study [8]. Keeping in mind the biological models it comes from, the hypothesis set on the parameters in (HJ) are done following [4], and assumptions are added for (HJ) to be in the framework of [8]. In what follows, we will then suppose that $\mathcal{H}$, $b$, $\mathcal{R}$, and $u^{in}$ satisfy the following assumptions:
• $\mathcal{H}: \mathbb{R}^d \to \mathbb{R}^d$, is of class $C^2$, and such that

$$ \mathcal{H} \text{ is strictly convex, and } \lim_{|p| \to +\infty} \frac{\mathcal{H}(p)}{|p|} = +\infty, \qquad (A1) $$

and moreover

$$ \mathcal{H}(0) = 0, \text{ and } \forall p \in \mathbb{R}^d, \mathcal{H}(p) \ge 0. \qquad (A2) $$

Remark 1.1. We assume here that $\mathcal{H}(0) = 0$, but it is in general not true in the models, see [5, 4, 35] for instance. However, this assumption is done for a sake of simplicity and can be relaxed. Indeed, remarking that

$$ \partial_t u + b\mathcal{H}(\nabla_x u) + \mathcal{R} = \partial_t u + b(\mathcal{H}(\nabla_x u) - \mathcal{H}(0)) + \mathcal{R} + b\mathcal{H}(0), $$

gives the key to adapt the results to more general convex Hamiltonians which are minimal at $p = 0$.

• $b: \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$, is a smooth positive and bounded function. We suppose that there exists $b_M > 1$ such that

$$ \forall x \in \mathbb{R}^d, \forall I \in \mathbb{R}, 0 < b_M^{-1} \le b(x, I) \le b_M, \text{ and } \forall (i, j) \in [1, d]^2, |\partial_{x_i x_j}^2 b(x, I)| \le b_M, \qquad (A3) $$

and such that it is also a Lipschitz constant of $b$ in both variables

$$ \forall x, y \in \mathbb{R}^d, \forall I, J \in \mathbb{R}, |b(x, I) - b(y, J)| \le b_M (|x - y| + |I - J|). \qquad (A4) $$

Moreover, we suppose that

$$ \forall x \in \mathbb{R}^d, I \mapsto b(x, I) \text{ is non-increasing.} \qquad (A5) $$

This hypothesis is relevant when considering the biological meaning of $b$ and $I$. Indeed, it implies that the larger the population is, the smaller its birth rate is, as the environment may not offer enough resources for the population.

• $\mathcal{R}: \mathbb{R}_+ \times \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$, is smooth, and decreasing in its last variable. More precisely, there exists a constant $K_1$ such that

$$ \forall t \ge 0, \forall x \in \mathbb{R}^d, \forall I \in \mathbb{R}, -K_1 \le \partial_I \mathcal{R}(t, x, I) \le -K_1^{-1} < 0, \qquad (A6) $$

and $I_M \ge I_m > 0$, such that

$$ \frac{I_M}{K_1} - b_M \ge 1, \quad \min_{t \ge 0, x \in \mathbb{R}^d} \mathcal{R}(t, x, I_M) \ge 0, \quad \max_{t \ge 0, x \in \mathbb{R}^d} \mathcal{R}(t, x, I_M) \le 0. \qquad (A7) $$

We emphasize on the fact, that the first condition is purely technical, and will appear naturally in what follows. Note that since $\mathcal{R}$ is decreasing in its last variable, $I_M$ can be chosen according to the third condition in (A7), and then be increased enough so that the first one is satisfied. Let us also define $K_2$ such that

$$ \forall t \ge 0, \forall I \in [\min(I_M/2, I_M - 1), 2I_M], \quad \|\mathcal{R}(t, \cdot, I)\|_{W^{2,\infty}(\mathbb{R}^d)} \le K_2, \qquad (A8) $$

and suppose in addition that $(t, I) \to \|\mathcal{R}(t, \cdot, I)\|_{L^\infty(\mathbb{R}^d)}$ is bounded on any bounded set of $\mathbb{R}_+ \times \mathbb{R}$. Finally, the dependence in time of $\mathcal{R}$ is supposed to be Lip$_{\mathcal{R}}$-Lipschitz,

$$ \forall x \in \mathbb{R}^d, \forall I \in \mathbb{R}, \forall t, s \in \mathbb{R}_+, |\mathcal{R}(t, x, I) - \mathcal{R}(s, x, I)| \le \text{Lip}_{\mathcal{R}}|t - s|, \qquad (A9) $$

and to simplify notations, we introduce

$$ K = \max\{K_1, K_2, \text{Lip}_{\mathcal{R}}\}. \qquad (1.3) $$

Remark 1.2. One can notice that $I \mapsto b(x, I)\mathcal{H}(p) + \mathcal{R}(t, x, I)$ is decreasing for all $t \in \mathbb{R}_+, x \in \mathbb{R}^d$ and $p \in \mathbb{R}^d$. This is indeed a consequence of (A2)-(A5)-(A6). It is an important condition for $I$ to be uniquely determined as the Lagrange multiplier of the viscosity solution of (HJ), see [8].

• $u^{\text{in}}$ is Lip$_{\text{in}}$-Lipschitz,

$$ \forall x, y \in \mathbb{R}^d, |u^{\text{in}}(x) - u^{\text{in}}(y)| \le \text{Lip}_{\text{in}}|x - y|, \qquad (A10) $$

is coercive

$$ \exists \bar{a} > \underline{a} > 0, \exists \underline{b}, \bar{b} \in \mathbb{R}, \forall x \in \mathbb{R}^d, \underline{a}|x| + \underline{b} \le u^{\text{in}}(x) \le \bar{a}|x| + \bar{b}, \qquad (A11) $$

and its minimum is equal to 0

$$ u^{\text{in}} \ge 0, \text{ and } \exists x_0 \in \mathbb{R}^d, u^{\text{in}}(x_0) = 0, \qquad (A12) $$

even though this minimum point is not necessarily unique. Note that $|\cdot|$ denotes here the Euclidean norm on $\mathbb{R}^d$.
Considering (HJ) with $\mathcal{R}(t, x, I) = \mathcal{R}(x, I)$, the following result holds

**Theorem 1.1** ([8]). Suppose that assumptions (A1) to (A12) are satisfied, and that $\mathcal{R}$ does not depend on $t$.

(i). Let $(u_1, I_1)$ and $(u_2, I_2)$ two solutions of (HJ) in $W_{\text{loc}}^1 \times BV_{\text{loc}}$ with the same initial data $u^{\text{in}}$. Then $u_1 = u_2$ and $I_1 = I_2$ almost everywhere.

(ii). Let $I \in BV(0, T)$ be given. Then the variational solution $u$ of

$$ \partial_t u(t, x) + b(x, I(t))\mathcal{H}(\nabla_x u(t, x)) + \mathcal{R}(x, I(t)) = 0, \quad t > 0, x \in \mathbb{R}^d, \qquad (1.4) $$

with initial data $u^{\text{in}}$, is the unique locally Lipschitz viscosity solution of (1.4) over $(0, T] \times \mathbb{R}^d$. Moreover, $u$ is independent of the choice of a representative of $I$ in $BV$. Namely, if (1.4) is considered with two source terms $I_1$ and $I_2$ in $BV(0, T)$ such that $I_1 = I_2$ almost everywhere, then $u_1 = u_2$.

**Remark 1.3.** The uniqueness result for the solution of (HJ) is a priori only true when $\mathcal{R}$ is independent of $t$, and when the previous assumptions are satisfied. However, we will assume that it is still valid for more general $\mathcal{R}$, with Lipschitz regularity with respect to $t$. From a modelization point of view, this enables the environment to change in time. Proving that Theorem 1.1 still holds when $\mathcal{R}$ depends on $t$, with Lipschitz regularity would deserve a dedicated study. As Lipschitz regularity is standard for time-dependency in Hamilton-Jacobi equations [3], we believe this conjecture likely to be true.

The analysis of a finite-differences scheme for a simplified (HJ), with $b \equiv 1$, $\mathcal{H}(p) = |p|^2$, and $\mathcal{R}$ independent of $t$, has been carried out in [7]. The goal of this paper is to generalize it, and to propose a framework for the numerical analysis of constrained Hamilton-Jacobi equations, such as (HJ). When the parameters are regular enough, and when dealing with bounded solutions, the approximation of viscosity solutions of Hamilton-Jacobi equations can be handled with finite-differences monotonic schemes, see [12, 38]. Here, the two main difficulties of the problem are hence the coercivity of $u^{\text{in}}$ (A11), and the lack of regularity of $I$. As in [7], we use Theorem 1.1 to overcome separately these concerns, since it enables decoupling the Hamilton-Jacobi equation from its constraint. Indeed, at the continuous level, one can consider $u$, the viscosity solution of (HJ) with $I$ given as a source term, show separately that $\min u(t, \cdot) = 0$ for all $t$, and conclude that $(u, I)$ solves (HJ). Following these considerations, we approximate the classical Hamilton-Jacobi part of the problem with a monotonic scheme which enjoys a discrete maximum principle [12, 38]. The approximation of $I$ is treated with a nonlinear problem to solve, so that the constraint is satisfied. In this step, the monotony of $I \mapsto b(x, I)\mathcal{H}(p) + \mathcal{R}(t, x, I)$ yields stability, and counterparts the lack of regularity of $I$.

However, the generalization of [7] to more general Hamiltonians brings additional difficulties, in particular when proving that the scheme preserves the semi-concavity of the solution of (HJ) (see [17]). This property is crucial for the stability of the scheme, and for the convergence of the approximation of $I$. It was easily satisfied in [7], but it may require much more precise estimates here. We hence distinguish two classes of schemes, depending on their behavior when investigating the discrete semi-concavity of the solutions. The scheme proposed in [7] belongs to the class designated as the flat setting. In this case, the discrete Hamiltonian is positive, as is $\mathcal{H}$, but it is not necessarily convex. As they do not represent a particular difficulty, the semi-concavity estimates for this class of schemes can be established in a general framework, so that the convergence of the schemes belonging to this class is true in any dimension, and even when adding a smooth dependency in $x$ to $\mathcal{H}$. On the contrary, discrete Hamiltonians are convex but not necessarily positive in the convex setting. This framework is natural, as the convexity of $\mathcal{H}$ is a crucial hypothesis for the well-posedness of the continuous problem (HJ), but the discrete semi-concavity estimates are more intricate. In this case, we establish the convergence of the scheme only in dimension 1, and we postpone the investigation of the generalization to a future work.

We then illustrate the construction of schemes for (HJ), by detailing the example where the Hamiltonian $\mathcal{H}$ is defined by (1.2). For this particular case, we propose five approximations of (HJ), and we discuss their properties. We also come back to the biological model, and the $\varepsilon$-dependent problem $(P_\varepsilon)$ that converges to (HJ) in the limit $\varepsilon \to 0$, and we propose and analyze an asymptotic-preserving scheme adapted to this asymptotics. Such schemes enjoy stability properties when $\varepsilon \to 0$, meaning that they are accurate in all the regimes, with no constraint on the discretization parameters. In particular, we show with numerical tests that they catch the asymptotics of the model, even when the population goes to an extinction, so that the asymptotic regime is no more described by (HJ).

The paper is organized as follows: the scheme for (HJ) is constructed and discussed in Section 2. In addition, examples of schemes for the integral-defined Hamiltonian (1.2) are presented in Section 3, where an asymptotic-preserving scheme for $(P_\varepsilon)$ is also constructed and analyzed. The convergence of the scheme for (HJ) is then proved in Section 4. Various properties of the schemes are illustrated and discussed via numerical tests in Section 5. Finally, technical points related to the asymptotic-preserving scheme of Section 3 are proved in Appendix A.
# Acknowledgements

The authors warmly thank Vincent Calvez for the many very interesting discussions they had about this problem. They also thank Pierre Navaro, for his expertise with Julia. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC consolidator grant WACONDY no 865711)

## 2 Construction of the scheme and main results

In this section, we present a general framework for the construction of schemes for (HJ) in dimension $d = 1$, which converge towards the viscosity solution of (HJ). This is a generalization of the scheme proposed in [7] for the special case $\mathcal{H}(p) = p^2$, and it is designed according to the following principles:

* The time-dependent Lagrange multiplier $I$ in (HJ) is implicated in the schemes, to make sure that the constraint $\min u = 0$ in (HJ) is preserved at the discrete level.

* Following [12, 38], $\mathcal{H}$ is discretized with a monotonous explicit scheme. However, in comparison to [12, 38], we relax the consistency hypothesis of the discrete Hamiltonian. Indeed, this assumption means that $\mathcal{H}(p)$ can be computed without approximation for any $p \in \mathbb{R}$, and this is for instance not relevant when $\mathcal{H}$ is defined as an integral, as in (1.2). To this end, we introduce a discretization parameter $\eta$, accounting for the quality of the approximation of $\mathcal{H}: \mathbb{R} \to \mathbb{R}$, and we denote $H_{\eta}: \mathbb{R}^2 \to \mathbb{R}$ the approximated Hamiltonian.

* It may happen, that the net growth rate $\mathcal{R}$ is not computed without approximation, especially when $\mathcal{H}$ is modified as in Remark 1.1. As previously, an approximation of $\mathcal{R}$ is then introduced. Still denoting $\eta$ the discretization parameter, we define $R_{\eta}: [0, T] \times \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$ the approximated net growth rate.

According to these considerations, the properties of $H_{\eta}$ and $R_{\eta}$ are stated as follows. First of all, as in [12, 38, 7], an appropriate discretization of $\mathcal{H}$ should be increasing in its first variable, and decreasing in the second one. However, this property can be slightly relaxed. Indeed, following the ideas of the Lax-Friedrichs scheme, one can add a viscosity to $H_{\eta}$ to force the monotony in a given bounded set of $\mathbb{R}^2$. To allow for these discretizations in the framework of our study, we let $\mathcal{L} > 0$ and assume that

$$ \forall (p, q) \in [-\mathcal{L}, \mathcal{L}], \forall e > 0, H_{\eta}(p, q + e) \le H_{\eta}(p, q) \le H_{\eta}(p + e, q). \quad (A13) $$

We also suppose that no discretization error is committed in (A2), meaning that

$$ H_{\eta}(0, 0) = 0. \quad (A14) $$

Note that this assumption is purely technical, but simplifies several expressions in what follows. However, one could also avoid this restriction by considering $H_{\eta} - H_{\eta}(0, 0)$. For the same reason, we also assume that the discrete Hamiltonian $H_{\eta}$ is non-negative on the diagonal

$$ \forall |p| \le \mathcal{L}, \forall \eta \in (0, 1], H_{\eta}(p, p) \ge 0. \quad (A15) $$

The quality of the approximation of $\mathcal{H}$ by $H_{\eta}$ is quantified as follows

$$ \exists C_{\mathcal{L}} \ge 0, \forall |p| \le \mathcal{L}, \forall \eta \in (0, 1], |H_{\eta}(p, p) - \mathcal{H}(p)| \le C_{\mathcal{L}} \eta, \quad (A16) $$

and the approximation of $\mathcal{R}$ by $R_{\eta}$ is also quantified with $\eta$

$$ \forall t \ge 0, \forall x \in \mathbb{R}, \forall I \ge 0, \forall \eta \in (0, 1], |\mathcal{R}(t, x, I) - R_{\eta}(t, x, I)| \le K\eta, \quad (A17) $$

where the constant $K$ is used again, as in (1.3), to simplify notations. In addition, we suppose that $R_{\eta}$ satisfies (A6)-(A7)-(A8) and (A9), as does $\mathcal{R}$. Finally, we assume that $H_{\eta}$ is smooth enough for a pseudo differential to be defined, and that the following function is well-defined, for all $L > 0$

$$ \mathrm{Lip}_{H_{\eta}}(L) = b_M \sup_{\eta \le 1} \sup_{|p|,|q| \le L} \sup_{(d_p, d_q) \in \partial H_{\eta}(p, q)} (|d_p| + |d_q|), \quad (2.1) $$

where $b_M$ is defined in (A3). This function will be used for various stability conditions and error estimates.

Suppose that $H_{\eta}$ and $R_{\eta}$ enjoy the above properties, and let $T > 0$, and the number $N_T$ of time steps be given. The time step is defined as $\Delta t = T/N_T$, and let $t^n = n\Delta t$ for $n \in [0, N_T]$. The trait step is denoted
$\Delta x > 0$, and the grid is defined with $x_i = x_0 + i\Delta x$ for all $i \in \mathbb{Z}$. For $n \in [0, N_T - 1]$ and $i \in \mathbb{Z}$, the scheme for (HJ) is given by

$$ \left\{ \begin{array}{l} \displaystyle \frac{u_i^{n+1} - u_i^n}{\Delta t} + b(x_i, I^{n+1}) H_{\eta} \left( \frac{u_i^n - u_{i-1}^n}{\Delta x}, \frac{u_{i+1}^n - u_i^n}{\Delta x} \right) + R_{\eta} (t^{n+1}, x_i, I^{n+1}) = 0 \\[1.5ex] \displaystyle \min_{j \in \mathbb{Z}} u_j^{n+1} = 0. \end{array} \right. \qquad (S) $$

It is initialized with

$$ u_i^0 = u^{\text{in}}(x_i), \qquad (2.2) $$

for all $i \in \mathbb{Z}$, and it is such that $\min_{j \in \mathbb{Z}} u_j^0 = 0$, thanks to (A12).

**Definition 2.1.** A scheme is $\mathcal{L}$-well chosen if $R_{\eta}$ satisfies (A6)-(A7)-(A8)-(A9), and (A17), if $H_{\eta}$ satisfies (A13), (A14), (A15), (A16), is smooth enough for (2.1) to be well-defined and finite for all $0 \le L \le \mathcal{L}$, and if one of the following conditions holds:

* Convex setting: $H_{\eta}$ is convex on $[-\mathcal{L}, \mathcal{L}]^2$,

* Flat setting: for all $\eta \in (0, 1]$ and $p, q$ such that $-\mathcal{L} \le p \le 0 \le q \le \mathcal{L}$, $H_{\eta}(p, q) = 0$.

**Remark 2.1.** It is possible for a scheme to be $\mathcal{L}$-well chosen for all $\mathcal{L}$. In that case, it will simply be said that it is well-chosen.

**Remark 2.2.** When in the flat setting, (A15) is automatically satisfied. Indeed, if $p \ge 0$ we have $0 = H_{\eta}(0, p) \le H_{\eta}(p, p)$, and similarly if $p < 0$.

**Remark 2.3.** One can notice that the flat and convex settings are not mutually exclusive. Moreover, it is always possible, when in the convex setting, to build a scheme that is also in the flat setting. Indeed, since $H(p) \ge 0$, one can consider $\max(H_{\eta}(p, q), 0)$.

Focus now on the properties of the trait-time mesh and the choice of $\eta$. Since an explicit discretization for the derivative of $u$ with respect to trait $x$ in (HJ) is chosen, stability conditions (CFL) must be satisfied. They allow the monotony of the numerical scheme in most results. It assumes

$$ \Delta t \text{Lip}_{H_{\eta}}(\text{Lip}_{\text{in}} + KT) \le \Delta x, \quad \text{and} \quad \Delta t \text{Lip}_{H_{\eta}}(\bar{a} + KT) \le \Delta x, \qquad (2.3) $$

where $\text{Lip}_{\text{in}}$ is defined in (A10), $\bar{a}$ in (A11), and $K$ in (1.3). Note that $\Delta t \text{Lip}_{H_{\eta}}(\bar{a} + KT) \le \Delta x$ is then automatically satisfied since $\bar{a} \le \bar{a}$, where $\bar{a}$ is also defined in (A11). However, a stronger condition is needed to conclude the proof of the convergence of $(S)$ towards the viscosity solution of (HJ),

$$ \Delta t \text{Lip}_{H_{\eta}}(14(\text{Lip}_{\text{in}} + KT) + 1) \le \Delta x, \qquad (2.4) $$

so that

$$ \text{L}_{\max} = 14(\text{Lip}_{\text{in}} + KT) + 1, \quad \text{and} \quad C_{H_{\eta}} = \text{Lip}_{H_{\eta}}(\text{L}_{\max}), \qquad (2.5) $$

are defined, and that the ratio $\Delta x/\Delta t$ is fixed, such that

$$ \Delta t \max \left( 2\text{Lip}_{H_{\eta}}(\bar{a} + KT), 2\text{Lip}_{H_{\eta}}(\text{Lip}_{\text{in}} + KT), C_{H_{\eta}} \right) \le \Delta x. \qquad (2.6) $$

Note that in most cases, the maximum of the two last items is $C_{H_{\eta}}$, and that this condition is restrictive. However, it can be relaxed in practice, see Section 4.3. Any monotonous function of $\Delta t$ that tends to 0 when $\Delta t \to 0$ can be chosen to define $\eta$. In what follows,

$$ \eta = \min(\sqrt{\Delta x}, 1), \qquad (2.7) $$

will be considered, as this choice simplifies the estimates in the convergence proof. The cut-off for $\Delta x > 1$ is introduced to use the definition of $\text{Lip}_{H_{\eta}}(L)$ in (2.1). Finally, a restriction on the discretization step is necessary, that is

$$ \Delta x \le \frac{1}{K w_T \text{Lip}_{H_{\eta}}(\text{Lip}_{\text{in}} + KT)}, \quad \text{and} \quad \Delta x \le \frac{4}{(e^T + 1)^2}, \qquad (2.8) $$

with

$$ w_T = e^{\beta T} \left( \gamma T + \sup_{x \in \mathbb{R}} u^{\text{in}''}(x) \right) \qquad (2.9) $$
$$ \gamma = K + b_M \max \{-H_{\eta}(-\text{Lip}_{\text{in}} - KT, \text{Lip}_{\text{in}} + KT), H_{\eta}(\text{Lip}_{\text{in}} + KT, -\text{Lip}_{\text{in}} - KT)\} \quad (2.10) $$

$$ \beta = 4\text{Lip}_{H_{\eta}}(\text{Lip}_{\text{in}} + KT). \quad (2.11) $$

Note that these assumptions are mostly technical. Indeed, removing the first bound requires slightly more regularity on $\mathcal{R}$, as $I_m - 1$ must then be lowered in (A7). The second one gives rise to the 1 in the definition of $L_{\max}$ in (2.5).

Let us now introduce the notion of an adapted discretization, to refer to a discretization $\Delta t$, $\Delta x$, $\eta$ satisfying these constraints:

**Definition 2.2.** A triplet $(\Delta t, \Delta x, \eta)$ is adapted to $H_{\eta}$ and $R_{\eta}$ if $(H_{\eta}, R_{\eta})$ is a $L_{\max}$-well chosen scheme, and (2.6), (2.7), (2.8) are satisfied.

A discretization $H_{\eta}, R_{\eta}, \Delta t, \Delta x, \eta$ is adapted if $(\Delta t, \Delta x, \eta)$ is adapted to $(H_{\eta}, R_{\eta})$.

*Remark 2.4.* $(H_{\eta}, R_{\eta})$ is assumed to be $L_{\max}$-well chosen to ensure that $C_{H_{\eta}}$ is well-defined.

As scheme $(\mathcal{S})$ is implicit, let us start by stating an existence result, along with qualitative properties on $\mathbf{u} = (u_i^n)_{n \in [0, N_T], i \in \mathbb{Z}}$ and on $\mathbf{I} = (I^{n+1})_{n \in [0, N_T-1]}$.

**Proposition 2.3** (Scheme $(\mathcal{S})$: existence of solutions and qualitative properties). Let $H_{\eta}, R_{\eta}, \Delta t, \Delta x, \eta$ be an adapted discretization (Definition 2.2). Then, scheme $(\mathcal{S})$ is well-posed: there exists $\mathbf{u} \in \mathbb{R}^{\mathbb{Z} \times [1, N_T]}$, and $\mathbf{I} \in \mathbb{R}^{[1, N_T]}$, satisfying $(\mathcal{S})$, where $\mathbf{u}^0 = (u_i^0)_{i \in \mathbb{Z}}$. Moreover $\mathbf{u}$ and $\mathbf{I}$ satisfy the following properties:

(i). Lipschitz property: $\mathbf{u}^n$ is $\text{Lip}_x(t^n)-\text{Lipschitz}$,

$$ \forall n \in [0, N_T], \forall (i, j) \in \mathbb{Z}^2, |u_i^n - u_j^n| \le \text{Lip}_x(t^n) \Delta x |i - j|, \quad (2.12) $$

where $\text{Lip}_x(t) = \text{Lip}_{\text{in}} + tK$, and with $\text{Lip}_{\text{in}}$ defined in (A10), and $K$ in (1.3).

(ii). Bounds for $\mathbf{u}^n$: let $b_t = \underline{b} - tb_M H_{\eta}(\underline{a}, -\underline{a}) - tK$, $\bar{b}_t = \bar{b} - tb_M H_{\eta}(-\bar{a}, \bar{a}) + tK$, where $\underline{a}, \bar{a}, \underline{b}, \bar{b}$ are defined in (A11) and $b_M$ in (A3). Then

$$ \forall n \in [0, N_T], \forall i \in \mathbb{Z}, \underline{a}|x_i| + \underline{b}^{n} \le u_i^n \le \bar{a}|x_i| + \bar{b}^{n}. \quad (2.13) $$

(iii). Let $n \in [1, N_T]$ and $j \in \mathbb{Z}$ such that $u_j^n = 0$, then

$$ R_{\eta}(t^n, x_j, I^n) \ge 0. \quad (2.14) $$

(iv). Bounds for $\mathbf{I}$: for all $n \in [1, N_T]$, $I_m - 1 \le I^n \le I_M$, where $I_m$ and $I_M$ are defined in (A7).

(v). Bound from below for the jumps of $\mathbf{I}$: there exists a constant $\kappa > 0$ such that for all $n \in [1, N_t-1]$,

$$ I^{n+1} - I^n \ge -\kappa \Delta t. \quad (2.15) $$

Moreover, in the flat setting, $\kappa = K^2$, with $K$ defined in (1.3). In the convex setting, $\kappa = K^2 + K \Delta x w(T)$ $\text{Lip}_{H_{\eta}}(\text{Lip}_x(T))$ suits, with $\text{Lip}_{H_{\eta}}$ defined in (2.6), $\text{Lip}_x$ in item (i), and $w$ in item (vi).

(vi). Semi-concavity: in the convex setting, there exists a function $w: [0, T] \to \mathbb{R}$ such that for all $t \in [0, T]$,

$$ w(t) \le e^{\beta t} \left( \gamma t + \sup_{x \in \mathbb{R}} u^{\text{in}''}(x) \right), \quad (2.16) $$

with $\beta, \gamma$ defined in (2.11)-(2.10), and such that for all $n \in [0, N_T]$, and for all $i \in \mathbb{Z}$

$$ \frac{u_{i-1}^n - 2u_i^n + u_{i+1}^n}{\Delta x^2} \le w(t^n). \quad (2.17) $$

(vii). Bound from below for $H_{\eta}$: for all $n \in [0, N_T]$, and for all $i \in \mathbb{Z}$,

$$ H_{\eta} \left( \frac{u_i^n - u_{i-1}^n}{\Delta x}, \frac{u_{i+1}^n - u_i^n}{\Delta x} \right) \ge -\Delta x \text{Lip}_{H_{\eta}}(\text{Lip}_x(t^n)) w(t^n)/b_M, \quad (2.18) $$

with $\text{Lip}_{H_{\eta}}$ defined in (2.1), $\text{Lip}_x$ in item (i), $w$ in (vi) and $b_M$ in (A3).
Moreover, **I** is uniquely determined if

$$b_M H_\eta(-\mathrm{Lip}_x(T), \mathrm{Lip}_x(T)) + K^{-1} \ge 0, \qquad (2.19)$$

where $b_M$ is defined in (A4), $\mathrm{Lip}_x$ in item (i) and $K$ in (1.3).

**Remark 2.5.** In the flat setting, (vi) is a priori not true. However (vii) holds, since $H_\eta(p, q) \ge 0$ for all $p$ and $q$. This property is a consequence of the flat setting definition (Definition 2.1), and of (A13).

**Remark 2.6.** As a consequence of Remark 2.5, it is worth noticing that (2.19) is automatically satisfied in the flat setting. In the convex setting, $H_\eta(-\mathrm{Lip}_x(T), \mathrm{Lip}_x(T))$ may be non-positive, so that condition (2.19) means that $R_\eta$ must be sufficiently decreasing w.r.t. $I$ in comparison to $b$, to have $I^{n+1}$ uniquely determined by (S). We emphasize on the fact that (2.19) is a sufficient condition, and that it seems that in practice (A5)-(A6) are enough, see Section 4.1.6.

**Remark 2.7.** Property (vi) is the discrete equivalent of a classical result for Hamilton-Jacobi equations: they preserve an upper bound on the second derivative, see [17]. As in the continuous case, the convexity of $H_\eta$ is a crucial assumption for this result.

**Remark 2.8.** When $\mathcal{R}$ does not depend on $t$, the time-dependent Lagrange multiplier $I$ in (HJ) is expected to be a non-decreasing function in $BV(0, T)$, see [8]. In the simpler case studied in [7], this property was conserved at the discrete level. Here, it is in fact still true in the flat setting, but the convex setting is not enough to have the monotony of $I$, see Section 5.1. However the bound (2.17) gives an estimate on the decreasing part.

**Remark 2.9.** Even when it is not expected to be non-decreasing, $I$ defined in (HJ) is in $BV(0, T)$. Moreover, no more regularity is to be expected, since it can have jumps, see Section (5.1). However, the bound from below for the jumps of $I$ in (v) yields that the decreasing jumps cannot be of order larger than $\Delta t$. It suggests that, in the continuous case (HJ), $I$ can only have increasing jumps, or, equivalently, that it can only decrease continuously.

**Prop. 2.3** strongly relies on monotony properties of scheme (S), it is proved in Section 4.1. Let us now focus on the convergence of these approximate solutions ($u, I$) towards the unique viscosity solution of (HJ). As in [7], the proof of this convergence uses compactness arguments, so that we need to pass from the discrete solutions ($u, I$) of (S), to the solutions ($u, I$) of (HJ). To that extend, let us introduce $I_{\Delta t}$ the constant by part reconstruction of $I$. Formally:

$$\forall n \in [0, N_T - 1], \forall s \in (0, \Delta t], I_{\Delta t}(t^n + s) := I^{n+1}, \text{ and } I_{\Delta t}(0) = I^1. \qquad (2.20)$$

Note that the choice of $I_{\Delta t}(0)$ is made only for technical purposes and does not bear any biological signification. For the extension of $u$, let us introduce for all $s \in (0, \Delta t]$, and for all $f: \mathbb{R} \to \mathbb{R}$, the operator $M_s^I$

$$\forall x \in \mathbb{R}, \forall I \in \mathbb{R}, M_s^I(f)(x) = f(x) - s b(x, I) H_\eta \left( \frac{f(x) - f(x - \Delta x)}{\Delta x}, \frac{f(x + \Delta x) - f(x)}{\Delta x} \right), \qquad (2.21)$$

and let us define $(t, x) \mapsto u_{\Delta t}(t, x)$ on $[0, T] \times \mathbb{R}$, such that for all $n \in [0, N_T - 1]$, $s \in (0, \Delta t]$ and $x \in \mathbb{R}$,

$$u_{\Delta t}(t^n + s, x) = M_s^{I_{\Delta t}(t^n + s)}(u_{\Delta t}(t^n, \cdot))(x) - s R_\eta(t^n + s, x, I_{\Delta t}(t^n + s)), \qquad (2.22)$$

supplemented by the initial condition $u_{\Delta t}(0, \cdot) = u^{\mathrm{in}}$. One can notice that for all $n \in [0, N_T]$, $i \in \mathbb{Z}$,

$$u_{\Delta t}(t^n, x_i) = u_i^n, \qquad (2.23)$$

that is, that the approximate solution corresponds with the scheme on the grid. If $R_\eta$ does not depend on $t$, one may also notice that, for all $x \in \mathbb{R}$, $u_{\Delta t}(\cdot, x)$ is piecewise linear. The following result states the convergence of $(u_{\Delta t}, I_{\Delta t})$ defined in (2.22)-(2.20) towards the viscosity solution of (HJ).

**Proposition 2.4** (Scheme (S): convergence). Let $(H_\eta, R_\eta)$ be a well chosen scheme (Def. 2.1), and $(\Delta t^n, \Delta x^n, \eta^n)$ be a sequence of discretizations adapted to $(H_\eta, R_\eta)$ (Def. 2.2) with $\Delta t^n \to 0$ when $n \to \infty$. Let $(u_{\Delta t^n}, I_{\Delta t^n})_n$ defined as in (2.22)-(2.20)-(2.23). Then, the following results hold:

* $u_{\Delta t^n}$ converges locally uniformly towards $u \in C^0([0, T], \mathbb{R})$, i.e. for all compact $\Omega$ of $\mathbb{R}$

$$\|u_{\Delta t^n} - u\|_{L^\infty([0,T] \times \Omega)} \xrightarrow{n \to +\infty} 0,$$

where $u$ is defined as the viscosity solution of (HJ).
• $I_{\Delta t^n}$ converges pointwise on $(0, T)$ towards $I_0 \in BV(0, T)$, which satisfies

$$I_m \le I_0 \le I_M,$$

and $I_0 = I$ almost everywhere on $(0, T)$, where $I$ is defined as the Lagrange-multiplier in (HJ).

Remark 2.10. Scheme (S) is introduced on an infinite trait grid $(x_i)_{i \in \mathbb{Z}}$, and Prop. 2.4 is stated in the same framework. In practice, Scheme (S) can be implemented on a grid that is reduced at each time step, so that no approximation has to be introduced on the boundaries, and that Prop. 2.4 holds. This represents however an increase in terms of computational cost, as grids larger than necessary are to be handled with. To work with a fixed grid, one can approximate the values that are outside of the grid. This has been proposed and numerically tested in [7], but the convergence cannot be properly established with these assumptions. Another issue, when dealing with grids larger than necessary, is that stability condition (2.6) may become much more restrictive on a larger grid, especially when trying to deal with data that are only Locally Lipschitz. This latter can be dealt with by linearizing $H_\eta$ where slopes are stronger than necessary, see Section 4.3.

Prop. 2.4 states the convergence of scheme (S), but does not give any clue on the convergence rate. This is due to the fact that the proof relies on compactness arguments for the sequences $(u_{\Delta t})_{\Delta t \ge 0}$, and $(I_{\Delta t})_{\Delta t \ge 0}$, thanks to the stability estimates of Prop. 2.3. Their limits are then identified with the viscosity solution of (HJ), using viscosity procedure. This goes through an appropriate regularization of $I_{\Delta t}$, as Lipschitz regularity is needed for time-dependency of source terms in the standard Hamilton-Jacobi framework [3, 11, 17]. We refer to Section 4.2 for details.

# 3 An example: Lotka-Volterra integral equations

In previous section the presentation of the hypothesis on $H_\eta$ was deliberately broad so that scheme (S) could be used in various applications, as soon as the equation satisfies the hypothesis of Theorem 1.1, and the discretization is adapted. In this section, we focus on a model of population dynamics described by $(P_e)$, and propose some schemes for (HJ)-(1.2). However, in order for (A2) to be satisfied, we will rather denote

$$\mathcal{H}(p) = \int_{\mathbb{R}} \mathcal{K}(z) e^{-zp} dz - 1, \qquad (3.1)$$

and modify $\mathcal{R}$ accordingly. Here $\mathcal{K}(z)$ is even, non negative and of integral 1. To make sure that $\mathcal{H}$ is properly defined, it is also supposed that $z \mapsto \mathcal{K}(z)e^{-zp}$, and $z \mapsto |z|\mathcal{K}(z)e^{-zp}$, are integrable for all $p \in \mathbb{R}$. This is for instance the case when $\mathcal{K}$ is a Gaussian or a compactly supported kernel.

## 3.1 Examples of schemes

In this part, we assume that $\mathcal{R}$ can be analytically computed, so that no approximation is needed. The choice $R_\eta = \mathcal{R}$ is natural in this context, and we discuss here the choice of the approximation $H_\eta$ of $\mathcal{H}$. Several classical schemes are covered by our analysis. Assume that we dispose of an approximation of the integral (3.1), and thus of $\mathcal{H}_\eta$ such that for all $|p| \le L_{\max}$:

$$|\mathcal{H}_\eta(p) - \mathcal{H}(p)| \le \eta, \quad \text{and} \quad \mathcal{H}_\eta(0) = \mathcal{H}(0) = 0.$$

Suppose also that $\mathcal{H}_\eta$ is convex, even, and that it is increasing for $p > 0$, as is $\mathcal{H}$. For instance, considering a quadrature of $\mathcal{H}$ with a symmetric grid, and renormalized for $\mathcal{H}_\eta$ to be 0 at $p = 0$ suits, although other choices may also be adapted. The so-called flat setting (see Def 2.1) includes the scheme (2.5) of [12]

$$\forall (p, q) \in \mathbb{R}, \quad H_{\eta}^{\text{CL}}(p, q) = \mathcal{H}_{\eta}(p^{+}) + \mathcal{H}_{\eta}(-q^{-}), \qquad (H_{\eta}^{\text{CL}})$$

where $p^{+} = \max(0, p)$ is the positive part of $p$ and $q^{-} = \max(0, -q)$ the negative part of $q$. One issue with this choice of discretization is that for a local maximum (i.e $p > 0$, $q < 0$), one can have $H_{\eta}^{\text{CL}}(p, q) > \eta + \sup_{d \in [q, p]} \mathcal{H}_{\eta}(d)$. This means that the value of the Hamiltonian at local maxima is overestimated by up to a factor two, when the local maximum enjoys symmetry. To avoid this issue, another scheme is for instance used in [22]

$$\forall (p, q) \in \mathbb{R}^2, \quad H_{\eta}^{\text{upw}}(p, q) = \max(\mathcal{H}_{\eta}(p^{+}), \mathcal{H}_{\eta}(-q^{-})). \qquad (H_{\eta}^{\text{upw}})$$

In this case, the issue is the lack of regularity when $\mathcal{H}_\eta(p^+)$ and $\mathcal{H}_\eta(-q^-)$ are equal. Both schemes $(H_\eta^{\text{CL}})-(H_\eta^{\text{upw}})$ are in the flat and convex settings, as defined in Definition 2.1. In [12], a Lax-Friedrichs scheme is also introduced

$$\forall (p, q) \in \mathbb{R}^2, \quad H_{\eta}^{\text{LF}}(p, q) = \mathcal{H}_{\eta}\left(\frac{p+q}{2}\right) - \theta(q-p), \qquad (H_{\eta}^{\text{LF}})$$
with $\theta$ large enough, so that $H_\eta$ increasing in its first variable, and decreasing in the second one. This choice no longer satisfies the flat hypothesis, but is still in the convex setting. Another choice stems from the limit of a natural asymptotic-preserving scheme introduced in [6] and analyzed in Section 3.2. Formally, introduce $h_\eta^+(p)$, $h_\eta^-(p)$, some monotonous and convex functions such that

$$ \left| h_{\eta}^{-}(p) - \int_{\mathbb{R}^{-}} \mathcal{K}(z)e^{-zq}dz + \frac{1}{2} \right| \leq \eta/2, \quad \left| h_{\eta}^{+}(q) - \int_{\mathbb{R}^{+}} \mathcal{K}(z)e^{-zq}dz + \frac{1}{2} \right| \leq \eta/2, $$

and

$$ h_{\eta}^{+}(0) = h_{\eta}^{-}(0) = 0. $$

Roughly speaking $h_\eta^+$ and $h_\eta^-$ respectively approximate the positive and negative part of the integral. This choice is intended to provide monotony properties (A13) to the approximation of $\mathcal{H}$. Using these notations, let

$$ \forall (p, q) \in \mathbb{R}^2, \quad H_{\eta}^{\mathrm{P1}}(p, q) = h_{\eta}^{-}(p) + h_{\eta}^{+}(q). \qquad (H_{\eta}^{\mathrm{P1}}) $$

As intended, $H_\eta$ is convex and enjoys the wished monotony properties (A13). However, when it is true for the solution of the continuous model (HJ), the monotony of $I$ is not true for the solution of scheme $(S)$, see section 5.1. This a drawback of $(H_\eta^{\mathrm{P1}})$. Indeed, when the net growth rate $R$ does not depend on $t$, the monotony of $I$ is a property that should be preserved.

We propose in what follows, a concave-convex-split scheme. It is intended to enforce the flat setting, and preserve the monotony of $I$ when it is relevant. Assume that the choice of $\mathcal{H}_\eta$ discretizing the whole integral is consistent with the choices of $h_\eta^+$ and $h_\eta^-$ made above, that is: $\mathcal{H}_\eta(p) = h_\eta^+(p) + h_\eta^-(p)$. We make use of the smoothness of the viscosity solutions in the convex area using $(H_\eta^{\mathrm{upw}})$, and the lack of regularity in the concave area using $(H_\eta^{\mathrm{P1}})$. Namely, we propose

$$ \forall (p, q) \in \mathbb{R}^2, \quad H_{\eta}^{\mathrm{CCS}}(p, q) = \begin{cases} H_{\eta}^{\mathrm{upw}}(p, q) & \text{if } p < q \\ H_{\eta}^{\mathrm{P1}}(p, q) & \text{else.} \end{cases} \qquad (H_{\eta}^{\mathrm{CCS}}) $$

It is worth noticing that $(H_\eta^{\mathrm{CCS}})$ is the maximum of $(H_\eta^{\mathrm{upw}})$ and $(H_\eta^{\mathrm{P1}})$. It is then a convex function of $(p, q)$.

Let us close this section with a remark about the discretization of the integral defining $\mathcal{H}$ in (3.1). In several cases, for instance the Gaussian kernel, the bilateral Laplace transform of $\mathcal{K}$ is known analytically, thus the addition of the parameter $\eta$ may seem artificial. Obviously, the construction and analysis of the schemes covered by this article is possible when the constant $C_\mathcal{L}$ in (A16) is zero.

## 3.2 Asymptotic-preserving scheme

The main goal of this paper is to provide a general framework for the numerical analysis of constrained Hamilton-Jacobi equations such as (HJ). The justification we have in mind for these models, is that they often appear as long-time and small mutations limits of population dynamics models. In this section, we focus more precisely on Lotka-Volterra integral equations $(P_\epsilon)$, supplemented with an initial condition $v_\epsilon^{\mathrm{in}}$ which satisfies (A10)-(A11). We also suppose that there exists two constants $c_m^{\mathrm{in}}$, $c_M^{\mathrm{in}}$ such that

$$ c_m^{\mathrm{in}} \varepsilon \leq \min v_\varepsilon^{\mathrm{in}} \leq c_M^{\mathrm{in}} \varepsilon, \qquad (A18) $$

and we assume that, $b$, $\mathcal{R}$, and $R = \mathcal{R} + b$ in $(P_\epsilon)$ satisfy the same assumptions as the ones in (HJ), and that $\mathcal{K}$ is defined as in (3.1). It is worth noticing that, in comparison to (HJ), problem $(P_\epsilon)$ is no longer a coupled problem where the unknowns are a solution of a PDE, and a Lagrange multiplier associated to a constraint. Instead, $J_\epsilon$ is determined here with $v_\epsilon$ with the second line of $(P_\epsilon)$. It accounts for a weighted measure of the size of the population, where the weight $\psi$ is such that

$$ \exists (\psi_m, \psi_M) \in \mathbb{R}^2, \forall x \in \mathbb{R}^d, 0 < \psi_m \leq \psi(x) \leq \psi_M < +\infty, \text{ and } \psi \in W^{2,\infty}(\mathbb{R}^d). \qquad (A19) $$

The following result holds

**Theorem 3.1.** ([4, 8]) Assume that the assumptions of Theorem 1.1 are satisfied, as well as (A18)-(A19), and that $\mathcal{R}$ does not depend on $t$. Let $v_\epsilon$ be the solution of $(P_\epsilon)$, and $J_\epsilon$ be defined in $(P_\epsilon)$. Suppose also that $(v_\epsilon^{\mathrm{in}})_\epsilon$ is a sequence of uniformly continuous functions which converges uniformly to $u^{\mathrm{in}}$. Then, $(v_\epsilon)_\epsilon$ converges locally uniformly to a function $u \in C^0([0, +\infty[ \times \mathbb{R}^d)$, and $(J_\epsilon)_\epsilon$ converges almost everywhere to a function $I$, such that $I \in BV(0, T)$ for all $T > 0$, and that $(u, I)$ is the unique viscosity solution of (HJ), with $\mathcal{H}$ defined in (3.1), and $\mathcal{R}$ replaced by $b + \mathcal{R}$.
Remark 3.1. As for Theorem 1.1, the above convergence result is a priori only true when $\mathcal{R}$ does not depend on $t$.

However, focusing on the asymptotic regime may be not completely relevant in terms of modeling. Considering the problem at some distance of the asymptotic regime may then give more information. Another issue when dealing only with asymptotic regimes such as (HJ) is, that it restricts the possibilities of biological situations that can be represented. Indeed, such asymptotic regimes are obtained under strong hypothesis on the birth rate $b$ and the net growth rate $\mathcal{R}$, that make sure that the population never extincts, nor grows uncontrolled. The question of the asymptotic behavior of the population in a regime of long time and small mutations, under more general hypothesis on the parameters of the model, has been addressed in the parabolic case for some changing environments in recent works [18, 19, 10, 39]. Up to our knowledge, it is an open problem for Lotka-Volterra integral equation ($P_\varepsilon$). If some of the assumptions on $b$ and $\mathcal{R}$ were not satisfied, Theorem 3.1 might not hold. Indeed, the population may vanish or explode, depending on the fact that the environment is poorly adapted or too advantageous. Very formally, suppose that it is possible to show that ($v_\varepsilon$)$_\varepsilon$ has a limit when $\varepsilon \to 0$. Then, these behaviors can be understood with the asymptotic behaviors of (min $v_\varepsilon$)$_\varepsilon$, and consequently of ($J_\varepsilon$)$_\varepsilon$:

* if min $v_\varepsilon$ converges to a positive limit, then $J_\varepsilon \to_{\varepsilon \to 0} 0$, meaning that the population vanishes,

* conversely, the population explodes if the limit of min $v_\varepsilon$ is negative, meaning that $J_\varepsilon \to_{\varepsilon \to 0} +\infty$.

The numerical resolution of ($P_\varepsilon$) requires specific attention, as ($P_\varepsilon$) becomes stiff when $\varepsilon$ is small. If no specific strategy was employed, the convergence of the numerical scheme for ($P_\varepsilon$) would fail in the asymptotic regime. Schemes specifically designed for such singular problems are said to be *Asymptotic Preserving (AP)*. Their properties are often summarized in the following diagram

$$ \begin{array}{ccc} \begin{array}{c} (P_\varepsilon) \\ \updownarrow \\ (S_\varepsilon^h) \end{array} & \xrightarrow{\varepsilon \to 0} & \begin{array}{c} (\text{HJ}) \\ \updownarrow \\ (S_0^h) \end{array} \\ h \to 0 & \longrightarrow & h \to 0. \end{array} \quad (3.2) $$

The first line represents the continuous level: the solution of ($P_\varepsilon$) converges when $\varepsilon \to 0$ to the viscosity solution of (HJ). The second line refers to the numerical schemes. The scheme ($S_\varepsilon^h$), where $h$ summarizes all the discretization parameters, is required to converge to the solution of ($P_\varepsilon$) when $\varepsilon > 0$ is fixed and $h \to 0$. Moreover, it must degenerate when $\varepsilon \to 0$ with fixed $h$ to another scheme ($S_0^h$), which approximates properly the viscosity (HJ) when $h \to 0$. An AP scheme can also enjoy the stronger property of being *Uniformly Accurate* (UA), meaning that its precision is independent on $\varepsilon$. AP schemes were introduced for kinetic equations [27, 28, 24], and various asymptotics have been considered [25, 16].

The design of AP schemes for Hamilton-Jacobi limits of biological models has been investigated more recently. An AP scheme for a Hamilton-Jacobi limit of a kinetic equation has been proposed in [23], and a model structured in age and phenotypic trait but in which no mutations are considered is treated in [1]. Regarding Lotka-Volterra models, an AP scheme for the parabolic case has been proposed and analyzed in [7]. A finite-differences based scheme which enjoys stability properties in the limit $\varepsilon \to 0$ has been proposed in [6], for a problem close to ($P_\varepsilon$) and in dimension $d=1$. It can be easily adapted in a scheme for ($P_\varepsilon$). Namely, let $\Delta z > 0$, and define for all $k \in \mathbb{Z}$, $z_k = k\Delta z$. For the other variables, we use the grids of scheme (S). To write the scheme in the spirit of the notation (2.21), we define, for all $J \in \mathbb{R}$, a function $M_\varepsilon^J: \mathbb{R}^J \to \mathbb{R}^J$, such that, for all $v = (v_i)_{i \in \mathbb{Z}}$, and for all $i \in \mathbb{Z}$,

$$ \begin{aligned} \mathcal{M}_\varepsilon^J(v)_i = v_i & - \Delta t \Delta z \sum_{k \in \mathbb{Z}, \varepsilon z_k \in [0, \Delta x]} \mathcal{K}(z_k) b(x_i + \varepsilon z_k, J) e^{-z_k(v_{i+1} - v_i)/\Delta x} & (3.3) \\ & - \Delta t \Delta z \sum_{k \in \mathbb{Z}, \varepsilon z_k \in [-\Delta x, 0)} \mathcal{K}(z_k) b(x_i + \varepsilon z_k, J) e^{-z_k(v_i - v_{i-1})/\Delta x} \\ & - \Delta t \Delta z \sum_{k \in \mathbb{Z}, |\varepsilon z_k| > \Delta x} \mathcal{K}(z_k) b(x_i + \varepsilon z_k, J) e^{-(\tilde{v}_{i,k} - v_i)/\varepsilon}, \end{aligned} $$

where $\tilde{v}_{i,k}$ denotes the linear interpolation of $(x, v)$ at abscissa $x_i + \varepsilon z_k$.

Remark 3.2. As $\tilde{v}_{i,k}$ is a linear interpolation of $v$, there is formally no difference between the difference quotients $(\tilde{v}_{i,k} - v_i)/\varepsilon$ and $z_k(v_{i+1} - v_i)/\Delta x$, if $0 < \varepsilon z_k \le \Delta x$ (or $z_k(v_i - v_{i-1})/\Delta x$, if $-\Delta x \le \varepsilon z_k < 0$).
In fact, they both are the slope of the linear interpolation of $(\boldsymbol{x}, \boldsymbol{v})$. The three lines of (3.3) may then seem artificial at first sight. However, when implemented, this reformulation ensures stability in the limit $\varepsilon \to 0$. Indeed, the expression $(\tilde{v}_{i,k} - v_i)/\varepsilon$ does not behave well when $\varepsilon \to 0$, as approximation errors in the linear interpolation are excessively increased, when divided by $\varepsilon \ll 1$, and injected in the exponential.

In order to state monotony properties for $\mathcal{M}_{\varepsilon}^J$, let us define for all $L \ge 0$,

$$ \mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}(L) = \Delta z b_M \sum_{k \in \mathbb{Z}} |z_k| \mathcal{K}(z_k) e^{L|z_k|}, \qquad (3.4) $$

where $b_M$ is defined in (A3). Note that $\mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}$ is well-defined since it was supposed that $z \mapsto |z|\mathcal{K}(z)e^{-zp}$, is integrable for all $p \in \mathbb{R}$. The monotony of $\mathcal{M}_{\varepsilon}^J$ is stated in the following lemma

**Lemma 3.2.** Let $L \ge 0$, $J \in \mathbb{R}$, and let $\boldsymbol{v}$ and $\boldsymbol{w} \in \mathbb{R}^Z$ such that for all $i \in \mathbb{Z}$, $|v_{i+1} - v_i| \le L\Delta x$, and $|w_{i+1} - w_i| \le L\Delta x$. Moreover, assume that $\Delta t \mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}(L) \le \Delta x$, with $\mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}$ defined in (3.4). We have,

(i). if $v_i \le w_i$ for all $i \in \mathbb{Z}$, then for all $i \in \mathbb{Z}$, $\mathcal{M}_{\varepsilon}^J(\boldsymbol{v})_i \le \mathcal{M}_{\varepsilon}^J(\boldsymbol{w})_i$,

(ii). If $\boldsymbol{v} - \boldsymbol{w} \in \ell^{\infty}(\mathbb{Z})$, then, $\mathcal{M}_{\varepsilon}^J(\boldsymbol{v}) - \mathcal{M}_{\varepsilon}^J(\boldsymbol{w}) \in \ell^{\infty}(\mathbb{Z})$, and $\|\mathcal{M}_{\varepsilon}^J(\boldsymbol{v}) - \mathcal{M}_{\varepsilon}^J(\boldsymbol{w})\|_{\infty} \le \|\boldsymbol{v} - \boldsymbol{w}\|_{\infty}$.

(iii). for all $i \in \mathbb{Z}$, $|\mathcal{M}_{\varepsilon}^J(\boldsymbol{v})_{i+1} - \mathcal{M}_{\varepsilon}^J(\boldsymbol{v})_i| \le L\Delta x$.

(iv). For all $i \in \mathbb{Z}$, $J \mapsto \mathcal{M}_{\varepsilon}^J(\boldsymbol{v})_i$ is non-decreasing.

*Proof.* The proof of this proposition is straightforward once one notices that the linear interpolation preserves inequalities. Namely, it means that for all $(i, k) \in \mathbb{Z}^2$, $\tilde{v}_{i,k} \le \tilde{w}_{i,k}$. Then, condition $\Delta t \mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}(L) \le \Delta x$ makes $\mathcal{M}_{\varepsilon}^J$ be a non-decreasing function of all its arguments, and the first point of the Lemma comes immediately. The second point is a consequence of the first one, of the inequalities $w_i - \|\boldsymbol{v} - \boldsymbol{w}\|_{\infty} \le v_i \le w_i + \|\boldsymbol{v} - \boldsymbol{w}\|_{\infty}$, and of

$$ \forall c \in \mathbb{R}, \forall i \in \mathbb{Z}, \mathcal{M}_{\varepsilon}^{J}(\boldsymbol{v} + c) = \mathcal{M}_{\varepsilon}^{J}(\boldsymbol{v}) + c. $$

The third point is done similarly with $w_i = v_{i+1} \pm L\Delta x$, while the last point comes from (A5), since $\mathcal{K}$ is non-negative. $\square$

Using notation (3.3), the scheme for $(P_{\varepsilon})$ is then defined for all $n \in [0, N_T - 1]$ and all $i \in \mathbb{Z}$, by

$$ \left\{ \begin{aligned} v_i^{n+1} &= \mathcal{M}_{\varepsilon}^{J^{n+1}}(\boldsymbol{v}^n)_i - \Delta t \mathcal{R}(t^{n+1}, x_i, J^{n+1}) \\ J^{n+1} &= \Delta x \sum_{i \in \mathbb{Z}} \psi(x_i) e^{-v_i^{n+1}/\varepsilon}. \end{aligned} \right. \qquad (S_{\varepsilon}) $$

It has been constructed with a quadrature in the integral kernel representing the births in $(P_{\varepsilon})$, and the stiffest term, $J$, is taken implicit to ensure stability in the regime $\varepsilon \to 0$. In what follows, it will be supposed that $\Delta z$ is not too large, so that the quadrature in $\mathcal{M}_{\varepsilon}^J$ is accurate enough. In particular, we will assume that

$$ \left| \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) - 1 \right| \le \frac{1}{2b_M}, \qquad (3.5) $$

with $b_M$ defined in (A3). Note that both $(\boldsymbol{v}^{n+1})_{n \in [0, N_T - 1]}$ and $\mathbf{J}$ depend on $\varepsilon$, although this dependency is omitted to simplify notations. The following proposition states the well-posedness of scheme $(S_{\varepsilon})$, and uniform w.r.t. $\varepsilon$ stability properties

**Proposition 3.3** (Scheme $(S_{\varepsilon})$: existence of solutions and stability properties). Suppose that $v_{\varepsilon}^{\mathrm{in}}$ satisfies (A10)-(A11)-(A12)-(A18), that $\psi$ satisfies (A19), that $b$ satisfies (A3)-(A4)-(A5), that $R$ and $\mathcal{R} = R + b$ satisfy (A6)-(A7)-(A8)-(A9), and that $\Delta t$ and $\Delta x$ are fixed such that

$$ \Delta x \mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}(\mathrm{Lip}_x(T)) \le \Delta t, \qquad (3.6) $$

with $\mathrm{Lip}_{H_{\eta}^{\mathrm{AP}}}$ defined in (3.4), $\mathrm{Lip}_x$ defined in Prop. 2.3-(i). Then, scheme $(S_{\varepsilon})$ is well-posed: there exists $\boldsymbol{v} \in \mathbb{R}^{\mathbb{Z} \times [0, N_T]}$ satisfying $(S_{\varepsilon})$. Moreover, there exists an $\varepsilon_0 > 0$, depending only on the constants arising in the assumptions, and on $\Delta t$ and $\Delta x$, such that for all $\varepsilon \in (0, \varepsilon_0)$, the sequence $(\boldsymbol{v}^n)_{n \in [0, N_T]}$ satisfies:
(i). Uniform Lipschitz continuity in trait: For all $n \in [0, N_T]$, $v^n$ enjoys Lip$_x(t^n)$-Lipschitz property:

$$ \forall i \in \mathbb{Z}, \left| \frac{v_{i+1}^n - v_i^n}{\Delta x} \right| \le \text{Lip}_x(t^n), $$

with Lip$_x$ defined in Prop. 2.3-(i).

(ii). Uniform bound from below for $v^n$: For all $n \in [0, N_T]$, for all $i \in \mathbb{Z}$,

$$ v_i^n \ge \underline{a}|x_i| + \underline{\beta}_{t^n}, $$

with $\underline{a}$ defined in (A11), and for all $t \in [0, T]$,

$$ \underline{\beta}_t = \underline{b} - b_M t \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) e^{a|x_i|} - tK, $$

with $\underline{b}$, $b_M$ and $K$ defined in (A11), (A3) and (1.3).

(iii). Uniform bounds for $J$: For all $n \in [0, N_T - 1]$,

$$ I_m/2 \le J^{n+1} \le 2I_M, $$

where $I_m$ and $I_M$ are defined in (A7).

(iv). Estimate for $\min v^n$: There exists $c_m$ and $c_M$ such that for all $n \in [0, N_T]$,

$$ c_m \varepsilon \le \min v^n \le c_M \varepsilon, $$

and $c_m \le c_m^{\text{in}}$, $c_M \ge c_M^{\text{in}}$ depend only on the constants defined in the assumptions, and on $\Delta x$ and $\Delta t$.

This proposition and its proof are very close to the stability properties proved in [7] for an AP scheme designed for parabolic Lotka-Volterra equations. The proof is adapted to the present case in Appendix A.1. Using the results of Prop. 3.3, the following proposition holds, that describes the asymptotic behaviour of $(S_\varepsilon)$ when $\varepsilon \to 0$. It is proved in Appendix A.2.

**Proposition 3.4** (Convergence of the scheme $(S_\varepsilon)$ when $\varepsilon \to 0$). Suppose that the assumptions of Theorem 3.1 are satisfied, and that $\Delta t$ and $\Delta x$ are fixed such that (3.6) holds. Let $v$ and $J$ be the $\varepsilon$-dependent sequences defined by $(S_\varepsilon)$. Then, for all $n \in [0, N_T - 1]$,

$$ \|v^{n+1} - u^{n+1}\|_{\infty} \underset{\varepsilon \to 0}{\longrightarrow} 0, \quad \text{and} \quad J^{n+1} \underset{\varepsilon \to 0}{\longrightarrow} I^{n+1}, $$

where $u$ and $I$ satisfy

$$ \left\{ \begin{array}{l} u_i^{n+1} = \mathcal{M}_0^{I^{n+1}}(u^n)_i - \Delta t \mathcal{R}(t^{n+1}, x_i, I^{n+1}) \\ \min u^{n+1} = 0, \end{array} \right. \qquad (S_0) $$

and for all $I \in \mathbb{R}$, for all $u = (u_i)_{i \in \mathbb{Z}}$ and for all $i \in \mathbb{Z}$, $\mathcal{M}_0^I : \mathbb{R}^\mathbb{Z} \to \mathbb{R}^\mathbb{Z}$ is defined by

$$ \mathcal{M}_0^I(\mathbf{u})_i = u_i - b(x_i, I) \Delta t \Delta z \left( \sum_{k \ge 0} \mathcal{K}(z_k) e^{-z_k(u_{i+1}-u_i)/\Delta x} + \sum_{k \ge -1} \mathcal{K}(z_k) e^{-z_k(u_i-u_{i-1})/\Delta x} \right). \quad (3.7) $$

**Remark 3.3.** The above proposition states the asymptotic behavior of scheme $(S_\varepsilon)$ when $\varepsilon \to 0$. It is easy to remark that the limit scheme $(S_0)$ is convergent towards the viscosity solution of the limit constrained Hamilton-Jacobi equation (HJ), with $\mathcal{H}$ defined in (3.1) and $\mathcal{R}$ replaced by $b+R$. Indeed, it can be rewritten with the formalism of $(S)$, with

$$ H_{\eta}^{\text{AP}}(p, q) = \Delta z \left( \mathcal{K}(0) + \sum_{k \le -1} \mathcal{K}(z_k) e^{-z_k p} + \sum_{k \ge 1} \mathcal{K}(z_k) e^{-z_k q} \right) - \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k), \quad (3.8) $$

and, if $R$ can be analytically computed,

$$ R_{\eta}^{\text{AP}}(t, x, I) = b(x, I) \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) + \mathcal{R}(t, x, I). \quad (3.9) $$

Notice then that $H_{\eta}^{\text{AP}}$ belongs to the class $(H_{\eta}^{\text{P1}})$ of choices proposed for $H_{\eta}$, in Section 3.1. Then, under the hypothesis of Prop. 3.4, the discretization $H_{\eta}^{\text{AP}}$, $R_{\eta}^{\text{AP}}$, $\Delta t$, $\Delta x$, $\eta$ is adapted if $\Delta z$ is small enough (Def. 2.2). In particular, $u$ and $I$ enjoy the properties of Prop. 2.3 and Prop. 2.4.
*Remark 3.4.* It is worth noticing that $M_0^I$, defined in (3.7), satisfies the properties of Lemma 3.2. As for $M_\epsilon^J$, defined in (3.3), this is a consequence of the monotony of $M_0^I$ in all its arguments, if the constraint (3.6) is satisfied.

*Remark 3.5.* Prop. 3.4 only holds if its assumptions are satisfied, that is if there is no extinction nor explosion of the population. However, numerical tests suggest that scheme $(S_\epsilon)$ is also stable when such behaviors are to be observed. In these cases, the asymptotics of $(P_\epsilon)$ can be conjectured by considering scheme $(S_\epsilon)$ with small $\epsilon$. We refer to section 5.3 for more details.

*Remark 3.6.* The scheme $(S_\epsilon)$ is defined for infinite grids in $x$ and $z$, that have to be truncated when implemented. The truncation in $z$ is done according to the decrease of $\mathcal{K}$. It is for instance easy when $\mathcal{K}$ is Gaussian or compactly supported. Once the grid in $z$ is truncated, scheme $(S_\epsilon)$ can be implemented on a grid in $x$ that is reduced at each time step, as in Remark 2.10. To reduce the computational cost, the values of $v$ outside of the grid may also be approximated, thanks to the coercivity of $v$, see Remark 2.10. In the numerical tests of Section 5, we use linear extrapolation of $v$ outside of the grid.

*Remark 3.7.* Prop. 3.4 gives the asymptotic behavior of scheme $(S_\epsilon)$ when $\epsilon \to 0$, and Prop. 2.4 gives the convergence of the limit scheme $(S_0)$. To complete the AP diagram (3.2), the convergence of scheme $(S_\epsilon)$ for fixed $\epsilon > 0$ is needed. Even though it yields tedious computations, this can be proved similarly as in [7], with stability estimates coming from the monotony properties (3.3), and truncation errors propagated with implicit function theorem. The main issue here is the truncation error, which is roughly defined as the error made by replacing derivatives by finite differences in the equation. If the solution $v_\epsilon$ of $(P_\epsilon)$ is smooth enough, it can be easily estimated with Taylor expansions. However, contrary to the parabolic case treated in [7], we do not have any result about the regularity of $v_\epsilon$. As a consequence, the convergence of scheme $(S_\epsilon)$ for fixed $\epsilon$ is true for smooth solutions with bounded derivatives, and remains an open question otherwise.

# 4 Convergence of scheme (S)

In this section, we come back to (HJ), to scheme (S), and to its reformulation (2.22), and we prove Prop. 2.3 and 2.4. These propositions strongly rely on monotony properties of scheme (S). Using the formalism of (2.22), this means that the operator $M_s^I$ in (2.21) is monotonic if a stability condition is satisfied. More precisely, the following result holds

**Lemma 4.1.** Let $I \in [I_m - 1, I_M]$, $s \in (0, \Delta t]$, $x \in \mathbb{R}$, $L > 0$, and assume that $\Delta t \text{Lip}_{H_\eta}(L) \le \Delta x$ we have:

(i). If $u_{\Delta t}$ and $v_{\Delta t}$ are such that $u_{\Delta t}(t^n, x) \le v_{\Delta t}(t^n, x)$, $u_{\Delta t}(t^n, x \pm \Delta x) \le v_{\Delta t}(t^n, x \pm \Delta x)$, and

$$|u_{\Delta t}(t^n, x \pm \Delta x) - u_{\Delta t}(t^n, x)| \le L\Delta x, \quad \text{and} \quad |v_{\Delta t}(t^n, x \pm \Delta x) - v_{\Delta t}(t^n, x)| \le L\Delta x,$$

$$\text{then } \mathcal{M}_s^I(u_{\Delta t}(t^n, \cdot))(x) \le \mathcal{M}_s^I(v_{\Delta t}(t^n, \cdot))(x), \text{ and } u_{\Delta t}(t^n + s, x) \le v_{\Delta t}(t^n + s, x).$$

(ii). If $u_{\Delta t}(t^n, \cdot)$ is L-Lipschitz, then $\mathcal{M}_s^I(u_{\Delta t}(t^n, \cdot))$ is also L-Lipschitz. If $I^{n+1} \in [I_m - 1, I_M]$, we have in addition that $u_{\Delta t}(t^n + s, \cdot)$ is $(L + sK)$-Lipschitz, with $K$ defined in (1.3).

*Proof.* Let $f$ be a L-Lipschitz function. Notice that $\mathcal{M}_s^I(f)(x)$ is increasing in $f(x)$, $f(x - \Delta x)$, $f(x + \Delta x)$ thanks to (2.1) so that the first point is straightforward. The second point, is a consequence of

$$\forall (x, y) \in \mathbb{R}^2, -L|x - y| + f(y) \le f(x) \le L|x - y| + f(y)$$

so that $\mathcal{M}_s^I(f)(x) - \mathcal{M}_s^I(f)(y) \le \mathcal{M}_s^I(L|\cdot - y| + f(y))(x) - \mathcal{M}_s^I(L|\cdot - y| + f(y))(y)$. Elementary computations give:

$$\mathcal{M}_s^I(L|\cdot - y| + f(y))(x) - \mathcal{M}_s^I(L|\cdot - y| + f(y))(y) \le L|x - y|,$$

and (A8) then yields the Lipschitz bound of $u_{\Delta t}(t^n + s, \cdot)$. $\square$

## 4.1 Proof of Prop. 2.3

In this section, the scheme is considered with formalism (S), and the proof of Prop. 2.3 is done by induction on the time step $n$. Thanks to (2.2), the properties of Prop. 2.3 are true for $n = 0$. Let us suppose that the hypothesis of Prop. 2.3 are satisfied and that $u^0, \dots, u^n \in \mathbb{R}^\mathbb{Z}$ satisfy (S)-(i)-(ii)-(iii) and (vi) if the convex setting is considered, and show that $u^{n+1}$ also satisfies these properties, that $I^{n+1}$ satisfies (iv) and (v), and that $u^n$ satisfies (vii).
### 4.1.1 $u^n$ satisfies property (vii)

First of all, remark that the fact that $u^n$ satisfies property (vii) is straightforward in the flat setting (see Remark 2.5), and that is a consequence of (vi) in the convex setting. Indeed, in the convex setting, for all $i \in \mathbb{Z}$, we have

$$ \frac{u_{i+1}^n - u_i^n}{\Delta x} = \frac{u_i^n - u_{i-1}^n}{\Delta x} + \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x} \le \frac{u_i^n - u_{i-1}^n}{\Delta x} + \Delta x w(t^n), $$

and the monotony of $H_\eta$, see (A13), and the stability condition (2.6) yield

$$ H_{\eta}\left(\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x},\frac{u_{i+1}^{n}-u_{i}^{n}}{\Delta x}\right) \geq H_{\eta}\left(\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x},\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x}+\Delta x w(t^{n})\right). $$

Property (vii) is then the consequence of an inequality of convexity for $H_\eta$, and of the definition of $\text{Lip}_{H_\eta}$ in (2.1), namely

$$ H_{\eta}\left(\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x},\frac{u_{i+1}^{n}-u_{i}^{n}}{\Delta x}\right) \geq H_{\eta}\left(\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x},\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x}\right)-\Delta x w(t^{n}) \text{ Lip}_{H_{\eta}}(\text{Lip}_{x}(t^{n}))/b_{M}. $$

Eventually, since $H_\eta(p, p) \ge 0$ for all $p$, see (A15), property (vii) holds.

### 4.1.2 Existence of $I^{n+1}$

The well posedness of (S) for a given $I^{n+1}$ is straightforward, since it is explicit in $u^n$. The only issue could stem from the implicit computation of $I$. Let us introduce, for any $I \in \mathbb{R}$, $u^{n+1, I}$ such that

$$ u_{i}^{n+1, I}=u_{i}^{n}-\Delta t b(x_{i}, I) H_{\eta}\left(\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x}, \frac{u_{i+1}^{n}-u_{i}^{n}}{\Delta x}\right)-\Delta t R_{\eta}\left(t^{n+1}, x_{i}, I\right), \quad (4.1) $$

and show that there exists $\mathcal{I} \in \mathbb{R}$, such that $\min_{i} u_{i}^{n+1, \mathcal{I}}=0$.

Remark 4.1. Coming back to Lemma 4.1, one can notice that if $u^n$, $v^n$ are such that $|u_{i+1}^n - u_i^n| \le L\Delta x$, $|v_{i+1}^n - v_i^n| \le L\Delta x$ for all $i \in \mathbb{Z}$, and if $\Delta t \text{Lip}_{H_\eta}(L) \le \Delta x$ then,

* If for all $i \in \mathbb{Z}$, $u_i^n \le v_i^n$, then $u_i^{n+1, I} \le v_i^{n+1, I}$ for all $i \in \mathbb{Z}$, and for all $I \in \mathbb{R}$,

* for all $i \in \mathbb{Z}$, for all $I \in [I_m - 1, I_M]$, $|u_{i+1}^{n+1, I} - u_i^{n+1, I}| \le (L + \Delta t K) \Delta x$.

Thanks to this remark, the bounds of $u^n$ in (2.13) can be propagated using the stability condition (2.6). It yields, for all $i \in \mathbb{Z}$, and for all $I \in \mathbb{R}$,

$$ \begin{aligned} & \underline{a}|x_i| + \bar{b}_{t^n} - \Delta t b_M H_\eta(\underline{a}, -\underline{a}) - \Delta t R_\eta(t^{n+1}, x_i, I) \\ & \le u_i^{n+1, I} \le \bar{a}|x_i| + \bar{b}_{t^n} - \Delta t b_M H_\eta(-\bar{a}, \bar{a}) - \Delta t R_\eta(t^{n+1}, x_i, I), \end{aligned} \quad (4.2) $$

where (A3) was used, considering that $H_\eta(-\bar{a}, \bar{a}) \le H_\eta(0, 0) = 0 \le H_\eta(\underline{a}, -\underline{a})$ thanks to (A13) and (A14). Then, (A8) yields the coercivity of $u^{n+1, I}$, so that

$$ \varphi: I \mapsto \min_{i \in \mathbb{Z}} u_i^{n+1, I}, \quad (4.3) $$

is well-defined for all $I \in \mathbb{R}$. Moreover, the minimum is in fact taken over a finite number of indices, thanks to (4.2) and $\varphi$ is hence continuous on all compact sets of $\mathbb{R}$. Let us now consider an index $j$ such that $u_j^n = \min u^n = 0$, then for all $I \in \mathbb{R}$,

$$ \begin{aligned} \varphi(I) &\le u_j^{n+1, I} = -\Delta t b(x_j, I) H_\eta\left(\frac{-u_{i-1}^n}{\Delta x}, \frac{u_{i+1}^n}{x}\right) - \Delta t R_\eta(t^{n+1}, x_i, I) \\ &\le -\Delta t b_M H_\eta(-\text{Lip}_x(T), \text{Lip}_x(T)) - \Delta t R_\eta(t^{n+1}, x_i, I), \end{aligned} $$

where the last inequality comes from (A3), (A13), and (2.12). Taylor expansion of $R_\eta(t^{n+1}, x_i, I)$ at $I_m$ and (A6) then yields for $I \le I_m$

$$ \varphi(I) \le -\Delta t b_M H_\eta(-\text{Lip}_x(T), \text{Lip}_x(T)) + \Delta t K^{-1}(I - I_m) - \Delta t R_\eta(t^{n+1}, x_i, I_m), $$
and we use (A7) to get

$$ \varphi(I) \leq -\Delta t b_M H_\eta (-\mathrm{Lip}_x(T), \mathrm{Lip}_x(T)) + \Delta t K^{-1} (I - I_m) \underset{I \to -\infty}{\longrightarrow} -\infty. $$

Hence, there exists $I^- \leq I_m$ such that $\varphi(I^-) < 0$. Coming back to (4.1), one has since $u_i^n \geq 0$ for all $i \in \mathbb{Z}$,

$$ u_i^{n+1, I} \geq -\Delta t b_M H_\eta (\mathrm{Lip}_x(T), -\mathrm{Lip}_x(T)) - \Delta t R_\eta (t^{n+1}, x_i, I), $$

where we used (A3), (A13) and (2.12). Then, as previously, (A6) yields for all $i \in \mathbb{Z}$, and for $I \geq I_M$

$$ u_i^{n+1, I} \geq -\Delta t b_M H_\eta (\mathrm{Lip}_x(T), \mathrm{Lip}_x(T)) + \Delta t (I - I_M) K^{-1} \underset{I \to +\infty}{\longrightarrow} +\infty, $$

so that there exists $I^+ \geq I_M$ such that $\varphi(I^+) > 0$. As $\varphi$ is continuous on $[I^-, I^+]$, there exists $\mathcal{I}$ such that $\varphi(\mathcal{I}) = 0$. Altough it may be not uniquely determined, $I^{n+1}$ is then well-defined, and so is $u^{n+1} = u^{n+1, I^{n+1}}$.

### 4.1.3 Properties (iii) and (iv)

The upper bound for $I^{n+1}$ is a consequence of (A7) and (A6). Indeed, since $u_i^n \geq 0$ for all $i \in \mathbb{Z}$, we use Remark 4.1, and (2.6)-(2.12)-(A14) to propagate this bound. It writes, for all $i \in \mathbb{Z}$

$$ u_i^{n+1} \geq -\Delta t R_\eta (t^{n+1}, x_i, I^{n+1}). $$

Considering now $j \in \mathbb{Z}$ such that $u_j^{n+1} = \min \mathbf{u}^{n+1} = 0$, then yields (iii). Moreover, remark that

$$ R_\eta (t^{n+1}, x_j, I^{n+1}) \geq 0 \geq R_\eta (t^{n+1}, x_j, I_M), $$

thanks to (A7), and use (A6) to conclude that $I^{n+1} \leq I_M$. To get the bound from below of $\mathcal{I}$, remark that as $\mathbf{u}^{n+1}$ satisfies (S), $u_i^{n+1} \geq 0$ for all $i \in \mathbb{Z}$. Let $j \in \mathbb{Z}$ such that $u_j^n = \min \mathbf{u}^n = 0$, we have

$$ u_j^{n+1} = -\Delta t b (x_j, I^{n+1}) H_\eta \left( \frac{u_j^n - u_{j-1}^n}{\Delta x}, \frac{u_{j+1}^n - u_j^n}{\Delta x} \right) - \Delta t R_\eta (t^{n+1}, x_j, I^{n+1}) \geq 0, $$

so that property (vii) yields

$$ R_\eta (t^{n+1}, x_j, I^{n+1}) \leq \Delta x w (t^n) \mathrm{Lip}_{H_\eta} (\mathrm{Lip}_x (t^n)). \quad (4.4) $$

Since $R_\eta (t^{n+1}, x_j, I_M) \geq 0$ (see (A7)), we obtain

$$ R_\eta (t^{n+1}, x_j, I^{n+1}) - R_\eta (t^{n+1}, x_j, I_M) \leq \Delta x w (t^n) \mathrm{Lip}_{H_\eta} (\mathrm{Lip}_x (t^n)). $$

Remark now that the result is straightforward if $I^{n+1} \geq I_m$, and suppose then that $I^{n+1} < I_m$. Since $R_\eta$ is decreasing w.r.t. $\mathcal{I}$, (A6) yields

$$ -K^{-1}(I^{n+1} - I_m) \leq \Delta x w(t^n) \mathrm{Lip}_{H_\eta}(\mathrm{Lip}_x(t^n)), \quad (4.5) $$

that is $I^{n+1} \geq I_m - 1$, thanks to (2.8).

Remark 4.2. It is worth noticing that the right-hand side of (4.4) can be replaced by 0 in the flat setting. In this case, the inequality $I^{n+1} \geq I_m$ holds, as in the continuous case, [4]. However, the convex setting does not preserve this qualitative property.

Remark 4.3. Note that (4.5) yields in fact a more precise bound from below for $I^{n+1}$,

$$ I^{n+1} \geq I_m - \Delta x w_T \mathrm{Lip}_{H_\eta}(\mathrm{Lip}_x(T)), \quad (4.6) $$

with $w_T$ defined by (2.9), and $\mathrm{Lip}_x$ in (i).
### 4.1.4 Property (v)

The proof of property (v) is very similar to the one of (iv). Starting from (4.4), we use property (iii) to get

$$R_{\eta}(t^{n+1}, x_j, I^{n+1}) - R_{\eta}(t^n, x_j, I^n) \leq \Delta x w(t^n) \text{ Lip}_{H_{\eta}}(\text{Lip}_x(t^n)),$$

where $j \in \mathbb{Z}$ is an index such that $u_j^n = \min \mathbf{u}^n = 0$. Introduce then $R_{\eta}(t^{n+1}, x_j, I^n)$ in the left-hand side, and use the Lipschitz-in-time regularity of $R_{\eta}$, see (A9), to get

$$R_{\eta}(t^{n+1}, x_j, I^{n+1}) - R_{\eta}(t^{n+1}, x_j, I^n) \leq K \Delta t + \Delta x w(t^n) \text{ Lip}_{H_{\eta}}(\text{Lip}_x(t^n)).$$

One can then notice as above that property (v) is straightforward if $I^{n+1} \geq I^n$, and use the fact that $R_{\eta}$ is decreasing w.r.t. $I$, (A6) to obtain if $I^{n+1} < I^n$

$$-K^{-1}(I^{n+1} - I^n) \leq K \Delta t + \Delta x w(t^n) \text{ Lip}_{H_{\eta}}(\text{Lip}_x(t^n)),$$

which is property (v).

*Remark 4.4.* The above result can be made more precise in the flat setting. Indeed, according to Remark 2.5, the right-hand side of (4.4) can be replaced by 0. It means that the approximation of the Hamiltonian $\mathcal{H}$ in (HJ) does not make $\mathbf{I}$ decrease, whereas it can happen in the convex case. This is the justification of the two expressions of $\kappa$ in property (v).

*Remark 4.5.* In the flat setting, and when $R_{\eta}$ does not depend on $t$, the previous computations yield that $I^{n+1} \geq I^n$. As a consequence, $\mathbf{I}$ is non-decreasing, as is $I$ defined in (HJ), see [8]. We emphasize on the fact that this qualitative property is not preserved in the convex setting, see Section 5.1.

To go on with the proof of Prop. 2.3, property (vi) must be proved in the convex setting (Def. 2.1). This is the purpose of the next section.

### 4.1.5 Property (vi) - convex setting

Suppose now that scheme ($\mathcal{S}$) is in the convex setting, so that $H_{\eta}$ is convex. To show that it preserves an upper bound on the discrete second derivative, let us introduce for all $i \in \mathbb{Z}$, $q_i^n = (u_{i+1}^n - u_i^n)/\Delta x$. Then $|q_i^n| \leq \text{Lip}_x(T)$ because of property (i), and

$$\begin{aligned} q_i^{n+1} = & q_i^n - \frac{\Delta t}{\Delta x} (b(x_{i+1}, I^{n+1}) H_{\eta}(q_i^n, q_{i+1}^n) - b(x_i, I^{n+1}) H_{\eta}(q_{i-1}^n, q_i^n)) \\ & - \frac{\Delta t}{\Delta x} (R_{\eta}(t^{n+1}, x_{i+1}, I^{n+1}) - R_{\eta}(t^{n+1}, x_i, I^{n+1})). \end{aligned} \quad (4.7)$$

In order to use the monotony of $H_{\eta}$ (see (A13)), remark that thanks to property (vii),

$$q_{i+1}^n \leq q_i^n + \Delta x w(t^n), \quad \text{and} \quad q_{i+2}^n \leq q_{i+1}^n + \Delta x w(t^n),$$

so that the following inequalities hold

$$H_{\eta}(q_{i-1}^n, q_i^n) \geq H_{\eta}(q_i^n - \Delta x w(t^n), q_i^n), \quad \text{and} \quad H_{\eta}(q_{i+1}^n, q_{i+2}^n) \geq H_{\eta}(q_{i+1}^n, q_{i+1}^n + \Delta x w(t^n)).$$

Using then the convexity of $H_{\eta}$, one has for all $(dp, dq) \in \partial H_{\eta}(q_i^n, q_{i+1}^n)$,

$$H_{\eta}(q_{i-1}^n, q_i^n) \geq H_{\eta}(q_i^n, q_{i+1}^n) - \Delta x w(t^n) dp - (q_{i+1}^n - q_i^n) dq \quad (4.8)$$

$$H_{\eta}(q_{i+1}^n, q_{i+2}^n) \geq H_{\eta}(q_i^n, q_{i+1}^n) + (q_{i+1}^n - q_i^n) dp + \Delta x w(t^n) dq, \quad (4.9)$$

and we inject these inequalities in $q_{i+1}^{n+1} - q_i^{n+1}$ expressed thanks to (4.7). We obtain

$$q_{i+1}^{n+1} - q_i^{n+1} \leq \mathfrak{A} + \mathfrak{B} + \mathfrak{C} + \mathfrak{D} + \mathfrak{E} \quad (4.10)$$

where

$$\mathfrak{A} = -\frac{\Delta t}{\Delta x} [R_{\eta}(t^{n+1}, x_{i+2}, I^{n+1}) - 2R_{\eta}(t^{n+1}, x_{i+1}, I^{n+1}) + R_{\eta}(t^{n+1}, x_i, I^{n+1})] \quad (4.11)$$

$$\mathfrak{B} = -\frac{\Delta t}{\Delta x} [b(x_{i+2}, I^{n+1}) - 2b(x_{i+1}, I^{n+1}) + b(x_i, I^{n+1})] H_{\eta}(q_i^n, q_{i+1}^n) \quad (4.12)$$
$$ \mathfrak{E} = -\frac{\Delta t}{\Delta x} [b(x_{i+2}, I^{n+1}) - b(x_i, I^{n+1})] \Delta x w(t^n) \, \mathrm{d}q $$

$$ \mathfrak{D} = -\frac{\Delta t}{\Delta x} [b(x_{i+2}, I^{n+1}) - b(x_i, I^{n+1})] (q_{i+1}^n - q_i^n) \, \mathrm{d}p $$

$$ \mathfrak{E} = \left[ \frac{\Delta t}{\Delta x} b(x_i, I^{n+1}) (dp - dq) \right] \Delta x w(t^n) + \left[ 1 - \frac{\Delta t}{\Delta x} b(x_i, I^{n+1}) (dp - dq) \right] (q_{i+1}^n - q_i^n), $$

and each term can be estimated separately. First of all, the $W^{2,\infty}$ estimate of $R_\eta$ in (A8), and the property $I_m - 1 \le I^{n+1} \le I_M$ yield

$$ |\mathfrak{A}| \le \Delta t \Delta x K, \qquad (4.13) $$

and similarly, (A4)-(2.1) give

$$ |\mathfrak{C}| \le 2\Delta t \Delta x \mathrm{Lip}_{H_\eta}(\mathrm{Lip}_x(T)) w(t^n). $$

The estimate of $\mathfrak{B}$ is a consequence of (A3) and of

$$ H_{\eta}(-\mathrm{Lip}_{x}(T), \mathrm{Lip}_{x}(T)) \le H_{\eta}(q_{i+1}^{n}, q_{i}^{n}) \le H_{\eta}(\mathrm{Lip}_{x}(T), -\mathrm{Lip}_{x}(T)), $$

so that

$$ |\mathfrak{B}| \le \Delta t \Delta x b_M \max(-H_{\eta}(-\mathrm{Lip}_x(T), \mathrm{Lip}_x(T)), H_{\eta}(\mathrm{Lip}_x(T), -\mathrm{Lip}_x(T))). \quad (4.14) $$

For the term $\mathfrak{D}$, it is worth noticing that since $H_\eta$ is increasing in its first variable, and decreasing in its second variable, then $dp - dq \ge 0$ for all $(dp, dq) \in \partial H_\eta (q_i^n, q_{i+1}^n)$. Moreover, using the stability condition (2.6) one has

$$ \forall (dp, dq) \in \partial H_{\eta}(q_{i}^{n}, q_{i+1}^{n}), \quad 0 \le \frac{\Delta t}{\Delta x} b(x_{i}, I^{n+1}) (dp - dq) \le 1, $$

meaning that $\mathfrak{E}$ is a convex combination of $\Delta x w(t^n)$ and $q_{i+1}^n - q_i^n$. It yields

$$ \mathfrak{E} \le \Delta x w(t^n). $$

However, $\mathfrak{D}$ cannot be estimated in all cases. Indeed, if $b(x_{i+2}, I^{n+1}) - b(x_i, I^{n+1}) \le 0$, it can be rewritten using (vii),

$$ \mathfrak{D} = \frac{\Delta t}{\Delta x} |b(x_{i+2}, I^{n+1}) - b(x_i, I^{n+1})| (q_{i+1}^n - q_i^n) dp \le 2\Delta t \Delta x \mathrm{Lip}_{H_\eta}(\mathrm{Lip}_x(T)) w(t^n), $$

but otherwise, the only estimate that can be used is $\mathfrak{D} \le 0$, if $q_{i+1}^n \ge q_i^n$. We can now gather these results to conclude that if $q_{i+1}^n \ge q_i^n$,

$$ \frac{q_{i+1}^{n+1} - q_i^{n+1}}{\Delta x} \le \Delta t \gamma + w(t^n) (1 + \beta \Delta t), \quad (4.15) $$

with $\gamma, \beta$ defined in (2.10)-(2.11). In the other case, namely when $q_{i+1}^n < q_i^n$, we reformulate (4.10) as

$$ q_{i+1}^{n+1} - q_i^{n+1} \le q_{i+1}^n - q_i^n + \mathfrak{A} + \mathfrak{B} + \mathfrak{F} + \mathfrak{G}, $$

where $\mathfrak{A}$ and $\mathfrak{B}$ have been defined in (4.11)-(4.12), and

$$ \mathfrak{F} = -\frac{\Delta t}{\Delta x} b(x_i, I^{n+1}) [H_{\eta}(q_{i-1}^n, q_i^n) - H_{\eta}(q_i^n, q_{i+1}^n)] $$

$$ \mathfrak{G} = -\frac{\Delta t}{\Delta x} b(x_{i+2}, I^{n+1}) [H_{\eta}(q_{i+1}^n, q_{i+2}^n) - H_{\eta}(q_i^n, q_{i+1}^n)]. $$

Use again (4.8)-(4.9), to get for all $(dp, dq) \in \partial H_\eta (q_i^n, q_{i+1}^n)$

$$ \mathfrak{F} \le \frac{\Delta t}{\Delta x} b(x_i, I^{n+1}) [\Delta x w(t^n) dp + (q_{i+1}^n - q_i^n) dq] $$

$$ \mathfrak{G} \le \frac{\Delta t}{\Delta x} b(x_{i+2}, I^{n+1}) [\Delta x w(t^n) (-dq) - (q_{i+1}^n - q_i^n) dp], $$

and use these estimates, the ones of $\mathfrak{A}, \mathfrak{B}$ in (4.13)-(4.14), the definition of $\gamma$ in (2.10) and the one of $\mathrm{Lip}_{H_\eta}$ in (2.1) to obtain

$$ q_{i+1}^{n+1} - q_i^{n+1} \le \left[ 1 - \frac{\Delta t}{\Delta x} \left( b(x_{i+2}, I^{n+1}) dp - b(x_i, I^{n+1}) dq \right) \right] (q_{i+1}^n - q_i^n) $$
$$+ \Delta t\gamma + 2\Delta t w(t^n) \mathrm{Lip}_{H_{\eta}}(\mathrm{Lip}_x(T)).$$

Eventually, the stability condition (2.6) yields that

$$1 - \frac{\Delta t}{\Delta x} (b(x_{i+2}, I^{n+1}) dp - b(x_i, I^{n+1}) dq) \geq 0,$$

and since $q_{i+1}^n - q_i^n < 0$, the above estimate can be simplified

$$\frac{q_{i+1}^{n+1} - q_i^{n+1}}{\Delta x} \leq \Delta t\gamma + w(t^n),$$

where the stability condition (2.6) also gave an upper bound of the last term. Hence, estimate (4.15) holds in fact in both cases, and

$$w(t^n) \leq e^{\beta t^n} (\gamma t^n + w(0)),$$

so that property (vi) is satisfied.

### 4.1.6 Conclusion

To conclude the induction, remark that properties (i) and (ii) are straightforward using Remark 4.1, item (iv) and (A8). Let us now discuss the fact that $I^{n+1}$ could be not uniquely determined as the solution of $\varphi(I) = 0$, where $\varphi$ is defined in (4.3). It is worth noticing that it is not a crucial issue in the proof above, since any suitable $I^{n+1}$ can be chosen in Section 4.1.2. Moreover, even if $I^{n+1}$ is not uniquely determined, Prop. 2.4 holds.

Actually, the fact that $I^{n+1}$ could be not uniquely determined is a purely numerical phenomenon, and it is only associated to the convex setting defined in Def. 2.1. Indeed, as $H_{\eta} \geq 0$, in the flat setting (see Remark 2.5), (A5) and (A6) yield that for all $i \in \mathbb{Z}$ and for all $n \in [0, N_t - 1]$,

$$I \mapsto u_i^{n+1, I}, \qquad (4.16)$$

with $u^{n+1, I}$ defined in (4.1), is increasing. As a consequence, if $I < J$ are such that $\min u^{n+1, I} = \min u^{n+1, J} = 0$, then for all $i \in \mathbb{Z}$,

$$u_i^{n+1, I} < u_i^{n+1, J}, \quad u_i^{n+1, I} \geq 0, \quad u_i^{n+1, J} \geq 0,$$

and considering an index $j$ such that $u_j^{n+1, J} = \min u^{n+1, J} = 0$ in the above expressions gives a contradiction. This remark also gives a sufficient condition for $I^{n+1}$ to be uniquely determined in the convex setting. Indeed, coming back to the definition of $u^{n+1, I}$ in (4.1), one can remark that the finite-differences in $H_{\eta}$ are bounded thanks to the Lipschitz bounds of $u^n$ stated in Prop. 2.3-(i). Hence, (4.16) is increasing if (2.19) is satisfied.

However, this condition is too restrictive in practice, since the bound of $H_{\eta}$ seems to be far too large around the minimal points of $u^{n+1, I}$. In fact, we observed in the numerical tests that $u^{n+1, I}$ is flat around its minimal points, although we were not able to quantify it precisely. As a consequence,

$$H_{\eta} \left( \frac{u_i^n - u_{i-1}^n}{\Delta x}, \frac{u_{i+1}^n - u_i^n}{\Delta x} \right),$$

is small around the minimal points of $u^n$, as are the finite differences in $H_{\eta}$. The condition (2.19) can then be relaxed in practice, and we could not exhibit a test case in which $I^{n+1}$ was not uniquely determined.

## 4.2 Proof of Prop. 2.4

In this section, we consider $u_{\Delta t}$ and $I_{\Delta t}$ defined by (2.20) and (S). As they are defined with a reformulation of scheme (S), such that $I_{\Delta t}$ is a constant by parts reconstruction of $I$, and that $u_{\Delta t}$ coincides with $u$ on the grid points, it is natural to expect $(u_{\Delta t}, I_{\Delta t})$ to satisfy similar properties as in Prop. 2.3. Indeed, the following result holds

**Lemma 4.2.** Let $H_{\eta}$, $R_{\eta}$, $\Delta t$, $\Delta x$, $\eta$ be an adapted discretization (Definition 2.2), and let $u_{\Delta t}$ and $I_{\Delta t}$ be defined by (2.22) and (2.20). We have:

(i). Lipschitz-in-space property: for all $t \in [0, T]$, $u_{\Delta t}(t, \cdot)$ is Lip$_x(t)$-Lipschitz, with Lip$_x$ defined in Prop. 2.3.

(ii). Lipschitz-in-time property: for all $x \in \mathbb{R}$, $u_{\Delta t}(\cdot, x)$ is Lip$_t$-Lipschitz with Lip$_t$ given by Lip$_t = b_M H_{\eta}(\mathrm{Lip}_x(T), -\mathrm{Lip}_x(T)) + K$, $b_M$ defined in (A3), and $K$ in (1.3).
(iii). Bounds for $u_{\Delta t}$: for all $(t, x) \in [0, T] \times \mathbb{R}$,

$$ \underline{a}|x| + \underline{b}_t \le u_{\Delta t}(t, x) \le \bar{a}|x| + \bar{b}_t, $$

with $\underline{a}$, $\bar{a}$ defined in (A11), and $\underline{b}_t$, $\bar{b}_t$ in Prop. 2.3.

(iv). Bounds for $I_{\Delta t}$: for all $t \in [0, T]$, $I_m - 1 \le I_{\Delta t}(t) \le I_M$, with $I_m$ and $I_M$ defined in (A7).

*Proof.* As $I_{\Delta t}$ is constant by parts and takes only the values of $I$, (iv) is a reformulation of Prop. 2.3-(iv). Items (i) and (iii) can then be proved by induction using Lemma 4.1, and similarly to what is done in the proof of Prop. 2.3. Item (ii) is a consequence of (i). More precisely, we only need to remark that for all $n \in [0, N_T - 1]$, $x \in \mathbb{R}$,

$$ \left| b(x, I^{n+1}) H_{\eta} \left( \frac{u_{\Delta t}(t^n, x) - u_{\Delta t}(t^n, x - \Delta x)}{\Delta x}, \frac{u_{\Delta t}(t^n, x + \Delta x) - u_{\Delta t}(t^n, x)}{\Delta x} \right) + R_{\eta}(t^{n+1}, x, I^{n+1}) \right| \le \text{Lip}_t. $$

This inequality is true thanks to (A3), (A13), (A8), and items (i) and (iv). $\square$

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

The proof of Prop. 2.4 then follows the lines of the proof of [7, Prop. 2.3], with technical adaptations to the case considered here. It is indeed a simple case of the general framework that is dealt with in this paper.

In Section 4.2.1, we show that $u_{\Delta t}$, and $I_{\Delta t}$ have limits when $\Delta t$, $\Delta x$, and $\eta$ go to 0, and some properties of these limits are stated. Section 4.2.2 is devoted to the regularization of $I$ and $I_{\Delta t}$, and associated numerical and viscosity solutions are introduced. Finally, Section 4.2.3 is devoted to the identification of the limits derived in Section 4.2.1, as the unique viscosity solution of (HJ).

### 4.2.1 Limits of $u_{\Delta t}$, and $I_{\Delta t}$

To start with, let us remark that Lemma 4.2 provides enough compactness to extract a limit. More precisely,

**Lemma 4.3.** Under the assumptions of Prop. 2.4, there exists $u_0 \in C^0([0, T], \mathbb{R})$, $I_0, I_0^+ \in BV(0, T)$, and a subsequence $(\Delta t^n)_n$ of the discretization such that $\Delta t^n \to_{n \to +\infty} 0$, and

(i). $u_{\Delta t^n}$ converges locally uniformly towards $u_0$: for all compact $\Omega$ of $\mathbb{R}$,

$$ \|u_{\Delta t^n} - u_0\|_{L^\infty([0,T] \times \Omega)} \xrightarrow{n \to +\infty} 0. $$

(ii). $u_0$ is Lip$_x(T)$-Lipschitz in trait, Lip$_t$-Lipschitz in time, with Lip$_x$ defined in Prop. 2.3-(i), and Lip$_t$ in Lemma 4.2-(ii). Moreover, $u_0$ satisfies:

$$ \forall t \in [0, T], \min_x u_0(t, x) = 0. $$

(iii). $(I_{\Delta t^n})_n$ converges almost everywhere on $(0, T)$ towards $I_0$ and $I_0^+$.

(iv). $I_0$ is lower semi-continuous, $I_0^+$ is upper semi-continuous, and both satisfy:

$$ I_m \le I_0 \le I_M, \quad \text{and} \quad I_m \le I_0^+ \le I_M. $$

*Proof.* This lemma is equivalent to [7, lemma 4.3]. The proof is similar, using Ascoli theorem for the convergence of $u_{\Delta t}$ in (i), and Helly's selection theorem for $I_{\Delta t}$ in (iii). Indeed, one can notice that Lemma 2.3-(iv)-(v) provides an uniform bound of $I_{\Delta t}$ in total variation. Actually, it is piecewise constant, equal to $I^{n+1}$ on $(t^n, t^{n+1})$ for all $n \in [0, N_T - 1]$, uniformly bounded on $(0, T)$, and with uniformly bounded total height of jumps down. Item (ii) is then a consequence of the local uniform convergence of the uniformly coercive and Lipschitz functions $(u_{\Delta t^n})_n$. Eventually, (iv) comes from (4.6). $\square$

*Remark 4.6.* In what follows, the mention $\Delta t \to 0$ will always refer to a subsequence for which the convergences of Lemma 4.3 hold true.

Note that neither $I_0$, nor $I_{\Delta t}$, are continuous on $(0, T)$, this behavior is showcased in Section 5.1. To be able to use the standard framework of Hamilton-Jacobi equations with Lipschitz Hamiltonian, we introduce in the next subsection a Lipschitz regularization of $I_0$, and $I_{\Delta t}$, and we regularize $u_0$ and $u_{\Delta t}$ accordingly.
### 4.2.2 Regularization of $I_{\Delta t}$ and $I_0$

Following [7, 8, 2] let us introduce, for all $\Delta t > 0$, $k > 0$, and $t \in [0, T]$:

$$I_{\Delta t}^{k,-} = \inf_{s \in [0,T]} I_{\Delta t}(s) + k|t-s|, \quad \text{and} \quad I_{\Delta t}^{k,+} = \sup_{s \in [0,T]} I_{\Delta t}(s) - k|t-s| \qquad (4.17)$$

Similarly, for $I_0$, we let for all $k > 0$, $t \in [0, T]$:

$$I_0^{k,-} = \inf_{s \in [0,T]} I_0(s) + k|t-s|, \quad \text{and} \quad I_0^{k,+} = \sup_{s \in [0,T]} I_0^+(s) - k|t-s|. \qquad (4.18)$$

These choices of regularization allow for stronger convergence properties than established in previous section. The needed results are provided by [7, Lemma 4.4], and they are recalled in the following lemma:

**Lemma 4.4.** Let $I_{\Delta t}^{k,-}$, $I_0^{k,-}$, $I_{\Delta t}^{k,+}$, $I_0^{k,+}$, and $I_{\Delta t}$ be defined by (4.17), (4.18) and (2.20). Suppose that the hypothesis of Lemma 4.3 are satisfied, so that $I_0$, $I_0^+$ are well-defined. Then, the following results hold:

(i). For all $\Delta t > 0$, for all $k > 0$,

$$\begin{aligned} I_m - 1 \le I_{\Delta t}^{k,-} &\le I_M, & I_m &\le I_0^{k,-} \le I_M, \\ I_m - 1 \le I_{\Delta t}^{k,+} &\le I_M, & I_m &\le I_0^{k,+} \le I_M, \end{aligned}$$

with $I_m$ and $I_M$ defined in (A7).

(ii). For all $\Delta t > 0$, for all $t \in [0, T]$, the following convergences hold when $k \to \infty$,

$$I_{\Delta t}^{k,-} \nearrow I_{\Delta t}, \quad I_0^{k,-} \nearrow I_0, \quad I_{\Delta t}^{k,+} \searrow I_{\Delta t}, \quad I_0^{k,+} \searrow I_0^+.$$

(iii). For all $\Delta t > 0$, for all $k > 0$, $I_{\Delta t}^{k,-}$, $I_0^{k,-}$, $I_{\Delta t}^{k,+}$ and $I_0^{k,+}$ are k-Lipschitz functions on $[0, T]$.

(iv). For all $k > 0$, $\|I_{\Delta t}^{k,-} - I_0^{k,-}\|_{\infty} \to_{\Delta t \to 0} 0$, and $\|I_{\Delta t}^{k,+} - I_0^{k,+}\|_{\infty} \to_{\Delta t \to 0} 0$.

The proof of this lemma is only computational, and partly done in [7]. The last point might seem surprising, but we emphasize on the fact that the convergences are not uniform in $k$.

Using these regularized versions of $I$ as parameters, one can then define associated viscosity solutions of (HJ) and solutions of scheme $(S)$. More precisely, for $t > 0$ and $x \in \mathbb{R}$, let us introduce $v^k$, and $w^k$ the viscosity solutions of

$$\partial_t v^k (t, x) + b(x, I_0^{k,-}(t)) \mathcal{H}(\nabla_x v^k)(t, x) + \mathcal{R}(t, x, I_0^{k,-}(t)) = 0 \qquad (4.19)$$

$$\partial_t w^k (t, x) + b(x, I_0^{k,+}(t)) \mathcal{H}(\nabla_x w^k)(t, x) + \mathcal{R}(t, x, I_0^{k,+}(t)) = 0, \qquad (4.20)$$

with initial condition $u^{\text{in}}$. Similarly, let $\Delta t > 0$, $n \in [0, N_T - 1]$, $s \in [0, \Delta t]$, and define $v_{\Delta t}^k$, and $w_{\Delta t}^k$ as in $(S)$, by

$$v_{\Delta t}^k (t^n + s, \cdot) = \mathcal{M}_s^{I_{\Delta t}^{k,-}(t^n+s)} (v_{\Delta t}^k (t^n, \cdot)) - sR_{\eta} (t^n + s, \cdot, I_{\Delta t}^{k,-} (t^n + s)) \qquad (4.21)$$

$$w_{\Delta t}^k (t^n + s, \cdot) = \mathcal{M}_s^{I_{\Delta t}^{k,+}(t^n+s)} (w_{\Delta t}^k (t^n, \cdot)) - sR_{\eta} (t^n + s, \cdot, I_{\Delta t}^{k,+} (t^n + s)), \qquad (4.22)$$

once again with initialization $u^{\text{in}}$. One can notice that, in the four cases, the constraint $\min u = 0$ is no longer enforced. Properties of $v^k$, $w^k$, $v_{\Delta t}^k$ and $w_{\Delta t}^k$ are stated in the following lemma,

**Lemma 4.5.** Suppose that the hypothesis of Prop. 2.3 are satisfied. Let $v^k$, $w^k$, $v_{\Delta t}^k$, $w_{\Delta t}^k$ and $u_{\Delta t}$ be defined by (4.19)-(4.20)-(4.21)-(4.22) and (2.22), and $I_0$, and $I_0^+$ defined in Lemma 4.3. Then,

(i). Uniform Lipschitz continuity in trait: for all $t \in [0, T]$, for all $k > 0$, for all $\Delta t > 0$, $v^k(t, \cdot)$, $w^k(t, \cdot)$, $v_{\Delta t}^k(t, \cdot)$, and $w_{\Delta t}^k(t, \cdot)$ are Lip$_x(T)$-Lipschitz, with Lip$_x$ defined in Prop. 2.3-(i).

(ii). Uniform Lipschitz continuity in time: for all $x \in \mathbb{R}$, for all $k > 0$, for all $\Delta t > 0$, $v^k(\cdot, x)$, $w^k(\cdot, x)$, $v_{\Delta t}^k(\cdot, x)$, and $w_{\Delta t}^k(\cdot, x)$ are Lip$_t$-Lipschitz, with Lip$_t$ defined in Lemmas 4.2-(ii).
(iii). Uniform bounds: for all $t \in [0, T]$, for all $x \in \mathbb{R}$, for all $k > 0$, for all $\Delta t > 0$,

$$
\begin{align*}
\underline{a}|x| + \underline{b}_t &\le v_{\Delta t}^k (t, x) \le \bar{a}|x| + \bar{b}_t \\
\underline{a}|x| + \underline{b}_t &\le w_{\Delta t}^k (t, x) \le \bar{a}|x| + \bar{b}_t \\
\underline{a}|x| + \underline{b}_t - t \left( K + \sup_{p \in [-\underline{a}, \underline{a}]} H(p) \right) &\le v^k (t, x) \le \bar{a}|x| + \bar{b}_t + tK \\
\underline{a}|x| + \underline{b}_t - t \left( K + \sup_{p \in [-\underline{a}, \underline{a}]} H(p) \right) &\le w^k (t, x) \le \bar{a}|x| + \bar{b}_t + tK,
\end{align*}
$$

with $\underline{a}, \bar{a}, \underline{b}, \bar{b}$ defined in (A3), $\underline{b}_t$ and $\bar{b}_t$ in Prop. 2.3-(ii), and $K$ in (1.3).

(iv). Monotony of the approximation: for $k \to \infty$ we have $v^k \nearrow v^\infty$, and $w^k \searrow w^\infty$, pointwise on $[0,T] \times \mathbb{R}$, where $v^\infty, w^\infty$ are viscosity solution of respectively

$$
\partial_t v^\infty + b(x, I_0) \mathcal{H}(\nabla_x v^\infty) + \mathcal{R}(t, x, I_0) = 0 \quad (4.23)
$$

$$
\partial_t w^\infty + b(x, I_0^+) \mathcal{H}(\nabla_x w^\infty) + \mathcal{R}(t, x, I_0^+) = 0, \quad (4.24)
$$

with initialization $u^{\text{in}}$.

(v). We have $v_{\Delta t}^k \le u_{\Delta t} \le w_{\Delta t}^k$.

Remark 4.7. Item (iv) and the fact that $I_0 = I_0^+$ a.e., yield that $v^\infty = w^\infty$ thanks to Theorem 1.1.

Proof. Items (i)-(ii) are classical properties of viscosity solutions, (iii) stems from a comparison principle, and (iv) from [8]. Eventually, (v) is a consequence of Lemma 4.1, and of the fact that $I \mapsto \mathcal{M}_s^I$ is non-decreasing, with $\mathcal{M}_s^I$ defined in (2.21). Indeed, if $v_{\Delta t}^k(t^n, \cdot) \le u_{\Delta t}(t^n, \cdot) \le w_{\Delta t}^k(t^n, \cdot)$, and since (2.6) is satisfied, then

$$
\mathcal{M}_s^{I_{\Delta t}^{k,-}(t^n)} v_{\Delta t}^k(t^n, \cdot) \le \mathcal{M}_s^{I_{\Delta t}^{k}(t^n)} v_{\Delta t}^k(t^n, \cdot) \le \mathcal{M}_s^{I_{\Delta t}^{k}(t^n)} u_{\Delta t}(t^n, \cdot),
$$

and

$$
\mathcal{M}_s^{I_{\Delta t}^{k}(t^n)} u_{\Delta t}(t^n, \cdot) \le \mathcal{M}_s^{I_{\Delta t}^{k}(t^n)} w_{\Delta t}^{k}(t^n, \cdot) \le \mathcal{M}_s^{I_{\Delta t}^{k,+}(t^n)} w_{\Delta t}^{k}(t^n, \cdot).
$$

These two chains of inequalities give the desired result thanks to (A6), and a straightforward induction. $\square$

4.2.3 (u₀, I₀) is the viscosity solution of (HJ)

In this section, we show that $v^\infty = w^\infty = u_0$, and that Prop. 2.4 holds. Following [12, 7], let us introduce for
$\alpha \in (0,1)$, $\sigma > 0$, and $(t, x, \tau, \xi) \in [0, T) \times \mathbb{R} \times [0, T] \times \mathbb{R}$,

$$
\begin{align}
\Psi_{\alpha, \sigma}^{-}(t, x, \tau, \xi) = v^{k}(t, x) - v_{\Delta t}^{k}(\tau, \xi) & - \frac{(x - \xi)^{2}}{2\sqrt{\Delta x}} - \frac{(t - \tau)^{2}}{2\sqrt{\Delta t}} \tag{4.25} \\
& - \alpha \frac{e^{t}}{2}(x^{2} + \xi^{2}) - \frac{\alpha}{T - t} - (\sigma + C_{H_{\eta}}^{2} \alpha e^{T}) t, \nonumber
\end{align}
$$

$$
\begin{align}
\Psi_{\alpha, \sigma}^{+}(t, x, \tau, \xi) ={}& w^{k}(t, x) - w_{\Delta t}^{k}(\tau, \xi) + \frac{(x-\xi)^{2}}{2\sqrt{\Delta x}} + \frac{(t-\tau)^{2}}{2\sqrt{\Delta t}} \tag{4.26} \\
& + \alpha \frac{e^{t}}{2}(x^{2}+\xi^{2}) + \frac{\alpha}{T-t} + (\sigma + C_{H_{\eta}}^{2}\alpha e^{T}) t, \nonumber
\end{align}
$$

with $C_{H_{\eta}}$ defined in (2.5). Thanks to Lemma 4.5, $\Psi_{\alpha, \sigma}^{-}$ and $\Psi_{\alpha, \sigma}^{+}$ respectively admit a maximum and minimum on $[0, T) \times \mathbb{R} \times [0, T] \times \mathbb{R}$, and the following lemma holds:

**Lemma 4.6.** Let $\Psi_{\alpha, \sigma}^{-}$ and $\Psi_{\alpha, \sigma}^{+}$ be defined respectively by (4.25), (4.26). We have the following properties:

* For all $\alpha \in (0,1)$, and for all $\sigma > 0$, $\Psi_{\alpha, \sigma}^{-}$ admits a maximum, and $\Psi_{\alpha, \sigma}^{+}$ a minimum on $[0, T) \times \mathbb{R} \times [0, T] \times \mathbb{R}$. They are reached respectively at $(t^{-}, x^{-}, \tau^{-}, \xi^{-})$ and $(t^{+}, x^{+}, \tau^{+}, \xi^{+})$.

* There exists $\sigma^{+} (\Delta t, k)$, $\sigma^{-} (\Delta t, k)$ positive, with $\sigma^{\pm} (\Delta t, k) \to_{\Delta t \to 0} 0$ for all $k > 0$ (though not uniformly in $k$), such that, for all $\alpha \in (0, 1)$, and $\Delta t$ small enough

$$
t^{\pm} \le 2\mathrm{Lip}_{t} \sqrt{\Delta t}. \qquad (4.27)
$$
Moreover,

$$
\begin{align}
\Psi_{\alpha, \sigma}^{-}(t^{-}, x^{-}, \tau^{-}, \xi^{-}) &\le 2\mathrm{Lip}_{t}^{2}\sqrt{\Delta t} + 10\mathrm{Lip}_{x}(T)^{2}\sqrt{\Delta x} \nonumber \\
\Psi_{\alpha, \sigma}^{+}(t^{+}, x^{+}, \tau^{+}, \xi^{+}) &\ge -2\mathrm{Lip}_{t}^{2}\sqrt{\Delta t} - 10\mathrm{Lip}_{x}(T)^{2}\sqrt{\Delta x}, \tag{4.28}
\end{align}
$$

where Lip$_x$ is defined in Prop. 2.3-(i), and Lip$_t$ in Prop. 4.2-(ii).

Proof. As it is the adaptation of the particular case [7, Lemma 4.6], the proofs of these two results are strongly related. It is however worth detailing it in what follows, since the approximations $H_\eta(p, p)$, and $R_\eta$ of $\mathcal{H}(p)$ and $\mathcal{R}$, and the general framework considered here brings technical difficulties. We focus here on the results for $\sigma^+(\Delta t, k)$ and $\Psi_{\alpha, \sigma}^+$. The others can be proved with straightforward adaptations.

The proof is divided in three steps, first show that

$$
\alpha e^{t^{+}} \max (|x^{+}|, |\xi^{+}|) \le 4 \mathrm{Lip}_{x}(T), \quad (4.29)
$$

and then that

$$
|t^+ - \tau^+| \le 2\mathrm{Lip}_t \sqrt{\Delta t}, \quad \text{and} \quad |x^+ - \xi^+| \le 10\mathrm{Lip}_x(T) \sqrt{\Delta x}. \tag{4.30}
$$

Finally, using these two results, define $\sigma^+(\tau, k)$, such that $\tau^+ = 0$. It yields $t^+ \le 2\mathrm{Lip}_t \sqrt{\Delta t}$, and deduce the desired bound on $\Psi_{\alpha, \sigma}^+$.

For the first step, notice that $\Psi_{\alpha, \sigma}^{+}(t^{+}, x^{+}, \tau^{+}, \xi^{+}) \le \Psi_{\alpha, \sigma}^{+}(t^{+}, 0, \tau^{+}, 0)$, and so that

$$
w^k(t^+, x^+) - w_{\Delta t}^k(\tau^+, \xi^+) + \alpha \frac{e^{t^+}}{2}(x^{+2} + \xi^{+2}) \le w^k(t^+, 0) - w_{\Delta t}^k(\tau^+, 0).
$$

Elementary computations and the Lipschitz-regularity of $w^k$, and $w_{\Delta t}^k$ provided by Lemma 4.5 then give (4.29)
Similarly, the bound on $(x^+ - \xi^+)$ is a consequence of $\Psi_{\alpha, \sigma}^+(t^+, x^+, \tau^+, \xi^+) \le \Psi_{\alpha, \sigma}^+(t^+, x^+, \tau^+, x^+)$. Thanks
to the non-negativity of $\xi^{+2}$, and after elementary computations, it yields

$$
\frac{(x^+ - \xi^+)^2}{2\sqrt{\Delta x}} \le w_{\Delta t}^k (\tau^+, \xi^+) - w_{\Delta t}^k (\tau^+, x^+) + \alpha \frac{e^{t^+}}{2} (x^{+2} - \xi^{+2}).
$$

Thus, (4.29), proved previously and Lemma 4.5, gives the desired result, since

$$
(x^+ - \xi^+)^2 \le 2\sqrt{\Delta x} \left[ \mathrm{Lip}_x(T) |x^+ - \xi^+| + \alpha \frac{e^{t^+}}{2} (x^+ - \xi^+) (x^+ + \xi^+) \right].
$$

For the other estimate in (4.30), use $\Psi_{\alpha, \sigma}^{+}(t^{+}, x^{+}, \tau^{+}, \xi^{+}) \le \Psi_{\alpha, \sigma}^{+}(t^{+}, x^{+}, t^{+}, \xi^{+})$ and Lemma 4.5 to get

$$
\frac{(t^+ - \tau^+)^2}{2\sqrt{\Delta t}} \le w_{\Delta t}^k (\tau^+, \xi^+) - w_{\Delta t}^k (t^+, \xi^+) \le \mathrm{Lip}_t |t^+ - \tau^+|.
$$

Let us now prove that $\tau^+ = 0$. We argue by contradiction, supposing that $\tau^+ > 0$, and defining $\sigma^+(\Delta t, k)$ such that this assumption cannot hold. Let us introduce for all $(t, x) \in [0, T) \times \mathbb{R}$,

$$
\varphi_{\alpha, \sigma}^{+}(t, x) = w_{\Delta t}^{k}(\tau^{+}, \xi^{+}) - \frac{(x - \xi^{+})^{2}}{2\sqrt{\Delta x}} - \frac{(t - \tau^{+})^{2}}{2\sqrt{\Delta t}} \\
- \alpha \frac{e^{t^{+}}}{2}(x^{+2} + \xi^{+2}) - \frac{\alpha}{T - t} - (\sigma + C_{H_{\eta}}^{2} \alpha e^{T}) t,
$$

so that $\Psi_{\alpha, \sigma}^{+}(t, x, \tau^{+}, \xi^{+}) = w^{k}(t, x) - \varphi_{\alpha, \sigma}^{+}(t, x)$. As $\varphi_{\alpha, \sigma}^{+}$ is minimal at $(t^{+}, x^{+})$, and since $w^{k}$ is the viscosity solution of (4.20), one has

$$
\partial_t \varphi_{\alpha, \sigma}^{+} (t^+, x^+) + b \left( x^+, I_0^{k,+} (t^+) \right) \mathcal{H} \left( \nabla_x \varphi_{\alpha, \sigma}^{+} (t^+, x^+) \right) + \mathcal{R} \left( t^+, x^+, I_0^{k,+} (t^+) \right) \ge 0, \quad (4.31)
$$

that can be reformulated as

$$
\begin{equation}
\begin{split}
& -\frac{t^+ - \tau^+}{\sqrt{\Delta t}} - \alpha \frac{e^{t^+}}{2} (x^{+2} + \xi^{+2}) - \frac{\alpha}{(T-t^+)^2} - \sigma - C_{H_\eta}^2 \alpha e^T \\
& \qquad + b \left(x^+, I_0^{k,+}(t^+)\right) \mathcal{H}\left(-\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+\right) + \mathcal{R}\left(t^+, x^+, I_0^{k,+}(t^+)\right) \ge 0.
\end{split}
\tag{4.32}
\end{equation}
$$
A similar inequality is obtained considering $w_{\Delta t}^k$ and scheme (4.22). As, for all $(\tau, \xi) \in [0, T] \times \mathbb{R}$, the inequality $\Psi_{\alpha, \sigma}^+(t^+, x^+, \tau, \xi) \ge \Psi_{\alpha, \sigma}^+(t^+, x^+, \tau^+, \xi^+)$ holds, one has

$$ \overline{w}_{\Delta t}^k (\tau, \xi) + k^+ \ge w_{\Delta t}^k (\tau, \xi), \qquad (4.33) $$

where for all $(\tau, \xi) \in [0, T] \times \mathbb{R}$,

$$ \overline{w}_{\Delta t}^k (\tau, \xi) = \frac{(x^+ - \xi)^2}{2\sqrt{\Delta x}} + \frac{(t^+ - \tau)^2}{2\sqrt{\Delta t}} + \alpha \frac{e^{t^+}}{2} \xi^2, \qquad (4.34) $$

and

$$ k^+ = w_{\Delta t}^k (\tau^+, \xi^+) - \frac{(x^+ - \xi^+)^2}{2\sqrt{\Delta x}} - \frac{(t^+ - \tau^+)^2}{2\sqrt{\Delta t}} - \alpha \frac{e^{t^+}}{2} \xi^2 = w_{\Delta t}^k (\tau^+, \xi^+) - \overline{w}_{\Delta t}^k (\tau^+, \xi^+). \qquad (4.35) $$

Recall now that $\tau^+ \in [0, T)$ is supposed such that $\tau^+ > 0$. Then, there exists $n^+ \in [0, N_T - 1]$, and $s^+ \in (0, \Delta t]$ such that $\tau^+ = t^{n^+} + s^+$. The next step consists in taking $\tau = t^{n^+}$ and $\xi = \xi^+$, and propagate the inequality (4.33) to time $\tau^+$ with Lemma 4.1. To that extend, we show that the arguments of $H_\eta$ in $M_{s^+}^{I_{\Delta t}^{k,+}(\tau^+)}(w_{\Delta t}^k(t^{n^+}), \cdot)(\xi^+)$ and $M_{s^+}^{I_{\Delta t}^{k,+}(\tau^+)}(\overline{w}_{\Delta t}^k(t^{n^+}), \cdot)(\xi^+)$ are bounded by $L_{\max}$. Concerning $w_{\Delta t}^k$, it is in fact a consequence of Lemma 4.5. Consider now $\overline{w}_{\Delta t}^k$, one has

$$ \frac{\overline{w}_{\Delta t}^k (\tau, \xi^+ \pm \Delta x) - \overline{w}_{\Delta t}^k (\tau, \xi^+)}{\Delta x} = \mp \frac{x^+ - \xi^+}{\sqrt{\Delta x}} + \frac{\sqrt{\Delta x}}{2} \pm \alpha e^{t^+} \xi^+ + \alpha e^{t^+} \frac{\Delta x}{2}, $$

so that (4.29)-(4.30), and (2.8) yield

$$ \left| \overline{w}_{\Delta t}^k (\tau, \xi^+ \pm \Delta x) - \overline{w}_{\Delta t}^k (\tau, \xi^+) \right| \le L_{\max} \Delta x. \qquad (4.36) $$

We may now compute the values of the image of each side of (4.33) by $M_{s^+}^{I_{\Delta t}^{k,+}}$. By definition of $w_{\Delta t}^k$, one has

$$ M_{s^+}^{I_{\Delta t}^{k,+}(\tau^+)} (w_{\Delta t}^k (t^{n^+}, \cdot)) (\xi^+) = w_{\Delta t}^k (\tau^+, \xi^+) + s^+ R_{\eta} (\tau^+, \xi^+, I_{\Delta t}^{k,+} (\tau^+)), \qquad (4.37) $$

and the expression of $M_s^I$ (see (2.21)) yields

$$ \begin{aligned} M_{s^+}^{I_{\Delta t}^{k,+}(\tau^+)} \left( \overline{w}_{\Delta t}^k (t^{n^+}, \cdot) + k^+ \right) (\xi^+) = & \overline{w}_{\Delta t}^k (t^{n^+}, \xi^+) + k^+ \\ -s^+ b (\xi^+, I_{\Delta t}^{k,+} (\tau^+)) H_{\eta} & \left( \frac{\overline{w}_{\Delta t}^k (t^{n^+}, \xi^+) - \overline{w}_{\Delta t}^k (t^{n^+}, \xi^+ - \Delta x)}{\Delta x}, \frac{\overline{w}_{\Delta t}^k (t^{n^+}, \xi^+ + \Delta x) - \overline{w}_{\Delta t}^k (t^{n^+}, \xi^+)}{\Delta x} \right). \end{aligned} \qquad (4.38) $$

The propagation of inequality (4.33) to time $\tau^+$ with Lemma 4.1 then gives

$$ \begin{aligned} & 0 \le \overline{w}_{\Delta t}^k (t^{n^+}, \xi^+) + k^+ - w_{\Delta t}^k (\tau^+, \xi^+) - s^+ R_{\eta} (\tau^+, \xi^+, I_{\Delta t}^{k,+} (\tau^+)) \\ & -s^+ b (\xi^+, I_{\Delta t}^{k,+} (\tau^+)) H_{\eta} \left( \frac{\overline{w}_{\Delta t}^k (t^{n^+}, \xi^+) - \overline{w}_{\Delta t}^k (t^{n^+}, \xi^+ - \Delta x)}{\Delta x}, \frac{\overline{w}_{\Delta t}^k (t^{n^+}, \xi^+ + \Delta x) - \overline{w}_{\Delta t}^k (t^{n^+}, \xi^+)}{\Delta x} \right), \end{aligned} $$

that is, using (4.35) and (4.34),

$$ \begin{aligned} & 0 \le \frac{t^+ - t^{n^+} - s^+/2}{\sqrt{\Delta t}} - R_{\eta} (\tau^+, \xi^+, I_{\Delta t}^{k,+} (\tau^+)) \\ & -b (\xi^+, I_{\Delta t}^{k,+} (\tau^+)) H_{\eta} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \frac{\sqrt{\Delta x}}{2} + \alpha e^{t^+} \xi^+ - \alpha e^{t^+} \frac{\Delta x}{2}, \right. \\ & \qquad \left. -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} + \frac{\sqrt{\Delta x}}{2} + \alpha e^{t^+} \xi^+ + \alpha e^{t^+} \frac{\Delta x}{2} \right). \end{aligned} \qquad (4.39) $$
The choice of $\sigma^+(\Delta t, k)$, and the conclusion of the argument are a consequence of (4.32) and the previous estimate. Indeed, adding the two expressions yields

$$ \mathcal{T}_1 + \mathcal{T}_2 + \mathcal{T}_3 \ge 0, \quad (4.40) $$

with

$$
\begin{align*}
\mathcal{T}_1 &= \frac{s^+}{2\sqrt{\Delta t}} - \alpha \frac{e^{t^+}}{2} (x^{+2} + \xi^{+2}) - \frac{\alpha}{(T-t^+)^2} - \sigma - C_{H_\eta}^2 \alpha e^T \\
\mathcal{T}_2 &= \mathcal{R} (t^+, x^+, I_0^{k,+} (t^+)) - R_\eta (\tau^+, \xi^+, I_{\Delta t}^{k,+} (\tau^+)) \\
\mathcal{T}_3 &= b (x^+, I_0^{k,+} (t^+)) \mathcal{H} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) - b (\xi^+, I_{\Delta t}^{k,+} (\tau^+)) \\
&\qquad H_\eta \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \frac{\sqrt{\Delta x}}{2} + \alpha e^{t^+} \xi^+ - \alpha e^{t^+} \frac{\Delta x}{2}, -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} + \frac{\sqrt{\Delta x}}{2} + \alpha e^{t^+} \xi^+ + \alpha e^{t^+} \frac{\Delta x}{2} \right),
\end{align*}
$$

and these three terms are then estimated independently. To begin with, as $s^+ \in (0, \Delta t]$, one has

$$ \mathcal{T}_1 \le \sqrt{\Delta t} - \alpha \frac{e^{t^+}}{2} (x^{+2} + \xi^{+2}) - \sigma - C_{H_\eta}^2 \alpha e^T. \quad (4.41) $$

The second one is a consequence of the K-Lipschitz regularity of $\mathcal{R}$ and $R_\eta$ in (A6)-(A8)-(A9), of the k-Lipschitz regularity of $I_0^{k,+}$ and $I_{\Delta t}^{k,+}$ (Lemma 4.4), and of the fact that $R_\eta$ is supposed to be a good approximation of $\mathcal{R}$,

$$ \mathcal{T}_2 \le K |t^+ - \tau^+| + K |x^+ - \xi^+| + K \left\| I_0^{k,+} - I_{\Delta t}^{k,+} \right\|_{\infty} + Kk |\tau^+ - t^+| + \|\mathcal{R} - R_\eta\|_{\infty}, $$

and thanks to (4.30), and to (A17)-(2.7), we obtain

$$ \mathcal{T}_2 \le C_2^k \left( \sqrt{\Delta t} + \sqrt{\Delta x} + \left\| I_0^{k,+} - I_{\Delta t}^{k,+} \right\|_{\infty} \right), \quad (4.42) $$

where $C_2^k$ depends only on $K$, $Lip_t$, $Lip_x(T)$ and on $k$. The estimate of $\mathcal{T}_3$ is obtained similarly, using the Lipschitz property of $b$ in (A4), the fact that $H_\eta(p,p)$ is an approximation of $\mathcal{H}(p)$ in (A16), and the $C_{H_\eta}$-Lipschitz property of $H_\eta$. Note that the latter is true because (4.36) holds. More precisely, $\mathcal{T}_3$ is reformulated as

$$ \mathcal{T}_3 = \mathcal{T}_3^a + \mathcal{T}_3^b + \mathcal{T}_3^c + \mathcal{T}_3^d, $$

with

$$
\begin{align*}
\mathcal{T}_3^a &= b \left( \xi^+, I_{\Delta t}^{k,+} (\tau^+) \right) H_\eta \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+, -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) \\
&\quad - b \left( \xi^+, I_{\Delta t}^{k,+} (\tau^+) \right) H_\eta \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \frac{\sqrt{\Delta x}}{2} + \alpha e^{t^+} \xi^+ - \alpha e^{t^+} \frac{\Delta x}{2}, -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} + \frac{\sqrt{\Delta x}}{2} + \alpha e^{t^+} \xi^+ + \alpha e^{t^+} \frac{\Delta x}{2} \right) \\
\mathcal{T}_3^b &= -b \left( \xi^+, I_{\Delta t}^{k,+} (\tau^+) \right) H_\eta \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+, -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) \\
&\quad + b \left( \xi^+, I_{\Delta t}^{k,+} (\tau^+) \right) \mathcal{H} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) \\
\mathcal{T}_3^c &= -b \left( \xi^+, I_{\Delta t}^{k,+} (\tau^+) \right) \mathcal{H} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) + b \left( x^+, I_0^{k,+} (\tau^+) \right) \mathcal{H} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) \\
\mathcal{T}_3^d &= -b \left( x^+, I_0^{k,+} (\tau^+) \right) \mathcal{H} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right) + b \left( x^+, I_0^{k,+} (t^+) \right) \mathcal{H} \left( -\frac{x^+ - \xi^+}{\sqrt{\Delta x}} - \alpha e^{t^+} x^+ \right).
\end{align*}
$$

Then, $\mathcal{T}_3^a$ is estimated with (2.1)-(2.5), $\mathcal{T}_3^b$ with (A16), $\mathcal{T}_3^c$ with the Lipschitz property of $b$ in (A4), and

$$ \sup_{|p| \le L_{\max}} |\mathcal{H}(p)| \le L_{\max} C_{H_\eta} + C_{L_{\max}}, $$

that holds true since $\mathcal{H}(0) = 0$, because $H_\eta$ enjoys $C_{H_\eta}$-Lipschitz regularity, and because it is an approximation of $\mathcal{H}$, with (A16)-(2.7). Eventually, the k-Lipschitz property of $I_0^{k,+}$ (Lemma 4.4) is used in $\mathcal{T}_3^d$, so that

$$ \mathcal{T}_3 \le C_3^k \left( \sqrt{\Delta x} + \sqrt{\Delta t} + \left\| I_{\Delta t}^{k,+} - I_0^{k,+} \right\|_{\infty} \right) + C_{H_\eta} \alpha e^{t^+} (|x^+| + |\xi^+|), \quad (4.43) $$
where $C_3^k$ depends only on the constants arising in the assumptions, namely $C_{H_\eta}$, $b_M$, $C_{L_{\max}}$, $Lip_x(T)$, $Lip_t$, and on $k$. We now come back to (4.40), and gather (4.41), (4.42) and (4.43), so that

$$ \sqrt{\Delta t} + (C_2^k + C_3^k) (\sqrt{\Delta x} + \sqrt{\Delta t} + \|I_{\Delta t}^{k,+} - I_0^{k,+}\|_{\infty}) \ge \sigma + \alpha C_{H_\eta}^2 (e^T - e^{t^+}) \\ + \alpha \frac{e^{t^+}}{2} (|x^+|^2 - 2C_{H_\eta}|x^+| + C_{H_\eta}^2) + \alpha \frac{e^{t^+}}{2} (|\xi^+|^2 - 2C_{H_\eta}|\xi^+| + C_{H_\eta}^2), $$

which finally gives the conclusion, as

$$ \sqrt{\Delta t} + (C_2^k + C_3^k) (\sqrt{\Delta x} + \sqrt{\Delta t} + \|I_{\Delta t}^{k,+} - I_0^{k,+}\|_{\infty}) \ge \sigma. \quad (4.44) $$

Indeed, as $\Delta t$ and $\Delta x$ depend on each other with (2.6), one can choose $\sigma = \sigma^+(\Delta t, k)$, with $\sigma^+(\Delta t, k) \to_{\Delta t \to 0} 0$ when $k$ is fixed, such that the inequality (4.44) is a contradiction. With this choice $\sigma^+(\Delta t, k)$, the assumption $\tau^+ > 0$ cannot hold, which yields the conclusion $\tau^+ = 0$. The estimate (4.27) is then a straightforward consequence of (4.30), and the inequality (4.28) comes from the expression of $\Psi_{\alpha,\sigma}^+$ in (4.26). Indeed, since $\tau^+ = 0$, and $w^k(0, \cdot) = w_{\Delta t}^k(0, \cdot)$, the following estimate holds true

$$ \Psi_{\alpha, \sigma}^{+} (t^{+}, x^{+}, \tau^{+}, \xi^{+}) \ge w^{k} (t^{+}, x^{+}) - w^{k} (0, x^{+}) + w_{\Delta t}^{k} (0, x^{+}) - w_{\Delta t}^{k} (0, \xi^{+}). $$

The Lipschitz regularity of $w^k$ and $w_{\Delta t}^k$ stated in Lemma 4.5, and (4.30) finally gives (4.28). $\square$

We may now show the convergence of the solution of scheme (S) towards the viscosity solution of (HJ). It is a consequence of

$$ v^{\infty} \le u_0 \le w^{\infty}, \quad (4.45) $$

where $v^\infty$, and $w^\infty$ are defined in (4.23)-(4.24) and $u_0$ in Lemma 4.3. As in the proof of Lemma 4.6, let us detail only the upper bound here, as the bound from below is very similar. For all $(t, x) \in [0, T] \times \mathbb{R}$, one has, using (4.28), and the fact that $(t^+, x^+, \tau^+, \xi^+)$ is a minimum point of $\Psi_{\alpha, \sigma}^+$,

$$ -2\mathrm{Lip}_t^2 \sqrt{\Delta t} - 10\mathrm{Lip}_x(T)^2 \sqrt{\Delta x} \le \Psi_{\alpha, \sigma}^+ (t^+, x^+, \tau^+, \xi^+) \le \Psi_{\alpha, \sigma}^+ (t, x, t, x), $$

that is,

$$ -2\mathrm{Lip}_t^2 \sqrt{\Delta t} - 10\mathrm{Lip}_x(T)^2 \sqrt{\Delta x} \le w^k(t, x) - w_{\Delta t}^k(t, x) + \alpha e^t x^2 + \frac{\alpha}{T-t} + (\sigma^+(\Delta t, k) + C_{H_\eta}^2 \alpha e^T)t. $$

Let now $\alpha \to 0$ in previous inequality, it gives

$$ u_{\Delta t}(t, x) - 2\mathrm{Lip}_t^2 \sqrt{\Delta t} - 10\mathrm{Lip}_x(T)^2 \sqrt{\Delta x} - \sigma^+(\Delta t, k) t \le w^k(t, x), $$

where we also used the estimate $u_{\Delta t} \le w_{\Delta t}^k$ stated in Lemma 4.5. Let now $\Delta t \to 0$, with $k > 0$ fixed. By definition of $\Delta x$, $\sigma^+(\Delta t, k)$ and $u_{\Delta t}$ the limit of the right-hand side is $u_0(t, x)$. Therefore, for all $k > 0$, $t \in [0, T]$, and $x \in \mathbb{R}$,

$$ u_0(t, x) \le w^k(t, x). $$

One can finally let $k \to \infty$ in the previous inequality, so that for all $t \in [0, T]$ and for all $x \in \mathbb{R}$, $u_0(t, x) \le w^\infty(t, x)$, which is the second estimate in (4.45).

To conclude, recall that $v^\infty = w^\infty$, thanks to Remark 4.7, and that $u_0$, $v^\infty$ and $w^\infty$ enjoy Lipschitz regularity. Then, inequalities in (4.45) are in fact equalities, and $u_0$ is the viscosity solution of (4.23) (or (4.24), equivalently). In addition, for all $t \in [0, T]$, $\min u_0(t, \cdot) = 0$, as proved in Lemma 4.3, so that Theorem 1.1 yields that $(u_0, I_0)$ is the viscosity solution of (HJ). Eventually, as Theorem 1.1 also states the uniqueness of the viscosity solution of (HJ), the restriction *up to a subsequence* of Remark 4.6, can then be removed, so that Prop. 2.4 holds.

## 4.3 Relaxation of the CFL condition

The convergence of scheme (S) is proved under strong assumptions on the relative size of the discretization parameters $\Delta t$ and $\Delta x$. Indeed, the CFL condition (2.6) must be satisfied, but the definition of $L_{\max}$ in (2.5) includes a security margin required by the proof of the convergence. It makes the CFL very restrictive, and refining time meshes, so that (2.6) holds, leads to very costly numerical computations. This is especially accurate when the Hamiltonian $\mathcal{H}$ contains an exponential, as in the example (3.1).

In practice, considering the CFL with the larger slope that appears in the computations, instead of $L_{\max}$, is enough to ensure the stability of the computations. One can even prove that the scheme converges, using this relaxed CFL condition. The key idea is to linearize $\mathcal{H}$ and $H_\eta$ for slopes stronger than what appears in the simulations. The process is the following:
(i). Make an educated guess of the steepest slope $S$ that could appear. In the worst-case scenario, $\text{Lip}_x(T)$ suits, but smaller values are targeted.

(ii). Let, for all $p \in \mathbb{R}$:

$$ \tilde{\mathcal{H}}(p) = \begin{cases} \mathcal{H}(S) + (p-S)\mathcal{H}'(S) & \text{if } p > S \\ \mathcal{H}(p) & \text{if } -S \le p \le S \\ \mathcal{H}(-S) + (p-S)\mathcal{H}'(-S) & \text{if } p < -S. \end{cases} $$

Note that if $S$ is greater than the stronger slope of the solution $u$ of (HJ), the values outside of $[-S, S]$ are never reached by $p = \nabla_x u$. Hence, replacing $\mathcal{H}$ by $\tilde{\mathcal{H}}$ in (HJ) does not change the equation. Then, linearize $H_\eta$ for slopes stronger than $S$. Let for all $(p, q) \in \mathbb{R}^2$:

$$ \tilde{H}_{\eta}(p, q) = \begin{cases} H_{\eta}(p, q) & \text{if } -S \le p \le S, -S \le q \le S \\ H_{\eta}(p, S) + (q - S)\mathcal{H}'(S) & \text{if } -S \le p \le S, q > S \\ H_{\eta}(p, -S) + (q + S)\mathcal{H}'(-S) & \text{if } -S \le p \le S, q < -S \\ H_{\eta}(S, q) + (p - S)\mathcal{H}'(S) & \text{if } p > S, -S \le q \le S \\ H_{\eta}(-S, q) + (p + S)\mathcal{H}'(-S) & \text{if } p < -S, -S \le q \le S \\ H_{\eta}(S, S) + (p - S)\mathcal{H}'(S) + (q - S)\mathcal{H}'(S) & \text{if } p > S, q > S \\ H_{\eta}(S, -S) + (p - S)\mathcal{H}'(S) + (q + S)\mathcal{H}'(-S) & \text{if } p > S, q < -S \\ H_{\eta}(-S, S) + (p + S)\mathcal{H}'(-S) + (q - S)\mathcal{H}'(S) & \text{if } p < -S, q > S \\ H_{\eta}(-S, -S) + (p + S)\mathcal{H}'(-S) + (q + S)\mathcal{H}'(-S) & \text{if } p < -S, q < -S. \end{cases} $$

Note that, if $S$ is large enough, replacing $H_\eta$ by $\tilde{H}_\eta$ does not change the scheme. However, in any case, $\tilde{H}_\eta$ satisfies the assumptions required for $H_\eta$. As the slope of $\tilde{H}_\eta$ does not increase for $|p| \ge S$ or $|q| \ge S$, no security margin is to be added in the proof of Section 4.2.3.

(iii). Choose $\Delta x$ according the CFL associated to $S$, for instance

$$ \Delta t \text{Lip}_{H_{\eta}}(S) = \frac{\Delta x}{2}, \qquad (4.46) $$

and $\eta = \sqrt{\Delta x}$ as in (2.7).

(iv). Compute the solutions of the new system, and check that for all $i$ and all $n$, $|u_i^n - u_{i-1}^n| \le S\Delta x$. If this does not hold, try again from step (i) with a more educated guess. For instance, if the first time with larger slopes is $0 < t < T$, one could try $S' = \text{Lip}_{\text{in}} + (S - \text{Lip}_{\text{in}})T/t$.

# 5 Numerical tests and discussion

In this section, we present some numerical tests to highlight the properties of schemes $(S)$ and $(S_\epsilon)$. Unless other choices are specified, the following birth rate, and net growth rate parameters will be used

$$ b(x, I) = \frac{8 - (x + 1)^2}{1 + (x + 1)^2} \frac{1}{2 + I}, \qquad (5.1) $$

$$ \mathcal{R}(t, x, I) = - \frac{(x+3)^2}{1+(x+3)^2} - \frac{x^2}{1+x^2}(t+1)I. \qquad (5.2) $$

Two initial data $u^{\text{in}}$ are considered, depending on the tests. The first one will be denoted $u^{\text{in}}_{\text{conv}}$. It is convex, and defined for all $x \in \mathbb{R}$ by,

$$ u_{\text{conv}}^{\text{in}} = \sqrt{1 + x^2} - 1. \qquad (5.3) $$

The second one has two local minima, is hence not convex, and will be denoted $u^{\text{in}}_{\text{not conv}}$. It is given for all $x \in \mathbb{R}$ by,

$$ u_{\text{not conv}}^{\text{in}}(x) = \left( \left(1 + (x+3)^2\right)^{1/4} - 1 \right) \left( \left(2 + (x-2)^2\right)^{1/4} - 1 \right). \qquad (5.4) $$

In all tests, the $x$ domain is $[-10, 10]$, but the final time $T$ varies and will always be specified. Regarding the discretization, $\Delta t$ is defined specifically for all tests, and $\Delta x$ is determined as a function of $\Delta t$ using the relation (4.46). The value of $S$ will always be indicated, and an order of magnitude of $\Delta x$ will be given for
information. Among the list of discretizations of $\mathcal{H}$ proposed in Section 3.1, we focus on $(H_{\eta}^{\text{P1}})$ and $(H_{\eta}^{\text{CCS}})$. We indeed favored, $(H_{\eta}^{\text{P1}})$ as it naturally appears as the limit of the AP scheme $(S_{\epsilon})$. It is in the convex setting but not in the flat setting, and $(H_{\eta}^{\text{CCS}})$ is its modification to a scheme in the flat setting. The choice will be made clear for each test.

## 5.1 Properties of $(S)$

We investigate the properties of scheme $(S)$, starting with its convergence. It is established in Proposition 2.4, but the compactness argument that is used in the proof does not yield any convergence rate. We test it numerically in what follows. Contrary to the particular case of the quadratic Hamiltonian $\mathcal{H}(p) = p^2$ that has been tested in [7], no analytical solution of (HJ) is available when $\mathcal{H}$ is defined by (3.1). We hence define a reference solution for scheme $(S)$, using the same scheme with a refined grid. Let $T = 0.5$, $S = 2$, and fix $b$, $\mathcal{R}$ and a convex initial data as in (5.1), (5.2) and (5.3). The reference solution $(u_{\text{ref}}, I_{\text{ref}})$ is computed with $\Delta t = 4 \cdot 10^{-5}$ (which yields $\Delta x \simeq 2.4 \cdot 10^{-3}$), and it is compared to solutions $(u_{\Delta t}, I_{\Delta t})$ computed for a sequence of $\Delta t$ varying from $10^{-2}$ to $10^{-4}$. The corresponding $\Delta x$ hence vary approximately from 0.6 to $6 \cdot 10^{-3}$. Both $(H_{\eta}^{\text{P1}})$ and $(H_{\eta}^{\text{CCS}})$ are tested, and the convergence errors

$$ \| u_{\text{ref}}^{N_{\text{ref}}^{\text{ref}}} - u_{\Delta t}^{N_{\Delta t}^{\Delta t}} \|_{\infty}, \quad \text{and} \quad \| I_{\text{ref}} - I_{\Delta t} \|_{L^1(0,T)}, \qquad (5.5) $$

are displayed as a function of $\Delta t$ in logarithmic scale in Figure 1. These errors are defined using suitable norms regarding the expected regularity of the solution of (HJ), see [7] for details. This suggests that the numerical order of scheme $S$ is 1, both for $(H_{\eta}^{\text{P1}})$ and $(H_{\eta}^{\text{CCS}})$.

<figure><img src="image_4.png" /><figcaption>Figure 1: Scheme $(S)$ with $T = 0.5$, $\Delta t_{\text{ref}} = 4 \cdot 10^{-5}$, $S = 2$, $b$, $\mathcal{R}$ and $u_{\text{conv}}^{\text{in}}$ defined in (5.1), (5.1) and (5.3). Convergence test for the solution of $(S)$ with $(H_{\eta}^{\text{P1}})$ and $(H_{\eta}^{\text{CCS}})$. Left: Error (5.5) on the component $u$, as a function of $\Delta t$ (log scale). Right: Error (5.5) on the component $I$, as a function of $\Delta t$ (log scale).</figcaption></figure>

Focus now on the qualitative properties of scheme $(S)$, and especially on the properties of $I$, or equivalently of $I_{\Delta t}$ defined by (2.20). When the dependency in $t$ of $\mathcal{R}$ is Lipschitz, (2.15) yields that the constraint $I$ in (HJ) can decrease, but only continuously. It is indeed the limit of $I_{\Delta t}$ when $\Delta t \to 0$, thanks to Proposition 2.4, and its finite differences are bounded from below. Roughly speaking, this means that $I$ has only increasing jumps. This property is illustrated on the left-hand side of Fig. 2, where $I_{\Delta t}$ is displayed for as a function of $t \in [0, 1/2]$. One can notice that $I_{\Delta t}$ starts with a smooth decay, before an increasing jump. The right-hand side of Fig. 2 presents $u$, or equivalently $u_{\Delta t}$ defined by (2.22), as a function of $t$ and $x$. Remark that the location of the minimum of $u_{\Delta t}$ jumps when $I_{\Delta t}$ jumps. These figures have been obtained with discretization $(H_{\eta}^{\text{P1}})$, parameters $b$ and $\mathcal{R}$ defined by (5.1) and (5.2), the initial data $u_{\text{not conv}}^{\text{in}}$ as in (5.4), $\Delta t = 2 \cdot 10^{-5}$ (hence $\Delta x \simeq 1.2 \cdot 10^{-3}$).

When $\mathcal{R}$ does not depend on $t$, it has been proved in [8], that $I$ is non-decreasing. It is worth noticing that this property is numerically preserved only in the flat setting. Fig. 3 highlights this behavior. We considered once again $b$ defined in (5.1) and the initial data $u_{\text{not conv}}^{\text{in}}$ (5.4). We used $\mathcal{R}(0, x, I)$ instead of $\mathcal{R}$ defined in (5.2), to make it independent of $t$. We represented $I_{\Delta t}$ as a function of $t$, when it is computed with a scheme in the flat setting as $(H_{\eta}^{\text{CCS}})$, or with a scheme in the convex setting as $(H_{\eta}^{\text{P1}})$. The left-hand side of Fig. 3 is computed with a very coarse grid, $\Delta t = 10^{-2}$ (that is, $\Delta x \simeq 0.6$). The numeric-induced decay of $I_{\Delta t}$ in the convex setting $(H_{\eta}^{\text{P1}})$ is clearly visible, while $I_{\Delta t}$ is non-decreasing in the flat setting $(H_{\eta}^{\text{CCS}})$. Of course,
<figure><img src="image_5.png" /><figcaption>Figure 2: Scheme (S) with $T = 0.5$, $\Delta t = 2 \cdot 10^{-5}$, $S = 2$, $b$, $\mathcal{R}$ and $u_{\text{not conv}}^{\text{in}}$ defined in (5.1), (5.2) and (5.4). Qualitative properties of (S) with $(H_{\eta}^{\text{P1}})$. Left: $I_{\Delta t}$ as a function of $t$. Right: $u_{\Delta t}$ as a function of $x$ (trait) and $t$ (time). The solid lines represent respectively the local minima of $u_{\Delta t}$, and the optimal trait. Light colors emphasize small values of $u_{\Delta t}$.</figcaption></figure>

this phenomenon is only due to numerical approximations, and the decay tends to 0 when $\Delta t \to 0$, as the limit $I$ of $I_{\Delta t}$ is non-decreasing. The right-hand side of Fig. 3 shows that the artificial decay of $(H_{\eta}^{\text{P1}})$ is much smaller when $\Delta t = 10^{-4}$ (and $\Delta x \simeq 6 \cdot 10^{-3}$).

<figure><img src="image_6.png" /><figcaption>Figure 3: Scheme (S) with $T = 1.5$, $S = 2$, $\mathcal{R}(0, x, I)$ defined in (5.2), and $b$ and $u_{\text{not conv}}^{\text{in}}$ defined in (5.1) and (5.4). $I_{\Delta t}$ as a function of $t$. Left: with $\Delta t = 10^{-2}$. Right: with $\Delta t = 10^{-4}$.</figcaption></figure>

## 5.2 Asymptotic-preserving property of scheme ($S_{\varepsilon}$)

Consider now the scheme ($S_{\varepsilon}$). As it is stated in Proposition 3.4, it enjoys the AP property. Indeed, its solution converges to the solution of scheme (S) when $\varepsilon$ goes to 0 with fixed discretization. The discretization being fixed, denote $(v^{\varepsilon}, J^{\varepsilon})$ the solution of scheme ($S_{\varepsilon}$) for a given $\varepsilon > 0$, $(u, I)$ the solution of scheme (S). Thanks to the AP property, $(v^{\varepsilon}, J^{\varepsilon})$ go to $(u, I)$ when $\varepsilon \to 0$. This property is tested in Fig. 4, where the errors

$$ \|v^{\varepsilon} - u\|_{\infty}, \quad \|J^{\varepsilon} - I\|_{L^{1}(0,T)}, \qquad (5.6) $$

are displayed as functions of $\varepsilon$, in logarithmic scale. Here again, the functional space have been chosen according to the expected regularity of the solutions, see [7]. As expected, these errors decrease with $\varepsilon$, a saturation for the small errors in the component $J_{\varepsilon}$ excepted. Moreover, according to this test the numerical order of convergence of the solution of ($P_{\varepsilon}$) to the solution of (HJ) with (3.1) is 1. The schemes (S) and ($S_{\varepsilon}$)
are run with $b, \mathcal{R}$ and $u_{\text{not conv}}^{\text{in}}$ defined in (5.1), (5.2) and (5.4), with $T = 0.5$, $S = 2$ and $\Delta t = 10^{-3}$ (hence, $\Delta x \simeq 6 \cdot 10^{-2}$). Scheme $(\mathcal{S})$ is implemented with discretization $(H_{\eta}^{\text{P1}})$.

<figure><img src="image_7.png" /><figcaption>Figure 4: With $T = 0.5$, $\Delta t = 10^{-3}$, $S = 2$, $b, \mathcal{R}$ and $u_{\text{not conv}}^{\text{in}}$ defined in (5.1), (5.2) and (5.4). AP property of scheme $(\mathcal{S}_{\varepsilon})$. Scheme $(\mathcal{S})$ with $(H_{\eta}^{\text{P1}})$. Left: Error (5.6) on the component $v$ as a function of $\varepsilon$ (log scale). Right: Error (5.6) on the component $J$ as a function of $\varepsilon$ (log scale).</figcaption></figure>

A stronger property of AP schemes is the fact that they can be Uniformly Accurate (UA). Indeed, in the numerical test above, the discretization was fixed, because AP scheme may present an order degeneracy for some values of $\varepsilon$, typically when their order of magnitude is comparable to the discretization parameters. On the contrary, the UA property certifies that the precision is independent of the scaling parameter $\varepsilon$. The proof of this property is far beyond the scope of this paper, as the estimated order of convergence of the limit scheme $(\mathcal{S})$ is not even theoretically determined. It can however be tested numerically. As for convergence numerical tests, a reference solution must be computed to perform the test. It is done with $b, \mathcal{R}$, and $u_{\text{not conv}}^{\text{in}}$ defined in (5.1), (5.2) and (5.4), and $T = 0.5$, $S = 2$. For all $\varepsilon > 0$, the reference solution is $(v_{\Delta t_{\text{ref}}}^{\varepsilon}, J_{\Delta t_{\text{ref}}}^{\varepsilon})$, computed with $(\mathcal{S}_{\varepsilon})$, and a refined grid $\Delta t_{\text{ref}} = 2 \cdot 10^{-5}$ (that is $\Delta x_{\text{ref}} \simeq 1.2 \cdot 10^{-3}$). Then, for $\varepsilon$ varying from $10^{-16}$ to 1, and for $\Delta t$ varying from $10^{-4}$ to $10^{-2}$ ($\Delta x$ from $6 \cdot 10^{-3}$ to 0.6), compute the $(\Delta t, \varepsilon)$-dependent errors,

$$ \|v_{\Delta t}^{\varepsilon} - v_{\Delta t_{\text{ref}}}^{\varepsilon}\|_{\infty}, \|J_{\Delta t}^{\varepsilon} - J_{\Delta t_{\text{ref}}}^{\varepsilon}\|_{L^1(0,T)}. \qquad (5.7) $$

They are displayed in Fig. 5, each value of $\Delta t$ being represented in logarithmic scale as a function of $\varepsilon$. The component in $v$ of the error is on the left-hand side of the figure, and the component in $J$ is on the right-hand side. Remark that the supremum in $\varepsilon$ of these errors decreases with $\Delta t$. It is in fact of the order of magnitude of $\Delta x$. This suggests that the scheme $(\mathcal{S}_{\varepsilon})$ enjoys the UA property, with uniform order of convergence 1.

<figure><img src="image_8.png" /><figcaption>Figure 5: With $T = 0.5$, $\Delta t = 10^{-3}$, $S = 2$, $b, \mathcal{R}$ and $u_{\text{not conv}}^{\text{in}}$ defined in (5.1), (5.2) and (5.4). UA property of scheme $(\mathcal{S}_{\varepsilon})$. Left: For a range of $\Delta t$, error (5.7) on the component $v$ as a function of $\varepsilon$ (log scale). Right: For a range of $\Delta t$, error (5.7) on the component $J$ as a function of $\varepsilon$ (log scale).</figcaption></figure>
## 5.3 Scheme ($S_\varepsilon$) when the population vanishes

Determining theoretically the asymptotic behaviour of ($P_\varepsilon$), for given $b$ and $\mathcal{R}$ not necessarily satisfying nice hypothesis, is -up to our knowledge- an open question. It can however be conjectured using ($S_\varepsilon$), as it is illustrated in this section. In this section, the following birth rate will be considered,

$$b(x, I) = 2, \qquad (5.8)$$

with the initial data $u_{\text{conv}}^{\text{in}}$ defined in (5.3). All the test cases are run with $T = 6$, $\Delta t = 2 \cdot 10^{-5}$ (that is, $\Delta x \approx 2.2 \cdot 10^{-2}$), and $S = 3$. We focus on changing environments, as it is done in [18, 19, 10] for the parabolic case.

Consider first a periodic varying environment, with

$$R(t, x, I) = -(x - 3 \sin(\alpha t))^2 - I. \qquad (5.9)$$

Note that $\mathcal{R}$ does not satisfy hypothesis (A7). However, with $\alpha = 0.5$, that is a slowly-varying environment, scheme ($S_\varepsilon$) suggests that the population does not extincts, and that the asymptotic behavior of ($P_\varepsilon$) is (HJ), see Fig. 6. The left-hand side of Fig. 6 displays $I_{\Delta t}$ (or, equivalently, $I$) computed with ($S$)-($H_{\eta}^{\text{P1}}$), and $J$ computed with ($S_\varepsilon$) and $\varepsilon = 10^{-8}$. Note that both coincide, and that they are bounded from above and below by positive constants, meaning that the total size of the population does not vanish nor grow uncontrolled. The function $u_{\Delta t}$ (or the sequence $u$) defined by ($S$)-($H_{\eta}^{\text{P1}}$) is on the right-hand side of Fig. 6 ($v$ defined by ($S_\varepsilon$) is not represented, as it cannot be distinguished from the one for $u_{\Delta t}$). The optimal trait and the minimum of $u_{\Delta t}$ are represented by solid lines on the figure. Note that the local minimum of $u_{\Delta t}$ follows the optimal trait with a delay, as expected thanks to [19]. Fig. 6 has been obtained with $T = 6$, $\Delta t = 2 \cdot 10^{-5}$ and $S = 3$.

<figure><img src="image_9.png" /><figcaption>Figure 6: With $T = 6$, $\Delta t = 2 \cdot 10^{-5}$, $S = 3$, $b$, $u_{\text{conv}}^{\text{in}}$ defined in (5.8), (5.3), $\mathcal{R}$ defined in (5.9) and $\alpha = 0.5$. Left: $I_{\Delta t}$ computed with ($S$) and ($H_{\eta}^{\text{P1}}$), and $J$ computed with ($S_\varepsilon$), as functions of $t$. Right: $u_{\Delta t}$ computed with ($S$)-($H_{\eta}^{\text{P1}}$) as a function of $x$ (trait) and $t$ (time). The solid lines represent respectively the local minima of $u_{\Delta t}$, and the optimal trait. Light colors emphasize small values of $u_{\Delta t}$.</figcaption></figure>

The parameter $\alpha$ in (5.9) drives the speed at which the environment changes. When it varies too quickly, the population does not have enough time to adapt and may vanish, see [18] for periodic-varying environments, or [19] for shifted environments. This behavior is illustrated in Fig. 7, where we took $\alpha = 3$ in (5.9), and where ($P_\varepsilon$) is approximated with ($S_\varepsilon$) for $\varepsilon = 10^{-8}$. With this value of $\alpha$, scheme ($S_\varepsilon$) suggests that the population vanishes in the asymptotics $\varepsilon \to 0$ of ($P_\varepsilon$). Indeed, $J$ defined by ($S_\varepsilon$) goes to 0 rather quickly. It is displayed on the left-hand side of Fig.7 as a function of $t$. The value of $v$ is on the right-hand side of Fig. 7, as a function of $t$ and $x$. As previously, note that the minimum of $v$ follows the optimal trait with a delay. The fact that the minimum of $v$ is not close to 0 after some time, informs on the extinction of the population.

Let us emphasize on the fact that the solutions of ($S_\varepsilon$) and ($S$) do not coincide when the population vanishes. Indeed, still considering $\mathcal{R}$ as in (5.9) with $\alpha = 3$, the solution of ($S$)-($H_{\eta}^{\text{P1}}$) is displayed on Fig. 8, with the plot of $u_{\Delta t}$ as a function of $t$ and $x$ on the right-hand side. Observe on the left-hand side that $I_{\Delta t}$ takes negative values. This is due to the fact that the constraint $\min u_{\Delta t} = 0$ is always satisfied in scheme ($S$), possibly at the cost of the positivity of $I_{\Delta t}$. This suggests that, in this case of fast oscillations of the environment, the $\varepsilon \to 0$ limit of ($P_\varepsilon$) is not given by (HJ), and that the asymptotic behavior of ($P_\varepsilon$) is better conjectured with ($S_\varepsilon$) than with ($S$).
<figure><img src="image_10.png" /><figcaption>Figure 7: With $T = 6$, $\Delta t = 2 \cdot 10^{-5}$, $S = 3$, $b$, $u_{\text{conv}}^{\text{in}}$ defined in (5.8), (5.3), $\mathcal{R}$ defined in (5.9) and $\alpha = 3$. Left: $J$ computed with $(S_{\varepsilon})$ as a function of $t$. Right: $v$ computed with $(S_{\varepsilon})$ as a function of $x$ (trait) and $t$ (time). The solid lines represent respectively the local minima of $v$, and the optimal trait. Light colors emphasize small values of $v$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: With $T = 6$, $\Delta t = 2 \cdot 10^{-5}$, $S = 3$, $b$, $u_{\text{conv}}^{\text{in}}$ defined in (5.8), (5.3), $\mathcal{R}$ defined in (5.9) and $\alpha = 3$. Left: $I_{\Delta t}$ computed with $(S)$ and $(H_{\eta}^{\text{P1}})$ as a function of $t$. Right: $u_{\Delta t}$ computed with $(S)-(H_{\eta}^{\text{P1}})$ as a function of $x$ (trait) and $t$ (time). The solid lines represent respectively the local minima of $u_{\Delta t}$, and the optimal trait. Light colors emphasize small values of $u_{\Delta t}$.</figcaption></figure>

Eventually, we consider a constant by parts environment, as in [10], in a regime of quick changes in the environment. In what follows, we will define $\mathcal{R}$ as

$$R(t, x, I) = -(x - 6 \lceil \sin(3t - 0.9) \rceil) \mathbb{1}_{t>0.3})^2 - I. \quad (5.10)$$

Following the previous test case, we conjecture the small $\varepsilon$ limit of $(P_{\varepsilon})$ to vanish, because of the too fast changes of the environment. Fig. 9 suggests that scheme $(S_{\varepsilon})$ is robust enough to deal with the lack of regularity of $\mathcal{R}$. Indeed, $J$ is represented as a function of $t$ on the left-hand side of Fig. 9. Once again, it vanishes quickly, as expected. The evolution of $v$ is displayed on the right-hand side of Fig. 9 as a function of $t$ and $x$. As previously, the minimum of $v$ is nonzero after some time, and its location follows the optimal trait with a delay, despite its lack of regularity. Fig. 9 has been obtained with $\varepsilon = 10^{-8}$ in $(S_{\varepsilon})$.

Note that in this test case, scheme $(S)$ again gives results far from the ones of $(S_{\varepsilon})$. They are displayed in Fig. 10, with $I_{\Delta t}$ as a function of $t$ on the left-hand side and $u_{\Delta t}$ as a function of $t$ and $x$ on the right-hand side. Once again, the constraint $\min u_{\Delta t} = 0$ leads to negative values for $I_{\Delta t}$. Moreover, observe that $I_{\Delta t}$ has decreasing jumps. This could not have been possible if $\mathcal{R}$ had Lipschitz regularity in $t$, see Prop. 2.3-(v). Up to our knowledge, the well-posedness of (HJ) is this context is an open question. Despite the robustness of $(S)-(H_{\eta}^{\text{P1}})$ regarding the lack of regularity of $\mathcal{R}$, it corroborates the fact that (HJ) is not the $\varepsilon \to 0$ limit of $(P_{\varepsilon})$ in this case.
<figure><img src="image_12.png" /><figcaption>Figure 9: With $T = 6$, $\Delta t = 2 \cdot 10^{-5}$, $S = 3$, $b$, $u_{\text{conv}}^{\text{in}}$ and $\mathcal{R}$ defined in (5.8), (5.3) and in (5.10). Left: $\mathbf{J}$ computed with $(S_{\varepsilon})$ as a function of $t$. Right: $\mathbf{v}$ computed with $(S_{\varepsilon})$ as a function of $x$ (trait) and $t$ (time). The solid lines represent respectively the local minima of $\mathbf{v}$, and the optimal trait. Light colors emphasize small values of $\mathbf{v}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: With $T = 6$, $\Delta t = 2 \cdot 10^{-5}$, $S = 3$, $b$, $u_{\text{conv}}^{\text{in}}$ and $\mathcal{R}$ defined in (5.8), (5.3) and in (5.10). Left: $I_{\Delta t}$ computed with $(S)-(H_{\eta}^{\text{P1}})$ as a function of $t$. Right: $u_{\Delta t}$ computed with $(S)-(H_{\eta}^{\text{P1}})$ as a function of $x$ (trait) and $t$ (time). The solid lines represent respectively the local minima of $u_{\Delta t}$, and the optimal trait. Light colors emphasize small values of $u_{\Delta t}$.</figcaption></figure>

## 5.4 Generalization

In this section, we investigate the generalization of the results of the paper, especially when the dimension is greater than 1. Let us start with Lotka-Volterra integral equation in the multi-dimensional case, with Hamiltonian (3.1) as in Section 3. Suppose that the mutation kernel $K_d$ satisfies

$$ \forall x \in \mathbb{R}^d, K_d(x) = \prod_{j=1}^{d} K(x_j), $$

with $K$ defined in (3.1). Then, the Hamiltonian $\mathcal{H}$ can be reformulated as

$$ \forall \mathbf{p} \in \mathbb{R}^d, \mathcal{H}(\mathbf{p}) = \prod_{j=1}^{d} \int_{\mathbb{R}} K(x_j) e^{-z_j \mathbf{p}} d\mathbf{z}_j - 1, $$

so that the different choices of $H_{\eta}$ proposed in Section 3.1 in dimension 1 can be plugged in the product, to define a numerical Hamiltonian in dimension greater than 1. Remark that one has to pay attention to the constant $-1$, that is not in the multiplication, but this is quickly addressed. The AP scheme of Section 3.2 is generalized similarly, and Prop. 3.4 holds. However, the convergence of the scheme for the limit constrained Hamilton-Jacobi equation cannot be established as easily. This is the purpose of this section.
We studied above a class of numerical schemes for constrained Hamilton-Jacobi equations such as (HJ), in dimension 1. We discuss here the design of a scheme for (HJ) when the dimension is greater than 1, and when the Hamiltonian $H$ also depends on the trait variable $x$. This latter case happens for instance in [35], where the Hamiltonian is implicitly defined. Let us consider $u$, the viscosity solution of

$$ \left\{ \begin{array}{l} \partial_t u(t, x) + b(x, I(t))\mathcal{H}(x, \nabla_x u(t, x)) + \mathcal{R}(t, x, I(t)) = 0, \quad x \in \mathbb{R}^d, t \ge 0, \\ \min_{x \in \mathbb{R}^d} u(t, x) = 0, \quad t \ge 0, \end{array} \right. \qquad (5.11) $$

where the initial condition $u^{\text{in}}$ satisfies (A10)-(A11)-(A12). As previously, we suppose that $b$ satisfies (A3)-(A4)-(A5), and that $\mathcal{R}$ satisfies (A6)-(A7)-(A8)-(A9). As $\mathcal{H}$ now depends also on $x$, the hypothesis must be adapted. We suppose that for all $x \in \mathbb{R}^d$, $\mathbf{p} \mapsto \mathcal{H}(x, \mathbf{p})$ is convex, and (A2) is replaced by

$$ \forall x \in \mathbb{R}^d, \mathcal{H}(x, 0) = 0, \quad \text{and} \quad \forall (x, \mathbf{p}) \in \mathbb{R}^d \times \mathbb{R}^d, \mathcal{H}(x, \mathbf{p}) \ge 0, \qquad (5.12) $$

and we emphasize once again on the fact that the condition $\mathcal{H}(x, 0) = 0$ for all $x \in \mathbb{R}^d$ can be relaxed by considering $(x, \mathbf{p}) \mapsto \mathcal{H}(x, \mathbf{p}) - \mathcal{H}(x, 0)$, and modifying $\mathcal{R}$ accordingly.

For a sake of simplicity, we define a uniform trait grid, with the same trait step in all the directions, but this could be easily relaxed. Let $x_0 \in \mathbb{R}^d$ and $\Delta x > 0$. For all $\mathbf{i} = (i_1, \dots, i_d) \in \mathbb{Z}^d$, define then $x_i = x_0 + \Delta x i$. The time grid is kept unchanged.

As in dimension 1, introduce $H_\eta: \mathbb{R}^d \times \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ and $R_\eta: \mathbb{R} \times \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$, such that $H_\eta(x, \mathbf{p}, \mathbf{p})$ is an approximation of the Hamiltonian $\mathcal{H}(x, \mathbf{p})$ for all $x \in \mathbb{R}^d$, and $\mathbf{p} \in \mathbb{R}^d$, and that $R_\eta$ approximates $\mathcal{R}(t, x, I)$ for all $t \ge 0$, all $x \in \mathbb{R}^d$ and all $I \in \mathbb{R}$. Assumptions (A13)-(A14)-(A15)-(A15)-(A16) and (A17) are adapted to this multi-dimensional context. Introduce $\mathcal{L} > 0$, and suppose that $H_\eta$ and $R_\eta$ satisfy the following properties:

* **Regularity in $x$:** $\forall (\mathbf{p}, \mathbf{q}) \in \mathbb{R}^d \times \mathbb{R}^d$, $x \mapsto H_\eta(x, \mathbf{p}, \mathbf{q})$ is smooth. Moreover, $H_\eta(\cdot, \mathbf{p}, \mathbf{q})$ and its first and second derivatives are bounded uniformly with respect to $(\mathbf{p}, \mathbf{q}) \in \mathbb{R}^d \times \mathbb{R}^d$ such that $\|\mathbf{p}\|_\infty \le \mathcal{L}$, and $\|\mathbf{q}\|_\infty \le \mathcal{L}$.

* **Monotony:** $\forall \|\mathbf{p}\|_\infty \le \mathcal{L}$, $\forall \|\mathbf{q}\|_\infty \le \mathcal{L}$, $\forall x \in \mathbb{R}^d$, $\forall e \in \mathbb{R}^d$ such that $e \ge 0$ (coordinate by coordinate),

$$ H_{\eta}(x, \mathbf{p}, \mathbf{q} + \mathbf{e}) \le H_{\eta}(x, \mathbf{p}, \mathbf{q}) \le H_{\eta}(x, \mathbf{p} + \mathbf{e}, \mathbf{q}), $$

* **Exact value at 0:** $\forall x \in \mathbb{R}^d$, $H_\eta(x, 0, 0) = 0$.

* **Positivity on the diagonal:** $\forall \mathbf{p} \in \mathbb{R}^d$ such that $\|\mathbf{p}\|_\infty \le \mathcal{L}$, $\forall x \in \mathbb{R}^d$, $H_\eta(x, \mathbf{p}, \mathbf{p}) \ge 0$.

* **Consistency for $H_\eta$:** $\exists C_\mathcal{L} > 0$, $\forall x \in \mathbb{R}^d$, $\forall \|\mathbf{p}\|_\infty \le \mathcal{L}$, $\forall \eta \in (0, 1]$,

$$ |H_{\eta}(x, \mathbf{p}, \mathbf{p}) - \mathcal{H}(x, \mathbf{p})| \le C_{\mathcal{L}}\eta. $$

* **Consistency for $R_\eta$:** $\forall t \ge 0$, $\forall x \in \mathbb{R}^d$, $\forall I \ge 0$, $\forall \eta \in (0, 1]$,

$$ |R_{\eta}(t, x, I) - \mathcal{R}(t, x, I)| \le K\eta, $$

for some constant $K$ as in (A17).

To adapt the Definition 2.1 of a $\mathcal{L}$-well chosen scheme, two classes of numerical Hamiltonians are considered:

* **Flat setting:** If $\|\mathbf{p}\|_\infty \le \mathcal{L}$, and $\|\mathbf{q}\|_\infty \le \mathcal{L}$ are such that $\exists \ell \in [1, d]$, $p_\ell < 0$ and $q_\ell > 0$, then for all $x \in \mathbb{R}^d$, and for all $\eta \in (0, 1]$,

$$ H_{\eta}(x, \mathbf{p}, \mathbf{q}) = H_{\eta} \left( x, \begin{pmatrix} p_1 \\ \vdots \\ p_{\ell} \\ \vdots \\ p_d \end{pmatrix}, \begin{pmatrix} q_1 \\ \vdots \\ q_{\ell} \\ \vdots \\ q_d \end{pmatrix} \right) = H_{\eta} \left( x, \begin{pmatrix} p_1 \\ \vdots \\ 0 \\ \vdots \\ p_d \end{pmatrix}, \begin{pmatrix} q_1 \\ \vdots \\ 0 \\ \vdots \\ q_d \end{pmatrix} \right), $$

* **Convex setting:** for all $x \in \mathbb{R}^d$, $H_\eta(x, \cdot, \cdot) : B_{\|\cdot\|_\infty}(0, \mathcal{L}) \times B_{\|\cdot\|_\infty}(0, \mathcal{L}) \to \mathbb{R}$ is convex.
The adaptation of Definition 2.2 is straightforward. Then, notations are introduced to write the adaptation of scheme (S) in dimension $d > 1$. For all $\ell \in [1, d]$, denote $T_{\ell}^{-} : \mathbb{R}^{Z^{d}} \to \mathbb{R}^{Z^{d}}$, and $T_{\ell}^{+} : \mathbb{R}^{Z^{d}} \to \mathbb{R}^{Z^{d}}$. They are defined such that for all $\mathbf{u} = (u_i)_{i \in \mathbb{Z}^d} \in \mathbb{R}^{Z^d}$, $T_{\ell}^{\pm}(\mathbf{u}) \in \mathbb{R}^{Z^d}$, and for all $\mathbf{i} \in \mathbb{Z}^d$,

$$T_{\ell}^{-}(\mathbf{u})_{\mathbf{i}} = u_{i_1, \dots, i_{\ell-1}, \dots, i_{\ell}}, \quad \text{and} \quad T_{\ell}^{+}(\mathbf{u})_{\mathbf{i}} = u_{i_1, \dots, i_{\ell+1}, \dots, i_{\ell}},$$

where we denoted $u_i = u_{i_1, \dots, i_{\ell}, \dots, i_{\ell}}$. Then, let $T^{\pm}: \mathbb{R}^{Z^d} \to (\mathbb{R}^{Z^d})^d$, such that,

$$\forall \mathbf{u} \in \mathbb{R}^{Z^d}, \forall \mathbf{i} \in \mathbb{Z}^d, T^{\pm}(\mathbf{u})_{\mathbf{i}} = \begin{pmatrix} T_1^{\pm}(\mathbf{u})_{\mathbf{i}} \\ \vdots \\ T_d^{\pm}(\mathbf{u})_{\mathbf{i}} \end{pmatrix},$$

and define the adaptation of scheme (S) for all $n \in [0, N_T - 1]$ and all $\mathbf{i} \in \mathbb{Z}^d$, by

$$\left\{ \begin{array}{l} \displaystyle \frac{u_{\mathbf{i}}^{n+1} - u_{\mathbf{i}}^{n}}{\Delta t} + b(x_{\mathbf{i}}, I^{n+1}) H_{\eta} \left( x_{\mathbf{i}}, \frac{I - T^{-}}{\Delta x}(\mathbf{u}^{n})_{\mathbf{i}}, \frac{T^{+} - I}{\Delta x}(\mathbf{u}^{n})_{\mathbf{i}} \right) + R_{\eta} (t^{n+1}, x_{\mathbf{i}}, I^{n+1}) = 0 \\ \min_{\mathbf{i} \in \mathbb{Z}^d} u_{\mathbf{i}}^{n} = 0, \end{array} \right. \qquad (5.13)$$

with initialization $u_0^i = u^{\text{in}}(x_i)$ for all $\mathbf{i} \in \mathbb{Z}^d$. Using the assumptions and notations above, one can follow the proofs of Prop. 2.3 and 2.4, and remark that everything can be rewritten (with a small change in (2.6)), except the semi-concavity of the numerical solution. It is in fact straightforward in the flat setting, as in dimension 1. However, in the convex setting, it seems that the adaptation of the proof of Prop. 2.3-(vi) (Section 4.1.5) in dimension greater than 1, or with $\mathcal{H}$ depending also on $x$, is more delicate. This question will be addressed in a dedicated study.

To conclude roughly, scheme (5.13) converges when $H_{\eta}$ is in the flat setting. More precisely, Prop. 2.3 and Prop. 2.4 hold, with the necessary adaptations of notations and hypothesis. In the convex setting, no such results hold. Coming back to Lotka-Volterra integral equation, with Hamiltonian 3.1, one can notice that the tensor product of non necessarily non-negative convex discrete Hamiltonian $H_{\eta}$ is not convex. Hence, the tensor product of the examples of discrete Hamiltonians of Section 3 does not give any example of multi-dimensional discrete Hamiltonian which would be in the *convex setting* but not in the *flat setting*. The investigation of the convergence of scheme (5.13) in the *convex setting* is postponed to a future work.

# A Asymptotic-preserving property of ($S_{\varepsilon}$)

## A.1 Proof of Prop. 3.3

In this section, we show that scheme ($S_{\varepsilon}$) is well-posed, and we prove stability estimates, which are uniform w.r.t. $\varepsilon$, if $\varepsilon$ is small enough. Suppose that the hypothesis of Prop. 3.3 are satisfied. The proof we propose above strongly relies on the monotony of (3.3), as for the proof of Prop. (2.3). Note also that, it is close with the proof of the stability Lemma in [7] although technical details are different. It is done by induction on the time step. Since, the items of Prop. 3.3 are true for the initial data, suppose now that $\mathbf{v}^n$ is constructed such that (i)-(ii)-(iv) hold. We show that there exists an unique pair $(\mathbf{v}^{n+1}, J^{n+1})$ satisfying ($S_{\varepsilon}$), and that it satisfies (i)-(ii)-(iii)-(iv).

The well-posedness of ($S_{\varepsilon}$) is immediate remarking that $J^{n+1}$ is solution of

$$J^{n+1} = \Delta x \sum_{\mathbf{i} \in \mathbb{Z}} \psi(x_{\mathbf{i}}) e^{-(M_{\varepsilon}^{J^{n+1}}(\mathbf{v}^n)_{\mathbf{i}} - \Delta t R(t^{n+1}, x_{\mathbf{i}}, J^{n+1})) / \varepsilon}, \qquad (A.1)$$

that the right-hand side of this equality is decreasing w.r.t $J^{n+1}$ thanks to (A6) and Lemma 3.2, and that the left-hand side in increasing. In other words, $J^{n+1}$ is the unique solution of $\varphi^{\text{AP}}(J) = 0$, with

$$\varphi^{\text{AP}}: J \mapsto J - \Delta x \sum_{\mathbf{i} \in \mathbb{Z}} \psi(x_{\mathbf{i}}) e^{-(M_{\varepsilon}^{J}(\mathbf{v}^n) - \Delta t R(t^{n+1}, x_{\mathbf{i}}, J)) / \varepsilon}. \qquad (A.2)$$

Since (A.1) has an unique solution, it can be plugged in ($S_{\varepsilon}$) so that $\mathbf{v}^{n+1}$ is also uniquely determined.

*Remark A.1.* In practice, $J^{n+1}$ is computed using Newton's method for the equation $\varphi^{\text{AP}}(J) = 0$. However, depending on the expressions of $R$ and $b$, the computation of the iterations may be tricky. Moreover, the dependency in $\varepsilon$ forbids the use of approximated Newton's method. Indeed, the problem becomes stiff when $\varepsilon \to 0$, and more computational time would be needed for small $\varepsilon$ (if the resolution does not simply break). To avoid this issue, automatic differentiation should be used.
The bound from below for $J^{n+1}$ in (iii) is obtained from

$$ \varphi^{\text{AP}}(J) \le J - \Delta x\psi(x_j)e^{-(M_{\varepsilon}^{J}(\mathbf{v}^n)_j - \Delta t R(t^{n+1}, x_j, J))/\varepsilon}, $$

which holds for all $j \in \mathbb{Z}$ and for all $J \in \mathbb{R}$. Moreover, as $b$ is positive, and thanks to (iv), the expression of $M_{\varepsilon}^J$ in (3.3) yields for $j \in \mathbb{Z}$, such that $v_j^n = \min \mathbf{v}^n$,

$$ M_{\varepsilon}^{J}(\mathbf{v}^{n})_{j} \le v_{j}^{n} \le c_{M}\varepsilon, $$

so that for all $J \le I_m/2$,

$$ \varphi^{\text{AP}}(J) \le \varphi^{\text{AP}}(I_m/2) \le \frac{I_m}{2} - \Delta x \psi_m e^{-c_M} e^{\Delta t R(t^{n+1}, x_j, I_m/2)/\varepsilon}, $$

where the fact that $\varphi^{\text{AP}}$ is increasing was used for the first inequality, and (A19) in the second. Then (A6)-(A7) yield,

$$ R(t^{n+1}, x_j, I_m/2) \ge R(t^{n+1}, x_i, I_m) + \frac{I_m}{2K} \ge \frac{I_m}{2K}, $$

and hence,

$$ \forall J \le I_m/2, \quad \varphi^{\text{AP}}(J) \le \frac{I_m}{2} - \Delta x \psi_m e^{-c_M} e^{\Delta t \frac{I_m}{2\varepsilon K}} \underset{\varepsilon \to 0}{\longrightarrow} -\infty. $$

One can conclude that there exists $\varepsilon_1 > 0$, which depends only on $I_m$, $\psi_m$, $c_M$, $K$, $\Delta t$ and $\Delta x$, such that for all $\varepsilon < \varepsilon_1$ and for all $J \le I_m/2$, $\varphi^{\text{AP}}(J) < 0$. As $\varphi^{\text{AP}}$ is increasing, it gives the bound from below for $J^{n+1}$

$$ \forall \varepsilon < \varepsilon_1, \quad J^{n+1} \ge I_m/2. \qquad (A.3) $$

*Remark A.2.* We emphasize on the fact that $\varepsilon_1$ determined above can be fixed once for all, and such that is does not depend on $n$.

The propagation of the bound from below of $\mathbf{v}^n$ in (ii) is a straightforward consequence of Lemma 3.2. Indeed, as (3.6) is satisfied, for all $i \in \mathbb{Z}$,

$$ M_{\varepsilon}^{J^{n+1}}(\mathbf{v}^n)_i \ge \underline{a}|x_i| + \underline{\beta}_{t^n} - b_M \Delta t \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) e^{\underline{a}|z_k|}, $$

and since $R$ is decreasing w.r.t. $J$, assumption (A8) gives

$$ R(t^{n+1}, x_i, J^{n+1}) \le R(t^{n+1}, x_i, I_m/2) \le K. $$

Hence, property (ii) is obtained with the choice

$$ \underline{\beta}_{t^{n+1}} = \underline{\beta}_{t^n} - b_M \Delta t \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) e^{\underline{a}|x_i|} - \Delta t K, $$

and together with (A.3), it yields the bound of $\min \mathbf{v}^{n+1}$. Indeed, we have

$$ I_m/2 \le J^{n+1} = \Delta x \sum_{i \in \mathbb{Z}} \psi(x_i) e^{-v_i^{n+1}/ep} \le \Delta x \psi_M \left( \sum_{|i| \le N} e^{-v_i^{n+1}/\varepsilon} + \sum_{|i| > N} e^{-(\underline{\beta}_T - \underline{a}|x_0| + \underline{a}|x_i - x_0|)/\varepsilon} \right), $$

for any $N \in \mathbb{N}$. As $x_i - x_0 = i\Delta x$, the above expression is simplified, such that

$$ I_m/2 \le (2N - 1)\Delta x \psi_M e^{-\min \mathbf{v}^{n+1}/\varepsilon} + \frac{2\Delta x \psi_M e^{-(\underline{\beta}_T - \underline{a}|x_0| + \underline{a}N\Delta x)/\varepsilon}}{1 - e^{-\underline{a}\Delta x/\varepsilon}}. $$

One can now fix $N$ such that $\underline{\beta}_T - \underline{a}|x_0| + \underline{a}N\Delta x > 0$, and determine $\varepsilon_2 > 0$ such that

$$ \forall \varepsilon \in (0, \varepsilon_2), \quad \frac{2\Delta x \psi_M e^{-(\underline{\beta}_T - \underline{a}|x_0| + \underline{a}N\Delta x)/\varepsilon}}{1 - e^{-\underline{a}\Delta x/\varepsilon}} \le I_m/4. $$

This gives $\min \mathbf{v}^{n+1} \le \varepsilon c_M$, with

$$ c_M = \max \left\{ -\ln \left( \frac{I_m}{4(2N-1)\Delta x \psi_M} \right), c_M^{\text{in}} \right\}. $$
Remark A.3. Note that $N$, $\varepsilon_2$ and $c_M$ only depend on the constants arising in the assumptions, and of $\Delta x$. They can hence be fixed once for all, and they are independent on $n$.

The upper bound for $J^{n+1}$ is a consequence of the monotony of $M_{\varepsilon}^J$, and of the previous result. Indeed, as (3.6) is satisfied, and since the inequality $v_i^n \ge c_m\varepsilon$, holds for all $i \in \mathbb{Z}$, it can be propagated so that

$$ \forall i \in \mathbb{Z}, \forall J \in \mathbb{R}, M_{\varepsilon}^{J}(\mathbf{v}^n)_i \ge c_m\varepsilon - \Delta t \Delta z b_M \sum_{k \in \mathbb{Z}} K(z_k), $$

where we used (A3). Consider now $\mathbf{v}^{n+1}$ at point $j$ such that $v_j^{n+1} = \min \mathbf{v}^{n+1}$. Such a minimum point exists, since $\mathbf{v}^{n+1}$ satisfies (ii). Coming back to the expression of $(S_{\varepsilon})$,

$$ \begin{aligned} \Delta t R(t^{n+1}, x_j, J^{n+1}) &= M_{\varepsilon}^{J^{n+1}}(\mathbf{v}^n)_j - v_j^{n+1} & (A.4) \\ &\ge (c_m - c_M)\varepsilon - \Delta t b_M \Delta z \sum_{k \in \mathbb{Z}} K(z_k) + \Delta t R(t^{n+1}, x_j, I_M), \end{aligned} $$

thanks to (A7). We conclude using (A6), to get another bound from below for $R(t^{n+1}, x_j, I_M)$,

$$ R(t^{n+1}, x_j, I_M) \ge R(t^{n+1}, x_j, 2I_M) + \frac{I_M}{K_1}, $$

which is plugged in (A.4), to get with (A7)-(3.5)

$$ \Delta t R(t^{n+1}, x_j, J^{n+1}) \ge \Delta t R(t^{n+1}, x_j, 2I_M) - \varepsilon(c_M - c_m) + \frac{\Delta t}{2}. $$

Then, we define $\varepsilon_3$ such that for all $\varepsilon \in (0, \varepsilon_3)$, $-\varepsilon(c_M - c_m) + \frac{\Delta t}{2} \ge 0$, and (A6) eventually yields $J^{n+1} \le 2I_M$, that is the second inequality in (iii).

Remark A.4. Once again, we emphasize on the fact that $\varepsilon_3$ depends only on the constants defined in the assumptions and on $\Delta t$. In particular, it does not depend on $n$.

The bound from below of $\min \mathbf{v}^{n+1}$ is straightforward using $J^{n+1} \le 2I_M$. Indeed, for all $i \in \mathbb{Z}$,

$$ \psi_m \Delta x e^{-v_i^{n+1}/\varepsilon} \le \Delta x \sum_{i \in \mathbb{Z}} \psi(x_i) e^{-v_i^{n+1}/\varepsilon} = J^{n+1} \le 2I_M, $$

so that (iv) holds, with $c_m$ independent of $n$. One can for instance consider

$$ c_m = \min \left\{ -\ln \left( \frac{2I_M}{\psi_m \Delta x} \right), c_m^{\text{in}} \right\}. $$

The Lipschitz bound (i) is then a consequence of (3.6), of Lemma 3.2, and of (iii). To conclude the proof, we eventually define $\varepsilon_0 = \min(\varepsilon_1, \varepsilon_2, \varepsilon_3)$, such that items (i)-(ii)-(iii)-(iv) hold.

## A.2 Proof of Prop. 3.4

As previously, the proof of Prop. 3.4 is done by induction. Thanks to the assumptions, there exists $\mathbf{u}^0$ such that $\|\mathbf{u}^0 - \mathbf{v}^0\|_{\infty} \to_{\varepsilon \to 0} 0$ and that $\min \mathbf{u}^0 = 0$. We suppose that it is true for a given $n \in [0, N_T - 1]$, and we show that there exists $\mathbf{u}^{n+1}$, and $I^{n+1}$ such that

$$ \|\mathbf{v}^{n+1} - \mathbf{u}^{n+1}\|_{\infty} \to_{\varepsilon \to 0} 0, \quad \text{and} \quad J^{n+1} \to_{\varepsilon \to 0} I^{n+1}, $$

and that $\mathbf{u}$ and $I$ satisfy scheme $(S_0)$. First of all, Prop. 3.3 gives that $(J^{n+1})_{\varepsilon \in (0, \varepsilon_0]}$ is uniformly bounded w.r.t. $\varepsilon$. Hence, there exists $I^{n+1}$ such that $J^{n+1} \to_{\varepsilon \to 0} I^{n+1}$, up to an extraction. Consider now this extraction and define $\mathbf{u}^{n+1, I^{n+1}}$ as in the first line of $(S_0)$. Let us show that

$$ \| \mathbf{v}^{n+1} - \mathbf{u}^{n+1, I^{n+1}} \|_{\infty} \underset{\varepsilon \to 0}{\longrightarrow} 0, \qquad (A.5) $$

and that $\min \mathbf{u}^{n+1, I^{n+1}} = 0$. Let $i \in \mathbb{Z}$, one has, with $(S_{\varepsilon})$ and $(S_0)$,

$$ |v_i^{n+1} - u_i^{n+1, I^{n+1}}| \le \Delta t |R(t^{n+1}, x_i, I^{n+1}) - R(t^{n+1}, x_i, J^{n+1})| + |M_0^{I^{n+1}}(\mathbf{u}^n)_i - M_0^{J^{n+1}}(\mathbf{u}^n)_i| $$
$$+ |\mathcal{M}_0^{J^{n+1}}(\mathbf{u}^n)_i - \mathcal{M}_0^{J^{n+1}}(\mathbf{v}^n)_i| + |\mathcal{M}_0^{J^{n+1}}(\mathbf{v}^n)_i - \mathcal{M}_\varepsilon^{J^{n+1}}(\mathbf{v}^n)_i|,$$

where the three first terms can be easily estimated. Indeed, using (A6), one has

$$|R(t^{n+1}, x_i, I^{n+1}) - R(t^{n+1}, x_i, J^{n+1})| \le K |I^{n+1} - J^{n+1}|,$$

while the definition of $\mathcal{M}_0^I$ in (3.7), the fact that $b$ is $b_M$-Lipschitz (see (A4)), and that $\mathbf{u}^n$ is Lip$_x(t^n)$-Lipschitz (Prop. 2.3-(i)) yield

$$|\mathcal{M}_0^{I^{n+1}}(\mathbf{u}^n)_i - \mathcal{M}_0^{I^{n+1}}(\mathbf{u}^n)_i| \le b_M |I^{n+1} - J^{n+1}| \Delta t \Delta z \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) e^{|z_k| \text{Lip}_x(t^n)}.$$

Eventually, as $\mathcal{M}_0^{J^{n+1}}$ satisfies Lemma 3.2,

$$\| \mathcal{M}_0^{J^{n+1}}(\mathbf{u}^n) - \mathcal{M}_0^{J^{n+1}}(\mathbf{v}^n) \|_{\infty} \le \| \mathbf{u}^n - \mathbf{v}^n \|_{\infty}.$$

The estimate of the last term is obtained thanks to its expression. It is rewritten to obtain

$$\frac{1}{\Delta t \Delta z} \left| \mathcal{M}_0^{J^{n+1}}(\mathbf{v}^n)_i - \mathcal{M}_\varepsilon^{J^{n+1}}(\mathbf{v}^n)_i \right| \\ \le \sum_{k \ge 0} \mathcal{K}(z_k) |b(x_i + \varepsilon z_k, J^{n+1}) - b(x_i, J^{n+1})| e^{-z_k(v_{i+1}^n - v_i^n)/\Delta x} \\ + \sum_{k \le -1} \mathcal{K}(z_k) |b(x_i + \varepsilon z_k, J^{n+1}) - b(x_i, J^{n+1})| e^{-z_k(v_i^n - v_{i-1}^n)/\Delta x} \\ + \sum_{\varepsilon z_k \ge \Delta x} \mathcal{K}(z_k) b(x_i + \varepsilon z_k, J^{n+1}) e^{-z_k(v_{i+1}^n - v_i^n)/\Delta x} \\ + \sum_{\varepsilon z_k \le -\Delta x} \mathcal{K}(z_k) b(x_i + \varepsilon z_k, J^{n+1}) e^{-z_k(v_i^n - v_{i-1}^n)/\Delta x} \\ + \sum_{\varepsilon |z_k| \ge \Delta x} \mathcal{K}(z_k) b(x_i + \varepsilon z_k, J^{n+1}) e^{(\tilde{v}_{i,k} - v_i^n)/\varepsilon}.$$

Then, Prop. 3.3-(i) and (A4) yield

$$\frac{1}{\Delta t \Delta z} \left| \mathcal{M}_0^{J^{n+1}}(\mathbf{v}^n)_i - \mathcal{M}_\varepsilon^{J^{n+1}}(\mathbf{v}^n)_i \right| \le \varepsilon b_M \sum_{k \in \mathbb{Z}} \mathcal{K}(z_k) |z_k| e^{|z_k| \text{Lip}_x(t^n)} \\ + 2b_M \sum_{\varepsilon |z_k| \ge \Delta x} \mathcal{K}(z_k) e^{|z_k| \text{Lip}_x(t^n)},$$

so that $\| \mathcal{M}_0^{J^{n+1}}(\mathbf{v}^n) - \mathcal{M}_\varepsilon^{J^{n+1}}(\mathbf{v}^n) \|_{\infty} \to_{\varepsilon \to 0} 0$. Equation (A.5) is a consequence of these estimates, and property $\min \mathbf{u}^{n+1, I^{n+1}} = 0$ comes from Prop. 3.3-(iv).

To conclude the proof, one has to show that all extractions of $(J^{n+1})_{\varepsilon \in (0, \varepsilon_0)}$ converge to the same limit. We argue by contradiction, supposing that there are two extractions converging respectively to $I_a^{n+1}$ and $I_b^{n+1}$, with $I_a^{n+1} < I_b^{n+1}$. With the same notations as above, it provides two extractions of $(v_i^{n+1})_{\varepsilon \in (0, \varepsilon_0)}$, which converge respectively to $u_i^{n+1, I_a^{n+1}}$ and $u_i^{n+1, I_b^{n+1}}$ when $\varepsilon \to 0$, with $\min \mathbf{u}^{n+1, I_a^{n+1}} = \min \mathbf{u}^{n+1, I_b^{n+1}} = 0$. However, considering

$$u_i^{n+1, I_a^{n+1}} - u_i^{n+1, I_b^{n+1}} = \mathcal{M}_0^{I_a^{n+1}}(\mathbf{u}^n)_i - \mathcal{M}_0^{I_b^{n+1}}(\mathbf{u}^n)_i + \Delta t (R(t^{n+1}, x_i, I_a^{n+1}) - R(t^{n+1}, x_i, I_b^{n+1})),$$

at the minimum point of $\mathbf{u}^{n+1, I_b^{n+1}}$ yields a contradiction, since $I \mapsto \mathcal{M}_0^I(\mathbf{u}^n)_i$ is non-decreasing (see Remark 3.4), and thanks to (A6).

## References

[1] L. Almeida, B. Perthame, and X. Ruan. An Asymptotic Preserving Scheme for Capturing Concentrations in Age-structured Models Arising in Adaptive Dynamics. HAL-02438316, 2020. working paper or preprint.
[2] L. Ambrosio, N. Fusco, and D. Pallara. *Functions of bounded variations and free discontinuities problems*. Oxford University Press, 2000.

[3] G. Barles. *Solutions de viscosité des équations de Hamilton-Jacobi*, volume 17 of *Mathématiques & Applications (Berlin) [Mathematics & Applications]*. Springer-Verlag, Paris, 1994.

[4] G. Barles, S. Mirrahimi, and B. Perthame. Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result. *Methods Appl. Anal.*, 16(3):321–340, 2009.

[5] G. Barles and B. Perthame. Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics. In *Recent developments in nonlinear partial differential equations*, volume 439 of *Contemp. Math.*, pages 57–68. Amer. Math. Soc., Providence, RI, 2007.

[6] V. Calvez, S. Figueroa Iglesias, H. Hivert, S. Méléard, A. Melnykova, and S. Nordmann. Horizontal gene transfer: numerical comparison between stochastic and deterministic approaches. In *CEMRACS 2018—numerical and mathematical modeling for biological and medical applications: deterministic, probabilistic and statistical descriptions*, volume 67 of *ESAIM Proc. Surveys*, pages 135–160. EDP Sci., Les Ulis, 2020.

[7] V. Calvez, H. Hivert, and H. Yoldaş. Concentration in Lotka-Volterra parabolic equations: an asymptotic-preserving scheme. working paper or preprint, April 2022.

[8] V. Calvez and K.-Y. Lam. Uniqueness of the viscosity solution of a constrained Hamilton-Jacobi equation. *Calc. Var. Partial Differential Equations*, 59(5):Paper No. 163, 22, 2020.

[9] J. Antonio Carrillo, S. Cuadrado, and B. Perthame. Adaptive dynamics via Hamilton-Jacobi approach and entropy methods for a juvenile-adult model. *Math. Biosci.*, 205(1):137–161, 2007.

[10] M. Costa, C. Etchegaray, and S. Mirrahimi. Survival criterion for a population subject to selection and mutations; application to temporally piecewise constant environments. *Nonlinear Anal. Real World Appl.*, 59:Paper No. 103239, 36, 2021.

[11] M. G. Crandall, H. Ishii, and P.-L. Lions. User's guide to viscosity solutions of second order partial differential equations. *Bull. Amer. Math. Soc. (N.S.)*, 27(1):1–67, 1992.

[12] M. G. Crandall and P. L. Lions. Two approximations of solutions of Hamilton-Jacobi equations. *Mathematics of Computation*, 43(167):1–19, 1984.

[13] L. Desvillettes, P.-E. Jabin, S. Mischler, and G. Raoul. On selection dynamics for continuous structured populations. *Commun. Math. Sci.*, 6(3):729–747, 2008.

[14] O. Diekmann. A beginner's guide to adaptive dynamics. In *Mathematical modelling of population dynamics*, volume 63 of *Banach Center Publ.*, pages 47–86. Polish Acad. Sci. Inst. Math., Warsaw, 2004.

[15] O. Diekmann, P.-E. Jabin, S. Mischler, and B. Perthame. The dynamics of adaptation: An illuminating example and a Hamilton-Jacobi approach. *Theoretical Population Biology*, 67(4):257–271, 2005.

[16] G. Dimarco and L. Pareschi. Numerical methods for kinetic equations. *Acta Numer.*, 23:369–520, 2014.

[17] L. C. Evans. *Partial differential equations*. American Mathematical Society, Providence, R.I., 2010.

[18] S. Figueroa Iglesias and S. Mirrahimi. Long time evolutionary dynamics of phenotypically structured populations in time-periodic environments. *SIAM J. Math. Anal.*, 50(5):5537–5568, 2018.

[19] S. Figueroa Iglesias and S. Mirrahimi. Selection and mutation in a shifting and fluctuating environment. *Commun. Math. Sci.*, 19(7):1761–1798, 2021.

[20] S. A. H. Geritz, E. Kisdi, G. Meszéna, and J. A. J. Metz. Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree. *Evolutionary Ecology*, 12:35–57, 1998.

[21] S. A. H. Geritz, J. A. J. Metz, E. Kisdi, and G. Meszéna. Dynamics of adaptation and evolutionary branching. *Phys. Rev. Lett.*, 78:2024–2027, Mar 1997.

[22] J. Guerand and M. Koumaiha. Error estimates for a finite difference scheme associated with Hamilton-Jacobi equations on a junction. *Numer. Math.*, 142(3):525–575, 2019.
[23] H. Hivert. A first-order asymptotic preserving scheme for front propagation in a one-dimensional kinetic reaction-transport equation. *J. Comput. Phys.*, 367:253–278, 2018.

[24] S. Jin. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. *SIAM J. Sci. Comput.*, 21(2):441–454, 1999.

[25] S. Jin. Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review. *Riv. Math. Univ. Parma (N.S.)*, 3(2):177–216, 2012.

[26] Y. Kim. On the uniqueness of solutions to one-dimensional constrained Hamilton-Jacobi equations. *Minimax Theory Appl.*, 6(1):145–154, 2021.

[27] A. Klar. An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit. *SIAM J. Numer. Anal.*, 35(3):1073–1094, 1998.

[28] A. Klar. An asymptotic preserving numerical scheme for kinetic equations in the low mach number limit. *SIAM J. Numer. Anal.*, 36(5):1507–1527, 1999.

[29] T. Lorenzi and C. Pouchol. Asymptotic analysis of selection-mutation models in the presence of multiple fitness peaks. *Nonlinearity*, 33(11):5791–5816, 2020.

[30] A. Lorz, T. Lorenzi, J. Clairambault, A. Escargueil, and B. Perthame. Modeling the effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors. *Bull. Math. Biol.*, 77(1):1–22, 2015.

[31] A. Lorz, S. Mirrahimi, and B. Perthame. Dirac mass dynamics in multidimensional nonlocal parabolic equations. *Comm. Partial Differential Equations*, 36(6):1071–1098, 2011.

[32] G. Meszéna, M. Gyllenberg, F. J. Jacobs, and J. A. J. Metz. Link between population dynamics and dynamics of darwinian evolution. *Phys. Rev. Lett.*, 95:078105, Aug 2005.

[33] J. A. J. Metz, S. A. H. Geritz, G. Meszéna, F. J. A. Jacobs, and J. S. van Heerwaarden. Adaptive dynamics, a geometrical study of the consequences of nearly faithful reproduction. In *Stochastic and spatial structures of dynamical systems (Amsterdam, 1995)*, volume 45 of *Konink. Nederl. Akad. Wetensch. Verh. Afd. Natuurk. Eerste Reeks*, pages 183–231. North-Holland, Amsterdam, 1996.

[34] S. Mirrahimi and J.-M. Roquejoffre. A class of Hamilton-Jacobi equations with constraint: uniqueness and constructive approach. *J. Differential Equations*, 260(5):4717–4738, 2016.

[35] S. Nordmann, B. Perthame, and C. Taing. Dynamics of concentration in a population model structured by age and a phenotypical trait. *Acta Appl. Math.*, 155:197–225, 2018.

[36] B. Perthame. *Transport Equations in Biology*. Frontiers in Mathematics. Birkhäuser Basel, 2007.

[37] B. Perthame and G. Barles. Dirac concentrations in Lotka-Volterra parabolic PDEs. *Indiana University Mathematics Journal*, 57(7):3275–3301, 2008.

[38] P. E. Souganidis. Approximation schemes for viscosity solutions of Hamilton-Jacobi equations. *Journal of Differential Equations*, 59(1):1–43, 1985.

[39] B. van Rensburg, F. Spill, and M. H. Duong. Adaptive dynamics of diverging fitness optima. *arXiv*-2311.00829, 2023. working paper or preprint.