# A Hamilton-Jacobi approach to nonlocal kinetic equations

Nadia Loy Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di
Torino, Corso Duca degli Abruzzi 24, 10129, Torino, Italy,
<EMAIL_ADDRESS>Benoit Perthame Sorbonne Université, CNRS, Université de
Paris, Inria, Laboratoire Jacques-Louis Lions, F-75005 Paris, France,
<EMAIL_ADDRESS>
###### Abstract

Highly concentrated patterns have been observed in a spatially heterogeneous,
nonlocal, model of BGK type implementing a velocity-jump process. We study
both a linear and a nonlinear case and describe the concentration profile. In
particular, we analyse a hyperbolic (or high frequency) regime that can be
interpreted both as a local (microscopic) or as a nonlocal (macroscopic)
rescaling. We consider a Hopf-Cole transform and derive a Hamilton-Jacobi
equation. The concentrations are then explained as a consequence of the
stationary points of the Hamiltonian that is spatially heterogeneous like the
velocity-jump process. After revising the classical hydrodynamic limits for
the aggregate quantities and the eikonal equation that can be derived from
those with a Hopf-Cole transform, we find that the Hamilton-Jacobi equation is
a second order approximation of the eikonal equation in the limit of small
diffusivity. For nonlinear turning kernels, the Hopf-Cole transform allows to
study the stability of the possible homogeneous configurations and of patterns
and the results of a linear stability analysis previously obtained are found
and extended to a nonlinear regime. In particular, it is shown that
instability (pattern formation) occurs when the Hamiltonian is convex-concave.

## Introduction

Kinetic equations have proved to be an effective mathematical framework for
modeling cell migration, both for bacteria [35, 10, 22, 34, 12, 19, 20] and
for cells in a tissue [21, 13, 30, 31, 32]. In fact, the typical migration
mode of a cell is the run and tumble, consisting in alternating runs over
straight lines and reorientations, that may be biased by the presence of
external signals affecting the choice of the direction, such as, for example,
chemicals and the cell population density itself. At the population, or
aggregate, level, the latter process may give rise to a tactic dynamics such
as chemotaxis and adhesion, respectively. The run and tumble process may be
modelled as a microscopic stochastic process named velocity-jump process [37].
It is a Markovian processes that prescribes a transition probability $T$ of
choosing a new velocity and a frequency of reorientation $\mu$. In particular,
the transition probability $T=T[\mathcal{S}]$ may be influenced by the
presence of an external signal $\mathcal{S}$, that may embody the presence of
a chemoattractant or of the cell population density. The kinetic equation that
implements a velocity-jump process of intensity $\mu>0$, that is a piecewise
deterministic Markov process in which we consider a transition probability
$T[\mathcal{S}]$, may be written as

$\partial_{t}f(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}f(t,x,v,\hat{v})=\mu\big{(}\rho(t,x)T[\mathcal{S}](v,\hat{v})-f(t,x,v,\hat{v})\big{)},$
(1)

where $f=f(t,x,v,\hat{v})$ is a probability density function describing the
distribution of the particle located at position
$x\in\Omega\in\mathbb{R}^{d}$, moving with speed $v\in[0,U]$ along direction
$\hat{v}\in\mathbb{S}^{d-1}$, for each time $t>0$. We also use the notation
${\bf v}$ for the microscopic velocity of the cells that is given by the
vector ${\bf v}=v\hat{v}$. As
$(v,\hat{v})\in\mathcal{V}:=[0,U]\times\mathbb{S}^{d-1}$, then ${\bf v}\in
B(0,U)$, that is compact in $\mathbb{R}^{d}$ and symmetric. The function
$\rho(t,x)$ denotes the number density of cells in position $x$ at time $t$:

$\rho(t,x)=\int_{0}^{U}\int_{\mathbb{S}^{d-1}}f(t,x,v,\hat{v})\,dvd\hat{v},\qquad(t,x)\in\mathbb{R}_{+}\times\Omega.$

Formally, (1) is a kinetic equation with linear relaxation operator of BGK
type.

Another important issue in modeling cell migration is the nonvanishing size of
the cell that gives rise to nonlocality in the physical space (see e.g. [1,
14] and references therein). Specifically at the kinetic level, a nonlocal
gradient of the chemoattractant $\mathcal{S}$ sensed with a sampling radius
has been introduced in [34, 23]. In [31, 32, 16, 15] the authors propose some
models in which $T$ depends on a fixed external signal $\mathcal{S}$ and also
consider the case in which $\mathcal{S}$ is the cell density $\rho$, thus
mimicking adhesion. In [31, 32] the authors derive the so called macroscopic
models for the aggregate quantities defined as averaged quantities (the
statistical moments of $f$), showing that keeping the nonlocality at the
aggregate description implies a strong nonvanishing advection term. In such
models, aggregation and concentrations have been observed, both in the case of
linear models (i.e. $T$ depends on $\mathcal{S}$), and nonlinear ones (i.e.
$T$ depends on $\rho$). In particular, in [33] a stability analysis of a model
with adhesion is performed and pattern formation is shown.

Pattern formation may be seen in this context as a formation of small
concentrations, typically persisting in time. As a matter of fact, this kind
of solution may be represented as the sum of Dirac masses, in the form

$\rho_{\varepsilon}(t,x)\approx\sum\bar{\rho}_{i}(t)\delta(x-\bar{x}_{i}(t)),$

where $\varepsilon$ is a small parameter determining a specific regime,
typically of high frequencies, in which there may be formation of
concentrations and patterns, $\bar{x}_{i}$ is the location of a concentration
point (that is a maximum point of the solution $\rho_{\varepsilon}$) and
$\bar{\rho}_{i}$ is the weight of the concentration. A popular tool for
analysing such concentration profiles is the real WKB ansatz or Hopf-Cole
transform. The leading idea is similar to approximate the Dirac mass by
Gaussians

$\delta(x-\bar{x})\approx\dfrac{1}{\sqrt{2\pi\varepsilon}}\exp^{\dfrac{-|x-\bar{x}|^{2}}{2\varepsilon}}=\exp^{\dfrac{-|x-\bar{x}|^{2}-\varepsilon
2\pi\ln(\varepsilon)}{2\varepsilon}}.$

Therefore, the assumption is to consider a real WKB ansatz for the high
frequency regime (or Hopf-Cole transform) in the form

$\rho_{\varepsilon}(t,x)=\exp^{-\dfrac{\varphi_{\varepsilon}(t,x)}{\varepsilon}},\qquad\lim_{\varepsilon\to
0}\varphi_{\varepsilon}\geq 0,$

where the concentration points are understood as the minima of the phase
$\varphi_{\varepsilon}$. This kind of analysis typically leads, in the limit
$\varepsilon\rightarrow 0$, to a constrained Hamilton-Jacobi equation for the
phase. Such an Hamilton-Jacobi equation enables a rigorous derivation of the
so-called canonical equation for the evolution of the maxima, which has been
formally proposed in the framework of adaptive dynamics in structured (by a
trait) population equations to describe the trait evolution. In the context of
adaptive dynamics, the maximum is interpreted as the ‘fittest trait’, [3, 29,
26, 28] and the Hamilton-Jacobi equation also allows to find an equation for
the location of the maxima, that is typically given by the roots of the
equation obtained by setting the growth parameter to zero.

In the context of kinetic equations, this kind of ansatz has been first
adopted by Bouin and Calvez in [6], where they study a BGK model with a
relaxation Maxwellian that is spatially homogeneous. In a high frequency
regime, they assume a Hopf-Cole transform of the distribution

$f_{\varepsilon}(t,x,v,\hat{v})=\exp^{-\dfrac{\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}$
(2)

and also assume a perturbed test function
$\phi_{\varepsilon}(t,x,v,\hat{v})=\phi(t,x)+\varepsilon\eta(t,x,v,\hat{v}),$
i.e. the leading order in $\varepsilon$ only depends on $(t,x)$. The authors
in [6] derive in the limit $\varepsilon\rightarrow 0^{+}$ a Hamilton-Jacobi
equation to which the potential is a viscosity solution. This has been
extended in [9] in order to take into account higher space dimensions, where
the Hamiltonian may lack $\mathcal{C}^{1}$ regularity. In [4, 8, 5] the
authors study front propagation in transport-reaction kinetic equations. In
[7] the authors study a BGK type equation in a high frequency regime and with
a Maxwellian with vanishing variance and derive a new constrained nonlocal
Hamilton-Jacobi equation.

In the present work, we want to analyse nonlocal kinetic equations of the same
class of models as introduced in [31, 33] with the real WKB ansatz and the
perturbed test function. The leading idea is to explain the concentration
profiles by means of an appropriate Hamilton-Jacobi equation. The latter is
derived from the kinetic equation after assuming a Hopf-Cole transform in the
same spirit as in [6] and interpreting the location of the maxima (that is a
spatial variable) as the fittest trait, like in adaptive dynamics.

In the first section, we will present the nonlocal kinetic equation under
study along with the evolution equations for the aggregate quantities that can
be derived in a regime of high frequencies and that can be obtained in a local
and in a nonlocal rescaling. In section 3, we shall perform the WKB analysis
of a linear nonlocal equation and present the concentration result as well as
the canonical equation for the evolution of the maxima and some examples along
with some numerical test. In section 4 we shall extend formally the analysis
to a nonlinear case. In section 5 we draw some conclusion.

##  1 Preliminaries

###  1.1 A nonlocal kinetic equation

In the same spirit as [31], we consider the kinetic equation (1) with
transition probability that depends on the external signal
$\mathcal{S}:\,\Omega\rightarrow\mathbb{R}_{+}$, that is measured nonlocally
in the physical space and that affects the choice of the direction of the
cells. The transition probability $T[\mathcal{S}](v,\hat{v})$, that is a
probability on $\mathcal{V}$ and depends on $x$ through $\mathcal{S}$, is in
general defined by

$T[\mathcal{S}](v,\hat{v})=c(x)\psi(v|\hat{v})b(\mathcal{S}(x+R\hat{v})),$ (3)

where $b(\cdot)$ is a function that weights the external field $\mathcal{S}$,
while the quantity $R$ is the sensing radius defining the neighborhood of the
particle where the field $\mathcal{S}$ is measured. The function
$\psi=\psi(v|\hat{v}):[0,U]\rightarrow\mathbb{R}_{+}$ is the probability
density function of the possible speeds $v\in[0,U]$ on a given direction
$\hat{v}$, satisfying

$\int_{0}^{U}\psi(v|\hat{v})\,dv=1\quad\forall\hat{v}\in\mathbb{S}^{d-1}.$

We denote its average speed (along direction $\hat{v}$) $V_{\psi}$ and second
statistical moment $D^{2}$ (that we assume to be independent of the
direction), i.e. they are defined by

$V_{\psi}(\hat{v}):=\int_{0}^{U}\psi(v|\hat{v})v\,dv,\qquad
D^{2}(\hat{v}):=\int_{0}^{U}\psi(v|\hat{v})v^{2}\,dv.$ (4)

The function $c(x)$ is a normalization function defined by

$c(x)^{-1}:=\int_{\mathbb{S}^{d-1}}b(\mathcal{S}(x+R\hat{v}))\,d\hat{v}.$

This ensures that $T[\mathcal{S}](v,\hat{v})$ is a probability density
function on $\mathcal{V}=[0,U]\times\mathbb{S}^{d-1}$ as

$\int_{0}^{U}\int_{\mathbb{S}^{d-1}}T[\mathcal{S}](v,\hat{v})\,d\hat{v}dv=1,$

in such a way that the number density is conserved at $(t,x)$. Then, we can
also define the average velocity of the transition probability
$T[\mathcal{S}]$ as

${\bf
U}_{\mathcal{S}}(x)=\int_{0}^{U}\int_{\mathbb{S}^{d-1}}T[\mathcal{S}](v,\hat{v}){\bf
v}\,d\hat{v}dv,$ (5)

and its variance-covariance matrix

$\mathbb{D}_{\mathcal{S}}(x)=\int_{0}^{U}\int_{\mathbb{S}^{d-1}}T[\mathcal{S}](v,\hat{v})({\bf
v}-{\bf U}_{\mathcal{S}})\otimes({\bf v}-{\bf U}_{\mathcal{S}})\,d\hat{v}dv.$
(6)

We remark that, when $\Omega$ is bounded, in order to deal with the boundary,
we must restrict the sensing radius using the formula

$R(x,\hat{v}):=\min\\{\lambda\in[0,R]:x+\lambda\hat{v}\in\Omega\\}.$ (7)

Eq. (1) needs to be coupled with initial and boundary conditions, defined by,
respectively

$f(0,x,v,\hat{v})=f^{0}(x,v,\hat{v}),\qquad(x,v,\hat{v})\in\Omega\times[0,U]\times\mathbb{S}^{d-1},$
(8) $f(t,x,v,\hat{v})=\mathcal{R}[f_{|\Gamma_{+}}](t,x,v,\hat{v})\qquad
x\in\partial\Omega,\,v\in[0,U],\,\hat{v}\in\Gamma_{-}(x),$ (9)

where

$\Gamma_{\pm}(x):=\\{\hat{v}\in\mathbb{S}^{d-1}:\hat{v}\cdot\boldsymbol{n}(x)\gtrless
0\\},$

with ${\boldsymbol{n}}(x)$ the outward normal to the boundary $\partial\Omega$
in the point $x$. As boundary conditions for the kinetic equation, we assume a
standard diffusive boundary condition [27, 36] called Maxwellian boundary
conditions, defined as

$\mathcal{R}[f_{+}](t,x,v,\hat{v})=\alpha(x)f(t,x,v,\mathcal{W}(\hat{v}))+(1-\alpha(x))M(x,v,\hat{v})\int_{0}^{U}\int_{\hat{v}^{*}\cdot\boldsymbol{n}\geq
0}f(t,x,v^{*},\hat{v}^{*})|\hat{v}^{*}\cdot\boldsymbol{n}|d\hat{v}^{*}dv^{*},$
(10)

where $\mathcal{W}(\hat{v})=-\hat{v}$ for the bounce back reflection condition
and
$\mathcal{W}(\hat{v})=\hat{v}-2(\hat{v}\cdot{\boldsymbol{n}}){\boldsymbol{n}}$
for the specular reflection. Diffusive boundary conditions are no-flux
boundary conditions at the macroscopic level [36], in the sense that the total
mass is conserved in $\Omega$. In fact, it may be proved that [36] if $f$ is a
solution to (1)-(8)-(9)-(10), then the following is satisfied

$\int_{0}^{U}\int_{\mathbb{S}^{d-1}}f(t,x,v,\hat{v}){\bf v}\cdot{\bf
n}(x)d\hat{v}\,dv=0,\quad\forall x\in\partial\Omega,\quad t>0.$ (11)

The equilibrium distribution of (1) is given by

$f^{\infty}(x,v,\hat{v})=\rho(x)T[\mathcal{S}](v,\hat{v}).$

With classical arguments (Jensen’s inequality) and assuming Maxwellian
boundary conditions (that are nonabsorbing boundary conditions [11]) it is
easy to see that given a convex function $\Phi$, then

$\dfrac{d}{dt}\int_{\Omega}\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\Phi\Big{(}f(t,x,v,\hat{v})\Big{)}\,dvd\hat{v}dx\leq
0,$

and the equality holds if and only if $f=f^{\infty}$. This equilibrium is
asymptotically stable and does not depend on the initial condition. As
$T[\mathcal{S}]$ depends on $x$ through $\mathcal{S}$, in order for
$f^{\infty}$ to be a stationary equilibrium, then the following must be
satisfied

$-\dfrac{{\bf v}\cdot\nabla_{x}\rho}{\rho}=\dfrac{{\bf
v}\cdot\nabla_{x}T[\mathcal{S}]}{T[\mathcal{S}]}.$ (12)

###  1.2 Rescaling

We now consider a regime in which reorientations occur at random exponential
times with rate $\mu=\dfrac{1}{\varepsilon}$, i.e., the dynamics is ruled by

$\partial_{t}f_{\varepsilon}(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}f_{\varepsilon}(t,x,v,\hat{v})=\dfrac{1}{\varepsilon}\left(\rho_{\varepsilon}(t,x)T[S]_{\varepsilon}(v,\hat{v})-f_{\varepsilon}(t,x,v,\hat{v})\right)$
(13)

where the limit $\varepsilon\rightarrow 0$ defines a high frequency regime.

On a one hand, a possible interpretation is to see Eq. (13) as the result of a
hyperbolic scaling of Eq. (1) defined by

$(t,v,x)\rightarrow(\dfrac{t}{\varepsilon},v,\dfrac{x}{\varepsilon}),$ (14)

that defines a long time scale (the equilibrium is reached fast) and a
macroscopic (or large) space scale in which the interactions are localized. In
fact, in this rescaling we also need to consider a scaling of the sensing
radius $R$, i.e.

$R\rightarrow\varepsilon R,$

that naturally leads to a localization of the interactions of the cells with
the background $\mathcal{S}$. In this sense the large scale limit of (13) for
$\varepsilon\rightarrow 0$ leads to the hydrodynamic (fluid) behavior of the
system on a macroscopic space scale that must be observed on a long time
scale.

However, a priori we can consider the perspective of the following
nondimensionalization

$x\rightarrow\dfrac{x}{L},\quad t\rightarrow\dfrac{t}{t_{0}},\quad
v\rightarrow\dfrac{v}{V},\quad\rho\rightarrow\dfrac{\rho}{\bar{\rho}},\quad
f\rightarrow\dfrac{f}{\bar{\rho}/V^{d}},\quad
T[\mathcal{S}]\rightarrow\dfrac{T[\mathcal{S}]}{V^{d}},$ (15)

where $t_{0}$ and $L$ are characteristic time and length scales of the system,
$V$ is the typical speed, while $\bar{\rho}$ is a reference density. Plugging
(15) in (1) we obtain

$\textrm{St}\,\partial_{t}f(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}f(t,x,v,\hat{v})=\dfrac{1}{\textrm{Kn}}\left(\rho(t,x)T[S](v,\hat{v})-f(t,x,v,\hat{v})\right),$
(16)

where the kinetic Strouhal number St and Knudsen number Kn are defined as

$\textrm{St}:=\dfrac{L}{Vt_{0}},\qquad\qquad\textrm{Kn}:=\dfrac{V}{L\mu}.$

The regime under consideration in (13) corresponds to having parameters given
by

$\textrm{St}\sim\mathcal{O}(1),\qquad\qquad\textrm{Kn}\sim\mathcal{O}(\varepsilon)\ll
1.$ (17)

By looking at (16) we observe that the parameters regime (17) corresponds to a
large time horizon $t_{0}$ satisfying

$t_{0}\sim\dfrac{\varepsilon^{-1}}{\mu}.$ (18)

This may be rephrased saying that we choose $\varepsilon$ such that

$\dfrac{V}{L}\sim\mathcal{O}(\varepsilon)$ (19)

and we choose a drift long time scale as

$t_{0}=\dfrac{L}{V},$

that satisfies (18) because of (19). We remark that if $V=\mathcal{O}(1)$,
then $\dfrac{1}{L}\sim\varepsilon$, that amounts to (14) (where we use again
$R\rightarrow\varepsilon R$).

On the other hand, if $L=\mathcal{O}(1)$, i.e., we observe the dynamics on the
microscopic space scale, then (19) amounts to a regime of very small speeds
$V\sim\varepsilon$ that must be observed, in order to balance the smallness of
the speed, on a long time scale. This can be seen as a scaling in the form

$(t,v,x)\rightarrow(\dfrac{t}{\varepsilon},\varepsilon v,x).$ (20)

The latter may also be seen as a nonlocal regime as the sensing radius $R$ is
not rescaled.

###  1.3 Aggregate behaviour

In the regime defined by (13), we may find limiting equations for the averaged
population quantities. The single conservation law induces that the aggregate
quantity is the mass and we obtain an evolution equations for the number
density $\rho$. Taking into account the equilibrium
$T[\mathcal{S}]_{\varepsilon}$, the formal expansion of (13) at order
$\mathcal{O}(\varepsilon)$ is a diffusion-advection equation with a dominating
drift term and small diffusivity

$\partial_{t}\rho_{\varepsilon}+\nabla_{x}\cdot(\rho_{\varepsilon}{\bf
U}_{\mathcal{S}}^{\varepsilon})=\varepsilon\nabla_{x}\cdot\left(\nabla_{x}\cdot(\mathbb{D}_{\mathcal{S}}^{\varepsilon}\rho_{\varepsilon})+\rho_{\varepsilon}{\bf
U}_{\mathcal{S}}^{\varepsilon}\nabla_{x}\cdot{\bf
U}_{\mathcal{S}}^{\varepsilon}\right),$ (21)

where ${\bf U}_{\mathcal{S}}^{\varepsilon}$ is the average of
$T[\mathcal{S}]_{\varepsilon}$ and $\mathbb{D}_{\mathcal{S}}^{\varepsilon}$
its variance-covariance matrix as defined by (5) and (6). The boundary
conditions can be found by imposing (11) to $f_{\varepsilon}$ [36] and this
gains

$(\rho_{\varepsilon}{\bf
U}_{\mathcal{S}}^{\varepsilon}-\varepsilon\big{(}(\mathbb{D}_{\mathcal{S}}^{\varepsilon}\nabla\rho_{\varepsilon})+\rho_{\varepsilon}(\nabla\cdot\mathbb{D}_{\mathcal{S}}^{\varepsilon}+{\bf
U}_{\mathcal{S}}^{\varepsilon}\nabla\cdot{\bf
U}_{\mathcal{S}}^{\varepsilon})\big{)}\cdot\boldsymbol{n}=0\qquad\text{on
}\;\partial\Omega.$ (22)

In the asymptotic limit $\varepsilon\rightarrow 0^{+}$, the dynamics is ruled
by the equilibrium of (13) at order zero in $\varepsilon$, that is defined by

$T[\mathcal{S}]_{0}:=\lim_{\varepsilon\rightarrow
0^{+}}T[\mathcal{S}]_{\varepsilon},$

in such a way that the evolution equation for $\rho=\rho_{0}$ is

$\partial_{t}\rho+\nabla_{x}\cdot(\rho{\bf U}_{\mathcal{S}}^{0})=0.$ (23)

The boundary conditions, that can be derived substituting $f_{\varepsilon}$ in
(11) and letting $\varepsilon\rightarrow 0$, are given by

$\rho{\bf U}_{\mathcal{S}}^{0}(x)\cdot\boldsymbol{n}(x)=0,\quad\text{for
$x\in\partial\Omega$},$ (24)

that are actually no-flux conditions for the conservation law (23). As
explained in Appendix Appendix. Boundary conditions, the latter are actually
only to be imposed on the entering region, but this can be in fact derived by
the underlying kinetic boundary conditions that are imposed on the entering
zone $\Gamma_{-}(x)$.

In conclusion, both rescalings, the local macroscopic one (14) or the nonlocal
microscopic one (20), are possible and lead to the same equation (23), but the
scales of the sensing radius differ and, as a consequence, ${\bf
U}_{\mathcal{S}}^{0}$ differs. For example, in the case
$b(\mathcal{S})=\mathcal{S}$, when the rescaling (14) is performed (assuming
$\mathcal{S}$ smooth), then
$T[\mathcal{S}]_{\varepsilon}(v,\hat{v})=c(x)\psi(v|\hat{v})\mathcal{S}_{\varepsilon}(x+\varepsilon
R\hat{v})$, so that

$T[\mathcal{S}]_{0}(v,\hat{v})=\psi(v|\hat{v})\dfrac{1}{|\mathbb{S}^{d-1}|},\quad{\bf
U}_{\mathcal{S}}^{0}=\int_{\mathbb{S}^{d-1}}V_{\psi}(\hat{v})\,\hat{v}d\hat{v},\quad\mathbb{D}_{\mathcal{S}}^{0}=D^{2}\mathbb{I}.$

If $V_{\psi}(\hat{v})$ is even (e.g., $V_{\psi}$ constant), then ${\bf
U}_{\mathcal{S}}^{0}=0$, so that from (21) $\partial_{t}\rho_{0}=0$. Anyway,
localization does not imply a vanishing drift in all cases. Let us, for
example, consider a comparative sensing [31], which means that we assume that
the turning rate depends on what is measured in $x+R\hat{v}$ and $x-R\hat{v}$,
i.e.

$b\big{(}\mathcal{S}(x+\lambda\hat{v}),\mathcal{S}(x-\lambda\hat{v})\big{)}=\alpha+\beta\dfrac{\mathcal{S}(x-\lambda\hat{v})-\mathcal{S}(x+\lambda\hat{v})}{2k+\mathcal{S}(x+\lambda\hat{v})+\mathcal{S}(x-\lambda\hat{v})}\,.$

Then, assuming fast adaptation $\beta\rightarrow\frac{\beta}{\varepsilon}$, we
find

$T[\mathcal{S}]_{0}(v,\hat{v})=\psi(v|\hat{v})\left(\alpha+\beta\dfrac{R\nabla\mathcal{S}(x)\cdot\hat{v}}{k+\mathcal{S}(x)}\right)\,,$

which means that, even in the localized interactions regime, if
$T[\mathcal{S}]$ is given by a comparative sensing, the equilibrium depends on
the directional derivative of the external field $\mathcal{S}$ along each
microscopic direction $\hat{v}$. Conversely, in the regime (20) we have,

$T[\mathcal{S}]_{0}=T[\mathcal{S}],\quad\text{and thus,}\quad{\bf
U}_{\mathcal{S}}^{0}=c(x)\int_{\mathbb{S}^{d-1}}V_{\psi}(\hat{v})\mathcal{S}(x+R\hat{v})\,\hat{v}d\hat{v},$

that is in general a nonvanishing quantity at the microscopic space scale. In
fact, we can remark that even in the case $V_{\psi}(\hat{v})=V_{\psi}$, then

${\bf
U}_{\mathcal{S}}^{0}=V_{\psi}c(x)\int_{\mathbb{S}^{d-1}}\mathcal{S}(x+R\hat{v})\,\hat{v}d\hat{v}$

is a nonvanishing quantity unless $\mathcal{S}$ is spatially homogeneous.
Therefore, the dominating drift term is due to the spatial heterogeneity that
is sensed nonlocally.

In conclusion, Eq. (23), if derived as a large scale limit in the regime (14),
has to be meant as a hydrodynamic limit on a macroscopic space scale in which
interactions are localized and the (now local) equilibrium is reached fast as
a longtime scale is observed. Conversely, it is derived in the regime (20), it
implies a high frequency and small speeds regime on the microscopic space
scale and slow time scale that is the same as the one of the original kinetic
equation (1). Therefore, when derived in this regime, Eq. (23) describes the
evolution of the average number density on the original (microscopic) phase
space.

####  1.3.1 Diffusive limit

When dealing with cell migration modeling, a typical rescaling is the
diffusive one. In the present framework, it corresponds to choosing in the
nondimensionalization a diffusive long time scale $t_{0}=\dfrac{L^{2}}{V^{2}}$
that actually satisfies $t_{0}=\dfrac{\varepsilon^{-2}}{\mu}$, in such a way
that $\textrm{St}=\varepsilon$. Therefore, the rescaled kinetic equation is in
the form

$\partial_{t}f_{\varepsilon}+\dfrac{1}{\varepsilon}{\bf
v}\cdot\nabla_{x}f_{\varepsilon}=\dfrac{1}{\varepsilon^{2}}(\rho_{\varepsilon}T[\mathcal{S}]_{\varepsilon}-f_{\varepsilon}).$
(25)

When we consider a nonlocal diffusive rescaling, i.e.

$(t,v,x)\rightarrow(\dfrac{t}{\varepsilon^{2}},\varepsilon v,x)$ (26)

and $T[\mathcal{S}]$ depends on $x$ and on $(v,\hat{v})$, then typically ${\bf
U}_{\mathcal{S}}$ is not a vanishing quantity and the aggregate equation for
$\rho$ is

$\partial_{t}\rho+\nabla_{x}\cdot({\bf
U}_{\mathcal{S}}\rho)=\nabla\cdot\nabla\cdot\left(\rho\mathbb{D}_{\mathcal{S}}\right).$
(27)

Conversely, when we consider a localized diffusive scaling, i.e.

$(t,v,x)\rightarrow(\dfrac{t}{\varepsilon^{2}},v,\dfrac{x}{\varepsilon})$ (28)

we can typically consider a Hilbert expansion for $\rho_{\varepsilon}$ and
$T[\mathcal{S}]_{\varepsilon}$, i.e.

$T[\mathcal{S}]_{\varepsilon}=T[\mathcal{S}]_{0}+\varepsilon
T[\mathcal{S}]_{1},\qquad\int_{\mathcal{V}}T[\mathcal{S}]_{0}dvd\hat{v}=1,\quad\int_{\mathcal{V}}T[\mathcal{S}]_{1}dvd\hat{v}=0.$

When the solvability condition ${\bf U}_{\mathcal{S}}^{0}=0$ is met, the
macroscopic equation for $\rho=\rho_{0}$ is

$\partial_{t}\rho+\nabla_{x}\cdot({\bf
U}_{\mathcal{S}}^{1}\rho)=\nabla\cdot\nabla\cdot\left(\rho\mathbb{D}_{\mathcal{S}}^{0}\right).$
(29)

Then, supposing that

$\mathbb{D}_{\mathcal{S}}^{0}=\nu^{2}\bar{\mathbb{D}}_{\mathcal{S}}^{0},$ (30)

when $\nu^{2}=\varepsilon$, i.e. for a small diffusivity, we essentially
recover (21) in the macroscopic limit.

###  1.4 Limit for small $R$

Let us define the characteristic length of variation of $\mathcal{S}$ as

$l_{\mathcal{S}}:=\dfrac{1}{\max\dfrac{|\nabla_{x}\mathcal{S}|}{\mathcal{S}}}.$
(31)

We remark that when $R\ll l_{\mathcal{S}}$, then we may consider the Taylor
expansion of $\mathcal{S}$ at first order

$\mathcal{S}(x+R\hat{v})=\mathcal{S}(x)+R\hat{v}\cdot\nabla_{x}\mathcal{S}(x)+\mathcal{O}(R^{2}),$

that is a positive quantity. Then, for example, in the case
$b(\mathcal{S})=\mathcal{S}$, we may approximate the probability density
function $T[\mathcal{S}]$ as

$T[\mathcal{S}](v,\hat{v})=\dfrac{\psi(v|\hat{v})}{|\mathbb{S}^{d-1}|}\left[1+R\hat{v}\cdot\dfrac{\nabla_{x}\mathcal{S}}{\mathcal{S}}\right].$

Then, choosing $L=l_{\mathcal{S}}$ the nondimensionalization of (1) leads to

$\textrm{St}\,\partial_{t}f_{\varepsilon}(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}f_{\varepsilon}(t,x,v,\hat{v})=\dfrac{1}{\textrm{Kn}}\left(\rho_{\varepsilon}(t,x)\dfrac{\psi(v|\hat{v})}{|\mathbb{S}^{d-1}|}(1+\eta\hat{v}\cdot\dfrac{\nabla_{x}\mathcal{S}(x)}{\mathcal{S}(x)})-f_{\varepsilon}(t,x,v,\hat{v})\right),$

where

$\eta:=\dfrac{R}{l_{\mathcal{S}}}.$

With the choice $\dfrac{1}{L}\sim\varepsilon$, that amounts to (14), then we
obtain $\eta=R\varepsilon$ and

$\partial_{t}f_{\varepsilon}(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}f_{\varepsilon}(t,x,v,\hat{v})=\dfrac{1}{\varepsilon}\left(\rho_{\varepsilon}(t,x)(T[S]_{0}+\varepsilon
T[\mathcal{S}]_{1})-f_{\varepsilon}(t,x,v,\hat{v})\right),$

where
$T[\mathcal{S}]_{0}(v,\hat{v})=\dfrac{\psi(v|\hat{v})}{|\mathbb{S}^{d-1}|}$
and
$T[\mathcal{S}]_{1}(v,\hat{v})=\dfrac{\psi(v|\hat{v})}{|\mathbb{S}^{d-1}|}R\hat{v}\cdot\dfrac{\nabla_{x}\mathcal{S}(x)}{\mathcal{S}(x)}$.
This corresponds to considering (1) with (14), as the leading order term is
local while the role of the sensing radius enters the dynamics as a higher
order term. In this case the evolution equations for $\rho_{\varepsilon}$
correspond in the macroscopic point of view (14) and in the high frequency
(microscopic) one (20), i.e.,

$\partial_{t}\rho_{\varepsilon}=\varepsilon\Delta\rho_{\varepsilon}.$ (32)

When choosing the diffusive scaling (choosing $V_{\psi}$ constant for
simplicity) we obtain (29) with $\mathbb{D}_{\mathcal{S}}^{0}=\mathbb{I}$ and
${\bf U}_{\mathcal{S}}^{1}=R\dfrac{\nabla_{x}\mathcal{S}}{\mathcal{S}}$, that
is the Keller and Segel model [31, 24].

##  2 Concentration profile and the Hamilton–Jacobi equation

We want to study the concentration profile of the solution $f_{\varepsilon}$
of Eq. (13) by studying the equation for a potential $\varphi_{\varepsilon}$
obtained through the Hopf-Cole transform (2). We expect that
$\varphi_{\varepsilon}(t,x,v,\hat{v})$ behaves like a quadratic and thus that
$f_{\varepsilon}(t,x,v,\hat{v})$ behaves like a Dirac mass near each
concentration point. For that reason, we study the limit of
$\varphi_{\varepsilon}$.

###  2.1 The Hamilton–Jacobi equation

At first we remark that, from Eq. (13), $\varphi_{\varepsilon}$ satisfies the
equation

$\partial_{t}\varphi_{\varepsilon}+{\bf
v}\cdot\nabla_{x}\varphi_{\varepsilon}=1-T[S]_{\varepsilon}(v,\hat{v})\int_{\mathcal{V}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})+\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}\,dwd\hat{w}.$
(33)

From this we get

$\left(1-\partial_{t}\varphi_{\varepsilon}-{\bf
v}\cdot\nabla_{x}\varphi_{\varepsilon}\right)=T[S]_{\varepsilon}(v,\hat{v})\int_{\mathcal{V}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})+\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}\,dwd\hat{w}.$
(34)

Following [6, 26], we may also look for $f_{\varepsilon}$ under the form

$f_{\varepsilon}(t,x,v,\hat{v})=Q_{\varepsilon}(t,x,v,\hat{v})\exp^{-\dfrac{\tilde{\varphi}_{\varepsilon}(t,x)}{\varepsilon}},\qquad\varphi_{\varepsilon}(t,x,v,\hat{v})=\tilde{\varphi}_{\varepsilon}(t,x)-\varepsilon\ln
Q_{\varepsilon}(t,x,v,\hat{v}),$ (35)

with $\tilde{\varphi}_{\varepsilon}$ and $Q_{\varepsilon}$ to be determined.
Setting

$p_{\varepsilon}=\nabla_{x}\tilde{\varphi}_{\varepsilon},\qquad
H_{\varepsilon}=-\partial_{t}\tilde{\varphi}_{\varepsilon},$ (36)

we can write Eq. (34) as

$\begin{array}[]{ll}\varepsilon[\partial_{t}Q_{\varepsilon}(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}Q_{\varepsilon}(t,x,v,\hat{v})]+\\!&(1+H_{\varepsilon}-v\hat{v}\cdot
p_{\varepsilon})Q_{\varepsilon}(t,x,v,\hat{v})\\\\[5.0pt] &=\displaystyle
T[\mathcal{S}]_{\varepsilon}(v,\hat{v})\int_{\mathcal{V}}Q_{\varepsilon}(t,x,w,\hat{w})d\hat{w}dw.\end{array}$
(37)

The formal limit as $\varepsilon\to 0$ gives us

$(1+H-v\hat{v}\cdot p)Q(t,x,v,\hat{v})=\displaystyle
T[\mathcal{S}]_{0}(v,\hat{v})\int_{\mathcal{V}}Q(t,x,w,\hat{w})d\hat{w}dw.$

This can be interpreted as the eigenvalue-eigenfunction problem in
$(v,\hat{v})$, with $(t,x)$ parameters, which is to find $(H,\mathcal{Q})$
such that

$(1+H-v\hat{v}\cdot
p)\mathcal{Q}(x,p,v,\hat{v})=T[\mathcal{S}]_{0}(v,\hat{v})\int_{\mathcal{V}}\mathcal{Q}(x,p,w,\hat{w})d\hat{w}dw.$
(38)

Thanks to the Krein-Rutman theory, see [25], with good properties of
$T[\mathcal{S}]_{0}$ to be discussed later, this eigenproblem has a unique
solution once normalized as

$\mathcal{Q}(x,p,v,\hat{v})>0,\qquad\int_{\mathcal{V}}\mathcal{Q}(x,p,v,\hat{v})d\hat{v}dv=1,\qquad\forall
p,\;x.$ (39)

The eigenvalue $H$ is solely determined by the parameters $p$ and
$\mathcal{S}(x)$ and we can write $H=H(x,p)$ which provides us with the
Hamilton-Jacobi equation for the dominant term in (36)

$\partial_{t}\varphi+H(x,\nabla_{x}\varphi)=0,$ (40)

with $\varphi$ the (formal) common limit of $\varphi_{\varepsilon}$ or
$\tilde{\varphi}_{\varepsilon}$

$\tilde{\varphi}:=\lim_{\varepsilon\rightarrow
0^{+}}\tilde{\varphi}_{\varepsilon}.$

Then, we recover the limiting corrector
$Q(t,x,v,\hat{v})=\mathcal{Q}\big{(}x,\nabla_{x}\varphi(t,x),v,\hat{v}\big{)}$.

Furthermore, adding the condition (assumed to hold initially)

$\int_{\mathcal{V}}Q_{\varepsilon}(t,x,v,\hat{v})d\hat{v}dv=1,\qquad\forall
t\geq 0,\;x\in\Omega,$

which in turn implies

$\int_{\Omega}\exp^{-\dfrac{\tilde{\varphi}_{\varepsilon}(t,x)}{\varepsilon}}dx=1,\qquad\forall
t,$

we can expect that the problem (37) itself has a particular solution
$Q_{\varepsilon}(t,x,v,\hat{v})$ similar to the principal bundle, see [25,
26], for parabolic equations. It is similar to a time dependent eigenvalue
problem. Up to our knowledge this notion has never been studied for kinetic
equations. We can expect it defines a time-dependent functional Hamiltonian
$H_{\varepsilon}(t,x,[\nabla_{x}\tilde{\varphi}_{\varepsilon}])\to
H(x,\nabla_{x}\varphi(t,x))$. This allows us to search for the solution of a
functional Hamilton-Jacobi equation

$\partial_{t}\tilde{\varphi}_{\varepsilon}(t,x)+H_{\varepsilon}\big{(}t,x,[\nabla_{x}\tilde{\varphi}_{\varepsilon}]\big{)}=0.$
(41)

Then $f_{\varepsilon}(t,x,v,\hat{v})$ in (35) is an exact particular solution
of Eq. (13) solely determined by the initial concentration profile
$\varphi_{\varepsilon}(0,x)$. As $\varepsilon$ vanishes, this particular
solution attracts all solutions with the same initial concentration profile
$\varphi(0,x)$.

These formal conclusions rely on the possibility to define a smooth
Hamiltonian $H(x,p)$, a question we analyse now.

###  2.2 The effective Hamiltonian

As in [6], one can characterize the eigenvalue $H(x,p)$ arising in Eq. (38)
which can be written as

$\mathcal{Q}(x,p,v,\hat{v})=\dfrac{T[\mathcal{S}]_{0}(v,\hat{v})\int_{\mathcal{V}}\mathcal{Q}(x,p,w,\hat{w})d\hat{w}dw}{1+H(x,p)-v\hat{v}\cdot
p}>0.$ (42)

We remind that $H(x,p)$ also depends explicitly on $x$ as $T[\mathcal{S}]_{0}$
depends on $x$. Integrating with respect to $v,\hat{v}$ and using (39), we
obtain the following problem: find $H$ such that

$1=\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})}{1+H(x,p)-v\hat{v}\cdot
p}\,dvd\hat{v}.$ (43)

In particular we obviously have

$1=\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})}{1+H(x,0)}\,dvd\hat{v},\qquad
H(x,0)=0.$

Eq. (43) can be uniquely solved by strict decay in $H$ and also gives that

$-U|p|<H(x,p)<U|p|,$

because when $H=U|p|$ the denominator is larger than $1$ for all $(v,\hat{v})$
and when $H=-U|p|$ the denominator is smaller than $1$.

However these bounds are not enough to compute from this Hamiltonian a
positive eigenfunction $\mathcal{Q}$. As observed in [9], it is necessary to
introduce some further assumption. We define the values of $H$ when the
denominator vanishes as

$\underline{H}(p)=-1+\max_{\mathcal{V}}\,[v\hat{v}\cdot p].$

We then need to assume

$\inf_{x}\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})}{1+\underline{H}(p)-v\hat{v}\cdot
p}\,dvd\hat{v}>1.$ (44)

The latter ensures that $H(x,p)>\underline{H}(p)$ and thus that the
denominator, and therefore $\mathcal{Q}$, are positive for the solution of Eq.
(43). This integral blows-up in 1D and the condition is always satisfied. In
higher dimension this restriction is needed. This is the so-called
‘dimensionality problem’ as mentioned in [9].

Additionally, differentiating (43) in $p$, we find (ignoring the dependence on
$x$ for simplicity)

$0=\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})(\nabla_{p}H(p)-v\hat{v})}{(1+H(p)-v\hat{v}\cdot
p)^{2}}\,dvd\hat{v}.$

Then, as we are interested in the minima points of $\varphi$ (that are the
maxima of $\rho$), we look for the values of the Hamiltonian $H$ in
$p=\nabla\varphi=0$, i.e., using the definition (5), we get

$\nabla_{p}H(x,0)={\bf U}_{\mathcal{S}}^{0}(x),$ (45)

which, in general, does not vanish, as already argued, as $T[\mathcal{S}]_{0}$
is not in principle symmetric as a function of $\hat{v}$, as instead assumed
in [6]. Moreover, differentiating twice, we find

$\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})D^{2}H(p)}{(1+H(p)-v\hat{v}\cdot
p)^{2}}\,dvd\hat{v}=2\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})(\nabla_{p}H(p)-v\hat{v})\otimes(\nabla_{p}H(p)-v\hat{v})}{(1+H(p)-v\hat{v}\cdot
p)^{3}}\,dvd\hat{v},$

and then we have that $D^{2}_{p}H$ is positive definite as (44) holds. We may
also compute

$D^{2}_{p}H(x,0)=2\mathbb{D}_{\mathcal{S}}(x).$ (46)

###  2.3 The concentration result

To simplify, we work in the full space, $\mathbb{R}^{d}\times\mathcal{V}$
instead of $\Omega\times\mathcal{V}$. We assume initially that uniformly in
$\varepsilon$

$\begin{cases}\varphi_{\varepsilon}^{0}\in
L_{\infty},\qquad|\nabla_{x}\varphi_{\varepsilon}^{0}|\in
L_{\infty},\qquad\varphi_{\varepsilon}^{0}(x,v,\hat{v})=\tilde{\varphi}^{0}(x)+O(\varepsilon),\\\
\partial_{t}\varphi_{\varepsilon}^{0}:=-{\bf
v}\cdot\nabla_{x}\varphi_{\varepsilon}^{0}+1-T[S]_{\varepsilon}(v,\hat{v}){\displaystyle\int_{\mathcal{V}}}\exp^{\dfrac{-\varphi_{\varepsilon}^{0}(x,w,\hat{w})+\varphi_{\varepsilon}^{0}(x,v,\hat{v})}{\varepsilon}}\,dwd\hat{w}\in
L_{\infty},\end{cases}$ (47)
$\int_{\mathbb{R}^{d}\times\mathcal{V}}f^{0}_{\varepsilon}(x,v,\hat{v})dxdvd\hat{v}=1,\qquad\int_{\mathbb{R}^{d}\times\mathcal{V}}|x|f^{0}_{\varepsilon}(x,v,\hat{v})dxdvd\hat{v}\quad\text{is
bounded},$ (48)

and for some constants $T_{-}>0$, $L_{M}>0$

$\begin{cases}T[\mathcal{S}]_{\varepsilon}\geq T_{-}>0,\qquad
T[\mathcal{S}]_{\varepsilon}+|\nabla_{x}T[\mathcal{S}]_{\varepsilon}|+|\nabla_{\bf
v}T[\mathcal{S}]_{\varepsilon}|\leq L_{M},\\\\[5.0pt]
T[\mathcal{S}]_{\varepsilon}\rightarrow
T[\mathcal{S}]_{0}\quad\text{uniformly}.\end{cases}$ (49)

Then we can prove the

######  Theorem 1.

We make the assumptions(44) and (47)–(49). Then, after extractions,
(i) $\varphi_{\varepsilon}$ is uniformly (in $\varepsilon$) bounded and
Lipschitz (locally in time),
(ii) $\varphi_{\varepsilon}$ converges locally uniformly on
$\mathbb{R}_{+}\times\mathbb{R}^{d}\times\mathcal{V}$ toward $\varphi$ where
$\varphi$ does not depend on $v,\hat{v}$. Moreover, $\varphi$ is the viscosity
solution of the Hamilton-Jacobi Eq. (40) with initial condition
$\varphi^{0}(x)$ and with a convex Hamiltonian $H(x,p)$ uniquely implicitly
determined by the formula (43),
(iii) $f_{\varepsilon}$ converges weakly to a measure $f$ supported by
$\\{\varphi=0\\}$.

Remarks. 1\. Compared to [6], the kernel $T[\mathcal{S}]_{\varepsilon}$
depends on $x$, which is an additional major technical difficulty. Also a
difference here is the $\varepsilon$ dependency which is not relevant with our
assumptions.
2\. The author in [9] faces the difficulty of gradient estimates as here. He
argues by limsup-liminf arguments which optimizes the assumptions. Here we do
not go to this elaborate method and use simpler arguments based on Lipschitz
estimates.
3\. When $\Omega$ is bounded then we impose no-flux boundary conditions [36]

$\int_{\mathcal{V}}f_{\varepsilon}v\hat{v}\cdot\boldsymbol{n}\,dvd\hat{v}=0,\qquad
x\in\partial\Omega.$

Therefore, considering (35), we have

$\exp^{-\dfrac{\tilde{\varphi}_{\varepsilon}}{\varepsilon}}\int_{\mathcal{V}}Q_{\varepsilon}(t,x,v,\hat{v})v\hat{v}\cdot\boldsymbol{n}\,dvd\hat{v}=0,\qquad
x\in\partial\Omega.$

and thus

$\int_{\mathcal{V}}Q_{\varepsilon}(t,x,v,\hat{v})v\hat{v}\cdot\boldsymbol{n}\,dvd\hat{v}=0,\qquad
x\in\partial\Omega,$

and then in the limit $\varepsilon\rightarrow 0^{+}$

$\int_{\mathcal{V}}\mathcal{Q}(x,p,v,\hat{v})v\hat{v}\cdot\boldsymbol{n}\,dvd\hat{v}=0,\qquad
x\in\partial\Omega.$ (50)

It would be interesting to investigate if this relation can be interpreted as
a Neumann boundary condition on $p=\nabla_{x}\varphi$. Remark that for small
$p$, Taylor expanding $\mathcal{Q}$ defined in (42) as a function of $p$ and
plugging the expression into (50), we find

$[\varphi{\bf
U}_{\mathcal{S}}^{0}+\mathbb{D}_{\mathcal{S}}^{0}\nabla\varphi]\cdot\boldsymbol{n}=0.$
(51)

4\. The corrector satisfies, according to the Hopf-Cole transform (35),
approximately
$\mathcal{Q}_{\varepsilon}=\dfrac{f_{\varepsilon}}{\rho_{\varepsilon}}$.
However, we also have that
$\dfrac{f_{\varepsilon}}{\rho_{\varepsilon}}\rightarrow T[\mathcal{S}]_{0}$
and $\mathcal{Q}_{\varepsilon}\rightarrow\mathcal{Q}$ when $\varepsilon$ tends
to zero. Actually, thanks to (38), we find
$\mathcal{Q}(x,0,v,\hat{v})=T[\mathcal{S}]_{0}(v,\hat{v})$.

###### Proof.

The proof uses standard arguments, see for instance [17, 2], and we only
sketch it. We begin with standard a priori estimates (i) for the solution of
Eq. (33). From assumption (47), we infer

$|\partial_{t}\varphi_{\varepsilon}(t,x,v,\hat{v})|\leq|\partial_{t}\varphi_{\varepsilon}^{0}(x,v,\hat{v})|\leq
C,$

therefore

$|\varphi_{\varepsilon}(t,x,v,\hat{v})|\leq|\varphi_{\varepsilon}^{0}(x,v,\hat{v})|+Ct.$

Also, still using the maximum principle for derivatives and the already proved
bounds (here Eq. (33) is used again for $x$-derivative of $T[S]$), we have
that

$\sum_{i=1}^{d}|\partial_{i}\varphi_{\varepsilon}(t,x,v,\hat{v})|\leq\sum_{i=1}^{d}|\partial_{i}\varphi^{0}_{\varepsilon}(x,v,\hat{v})|e^{Ct}+Ct.$

Since this estimate is more elaborate, we prove it. Differentiating Eq. (33)
in $x_{i}$, and setting
$\psi_{i}=\partial_{i}\varphi_{\varepsilon}(t,x,v,\hat{v})$, we find

$\partial_{t}\psi_{i}+{\bf
v}\cdot\nabla\psi_{i}=-T[S]_{\varepsilon}\int_{\mathcal{V}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})+\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}\frac{\psi_{i}(v,\hat{v})-\psi_{i}(w,\hat{w})}{\varepsilon}\,dwd\hat{w}+RHS$

where the RHS term is

$RHS=\partial_{i}T[S]_{\varepsilon}(v,\hat{v})\int_{\mathcal{V}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})+\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}\,dwd\hat{w}.$

Using again Eq. (33), it can be estimated as

$\begin{array}[]{ll}|RHS|&\leq\frac{|\partial_{i}T[S]_{\varepsilon}(v,\hat{v})|}{T[S]_{\varepsilon}(v,\hat{v})}T[S]_{\varepsilon}(v,\hat{v})\int_{\mathcal{V}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})+\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}\,dwd\hat{w}\\\\[5.0pt]
&=\frac{|\partial_{i}T[S]_{\varepsilon}(v,\hat{v})|}{T[S]_{\varepsilon}(v,\hat{v})}[1-\partial_{t}\varphi_{\varepsilon}-{\bf
v}\cdot\nabla\varphi_{\varepsilon}],\end{array}$

and, thus, using the time derivative estimate, we conclude that

$|RHS|\leq Ct+C\sum_{i=1}^{d}|\psi_{i}|.$

With this observation, we can use the maximum principle for
$\psi_{i}e^{-Ct}+C$ and conclude the bounds on the $x$ derivatives.

With these estimates, we conclude that for $t\leq T$, we have

$\int_{\mathcal{V}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})+\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}\,dwd\hat{w}\leq
C(T),$

which tells us that a limit of $\varphi_{\varepsilon}$ depends only on
$(t,x)$. As in [6], it also gives directly the last estimate of (i), that is

$|\partial_{{\bf
v}}\varphi_{\varepsilon}(t,x,v,\hat{v})|\leq\big{[}|\partial_{x}\varphi^{0}_{\varepsilon}(x,v,\hat{v})|+|\partial_{{\bf
v}}\varphi^{0}_{\varepsilon}(x,v,\hat{v})|\big{]}e^{Ct}+Ct.$

We are now in the same situation as [6] and the rest of the argument follows
in a similar way. Using the perturbed test function method, [18], we obtain
the statement (ii) thanks to the assumption (44) which allows us to handle
$\mathcal{Q}.$

Finally, for the statement (iii), we notice that the mass conservation is
immediate. Then, we observe that

$\frac{d}{dt}\int_{\mathbb{R}^{d}\times\mathcal{V}}|x|f_{\varepsilon}(t,x,v,\hat{v})dxdvd\hat{v}\leq\int_{\mathbb{R}^{d}\times\mathcal{V}}|{\bf
v}|f_{\varepsilon}(t,x,v,\hat{v})dxdvd\hat{v}\leq U.$

Therefore $f_{\varepsilon}$ is a tight probability measure and, after
extraction, it converges weakly to a probability measure and the only possible
concentration points are when $\varphi(t,x)$ is zero (see [3, 29] for details
and consequences). ∎

###  2.4 An eikonal equation

Given (21), in the limit $\varepsilon\rightarrow 0^{+}$, the phase
$\phi_{\varepsilon}=-\varepsilon\log\rho_{\varepsilon}$ satisfies the
Hamilton-Jacobi equation

$\partial_{t}\phi+{\bf
U}_{\mathcal{S}}^{0}\cdot\nabla\phi+\textrm{tr}\left[\mathbb{D}_{\mathcal{S}}^{0}(\nabla_{x}\phi\otimes\nabla_{x}\phi)\right]=0.$
(52)

When working in a bounded domain, from (22), we additionally obtain the
boundary condition

$[\phi{\bf
U}_{\mathcal{S}}^{0}+\mathbb{D}_{\mathcal{S}}^{0}\nabla\phi]\cdot\boldsymbol{n}=0,$

that is the same as (51). As $\varphi$ satisfies (40) and $\phi$ satisfies
(52), we should be lead to conclude, as observed in [6], that the two
procedures (aggregate quantities limit and WKB analysis) do not commute in
general, in particular because the Hilbert expansion is additive, while the
Hopf-Cole one is multiplicative. However, Eq. (52) may be seen as (40) where
the quadratic expansion of the Hamiltonian $H$ in a neighborhood of
$\nabla\phi=0$, that characterizes the minima points, is considered,
remembering (45)-(46). By exploiting this observation, it is possible to
detect a regime in which the two procedures may commute.

Let us consider the regime (27) or (29) and the assumption of small
diffusivity (30). Then, considering $\rho_{\nu}=\exp^{-\dfrac{\phi}{\nu}}$ and
letting $\nu\rightarrow 0$, Eq. (27) becomes

$\partial_{t}\phi+{\bf
U}_{\mathcal{S}}\cdot\nabla\phi+\textrm{tr}\left[\bar{\mathbb{D}}_{\mathcal{S}}(\nabla_{x}\phi\otimes\nabla_{x}\phi)\right]=0,$
(53)

while Eq. (29) becomes

$\partial_{t}\phi++{\bf
U}_{\mathcal{S}}^{1}\cdot\nabla\phi+\textrm{tr}\left[\bar{\mathbb{D}}_{\mathcal{S}}^{0}(\nabla_{x}\phi\otimes\nabla_{x}\phi)\right]=0.$
(54)

We now consider

$\varepsilon\rightarrow\varepsilon\nu$

and plug (2) in (25); in the limit $\varepsilon\rightarrow 0$, we obtain

$\partial_{t}(\left[\dfrac{\varphi}{\nu}\right])+H_{\nu}(\nabla\left[\dfrac{\varphi}{\nu}\right])=0,$
(55)

where $H_{\nu}$ is implicitly defined by

$1=\nu\int_{\mathcal{V}}\dfrac{T[S]_{0}(v,\hat{v})}{\nu+H(x,p)-v\hat{v}\cdot
p}\,dvd\hat{v},\quad p=\nabla\left[\dfrac{\varphi}{\nu}\right].$

Then, as $D^{2}H(0)=\dfrac{2}{\nu^{2}}\mathbb{D}_{\mathcal{S}}^{0}$,
considering $\nu$ small and assuming the small diffusivity (30), we obtain

$H(x,\nabla\varphi)=H(x,0)+\nabla\varphi\cdot\nabla_{p}H(x,0)+\dfrac{1}{2}\textrm{tr}\left[\nabla\varphi\otimes\nabla\varphi
D_{p}^{2}H(x,0)\right]={\bf
U}_{\mathcal{S}}^{0}\cdot\nabla\varphi+\textrm{tr}\left[\bar{\mathbb{D}}_{\mathcal{S}}^{0}\nabla\varphi\otimes\nabla\varphi\right].$

Plugging the latter in (55) allows to obtain the equivalent form of Eq. (52).

Let us now consider a spatially homogeneous $T$ such that

$\int_{\mathcal{V}}Tdvd\hat{v}=1,\qquad\int_{\mathcal{V}}Tvdvd\hat{v}=0,\qquad\int_{\mathcal{V}}Tv^{2}\hat{v}\otimes\hat{v}\,dvd\hat{v}^{2}=\nu^{2}\mathbb{I}.$

It is the case for example, in one-dimension, where $\mathcal{V}=[-U,U]$, and
we choose $T=\exp^{-\dfrac{v^{2}}{2}}$, $T=\dfrac{1}{U}$, or
$T=\dfrac{1}{2}\left(\delta(v-U)+\delta(v+U)\right)$, see [6]. Then,
performing the WKB analysis leads to (55) and, then, to
$\partial_{t}\varphi+|\nabla^{2}\varphi|=0$, while, starting from (32) we
obtain $\partial_{t}\phi+|\nabla^{2}\phi|=0$, that is (54) with ${\bf
U}_{\mathcal{S}}^{1}=0$. We highlight, in fact, that in these cases there is
no correction term $T[\mathcal{S}]_{1}$ (and, then no ${\bf
U}_{\mathcal{S}}^{1}$). In conclusion, the two procedures commute in the
regime of small $\nu$.

Interestingly, in the case

$T[\mathcal{S}]=c(x)\mathcal{S}(x+R\hat{v}),$

considering (20) or (26) and a large $R$, the two procedures lead to (27) and
(53), respectively. Conversely, if we consider a small $R$ and
$T[\mathcal{S}]_{0},T[\mathcal{S}]_{1}$ as defined as a consequence of a
localized scaling (28), we obtain (29) with
$\mathbb{D}_{\mathcal{S}}^{0}=\mathbb{I}$, ${\bf
U}_{\mathcal{S}}^{1}=R\dfrac{\nabla\mathcal{S}}{\mathcal{S}}$. If in the WKB
analysis we consider (14), we obtain the eikonal equation
$\partial_{t}\varphi+|\nabla^{2}\varphi|=0$ and the higher order effect is
naturally lost in the localized hyperbolic scaling.

###  2.5 Dynamics of the concentration points

Not only the Hamilton-Jacobi and eikonal equations give the microscopic shape
of the Dirac concentration of solutions, when they occur, but it also allows
to recover their dynamics. For this we look for the trajectory $\bar{x}(t)$ of
the maxima of $\rho$. In the context of adaptive dynamics, this is interpreted
as the ‘fittest trait’, [3, 29, 26]. First of all, we remark that in the
regime (13) for $\varepsilon\rightarrow 0$, each trajectory in the physical
space (see (23)) follows the differential equation

$\dot{x}(t)={\bf U}_{\mathcal{S}}^{0}(x(t)).$ (56)

For example, in the case $b(\mathcal{S})=\mathcal{S}$ and specifically in the
regime (14), we have that ${\bf U}_{\mathcal{S}}^{0}(x(t))=0\,\forall x$,
while in the regime (20) in general ${\bf U}_{\mathcal{S}}^{0}(x(t))$ does not
vanish and $\dot{x}_{i}(t)=0$ if $\exists x_{i}(t)\,:\,{\bf
U}_{\mathcal{S}}(x_{i}(t))=0$. This does not tell us the position of the
maxima, but as it is satisfied by all points, then it will be also satisfied
by the maxima. In order to find the trajectory of the maxima and their
position, we should look for the points $\bar{x}_{i}$ s.t.

$\varphi(t,\bar{x}_{i})=0$ (57)

that are the points of minimum of $\varphi$, i.e.

$\nabla_{x}\varphi(t,\bar{x}_{i})=0$ (58)

and

$D^{2}_{x}\varphi(t,\bar{x}_{i})\geq 0.$

Then, differentiating (57) with respect to time along the trajectories we find

$0=\dfrac{d}{dt}\varphi(t,\bar{x}_{i})=\partial_{t}\varphi(t,\bar{x}_{i})+\dot{\bar{x}}_{i}\cdot\nabla_{x}\varphi(t,\bar{x}_{i}).$

Because of (58), then we have that

$\partial_{t}\varphi(t,\bar{x}_{i})=0$

and, therefore, from (40) the minima satisfy $H(\bar{x}_{i},0)=0$. Then we
compute

$0=\dfrac{d}{dt}\nabla_{x}\varphi(t,\bar{x}_{i})=\partial_{t}\nabla_{x}\varphi(t,\bar{x}_{i})+\dot{\bar{x}}_{i}D^{2}_{x}\varphi(t,\bar{x}_{i}).$

Differentiating (40) with respect to $x$, we get

$\partial_{t}\nabla_{x}\varphi(t,\bar{x}_{i})=-\nabla_{x}H(x,\nabla_{x}\varphi(t,x))-\nabla_{p}H(x,\nabla_{x}\varphi(t,x))D^{2}_{x}\varphi(t,x),$

and, specializing it in $\bar{x}_{i}$, and remembering (45) and $H(x,0)=0$, we
get

$\partial_{t}\nabla_{x}\varphi(t,\bar{x}_{i})=-\nabla_{p}H(x,\nabla_{x}\varphi(t,\bar{x}_{i}))\cdot
D^{2}_{x}\varphi(t,\bar{x}_{i})=-{\bf
U}_{\mathcal{S}}^{0}(\bar{x}_{i}(t))\cdot D^{2}_{x}\varphi(t,\bar{x}_{i}).$

Therefore we obtain

$\dot{\bar{x}}_{i}(t)={\bf U}_{\mathcal{S}}^{0}(\bar{x}_{i}(t)),$

that is the same as (56). In particular the long term limit is solely
determined by ${\bf U}_{\mathcal{S}}^{0}$ and does not depend on the initial
condition of $\rho$ as it is usual in adaptive dynamics and as it follows from
the H-Theorem in kinetic theory, that establishes in this linear case that the
equilibrium is asymptotically stable and does not depend on the initial
condition.From the eikonal equation we may expect that

$\rho_{0}=\sum\dfrac{\rho^{\infty}}{\textrm{tr}\mathbb{D}_{\mathcal{S}}^{0}(\bar{x}_{i})}\delta(x-\bar{x}_{i}),\qquad\dot{\bar{x}}_{i}={\bf
U}_{\mathcal{S}}(\bar{x}_{i})=0.$

###  2.6 Examples

We illustrate the results with two examples in one dimension. We first choose
the signal

$\mathcal{S}(x)=\mathcal{S}_{0}\exp^{-\dfrac{(x-\bar{x})^{2}}{2\sigma^{2}}},$
(59)

with $\bar{x}$ a given point in $\Omega$. When $b(\mathcal{S})=\mathcal{S}$,
it generates a transition probability given by

$T[\mathcal{S}](v,\hat{v})=\psi(v|\hat{v})\dfrac{\exp^{-\dfrac{R(x-\bar{x})\cdot\hat{v}}{2\sigma^{2}}}}{\int_{\mathbb{S}^{d-1}}\exp^{-\dfrac{R(x-\bar{x})\cdot\hat{v}}{2\sigma^{2}}}\,d\hat{v}}.$

Firstly, we remark that

$l_{\mathcal{S}}=\dfrac{\sigma^{2}}{\max|x-\bar{x}|}.$

In the regime (20), we expect that the singular point ${\bf
U}_{\mathcal{S}}^{0}(x(t))=0$ only occurs when $x(t)=\bar{x}$. Moreover, we
expect a unique nonhomogeneous stationary state whose profile needs to satisfy
(12). In the regime (14), conversely, as the limiting $T[\mathcal{S}]_{0}$
does not depend on $x$, we expect that the homogeneous configuration is the
stationary equilibrium, as $\partial_{t}\rho=0$. We now consider the Hopf-Cole
analysis. In 1D we have that $\hat{v}=\pm 1$ and choosing
$\psi(v|\hat{v})=\delta(v-\hat{v}V)$ with $V,R=\mathcal{O}(1)$ we get

$H(x,p)=\dfrac{V^{2}p^{2}+Vp\mathcal{D}_{\mathcal{S}}(x)}{1+\sqrt{1+4V^{2}p^{2}+4Vp\mathcal{D}_{\mathcal{S}}(x)}},\qquad\mathcal{D}_{\mathcal{S}}(x)=\dfrac{\exp^{\dfrac{R(x-\bar{x})}{2\sigma^{2}}}-\exp^{-\dfrac{R(x-\bar{x})}{2\sigma^{2}}}}{\exp^{\dfrac{R(x-\bar{x})}{2\sigma^{2}}}+\exp^{-\dfrac{R(x-\bar{x})}{2\sigma^{2}}}}.$
(60)

We remark that, therefore, in the regime (20), the Hamiltonian has two
different zeros $p(x)=0,-\dfrac{\mathcal{D}_{\mathcal{S}}(x)}{2V}$.
Conversely, in the regime (14), in (2) (and, then, in (60)) $R$ is to be
replaced with $\varepsilon R$ and, hence, the hamiltonian vanishes only in
$p=0$ and every concentration disappears as $\nabla_{p}H(x,0)\equiv 0$.

We numerically solve equation (1). We consider $\Omega=[0,1]$, $\mu=1$ and
$V=5\cdot 10^{-5}$ and $\psi(v)=\dfrac{1}{V}$. In all simulations the space
grid has a uniform mesh defined by $dx=10^{-3}$. In Fig. 1 we use
$\mathcal{S}$ given by (59) with $\sigma=0.05$ and $\bar{x}=1$, $R=0.01$.
Therefore $l_{\mathcal{S}}=5\cdot 10^{-3}$. Let us consider $L=1$. We are then
in regime (20) with $\varepsilon=10^{-5}$. We consider two different initial
conditions $\rho^{0}=0.1$ and $\rho^{0}$ Gaussian centered in $1.5$. As
already mentioned the stationary state is unique and does not depend on the
initial condition. In the second line of Fig 1, second and third panel, we
plot the Hamiltonian (60). We remark that the Hamiltonian is not always
positive and there is a concentration profile.

Figure 1: First line. First panel: $\mathcal{S}$ (red curve) given by (59) and
two different initial conditions $\rho^{0}$: the constant one (blue) and an
asymmetric Gaussian (green). Second panel: temporal evolution of $\rho(t,x)$
in the case of constant $\rho^{0}$. Third panel: temporal evolution of
$\rho(t,x)$ in the case of an asymmetric $\rho^{0}$. Second line. First panel:
$H(x,p)$ in the local regime (14), second and third panel: $H(x,p)$ in the
nonlocal regime (20). In the third panel the red dashed lines correspond to
$\pm U|p|$, while the black horizontal dashed line indicates the level zero.

Conversely, in a regime defined by $V=5\cdot 10^{-5}$, $L=l_{\mathcal{S}}$
(i.e. regime (14) with $\varepsilon=10^{-2}$) the stationary state
$\rho_{\infty}$ is the stationary homogeneous configuration even for a
nonhomogeneous initial condition (not shown), and this is true in both regimes
$R\lessgtr l_{\mathcal{S}}$. The corresponding Hamiltonian is plotted in
Fig.1, second line, first panel.

Second, we choose a bimodal signal

$\mathcal{S}(x)=\mathcal{S}_{1}\exp^{-\dfrac{(x-\bar{x}_{1})^{2}}{2\sigma_{1}^{2}}}+\mathcal{S}_{2}\exp^{-\dfrac{(x-\bar{x}_{2})^{2}}{2\sigma_{2}^{2}}}.$
(61)

Then the number of singular points $x_{i}$ satisfying ${\bf
U}_{\mathcal{S}}(x_{i})=0$ depends on $R$, $\bar{x}_{1}-\bar{x}_{2}$ and on
$\sigma^{2}$.

In Fig. 2 we consider the bimodal signal $\mathcal{S}$ given by (61) with
$\sigma_{1}=\sigma_{2}=0.03$ and three different couples of
$\bar{x}_{1},\bar{x}_{2}$ according to the value of their distance with
respect to $R=0.4$. Here again $\varepsilon=10^{-5}$ (same values of
$V,L,\mu$). We remark that, when $\bar{x}_{2}-\bar{x}_{1}\leq R$ (see Fig.
2(b) for the case $\bar{x}_{2}-\bar{x}_{1}=R$), then there is a single peak as
$\exists!$ $x_{1}(=0)$ such that ${\bf U}_{\mathcal{S}}(x_{1})=0$. The case
$\bar{x}_{2}-\bar{x}_{1}<R$ behaves like $\bar{x}_{2}-\bar{x}_{1}=R$ (not
shown). When $\bar{x}_{2}-\bar{x}_{1}>R$ (see Fig. 2(c)) then $\exists$
$\bar{x}_{1},\bar{x}_{2}$ such that ${\bf
U}_{\mathcal{S}}(\bar{x}_{i})=0,\,i=1,2$. In Fig. 2(d) we plot the Hamiltonian
in the case $\bar{x}_{2}-\bar{x}_{1}>R$.

(a) $\mathcal{S}$ (b) $\bar{x}_{2}-\bar{x}_{1}=R$

(c) $\bar{x}_{2}-\bar{x}_{1}>R$ (d) $\bar{x}_{2}-\bar{x}_{1}>R$

Figure 2: Temporal evolution of $\rho(t,x)$ in case of $\mathcal{S}$ given by
(61). Here $R=0.4$, the initial condition $\rho^{0}$ is constant. In (a):
three different $\mathcal{S}$ as given by (61) with
$\sigma_{1}=\sigma_{2}=0.03$ and three different couples of
$\bar{x}_{2},\bar{x}_{1}$. In (b)-(c) temporal evolution of $\rho(t,x)$ for
the two different $\mathcal{S}$: in (b) $\bar{x}_{2}-\bar{x}_{1}=R$, in (c)
$\bar{x}_{2}-\bar{x}_{1}>R$. In (d) we plot the Hamiltonian for the case
$\bar{x}_{2}-\bar{x}_{1}>R$.

In Fig. 3 we consider $\mathcal{S}$ given by (61) with
$\sigma_{1}=2\sigma_{2},\sigma_{2}=0.03$ while for the sensing radius we have
again $R=0.4$. We consider two different couples of $\bar{x}_{1},\bar{x}_{2}$
as shown in Fig. 3(a). We remark that the peaks of $\rho$ do not coincide with
the peaks of $\mathcal{S}$ and this is due to the nonlocality ($R>0$). In
particular, in the first case (b) the distance between the maxima of
$\mathcal{S}$ is larger than the sensing radius and, thus, the stationary
solution has two peaks even though with different convexity due to the
configuration of $\mathcal{S}$. In the second case (c), the distance between
the peaks of $\mathcal{S}$ is exactly $R/2$, so that the peak of the
stationary solution is unique, even though it is asymmetric, because of the
asymmetry of $\mathcal{S}$.

(a) $\mathcal{S}$ (b) $\bar{x}_{2}-\bar{x}_{1}=2R$ (c)
$\bar{x}_{2}-\bar{x}_{1}=R/2$

Figure 3: Temporal evolution of $\rho(t,x)$ in case of $\mathcal{S}$ given by
(61). Here $R=0.4$ and $\sigma_{1}=0.06,\sigma_{2}=0.03$. In (a): two
different $\mathcal{S}$ given by (61) with two different couples
$\bar{x}_{2},\bar{x}_{1}$. In (b) and (c) the corresponding time evolution of
the densities $\rho$.

In conclusion, in this linear case the analysis of the kinetic equations and
of the aggregate limits give almost a complete set of information concerning
the dynamics of the maxima, except the concentration result, for which the WKB
analysis is needed. Therefore, we now consider a nonlinear case in which the
study of the kinetic and aggregate equations may not be able to convey all the
necessary information regarding the dynamics of the maxima points, while the
Hamilton-Jacobi formalism offers promising tools in order to describe the
concentration profiles.

##  3 A nonlinear case

When the external field affecting the choice of the reorientation is the cell
density itself, i.e., $\mathcal{S}=\rho$, then Eq. (1) becomes nonlinear. The
kinetic equation is

$\partial_{t}f(t,x,v,\hat{v})+{\bf
v}\cdot\nabla_{x}f(t,x,v,\hat{v})=\mu\left(\rho(t,x)T[\rho](v,\hat{v})-f(t,x,v,\hat{v})\right),$
(62)

where

$T[\rho]=\dfrac{\rho(t,x+R\hat{v})}{\int_{\mathbb{S}^{d-1}}\rho(t,x+R\hat{v})\,d\hat{v}}\psi(v|\hat{v}).$
(63)

In [33] the authors perform a linear stability analysis. In 1D and with the
choice $\psi(v|\hat{v})=\delta(v-V_{\psi}(\hat{v}))$, they show that the
uniform homogeneous configuration is stable if, using the notation (4),

$\dfrac{V}{R\mu}>1,\qquad V=\dfrac{V_{\psi}(+1)+V_{\psi}(-1)}{2}.$ (64)

We now consider regime (14), and in this nondimensionalized regime
$V,R,\mu=\mathcal{O}(1)$, then Eq. (62) reads

$\dfrac{\partial f_{\varepsilon}}{\partial t}(t,x,v,\hat{v})+{\bf
v}\cdot\nabla
f_{\varepsilon}(t,x,v,\hat{v})=\dfrac{\mu}{\varepsilon}\,\Big{(}\rho_{\varepsilon}T[\rho]_{\varepsilon}-f_{\varepsilon}(t,x,v,\hat{v})\Big{)},$
(65)

where

$T[\rho]_{\varepsilon}=\dfrac{\rho_{\varepsilon}(t,x+\varepsilon
R\hat{v})}{\int_{\mathbb{S}^{d-1}}\rho_{\varepsilon}(t,x+\varepsilon
R\hat{v})\,d\hat{v}}\psi(v|\hat{v}).$

In particular, in the rescaled regime (14), relation (64) is unchanged since

$\dfrac{\varepsilon V}{R\varepsilon\mu}=\dfrac{V}{R\mu}.$

At the macroscopic level we have that

$\partial_{t}\rho_{\varepsilon}+\nabla_{x}\cdot({\bf
U}_{\rho}^{\varepsilon}\rho_{\varepsilon})=\varepsilon\nabla_{x}\cdot\nabla_{x}\cdot\left(\mathbb{D}_{\rho}^{\varepsilon}\rho\right)$
(66)

with ${\bf U}_{\rho}^{\varepsilon}$ and $\mathbb{D}_{\rho}^{\varepsilon}$
defined thanks to (5), (6).

When $\rho_{\varepsilon}$ is smooth enough, the limiting transition
probability becomes

$\lim_{\varepsilon\rightarrow
0}T[\rho]_{\varepsilon}=T[\rho]_{0}:=\dfrac{1}{|\mathbb{S}^{d-1}|}\psi(v|\hat{v}).$

###  3.1 Concentration profile and the Hamilton–Jacobi equation

However, when $\rho_{\varepsilon}$ undergoes concentration, we may consider
(2) and we have

$\displaystyle\begin{aligned}
&T[\rho]_{\varepsilon}=\psi(v|\hat{v})\dfrac{\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\exp^{-\dfrac{\varphi_{\varepsilon}(t,x+\varepsilon
R\hat{v},w,\hat{w})}{\varepsilon}}dwd\hat{w}}{\int_{\mathbb{S}^{d-1}}\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\exp^{-\dfrac{\varphi_{\varepsilon}(t,x+\varepsilon
R\hat{v},w,\hat{w})}{\varepsilon}}\,dwd\hat{w}d\hat{v}}\\\
&\approx\psi(v|\hat{v})\dfrac{\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\exp^{\dfrac{-\varphi_{\varepsilon}(t,x,w,\hat{w})-\varepsilon
R\hat{v}\cdot\nabla\varphi_{\varepsilon}(t,x,w,\hat{w})}{\varepsilon}}\,dwd\hat{w}}{\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\int_{\mathbb{S}^{d-1}}\exp^{-\dfrac{\varphi_{\varepsilon}(t,x,w,\hat{w})+\varepsilon
R\hat{v}\cdot\nabla_{x}\varphi}{\varepsilon}}dwd\hat{w}\,d\hat{v}}\end{aligned}$

where we have used

$\varphi_{\varepsilon}(t,x+\varepsilon
R\hat{v},w,\hat{w})=\varphi_{\varepsilon}(t,x,w,\hat{w})+\varepsilon
R\hat{v}\cdot\nabla_{x}\varphi_{\varepsilon}(t,x,w,\hat{w}).$

Then, assuming (35) and, remembering that $\tilde{\varphi}_{\varepsilon}$ does
not depend on $v,\hat{v}$, therefore

$T[\rho]_{\varepsilon}\approx\psi(v|\hat{v})\dfrac{\exp^{\dfrac{-\tilde{\varphi}_{\varepsilon}(t,x)}{\varepsilon}}\exp^{-R\hat{v}\cdot\nabla\tilde{\varphi}_{\varepsilon}(t,x)}\int_{0}^{U}\int_{\mathbb{S}^{d-1}}Q_{\varepsilon}\,dwd\hat{w}}{\exp^{\dfrac{-\tilde{\varphi}_{\varepsilon}(t,x)}{\varepsilon}}\int_{\mathbb{S}^{d-1}}\int_{0}^{U}\int_{\mathbb{S}^{d-1}}Q_{\varepsilon}dwd\hat{w}\exp^{-R\hat{v}\cdot\nabla\tilde{\varphi}_{\varepsilon}}\,d\hat{v}},$

where we have used the Fubini-Tonelli theorem. Then, letting
$\varepsilon\rightarrow 0^{+}$, we obtain

$T[\rho]_{\varepsilon}\rightarrow
G_{R}(v,\hat{v},\nabla_{x}\varphi)=\psi(v|\hat{v})\dfrac{\exp^{-R\hat{v}\cdot\nabla_{x}\varphi}}{\displaystyle\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot\nabla_{x}\varphi}\,d\hat{v}}.$
(67)

Therefore, plugging (2) in (65) we obtain

$\displaystyle\mu-\partial_{t}\varphi_{\varepsilon}-{\bf
v}\cdot\nabla_{x}\varphi_{\varepsilon}=$
$\displaystyle\exp^{\dfrac{\varphi_{\varepsilon}(t,x,v,\hat{v})}{\varepsilon}}$
(68)
$\displaystyle\dfrac{\mu}{\varepsilon}\left[\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\exp^{-\dfrac{\varphi_{\varepsilon}(t,x,w,\hat{w})}{\varepsilon}}\\!\\!\\!\\!dwd\hat{w}\;T[\rho]_{\varepsilon}\right].$

Furthermore, considering the expansion (35) and by letting (formally)
$\varepsilon$ go to $0^{+}$, we obtain

$\displaystyle\mu-\partial_{t}\varphi$ $\displaystyle-{\bf
v}\cdot\nabla_{x}\varphi$ (69)
$\displaystyle=\mu\left[Q^{-1}(x,v,\hat{v})\int_{0}^{U}\int_{\mathbb{S}^{d-1}}Q(x,w,\hat{w})\,dwd\hat{w}G_{R}(v,\hat{v},\nabla_{x}\varphi)\right].$

Again, we assume that $Q(x,w,\hat{w})$ is positive according to assumption
(44). Then, like in the previous section, we define an eigenvalue-eigenvector
problem

$\displaystyle(\mu+H(p)-v\hat{v}\cdot p){\mathcal{Q}}(p,v,\hat{v})=\mu
G_{R}(v,\hat{v},p)\int_{0}^{U}\int_{\mathbb{S}^{d-1}}{\mathcal{Q}}(p,w,\hat{w})\,dwd\hat{w}.$

We remark that the term $G_{R}$ results from the interaction kernel and arises
due to the nonlocal sensing of $\rho$. As such, its role is to drive the
dynamics of $f$ toward the equilibrium in $v,\hat{v}$. We remark that it
satisfies:

$\nabla_{p}G_{R}(v,\hat{v},p)=\psi(v|\hat{v})\left[\dfrac{-R\hat{v}\exp^{-R\hat{v}\cdot
p}}{\displaystyle\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{w}\cdot
p}\,d\hat{w}}+\dfrac{\exp^{-R\hat{v}\cdot
p}R\displaystyle\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{w}\cdot
p}\hat{w}\,d\hat{w}}{\left(\displaystyle\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{w}\cdot
p}\,d\hat{w}\right)^{2}}\right],$

and, then

$G_{R}(v,\hat{v},0)=\dfrac{\psi(v|\hat{v})}{|\mathbb{S}^{d-1}|},\qquad\nabla_{p}G_{R}(v,\hat{v},0)=-\psi(v|\hat{v})\dfrac{R}{|\mathbb{S}^{d-1}|}\hat{v}.$

From (68), (69), we find that the (formal) limit $\varphi$ is the solution of

$\partial_{t}\varphi+H(\nabla_{x}\varphi)=0,$ (70)

where the Hamiltonian is implicitly defined by

$1=\mu\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\dfrac{G_{R}(v,\hat{v},p)}{\mu+H(p)-v\hat{v}\cdot
p}\,dvd\hat{v}.$ (71)

We remark that now $H$ only depends on $p$ and not on $x$. It is easy to see
that $H(0)=0$. Then, by differentiating (71) with respect to $p$ we obtain

$0=\mu\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\dfrac{\nabla_{p}G_{R}(v,\hat{v},p)}{(\mu+H(p)-v\hat{v}\cdot
p)}\,dvd\hat{v}-\mu\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\dfrac{G_{R}(v,\hat{v},p)(\nabla_{p}H(p)-v\hat{v})}{(\mu+H(p)-v\hat{v}\cdot
p)^{2}}\,dvd\hat{v}.$

Therefore $\nabla_{p}H(p)={\bf U}_{R}^{G}(p)$, where

${\bf
U}_{R}^{G}(p)=\dfrac{\int_{\mathbb{S}^{d-1}}V_{\psi}(\hat{v})\exp^{-R\hat{v}\cdot\nabla_{x}\varphi}\,\hat{v}d\hat{v}}{\int_{\mathbb{S}^{d-1}}V_{\psi}(\hat{v})\exp^{-R\hat{v}\cdot\nabla_{x}\varphi}\,d\hat{v}},$
(72)

so that

$\nabla_{p}H(0)=\int_{\mathbb{S}^{d-1}}V_{\psi}(\hat{v})\hat{v}\,d\hat{v}={\bf
U}_{R}^{G}(0).$ (73)

We remark that $\nabla_{p}H(0)$ vanishes in the case in which $V_{\psi}$ is
even.

Differentiating further, the Hessian of $H$ satisfies

$\displaystyle\int_{0}^{U}\int_{\mathbb{S}^{d-1}}$
$\displaystyle\dfrac{G_{R}(v,\hat{v},p)D^{2}H(p)}{(\mu+H-v\hat{v}p)^{2}}\,dvd\hat{v}=2\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\dfrac{G_{R}(\nabla_{p}H-v\hat{v})\otimes(\nabla_{p}H-v\hat{v})}{(\mu+H-v\hat{v}p)^{3}}\,dvd\hat{v}$
$\displaystyle+\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\dfrac{D_{p}^{2}G_{R}}{\mu+H-v\hat{v}p}\,dvd\hat{v}-2\int_{0}^{U}\int_{\mathbb{S}^{d-1}}\dfrac{\nabla_{p}G_{R}(\nabla_{p}H-v\hat{v})}{(\mu+H-v\hat{v}p)^{2}}\,dvd\hat{v}.$

We have $D_{p}^{2}G_{R}(v,\hat{v},0)=0$ because

$\displaystyle D_{p}^{2}G_{R}(v,\hat{v},p)=$
$\displaystyle\psi(v|\hat{v})R^{2}\exp^{-R\hat{v}\cdot
p}\left[\dfrac{\hat{v}\otimes\hat{v}}{\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\,d\hat{v}}-2\dfrac{\hat{v}\otimes\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\hat{v}\,d\hat{v}}{(\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\,d\hat{v})^{2}}\right.$
$\displaystyle\left.-\dfrac{\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\hat{v}\otimes\hat{v}\,d\hat{v}}{(\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\,d\hat{v})^{2}}+2\dfrac{\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\,d\hat{v}\otimes\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\,d\hat{v}}{(\int_{\mathbb{S}^{d-1}}\exp^{-R\hat{v}\cdot
p}\,d\hat{v})^{3}}\right].$

Then, in the case in which $\psi$ does not depend on $\hat{v}$, we have
$V=V_{\psi}$ and

$\displaystyle
D_{p}^{2}H(0)=\dfrac{2dD^{2}}{\mu}\mathbb{I}-2dVR\mathbb{I}=2d\big{(}\dfrac{D^{2}}{\mu}-VR\big{)}\;\mathbb{I}$
(74)

that is positive definite when

$\dfrac{D^{2}}{V}>R\mu.$ (75)

Notice that $D^{2}=V^{2}$ when $\psi$ is a Dirac delta and we recover, here in
any dimension $d$, the linear stability region determined in 1D in [33]. For
any other choice of $\psi$, we have that $D^{2}=V^{2}+e$, where $e\geq 0$ is
the variance of $\psi$. Then, when (64) is satisfied, the condition (75) with
$D^{2}=V^{2}+e$ is automatically satisfied. In conclusion, the choice of
$\psi$ being a Dirac delta is the most unstable one and, therefore, this
choice makes computations feasible and allows to predict a wider instability
region.

As in the linear case, we may also use the phase
$\phi_{\varepsilon}=-\varepsilon\log\rho_{\varepsilon}$ in the aggregate
equation (66). Letting $\varepsilon\rightarrow 0$ and assuming
$\phi_{\varepsilon}\rightarrow\phi$, we formally obtain

$\partial_{t}\phi+{\bf
U}_{R}^{G}\cdot\nabla\phi+\textrm{tr}\left[\mathbb{D}_{R}^{G}(\nabla_{x}\phi\otimes\nabla_{x}\phi)\right]=0,$
(76)

because ${\bf U}_{\rho}^{\varepsilon}\to{\bf U}_{R}^{G}(\nabla\phi)$ (defined
in (72) and
$\mathbb{D}_{\rho}^{\varepsilon}\to\mathbb{D}_{R}^{G}(\nabla\phi)$, which is
the variance covariance matrix of (67). Therefore, Eq. (76), in opposition to
the linear case, is not the quadratic expansion near $\nabla\phi=0$ of
equation (70) except in the regime $R$ small.

###  3.2 The regime $R$ small

Let us now consider $R$ small in the sense of (31) with $\mathcal{S}=\rho$,
i.e.,

$0<R\ll l_{\rho}:=\dfrac{1}{\max\dfrac{|\nabla\rho|}{\rho}},$

we may expand $\rho$ as

$\rho(t,x+R\hat{v})\approx\rho(t,x)+R\hat{v}\cdot\nabla\rho(t,x),$

so that the normalization function of (63) is

$c(t,x)\approx\rho(t,x)|\mathbb{S}^{d-1}|.$

Therefore Eq. (62) becomes

$\dfrac{\partial f}{\partial t}(t,x,v,\hat{v})+{\bf v}\cdot\nabla
f(t,x,v,\hat{v})=\mu\,\Big{(}\rho(t,x)(\dfrac{1}{|\mathbb{S}^{d-1}|}+\dfrac{R}{|\mathbb{S}^{d-1}|}\hat{v}\cdot\dfrac{\nabla\rho(t,x)}{\rho(t,x)})-f(t,x,v,\hat{v})\Big{)}.$
(77)

The latter may be assimilated to a case in which $T[\rho]$ is evaluated in a
small perturbation of a nondimensionalized homogeneous configuration set equal
to $1$ as $\dfrac{R\nabla\rho}{\rho}$ is small being $R\ll l_{\rho}$. This is
exactly the regime in which the linear stability analysis is performed in
[33]. Considering (77) with (14), we have a linearized kinetic equation. Then
we obtain $G_{R}=1-R\hat{v}\cdot p$ and we recover that $H(p)>0$ is equivalent
to the condition $V/\mu R>1$.

In the case of (77) with (14), we also remark that

$\partial_{t}\rho^{0}+\nabla_{x}\cdot(\rho^{0}{\bf U}_{\rho}^{0})=0$ (78)

and that ${\bf U}_{\rho}^{0}={\bf U}_{R}^{G}(0)$ and
$\mathbb{D}_{\rho}^{0}=D_{p}^{2}H(0)$, so that like in the linear case, in
this linearized regime obtained for $R$ small, the eikonal equation (76) is a
quadratic expansion of (70).

###  3.3 Concentration profile

Concerning the concentration points dynamics, we remark that, like in the
linear case, if $\bar{x}_{i}$ is a maximum of $\rho$, then

$\dot{\bar{x}}_{i}(t)={\bf U}_{R}^{G}(\bar{x}_{i}(t),p)$

and, as $p=\nabla_{x}\varphi(\bar{x}_{i})=0$, then, as a consequence of (73)

$\dot{\bar{x}}(t)=\int_{\mathbb{S}^{d-1}}V_{\psi}(\hat{v})\hat{v}\,d\hat{v},$

that is a nonvanishing quantity in the case where $V_{\psi}(\hat{v})$ is not
even as a function of $\hat{v}$. In such a case it is possible to observe
moving patterns, as showed in [33]. Again, in the regime of small $R$, that is
a linearized case, when the formal aggregate limit is (78) and ${\bf
U}_{\rho}^{0}={\bf U}_{R}^{G}(0)$, the aggregate limit procedure and the WKB
analysis give exactly the same amount of information about the evolution of
the maxima.

However, the Hamilton-Jacobi equation (70) also gives the microscopic
concentration profile. In the stability regime, i.e. when (64) is satisfied,
we actually have that $p=0$ is a minimum of $H$ as $H(0)=0$,
$\nabla_{p}H(0)=0$ when $V_{\psi}$ is even and $D^{2}_{p}(0)>0$. Therefore, we
have $H>0$ which implies that $\varphi(t,x)$ decreases and possible initial
concentrations will disappear. In the instability regime, the situation is
more interesting and the prototype of the shape of the Hamiltonian is depicted
in Fig. 4. In fact, being $D^{2}_{p}(0)<0$, then $p=0$ is a maxima. Therefore,
there will be a range of values of $p$ where the Hamiltonian is negative. As
for $p\rightarrow\pm\infty$, $H\rightarrow\infty$, and if $H$ is $C^{1}$,
there will be a value $\bar{p}>0$ where $H(\pm\bar{p})=0$, and the slopes
$\pm\bar{p}$ determine a saw tooth stationary state since $p=0$ is an unstable
state for (70). This explains the numerical profile obtained in Fig. 4. When
$\nabla\varphi$ is small (which is the case near concentration points), then
$H(\nabla\varphi)<0$, meaning that $\varphi$ will increase and the
concentration will get stronger.

Figure 4:  The convex-concave Hamiltonian $H$ in the one-dimensional case as
computed explicitly in formula (79). The values of $\mu,V,R$ are those
corresponding to the example in section 3.4. The red dashed lines correspond
to the asymptotes $\pm U|p|$, the black horizontal dashed line is the zero-
level. The values $\pm\bar{p}\neq 0$ with $H(\bar{p})=0$ determines the slopes
of the saw tooth solutions in Fig. 5

###  3.4 An example

As an example, we consider the 1D eigenproblem (71) (when $\hat{v}=\pm 1$).
Then $H$ is defined by

$1=\dfrac{\mu}{(\exp^{-Rp}+\exp^{Rp})}\int_{0}^{U}\psi(v)\left(\dfrac{\exp^{-Rp}}{\mu+H-vp}+\dfrac{\exp^{Rp}}{\mu+H+vp}\right)\,dv.$

Therefore

$1=\dfrac{\mu}{2(\exp^{-Rp}+\exp^{Rp})}\int_{0}^{U}\psi(v)\left(\dfrac{\exp^{-Rp}(\mu+H+vp)+\exp^{Rp}(\mu+H-vp)}{(\mu+H)^{2}-v^{2}p^{2}}\right)\,dv.$

When $\psi(v|\hat{v})=\delta(v-\hat{v}V)$, this reduces to

$1=\dfrac{\mu}{(\exp^{-Rp}+\exp^{Rp})}\left(\dfrac{\exp^{-Rp}(\mu+H+Vp)+\exp^{Rp}(\mu+H-Vp)}{(\mu+H)^{2}-V^{2}p^{2}}\right),$

and, therefore

$H^{2}+H\mu-V^{2}p^{2}+\mu VpD_{R}(p)=0,\qquad
D_{R}(p)=\dfrac{\exp^{Rp}-\exp^{-Rp}}{\exp^{-Rp}+\exp^{Rp}}=\tanh(Rp),$

in such a way that

$H(p)=\dfrac{-\mu+\sqrt{\mu^{2}+4V^{2}p^{2}-4\mu VpD_{R}(p)}}{2}.$ (79)

Consequently, the sign of $H$ is determined by the sign of
$V|p|-\mu|D_{R}(p)|$ and we obtain

$H(p)>0\qquad\textrm{iff}\qquad\dfrac{V}{\mu}>\dfrac{D_{R}(Rp)}{p}.$

For $|p|$ small this is in accordance with the sign of the second derivative
in formula (74). In the regime when $R|p|$ is small then
$\dfrac{D_{R}(Rp)}{p}=\dfrac{\tanh(Rp)}{p}\sim R$ and the latter condition
becomes (64). As $p=\nabla\varphi$ and
$|\nabla\varphi|\approx\dfrac{|\nabla\rho|}{\rho}$, then it coherently
corresponds to the analysis performed in the regime of $R$ small.

We now show some numerical tests. We solve numerically the kinetic equation
(62) in the regime (14). In particular we choose the following parameter
values: $V=1,\mu=100,R=5\cdot 10^{-2}$. Therefore $\varepsilon=10^{-2}$ and we
are in the regime of linear instability as $\dfrac{V}{\mu R}=0.2$. We consider
three different initial conditions: (a) a perturbation of the homogeneous
configuration, (b) a bimodal gaussian (i.e. $\rho^{0}$ as in (61)) centered in
$\bar{x}_{1}=2.3,\bar{x}_{2}=2.7$, (c) a bimodal gaussian (i.e. $\rho^{0}$ as
in (61)) with $\bar{x}_{1}=2.4,\bar{x}_{2}=2.6$.

As we are in a regime of linear instability, in Fig.5(a) we observe pattern
formation, while in figure (b), as the two initial peaks are far enough, they
stay so along the dynamics. In Fig. 5(c) we have that the two peaks merge,
because the sensing radius is large enough. In the second line (Fig.
5(d)-(e)-(f), respectively), we plot the corresponding $-\log(\rho)$.

Figure 5:  One dimensional example with parameters $V=1,\mu=100,R=5\cdot
10^{-2}$. In (a) the initial condition is the perturbation of the homogeneous
steady profile, in (b) the initial condition is a bimodal gaussian centered in
$\bar{x}_{1}=2.3,\bar{x}_{2}=2.7$, in (c) the initial condition is a bimodal
gaussian (i.e. $\rho_{0}$ as in (61)) centered in
$\bar{x}_{1}=2.4,\bar{x}_{2}=2.6$. In the second line we report the
corresponding profiles of $-\log(\rho(t,x))$. These saw tooth curves result
from the convex-concave Hamiltonian in fig. 4

##  4 Conclusion

We have considered a kinetic equation with a BGK relaxation operator in which
the transition probability is nonlocal in the physical space which can be
linear or nonlinear. For high frequencies, both in the localized and
nonlocalized regime, highly concentrated patterns may occur. We analyze them
thanks to the WKB ansatz, following [6], and obtain a Hamilton-Jacobi equation
in the limit. This method, in the spirit of adaptive dynamics, provides us
with the position evolution of the concentration points and with the
concentration profile. We can conclude that

  * •

In the linear case, the dynamics is driven by an external field
$\mathcal{S}(x)$. At the leading (hyperbolic) order, the kinetic equation and
corresponding aggregate limit almost give a complete and coinciding
information as ell as the WKB method through the canonical equation for the
maxima point. Indeed, we know the position of the maxima; the concentration
points do not depend on the initial condition but only on the field
$\mathcal{S}$ as the asymptotic equilibrium is determined by the H theorem and
it is independent on the initial condition. The WKB analysis allows also to
state the concentration profile around the maxima, in particular when the
leading order average velocity of $T$ does not vanish.

  * •

In the nonlinear case, the same link holds between the aggregate limit
analysis and the WKB analysis in the linearized regime ($R$ small). When $R$
is not small (and $\rho$ is not smooth) and we expect concentrations, then
stating (even formally) the macroscopic limit is not banal. Then, the WKB
method gives rise to an unusual convex-concave Hamiltonian, explaining saw
tooth patterns which are obtained numerically. It also tells us more with
respect to the aggregate equation through the canonical equation for the
maxima. A difference lies in the fact that it is not possible to find a
condition for the location of the concentration points. Moreover, it is
possible to recover the linear stability condition found in [33] in the
special case of one dimension and $\psi$ a Dirac delta in the limit of $R$
small where $\psi$ being the Dirac delta was chosen in order to determine
explicitly the instability condition. Furthermore, the present analysis
actually shows that it is possible to extend the results to other
distributions $\psi$ that have a larger second moment with respect to the
Dirac delta, that is the most unstable one, in the sense that it prescribes a
larger parameters region where we have linear instability. Moreover, the
analysis can be done in any dimension and the study of the positivity of the
Hessian matrix of the Hamiltonian allows to state the same result in any
dimension.

Both in the linear and nonlinear cases the aggregate and the WKB analysis give
compatible eikonal equations in the suitable limit. However, in the nonlinear
case the location of the concentration points cannot be explicitly determined.
The WKB analysis, as a further contribution, allows to determine a
concentration profile. In the fully nonlinear case ($R$ not small), the WKB
analysis can be used in order to understand the dynamics as the analysis of
the aggregate limits is not banal and as it goes beyond the regime of the
linear stability analysis. In the particular case analysed here we obtain the
same results, but we can expect that in other choices of transition
probability $T$ the WKB analysis predicts a different region. Moreover, the
WKB can be performed only in the (14) regime, but $R$ in the dimensional
system is not needed to be small like in the linear stability analysis. These
considerations suggest that the WKB analysis could be used in order to extend
the results of a linear stability analysis to other transition probabilities.

In the context of the study of kinetic eikonal equations in the same spirit as
[6], this work has allowed to make some steps further as $i)$ the Markovian
probability in the relaxation operator depends on the spatial variable as it
is nonlocal in the physical space, $ii)$ it was applied in order to study the
space dependent equilibrium in a regime in which concentrations are shown,
$iii)$ a regime in which the WKB and aggregate limit procedure may commute was
detected.

Moreover, as $T$ depends on $x$ and on the small parameter $\varepsilon$, in
the linear case $H$ depends on both $x$ and the $\nabla\varphi$, then we
obtain a time an evolution equation (41) for $H_{\varepsilon}$ and, then, a
time dependent eigenvalue problem similar to the principal bundle for
parabolic equations. Another open problem is to determine the boundary
conditions for the Hamilton-Jacobi equations.

### Acknowledgements

N.L. is member of INdAM-GNFM. N.L. acknowledges support by the Italian
Ministry for Education, University and Research (MUR) through the
“Dipartimenti di Eccellenza” Programme (2018- 2023) of the Department of
Mathematical Sciences, G. L. Lagrange, Politecnico di Torino (CUP:
E11G18000350001). N.L. gratefully acknowledges support from the Italian
Ministry of University and Research (MUR) through the grant PRIN2022-PNRR
project (No. P2022Z7ZAJ) “A Unitary Mathematical Framework for Modelling
Muscular Dystrophies” (CUP: E53D23018070001). N.L. gratefully acknowledges
support from the CNRS International Research Project ‘Modélisation de la
biomécanique cellulaire et tissulaire’ (MOCETIBI).

## Appendix. Boundary conditions

In Section 1.2 we have shown the boundary conditions for the hyperbolic limit
(23) of the kinetic equation (1) that are given by (24). In particular, both
(23) and (24) are derived from the kinetic equation: (23) is derived from (13)
and (24) are derived from (11) that is satisfied by any $f$ that obeys kinetic
boundary conditions in the form (10). The kinetic boundary conditions (10) are
actually imposed on the entering boundary, i.e. on $\Gamma_{-}(x)$. The
derived aggregate boundary conditions are noflux boundary conditions for the
conservation law (23). Actually, we want to verify that those boundary
conditions are to be imposed on the entering zone only, i.e. for
$x\in\partial\Omega$ such that ${\bf
U}_{\mathcal{S}}^{0}(x)\cdot\boldsymbol{n}(x)<0$ as in the outgoing region,
i.e. for $x\in\partial\Omega$ such that ${\bf
U}_{\mathcal{S}}^{0}\cdot\boldsymbol{n}>0$ they are granted by the underlying
kinetic boundary conditions.

Therefore, we need to compute the average ${\bf U}_{\mathcal{S}}^{0}(x)$ for
$x\in\partial\Omega$ that we denote as ${\bf
U}_{\mathcal{S}_{|\partial\Omega}}^{0}$. First of all we need to define
$T[\mathcal{S}]_{0_{|\partial\Omega}}$. Working in the regime (20), we define
it as

$T[\mathcal{S}]_{0_{|\partial\Omega}}(v,\hat{v})=c(x)b(\mathcal{S}(x+R(x,\hat{v})\hat{v}))\psi_{|\partial\Omega}(v|\hat{v}),$

where $R(x,\hat{v})$ is defined in (7) and $\psi_{|\partial\Omega}(v|\hat{v})$
is to be dependent on $x\in\partial\Omega$, as for
$\hat{v}\cdot\boldsymbol{n}(x)>0$, then we should set $V_{\psi}=0.$ Therefore,
we have that

$T[\mathcal{S}]_{0_{|\partial\Omega}}(v,\hat{v})=\begin{cases}\dfrac{\psi_{|\partial\Omega}(v|\hat{v})}{|\mathbb{S}^{d-1}|}\qquad&\textrm{if}\quad\hat{v}\cdot\boldsymbol{n}>0,\\\
c(x)b(\mathcal{S}(x+R\hat{v}))\psi_{|\partial\Omega}(v|\hat{v})\qquad&\textrm{if}\quad\hat{v}\cdot\boldsymbol{n}<0,\end{cases}$

as when $\hat{v}\cdot\boldsymbol{n(x)}>0$ then $R(x,\hat{v})=0$, while when
$\hat{v}\cdot\boldsymbol{n}<0$, $x+R\hat{v}\in\Omega$ if, for example $\Omega$
is convex.

Now, as $T[\mathcal{S}]_{0}$ is in fact the equilibrium, it must satisfy the
boundary conditions (10). We analyse the two cases $\alpha=0$ (purely
Maxwellian) and $\alpha=1$ (pure reflection), any case in between follows as a
convex combination. If we consider the Maxwellian boundary conditions, then we
have

$M(x,v,\hat{v})=c(x)b(\mathcal{S}(x+R(x,\hat{v})\hat{v})\psi(v|\hat{v}),\qquad\hat{v}\in\Gamma_{-}(x).$

and

$\displaystyle{\bf U}_{\mathcal{S}_{\partial\Omega}}^{0}$
$\displaystyle=\left[\int_{0}^{U}v\left(\int_{\hat{v}\cdot\boldsymbol{n}<0}M(x,v,\hat{v})\hat{v}d\hat{v}+\int_{\hat{v}\cdot\boldsymbol{n}>0}\dfrac{1}{|\mathbb{S}^{d-1}|}\psi(v|\hat{v})\hat{v}d\hat{v}\right)dv\right]$
$\displaystyle=\int_{0}^{U}v\int_{\hat{v}\cdot\boldsymbol{n}<0}M(x,v,\hat{v})\hat{v}d\hat{v}dv.$

Therefore

${\bf
U}_{\mathcal{S}_{\partial\Omega}}^{0}\cdot\boldsymbol{n}=\int_{0}^{U}v\int_{\hat{v}\cdot\boldsymbol{n}<0}M(x,v,\hat{v})\hat{v}\cdot\boldsymbol{n}d\hat{v}dv<0.$

In conclusion the whole boundary is an entering zone and then we need to
impose (22). If $\alpha=1$, then $T[\mathcal{S}]_{0}$ must satisfy the
following boundary conditions, if $\hat{v}\in\Gamma_{-}(x)$

$T[\mathcal{S}](v,\hat{v})_{\Gamma_{-}}=T[\mathcal{S}](v,\mathcal{W}(\hat{v})),$

and here $\mathcal{W}(\hat{v})\cdot\boldsymbol{n}>0$. Therefore
$T[\mathcal{S}]_{|\Gamma_{-}}=\dfrac{\psi(v|\hat{v})}{|\mathbb{S}^{d-1}|}$. In
conclusion, following the same computations as for the case $\alpha=0$, we
find ${\bf U}_{\mathcal{S}}^{0}=0,$ i.e. the velocity vector vanishes on the
whole boundary and (22) is satisfied.

Conversely, in the regime (14), we have that if $T[\mathcal{S}]_{0}$ does not
depend on $\hat{v}$ (and on x) because of the localization, then ${\bf
U}_{\mathcal{S}}^{0}=0$ on $\partial\Omega$.

## References

  * [1] N. J. Armstrong, K. J. Painter, and J. A. Sherratt.  A continuum approach to modelling cell-cell adhesion.  Journal of theoretical biology, 243 1:98–113, 2006.
  * [2] G. Barles.  Solutions de viscosité des équations de Hamilton-Jacobi.  Springer-Verlag Berlin Heidelberg, 1994.
  * [3] G. Barles and B. Perthame.  Concentrations and constrained Hamilton-Jacobi equations arising in adpative dynamics.  Contemporary Mathematics, 439:57, 2007.
  * [4] E. Bouin.  A hamilton-jacobi approach for front propagation in kinetic equations, 2015.
  * [5] E. Bouin and N. Caillerie.  Spreading in kinetic reaction–transport equations in higher velocity dimensions.  European Journal of Applied Mathematics, 30(2):219–247, 2019.
  * [6] E. Bouin and V. Calvez.  A kinetic eikonal equation.  Comptes Rendus Mathematique, 350(5):243–248, 2012.
  * [7] E. Bouin, V. Calvez, E. Grenier, and G. Nadin.  Large-scale asymptotics of velocity-jump processes and nonlocal hamilton–jacobi equations.  Journal of the London Mathematical Society, 108(1):141–189, 2023\.
  * [8] E. Bouin, V. Calvez, and G. Nadin.  Propagation in a kinetic reaction-transport equation: Travelling waves and accelerating fronts.  Archive for Rational Mechanics and Analysis, 217, 08 2015.
  * [9] N. Caillerie.  Large deviations of a forced velocity-jump process with a Hamilton–Jacobi approach.  Annales de l’Institut Fourier, 71(4):1733–1755, 2021.
  * [10] V. Calvez, G. Raoul, and C. Schmeiser.  Confinement by biased velocity jumps: Aggregation of escherichia coli.  Kinetic and Related Models, 8(4):651–666, 2015.
  * [11] C. Cercignani.  The Boltzmann Equation and its Applications.  Springer, New York, 1987.
  * [12] F. A. C. C. Chalub, P. A. Markowich, B. Perthame, and C. Schmeiser.  Kinetic models for chemotaxis and their drift-diffusion limits.  Monatshefte für Mathematik, 142(1):123–141, Jun 2004.
  * [13] A. Chauviere, T. Hillen, and L. Preziosi.  Modeling cell movement in anisotropic and heterogeneous network tissues.  Networks $\&$ Heterogeneous Media, 2(2):333, 2007.
  * [14] L. Chen, K. Painter, C. Surulescu, and A. Zhigun.  Mathematical models for cell migration: a non-local perspective.  Philosophical Transactions of the Royal Society B: Biological Sciences, 375(1807):20190379, 2020.
  * [15] M. Conte and N. Loy.  A non-local kinetic model for cell migration: A study of the interplay between contact guidance and steric hindrance.  SIAM Journal on Applied Mathematics, 0(0):S429–S451, 0.
  * [16] M. Conte and N. Loy.  Multi-cue kinetic model with non-local sensing for cell migration on a fiber network with chemotaxis.  Bull Math Biol., 84(3), 2022.
  * [17] M. G. Crandall, H. Ishii, and P.-L. Lions.  User’s guide to viscosity solutions of second order partial differential equations.  Bulletin of the American Mathematical Society, 27(1):1–67, 1992\.
  * [18] L. Evans.  The perturbed test function method for viscosity solutions of nonlinear PDE.  Proc. Roy. Soc. Edinburgh Sect. A, 111(3-4):359–375, 1989.
  * [19] F. Filbet, P. Laurencot, and B. Perthame.  Derivation of hyperbolic models for chemosensitive movement.  Journal of Mathematical Biology, 50:189–207, 03 2005.
  * [20] F. Filbet and N. Vauchelet.  Numerical simulation of a kinetic model for chemotaxis.  Kinetic and Related Models, 3:B348–B366, 09 2010.
  * [21] T. Hillen.  M5 mesoscopic and macroscopic models for mesenchymal motion.  Journal of mathematical biology, 53:585–616, 11 2006.
  * [22] T. Hillen and H. G. Othmer.  The diffusion limit of transport equations derived from velocity-jump processes.  SIAM Journal of Applied Mathematics, 61:751–775, 2000.
  * [23] T. Hillen, K. Painter, and C. Schmeiser.  Global existence for chemotaxis with finite sampling radius, 2007.
  * [24] E. F. Keller and L. A. Segel.  Initiation of slime mold aggregation viewed as an instability.  Journal of Theoretical Biology, 26(3):399 – 415, 1970.
  * [25] K. Lam and Y. Lou.  Introduction to Reaction-Diffusion Equations.  Lecture Notes on Mathematical Modelling in the Life Sciences. Springer International Publishing, 2022.
  * [26] K. Lam, Y. Lou, and B. Perthame.  A Hamilton-Jacobi approach to evolution of dispersal.  Communications in Partial Differential Equations, 48(1):86–118, 2023.
  * [27] B. Lods.  Semigroup generation properties of streaming operators with noncontractive boundary conditions.  Mathematical and Computer Modelling, 42:1441–1462, 12 2005.
  * [28] Tommaso Lorenzi and Camille Pouchol.  Asymptotic analysis of selection-mutation models in the presence of multiple fitness peaks.  Nonlinearity, 33(11):5791, oct 2020.
  * [29] A. Lorz, S. Mirrahimi, and B. Perthame.  Dirac mass dynamics in multidimensional nonlocal parabolic equations.  Communications in Partial Differential Equations, 36(6):1071–1098, 2011.
  * [30] N. Loy, T. Hillen, and K. Painter.  Direction dependent turning leads to anisotropic diffusion and persistence.  European Journal of Applied Mathematics, 33(4):729–765, 2022.
  * [31] N. Loy and L. Preziosi.  Kinetic models with non-local sensing determining cell polarization and speed according to independent cues.  Journal of Mathematical Biology, 80:373–421, 2020.
  * [32] N. Loy and L. Preziosi.  Modelling physical limits of migration by a kinetic model with non-local sensing.  Journal of Mathematical Biology, 80, 2020.
  * [33] N. Loy and L. Preziosi.  Stability of a non-local kinetic model for cell migration with density dependent orientation bias.  Kinetic and Related Models, 13(5):1007–1027, 2020.
  * [34] H. Othmer and T. Hillen.  The diffusion limit of transport equations ii: Chemotaxis equations.  SIAM Journal of Applied Mathematics, 62:1222–1250, 04 2002.
  * [35] H. G. Othmer, S. R. Dunbar, and W. Alt.  Models of dispersal in biological systems.  Journal of Mathematical Biology, 26(3):263–298, Jun 1988.
  * [36] R.G. Plaza.  Derivation of a bacterial nutrient-taxis system with doubly degenerate cross-diffusion as the parabolic limit of a velocity-jump process.  Journal of mathematical biology, 2019.
  * [37] D. W. Stroock.  Some stochastic processes which arise from a model of the motion of a bacterium.  Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 28(4):305–315, Dec 1974.