# A particle method for non-local advection-selection-mutation equations

Frank Ernesto Alvarez and Jules Guilberteau  CEREMADE (CNRS UMR no. 7534), PSL
University, Université Paris - Dauphine, Place du Maréchal De Lattre De
Tassigny, Paris CEDEX 16, 75775, France<EMAIL_ADDRESS>GMM, INSA Toulouse, 135 Avenue de Rangueil, Toulouse, 31000, France. alvarez-
<EMAIL_ADDRESS>Sorbonne Université, CNRS, Université Paris Cité, Inria,
Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France.
<EMAIL_ADDRESS>
###### Abstract

The well-posedness of a non-local advection-selection-mutation problem
deriving from adaptive dynamics models is shown for a wide family of initial
data. A particle method is then developed, in order to approximate the
solution of such problem by a regularised sum of weighted Dirac masses whose
characteristics solve a suitably defined ODE system. The convergence of the
particle method over any finite interval is shown and an explicit rate of
convergence is given. Furthermore, we investigate the asymptotic-preserving
properties of the method in large times, providing sufficient conditions for
it to hold true as well as examples and counter-examples. Finally, we
illustrate the method in two cases taken from the literature.

Keywords— Adaptive dynamics, non-local advection, Particle methods, Asymptotic
preservation

##  1 Introduction

### Presentation of the model

The goal of this paper is to develop a numerical method allowing to
approximate the solutions of equations of the form

$\displaystyle\begin{cases}\partial_{t}v(t,x)+\nabla_{x}\cdot(a(t,x,I_{a}v(t,x))v(t,x))=R(t,x,I_{g}v(t,x))v(t,x)+\displaystyle{\int\limits_{\mathbb{R}^{d}}}{m(t,x,y,I_{d}v(t,x))v(t,y)dy},\\\
v\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d})),\\\ v(0,\cdot)=v^{0}(\cdot)\in
W^{1,1}(\mathbb{R}^{d}),\end{cases}$ (1)

where

$(I_{l}u)(t,x)=\int\limits_{\mathbb{R}^{d}}\psi_{l}(t,x,y)u(t,y)dy,\quad
l=a,g,d$

are non-local terms and $a$, $R$, $m$ and $\psi_{l}$ are smooth functions.

This general formulation aims to bring together a wide family of PDE models
typically used in the field of adaptive dynamics. In this context, $x$
represents a phenotypic trait (usually simply called ‘phenotype’ or ‘trait’)
which is a characteristic inherent to individuals, and $x\mapsto v(t,x)$
represents the density of the studied population at time $t\geq 0$. One
purpose of adaptive dynamics is to understand the combined effect of selection
and mutations (which are usually assumed to be rare and small [1, 2, 3]) on
living populations [4]. The literature concerning phenotype-structured
equations is abundant [5, 6, 7, 8, 9, 10, 11]. The model proposed in this
paper (which includes, among others, the equations studied in [12, 13, 14, 15,
16, 17]) takes into account

  * •

Selection and growth, via the term ‘$R(t,x,I_{g}v(t,x))v(t,x)$’, where $R$ can
be interpreted as the instantaneous growth rate, which depends on the trait
$x$ and the whole population.

  * •

Mutation, via the term
‘$\int_{\mathbb{R}^{d}}{m(t,x,y,I_{d}v(t,x))v(t,y)dy}$’, where the function
$m$ can be seen as the probability density for a cell of trait $y$ to mutate
into a cell of trait $x$.

  * •

Advection, via the term ‘$\nabla_{x}\cdot(a(t,x,I_{a}v(t,x))v(t,x))$’. This
term models how the environment drives the individuals towards specific
regions, as opposed to random mutations. Among others, this term can be used
in order to model a cell differentiation phenomenon.

Mutations can also be modelled through a second order differential operator
such as in [18] and [19]. Laplacian-like terms can be approximated by integral
operators, as shown in [20], which means that, after choosing an appropriate
integral approximation, our analysis could be extended to deal with second
order equations. The other two non local-terms ($I_{a}$ and $I_{g}$) allow to
take into account the influence of the environment, created by the whole
population, over the behaviour of the individuals [21, 16], or competition
between individuals [7].

The long time behaviour of models considering only one phenomenon among
selection, advection and mutation is well-known: broadly speaking, it has been
shown that considering selection alone, or advection alone, leads to
concentration phenomena (towards a finite number of traits) [12, 7, 22],
meaning that the density converges to a sum of Dirac masses, while mutations
by themselves have a smoothing effect [23]. Nevertheless, the combined effects
of these terms remains unclear, and may lead to different and non-intuitive
behaviours. As an example, considering both selection and advection can lead
to convergence either to a Dirac mass or to a continuous function [15, 24],
and considering mutation and selection leads either to convergence to a non-
smooth measure or to a continuous function [14]. Note that this model also
includes the equation studied in [25], which was also approximated with a
particle method.

Upon establishing the well-posedness of (1), this paper is concerned with the
derivation of a particle method inspired by [20], the analysis of its
convergence and asymptotic-preserving properties. However, we must emphasise
the two main novelties with respect to that work: First, the use of non local
terms, which as we will show, poses technical difficulties and affects the
existence of smooth solutions in certain cases. Secondly, the study of the
asymptotic preserving property, which guarantees that, under certain
hypotheses, the long time behaviour of the solution is conserved. As will be
seen, these equations are naturally posed in the space of Radon measures,
making particle methods a natural tool to approximate then. Compared to finite
volume or finite element methods, they are more easily implemented.
Furthermore, a change in model leads to very few changes in the corresponding
code, a clear advantage over other methods.

### Particle method

Particle methods use ODE resolution in order to approximate the solution of
PDEs. This makes them particularly easy to implement, as they only require a
classical ODE solver. The main idea is to seek a sum of weighted Dirac masses,
called particle solution, which is denoted

$v^{N}(t)=\sum\limits_{i=1}^{N}{\alpha_{i}(t)\delta_{x_{i}(t)}},$ (2)

where the weights $\alpha_{i}$ and the points $x_{i}$ are solutions of a
suitable ODE system.

In order to recover a smooth function close to the solution of the studied
PDE, the particle solution needs to be regularised: this is usually done by
means of a convolution with a so-called ‘cut-off function’
$\varphi_{\varepsilon}$ which must satisfy some specific properties. We denote
this regular solution

$v^{N}_{\varepsilon}(t,x)=\sum\limits_{i=1}^{N}{\alpha_{i}(t)\varphi_{\varepsilon}(x-x_{i}(t))},$
(3)

where the scaling parameter $\varepsilon$ is a function of $N$.

This method is especially adapted for the linear advection equation
‘$\partial_{t}v(t,x)+\nabla\cdot\left(a(x)v(t,x)\right)=0$’, but has been
generalised to many other kinds of equations which mostly come from physics
[26], such as diffuson equations [27, 28, 29, 30, 31], advection-diffusion
equations [32, 33], convection-diffusion equations [20], the Navier-Stokes
equation [34, 35] or the Vlasov-Poisson equation [36, 37].

We apply the particle method by following its main three steps, as described
in [38]:

  1. (i)

Particle approximation of the initial data. This first step consists of
approaching the initial condition of $v^{0}$ with a sum of weighted Dirac
masses, $i.e.$ choosing $N\in\mathbb{N}$,
$x_{1}^{0},...,x_{N}^{0}\in\mathbb{R}^{d}$,
$\alpha_{1}^{0},...,\alpha_{N}^{0}\in\mathbb{R}$ such that

$v^{N}_{0}:=\sum\limits_{i=1}^{N}{\alpha_{i}^{0}\delta_{x_{i}^{0}}}\sim
v^{0},$

in the sense of Radon measures, which means that, for any
$\phi\in\mathcal{C}^{0}_{c}(\mathbb{R}^{d})$,

$\sum\limits_{i=1}^{N}{\alpha_{i}^{0}\phi(x_{i}^{0})}\underset{N\to+\infty}{\longrightarrow}\int_{\mathbb{R}^{d}}{\phi(x)v^{0}(x)dx}.$

Assuming that $v^{0}$ has a compact support, a canonical way of choosing these
values is to choose a finite collection of subsets
$\Omega^{0}_{i}\subset\mathrm{supp}\left(v^{0}\right)$ satisfying

$\Omega^{0}_{i}\cap\Omega^{0}_{j}=\emptyset,\mbox{ if }i\neq j,\mbox{ and
}\bigcup\limits_{i\in\\{1,...,N\\}}\Omega^{0}_{i}=\mathrm{supp}\left(v^{0}\right),$

and to take, for any $i\in\\{1,...,N\\}$

$x_{i}^{0}\in\Omega^{0}_{i},\quad
w_{i}^{0}=\lvert\Omega^{0}_{i}\rvert,\quad\nu_{i}^{0}=v^{0}(x_{i}^{0})\quad\alpha_{i}=\nu_{i}w_{i}.$

  2. (ii)

Time-evolution of the particles. By using a weak formulation of the PDE, we
determine the ODE satisfied by the positions (denoted $x_{i}$), the volumes
(denoted $w_{i}$) and the weights (denoted $\nu_{i}$) associated to each of
the $N$ particles, with initial conditions ($x_{i}^{0}$, $w_{i}^{0}$,
$\nu_{i}^{0}$) given at the previous step. The exact ODE, and the way the
particles are correlated with each other depends on the complexity of the PDE.
In the case where the advection term is local ($a(t,x,I)=a(t,x)$), the
positions and the volumes satisfy

$\dot{x}_{i}(t)=a(t,x_{i}(t)),\quad\dot{w}_{i}(t)=\nabla_{x}\cdot(a(t,x_{i}(t)))w_{i}(t)$
(4)

and the formula for the $\nu_{i}$ depends on the selection and the mutation
terms. The method described in the core of this article also allows to use
this method in the case of a non-local advection, which modifies the ODE
satisfied by the positions of the particles. Full formulas are given in
Section 3.

Rewriting $\alpha_{i}=\nu_{i}w_{i}$, with the volumes $w_{i}$ which satisfy
(4) is required for the approximation of the different integral terms. Indeed,
by Liouville’s formula [39], for any $f$ smooth enough,

$\int_{\mathbb{R}^{d}}{f(x)dx}\sim\sum\limits_{i=1}^{N}{f(x_{i}(t))w_{i}(t)}.$

Formally, $v$ and $v^{N}$ (as defined by (2)) are both solution of (1) in the
weak sense, with $v^{N}(0)\sim v^{0}$, which implies that, assuming that the
parameters of the PDE are smooth enough, for any given time $T>0$,
$v^{N}(T)\sim v(T,\cdot)$.

  3. (iii)

Regularisation. In order to transform the discrete measure $v^{N}$ into a
smooth function, we use a regularisation process based on convolution, which
writes as a sum, shown in (3), since the convoluted measure is a sum of Dirac
masses. The function

$\varphi_{\varepsilon}:=\frac{1}{\varepsilon}\varphi\left(\frac{\cdot}{\varepsilon}\right),$

depending on a parameter $\varepsilon>0$, is a scaling of the so-called cut-
off function $\varphi$, which must satisfy some regularity and symmetry
properties (which we specify in Section 4.1). The choice of $\varepsilon$,
which intimately depends on the choice of $N$, is intricate: if $\varepsilon$
is too large, then the solution is ‘over-regularised’, and the scheme loses
its accuracy. Conversely, if $\varepsilon$ is to small, then some of the
particles will be neglected, and the scheme does not converge towards the
solution. Choosing the optimal $\varepsilon$ as a function of $N$ is thus not
a trivial question, and it is possible in some cases to optimise the
convergence rate by improving this regularisation step [40, 41].

### Main results

##### Well-posedness.

We first prove that problem (1) is well-posed, i.e. that for any family of
parameters satisfying some regularity properties, defined in
$\mathcal{C}([0,+\infty),L^{1}(\mathbb{R}^{d}))$. The proof heavily relies on
the use of the characteristic curves $X_{u}(t,y)$, solution to the equation

$\displaystyle\begin{cases}\dot{X}_{u}(t,y)=a(t,X_{u}(t,y),(I_{a}u)(t,X_{u}(t,y)))=:\mathcal{A}_{u}(t,X_{u}(t,y)),\quad
t\in[0,T],\\\ X_{u}(0,y)=y,\end{cases}$

for $y\in\mathbb{R}^{d}$ and $u\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$.
The main difference with respect to the approach taken for similar problems,
as in [20], is the need for continuity results for $X_{u}(t,y)$, not only with
respect to the trait variable $y$, but with respect to $u$ as well. The
required results are stated in Section 2.1 and proved in Appendix A.

The well posedness of the problem for smooth initial data is proved in Section
2.2 using a standard fixed point argument. Moreover, in Section 2.3, we
consider a more general family of initial conditions and we prove that the
regularity of $v(t,\cdot)$ is linked to that of $v^{0}$: more precisely, if
$v^{0}\in W^{k,\infty}(\mathbb{R}^{d})$ with compact support, then
$v\in\mathcal{C}([0,+\infty),W^{k-1,1}(\mathbb{R}^{d}))$. This result can be
improved if the advection is local, i.e. $a(t,x,I)=a(t,x)$: In this case, for
any initial condition in $W^{k,1}(\mathbb{R}^{d})$, the solution $v$ is in
$\mathcal{C}([0,+\infty),W^{k,1}(\mathbb{R}^{d}))$.

##### Particle method: definition and well-posedness.

In Section 3, we define the particle method corresponding to this PDE, by
deriving the ODE system satisfied by the particles ($x_{i},w_{i}$ and
$\nu_{i}$). For the non-local case, the equation we obtain is a coupled system
with infinitely many equations and unknowns. We generalise some classical
results from the Cauchy-Lipschitz theory in order to deal with this problem,
and prove that the ODE system is well posed.

##### Convergence of the particle method.

Section 4 is structured as follows: in Section 4.1 we prove the following
estimate, detailed in Theorem 6:

$\lVert v-v_{\varepsilon}^{h}\rVert_{L^{1}(\mathbb{R}^{d})}\leqslant
C\bigl{(}\varepsilon^{r}+\big{(}\frac{h}{\varepsilon}\big{)}^{\kappa}+h^{\kappa}\bigr{)}\lVert
v^{0}\rVert_{W^{\mu,1}(\mathbb{R}^{d})},\quad\mbox{ for all }\;0\leqslant
t\leqslant T,$

with $h=\frac{1}{N}$, $C>0$, $r\in\mathbb{N}^{*}$ (which depends on the
regularity of $\varphi$), $\kappa\in\mathbb{N}^{*}$ (which depends on $k$ and
the local or non-local nature of $a$) and $\mu\in\mathbb{N}^{*}$ (which
depends on $\kappa$ and $r$). This result ensures the convergence of the
particle method in finite time, upon choosing $\varepsilon$ as an suitable
function of $h$ which allows to obtain a convergence in
$(\varepsilon(h)^{r}+\big{(}\frac{h}{\varepsilon(h)}\big{)}^{\kappa}+h^{\kappa})$
which yields a convergence rate of $h^{\frac{\kappa r}{\kappa+r}}$ after
choosing the optima value of $\varepsilon(h)\sim h^{\frac{\kappa}{\kappa+r}}$.
Nevertheless, in general this scheme is not asymptotic preserving. Therefore,
in subsection 4.2, we show examples for which the scheme is asymptotic
preserving, and others for which it is not. In general, the asymptotic
behaviour of a solution is preserved when it converges to a sum of Dirac
masses, and is not when it converges to a smooth solution.

### Perspectives and open problems

Although a loss of regularity appears to be taking place when advection is non
local, we do not know whether such a loss of regularity does happen in certain
cases. The construction of such an example or, on the contrary, the
improvement of our results in order to prove that, in fact, no regularity is
lost could be a first way to extend our work.
Another open problem is the optimisation of the order of convergence for the
numerical solution, improving upon the order $\frac{\kappa r}{\kappa+r}$
obtained in the present work. The approach from [40], where the local averages
are viewed as point values of an approximation of the solution, and the
regularisation of the solution at time $t>0$ is performed by interpolation
rather than convolution, could be a suitable choice.
Lastly, as mentioned before, another direction could be the extension of our
results in order to deal with second order equations, as done in [20], where
the Laplacian operator is approximated by an appropriate sequence of mutation
kernels.

##  2 The problem

For $T>0$ and $k\in\mathbb{N}$ we consider the functions

$\displaystyle(t,x,I)\mapsto a(t,x,I)$ $\displaystyle\in
W^{1,\infty}\left([0,T],(W^{k+1,\infty}(\mathbb{R}^{d+1}))^{d}\right),$ (5)
$\displaystyle(t,x,I)\mapsto R(t,x,I)$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x},W^{k+1,\infty}_{loc}(\mathbb{R}_{I})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{I},W^{k+1,\infty}(\mathbb{R}^{d}_{x})\right).$
(6)

We consider as well $(t,x,y,I)\mapsto m(t,x,y,I)$ such that

$0\leqslant
m\in\mathcal{C}\left([0,T]\times\mathbb{R}_{x}^{d}\times\mathbb{R}_{y}^{d},W^{k+1,\infty}(\mathbb{R}_{I})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{x}^{d}\times\mathbb{R}_{I},L^{\infty}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{y}^{d}\times\mathbb{R}_{I},\mathcal{C}^{k}_{c}(\mathbb{R}^{d}_{x})\right),$

which is globally Lipschitz with respect to the non local variables and
uniformly compactly supported with respect to the $x$ variable. That is, we
suppose $m$ to satisfy the following hypotheses:

  * •

There exists $\mu>0$ such that

$\sup\limits_{t,x,y}\sum\limits_{i=1}^{k}\sum\limits_{|\alpha|=i}\sum\limits_{j=1}^{k}|\partial^{\alpha}_{x}\partial^{j}_{I}m(t,x,y,I)-\partial^{\alpha}_{x}\partial^{j}_{I}m(t,x,y,J)|\leqslant\mu|I-J|.$
(7)

  * •

There exists a compact set $\mathcal{K}$ such that the function

$M(x):=\sup\limits_{t,y,I}\sum\limits_{i=1}^{k}\sum\limits_{|\alpha|=i}\sum\limits_{j=1}^{k}|\partial^{\alpha}_{x}\partial^{j}_{I}m(t,x,y,I)|,$
(8)

satisfies

$\mbox{supp }M(x)\subset\mathcal{K}.$ (9)

Furthermore, we assume that

$\displaystyle\|M\|_{L^{\infty}(\mathbb{R}^{d})}\leqslant\overline{M}<\infty.$
(10)

From a modelling point of view, assuming $m$ to be uniformly compactly
supported reflects the fact that some traits are realistically out of reach
for a given population, and that, in general, mutations are rare and small.

We remark that hypotheses (9) and (10) imply that

$\|M\|_{L^{1}(\mathbb{R}^{d})}\leqslant|\mathcal{K}|\overline{M}=:K<\infty.$
(11)

For all functions $u\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$ we consider
the linear mappings $I_{a}$, $I_{g}$ and $I_{d}$ which satisfy, for all
$t\in[0,T]$, $x\in\mathbb{R}^{d}$,

$\displaystyle(I_{a}u)(t,x)$
$\displaystyle:=\int\limits_{\mathbb{R}^{d}}\psi_{a}(t,x,y)u(t,y)dy,$ (12)
$\displaystyle(I_{g}u)(t,x)$
$\displaystyle:=\int\limits_{\mathbb{R}^{d}}\psi_{g}(t,x,y)u(t,y)dy,$ (13)
$\displaystyle(I_{d}u)(t,x)$
$\displaystyle:=\int\limits_{\mathbb{R}^{d}}\psi_{d}(t,x,y)u(t,y)dy,$ (14)

where

$\displaystyle\psi_{a}$ $\displaystyle\in
W^{1,\infty}\left([0,T]\times\mathbb{R}^{d}_{x},L^{\infty}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y},W^{k+1,\infty}(\mathbb{R}^{d}_{x})\right),$
(15) $\displaystyle 0<\underline{\psi_{g}}\leqslant\psi_{g}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x},L^{\infty}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y},W^{k+1,\infty}(\mathbb{R}^{d})\right),$
(16) $\displaystyle\psi_{d}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x},L^{\infty}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y},W^{k+1,\infty}(\mathbb{R}^{d})\right),$
(17)

for a certain $\underline{\psi_{g}}>0$. We remark that
$I_{a}u,I_{g}u,I_{d}u\in\mathcal{C}([0,T],L^{\infty}(\mathbb{R}^{d}))$. The
functions $\psi_{a}$ and $\psi_{d}$ do not need to be positive, reflecting
this way how different traits have different impacts (which are not always
beneficial) on the environment, and ultimately on the population itself. On
the other hand, $\psi_{g}$ has to be bounded away from zero. This hypothesis
reflects the fact that, at least for the growth term, all interactions between
individuals are of the same type. These interactions may be interpreted as
either strictly competitive or strictly cooperative. In particular, this means
that only “very” non-local dependence with respect to u are allowed. This
excludes partial densities, for instance on part of the traits.
Lastly, we assume that there exist non-negative constants $I^{*}$ and $r^{*}$
such that, for all $t\in[0,T]$, $x\in\mathbb{R}^{d}$, and $I\geqslant I^{*}$,

$R(t,x,I)+K<-r^{*},$ (18)

uniformly on $t$ and $x$. It is somewhat natural to assume that $R$ is
negative for a large population: for example, if a carrying capacity is
assumed to exist, and the population size is approaching such value, then the
growth rate will inevitably drop to levels where no amount of mutations will
be able to compensate for it.
For a given function $v^{0}\in L^{1}(\mathbb{R}^{d})$, we will study the
existence and uniqueness of solution for the problem

$\displaystyle\begin{cases}\partial_{t}v(t,x)+\nabla_{x}\cdot(a(t,x,I_{a}v(t,x))v(t,x))=R(t,x,I_{g}v(t,x))v(t,x)+\displaystyle{\int\limits_{\mathbb{R}^{d}}}{m(t,x,y,I_{d}v(t,x))v(t,y)dy},\\\
v\in\mathcal{C}([0,T),L^{1}(\mathbb{R}^{d})),\\\
v(0,\cdot)=v^{0}(\cdot).\end{cases}$ (19)

We will show that, under additional hypotheses, either over $a$ or $v^{0}$, we
can guarantee the well-posedness of this problem. In particular, we provide
the results regarding the cases of local advection ($\partial_{I}a=0$) and
non-local advection ($\partial_{I}a\neq 0$). We will see that this distinction
directly affects the set of initial data $v^{0}$ for which the existence of
solutions is guaranteed.

###  2.1 Some bounds over the characteristics

Consider $a$ satisfying (5) and $\psi_{a}$ satisfying (15) for some
$k\geqslant 1$. For all $y\in\mathbb{R}^{d}$ and
$u\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$ we define the characteristic
curve $t\mapsto X_{u}(t,y)$ as the unique solution the following ODE

$\displaystyle\begin{cases}\dot{X}_{u}(t,y)=a(t,X_{u}(t,y),(I_{a}u)(t,X_{u}(t,y)))=:\mathcal{A}_{u}(t,X_{u}(t,y)),\quad
t\in[0,T],\\\ X_{u}(0,y)=y,\end{cases}$ (20)

where $(I_{a}u)(t,x)$ is defined in (12). Since the function $\psi_{a}$
belongs to
$W^{1,\infty}\left([0,T]\times\mathbb{R}^{d}_{x},L^{\infty}(\mathbb{R}^{d}_{y})\right)$,
then $(I_{a}u)(t,x)$ belongs to
$W^{1,\infty}\left([0,T]\times\mathbb{R}^{d}_{x}\right)\subset\mathcal{C}^{0,1}([0,T]\times\mathbb{R}^{d}_{x})$.
The regularity of $a$ then implies that $\mathcal{A}_{u}(t,x)$ is a Lipschitz
function with respect to the $x$ variable, uniformly with respect to $t$,
guaranteeing this way the global existence of solution for (20).

For all $u\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$, we define the norms

$\|u\|_{1}:=\|u\|_{L^{1}([0,T]\times\mathbb{R}^{d})}=\int\limits_{0}^{T}\int\limits_{\mathbb{R}^{d}}|u(t,x)|dxdt\mbox{
and }\|u\|:=\sup\limits_{t\in[0,T]}\|u(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}.$

We present some results involving the characteristics. The proofs for such
results are given in Appendix A.

The first property we describe is the continuity of the family of
characteristics with respect to the spatial variables and the function $u$.

######  Lemma 1.

Let $\psi_{a}$ satisfy (15) and $a$ satisfy (5) for $k=0$. Consider
$y_{1},y_{2}\in\mathbb{R}^{d}$ and
$u_{1},u_{2}\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$. Then there exists a
positive constant $C(T,\|u_{1}\|,\|u_{2}\|)$, such that the solutions
$X_{u_{1}}$ and $X_{u_{2}}$ of (20) satisfy for any $t\in[0,T]$,

$\sum\limits_{j=1}^{d}|X^{j}_{u_{1}}(t,y_{1})-X^{j}_{u_{2}}(t,y_{2})|\leqslant
C(T,\|u_{1}\|,\|u_{2}\|)\left(|y_{1}-y_{2}|+\|u_{1}-u_{2}\|_{1}\right).$

Secondly, we claim that the spatial derivatives of the characteristics remain
bounded by a constant only depending on $T$ and $\|u\|$. We also claim that
the spatial derivatives are continuous with respect to the spatial variables
and the function $u$.

######  Lemma 2.

Let $a$ satisfy (5) and $\psi_{a}$ satisfy (15) for some $k\geqslant 1$.
Consider $u\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$. Then there exists a
positive constant $C(T,\|u\|)$, such that the solution $X_{u}(t,y)$ of (20)
satisfies, for all $t\in[0,T]$ and $y\in\mathbb{R}^{d}$,

$\sum\limits_{i=1}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}X^{j}_{u}(t,y)|\leqslant
C(T,\|u\|).$ (21)

Furthermore, for any two points $y_{1},y_{2}\in\mathbb{R}^{d}$, and any two
functions $u_{1},u_{2}\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$, there
exists a positive constant $C_{2}(T,\|u_{1}\|,\|u_{2}\|)$, such that the
solutions $X_{u_{1}}$ and $X_{u_{2}}$ of (20) satisfy for any $t\in[0,T]$

$\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}X^{j}_{u_{1}}(t,y_{1})-\partial^{\alpha}_{y}X^{j}_{u_{2}}(t,y_{2})|\leqslant
C(T,\|u_{1}\|,\|u_{2}\|)(|y_{1}-y_{2}|+\|u_{1}-u_{2}\|_{1}).$ (22)

Since for all $t\in[0,T]$, $y\mapsto X_{u}(t,y)$ is a
$\mathcal{C}^{1}$-diffeomorphism from $\mathbb{R}^{d}$ onto itself, we may
define its inverse as the function satisfying $X_{u}(t,X^{-1}_{u}(t,x))=x$ for
all $(t,x)\in[0,T]\times\mathbb{R}^{d}$. We have the following results for
$X^{-1}_{u}$.

######  Lemma 3.

Let $a$ satisfy (5) and $\psi_{a}$ satisfy (15) for some $k\geqslant 1$.
Consider $u\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$. Then there exists a
positive constant $\widetilde{C}(T,\|u\|)$, such that the inverse of the
solution $X_{u}(t,y)$ of (20) satisfies, for all $t\in[0,T]$,
$x\in\mathbb{R}^{d}$

$\sum\limits_{i=1}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{x}\left(X^{-1}_{u}\right)^{j}(t,x)|\leqslant\widetilde{C}(T,\|u\|).$
(23)

######  Lemma 4.

Let $a$ satisfy (5) and $\psi_{a}$ satisfy (15) for some $k\geqslant 1$.
Consider any two functions
$u_{1},u_{2}\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$. Then there exists a
positive constant $\widetilde{C}(T,\|u_{1}\|,\|u_{2}\|)$, which satisfies
$\lim\limits_{T\rightarrow 0}\widetilde{C}(T,\|u_{1}\|,\|u_{2}\|)=0$ and such
that the inverses of the solutions $X_{u_{1}}$ and $X_{u_{2}}$ of (20)
satisfy, for all $t\in[0,T]$ and $x\in\mathbb{R}^{d}$

$\sum\limits_{i=0}^{k-1}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\left(\partial^{\alpha}_{x}X^{-1}_{u_{1}}\right)^{j}(t,x)-\partial^{\alpha}_{x}\left(X^{-1}_{u_{2}}\right)^{j}(t,x)|\leqslant\widetilde{C}(T,\|u_{1}\|,\|u_{2}\|)\|u_{1}-u_{2}\|_{1}.$
(24)

We remark that thanks to the relation $\|u\|_{1}\leqslant T\|u\|$, the
relations (22) and (24) also hold true when replacing $\|u_{1}-u_{2}\|_{1}$ by
$\|u_{1}-u_{2}\|$. Lastly, we give a result regarding the regularity of
$X_{u}(t,x)$ with respect to $t$.

######  Lemma 5.

Let $a$ satisfy (5) and $\psi_{a}$ satisfy (15) for some $k\geqslant 1$.
Consider $u\in\mathcal{C}^{1}([0,T_{1})\times\mathbb{R}^{d})$ such that
$\sup\limits_{t\in[0,T_{1})}(\|u(t,\cdot)\|_{L^{1}(R^{d})}+\|\partial_{t}u(t,\cdot)\|_{L^{1}(R^{d})})<+\infty$.
Then
$X_{u}(t,y)\in\mathcal{C}^{1}([0,T_{1}],\mathcal{C}^{k}(\mathbb{R}^{d}))$. As
a consequence,
$X^{-1}_{u}(t,x)\in\mathcal{C}^{1}([0,T_{1}],\mathcal{C}^{k}(\mathbb{R}^{d}))$.

###### Proof.

Thanks to Lemma 2, we know that under these hypotheses, for all $T<T_{1}$,
$X_{u}(t,y)$ exists and belongs to
$\mathcal{C}^{1}([0,T],\mathcal{C}^{k}(\mathbb{R}^{d}))$. Consider
$0<t_{1},t_{2}<T_{1}$, then

$\displaystyle\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}X^{j}_{u}(t_{1},y)-\partial^{\alpha}_{y}X^{j}_{u}(t_{2},y)|$
$\displaystyle=\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\int_{t_{1}}^{t_{2}}\partial^{\alpha}_{y}\dot{X}^{j}_{u}(s,y)ds|$
$\displaystyle=\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\int_{t_{1}}^{t_{2}}\partial^{\alpha}_{y}a_{j}(s,X_{u}(s,y),(I_{a}u)(s,X_{u}(s,y)))ds|.$

Thanks to the regularity of $a$ and $\psi_{a}$, the bounds given in Lemma 2
for the derivatives of $X_{u}(t,y)$ and the uniform bound for
$\|u\|_{L^{1}(\mathbb{R}^{d})}$ we conclude that there exists a positive
constant such that

$\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}X^{j}_{u}(t_{1},y)-\partial^{\alpha}_{y}X^{j}_{u}(t_{2},y)|\leqslant
C|t_{1}-t_{2}|.$

Similarly, we have

$\displaystyle\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}\dot{X}^{j}_{u}(t_{1},y)-\partial^{\alpha}_{y}\dot{X}^{j}_{u}(t_{2},y)|$
$\displaystyle=\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}(\mathcal{A}_{u})_{j}(t_{1},y)-\partial^{\alpha}_{y}(\mathcal{A}_{u})_{j}(t_{2},y)|,$

where

$(\mathcal{A}_{u})_{j}(t,y):=a_{j}(t,X_{u}(t,y),(I_{a}u)(t_{1},X_{u}(t,y))).$

Again, the regularity up to order $k+1$ of the involved coefficients allow us
to conclude that

$\sum\limits_{i=0}^{k}\sum\limits_{|\alpha|\leqslant
i}\sum\limits_{j=1}^{d}|\partial^{\alpha}_{y}\dot{X}^{j}_{u}(t_{1},y)-\dot{X}^{j}_{u}(t_{1},y)|\leqslant
C|t_{1}-t_{2}|.$

We have shown that $X_{u}(t,\cdot)$ is a Cauchy sequence in
$\mathcal{C}^{k}(\mathbb{R}^{d})$ when $t$ goes to $T_{M}$, therefore, there
exists $X^{*}(x)\in\mathcal{C}^{k}(\mathbb{R}^{d})$ and
$Y^{*}(x)\in\mathcal{C}^{k}(\mathbb{R}^{d})$ such that

$\lim\limits_{t\rightarrow
T_{M}}\|X_{u}(t,\cdot)-X^{*}\|_{\mathcal{C}^{k}(\mathbb{R}^{d})}+\|\dot{X}_{u}(t,\cdot)-Y^{*}\|_{\mathcal{C}^{k}(\mathbb{R}^{d})},$

which is the desired result. ∎

###  2.2 Existence of solution for smooth initial data

We first provide the proof of existence and uniqueness of solution for problem
(19) when the initial condition $v^{0}$ is a smooth enough function. We still
assume hypotheses (5) through (18) hold. For a smooth initial condition
$v^{0}$ we denote by _classical solution_ any function
$v\in\mathcal{C}^{1}([0,T]\times\mathbb{R}^{d})$ which satisfies problem (19).
The following _a priori_ estimate will allow us to guarantee the global
existence of a classical solution given that such solution exists over a
certain interval $[0,T_{1}]$.

######  Lemma 6.

Let $v^{0}\in\mathcal{C}^{1}_{c}(\mathbb{R}^{d})$ and $T_{1}>0$ be such that a
classical solution $v\in\mathcal{C}^{1}([0,T]\times\mathbb{R}^{d})$ exists for
problem (19) for all $T<T_{1}$, which is positive and has compact support with
respect to the $x$ variable. Then, such solution satisfies the estimate

$\sup\limits_{t\in[0,T_{1}]}\|v(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}\leqslant\max\\{\|v^{0}\|_{L^{1}(\mathbb{R}^{d})},\frac{I^{*}}{\underline{\psi_{g}}}\\}.$
(25)

###### Proof.

Let $v$ be the aforementioned positive solution. Denoting
$\rho(t):=\|v(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}$, we see directly from
equation (19) that

$\dot{\rho}(t)=\int\limits_{\mathbb{R}^{d}}\left(R(t,y,(I_{g}v)(t,y))+\int_{\mathbb{R}^{d}}m(t,x,y,(I_{d}v)(t,x))dx\right)v(t,y)dy.$

If there exists $t$ such that $\rho(t)>\frac{I^{*}}{\underline{\psi_{g}}}$,
then $(I_{g}v)(t,y)\geqslant\underline{\psi_{g}}\rho(t)>I^{*}$, which allows
us to use hypothesis (18) in order to conclude

$\dot{\rho}(t)<-r^{*}\rho(t)<0.$

This way, we see that either $\rho(t)$ is smaller than
$\frac{I^{*}}{\underline{\psi_{g}}}$ or $\rho(t)$ is decreasing, which in turn
implies the bound (25). ∎

A fixed point argument together with estimate (25) will allow us to conclude
the existence of solution for problem (19) for smooth initial conditions.

######  Theorem 1.

Consider $k\geqslant 1$ and $T>0$. Under hypotheses (15) through (18), for all
non-negative functions $v^{0}\in\mathcal{C}^{k}_{c}(\mathbb{R}^{d})$, there
exists a unique non-negative classical solution
$v\in\mathcal{C}^{1}([0,T],\mathcal{C}^{k}_{c}(\mathbb{R}^{d}))$ to problem
(19). Furthermore, such solution satisfies

$\displaystyle\sup\limits_{t\in[0,T]}\|v(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}\leqslant\max\\{\|v^{0}\|_{L^{1}(\mathbb{R}^{d})},\frac{I^{*}}{\underline{\psi_{g}}}\\},$
(26)
$\displaystyle\sup\limits_{t\in[0,T]}\|v(t,\cdot)\|_{W^{k,1}(\mathbb{R}^{d})}\leqslant
C_{T}\|v^{0}\|_{W^{k,1}(\mathbb{R}^{d})}.$ (27)

###### Proof.

Consider $v^{0}\in\mathcal{C}^{k}_{c}(\mathbb{R}^{d})$. For $t\geq 0$ and the
function $M$ introduced in (8) we define $r_{t}:=2\|a\|_{L^{\infty}}t$,
$B_{r_{t}}$ the open ball centred at $0$ and of radius $r_{t}$,

$O_{t}:=\mbox{supp }(v^{0})\cup\mbox{supp }(M)+B_{r_{t}},$

and for $\alpha>1$ we define
$\overline{\rho}_{\alpha}:=\max\\{\alpha\|v^{0}\|_{L^{1}(\mathbb{R}^{d})},\frac{I^{*}}{\underline{\psi_{g}}}\\}$.
We consider

$u\in
D^{T}_{\alpha}:=\left\\{u\in\mathcal{C}([0,T]\times\mathbb{R}^{d}):u\geqslant
0,\mbox{supp }(u(t,\cdot))\subset
O_{t},\int_{\mathbb{R}^{d}}u(t,x)dx\leqslant\overline{\rho}_{\alpha},\ \forall
t\in[0,T]\right\\}.$

We denote as $\Phi u$ the mapping defined by $v=\Phi u$, where $v$ is the
solution of

$\displaystyle\begin{cases}\partial_{t}v(t,x)+\nabla_{x}\cdot(a(t,x,(I_{a}u)(t,x))v(t,x))-R(t,x,(I_{g}u)(t,x))v(t,x)=\int_{\mathbb{R}^{d}}{m(t,x,y,(I_{d}u)(t,x))u(t,y)dy}\\\
v(0,\cdot)=v^{0}(\cdot).\end{cases}$

Let us denote, for any $y\in\mathbb{R}^{d}$, as $X_{u}(\cdot,y)$ the unique
solution of

$\displaystyle\begin{cases}\dot{X}_{u}(t,y)=a(t,X_{u}(t,y),(I_{a}u)(t,X_{u}(t,y)))=:\mathcal{A}_{u}(t,X_{u}(t,y))\quad
t\geq 0,\\\ X(0,y)=y.\end{cases}$

As stated before, for any $t\in[0,T]$, $x\mapsto X_{u}(t,x)$ is a
$\mathcal{C}^{1}$-diffeomorphism, therefore, for all $t\geq 0$ and all
$x\in\mathbb{R}^{d}$ there exists a unique $y\in\mathbb{R}^{d}$ such that
$x=X_{u}(t,y)$, which we denote $y=X^{-1}_{u}(t,x)$.
We see that,

$\displaystyle\frac{d}{dt}v(t,X_{u}(t,y))=$
$\displaystyle\big{[}R(t,X_{u}(t,y),(I_{g}u)(t,X_{u}(t,y)))-\text{div
}\mathcal{A}_{u}(t,X_{u}(t,y))\big{]}v(t,X_{u}(t,y))$
$\displaystyle+\int_{\mathbb{R}^{d}}{m(t,X_{u}(t,y),z,(I_{d}u)(t,X_{u}(t,y)))u(t,z)dz},$

and thus, denoting

$\mathcal{G}_{u}(t,y):=R(t,X_{u}(t,y),(I_{g}u)(t,X_{u}(t,y)))-\text{div
}\mathcal{A}_{u}(t,X_{u}(t,y)),$

we get

$\displaystyle\Phi u(t,X_{u}(t,y))=$ $\displaystyle v(t,X_{u}(t,y))$
$\displaystyle=$ $\displaystyle
v^{0}(y)\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{u}(s,y)ds\biggr{)}$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}m(s,X_{u}(s,y),z,(I_{d}u)(s,X_{u}(s,y)))u(s,z)dz\exp\left(\int_{s}^{t}\mathcal{G}_{u}(\tau,y)d\tau\right)ds$

or, equivalently

$\displaystyle\Phi u(t,x)=$ $\displaystyle v(t,x)$ $\displaystyle=$
$\displaystyle
v^{0}(X^{-1}_{u}(t,x))\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{u}(s,X^{-1}_{u}(t,x))ds\biggr{)}$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}m(s,X_{u}(s,X^{-1}_{u}(t,x)),z,(I_{d}u)(s,X_{u}(s,X^{-1}_{u}(t,x))))u(s,z)dz\exp\left(\int_{s}^{t}\mathcal{G}_{u}(\tau,X^{-1}_{u}(t,x))d\tau\right)ds.$

The solution $v$ is thus non-negative, according to the non-negativity of
$v^{0}$, $u$ and $m$. Thanks to Lemma 3, $v$ belongs to
$\mathcal{C}^{1}([0,T],\mathcal{C}^{\kappa}(\mathbb{R}^{d}))\subset\mathcal{C}([0,T]\times\mathbb{R}^{d})$.
Furthermore, the fact that for all $t\in[0,T]$,
$|X(t,y)-y|\leqslant\|a\|_{L^{\infty}}t$ implies that $\mbox{supp
}(v(t,\cdot))\subset O_{t}$.
Additionally, directly from its definition, we see that

$\mathcal{G}_{u}(s,X^{-1}_{u}(t,x))\leqslant\gamma:=\|R\|_{L^{\infty}_{t,x,I}}+\|a\|_{W^{1,\infty}_{x}L^{\infty}_{t,I}}+\|a\|_{W^{1,\infty}_{I}L^{\infty}_{t,x}}\|\psi_{a}\|_{W^{1,\infty}_{x}L^{\infty}_{y}}\overline{\rho}_{\alpha},$

and consequently, for all $u\in D^{T}_{\alpha}$, we have

$\displaystyle\|\Phi u(t,\cdot))\|_{L^{1}(\mathbb{R}^{d})}\leqslant$
$\displaystyle e^{\gamma
T}\left(\int\limits_{\mathbb{R}^{d}}v^{0}(X^{-1}_{u}(t,x))dx\right.$
$\displaystyle\left.+\int_{0}^{t}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}m(s,X_{u}(s,X^{-1}_{u}(t,x)),z,(I_{d}u)(s,X_{u}(s,X^{-1}_{u}(t,x))))u(s,z)dzdxds\right)$
$\displaystyle\leqslant$ $\displaystyle e^{\gamma
T}\left(\int\limits_{\mathbb{R}^{d}}v^{0}(X^{-1}_{u}(t,x))dx+\int_{0}^{t}\int_{\mathbb{R}^{d}}M(s,X_{u}(s,X^{-1}_{u}(t,x)))dx\int_{\mathbb{R}^{d}}u(s,z)dzds\right).$

Making the changes of variables $y=X^{-1}_{u}(t,x)$ and
$y=X_{u}(s,X^{-1}_{u}(t,x))$ respectively on each of the integrals on the last
expression, recalling that, according to Liouville’s formula

$\displaystyle|J_{X_{u}(t,y)}|$ $\displaystyle=e^{\int_{0}^{t}\text{div
}\mathcal{A}_{u}(s,X_{u}(t,y))ds},$
$\displaystyle|J_{X_{u}(t,X^{-1}_{u}(s,y))}|$
$\displaystyle=e^{\int_{0}^{t}\text{div
}\mathcal{A}_{u}(\tau,X_{u}(t,X^{-1}_{u}(s,y)))d\tau-\int_{0}^{s}\text{div
}\mathcal{A}_{u}(\tau,y)d\tau},$

and using the hypotheses over $a$ and $m$ we obtain that for all $t\in[0,T]$,

$\displaystyle\|\Phi u(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}$
$\displaystyle\leqslant
e^{(\gamma+2\tilde{a})T}\left(\|v^{0}\|_{L^{1}(\mathbb{R}^{d})}+KT\|u\|_{L^{1}(\mathbb{R}^{d})}\right),$

where
$\tilde{a}:=\|a\|_{W^{1,\infty}_{x}L^{\infty}_{t,I}}+\|a\|_{W^{1,\infty}_{I}L^{\infty}_{t,x}}\|\psi_{a}\|_{W^{1,\infty}_{x}L^{\infty}_{y}}\overline{\rho}_{\alpha}$.
Finally, using the hypothesis over $v^{0}$ and $u$ we see that

$\|\Phi u(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}\leqslant
e^{(\gamma+2\tilde{a})T}\left(\frac{1}{\alpha}+KT\right)\overline{\rho}_{\alpha}.$

Thanks to the condition $\alpha>1$ there exists $T_{\alpha}$ (only depending
on $\alpha$ and on the coefficients of the problem) such that $\|\Phi
u\|_{L^{1}(\mathbb{R}^{d})}\leqslant\overline{\rho}_{\alpha}$ for all
$t\in[0,T_{\alpha}]$. In other words, $\Phi:D^{T_{\alpha}}_{\alpha}\rightarrow
D^{T_{\alpha}}_{\alpha}$.

We now claim that the mapping $\Phi u$ is a contraction on some
$D^{T_{1}}_{\alpha}$, $0<T_{1}\leqslant T_{\alpha}$, with respect to the usual
norm in $\mathcal{C}([0,T]\times\mathbb{R}^{d})$. For two functions
$u_{1},u_{2}\in D^{T_{\alpha}}_{\alpha}$ and any $t\in[0,T_{\alpha}]$, thanks
to Lemma 1 and Lemma 4, we have

$\displaystyle|\Phi u_{1}-\Phi u_{2}|\leqslant$ $\displaystyle Ce^{\gamma
t}|X^{-1}_{u_{1}}(t,x)-X^{-1}_{u_{2}}(t,x)|$ $\displaystyle\leqslant$
$\displaystyle Ce^{\gamma t}\|u_{1}-u_{2}\|_{1}$ $\displaystyle\leqslant$
$\displaystyle Ce^{\gamma
t}|O_{t}|t\|u_{1}-u_{2}\|_{\mathcal{C}([0,t]\times\mathbb{R}^{d})}.$

Clearly, for $t=T_{1}$ small enough, $\Phi u$ is a contraction, and therefore,
thanks to the Banach fixed point theorem, there exists a unique $v\in
D^{T_{1}}_{\alpha}$ such that $\Phi v=v$. Such $v$ is a solution of (19) over
$[0,T_{1}]$. Furthermore, directly from the relation $v=\Phi v$ we see that
$v\in\mathcal{C}^{1}([0,T],\mathcal{C}^{\kappa}_{c}(\mathbb{R}^{d}))$.

Let us now assume that there exists $T_{M}$, a finite maximal time such that a
solution exists in $B^{T}_{\alpha}$ for all $T<T_{M}$, and let $v$ be such
solution. Directly from Lemma 6, the solution $v$ satisfies estimate (25) over
$[0,T_{M})$. Furthermore, thanks to the relation $v=\Phi v$, we are able to
show that

$\sup\limits_{[0,T_{M})}(\|v(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}+\|\dot{v}(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})})<+\infty.$

From Lemma 5 we get then that
$X^{-1}_{v}(t,x)\in\mathcal{C}^{1}([0,T_{M}],\mathcal{C}^{k}(\mathbb{R}^{d}))$,
and by composition of functions, so is $v=\Phi v$. We can then iterate the
previous ideas using $v(T_{M})\in\mathcal{C}^{k}(\mathbb{R}^{d})$ as a
starting point in order to obtain the existence of solution over a certain
interval $[T_{M},T_{M}+\delta)$, contradicting this way the maximal character
of $T_{M}$. Hence, there exists a classical solution of (19) for all $t>0$.

The uniqueness on
$\mathcal{C}([0,T],\mathcal{C}^{\kappa}_{c}(\mathbb{R}^{d}))$ comes from the
fact that every other solution on $D^{\infty}_{\alpha}$ will coincide with $v$
at least over a small interval $(0,t_{0})$ and then, by continuity, the same
would hold for all $t$.

To obtain the $W^{k,1}(\mathbb{R}^{d})$ estimates, we differentiate the
relation $v=\Phi v$ and notice that for all multi-index $\beta$ such that
$|\beta|\leqslant k$

$\displaystyle\partial^{\beta}_{x}v=$
$\displaystyle\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{v}(s,X^{-1}_{v}(t,x))ds\biggr{)}\sum\limits_{|\gamma|\leqslant|\beta|}\partial^{\gamma}_{x}v^{0}(X^{-1}_{v}(t,x))F^{\gamma}_{1}(t,x)$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}F_{2}(s,x,y)v(s,y)dy\exp\left(\int_{s}^{t}\mathcal{G}_{v}(\tau,x)d\tau\right)ds$

where the functions $F^{\gamma}_{1}$ and $F_{2}$ are combinations of sums and
multiplications of the derivatives up to order $|\beta|$ of $X^{-1}_{v}$, $R$,
$m$, $I_{g}$ and $I_{d}$. Taking absolute values, integrating over
$\mathbb{R}^{d}$ and using the boundedness of all the involved coefficients we
arrive at

$\|v(t,\cdot)\|_{W^{k,1}(\mathbb{R}^{d})}\leqslant
C^{1}_{T}\|v^{0}\|_{W^{k,1}(\mathbb{R}^{d})}+C^{2}_{T}\int_{0}^{t}\|v(s,\cdot)\|_{W^{k,1}(\mathbb{R}^{d})}ds$

and thanks to Grönwall’s lemma we get (27). ∎

###  2.3 Existence of solution for more general initial data

Depending on whether $\partial_{I}a=0$ or $\partial_{I}a\neq 0$, we will have
a different class of initial data for which we are able to guarantee existence
of solution for problem (19). Furthermore, the regularity of such solution
might also be affected.
We first prove that, when $\partial_{I}a\neq 0$, a solution exists (in a sense
that will be defined below) for any initial condition $v^{0}\in
W^{k,\infty}(\mathbb{R}^{d})$, with compact support. However, we do not prove
that the regularity of the solution is preserved over time, even if we did not
manage to highlight the existence of cases where a loss of regularity is
observed. Secondly, we will show that, when $\partial_{I}a=0$, not only the
set of initial conditions for which we can claim existence of solution is more
general ($v^{0}\in W^{k,1}(\mathbb{R}^{d})$), but the regularity of such
solution is preserved for all $t>0$.
We introduce the definition of _weak solution_ for problem (19). We say that
$v$ is a _weak solution_ of problem (19) associated to $v^{0}\in
L^{p}(\mathbb{R}^{d})$ if

$v\in L^{\infty}([0,T],L^{p}(\mathbb{R}^{d})),$

and it satisfies the equation in the following weak sense

$\int_{0}^{T}\int_{\mathbb{R}^{d}}vL^{*}_{v}\varphi
dxdt=\int_{\mathbb{R}^{d}}v^{0}\varphi dx,$

for any $\varphi\in\mathcal{C}^{1}_{c}([0,T)\times\mathbb{R}^{d})$, where we
define the operator $L^{*}_{v}$ by

$L^{*}_{v}\varphi(t,x)=-\partial_{t}\varphi(t,x)-a(t,x,(I_{a}v)(t,x))\cdot\nabla\varphi(t,x)-R(t,x,(I_{g}v)(t,x))\varphi(t,x)-\int_{\mathbb{R}^{d}}m(t,y,x,(I_{d}v)(t,y))\varphi(t,y)dy,$

for all $t\in[0,T]$, $y\in\mathbb{R}^{d}$. We remark that a _classical
solution_ is always a _weak solution_.

######  Theorem 2.

Under hypotheses (5) through (18), for all $k\geqslant 1$ and any non-negative
functions $v^{0}\in W^{k,\infty}(\mathbb{R}^{d})$ with compact support, there
exists a unique non-negative weak solution
$v\in\mathcal{C}([0,T],\mathcal{C}^{k-1}_{c}(\mathbb{R}^{d}))$ to problem
(19). Furthermore, such solution satisfies

$\displaystyle\sup\limits_{t\in[0,T]}\|v(t,\cdot)\|_{L^{1}(\mathbb{R}^{d})}\leqslant\max\\{\|v^{0}\|_{L^{1}(\mathbb{R}^{d})},\frac{I^{*}}{\underline{\psi_{g}}}\\},$
(28)
$\displaystyle\sup\limits_{t\in[0,T]}\|v(t,\cdot)\|_{W^{k-1,1}(\mathbb{R}^{d})}\leqslant
C_{T}\|v^{0}\|_{W^{k-1,1}(\mathbb{R}^{d})},$ (29)

and, for $k\geqslant 2$,
$v\in\mathcal{C}^{1}([0,T],\mathcal{C}^{k-1}_{c}(\mathbb{R}^{d}))$.

###### Proof.

Directly from Morrey’s inequality, we get the relation $v^{0}\in
W^{k,\infty}(\mathbb{R}^{d})\subset\mathcal{C}^{k-1,1}(\mathbb{R}^{d})$. If
$k\geqslant 2$, we are able to apply Theorem 1 in order to get the desired
result.
Consider now $k=1$. The compact support of $v^{0}$ and the
$W^{1,\infty}(\mathbb{R}^{d})$ regularity imply that $v^{0}\in
W^{k,1}(\mathbb{R}^{d})$. This means that, there exists a sequence of
compactly supported functions
$v^{0}_{\varepsilon}\in\mathcal{C}^{1}_{c}(\mathbb{R}^{d})$ such that

$\displaystyle\mbox{supp }(v^{0})_{\varepsilon}\subset\mbox{ supp }(v^{0}),$
$\displaystyle\|v^{0}_{\varepsilon}\|_{W^{1,\infty}(\mathbb{R}^{d})}\leqslant\|v^{0}\|_{W^{1,\infty}(\mathbb{R}^{d})},$
$\displaystyle\lim\limits_{\varepsilon\rightarrow
0}\|v^{0}-v^{0}_{\varepsilon}\|_{W^{1,1}(\mathbb{R}^{d})}=0.$

We denote as $v_{\varepsilon}$ the solution of problem (19) associated to
$v^{0}_{\varepsilon}$, and we claim that $v_{\varepsilon}$ is a Cauchy
sequence in $\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$.
We recall that for all $\varepsilon$,

$\underset{t\in[0,T]}{\sup}\;\|v_{\varepsilon}\|_{L^{1}(\mathbb{R}^{d})}\leqslant\overline{\rho}_{\varepsilon}:=\max\\{\|v^{0}_{\varepsilon}\|_{L^{1}(\mathbb{R}^{d})},\frac{I^{*}}{\underline{\psi_{g}}}\\}.$

Furthermore, the equality $v_{\varepsilon}=\Phi v_{\varepsilon}$ holds true,
where $\Phi$ was defined on the proof of Theorem 1. Consequently, for
$\varepsilon_{1},\varepsilon_{2}>0$ we have the relation

$\displaystyle\Delta V_{\varepsilon_{1}\varepsilon_{2}}:=$ $\displaystyle
v_{\varepsilon_{1}}-v_{\varepsilon_{2}}$ $\displaystyle=$ $\displaystyle\Phi
v_{\varepsilon_{1}}-\Phi v_{\varepsilon_{2}}$ $\displaystyle=$
$\displaystyle\biggl{(}v^{0}_{\varepsilon_{1}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{2}}}(t,x))\biggr{)}\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\biggr{)}-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{2}}}(t,x))\biggl{(}\Delta
E(0,t,x)\biggr{)}$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}\biggl{(}\mathcal{M}_{1}(s,t,x,z)-\mathcal{M}_{2}(s,t,x,z)\biggr{)}v_{\varepsilon_{1}}(s,z)dz\exp\left(\int_{s}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\right)ds$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}\mathcal{M}_{2}(s,t,x,z)\biggl{(}\Delta
V_{\varepsilon_{1}\varepsilon_{2}}(s,z)\biggr{)}dz\exp\left(\int_{s}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\right)ds$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}\mathcal{M}_{2}(s,t,x,z)v_{\varepsilon_{j}}(s,z)dz\biggl{(}\Delta
E(s,t,x)\biggr{)}ds,$

where

$\displaystyle\mathcal{G}_{i}(\tau,t,x)$
$\displaystyle=\mathcal{G}_{v_{\varepsilon_{i}}}(\tau,X^{-1}_{v_{\varepsilon_{i}}}(t,x)),$
$\displaystyle\Delta E(s,t,x)$
$\displaystyle:=\exp\left(\int_{s}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\right)-\exp\left(\int_{s}^{t}\mathcal{G}_{2}(\tau,t,x)d\tau\right),$
$\displaystyle\mathcal{M}_{i}(s,t,x,z)$
$\displaystyle:=m(s,X_{v_{\varepsilon_{i}}}(s,X^{-1}_{v_{\varepsilon_{i}}}(t,x)),z,(I_{d}v_{\varepsilon_{i}})(s,X_{v_{\varepsilon_{i}}}(s,X^{-1}_{v_{\varepsilon_{i}}}(t,x)))).$

We write

$\begin{matrix}v^{0}_{\varepsilon_{1}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{2}}}(t,x))&=&v^{0}_{\varepsilon_{1}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))\\\
&&\\\
&&+v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{2}}}(t,x)).\end{matrix}$
(30)

Thanks to the change of variables $y=X^{-1}_{v_{\varepsilon_{1}}}(t,x)$ and
the relations

$\displaystyle|J_{X_{v_{\varepsilon_{1}}}(t,y)}|=$ $\displaystyle
e^{\int_{0}^{t}\text{div
}\mathcal{A}_{v_{\varepsilon_{1}}}(s,X_{v_{\varepsilon_{1}}}(t,y))ds},$
$\displaystyle\mathcal{G}(s,X^{-1}_{v_{\varepsilon_{1}}}(t,x))\leqslant$
$\displaystyle\gamma:=\|R\|_{L^{\infty}_{t,x,I}}+\|a\|_{W^{1,\infty}_{x}L^{\infty}_{t,I}}+\|a\|_{W^{1,\infty}_{I}L^{\infty}_{t,x}}\|\psi_{a}\|_{W^{1,\infty}_{x}L^{\infty}_{y}}\overline{\rho},$
$\displaystyle\text{div
}\mathcal{A}_{v_{\varepsilon_{1}}}(s,X_{v_{\varepsilon_{1}}}(t,y))\leqslant$
$\displaystyle\tilde{a}:=\|a\|_{W^{1,\infty}_{x}L^{\infty}_{t,I}}+\|a\|_{W^{1,\infty}_{I}L^{\infty}_{t,x}}\|\psi_{a}\|_{W^{1,\infty}_{x}L^{\infty}_{y}}\overline{\rho},$

we conclude that

$\displaystyle\int\limits_{\mathbb{R}^{d}}\biggl{(}v^{0}_{\varepsilon_{1}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))\biggr{)}\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\biggr{)}dx$
$\displaystyle=$
$\displaystyle\int\limits_{\mathbb{R}^{d}}\biggl{(}v^{0}_{\varepsilon_{1}}(y)-v^{0}_{\varepsilon_{2}}(y)\biggr{)}\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{v_{\varepsilon_{i}}}(\tau,y)d\tau\biggr{)}|J^{-1}_{X_{v_{\varepsilon_{1}}}(t,y)}|dy$
$\displaystyle\leqslant$ $\displaystyle
e^{(\gamma+\tilde{\alpha})T}\|v^{0}_{\varepsilon_{1}}-v^{0}_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R})}.$
(31)

On the other hand, the compactness of the support of
$v^{0}_{\varepsilon_{2}}$, together with the relation

$|X_{v_{\varepsilon}}(t,y)-y|\leqslant\|a\|_{L^{\infty}}t,$

implies that

$v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{2}}}(t,x))=0,\mbox{
for all }x\not\in O_{t}:=\mbox{supp }v^{0}+B_{r_{t}},$

where $B_{r_{t}}$ is the ball of radius $\|a\|_{L^{\infty}}t$. Hence, we
obtain

$\displaystyle\int\limits_{\mathbb{R}^{d}}\biggl{(}v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{1}}}(t,x))-v^{0}_{\varepsilon_{2}}(X^{-1}_{v_{\varepsilon_{2}}}(t,x))\biggr{)}\exp\biggl{(}\int_{0}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\biggr{)}dx$
$\displaystyle\leqslant$
$\displaystyle|O_{t}|\|v^{0}_{\varepsilon_{2}}\|_{W^{1,\infty}(\mathbb{R}^{d})}e^{\gamma
T}|X^{-1}_{v_{\varepsilon_{1}}}(t,x)-X^{-1}_{v_{\varepsilon_{2}}}(t,x)|$
$\displaystyle\leqslant$
$\displaystyle\widetilde{C}_{T}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1},$
(32)

where we have used (24) on the second line111This term is responsible for the
possible loss of regularity for $t>0$: In order to prove that
$v_{\varepsilon}$ is a Cauchy sequence in
$\mathcal{C}([0,T],W^{1,1}(\mathbb{R}^{d}))$, we would need a
$W^{2,\infty}(\mathbb{R}^{d})$ estimate over $v^{0}_{\varepsilon}$, which we
do not have..
From the definition of $\mathcal{G}_{u}(t,x)$, we observe that

$\displaystyle|\Delta E(s,t,x)|\leqslant$ $\displaystyle e^{\gamma
T}\int_{s}^{t}|\mathcal{G}_{1}(\tau,t,x)-\mathcal{G}_{2}(\tau,t,x)|d\tau$
$\displaystyle\leqslant$
$\displaystyle\widetilde{C}_{T}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1},$
(33)

where $\widetilde{C}_{T}$ depends on $\overline{\rho}$, $T$, the derivatives
of $R$, $a$ and $\psi_{a}$, and on the constant appearing in (24).
Using the change of variables
$y=X_{v_{\varepsilon_{2}}}(s,X^{-1}_{v_{\varepsilon_{2}}}(t,x))$ and recalling
that

$\displaystyle|J_{X_{v_{\varepsilon_{2}}}(t,X^{-1}_{v_{\varepsilon_{2}}}(s,x))}|$
$\displaystyle=e^{\int_{0}^{t}\text{div
}\mathcal{A}_{u}(\tau,X_{v_{\varepsilon_{2}}}(t,X^{-1}_{v_{\varepsilon_{2}}}(s,x)))d\tau-\int_{0}^{s}\text{div
}\mathcal{A}_{v_{\varepsilon_{2}}}(\tau,x)d\tau},$

we see that

$\displaystyle\int\limits_{\mathbb{R}^{d}}\mathcal{M}_{2}(s,t,x,z)dx$
$\displaystyle=\int\limits_{\mathbb{R}^{d}}m(s,y,z,(I_{d}v_{\varepsilon_{2}})(s,y))|J_{X_{v_{\varepsilon_{2}}}(t,X^{-1}_{v_{\varepsilon_{2}}}(s,y))}|^{-1}dy$
$\displaystyle\leqslant
e^{2\tilde{a}T}\int\limits_{\mathbb{R}^{d}}\sup\limits_{s,z,I}m_{\varepsilon_{i}}(s,y,z,I)dy\leqslant
e^{2\tilde{a}T}|\mathcal{K}|\overline{M}.$

Therefore, we have the bounds

$\displaystyle\int\limits_{\mathbb{R}^{d}}v^{0}(X^{-1}_{v_{\varepsilon_{2}}}(t,x))|\Delta
E(0,t,x)|dx$ $\displaystyle\leqslant
e^{\tilde{\alpha}T}\|v^{0}\|_{L^{1}(\mathbb{R}^{d})}\widetilde{C}_{T}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1},$
(34)
$\displaystyle\int\limits_{0}^{t}\int\limits_{\mathbb{R}^{d}}\int\limits_{\mathbb{R}^{d}}\mathcal{M}_{2}(s,t,x,z)\biggl{(}\Delta
V_{\varepsilon_{1}\varepsilon_{2}}(s,z)\biggr{)}dz\exp\left(\int\limits_{s}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\right)dxds$
$\displaystyle\leqslant
e^{(\gamma+2\tilde{a})T}|\mathcal{K}|\overline{M}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1},$
(35)
$\displaystyle\int\limits_{0}^{t}\int\limits_{\mathbb{R}^{d}}\int\limits_{\mathbb{R}^{d}}\mathcal{M}_{2}(t,x,z)v_{\varepsilon_{j}}(s,z)dz\biggl{(}|\Delta
E(s,t,x)|\biggr{)}dxds$ $\displaystyle\leqslant
e^{2\tilde{a}T}|\mathcal{K}|\overline{M}\overline{\rho}T\widetilde{C}_{T}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1}.$
(36)

The function $m(t,x,z,I)$ having a compact support on the $x$ variable, leads
to

$\mathcal{M}_{1}(s,t,x,z)-\mathcal{M}_{2}(s,t,x,z)=0,\mbox{ for all
}x\not\in\mathcal{K}+B_{2r_{t}}.$

On the other hand, the function $m$ being differentiable and Lipschitz,
implies that, for all $x\in\mathcal{K}+B_{2r_{t}}$

$\displaystyle|\mathcal{M}_{1}(s,t,x,z)-\mathcal{M}_{2}(s,t,x,z)|\leqslant$
$\displaystyle\|m\|_{W^{1,\infty}_{x}}|X_{v_{\varepsilon_{1}}}(s,X^{-1}_{v_{\varepsilon_{1}}}(t,x))-X_{v_{\varepsilon_{2}}}(s,X^{-1}_{v_{\varepsilon_{2}}}(t,x))|$
$\displaystyle+\mu|(I_{d}v_{\varepsilon_{1}})(s,X_{v_{\varepsilon_{1}}}(s,X^{-1}_{v_{\varepsilon_{1}}}(t,x)))-(I_{d}v_{\varepsilon_{2}})(s,X_{v_{\varepsilon_{2}}}(s,X^{-1}_{v_{\varepsilon_{2}}}(t,x)))|$
$\displaystyle\leqslant$
$\displaystyle(\|m\|_{W^{1,\infty}_{x}}+\|\psi_{d}\|_{W^{1,\infty}_{x}}\overline{\rho})|X_{v_{\varepsilon_{1}}}(s,X^{-1}_{v_{\varepsilon_{1}}}(t,x))-X_{v_{\varepsilon_{2}}}(s,X^{-1}_{v_{\varepsilon_{2}}}(t,x))|$
$\displaystyle+\|\psi_{d}\|_{L^{\infty}}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}.$

Using first Lemma 1 and then Lemma 4, we conclude that there exists a constant
$\widetilde{C}_{T}$ such that

$\displaystyle|\mathcal{M}_{1}(s,t,x,z)-\mathcal{M}_{2}(s,t,x,z)|\leqslant$
$\displaystyle\widetilde{C}_{T}\biggl{(}|X^{-1}_{v_{\varepsilon_{1}}}(t,x)-X^{-1}_{v_{\varepsilon_{2}}}(t,x)|+\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1}+\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}\biggr{)}$
$\displaystyle\leqslant$
$\displaystyle\widetilde{C}_{T}\biggl{(}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1}+\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}\biggr{)}.$

Therefore, we have the bound

$\displaystyle\int\limits_{0}^{t}\int\limits_{\mathbb{R}^{d}}\int\limits_{\mathbb{R}^{d}}\biggl{(}\mathcal{M}_{1}(s,t,x,z)-\mathcal{M}_{2}(s,t,x,z)\biggr{)}v_{\varepsilon_{1}}(s,z)dz\exp\left(\int_{s}^{t}\mathcal{G}_{1}(\tau,t,x)d\tau\right)dxds$
$\displaystyle\leqslant$
$\displaystyle|\mathcal{K}+B_{2r_{T}}|\overline{\rho}\widetilde{C}_{T}e^{\gamma
T}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{1}.$ (37)

Putting together the bounds (31) through (37), we get

$\displaystyle\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}\leqslant\widetilde{C}_{T}\left(\|v^{0}_{\varepsilon_{1}}-v^{0}_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R})}+\int_{0}^{T}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}ds\right).$

Thanks to Grönwall’s lemma, we have then the relation

$\sup\limits_{t\in[0,T]}\|v_{\varepsilon_{1}}-v_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}\leqslant\widetilde{C}_{T}\|v^{0}_{\varepsilon_{1}}-v^{0}_{\varepsilon_{2}}\|_{L^{1}(\mathbb{R}^{d})}$

for some $\widetilde{C}_{T}$ independent of $\varepsilon_{1}$ and
$\varepsilon_{2}$, which proves that, up to the extraction of a sub-sequence,
$v_{\varepsilon}$ is a Cauchy sequence in
$\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$. Therefore, there exists
$v\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$ such that

$\lim\limits_{\varepsilon\rightarrow
0}\sup\limits_{t\in[0,T]}\|v_{\varepsilon}-v\|_{L^{1}(\mathbb{R}^{d})}=0.$

Furthermore, such function satisfies the bounds (25) and (27).
We claim now that the sequence $L^{*}_{v_{\varepsilon}}\varphi$ converges to
$L^{*}_{v}\varphi$ in $L^{\infty}([0,T]\times\mathbb{R}^{d})$ for all
$\varphi\in\mathcal{C}^{1}_{c}([0,T)\times\mathbb{R}^{d})$. This is a direct
consequence of the relation

$|L^{*}v_{\varepsilon}\varphi-L^{*}v_{\varepsilon}\varphi|\leqslant(L_{r}\|\psi_{g}\|_{L^{\infty}}+\mu\|\psi_{d}\|_{L^{\infty}})\|\varphi\|_{L^{\infty}(\mathbb{R}^{d})}\|v^{\varepsilon}-v\|_{L^{1}(\mathbb{R}^{d})}.$

In order to conclude, we recall that all classical solutions are weak
solutions, and therefore, for all $\varepsilon>0$

$\int_{0}^{T}\int_{\mathbb{R}^{d}}v_{\varepsilon}L^{*}_{v_{\varepsilon}}\varphi
dxdt=\int_{\mathbb{R}^{d}}v^{0}_{\varepsilon}\varphi dx,$

and taking the limit when $\varepsilon$ goes to $0$ we see that $v$ is a weak
solution of problem (19). ∎

Considering initial data with compact support might be enough in order to
model most of the biological scenarios found in nature. However, the
hypothesis $v^{0}\in W^{k,\infty}(\mathbb{R}^{d})$ might be too restrictive
for some real life scenarios. Furthermore, the study of the problem when more
general initial conditions are present, is of theoretical interest. We show
below that, when $\partial_{I}a=0$, a solution exists for any initial data
$v^{0}\in W^{k,1}(\mathbb{R}^{d})$, $k\geqslant 1$.

######  Theorem 3.

Under hypothesis (15) through (18), if $\partial_{I}a=0$, for all non-negative
functions $v^{0}\in W^{k,1}(\mathbb{R}^{d})$, there exists a unique non-
negative weak solution $v\in\mathcal{C}([0,T],W^{k,1}(\mathbb{R}^{d}))$ of
problem (19). Furthermore, such a solution satisfies

$\displaystyle\sup\limits_{t\in[0,T]}\|v\|_{L^{1}(\mathbb{R}^{d})}\leqslant\max\\{\|v^{0}\|_{L^{1}(\mathbb{R}^{d})},\frac{I^{*}}{\underline{\psi_{g}}}\\},$
(38)
$\displaystyle\sup\limits_{t\in[0,T]}\|v\|_{W^{k,\infty}(\mathbb{R}^{d})}\leqslant
C_{T}\|v^{0}\|_{W^{k,1}(\mathbb{R}^{d})}.$ (39)

###### Proof.

As in the proof for $k=1$ when $\partial_{I}a\neq 0$, we can approximate any
function $v^{0}\in W^{k,1}$ by a smooth, compactly supported sequence of
functions $v^{0}_{\varepsilon}$. The same arguments as in the previous proof
will show that $v_{\varepsilon}$, the sequence of solutions associated to
$v^{0}_{\varepsilon}$, is a Cauchy sequence in
$\mathcal{C}([0,T],L^{1}(\mathbb{R}^{d}))$. Furthermore, given that the second
term in (30), which is responsible for the possible loss of regularity in the
previous case, is equal $0$ when $\partial_{I}a=0$, we show that
$v_{\varepsilon}$ is a Cauchy sequence in
$\mathcal{C}([0,T],W^{k,1}(\mathbb{R}^{d}))$ as well. We prove as in Theorem 2
that the limit of $v_{\varepsilon}$ is the required weak solution. ∎

Given that the regularity of the solution varies depending on whether
$\partial_{I}a=0$ or $\partial_{I}a\neq 0$, and that such regularity will be
of importance in the upcoming sections, we define the parameter

$\kappa:=\left\\{\begin{matrix}k-1,&\mbox{ if }\partial_{I}a\neq 0,\\\ &\\\
k,&\mbox{ if }\partial_{I}a=0,\end{matrix}\right.$

which encompasses the information over said regularity.

We remark that if we had $\partial_{I}m=0$, we might obtain existence for a
larger class of mutation functions $m$. For conciseness, we will however not
consider such cases in the present work.

##  3 Particle Method

The particle method basically consists in searching for an approximate
solution of problem (19) which is a sum of weighted Dirac masses.
Throughout the following section we suppose

$\displaystyle\psi_{a}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}_{x}^{d},W^{1,\infty}(\mathbb{R}^{d})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{y},\mathcal{C}^{2}(\mathbb{R}^{d}_{x})\cap
W^{2,\infty}(\mathbb{R}^{d}_{x})\right),$ (40) $\displaystyle
0<\underline{\psi_{g}}\leqslant\psi_{g}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}_{x}^{d},W^{1,\infty}(\mathbb{R}_{y}^{d})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{y}^{d},\mathcal{C}^{1}(\mathbb{R}_{x}^{d})\cap
W^{1,\infty}(\mathbb{R}_{x}^{d})\right)$ (41) $\displaystyle\psi_{d}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x},W^{1,\infty}(\mathbb{R}_{y}^{d})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y},\mathcal{C}^{1}(\mathbb{R}_{x}^{d})\cap
W^{1,\infty}(\mathbb{R}_{x}^{d}).\right)$ (42)

Notice that, unlike the set of hypotheses (15)-(17), we have imposed
$W^{1,\infty}(\mathbb{R}^{d})$ regularity for the $y$ variable, which is
needed in order to approximate the integral terms by sums over a countable
set.
Consider as well

$0\leqslant
m\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x}\times\mathbb{R}^{d}_{y},W^{1,\infty}(\mathbb{R}_{I})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x}\times\mathbb{R}_{I},W^{1,\infty}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y}\times\mathbb{R}_{I},\mathcal{C}^{1}_{c}(\mathbb{R}_{x}^{d})\right)$
(43)

satisfying hypotheses (7) through (18).
For $h>0$, consider a countable set of indices
$\mathcal{J}_{h}\in\mathbb{Z}^{d}$, points $x^{0}_{i}\in\mathbb{R}^{d}$ and
weights $w^{0}_{i}$ for $i\in\mathcal{J}_{h}$. The weights $w^{0}_{i}$ can be
regarded as the respective masses of a collection of subsets
$\Omega^{0}_{i}\subset\mathbb{R}^{d}$ satisfying

$\Omega^{0}_{i}\cap\Omega^{0}_{j}=\emptyset,\mbox{ if }i\neq j,\mbox{ and
}\bigcup\limits_{i\in\mathcal{J}_{h}}\Omega^{0}_{i}=\mathbb{R}^{d}.$ (44)

For example, we may choose the $\Omega^{0}_{i}$ as the set of all non
intersecting cubes of side length equal $h$ having the points $hi$ as centers
, $i\in\mathbb{Z}^{d}$. This way, $w^{0}_{i}=h^{d}$, with each of the
$x^{0}_{i}$ being a point in $\Omega^{0}_{j}$. In general we assume that there
exist positive constants $c$ and $C$ such that

$\displaystyle ch\leqslant|x^{0}_{i}-x^{0}_{j}|\leqslant Ch,\ \forall i\neq
j,$ (45) $\displaystyle ch^{d}\leqslant w^{0}_{i}\leqslant Ch^{d},\ \forall
i\in\mathcal{J}_{h}.$ (46)

Following [20], the particle method then consists in looking for a measure
$\nu_{h}$ of the form

$\nu_{h}(t)=\sum\limits_{i\in\mathcal{J}_{h}}\nu_{i}(t)w_{i}(t)\delta_{x_{i}(t)},$

where
$(\nu:=\\{v_{i}(t)\\}_{i\in\mathcal{J}_{h}},w:=\\{w_{i}(t)\\}_{i\in\mathcal{J}_{h}},\overline{x}:=\\{x_{i}(t)\\}_{i\in\mathcal{J}_{h}})$,
is the solution of the following system

$\left\\{\begin{matrix}[l]\dot{x}_{i}(t)&=&A_{\nu,w}(t,x_{i}),\\\ &&\\\
\dot{w}_{i}(t)&=&\text{div }A_{\nu,w}(t,x_{i}(t))w_{i}(t),\\\ &&\\\
\dot{\nu}_{i}(t)&=&\big{(}-\text{div
}A_{\nu,w}(t,x_{i}(t))+R(t,x_{i}(t),I_{g}(t,x_{i}(t),\nu,w))\big{)}\nu_{i}(t)\\\
&&\\\
&&+\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)\nu_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t),\nu,w)),\\\
&&\\\ x_{i}(0)&=&x^{0}_{i},\ w_{i}(0)\ =\ w^{0}_{i},\ \nu_{i}(0)\ =\
v^{0}(x^{0}_{i}),\end{matrix}\right.$ (47)

where

$A_{\nu,w}(t,x)=a(t,x,I_{a}(t,x,\nu,w)),$

and

$\displaystyle I_{l}(t,x,\nu,w)$
$\displaystyle:=\sum\limits_{j\in\mathcal{J}_{h}}\nu_{j}(t)w_{j}(t)\psi_{l}(t,x,x_{j}(t)),$

with $l\in\\{a,g,d\\}$.
In what follows we assume that $h$ and $x^{0}_{k}$ are chosen in such a way
that

$\|v^{0}\|_{1,h}:=\sum\limits_{i\in\mathcal{J}_{h}}v^{0}(x^{0}_{i})w^{0}_{i}<\infty.$
(48)

We define the subset of indices

$\mathcal{J}^{m}_{h}:=\\{i\in\mathcal{J}_{h}:x^{0}_{i}\in\mbox{supp
}m+B_{\|a\|_{L^{\infty}}T}\\},$

where $B_{\|a\|_{L^{\infty}}T}$ is the ball of radius $\|a\|_{L^{\infty}}T$.
The compact support of $m$ implies that $|\mathcal{J}^{m}_{h}|<\infty$.

For a positive value of $h$ and a set of indexes $\mathcal{J}_{h}$ we define
the functional spaces

$\displaystyle\mathcal{l}^{1}(\mathcal{J}_{h})$
$\displaystyle:=\big{\\{}u=\\{u_{i}\\}_{i\in\mathcal{J}_{h}}:\sum\limits_{i\in\mathcal{J}_{h}}|u_{i}|<+\infty\big{\\}},$
$\displaystyle\mathcal{l}^{\infty}(\mathcal{J}_{h})$
$\displaystyle:=\big{\\{}w=\\{w_{i}\\}_{i\in\mathcal{J}_{h}}:\sup\limits_{i\in\mathcal{J}_{h}}|w_{i}|<+\infty\big{\\}}.$

We equip these spaces with the norms

$\|u\|_{\mathcal{l}^{1}}:=\sum\limits_{i\in\mathcal{J}_{h}}|u_{i}|\mbox{ and
}\|w\|_{\mathcal{l}^{\infty}}:=\sup\limits_{i\in\mathcal{J}_{h}}|w_{i}|$

respectively. It is clear that for all $u\in\mathcal{l}^{1}(\mathcal{J}_{h})$
and $v\in\mathcal{l}^{\infty}(\mathcal{J}_{h})$, then
$uv\in\mathcal{l}^{1}(\mathcal{J}_{h})$. For $T>0$ we define as well the
spaces

$X^{T}_{h}:=\mathcal{C}([0,T],\mathcal{l}^{1}(\mathcal{J}_{h}))\mbox{ and
}Y^{T}_{h}:=\mathcal{C}([0,T],\mathcal{l}^{\infty}(\mathcal{J}_{h})),$

equipped with the norms

$\|u\|_{1,h}:=\sup\limits_{t\in[0,T]}\|u(t)\|_{\mathcal{l}^{1}}\mbox{ and
}\|w\|_{\infty,h}:=\sup\limits_{t\in[0,T]}\|w(t)\|_{\mathcal{l}^{\infty}}.$

Problem (47) is a strongly coupled system of ODEs, with an infinite number of
unknowns and equations. In some cases the system becomes uncoupled (for
example if $\partial_{I}a=0$) or with a finite number of equations and
unknowns (for example if $v^{0}$, $a$ and $m$ have compact support), however,
for the sake of generality, we present below the proof of existence of
solution in the general case, and later discuss briefly these particular
scenarios.
We start by giving two results that will be of great use for the proof of
existence for problem (47). First, we deal with the existence of solution for
a simpler system of infinite equations with infinitely many unknowns:

######  Lemma 7.

Consider $a\in\mathcal{C}([0,T],(W^{1,\infty}(\mathbb{R}^{d+1})))^{d}$, $u\in
X^{T}_{h}$, $w\in Y^{T}_{h}$ and $\psi_{a}$ satisfying hypothesis (40). Then
there exists a unique family of functions
$x:=\\{x_{i}\\}_{i\in\mathcal{J}_{h}}$, $x_{i}\in\mathcal{C}^{1}([0,T])$ for
all $i\in\mathcal{J}_{h}$ which is solution of the system of equations

$\dot{x_{i}}(t)=A_{u,w}(t,x_{i}),\quad t\in[0,T],\ x_{i}(0)=x^{0}_{i}.$ (49)

When $\partial_{I}a=0$, system (49) becomes uncoupled, each individual
equation has a solution, thanks to the classic Cauchy-Lipschitz theory. The
proof of the general case is given in Appendix B.
The second auxiliary result comes from approximation theory, and it will also
be of great use in Section 4:

######  Lemma 8.

$\displaystyle\forall\varphi\in
W^{k,1}(\mathbb{R}^{d}),\quad\bigg{\lvert}\int_{\mathbb{R}^{d}}{\varphi(x)dx}-\underset{i\in\mathcal{J}_{h}}{\sum}{w_{i}(t)\varphi(x_{i}(t))}\bigg{\rvert}\leqslant
Ch^{k}\lVert\varphi\rVert_{k,1},$

where $C$ is a constant which depends on $a$, $\psi_{a}$, $\|vw\|_{1,h}$ and
$T$.

This result is a direct corollary of Lemma 8 in [42]. More details regarding
its proof are given in Appendix C.
From now on, we suppose $h$ to be small enough so that for any $t\in[0,T]$,

$\underset{i\in\mathcal{J}_{h}}{\sum}{w_{i}(t)m(t,x_{i}(t),y,I_{d}(t,x_{i}(t),\nu,w))}<K+\frac{r^{*}}{2},$
(50)

where the values of $K$ and $r^{*}$ are given in (11) and (18) respectively.
Such a choice is always possible thanks to Lemma 8.

######  Theorem 4.

Under hypothesis (5) through (18) and (40) through (46), for all $T>0$ and all
non-negative initial data
$v^{0}\in\mathcal{l}^{1}(\mathcal{J}_{h},\Omega^{0})$ there exists a unique
solution $x_{i}\in\mathcal{C}^{1}([0,T])$, for all $i\in\mathcal{J}_{h}$,
$w:=\\{w_{i}(\cdot)\\}_{i\in\mathcal{J}_{h}}\in\mathcal{C}([0,T],\mathcal{l}^{\infty}(\mathcal{J}_{h}))$
and
$0\leqslant\nu:=\\{\nu_{i}(\cdot)\\}_{i\in\mathcal{J}_{h}}\in\mathcal{C}([0,T],\mathcal{l}^{1}(\mathcal{J}_{h}))$
of problem (47). Furthermore, there exist positive constants $c_{T}$ and
$C_{T}$ such that the solution satisfies, for all $t\in[0,T]$

$\displaystyle c_{T}h\leqslant|x_{i}(t)-x_{j}(t)|\leqslant C_{T}h,\ \forall
i,j\in\mathcal{J}_{h},\ i\neq j,$ (51) $\displaystyle c_{T}h^{d}\leqslant
w_{i}(t)\leqslant C_{T}h^{d},\ \forall i\in\mathcal{J}_{h},$ (52)
$\displaystyle\|\nu
w\|_{1,h}\leqslant\max\\{\|v^{0}h^{d}\|_{\mathcal{l}^{1}},\frac{I^{*}}{\underline{\psi_{g}}}\\}.$
(53)

###### Proof.

Consider $v^{0}\in\mathcal{l}^{1}(\mathcal{J}_{h},\Omega^{0})$, satisfying
$v^{0}\geqslant 0$. Consider as well $\alpha>1$, and define

$\displaystyle\overline{\rho}_{\alpha}$ $\displaystyle:=\max\\{\alpha
h^{d}\|v^{0}\|_{\mathcal{l}^{1}},\frac{I^{*}}{\underline{\psi_{g}}}\\},$
$\displaystyle\tilde{a}$
$\displaystyle:=\|a\|_{W^{1,\infty}_{x}L^{\infty}_{I,t}}+\|a\|_{W^{1,\infty}_{I}L^{\infty}_{t,x}}\|\psi_{a}\|_{W^{1,\infty}_{x}}C_{\alpha}.$

For $T>0$ we define the set

$D^{T}_{\alpha}:=\\{(u,w)\in X^{T}_{h}\times
Y^{T}_{h}:\|uw\|_{1,h}\leqslant\overline{\rho}_{\alpha},\forall
t\in[0,T],u(t)\geqslant 0,w(t)\geqslant 0,h^{d}e^{-\tilde{a}t}\leqslant
w_{k}(t)\leqslant h^{d}e^{\tilde{a}t}\\}.$

For any $(u,w)\in D^{T}_{\alpha}$ we introduce the problem, for $t\in[0,T]$,

$\left\\{\begin{matrix}[l]\dot{x}_{i}(t)&=&A_{u,w}(t,x_{i}),\\\ &&\\\
\dot{\omega}_{i}(t)&=&\text{div }A_{u,w}(t,x_{i}(t))\omega_{i}(t),\\\ &&\\\
\dot{\nu}_{i}(t)&=&\big{(}-\text{div
}A_{u,w}(t,x_{i}(t))+R(t,x_{i}(t),I_{g}(t,x_{i}(t),u,w))\big{)}\nu_{i}(t)\\\
&&\\\
&&+\sum\limits_{j\in\mathcal{J}_{h}}\omega_{j}(t)u_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t),u,w)),\\\
&&\\\ x_{i}(0)&=&x^{0}_{i},\ \omega_{i}(0)\ =\ w^{0}_{i},\ \nu_{i}(0)\ =\
v^{0}(x^{0}_{i}).\end{matrix}\right.$ (54)

We denote $(\nu,\omega)=\Phi(u,w)$.

For each pair $(u,w)$, the existence and uniqueness of $x_{i}$ is immediate
from Lemma 7. Furthermore, for all values of $i$, we have the following
explicit expression for $\omega_{i}$

$\omega_{i}(t)=w^{0}_{i}e^{\int_{0}^{t}\text{div }A_{u,w}(s,x_{i}(s))ds},$

which satisfies, for any $t\in[0,T]$

$h^{d}e^{-\tilde{a}t}\leqslant\omega_{i}(t)\leqslant h^{d}e^{\tilde{a}t}.$

On the other hand, for all $(u,w)\in D^{T}_{\alpha}$ and all values of $i$,
the right-hand side of the differential equation in (54) is well defined, as
we have for all $t\in[0,T]$

$\displaystyle\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)u_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t),u,w))$
$\displaystyle\leqslant\overline{M}\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)u_{j}(t)\leqslant\overline{M}\overline{\rho}_{\alpha}.$

Therefore, the expression for $\nu_{i}$ is given by

$\nu_{i}(t)=v^{0}(x^{0}_{i})e^{\int_{0}^{t}\mathcal{G}_{i}(s)ds}+\int\limits_{0}^{t}\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(s)u_{j}(s)m(s,x_{i}(s),x_{j}(s),I_{d}(s,x_{i}(s),u,w))e^{\int_{s}^{t}\mathcal{G}_{i}(\tau)d\tau}ds,$
(55)

where

$\mathcal{G}_{i}(t):=-\text{div
}A_{u,w}(s,x_{i}(s))+R(s,x_{i}(s),I_{g}(s,x_{i}(s),u(s),w(s))),$

satisfies

$\underset{t\in[0,T]}{\sup}\,|\mathcal{G}_{i}(t)|\leqslant\gamma:=\|R\|_{L^{\infty}_{t,x,I}}+\|a\|_{W^{1,\infty}_{x}L^{\infty}_{t,I}}+\|a\|_{W^{1,\infty}_{I}L^{\infty}_{t,x}}\|\psi_{a}\|_{W^{1,\infty}_{x}L^{\infty}_{y}}\overline{\rho}_{\alpha}.$

The positiveness of $\nu$ is immediate from the positiveness of $v^{0}$ and
$m$.

Furthermore, given that $|x_{i}(t)-x^{0}_{i}|\leqslant\|a\|_{L^{\infty}}T$,
for all $k\in\mathcal{J}_{h}$ and $t\in[0,T]$, we have

$m(t,x_{i}(s),y,I)=0,$

for all $i\not\in\mathcal{J}^{m}_{h}$, $t\in[0,T]$, $y\in\mathbb{R}^{d}$ and
$I\in\mathbb{R}$. As a result, multiplying (55) by $\omega_{i}(t)$ for each
$i$ and adding for all values of $i$, we obtain

$\displaystyle\sum\limits_{i\in\mathcal{J}_{h}}\nu_{i}(t)\omega_{i}(t)$
$\displaystyle\leqslant
e^{(\gamma+\tilde{\alpha})T}\left(\sum\limits_{i\in\mathcal{J}_{h}}v^{0}(x^{0}_{i})h^{d}+h^{d}\int\limits_{0}^{t}\sum\limits_{i\in\mathcal{J}_{h}}\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(s)u_{j}(s)m(s,x_{i}(s),x_{j}(s),I_{d}(s,x_{i}(s),u,w))ds\right)$
$\displaystyle=e^{(\gamma+\tilde{\alpha})T}\left(\sum\limits_{i\in\mathcal{J}_{h}}v^{0}(x^{0}_{i})h^{d}+h^{d}\int\limits_{0}^{t}\sum\limits_{i\in\mathcal{J}^{m}_{h}}\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(s)u_{j}(s)m(s,x_{i}(s),x_{j}(s),I_{d}(s,x_{i}(s),u,w))ds\right)$
$\displaystyle\leqslant
e^{(\gamma+\tilde{\alpha})T}\left(\sum\limits_{i\in\mathcal{J}_{h}}v^{0}(x^{0}_{i})h^{d}+h^{d}|\mathcal{J}^{m}_{h}|\overline{M}\int\limits_{0}^{t}\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(s)u_{j}(s)ds\right)$
$\displaystyle\leqslant
e^{(\gamma+\tilde{\alpha})T}(\frac{1}{\alpha}+TK_{h})\overline{\rho}_{\alpha},$

where $K_{h}:=h^{d}|\mathcal{J}^{m}_{h}|\overline{M}$222Notice that
$h^{d}|\mathcal{J}^{m}_{h}|\approx|\mbox{supp }m+B_{\|a\|_{L^{\infty}}T}|$.
Thanks to the condition $\alpha>1$ there exists $T_{\alpha}$ (only depending
on $\alpha$ and on the coefficients of the problem) such that
$\Phi:D^{T}_{\alpha}\rightarrow D^{T}_{\alpha}$, for all $T\leqslant
T_{\alpha}$.
We now prove that there exists $T\in(0,T_{\alpha})$ such that $\Phi$ is a
contraction over $D^{T}_{\alpha}$.
Step 1: Bounds over $x=\\{x_{i}\\}_{i\in\mathcal{J}_{h}}$
Let $(u^{1},w^{1})$ and $(u^{2},w^{2})$ be two pairs in $D^{T}_{\alpha}$, and
let $x^{1},x^{2}$ be the respective solutions of

$\dot{x^{j}_{i}}=A_{u^{j},w^{j}}(t,x^{j}_{i}).$

By following the same ideas as in the proof of Lemma 1 (see Appendix A), we
obtain that for all $t\in[0,T]$,

$\|x^{1}(t)-x^{2}(t)\|_{\infty,h}\leqslant
C(T,h)\left(\|w^{1}-w^{2}\|_{\infty,h}+\|u^{1}-u^{2}\|_{1,h}\right),$

where the constant $C(T,h)$ satisfies $\lim\limits_{T\rightarrow 0}C(T,h)=0$.
Step 2: Bounds over $\omega=\\{\omega_{i}\\}_{i\in\mathcal{J}_{h}}$
From the expression for $\omega$, we get, for all $t\in[0,T]$

$\displaystyle|\omega^{1}_{i}(t)-\omega^{2}_{i}(t)|\leqslant$ $\displaystyle
h^{d}Te^{\tilde{a}T}\left(\|a\|_{W^{2,\infty}_{x,I}}(1+|\partial_{x}I_{a}(t,x^{2}_{i},u^{2},w^{2})|)\|x^{1}-x^{2}\|_{\infty,h}\right.$
$\displaystyle+\|a\|_{W^{2,\infty}_{x,I}}(1+|\partial_{x}I_{a}(t,x^{2}_{i},u^{2},w^{2})|)|I_{a}(t,x^{1}_{i},u^{1},w^{1})-I_{a}(t,x^{2}_{i},u^{2},w^{2})|$
$\displaystyle\left.+\|a\|_{W^{1,\infty}_{I}}|\partial_{x}I_{a}(t,x^{1}_{i},u^{1},w^{1})-\partial_{x}I_{a}(t,x^{2}_{i},u^{2},w^{2})|\right).$

On the other hand we have

$\displaystyle|\partial_{x}I_{a}(t,x^{2}_{i},u^{2},w^{2})|\leqslant$
$\displaystyle\|\psi_{a}\|_{W^{1,\infty}_{x}}\overline{\rho}_{\alpha},$
$\displaystyle|I_{a}(t,x^{1}_{i},u^{1},w^{1})-I_{a}(t,x^{2}_{i},u^{2},w^{2})|\leqslant$
$\displaystyle\overline{\rho}_{\alpha}\|\psi_{a}\|_{W^{1,\infty}_{x,y}}\|x^{1}-x^{2}\|_{\infty,h},$
$\displaystyle+\|\psi_{a}\|_{L^{\infty}_{x,y}}e^{\tilde{a}T}(h^{d}\|u^{1}-u^{2}\|_{1,\infty}+\frac{\overline{\rho}_{\alpha}}{h^{d}}\|w^{1}-w^{2}\|_{\infty,h}),$

and

$\displaystyle|\partial_{x}I_{a}(t,x^{1}_{i},u^{1},w^{1})-\partial_{x}I_{a}(t,x^{2}_{i},u^{2},w^{2})|\leqslant$
$\displaystyle\overline{\rho}_{\alpha}\|\psi_{a}\|_{W^{2,\infty}_{x,y}}\|x^{1}-x^{2}\|_{\infty,h}$
$\displaystyle+\|\psi_{a}\|_{W^{1,\infty}_{x}}e^{\tilde{a}T}(h^{d}\|u^{1}-u^{2}\|_{1,\infty}+\frac{\overline{\rho}_{\alpha}}{h^{d}}\|w^{1}-w^{2}\|_{\infty,h}).$

In conclusion, there exists a constant $C(T,h)$, satisfying
$\lim\limits_{T\rightarrow 0}C(T,h)=0$ such that

$\|\omega^{1}-\omega^{2}\|_{\infty,h}\leqslant
C(T,h)\left(\|w^{1}-w^{2}\|_{\infty,h}+\|u^{1}-u^{2}\|_{1,h}\right).$

Step 3: Bounds over $\nu=\\{\nu_{i}\\}_{i\in\mathcal{J}_{h}}$
Using the expression for $\nu$, the regularity of $m$ and bounds similar to
those used for $\omega$, we see that there exists a constant $C(T,h)$
satisfying $\lim\limits_{T\rightarrow 0}C(T,h)=0$ such that

$\displaystyle\|\nu^{1}-\nu^{2}\|_{1,h}\leqslant
C(T,h)\left(\|w^{1}-w^{2}\|_{\infty,h}+\|u^{1}-u^{2}\|_{1,h}\right).$

Consequently, there exists a constant $C(T,h)$ satisfying
$\lim\limits_{T\rightarrow 0}C(T,h)=0$, such that

$\|\omega^{1}-\omega^{2}\|_{\infty,h}+\|\nu^{1}-\nu^{2}\|_{1,h}\leqslant
C(T,h)\left(\|w^{1}-w^{2}\|_{\infty,h}+\|u^{1}-u^{2}\|_{1,h}\right),$

which implies that, for $0<T_{1}\leqslant T_{\alpha}$ small enough, $\Phi$ is
a contraction over $D^{T_{1}}_{\alpha}$, and therefore it has a unique fixed
point. Such fixed point is a solution of problem (47) over $[0,T_{1})$.
We now claim that the solution exists for $T$ arbitrary, and furthermore, it
satisfies the relation (53). Let $T_{f}$ be the maximal time of existence of
solution. Suppose that there exists $t_{0}\in(0,T_{f}]$ such that $\|\nu
w\|_{1,h}>\overline{\rho}_{\alpha}$. This implies that there exist
$\delta\geqslant 0$ and $t^{*}>0$ such that for a certain finite subset of
$\mathcal{J}_{h}$, that we denote as $\mathcal{K}_{h}$, the following
statements are true:

$\displaystyle\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}(t)w_{i}(t)\leqslant\overline{\rho}_{\alpha},\
\forall\ t\in[t^{*}-\delta,t^{*}],$
$\displaystyle\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}(t)w_{i}(t)>\overline{\rho}_{\alpha},\
\forall\ t\in(t^{*},t^{*}+\delta].$

This implies the existence of $t_{1}\in[t^{*},t^{*}+\delta]$ such that the
following properties are satisfied simultaneously

$\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}(t_{1})w_{i}(t_{1})\geqslant\overline{\rho}_{\alpha}\mbox{
and
}\left(\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}w_{i}\right)^{\prime}(t_{1})\geqslant
0.$ (56)

Multiplying the equation satisfied by $\nu_{i}(t)$ by $w_{i}(t)$ we get the
relation

$\displaystyle\dot{\nu}_{i}(t)w_{i}(t)=$ $\displaystyle\big{(}-\text{div
}A_{\nu,w}(t,x_{i}(t))+R(t,x_{i}(t),I_{g}(t,x_{i}(t),\nu,w))\big{)}\nu_{i}(t)w_{i}(t)$
$\displaystyle+w_{i}(t)\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)\nu_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{g}(t,x_{i}(t),\nu,w)),$

while, directly from the equation for $w_{i}(t)$ we deduce

$\nu_{i}(t)\dot{w}_{i}(t)=\text{div }A_{\nu,w}(t,x_{i}(t))\nu_{i}(t)w_{i}(t).$

Therefore, adding both relations for $i\in\mathcal{K}_{h}$ and using (50), we
get

$\displaystyle\left(\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}w_{i}\right)^{\prime}(t)=$
$\displaystyle\left(\sum\limits_{i\in\mathcal{K}_{h}}R(t,x_{i}(t),I_{g}(t,x_{i}(t),\nu,w))\nu_{i}(t)w_{i}(t)\right)$
$\displaystyle+\sum\limits_{j\in\mathcal{J}}w_{j}(t)\nu_{j}(t)\sum_{i\in\mathcal{K}_{h}\cap\mathcal{J}^{m}_{h}}m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t),\nu,w))w_{i}(t)$
$\displaystyle\leqslant$
$\displaystyle\left(\sum\limits_{i\in\mathcal{K}_{h}}R(t,x_{i}(t),I_{g}(t,x_{i}(t)))\nu_{i}(t)w_{i}(t)\right)+(K+\frac{r^{*}}{2})\sum\limits_{j\in\mathcal{K}_{h}}w_{j}(t)\nu_{j}(t).$

Given that
$\|\nu(t_{1})w(t_{1})\|_{\mathcal{l}^{1}}\geqslant\overline{\rho}_{\alpha}\geqslant\frac{I^{*}}{\underline{\psi}_{g}}$,
then $I_{g}(t_{1},x_{i}(t_{1}),\nu(t_{1}),w(t_{1}))\geqslant I^{*}$, and
consequently

$R(t_{1},x_{i}(t_{1}),I_{g}(t_{1},x_{i}(t_{1}),\nu(t_{1}),w(t_{1})))<-r^{*}-K,$

for all values of $i$, which in turn implies that

$\displaystyle\left(\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}w_{i}\right)^{\prime}(t_{1})\leqslant-\frac{r^{2}}{2}\sum\limits_{i\in\mathcal{K}_{h}}\nu_{i}(t_{1})w_{i}(t_{1})<0,$

which contradicts (56). Therefore, $\|\nu
w\|_{1,h}\leqslant\overline{\rho}_{\alpha}$ for all values of $\alpha>1$ and
for all $t\in(0,T_{f})$. We can then iterate the arguments used to prove
existence of a solution, and conclude that the solution can be extended to any
interval $[0,T]$. As $\|\nu w\|_{1,h}\leqslant\overline{\rho}_{\alpha}$ for
all $t$, independently of $\alpha$, taking the limit when $\alpha$ goes to
$1$, we obtain (53). ∎

##  4 Convergence of the numerical solution towards a weak solution

We study now the conditions under which a solution of problem (54) converges
towards a solution of problem (19), in a certain sense that will be defined
later. We split our analysis in two cases: first, the study of convergence on
a finite interval of time $[0,T]$. We will see that for any $T>0$, the
solution obtained through the particle method converges towards the solution
of the PDE (19). However, the speed of convergence might be affected by the
value of $T$. The second case we study is the asymptotic proximity of both
solutions when $t$ goes to $\infty$. This is a far more complex and
interesting issue, and we show different examples exposing some of the
behaviours that can be observed.

Directly from the study of existence of solutions for each problem, we notice
that the sets of hypotheses we have used, do not coincide. We give a set of
hypotheses which simultaneously guarantees the existence of solution for both
problems, while taking into account the distinction of cases involved in the
definition of $\kappa$.
For a certain $T>0$ and $k>0$, we consider the functions

$\displaystyle\psi_{a}$
$\displaystyle\in\mathcal{C}([0,T]\times\mathbb{R}_{x}^{d},W^{1,\infty}(\mathbb{R}_{y}^{d}))\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{y}^{d},\mathcal{C}^{k+1}(\mathbb{R}^{d}_{x})\cap
W^{k+1,\infty}(\mathbb{R}^{d}_{x})\right),$ (57) $\displaystyle
0<\underline{\psi_{g}}\leqslant\psi_{g}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}_{x},W^{1,\infty}(\mathbb{R}_{y}^{d})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{y},\mathcal{C}^{\kappa}(\mathbb{R}^{d}_{x})\cap
W^{\kappa,\infty}(\mathbb{R}^{d}_{x})\right),$ (58) $\displaystyle\psi_{d}$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x},W^{1,\infty}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y},\mathcal{C}^{\kappa}(\mathbb{R}^{d}_{x})\cap
W^{\kappa,\infty}(\mathbb{R}^{d}_{x})\right).$ (59)

As in Section 2 we introduce

$\displaystyle a$
$\displaystyle\in\mathcal{C}([0,T],W^{k+1,\infty}(\mathbb{R}^{d+1})),$ (60)
$\displaystyle R$
$\displaystyle\in\mathcal{C}\left([0,T]\times\mathbb{R}_{x},W^{\kappa,\infty}_{loc}(\mathbb{R}^{d}_{y})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}_{y},\mathcal{C}^{\kappa}(\mathbb{R}^{d}_{x})\cap
W^{\kappa,\infty}(\mathbb{R}^{d}_{x}))\right)$ (61)

We consider as well

$0\leqslant
m\in\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x}\times\mathbb{R}^{d}_{y},W^{\kappa,\infty}(\mathbb{R}_{I})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{x}\times\mathbb{R}^{d}_{I},W^{\kappa,\infty}(\mathbb{R}_{y}^{d})\right)\cap\mathcal{C}\left([0,T]\times\mathbb{R}^{d}_{y}\times\mathbb{R}^{d}_{I},\mathcal{C}^{\kappa}_{c}(\mathbb{R}^{d}_{x})\right).$
(62)

satisfying hypothesis (7) through (18).
Finally, we consider $v^{0}\in W^{k,1}(\mathbb{R}^{d})$ if $\partial_{I}a=0$
and $v^{0}\in W^{k,1}(\mathbb{R}^{d})\cap W^{k,\infty}(\mathbb{R}^{d})$ with
compact support otherwise.

###  4.1 Convergence on a finite time interval

We recall that the function $v$ represents the solution of problem (19) while
$x_{i}$, $w_{i}$ and $\nu_{i}$, $i\in\mathcal{J}_{h}$ represents that of
problem (54). We recall as well that

$\max\\{\|v\|,\|\nu
w\|_{1,h}\\}\leqslant\max\\{\|v^{0}\|_{L^{1}(\mathbb{R}^{d})},\|v^{0}h^{d}\|_{\mathcal{l}^{1}},\frac{I^{*}}{\underline{\psi_{g}}}\\}=:\overline{\rho}.$

Let $\varepsilon>0$, $r\in\mathbb{R}$ and
$\varphi\in\mathcal{C}_{c}(\mathbb{R}^{d})$ satisfy the following conditions

$\displaystyle\int_{\mathbb{R}^{d}}{\varphi(x)dx}$ $\displaystyle=1,$ (63)
$\displaystyle\int_{\mathbb{R}^{d}}{x^{\alpha}\varphi(x)dx}$
$\displaystyle=0,\quad\forall\;\alpha\in\mathbb{N}^{n},\lvert\alpha\rvert\leqslant
r-1.$ (64)

. We define, for all $t\in[0,T]$, $x\in\mathbb{R}^{d}$

  1. i)

$v^{h}(t,x)=\sum\limits_{i\in\mathcal{J}_{h}}{\nu_{i}(t)w_{i}(t)\delta(x-x_{i}(t))},$
(65)

a time dependent measure obtained as a sum of weighted Dirac deltas at
$x_{i}(t)$,

  2. ii)

$v_{\varepsilon}^{h}(t,x)=(v^{h}(t)*\varphi_{\varepsilon})(x)=\sum\limits_{i\in\mathcal{J}_{h}}{\nu_{i}(t)w_{i}(t)\varphi_{\varepsilon}(x-x_{i}(t))},$
(66)

a regular function obtained as the space convolution of $v^{h}(t,x)$ and
$\varphi_{\varepsilon}(x)$, where

$\varphi_{\varepsilon}:=\frac{1}{\varepsilon^{d}}\varphi(\frac{\cdot}{\varepsilon}).$

We also introduce the following operator, for any function $v\in
L^{\infty}(0,T;L^{1}(\mathbb{R}^{d}))$:

$\biggl{(}\Pi_{\varepsilon}^{h}(t)v\biggr{)}(x)=\sum\limits_{i\in\mathcal{J}_{h}}{w_{i}(t)v(t,x_{i}(t))\varphi_{\varepsilon}(x-x_{i}(t))}.$

We recall a direct corollary of the Theorem 3 in [42]:

######  Proposition 5.

Let $k$, $r$ be two integers, with $k>d$, and let us assume that $a\in
L^{\infty}\left(0,T;W^{k+1,\infty}(\mathbb{R}^{d})^{d}\right)$, and that
$\varphi\in\mathcal{C}_{c}^{1}(\mathbb{R}^{d})\cap W^{k+1,1}(\mathbb{R}^{d})$
satisfies conditions (63) and (64). Then, for any $p\in[1,+\infty]$, there
exists $C=C(T)>0$ such that, for any $u\in W^{\mu,1}(\mathbb{R}^{d})$
($\mu=\max(r,k)$),

$\displaystyle\lVert
u-{\Pi}_{\varepsilon}^{h}(t)u\rVert_{L^{p}(\mathbb{R}^{d})}\leqslant
C\big{(}\varepsilon^{r}\lVert
u\rVert_{W^{r,p}(\mathbb{R}^{d})}+\big{(}\frac{\varepsilon}{h}\big{)}^{k}\lVert
u\rVert_{W^{k,p}(\mathbb{R}^{d})}\big{)}.$

We seek to prove the following approximation result between
$v_{\varepsilon}^{h}$ and $v$, the solution of problem (19).

######  Theorem 6.

Assume that hypotheses (57) through (62) are satisfied, and that
$\varphi\in\mathcal{C}_{c}^{1}(\mathbb{R}^{d})\cap W^{k+1,1}(\mathbb{R}^{d})$
satisfies (63) and (64). Then, there exists $C=C(T,a,R,m,\overline{\rho})>0$,
a positive constant which depends on $T$, $a$, $R$, $m$ and $\overline{\rho}$
such that

$\lVert v-v_{\varepsilon}^{h}\rVert_{L^{1}(\mathbb{R}^{d})}\leqslant
C\bigl{(}\varepsilon^{r}+\big{(}\frac{h}{\varepsilon}\big{)}^{\kappa}+h^{\kappa}\bigr{)}\lVert
v^{0}\rVert_{W^{\mu,1}(\mathbb{R}^{d})},\quad\forall\;0\leqslant t\leqslant
T,$

where $\mu=\max(r,\kappa)$.

The proof of Theorem 6 strongly relies on Proposition 5 and the following
result:

######  Proposition 7.

Under hypotheses (57) through (62), there exists a constant $C_{T}>0$,
depending only on $T$ and on the parameters of problems (19) and (54), such
that their respective solutions satisfy, for all $t\in[0,T]$,

$\sum\limits_{i\in\mathcal{J}_{h}}|v(t,x_{i}(t))-\nu_{i}(t)|w_{i}(t)\leqslant
C_{T}h^{k-1}\lVert v^{0}\rVert_{W^{\mu,1}(\mathbb{R}^{d})}.$ (67)

###### Proof.

Consider $\beta_{\varepsilon}$ as in (88). We define
$e=\\{e_{i}(\cdot)\\}_{i\in\mathcal{J}}$ where for all $i\in\mathcal{J}$ and
$t\in[0,T]$,

$\displaystyle e_{i}(t)$ $\displaystyle:=v(t,x_{i}(t))-\nu_{i}(t),$
$\displaystyle e_{\varepsilon,i}(t)$
$\displaystyle:=\beta_{\varepsilon}(e_{i}(t))w_{i}(t),$

and compute

$\dot{e}_{\varepsilon,i}(t)=\beta^{\prime}_{\varepsilon}(e_{i}(t))\dot{e}_{i}(t)w_{i}(t)+\beta_{\varepsilon}(e_{i}(t))\dot{w}_{i}(t).$
(68)

We recall that

$\displaystyle\dot{e}_{i}(t)=$
$\displaystyle\biggl{(}G_{\nu,w}(t,x_{i}(t))-\mathcal{G}_{v}(t,x_{i}(t))\biggr{)}\nu_{i}(t)+\biggl{(}A_{\nu,w}(t,x_{i}(t))-\mathcal{A}_{v}(t,x_{i}(t))\biggr{)}\nabla
v(t,x_{i}(t))$
$\displaystyle-\mathcal{G}_{v}(t,x_{i}(t))e_{i}(t)+\Delta\mathcal{M}(t,x_{i}(t)),$

where

$\displaystyle G_{\nu,w}(t,x_{i}(t))$ $\displaystyle:=\text{div
}A_{\nu,w}(t,x_{i}(t))-R(t,x_{i}(t),I_{g}(t,x_{i}(t),\nu,w)),$
$\displaystyle\mathcal{G}_{v}(t,x_{i}(t))$ $\displaystyle:=\text{div
}\mathcal{A}_{v}(t,x_{i}(t))+R(t,x_{i}(t),(I_{g}v)(t,x_{i}(t))),$
$\displaystyle\Delta\mathcal{M}(t,x_{i}(t))$
$\displaystyle:=\int\limits_{\mathbb{R}^{d}}{m(t,x_{i}(t),y,(I_{d}v)(t,x_{i}))v(t,y)dy}-\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)\nu_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t),\nu,w)).$

The functions $a$ and $R$ being Lipschitz, there exists a constant $C$
depending on the parameters of the problem and the value $\overline{\rho}$,
such that

$\displaystyle\biggl{|}G_{\nu,w}(t,x_{i}(t))-\mathcal{G}_{v}(t,x_{i}(t))\biggr{|}\leqslant$
$\displaystyle C\biggl{(}|(I_{a}v)(t,x_{i}(t))-I_{a}(t,x_{i}(t),\nu,w)|$
$\displaystyle+|\partial_{x}(I_{a}v)(t,x_{i}(t))-\partial_{x}I_{a}(t,x_{i}(t),\nu,w)|$
$\displaystyle+|(I_{g}v)(t,x_{i}(t))-I_{g}(t,x_{i}(t),\nu,w)|\biggr{)}.$

Notice that

$\displaystyle|(I_{a}v)(t,x_{i}(t))-I_{a}(t,x_{i}(t),\nu,w)|=$
$\displaystyle\biggl{|}\int\limits_{\mathbb{R}^{d}}\psi_{a}(t,x_{i}(t),y)v(y)dy-\sum\limits_{j\in\mathcal{J}_{h}}\psi_{a}(t,x_{i}(t),x_{j}(t))\nu_{j}(t)w_{j}(t)\biggr{|}$
$\displaystyle\leqslant$
$\displaystyle\biggl{|}\int\limits_{\mathbb{R}^{d}}\psi_{a}(t,x_{i}(t),y)v(y)dy-\sum\limits_{j\in\mathcal{J}_{h}}\psi_{a}(t,x_{i}(t),x_{j}(t))v(t,x_{j}(t))w_{j}(t)\biggr{|}$
$\displaystyle+\biggl{|}\sum\limits_{j\in\mathcal{J}_{h}}\psi_{a}(t,x_{i}(t),x_{j}(t))\biggl{(}v(t,x_{j}(t)-\nu_{j}(t))\biggr{)}w_{j}(t)\biggr{|}$
$\displaystyle\leqslant$ $\displaystyle
C\biggl{(}h^{\kappa}\|v\|_{W^{\kappa,1}(\mathbb{R}^{d})}+\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)},$

where in the last line we have used Lemma 8 and the
$W^{\kappa,1}(\mathbb{R}^{d})$ regularity of $v$. Similar results are true for
$|\partial_{x}(I_{a}v)(t,x_{i}(t))-\partial_{x}I_{a}(t,x_{i}(t),\nu,w)|$ and
$|(I_{g}v)(t,x_{i}(t))-I_{g}(t,x_{i}(t),\nu,w)|$. In conclusion, thanks to
(27), there exists a constant $C_{T}$, only depending on $T$, the parameters
of the problem and the value $\overline{\rho}$, such that

$\displaystyle\biggl{|}G_{\nu,w}(t,x_{i}(t))-\mathcal{G}_{v}(t,x_{i}(t))\biggr{|}\leqslant
C_{T}\biggl{(}h^{\kappa}\|v^{0}\|_{W^{\kappa,1}(\mathbb{R}^{d})}+\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)}.$
(69)

Again, using the Lipschitz regularity of $a$, we see that

$\displaystyle\biggl{|}A_{\nu,w}(t,x_{i}(t))-\mathcal{A}_{v}(t,x_{i}(t))\biggr{|}\leqslant$
$\displaystyle C|(I_{a}v)(t,x_{i}(t))-I_{a}(t,x_{i}(t),\nu,w)|$
$\displaystyle\leqslant$ $\displaystyle
C_{T}\biggl{(}h^{\kappa}\|v^{0}\|_{W^{\kappa,1}(\mathbb{R}^{d})}+\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)}.$
(70)

The boundedness of $a$ and $R$ implies the existence of a constant
$\overline{G}$ such that for all $i\in\mathcal{J}$ and $t\in[0,T]$,

$|\mathcal{G}_{v}(t,x_{i}(t))|\leqslant\overline{G}.$ (71)

We recall from the previous section that

$\Delta\mathcal{M}(t,x_{i}(t))=0,\ \forall i\not\in\mathcal{J}^{m}_{h},$

where the set of indexes $\mathcal{J}^{m}_{h}$ has a finite number of
elements, which depends on $T$. On the other hand, for those
$i\in\mathcal{J}^{m}_{h}$, we have

$\displaystyle|\Delta\mathcal{M}(t,x_{i}(t))|$ $\displaystyle\leqslant$
$\displaystyle\biggl{|}\int\limits_{\mathbb{R}^{d}}{m(t,x_{i}(t),y,(I_{d}v)(t,x_{i}(t)))v(t,y)dy}-\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)v(t,x_{j}(t))m(t,x_{i}(t),x_{j}(t),(I_{d}v)(t,x_{i}(t)))\biggr{|}$
$\displaystyle+\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)v(t,x_{j}(t))\biggl{|}m(t,x_{i}(t),x_{j}(t),(I_{d}v)(t,x_{i}(t)))-m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t),\nu,w))\biggr{|}$
$\displaystyle+\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)|e_{j}(t)|m(t,x_{i}(t),x_{j}(t),(I_{d}v)(t,x_{i}(t)))$
$\displaystyle\leqslant$ $\displaystyle
C\biggl{(}h^{\kappa}\|v\|_{W^{\kappa,1}(\mathbb{R}^{d})}+\mu|(I_{d}v)(t,x_{i}(t))-I_{d}(t,x_{i}(t),\nu,w)|\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)v(t,x_{j}(t))+\overline{M}\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)},$

where we have used again Lemma 8, the $W^{\kappa,1}(\mathbb{R}^{d})$
regularity of $m(t,x,y,I)v(t,y)$ with respect to the $y$ variable, and the
Lipschitz regularity of $m$. Furthermore, Lemma 8 gives us the bound

$\sum\limits_{j\in\mathcal{J}_{h}}w_{j}(t)v(t,x_{j}(t))\leqslant\|v\|_{L^{1}(\mathbb{R}^{d})}+Ch^{\kappa}\|v\|_{W^{\kappa,1}(\mathbb{R}^{d})},$

which together with manipulations similar to those made for
$|(I_{a}v)(t,x_{i}(t))-I_{a}(t,x_{i}(t),\nu,w)|$, and the bound (27), gives

$|\Delta\mathcal{M}(t,x_{i}(t))|\leqslant
C_{T}\biggl{(}h^{\kappa}\|v^{0}\|_{W^{\kappa,1}(\mathbb{R}^{d})}+\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)}.$
(72)

We denote as $\mathcal{K}$ an arbitrary finite subset of $\mathcal{J}_{h}$. If
we add (68) for all values of $i\in\mathcal{K}$, and use bounds (69) through
(72), together with the equation for $w_{i}$, we get

$\displaystyle\left(\sum\limits_{i\in\mathcal{K}}e_{\varepsilon,i}(t)\right)^{\prime}\leqslant$
$\displaystyle
C_{T}B(t)\biggl{(}h^{\kappa}\|v^{0}\|_{W^{\kappa,1}(\mathbb{R}^{d})}+\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)},$

where

$B(t):=\sum\limits_{i\in\mathcal{K}}(\nu_{i}(t)+|\nabla
v(t,x_{i}(t))|)w_{i}(t)+\overline{G}+|\mathcal{K}\cap\mathcal{J}^{m}_{h}|h^{d}e^{\tilde{a}T}+\tilde{a}.$

Given that

$\sum\limits_{i\in\mathcal{K}}\nu_{i}(t)w_{i}(t)\leqslant\overline{\rho},$

and

$\sum\limits_{i\in\mathcal{K}}w_{j}(t)|\nabla v(t,x_{j}(t))|\leqslant\|\nabla
v\|_{L^{1}(\mathbb{R}^{d})}+Ch^{\kappa-1}\|\nabla
v\|_{W^{\kappa-1,1}(\mathbb{R}^{d})},$

thanks to (27), we conclude that there exists a constant $B_{T}$, independent
of the choice of $\mathcal{K}$ and $h$, such that $B(t)\leqslant B_{T}$.
Consequently, for all values of $t\in[0,T]$, $h$ small enough and any finite
subset of $\mathcal{J}_{h}$, we have the relation

$\left(\sum\limits_{i\in\mathcal{K}}e_{\varepsilon,i}(t)\right)^{\prime}\leqslant
C_{T}\biggl{(}h^{\kappa}+\sum\limits_{j\in\mathcal{J}_{h}}|e_{j}(t)|w_{j}(t)\biggr{)}.$

Integrating between $0$ and $t$, taking the limit when $\varepsilon$ goes to
zero, and using Grönwall’s lemma, we obtain that there exists a constant
$C_{T}$, independent of $\mathcal{K}$ such that

$\sum\limits_{i\in\mathcal{K}}|e_{i}(t)|w_{i}(t)\leqslant
C_{T}h^{\kappa}\|v^{0}\|_{W^{\kappa,1}(\mathbb{R}^{d})}.$

Being $C_{T}$ independent of $\mathcal{K}$ and $h$, (67) is immediate. ∎

In other words, we proved in Proposition 7 that the piece-wise constant
functions that take values $v(t,x_{k}(t))$ and $v_{k}(t)$ respectively over
the intervals $\Omega_{k}(t)$ are close in $L^{1}(\mathbb{R}^{d})$.

###### Proof of Theorem 6..

According to the triangle inequality,

$\displaystyle\lVert
v-v_{\varepsilon}^{h}\rVert_{L^{1}(\mathbb{R}^{d})}\leqslant\lVert
v-{\Pi}_{\varepsilon}^{h}(t)v\rVert_{L^{1}(\mathbb{R}^{d})}+\lVert{\Pi}_{\varepsilon}^{h}(t)v-v_{\varepsilon}^{h}\rVert_{L^{1}(\mathbb{R}^{d})},$
(73)

it only remains to bound both terms on the right hand side.

  * _i_)

According to Proposition 5 with $p=1$, and bound (27)

$\displaystyle\lVert
v-{\Pi}_{\varepsilon}^{h}(t)v\rVert_{L^{1}(\mathbb{R}^{d})}\leqslant
C(\varepsilon^{r}+\big{(}\frac{\varepsilon}{h}\big{)}^{\kappa})\lVert
v\rVert_{\mu,p}\leqslant
C_{T}(\varepsilon^{r}+\big{(}\frac{\varepsilon}{h}\big{)}^{\kappa})\lVert
v^{0}\rVert_{\mu,p}.$

  * _i_ i)

On the other hand, one computes

$\displaystyle\lVert{\Pi}_{\varepsilon}^{h}(t)v-v_{\varepsilon}^{h}\rVert_{L^{1}(\mathbb{R}^{d})}$
$\displaystyle=\int_{\mathbb{R}^{d}}{\big{\lvert}\sum\limits_{i\in\mathcal{J}_{h}}{w_{i}(t)\varphi_{\varepsilon}(x-x_{i}(t))\big{(}v(t,x_{i}(t))-\nu_{i}(t)\big{)}}\big{\rvert}dx}$
$\displaystyle\leqslant\sum\limits_{i\in\mathcal{J}_{h}}{\biggl{(}w_{i}(t)\lvert
v(t,x_{i}(t))-\nu_{i}(t)\rvert\int_{\mathbb{R}^{d}}{\lvert\varphi_{\varepsilon}(x-x_{i}(t))\rvert
dx}\biggr{)}}.$

According to the definition of $\varphi_{\varepsilon}$, with the change of
variable $x^{\prime}=\frac{x-x_{k}(t)}{\varepsilon}$, we note that

$\int_{\mathbb{R}^{d}}{\lvert\varphi_{\varepsilon}(x-x_{k}(t))\rvert
dx}=\int_{\mathbb{R}^{d}}{\lvert\varphi(x)\rvert dx}<+\infty,$

by hypothesis on $\varphi$. We have then, according to Proposition 67, that

$\lVert{\Pi}_{\varepsilon}^{h}(t)v-v_{\varepsilon}^{h}\rVert_{L^{1}(\mathbb{R}^{d})}\leqslant
C_{T}h^{\kappa}\lVert v^{0}\rVert_{\kappa,1}\leqslant Ch^{\kappa}\lVert
v^{0}\rVert_{\mu,1},$

which concludes the proof of Theorem 6.

∎

###  4.2 Asymptotic preserving properties

The study of the asymptotic behaviour of the solution for adaptive dynamics
models, such as (19), is one of the main interests often treated in the
literature (see [5, 14, 13, 24, 6, 17, 16, 15, 11, 12]). For this reason, the
design of numerical methods which preserve the asymptotic behaviour, or at
least, the identification of the problems for which the asymptotics are
preserved under a certain numerical scheme, is a priority. In other words,
given that $v(t,\cdot)$ converges to a measure $\mu$ when $t$ goes to
infinity, we expect to identify the conditions under which $\lim\limits_{h\to
0}\lim\limits_{t\to+\infty}v^{h}_{\varepsilon(h)}(t,\cdot)=\mu$, that is,
ensuring the commutativity of diagram (74):

$\begin{matrix}v(t,\cdot)&\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern
3.0pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern
0.0pt\offinterlineskip\halign{\entry@#!@&&\entry<EMAIL_ADDRESS>0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{{}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces{\hbox{\kern
12.86252pt\raise 0.0pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{$\scriptstyle{}$}}}\kern
3.0pt}}}}}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern
7.00002pt\raise-5.25049pt\hbox{{}\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern
0.0pt\raise-2.2505pt\hbox{$\scriptstyle{t\rightarrow\infty}$}}}\kern
3.0pt}}}}}}\ignorespaces{\hbox{\kern 28.72504pt\raise 0.0pt\hbox{\hbox{\kern
0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern
28.72504pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}\ignorespaces}}}}\ignorespaces&\mu\\\ &&\\\
{\scriptstyle h\rightarrow 0}{\left\uparrow\vbox
to14.22636pt{}\right.\kern-1.2pt}&&\hskip 17.07182pt{\left\uparrow\vbox
to14.22636pt{}\right.\kern-1.2pt}{\scriptstyle h\rightarrow 0}\\\ &&\\\
v^{h}_{\varepsilon(h)(t)}&\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern
3.0pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern
0.0pt\offinterlineskip\halign{\entry@#!@&&\entry<EMAIL_ADDRESS>0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise
0.0pt\hbox{$\textstyle{{}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces{\hbox{\kern
12.86252pt\raise 0.0pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{$\scriptstyle{}$}}}\kern
3.0pt}}}}}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern
7.00002pt\raise-5.25049pt\hbox{{}\hbox{\kern 0.0pt\raise
0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern
0.0pt\raise-2.2505pt\hbox{$\scriptstyle{t\rightarrow\infty}$}}}\kern
3.0pt}}}}}}\ignorespaces{\hbox{\kern 28.72504pt\raise 0.0pt\hbox{\hbox{\kern
0.0pt\raise
0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern
28.72504pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern
3.0pt\raise
0.0pt\hbox{$\textstyle{}$}}}}}}}\ignorespaces}}}}\ignorespaces&\mu_{h}\end{matrix}$
(74)

In what follows, we formally define the concept of an asymptotic preserving
approximation, and give examples and counter-examples of this concept.

We recall that, according to the Riesz representation theorem, the space of
finite Radon measures can be identified with the topological dual space of
$\mathcal{C}_{c}(\mathbb{R}^{d})$. Hence, we say that a sequence of finite
Radon measures $\\{\mu_{n}\\}_{n\in\mathbb{N}}$ converges weakly to a finite
Radon measure $\mu$ (denoted
$\mu_{n}\underset{n\to+\infty}{\rightharpoonup}\mu$) if for all
$\phi\in\mathcal{C}_{c}(\mathbb{R}^{d})$,

$\int_{\mathbb{R}^{d}}{\phi(x)d\mu_{n}(x)}\underset{n\to+\infty}{\longrightarrow}\int_{\mathbb{R}^{d}}{\phi(x)d\mu(x)}.$

This leads us to introduce the following definition:

######  Definition 1.

We say that the particle solution $v^{h}_{\varepsilon}$ defined in (66) is an
_asymptotic preserving approximation_ of $v$, the solution to (19), if for all
$\varepsilon:(0,1]\to\mathbb{R}_{+}^{*}$ which converges to $0$ when $h$ goes
to $0$, and all $\phi\in\mathcal{C}_{c}(\mathbb{R}^{d})$,

$\underset{t\to+\infty}{\limsup}\left\lvert\int_{\mathbb{R}^{d}}{\phi(x)v_{\varepsilon(h)}^{h}(t,x)dx}-\int_{\mathbb{R}^{d}}{\phi(x)v(t,x)dx}\right\rvert\underset{h\to
0}{\longrightarrow}0.$

The following lemma ensures that, in the previous definition,
$v_{\varepsilon(h)}^{h}$, introduced in (66), can be replaced by $v^{h}$,
introduced in (65).

######  Lemma 9.

The function $v_{\varepsilon}^{h}$ is an asymptotic preserving approximation
of $v$ if and only if

$\underset{t\to+\infty}{\limsup}\left\lvert\int_{\mathbb{R}^{d}}{\phi(x)dv^{h}(t,x)}-\int_{\mathbb{R}^{d}}{\phi(x)v(t,x)dx}\right\rvert\underset{h\to
0}{\longrightarrow}0.$

###### Proof.

Let us prove that for all $\varepsilon:(0,1]\to\mathbb{R}_{+}^{*}$ which
converges to $0$ as $h$ goes to $0$, and all
$\phi\in\mathcal{C}_{c}(\mathbb{R}^{d})$,

$\underset{t\to+\infty}{\limsup}\left\lvert\int_{\mathbb{R}^{d}}{\phi(x)v^{h}_{\varepsilon(h)}(t,x)dx}-\int_{\mathbb{R}^{d}}{\phi(x)dv^{h}(t,x)}\right\rvert\underset{h\to
0}{\longrightarrow}0.$

Let $h>0$. According to the definitions of $v^{h}$ and $v^{h}_{\varepsilon}$,
and since $\int_{\mathbb{R}^{d}}{\varphi_{\varepsilon(h)}(x-x_{i}(t))dx}=1$
for all $i\in\mathcal{J}_{h}$ we get, for all $t\geqslant 0$,

$\displaystyle\left\lvert\int_{\mathbb{R}^{d}}{\phi(x)v^{h}_{\varepsilon(h)}(t,x)dx}-\int_{\mathbb{R}^{d}}{\phi(x)dv^{h}(t,x)}\right\rvert=$
$\displaystyle\left\lvert\int_{\mathbb{R}^{d}}{\sum\limits_{i\in\mathcal{J}_{h}}{\nu_{i}(t)w_{i}(t)\big{(}\phi(x)-\phi(x_{i}(t))\big{)}\varphi_{\varepsilon(h)}(x-x_{i}(t))}dx}\right\rvert$
$\displaystyle\leqslant\sum\limits_{i\in\mathcal{J}_{h}}{\lvert\nu_{i}(t)w_{i}(t)\rvert\int_{\mathbb{R}^{d}}{\big{\lvert}\big{(}\phi(x)-\phi(x_{i}(t))\big{)}\varphi_{\varepsilon(h)}(x-x_{i}(t))\big{\rvert}dx}}.$

With the change of variable ‘$y=\frac{x-x_{i}(t)}{\varepsilon(h)}$’, we get,
for all $i\in\mathcal{J}_{h}$,

$\int_{\mathbb{R}^{d}}{\big{\lvert}\big{(}\phi(x)-\phi(x_{i}(t))\big{)}\varphi_{\varepsilon(h)}(x-x_{i}(t))\big{\rvert}dx}=\int_{K}{\big{\lvert}\big{(}\phi(\varepsilon(h)y+x_{i}(t))-\phi(x_{i}(t))\big{)}\varphi(y)\big{\rvert}dy},$

where $K$ is the support of $\varphi$. Let $\eta>0$. Since $\phi$ is
continuous with a compact support, and thus uniformly continuous, then
$\lvert\phi(\varepsilon(h)x+x_{i}(t))-\phi(x_{i}(t))|\leqslant\eta$ for all
$i\in\mathcal{J}_{h}$, $x\in K$, $t\geqslant 0$ and any $h$ small enough.
Therefore, for any $h$ small enough,

$\left\lvert\int_{\mathbb{R}^{d}}{\phi(x)v^{h}_{\varepsilon(h)}(t,x)dx}-\int_{\mathbb{R}^{d}}{\phi(x)dv^{h}(t,x)}\right\rvert\leqslant\eta\sum\limits_{i\in\mathcal{J}_{h}}{\nu_{i}(t)w_{i}(t)},$

which concludes the proof, since there exists $\overline{\rho}>0$ such that
$0\leqslant\sum\limits_{i\in\mathcal{J}_{h}}{\nu_{i}(t)w_{i}(t)}\leqslant\overline{\rho}$
for all $h>0$, $t\geqslant 0$, as proved in Theorem 4. ∎

The problem of determining if $v_{\varepsilon}^{h}$ is an asymptotic
preserving approximation of $v$ is generally a difficult question. In what
follows, we deal with cases where we are able to determine the asymptotic
behaviour of both $v$ and $v^{h}$, and we check if the necessary and
sufficient condition from Lemma 9 holds. From now on, we assume that $a$ is
local and not time dependent, i.e. $a(t,x,I)=a(x)$ and that that the functions
$m$, $R$, $\psi_{g}$ and $\psi_{d}$ are not time-dependent. We assume as well
that the function $(x,y,I)\mapsto m(x,y,I)$ is not only uniformly compactly
supported as a function of the $x$ variable, but relative to the $y$ variable
as well. That is, there exist two compact sets $K_{x}$ and $K_{y}$ such that
$\sup\limits_{y,I}m(x,y,I)=0$ for all $y$ outside of $K_{x}$ and
$\sup\limits_{x,I}m(x,y,I)=0$ for all $x$ outside of $K_{y}$. We denote
$K_{xy}:=K_{x}\cup K_{y}$. Finally, we assume $\psi_{g}$ and $\psi_{d}$ to be
compactly supported as functions of the $y$ variable.

####  4.2.1 Necessary conditions of convergence towards a Radon measure

With the help of necessary conditions, we would be able to rule out those
cases where $v$ does not converge towards certain types of Radon measures,
which are the object of our interest. We start by giving a general result,
involving the necessary conditions of convergence towards any Radon measure.

######  Lemma 10.

Let us assume that $v(t,\cdot)\underset{t\to+\infty}{\rightharpoonup}\mu$ in
the weak sense in the space of finite Radon measures. Then, for
$\gamma\in\\{g,d\\}$,

$I_{\alpha}(t,x)\underset{t\to+\infty}{\longrightarrow}\int_{\mathbb{R}^{d}}{\psi_{\alpha}(x,y)d\mu(y)}=:\overline{I_{\alpha}}(x)\quad\forall
x\in\mathbb{R}^{d},$

and for all $\phi\in\mathcal{C}_{0}^{1}(\mathbb{R}^{d})$,

$\int_{\mathbb{R}^{d}}{\bigg{(}a(x)\cdot\nabla\phi(x)+R\left(x,\overline{I_{g}}(x)\right)\phi(x)\bigg{)}d\mu(x)}+\int_{\mathbb{R}^{d}}{\left(\int_{\mathbb{R}^{d}}{m\left(x,y,\overline{I_{d}}(x)\right)}\phi(x)dx\right)}d\mu(y)=0.$

###### Proof.

Let us assume that $v(t,\cdot)\underset{t\to+\infty}{\rightharpoonup}\mu$. By
definition of the weak convergence in the space of finite Radon measure, for
any $\phi\in\mathcal{C}_{c}^{1}(\mathbb{R}^{d})$,

$\int_{\mathbb{R}^{d}}{\phi(x)v(t,x)dx}\underset{t\to+\infty}{\rightarrow}\int_{\mathbb{R}^{d}}{\phi(x)d\mu(x)}.$

The first identity is thus a direct consequence of the definition of the weak
convergence, applied with $\phi=\psi_{\gamma}(x,\cdot)$, $\gamma\in\\{g,d\\}$.

Let $\phi\in\mathcal{C}_{c}(\mathbb{R}^{d})$. One computes

$\displaystyle\frac{d}{dt}\int_{\mathbb{R}^{d}}{\phi(x)v(t,x)dx}=$
$\displaystyle-\int_{\mathbb{R}^{d}}{\phi(x)\nabla\cdot\left(a(x)v(t,x)\right)}+\int_{\mathbb{R}^{d}}{R(x,I_{g}(t,x))\phi(x)v(t,x)dx}$
$\displaystyle+\int_{\mathbb{R}^{d}}{\left(\int_{\mathbb{R}^{d}}{m(x,y,I_{d}(t,x))v(t,y)}\right)\phi(x)dx}$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{d}}{\bigg{(}a(x)\cdot\nabla\phi(x)+R(x,I_{g}(t,x))\phi(x)\bigg{)}v(t,x)dx}$
$\displaystyle+\int_{\mathbb{R}^{d}}{\bigg{(}\int_{\mathbb{R}^{d}}{m(x,y,I_{d}(t,x))\phi(x)dx}\bigg{)}v(t,y)dy}$
$\displaystyle\underset{t\to+\infty}{\longrightarrow}$
$\displaystyle\int_{\mathbb{R}^{d}}{\bigg{(}a(x)\cdot\nabla\phi(x)+R\left(x,\overline{I_{g}}(x)\right)\phi(x)\bigg{)}d\mu(x)}$
$\displaystyle+\int_{\mathbb{R}^{d}}{\left(\int_{\mathbb{R}^{d}}{m\left(x,y,\overline{I_{d}}(x)\right)\phi(x)dx}\right)d\mu(y)}.$

Thus, $t\mapsto\int_{\mathbb{R}^{d}}{\phi(x)v(t,x)dx}$ is a convergent
function with a convergent derivative, which ensures that the limit of its
derivative is zero, which concludes the proof. ∎

The following proposition provides a necessary condition for the convergence
to a sum of Dirac masses.

######  Proposition 8.

Let us assume that
$v(t,\cdot)\underset{t\to+\infty}{\rightharpoonup}\sum\limits_{i=1}^{N}{C_{i}\delta_{x_{i}}}$,
with $x_{1},\ldots,x_{N}\in\mathbb{R}^{d}$, $C_{1},\ldots,C_{N}>0$. Then, for
$\alpha=g,d$,
$I_{\alpha}(t,x)\underset{t\to+\infty}{\rightarrow}\sum\limits_{i=1}^{N}{C_{i}\psi_{\alpha}(x,x_{i})}=:\overline{I_{\alpha}}(x)$.
Moreover, for all $i\in\\{1,\ldots,N\\}$, $a(x_{i})=0$,
$R\left(x_{i},\overline{I_{g}}(x_{i})\right)=0$, and for all
$x\in\mathbb{R}^{d},$ $m(x,x_{i},\overline{I_{d}}(x))=0$.

###### Proof.

According to the previous lemma, for all
$\phi\in\mathcal{C}_{c}^{1}(\mathbb{R}^{d})$,

$\sum\limits_{i=1}^{N}{C_{i}a(x_{i})\cdot\nabla\phi(x_{i})}+\sum\limits_{i=1}^{N}{C_{i}R\left(x_{i},\overline{I_{g}}(x_{i})\right)\phi(x_{i})}+\sum\limits_{i=1}^{N}{C_{i}\int_{\mathbb{R}^{d}}{m\left(x,x_{i},\overline{I_{d}}(x)\right)\phi(x)dx}}=0.$

Let $\varepsilon>0$. For any non-negative function
$\phi_{\varepsilon}\in\mathcal{C}_{c}^{1}(\mathbb{R}^{d})$ with a support on
$\mathbb{R}^{d}\backslash\bigcup_{i=1}^{N}{B(x_{i},\varepsilon)}$, we have
that

$\int_{\mathbb{R}^{d}}{\sum\limits_{i=1}^{N}{C_{i}m\left(x,x_{i},\overline{I_{d}}(x)\right)\phi_{\varepsilon}(x)}dx}=0,$

which proves that
$x\mapsto\sum\limits_{i=1}^{N}{C_{i}m\left(x,x_{i},\overline{I_{d}}(x)\right)}$
is $0$ on $\mathbb{R}^{d}\backslash\bigcup_{i=1}^{N}{B(x_{i},\varepsilon)}$.
Since $m$ is a non-negative function, for all $i\in\\{1,...N\\}$, $x\mapsto
m\left(x,x_{i},\overline{I_{d}}(x)\right)$ is $0$ as well on
$\mathbb{R}^{d}\backslash\bigcup_{i=1}^{N}{B(x_{i},\varepsilon)}$, and
therefore on $\mathbb{R}^{d}$, since the result holds for any $\varepsilon>0$.
Hence,

$\sum\limits_{i=1}^{N}{C_{i}a(x_{i})\cdot\nabla\phi(x_{i})}+\sum\limits_{i=1}^{N}{C_{i}R\left(x_{i},\overline{I_{g}}(x_{i})\right)\phi(x_{i})}=0,$

for any $\phi\in\mathcal{C}_{c}^{1}(\mathbb{R}^{d})$. By choosing
$\phi^{1}_{j}$ such that

$\phi^{1}_{j}(x_{i})=\delta_{ij}\mbox{ and }\nabla\phi^{1}_{j}(x_{i})=0,\
i,j\in\\{1,...N\\},$

where $\delta_{ij}$ represents the Kronecker delta, one proves that
$R\left(x_{i},\overline{I_{g}}(x_{i})\right)=0$ for all $i\in\\{1,...N\\}$.
Finally, by choosing $\phi^{2}_{ij}$ such that

$\nabla\phi^{2}_{jl}(x_{i})=\delta_{ij}e_{l},\ i,j\in\\{1,...N\\},\
l\in\\{1,...d\\},$

where $\\{e_{l}\\}_{l=1}^{d}$ represents the euclidean basis in
$\mathbb{R}^{d}$, one proves that $a(x_{i})=0$, for any
$i\in\\{1,\ldots,N\\}$. ∎

####  4.2.2 Limit identification and asymptotic preserving approximations

In some cases, it is possible to guarantee the existence of a limit for
$v^{h}$ and identify it.
Assume that there exists $\hat{x}\in\mathbb{R}^{d}$ an asymptotically stable
equilibrium for the ODE ‘$\dot{x}=a(x)$’ and that there exists $C,\delta>0$
such that

$\forall y\in\mbox{supp}(n^{0}),t\geq 0,\quad\lVert X(t,y)-\hat{x}\rVert\leq
Ce^{-\delta t}.$ (75)

Moreover, let us assume that there exist positive values $D$, $I^{m}$ and
$I^{M}$ such that

$R(x,I^{m})\geqslant 0,\quad R(x,I^{M})\leqslant 0\mbox{ and
}\partial_{I}R(x,I)\leqslant-D,\quad\forall x\in\textrm{supp}(v^{0}).$ (76)

Then, we can compute the limit of $v^{h}$ when $t$ goes to $+\infty$, whatever
the value of $m$, as stated in the following proposition.

######  Proposition 9.

Let us assume that $\textrm{supp}\left(v^{0}\right)$ is a compact set such
that hypotheses (75) and (76) hold. Then, $v^{h}$ converges to
$\hat{\rho}_{h}\>\delta_{\hat{x}}$ in the weak sense in the space of Radon
measures, where $\hat{\rho}_{h}$ is the unique solution of

$R(\hat{x},\psi_{g}(\hat{x},\hat{x})\hat{\rho}_{h})=0.$

The following lemma, proved in Appendix D, is required in the proof of this
result.

######  Lemma 11.

Let $u\in\mathcal{C}^{2}\left(\mathbb{R}_{+},\mathbb{R}\right)$ be a bounded
function, and let us assume that there exist $p_{0}>0$,
$p:\mathbb{R}_{+}\to\mathbb{R}_{+}$ a function which satisfies $p\geqslant
p_{0}$ and $B\in L^{1}(\mathbb{R}_{+})$ an integrable function such that

$\ddot{u}(t)\geqslant-p(t)\dot{u}(t)+B(t).$

Then, there exists $u_{\infty}\in\mathbb{R}$ such that
$\lim\limits_{t\rightarrow+\infty}u(t)=u_{\infty}$.

###### Proof of Proposition 9.

Given that $\textrm{supp}\left(v^{0}\right)\cup K_{xy}$ is a compact set which
is strictly contained in the basin of attraction of $\hat{x}$ we will have the
existence of $\mathcal{J}^{0}_{h}\subset\mathcal{J}_{h}$, with
$\lvert\mathcal{J}^{0}_{h}\rvert<+\infty$ such that $\nu^{h}(t)\not\equiv 0$
only for $i\in\mathcal{J}^{0}_{h}$.
Let us denote, for all $i\in\mathcal{J}^{0}_{h}$,
$\alpha_{i}(t):=\nu_{i}(t)w_{i}(t)$, and

$\rho_{h}(t):=\sum\limits_{i\in\mathcal{J}_{h}}{\alpha_{i}(t)}=\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\alpha_{i}(t)}.$

Let us note that, according to the hypotheses on $a$, for all
$i\in\mathcal{J}^{0}_{h}$, $x_{i}(t)$ converges to $\hat{x}$. Thus,

$v^{h}(t)-\rho_{h}(t)\delta_{\hat{x}}\underset{t\to+\infty}{\rightharpoonup}0.$

Hence, it only remains to prove that $\rho_{h}$ converges to the expected
limit. According to the definition of $\rho_{h}$,

$\displaystyle\dot{\rho}_{h}(t)=\sum\limits_{i\in\mathcal{J}^{0}_{h}}{R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}+\underbrace{\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t)m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)}}_{:=\varepsilon(t)}.$
(77)

According to hypothesis 75, there exist $C,\delta>0$ such that

$\left\lVert x_{i}(t)-\hat{x}\right\rVert\leqslant Ce^{-\delta t},\quad\forall
i\in\mathcal{J}^{0}_{h},\quad\forall t\geqslant 0.$

Thus, for all $i\in\mathcal{J}^{0}_{h}$,

$w_{i}(t)=e^{\int_{0}^{t}{\text{div
}a(x_{i}(s))ds}}w_{i}(0)\leqslant\underbrace{e^{\int_{0}^{t}{\text{div
}a(x_{i}(s))-\text{div }a(\hat{x})ds}}}_{\leqslant\tilde{C}}e^{\text{div
}a(\hat{x})t}w_{i}(0),$

which proves, since $\text{div }a(\hat{x})<0$, that there exist
$C^{\prime},\delta^{\prime}$ such that for all $t\geqslant 0,$

$\displaystyle
0\leqslant\underset{i\in\mathcal{J}^{0}_{h}}{\max}\;w_{i}(t)\leqslant
C^{\prime}e^{-\delta^{\prime}t}.$ (78)

In particular, it proves that $t\mapsto\varepsilon(t)$ defined in (77),
converges to zero with an exponential speed, since

$\displaystyle\lvert\varepsilon(t)\rvert\leqslant\underset{i\in\mathcal{J}_{h}}{\max}\;w_{i}(t)\lvert\mathcal{J}^{0}_{h}\rvert\lVert
m\rVert_{L^{\infty}}\rho_{h}(t),$

and $\rho_{h}$ is bounded, according to Theorem 4.

Moreover, according to the hypothesis on $\psi_{g}$,
$\underline{\psi_{g}}\,\rho(t)\leqslant
I_{g}(t,x_{i}(t))\leqslant\lVert\psi_{g}\rVert_{\infty}\,\rho(t)$. Thus,
according to the hypotheses (76) on $R$, the relation

$\dot{\rho}_{h}(t)\geqslant\left(\underset{i\in\mathcal{J}^{0}_{h}}{\min}\left(R(x_{i}(t),I_{g}(t,x_{i}(t)))\right)\right)\,\rho_{h}(t)$

implies that, as soon as $\rho_{h}$ becomes small, and so does $I_{g}$,
$\rho_{h}$ becomes increasing, which proves that $\rho_{h}$ is lower bounded
by a positive constant. Moreover, since $\rho_{h}$ and $\varepsilon$ are
bounded, and

$\displaystyle\|\dot{\alpha}\|_{\mathcal{l}_{1}}=\sum\limits_{i\in\mathcal{J}^{0}_{h}}|\dot{\alpha}(t)|\leqslant$
$\displaystyle\sum\limits_{i\in\mathcal{J}^{0}_{h}}\left(2|\text{div
}a(x_{i}(t))|+|R(x_{i}(t),I_{g}(t,x_{i}(t)))|\right)\alpha_{i}(t)+\varepsilon(t)$
$\displaystyle\leqslant$
$\displaystyle\left(2\|a\|_{W^{1,\infty}(\mathbb{R}^{d})}+\overline{R}\right)\rho_{h}(t)+\varepsilon(t),$

where

$\overline{R}:=\max\limits_{t\in\mathbb{R},i\in\mathcal{J}^{0}_{h}}|R(x_{i}(t),I_{g}(t,x_{i}(t)))|,$

then $\|\dot{\alpha}\|_{\mathcal{l}_{1}}$ is also bounded.

Now, let us prove that $\rho_{h}$ satisfies the equality of Lemma 11. First,
for $\gamma\in\\{g,d\\}$, we compute

$\displaystyle\left\lvert\frac{d}{dt}I_{\gamma}(t,x_{i}(t))-\psi_{\gamma}(\hat{x},\hat{x})\dot{\rho}_{h}(t)\right\lvert\leqslant$
$\displaystyle\left\lvert\sum\limits_{j\in\mathcal{J}^{0}_{h}}{\bigg{(}a(x_{i}(t))\partial_{x}\psi_{\gamma}(x_{i}(t),x_{j}(t))+a(x_{j}(t))\partial_{y}\psi_{\gamma}(x_{i}(t),x_{j}(t))\bigg{)}\alpha_{j}(t)}\right\rvert$
$\displaystyle+\left\lvert\sum\limits_{j\in\mathcal{J}^{0}_{h}}{\big{(}\psi_{\gamma}(x_{i}(t),x_{j}(t))-\psi_{\gamma}(\hat{x},\hat{x})\big{)}\dot{\alpha}_{j}(t)}\right\rvert$
$\displaystyle\leqslant$ $\displaystyle
2\max\limits_{i\in\mathcal{J}^{0}_{h}}\lvert
a(x_{i}(t))\rvert\lVert\psi_{\gamma}\rVert_{W^{1,\infty}(\mathbb{R}^{d})}\rho_{h}(t)+\max\limits_{i,j\in\mathcal{J}^{0}_{h}}\lvert\psi_{\gamma}(x_{i}(t),x_{j}(t))-\psi_{\gamma}(\hat{x},\hat{x})\rvert\|\dot{\alpha}\|_{\mathcal{l}_{1}}.$

Since hypothesis (75), is satisfied, the functions
$t\mapsto\max\limits_{i\in\mathcal{J}^{0}_{h}}\lvert a(x_{i}(t))\rvert$ and
$t\mapsto\max\limits_{i,j\in\mathcal{J}^{0}_{h}}\lvert\psi_{\gamma}(x_{i}(t),x_{j}(t))-\psi_{\gamma}(\hat{x},\hat{x})\rvert$
converge to zero with an exponential speed. Since
$|\mathcal{J}^{0}_{h}\rvert<+\infty$, this proves that, for
$\gamma\in\\{g,d\\}$,

$\displaystyle\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\left\lvert\frac{d}{dt}I_{\gamma}(t,x_{i}(t))-\psi_{\gamma}(\hat{x},\hat{x})\dot{\rho}_{h}(t)\right\rvert}=O(e^{-\delta
t}),$ (79)

for a certain $\delta>0$.

Thus, by differentiating (77), we get

$\displaystyle\ddot{\rho}_{h}(t)=$
$\displaystyle\sum\limits_{i\in\mathcal{J}^{0}_{h}}{a(x_{i}(t))\partial_{x}R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}+\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\left(\frac{d}{dt}I_{g}(t,x_{i}(t))\right)\partial_{I}R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}$
$\displaystyle+\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\dot{\alpha}_{i}(t)R(x_{i}(t),I_{g}(t,x_{i}(t)))}+\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t){m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\dot{\alpha}_{j}(t)}}$
$\displaystyle+\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{\dot{w}_{i}(t){m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)}}+\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t)\frac{d}{dt}m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)},$

where

  * _i)_

$\sum\limits_{i\in\mathcal{J}^{0}_{h}}{a(x_{i}(t))\partial_{x}R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}=O(e^{-\delta
t})$, since $\max\limits_{i\in\mathcal{J}_{h}}|a(x_{i}(t))|$ converges to zero
with exponential speed, and $\partial_{x}R$ and $\rho_{h}$ are bounded,

  * _ii)_

according to
(79)$,\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\left(\frac{d}{dt}I_{g}(t,x_{i}(t))\right)\partial_{I}R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}=-p(t)\dot{\rho}_{h}(t)+O(e^{-\delta
t})$, where

$p(t)=-\psi_{g}(\hat{x},\hat{x})\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\partial_{I}R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}\geqslant
D\underline{\psi_{g}}\min\limits_{t\geqslant 0}\,\rho_{h}(t)>0,$

  * _iii)_

$\displaystyle\sum\limits_{i\in\mathcal{J}^{0}_{h}}{\dot{\alpha}_{i}(t)R(x_{i}(t),I_{g}(t,x_{i}(t)))}=$
$\displaystyle\underbrace{\sum\limits_{i\in\mathcal{J}^{0}_{h}}{R(x_{i}(t),I_{g}(t,x_{i}(t)))^{2}}\alpha_{i}(t)}_{:=P(t)\geqslant
0}$
$\displaystyle+\underbrace{\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t)R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))}}_{=O(e^{-\delta
t})},$

where the relation for the second term was proved thanks to the bound

$\left\lvert\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t)R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{j}(t)m(t,x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))}\right\rvert\leqslant\overline{M}\,\overline{R}\rho_{h}(t)\lvert\mathcal{J}^{0}_{h}\rvert\max\limits_{i\in\mathcal{J}_{h}}{w_{i}(t)},$

the boundedness of $\rho_{h}$ and inequality (78),

  * _iv)_

$\displaystyle\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t){m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\dot{\alpha}_{j}(t)}}+\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{\dot{w}_{i}(t){m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)}}$
$\displaystyle+\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t)\frac{d}{dt}m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)}$
$\displaystyle=O(e^{-\delta t}),$

since for all $t\geq 0$,

$\displaystyle\left\lvert\sum\limits_{i,j\in\mathcal{J}^{0}_{h}}{w_{i}(t){m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\dot{\alpha}_{j}(t)}}\right\rvert\leqslant\overline{M}\lvert\mathcal{J}^{0}_{h}\rvert\|\dot{\alpha}\|_{\mathcal{l}^{1}}\max\limits_{i\in\mathcal{J}^{0}_{h}}w_{i}(t),$
$\displaystyle\left\lvert\sum\limits_{i,j\in\mathcal{J}_{h}}{\dot{w}_{i}(t){m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)}}\right\rvert\leqslant\overline{M}\rho_{h}(t)\lvert\mathcal{J}^{0}_{h}\rvert\lVert
a\rVert_{W^{1,\infty}(\mathbb{R}^{d})}\max\limits_{i\in\mathcal{J}^{0}_{h}}w_{i}(t),$

and

$\displaystyle\left\lvert\sum\limits_{i,j\in\mathcal{J}_{h}}{w_{i}(t)\frac{d}{dt}m(x_{i}(t),x_{j}(t),I_{d}(t,x_{i}(t)))\alpha_{j}(t)}\right\rvert\leqslant\overline{M}\left(\|a\|_{L^{\infty}(\mathbb{R})}+\big{\lvert}\frac{d}{dt}I_{g}(t,x_{i}(t))\big{\rvert}\right)\rho_{h}(t)\lvert\mathcal{J}^{0}_{h}\rvert\max\limits_{i\in\mathcal{J}^{0}_{h}}w_{i}(t),$

where $\frac{d}{dt}I_{d}(t,x_{i}(t))$ is bounded thanks to (79). For these
three inequalities, we conclude with (78).

Hence, $\ddot{\rho}_{h}(t)\geqslant-p(t)\dot{\rho}_{h}(t)+O(e^{-\delta t})$.
According to Lemma 11, $\rho_{h}$ has a limit when $t$ goes to $\infty$, which
we denote $\hat{\rho}_{h}$. Since
$\dot{\rho}_{h}(t)=\sum\limits_{i\in\mathcal{J}^{0}_{h}}{R(x_{i}(t),I_{g}(t,x_{i}(t)))\alpha_{i}(t)}+\varepsilon(t)$,
which converges to
$R(\hat{x},\psi_{g}(\hat{x},\hat{x})\hat{\rho}_{h})\hat{\rho}_{h}$, we deduce
that $R(\hat{x},\psi_{g}(\hat{x},\hat{x})\hat{\rho_{h}})=0$, since
$\hat{\rho}_{h}>0$. ∎

When $m\equiv 0$ and under the same assumptions for the remaining coefficients
of the problem as in Proposition 9, we are able to identify the limit of $v$
and prove that it coincides with the limit of $v^{h}$. According to Lemma 9,
this ensures that $v^{h}_{\varepsilon}$ is an asymptotic preserving
approximation of $v$.

######  Theorem 10.

Let us assume that there exists $\hat{x}\in\mathbb{R}^{d}$ which is an
asymptotically stable equilibrium for the ODE $\dot{x}=a(x)$ such that
hypothesis (75) holds. We assume as well that $m\equiv 0$ and that hypotheses
(76) hold. Then, $v$ converges to $\hat{\rho}\>\delta_{\hat{x}}$ in the weak
sense in the space of Radon measures, where $\hat{\rho}$ is the unique
solution of

$R(\hat{x},\psi_{g}(\hat{x},\hat{x})\hat{\rho})=0.$

Consequently, $v^{h}_{\varepsilon}$ is an asymptotic preserving approximation
of $v$.

###### Proof.

Let us recall that

$\rho(t)=\int_{\mathbb{R}^{d}}{v(t,x)dx}.$

We recall as well that, the function $a$ being only dependent of $x$, the
characteristic lines $X_{t}(x):=X(t,x)$ satisfy
$X^{-1}_{t}(x):=X^{-1}(t,x)=X_{-t}(x)$ and $X_{t}(X_{s}(x))=X_{t+s}(x)$. By
using the fact that, for any $x\in\mathbb{R}^{d}$,

$v(t,x)=v^{0}(X_{-t}(x))e^{\mathcal{G}(0,t,x)},$

where

$\mathcal{G}(s,t,x):=\int_{s}^{t}{R(X_{\tau-t}(x),(I_{g}v)(\tau,X_{\tau-t}(x)))-\text{div
}a(X_{\tau-t}(x))d\tau},$

and that, by hypothesis, $K:=\text{supp}(v^{0})$ is a compact set included in
the basin of attraction of $\hat{x}$, one proves that
$\text{supp}(v(t,\cdot))$ is the image of $\text{supp}(v^{0})$ by
$X_{t}(\cdot)$. Since $X_{t}(y)$ converges to $\hat{x}$ for all
$y\in\text{supp}(v^{0})$, we prove that
$K_{t}:=\text{supp}(v(t,\cdot))=X_{t}(\text{supp}(v^{0}))=X_{t}(K)$ is a
compact set included in the basin of attraction of $\hat{x}$, for all
$t\geqslant 0$. By (75), there exist $C>0$ and $\delta>0$ such that

$\displaystyle\left\lVert X_{t}(y)-\hat{x}\right\rVert\leqslant Ce^{-\delta
t},\quad\forall y\in K,\quad\forall t\geqslant 0.$

Let $\phi\in\mathcal{C}_{c}(\mathbb{R}^{d})$. By definition of $\rho$ and by
using the change of variable ‘$x=X_{t}(y)$’ we get

$\displaystyle\bigg{\lvert}\int_{\mathbb{R}^{d}}{\phi(x)v(t,x)}-\rho(t)\phi(\hat{x})\bigg{\rvert}$
$\displaystyle\leqslant\int_{K_{t}}{\lvert\phi(x)-\phi(\hat{x})\rvert
v(t,x)dx}$ $\displaystyle=\int_{K}{\lvert\phi(X_{t}(y))-\phi(\hat{x})\rvert
v(t,X_{t}(y))e^{\int_{0}^{t}\text{div }a(X_{t}(y))}dy}$
$\displaystyle\leqslant\max\limits_{y\in
K}\lvert\phi(X_{t}(y))-\phi(\hat{x})\rvert\int_{K}{v(t,X_{t}(y))e^{\int_{0}^{t}\text{div
}a(X_{t}(y))}dy}$ $\displaystyle=\max\limits_{y\in
K}\lvert\phi(X_{t}(y))-\phi(\hat{x})\rvert\rho(t).$

Since $\rho$ is bounded, this proves that

$\displaystyle
v(t,\cdot)-\rho(t)\delta_{\hat{x}}\underset{t\to+\infty}{\rightharpoonup}0.$
(80)

Hence, it only remains to prove that $\rho$ converges to the expected limit.
We see that

$\displaystyle\dot{\rho}(t)=\int_{\mathbb{R}^{d}}{R(x,(I_{g}v)(t,x))v(t,x)}dx.$
(81)

First, let us note that $\rho$ has a positive lower bound. Indeed, according
to the hypothesis on $\psi_{g}$, for all $x\in\mathbb{R}^{d}$,
$\underline{\psi_{g}}\,\rho(t)\leqslant(I_{g}v)(t,x)\leqslant\lVert\psi_{g}\rVert_{L^{\infty}}\,\rho(t)$.
Thus,

$\dot{\rho}(t)\geqslant\left(\underset{x\in
K_{t}}{\min}\left(R(x,(I_{g}v)(t,x)\right)\right)\,\rho(t).$

Hence, as soon as $\rho(t)$ becomes small, and so does $I_{g}(t,x)$, $\rho$
becomes increasing, which proves that $\rho$ is lower bounded by a positive
constant. We get an upper bound for $|\dot{\rho}(t)|$ from the relation

$|\dot{\rho}(t)|\leqslant\left(\underset{x\in
K_{t}}{\max}\left(R(x,(I_{g}v)(t,x)\right)\right)\rho(t),$

and the boundedness of $\rho(t)$.
We introduce now the function

$\tilde{v}(t,y)=v(t,X_{t}(y))e^{\int_{0}^{t}\text{div }a(X_{s}(y))ds},$

which satisfies

$\int_{K}\tilde{v}(t,y)dy=\int_{K_{t}}v(t,x)dx=\rho(t).$

Moreover,

$\displaystyle\partial_{t}\tilde{v}(t,y)=$ $\displaystyle
R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\tilde{v}(t,y).$ (82)

Before proving that $\rho$ satisfies the equality of Lemma 11, we observe that
the following relation holds for all $y\in K$:

$\displaystyle\left\lvert\frac{d}{dt}(I_{g}v)(t,X_{t}(y))-\psi_{g}(\hat{x},\hat{x})\dot{\rho}(t)\right\rvert\leqslant$
$\displaystyle\left\lvert\int_{K_{t}}{a(X_{t}(y))\cdot\nabla_{x}\psi_{g}(X_{t}(y),z)v(t,z)dz}\right\rvert$
$\displaystyle+\left\rvert\int_{K_{t}}{\big{(}\psi_{g}(X_{t}(y),z)-\psi_{g}(\hat{x},\hat{x})\big{)}\partial_{t}v(t,z)dz}\right\rvert.$
(83)

From the equation satisfied by $v$, we see that

$\displaystyle\int_{K_{t}}{\big{(}\psi_{g}(X_{t}(y),z)-\psi_{g}(\hat{x},\hat{x})\big{)}\partial_{t}v(t,z)dz}$
$\displaystyle=$
$\displaystyle\int_{K_{t}}{\nabla_{z}\psi_{g}(X_{t}(y),z)a(z)v(t,z)dz}$
$\displaystyle+\int_{K_{t}}{\big{(}\psi_{g}(X_{t}(y),z)-\psi_{g}(\hat{x},\hat{x})\big{)}R(z,(I_{g}v)(t,z))v(t,z)dz}$
$\displaystyle=$
$\displaystyle\int_{K}{\nabla_{z}\psi_{g}(X_{t}(y),X_{t}(\bar{z}))a(X_{t}(\bar{z}))\tilde{v}(t,\bar{z})d\bar{z}}$
$\displaystyle+\int_{K}{\big{(}\psi_{g}(X_{t}(y),X_{t}(\bar{z}))-\psi_{g}(\hat{x},\hat{x})\big{)}R(X_{t}(\bar{z}),(I_{g}v)(t,X_{t}(\bar{z})))\tilde{v}(t,\bar{z})d\bar{z}},$

which allows us to conclude that

$\displaystyle\left\lvert\int_{K_{t}}{\big{(}\psi_{g}(X_{t}(y),z)-\psi_{g}(\hat{x},\hat{x})\big{)}\partial_{t}v(t,z)dz}\right\rvert\leqslant$
$\displaystyle\max\limits_{\bar{z}\in
K}\|a(X_{t}(\bar{z}))\|\|\psi_{g}\|_{W^{1,\infty}(\mathbb{R}^{d})}\rho(t)$
$\displaystyle+\max\limits_{y,\bar{z}\in
K}\left\lvert\psi_{g}(X_{t}(y),X_{t}(\bar{z}))-\psi_{g}(\hat{x},\hat{x})\right\rvert\overline{R}\rho(t).$

Using this relation in (83) gives

$\displaystyle\left\lvert\frac{d}{dt}(I_{\gamma}v)(t,X_{t}(y))-\psi_{\gamma}(\hat{x},\hat{x})\dot{\rho}(t)\right\rvert\leqslant$
$\displaystyle 2\max\limits_{\bar{z}\in
K}\|a(X_{t}(\bar{z}))\|\|\psi_{\gamma}\|_{W^{1,\infty}(\mathbb{R}^{d})}\rho(t)$
$\displaystyle+\max\limits_{y,\bar{z}\in
K}\left\lvert\psi_{\gamma}(X_{t}(y),X_{t}(\bar{z}))-\psi_{\gamma}(\hat{x},\hat{x})\right\rvert\overline{R}\rho(t),$

which proves, according to hypothesis (75), that for all $y\in K$

$\displaystyle\frac{d}{dt}I_{g}(t,X(t,y))=\psi_{g}(\hat{x},\hat{x})\dot{\rho}(t)+O(e^{-\delta
t})),$ (84)

for a certain $\delta>0$.
Differentiating (81) we get

$\displaystyle\ddot{\rho}(t)=$
$\displaystyle\frac{d}{dt}\left(\int_{K_{t}}{R(x,(I_{g}v)(t,x))v(t,x)}dx\right)$
$\displaystyle=$
$\displaystyle\frac{d}{dt}\left(\int_{K}{R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\tilde{v}(t,y)}dy\right)$
$\displaystyle=$
$\displaystyle\int_{K}{a(X_{t}(y))\cdot\nabla_{x}R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\tilde{v}(t,y)}dy$
$\displaystyle+\int_{K}{\frac{d}{dt}(I_{g}v)(t,X_{t}(y))\partial_{I}R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\tilde{v}(t,y)}dy$
$\displaystyle+\int_{K}{R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\partial_{t}\tilde{v}(t,y)dy}.$

Let us note that

  * _i)_

Since $a(X_{t}(y))$ converges uniformly to zero with an exponential speed,
then

$\int_{\mathbb{R}^{d}}{a(X_{t}(y))\cdot\nabla_{x}R\left(X_{t}(y),(I_{g}v)(t,X_{t}(y))\right)\tilde{v}(t,y))dy}=O(e^{-\delta
t}),$

thanks to the boundedness of $\nabla_{x}R$ and $\rho(t)$.

  * _ii)_

According to (84),

$\int_{K}{\frac{d}{dt}(I_{g}v)(t,X_{t}(y))\partial_{I}R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\tilde{v}(t,y)}dy=-p(t)\dot{\rho}(t)+O(e^{-\delta
t}),$

with

$\displaystyle p(t)$
$\displaystyle:=-\psi_{g}(\hat{x},\hat{x})\int_{K}{\partial_{I}R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\tilde{v}(t,y)}dy\geqslant\psi_{g}(\hat{x},\hat{x})D\min\limits_{t\geqslant
0}\rho(t)>0.$

  * _iii)_

Directly from (82),

$\displaystyle\int_{K}{R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\partial_{t}\tilde{v}(t,y)dy}$
$\displaystyle=\int_{K}{\Big{(}R(X_{t}(y),(I_{g}v)(t,X_{t}(y)))\Big{)}^{2}\tilde{v}(t,y)dy}$
$\displaystyle=:P(t)\geqslant 0.$

Thus, $\ddot{\rho}(t)\geqslant-p(t)\dot{\rho}(t)+O(e^{-\delta t})$, hence,
$\rho$ converges, thanks to Lemma 11.

Recalling (80), $v(t,\cdot)$ thus converges to $\hat{\rho}\delta_{\hat{x}}$,
where $\hat{\rho}$ is the limit of $\rho$. We conclude, according to
Proposition 8, that $\hat{\rho}$ satisfies the expected equality.
Having proved that $v$ and $v^{h}$ share the same limit is enough then to
conclude, thanks to Lemma 9, that $v_{\varepsilon}^{h}$ is an asymptotic
preserving approximation. ∎

If $m$ is not $0$, under very specific hypotheses over its support, we can
extend the result of Theorem 10. The explanation behind this is simple: as
long as the population is composed of traits that are not prone to
mutatations, it will evolve as in the case where mutations are not possible at
all.

######  Theorem 11.

Let us assume that there exists $\hat{x}\in\mathbb{R}^{d}$ which is an
asymptotically stable equilibrium for the ODE $\dot{x}=a(x)$ such that (75)
holds. Moreover, let us assume that $\textrm{supp}\left(v^{0}\right)\cup
K_{x}$ is a compact set such that there exist $C^{\prime},\delta^{\prime}>0$
such that

$\forall y\in\textrm{supp}\left(v^{0}\right)\cup K_{x},\quad t\geq
0,\quad\lVert X(t,y)-\hat{x}\rVert\leq C^{\prime}e^{-\delta^{\prime}t},$

that $\bigcup\limits_{s\geqslant 0}\left(X_{s}(\text{supp}(v^{0})\cup
K_{x})\right)\cap K_{y}=\emptyset$ and that hypothesis (76) holds. Then, $v$
converges to $\hat{\rho}\>\delta_{\hat{x}}$ in the weak sense in the space of
Radon measures, where $\hat{\rho}$ is the unique solution of

$R(\hat{x},\psi_{g}(\hat{x},\hat{x})\hat{\rho})=0.$

Consequently, $v^{h}_{\varepsilon}$ is an asymptotic preserving approximation
of $v$.

###### Proof.

By using the fact that, for any $x\in\mathbb{R}^{d}$,

$v(t,x)=v^{0}(X_{-t}(x))e^{\mathcal{G}(0,t,x)}+\int_{0}^{t}\int_{\mathbb{R}^{d}}m(X_{s-t}(x),z,(I_{d}v)(s,X_{s-t}(x)))v(s,z)dze^{\mathcal{G}(s,t,x)}ds,$

where

$\mathcal{G}(s,t,x):=\int_{s}^{t}{R(X_{\tau-t}(x),(I_{g}v)(\tau,X_{\tau-t}(x)))-\text{div
}a(X_{\tau-t}(x))d\tau},$

we observe that $\text{supp}(v(t,x))\subset
X_{t}(\text{supp}(v^{0}))\cup\bigcup\limits_{0\leqslant s\leqslant
t}X_{s}(K_{x})\subset\bigcup\limits_{s\geqslant
0}\left(X_{s}(\text{supp}(v^{0})\cup K_{x})\right)$, therefore

$\displaystyle v(t,x)=$ $\displaystyle
v^{0}(X_{-t}(x))e^{\mathcal{G}(0,t,x)}+\int_{0}^{t}\int_{\mathbb{R}^{d}}m(X_{s-t}(x),z,(I_{d}v)(s,X_{s-t}(x)))v(s,z)dze^{\mathcal{G}(s,t,x)}ds$
$\displaystyle=$ $\displaystyle
v^{0}(X_{-t}(x))e^{\mathcal{G}(0,t,x)}+\int_{0}^{t}\int_{\text{supp}(v(t,x))\cap
K_{y}}m(X_{s-t}(x),z,(I_{d}v)(s,X_{s-t}(x)))v(s,z)dze^{\mathcal{G}(s,t,x)}ds$
$\displaystyle=$ $\displaystyle v^{0}(X_{-t}(x))e^{\mathcal{G}(0,t,x)}.$