# On the Visual-based Safe Landing of UAVs in Populated Areas: a Crucial
Aspect for Urban Deployment.

Javier González-Trejo, Diego Mercado-Ravell*, Israel Becerra and Rafael
Murrieta-Cid Manuscript received: February, 24, 2021; Revised June, 5, 2021;
Accepted July, 13, 2021.This paper was recommended for publication by Editor
Pauline Pounds upon evaluation of the Associate Editor and Reviewers’
comments. This work was supported by the Mexican National Council of Science
and Technology CONACyT, the FORDECyT project 296737 “Consorcio en Inteligencia
Artificial”, and Intel Corporation. Authors are with the Center for Research
in Mathematics CIMAT AC, Mexico. D. Mercado-Ravell & I. Becerra are also with
Cátedras CONACyT (*corresponding author email: diego.mercado@cimat.mx).Digital
Object Identifier (DOI): see top of this page.

###### Abstract

Autonomous landing of Unmanned Aerial Vehicles (UAVs) in crowded scenarios is
crucial for successful deployment of UAVs in populated areas, particularly in
emergency landing situations where the highest priority is to avoid hurting
people. In this work, a new visual-based algorithm for identifying Safe
Landing Zones (SLZ) in crowded scenarios is proposed, considering a camera
mounted on a UAV, where the people in the scene move with unknown dynamics. To
do so, a density map is generated for each image frame using a Deep Neural
Network, from where a binary occupancy map is obtained aiming to overestimate
the people’s location for security reasons. Then, the occupancy map is
projected to the head’s plane, and the SLZ candidates are obtained as circular
regions with a minimum security radius, within the head’s plane. Finally, to
keep track of the SLZ candidates, a multiple instance tracking algorithm is
implemented using Kalman Filters along with the Hungarian algorithm for data
association. Several scenarios were studied to prove the validity of the
proposed strategy, including public datasets and real uncontrolled scenarios
with people moving in public squares, taken from a UAV in flight. The study
showed promising results in the search of preventing the UAV from hurting
people during emergency landing.

###### Index Terms:

UAVs, Autonomous Landing, Crowds Detection, Visual-based Landing, Density
Maps.

##  I Introduction

Over the last decades, the world have witnessed a great development and
production of Unmanned Aerial Vehicles (UAVs), flooding the market with all
kinds of drones for different budgets and with continuously increasing
capabilities. As a direct consequence, the number of successful civilian
applications has been steadily increasing, and it is more common to see UAVs
integrated in our daily lives, interacting closer to people in a wide variety
of tasks ranging from surveillance, photography, structure inspection, load
transportation, fast response to disasters, to simple entertainment, among
others [nex13_uav_mappin_applic], [lucieer13_mappin_lands_displ_using_struc],
[menouar17_uav_enabl_intel_trans_system_smart_city].

Nevertheless, the increased use of UAVs has also raised several well founded
legal and public security concerns, in part due to the threat of violating
privacy, but mainly derived from the inherent risk of an accident produced by
a system failure or human error. Although big efforts have been devoted to
increase the robustness of these systems, e.g., safety protocols for returning
home have helped to attenuate security issues, it is not enough to warrant
people’s safety on the streets in case of a system failure. Consequently, this
has restricted the use of drones to rural environments, or wide urban areas
clear of people, precluding some very interesting applications in populated
zones such as surveillance of crowds, monitoring mass events, fast response in
health emergencies, criminals chasing and law enforcing, etc.
[motlagh17_uav_based_iot_platf].

Accordingly, the legislation generally prohibits any UAV activity over a crowd
[lawDrones], but in practical terms, it is unavoidable that UAVs will flight
above an unaware crowd putting at risk the people’s safety in case of any UAV
malfunctioning, such as communication loss or power shortage, or even a human
error. Henceforth, it is of crucial relevance to further provide every UAV
with safety protocols to at least avoid hurting people in case of an
emergency, being able to autonomously find landing zones without putting
people at risk. If every UAV were provided with such an emergency autonomous
landing system, this would not only help to prevent hurting, or even killing
people during drones’ crashes, but it also would considerably boost all their
potential for successful urban deployment, even in crowded scenarios.

Figure 1: Autonomous visual-based Safe Landing Zone (SLZ) proposals in
populated areas. In the case of a system failure, UAVs must at least avoid
hurting people. The red regions are places where people is detected using
density maps. The green circles represent the regions free of people where the
drone can land. Finally, the blue circles indicate landing zones tracked by
KFs to ensure time consistency.

The idea of visual-based autonomous landing of flying machines has been
revisited in the past, and some examples are to be found mainly to locate
predefined landing tags [saripalliil_vision], [garcia2002towards], or to
identify and avoid hazardous terrain
[johnsonil_vision_guided_landin_auton_helic_hazar_terrain]. Nevertheless,
autonomous safe landing in crowded scenarios is an almost unexplored and very
challenging task, due to the implied security constraints and the inherent
difficulty of the problem, where the camera moves, and people may cluster in
dense groups producing severe occlusions, while freely moving in the scene
with unknown dynamics. It is only thanks to recent advancements in modern deep
learning techniques for computer vision that we can attempt to solve this
complicated task, where just a few related works are to be found in the
literature.

In [tzelepi19_graph_embed_convol_neural_networ], the authors propose a new
neural network and regularization technique based on the graph embedding
network to detect what they call no-fly zones, or zones where a crowd is
present in the scenario. Moreover, in [liu19_geomet_physic_const_drone_based]
a new neural network is proposed for time consistency across the frames in a
video, taking into account the perspective distortion caused by the projection
transformation from the real world to the image. The Safe Landing Zones (SLZ)
are a byproduct of the density consistency achieved in the work. More
recently, in [castellano20_crowd_detec_aerial_images_using], a model with a
MobileNet back-end is presented, specifically tailored for images coming from
a UAV, and fast enough for real-time applications. However, all these works
are limited to segment the scene in areas free of people, although this is a
necessary and important first step for the autonomous safe landing in crowded
areas, none of them attempts to solve the problem of indicating which regions
are the most appropriate to safely land the vehicle. It is important to stress
out that we do not use classical methods to avoid hazardous areas, because in
this work the priority is to avoid landing over people and not to ensure the
drone’s integrity. Nevertheless, the best overall solution for autonomous
landing during urban deployment would be a combination of both, namely,
approaches to detect flat surfaces
[johnsonil_vision_guided_landin_auton_helic_hazar_terrain] whenever possible
and crowd detection to prevent the vehicle from harming people.

In that regard, this paper presents an algorithm to detect Safe Landing Zones
(SLZ) in real populated and uncontrolled scenarios, considering a camera
mounted on a UAV, and where people are moving with unknown dynamics,
prioritizing their integrity in emergency situations, as depicted in Fig. 1.
To do so, we consider the use of a Deep Neural Network (DNN) density map
generator to infer people’s location on the image, then, an occupancy binary
mask is obtained aiming at overestimating the people’s location for security
reasons. In our previous work in [gonzalez-
trejo21_light_densit_map_archit_uavs], we proposed and implemented on a mobile
phone a lightweight architecture for real-time density map generation. A
complete study on the performance of the density map generator can be found
there. Moreover, as a first approach to detect SLZ, the largest circle free of
people was found directly on the image plane using the Polylabel algorithm
[polylabel], without considering the camera movement. In contrast, in this
work, the occupancy map is projected to the head’s plane in the real world,
which is more appropriate for the landing task, where multiple SLZ candidates
are obtained as circular regions. Note that this also contrasts with the
solution in [gonzalez-trejo21_light_densit_map_archit_uavs], where only a
single SLZ is computed. In addition to this, we have also considerably
improved the time response of SLZ detection by changing the algorithm to find
a suitable SLZ in the person-free space in the heads’ plane, that is, by
replacing the Polylabel algorithm [polylabel] with the Euclidean distance
transform [felzenszwalb2012distance].

Furthermore, by means of Kalman Filters (KF) and the Hungarian algorithm for
data association [bruff2005assignment], a multiple instances tracking
technique is implemented using the SLZ candidates as measurements, smoothing
the SLZ proposals movement and providing time consistency along frames.
Between two time instants $k$ and $k+1$, a KF for each landing zone gives an
estimate of the location and radius of a landing zone in the head’s plane. The
Hungarian algorithm finds the most similar landing zone at time $k+1$ for a
given zone at time $k$; that landing zone will be the one that maximizes the
area of Intersection over Union (IoU) of the landing zones at times $k$ and
$k+1$. Indeed, the Hungarian algorithm orders the similarity of all possible
matching using that criterion. Having multiple suitable SLZ allows to study
different criteria to choose “the best” SLZ, based, for instance, on their
size, their time consistency or their distance to the UAV.

Several scenarios were studied from available crowd datasets. Some metrics for
the SLZ detection task were defined and three different solutions were
compared. Lastly, real-time tests were carried out during tele-operated
landing missions in real uncontrolled and challenging scenarios, from a flying
drone in public squares with people moving at will. Although there is still
room for improvement, the proposed strategy showed good performance in spite
of the challenging scenarios, proving to be helpful to prevent hurting people
in case of emergency landing.

The reminder of the paper is structured as follows. In Section II, we briefly
introduce the density map generators and their utility for detecting crowds in
images. Section III describes the SLZ proposal generator in the head’s plane.
Section IV presents the multiple instances tracking algorithm. In Section V,
the experimental results are presented. Finally, in Section VI, conclusions
and future work are discussed.

##  II Crowd Detection using Density Maps

In order to find suitable landing zones where no person could be harmed, it is
necessary to distinguish the crowd from a SLZ. One of the most popular methods
for crowd detection are the density map generators [kang2018beyond]. Density
maps are obtained by means of DNN trained utilizing crowd images containing
head annotations, and provide spatial information on the crowd location even
subject to large occlusions, furthermore, they provide an estimation of the
number of persons in an image or scene as the integral of the corresponding
density map [zhang16_singl_image_crowd_count_multi]. These density maps have
excellent results in the counting task for a wide variety of crowded scenes,
ranging from sparse to highly dense crowds, where person-to-person occlusion
could be so severe to the point that the only visible part of a person is a
portion of the head.

In contrast to the classical detect-and-count approaches
[bochkovskiy2020yolov4], [tan2019efficientdet], which try to detect the whole
body for each person, the density map generators are trained to detect the
heads of persons in a crowd [zhang16_singl_image_crowd_count_multi],
increasing the robustness against occlusion, at the cost of being more prone
to overestimate the people’s locations under rich textured backgrounds. The
first works using density maps with deep learning algorithms are found in
[zhang16_singl_image_crowd_count_multi], with the algorithm Multi Column
Convolutional Neural Network (MCNN) that used multiple columns to capture
different perspectives and heads’ scales. The latest works focus not only on
obtaining more accurate architectures, but also on improving the manual
annotations through learning [bai2020adaptive], or creating new loss functions
to better learn how to count and locate the crowd in an image
[wang20_distr_match_crowd_count].

The accuracy of people’s location is of great relevance for the SLZ detection
task in order to avoid human accidents. In that sense, due to the inherent
security constraints, it is preferred to overestimate the people’s location,
hence, it is preferred to have false positive detections, but it is
unacceptable to miss detection that can result in hurting people. Henceforth,
the employed density map generator should be explicitly trained to
overestimate the people’s locations to increment the security.

Another important aspect to take into account when selecting the density map
generator is its computational complexity and time response. It is of crucial
importance to have the information at a high rate in order to perform the
landing maneuver in real-time, preferably embedded on the vehicle. Current
state-of-the-art algorithms, especially the works using pre-trained DNN back-
ends like VGG (Visual Geometry Group), are slow, mostly in embedded devices,
which hinders the objective of implementing them in most UAVs. Accordingly, it
is required to have a good trade-off between accuracy and computational cost
when selecting the density map generator, where recent advancements in
lightweight architectures appear as a good alternative for the SLZ detection
task. These lightweight architectures could be defined as DNN below the 1
million of parameters, being suitable to be implemented in real-time solutions
for both high-end computers and embedded devices [gonzalez-
trejo21_light_densit_map_archit_uavs]. For this work, we have selected our own
custom-designed lightweight density map generator called Pruned Bayes Loss
Compact Convolutional Neural Network (Pruned BL CCNN) [gonzalez-
trejo21_light_densit_map_archit_uavs], but any other state-of-the art density
map generator can be used instead depending on the available hardware.

The Pruned BL CCNN has a lightweight density map architecture developed from a
pruned density map generator to reduce the number of parameters, and trained
using the Bayesian Loss to improve its accuracy, obtaining state-of-the-art
Mean Square Error accuracy, in addition to maintaining a reasonable number of
parameters. It obtained good results in the crowd counting and detection
tasks, while being light enough to be implemented for real-time application on
a mobile phone app. More specifically, that algorithm achieves $154.07$ as
Mean Absolute Error and $241.77$ for Mean Square Error, well within the state-
of-the-art for crowd detection and counting, while maintaining a number of
$0.067M$ parameters, which translates to $12.868$ giga flops to process an
image of resolution $856\times 480$. This network was explicitly trained to
overestimate the crowd for this particular task, in order to increase the
safety of the algorithm.

##  III Safe Landing Zone Proposals

Figure 2: UAV moving in a populated environment. Three reference frames are
defined, an inertial one fixed to the ground, a body fixed frame attached to
the UAV and a camera frame. Then, the camera pose can be retrieved from the
UAV on-board sensors.

The objective of our algorithm is to propose SLZ in populated environments to
prevent the drone from hurting people in case of a system failure. For that
matter, we consider a camera attached to a UAV. Moreover, people may freely
move in the scene with unknown dynamics. For now, we restrict our study to
scenes where people are primarily moving on a horizontal main plain.

In the present modeling, three reference frames are used. The first one is the
inertial or global reference frame, which is denoted as $\mathbf{F_{W}}$. The
body reference frame is attached to the UAV’s center of mass and is referred
as $\mathbf{F_{B}}$. The homogeneous transformation from the UAV’s body
reference frame $\mathbf{F_{B}}$ to the global reference frame
$\mathbf{F_{W}}$ is denoted by $\mathbf{T}^{W}_{B}\in SE(3)$ and is assumed to
be retrievable by means of the UAV’s on-board sensors. The third reference
frame uses the camera center as its origin and is denoted as $\mathbf{F_{C}}$,
and its transformation with respect to the body fixed frame
$\mathbf{T}^{B}_{C}\in SE(3)$ is known a priori. Fig. 2 shows the three
reference frames and their transformations.

An important element in the proposed solution is the so-called head plane,
$P$, which is a plane in the scene parallel to the floor that crosses the
average height of a person at $h_{h}$. More precisely, let $P$ be the plane
defined in $\mathbf{F_{W}}$, with $P_{0}=(0,0,h_{h})$ and
$\vec{\bf{n}}=(0,0,1)$ defining its point-normal form. The head plane $P$ is
crucial since computations of SLZ proposals are performed on it. Additionally,
let $I$ be the image plane as customary defined in $\mathbf{F_{C}}$.

The algorithm is divided into two main stages: the SLZ candidate proposal and
the multiple SLZ tracking that ensures time consistency between frames. The
former is described in this section, while the latter in Section IV. In the
SLZ proposal stage, we obtain at most $N_{p}$ SLZ candidates in the $P$ plane,
with a safety radius of at least $r_{0}$ meters. Towards the generation of
these SLZ proposals, we first calculate the density map $D$ in the image plane
$I$ using a density map generator. For our experiments, we use the lightweight
deep learning architecture Pruned BL CCNN [gonzalez-
trejo21_light_densit_map_archit_uavs].

The next step is to obtain an occupancy map $O$ in the $I$ plane, based on the
density map $D$. The map $O$ has the same size as $D$, where each pixel has
either a value of $255$ or $0$ denoting whether there is a person or not
according to the density values from $D$. For that, all occupancy map pixel
values $O_{i,j}$ are set according to the next expression:

$O_{i,j}=\begin{cases}255&D_{i,j}-\min\\{D_{i,j}\\}>0,\\\
0&D_{i,j}-\min\\{D_{i,j}\\}=0.\end{cases}\\\ $ (1)

Note that in the unlikely case where all the image is detected as a crowd
($\min\\{D_{i,j}\\}=255$), the SLZ algorithms would not be required and the
best strategy would be to fly away from that area. After that, the zones with
values equal to $255$ are dilated two times using a kernel of ones with size
$\mathbb{R}^{5\times 5}$ to further overestimate the crowd’s location for
security reasons. Finally, we perform a bitwise-not operation to set all
values equal to $0$ to $255$ and vice versa.

Thereafter, the occupancy map $O$ is projected to the head’s plane $P$ based
on the rigid transformations that relate reference frames $\mathbf{F_{W}}$,
$\mathbf{F_{B}}$ and $\mathbf{F_{C}}$, which are obtained from the UAV’s
sensors and the a priori known camera’s pose with respect to the UAV’s body.
These assumptions are consistent with the information commonly available from
most UAVs. More precisely, let $\mathbf{T}^{B}_{W}\in SE(3)$ be the rigid
transformation from inertial frame coordinates to body frame coordinates, and
$\mathbf{T}^{C}_{B}\in SE(3)$ the transform from body frame coordinates to
camera frame coordinates (see Fig. 2). Thus, denoting the intrinsic camera
matrix by $\mathbf{K}$, the projection from a point
$(x_{W},y_{W},h_{h},1)^{T}$ in the $P$ plane to a point $(x_{I},y_{I},1)^{T}$
in the image plane $I$ is given as:

$\lambda\left[\begin{array}[]{c}x_{I}\\\ y_{I}\\\
1\end{array}\right]=\mathbf{K}\mathbf{T}^{C}_{B}\mathbf{T}^{B}_{W}\left[\begin{array}[]{c}x_{W}\\\
y_{W}\\\ h_{h}\\\ 1\end{array}\right],$ (2)

with $\lambda$ a scale factor, which can be easily retrieved from the UAV
altitude and the camera pose.

Additionally, to perform the projection of the occupancy map $O$ to the head
plane $P$, denote it as $O^{W}$, a grid $G\in P$ is defined. Each element
$G_{i,j}$ of the grid stores a coordinate $(x_{W},y_{W})$ on the $P$ plane;
refer to them as $(G^{x}_{i,j},G^{y}_{i,j})$. Then, the corresponding
occupancy values for $O^{W}$ are retrieved as follows. Utilizing Eq. (2) and
given a grid element $G_{i,j}$, its stored inertial coordinates
$(G^{x}_{i,j},G^{y}_{i,j})$ are mapped to an image coordinate $(x_{I},y_{I})$
as:

$\lambda\cdot(x_{I},y_{I},1)^{T}=\mathbf{K}\mathbf{T}^{C}_{B}\mathbf{T}^{B}_{W}\cdot(G^{x}_{i,j},G^{y}_{i,j},h_{h},1)^{T}.$
(3)

The image coordinate $(x_{I},y_{I})$ is employed to query the occupancy map
$O$. Thus, the respective occupancy value in $O$ is used as the occupancy
value for $O^{W}_{i,j}$. The implementation of the past procedure is based on
the sampler module of the Spatial Transformer Network [jaderberg2015spatial].

In the final step, we obtain a distance map $C$ from the occupancy map $O^{W}$
in the head’s plane. The distance map $C$ considers the binary map encoded in
$O^{W}$, and stores the distance of each pixel with nonzero value (people
free) to the the nearest zero value pixel (occupied by people). By finding the
position and the value of the pixel, we effectively find the biggest circular
landing zone center and its radius in the head plane $P$. The algorithm used
to calculate the distance map is the Euclidean distance transform
[felzenszwalb2012distance]. Its algorithmic complexity is $O(aB)$, where $a$
is the dimension of the map, in this case $2$, and $B$ is the number of
elements on the map. Next, we mark the landing zone in the projected
occupation map $O^{W}$ as occupied, and repeat the process until a number of
$N_{p}$ SLZ candidates are encountered, or until a small landing zone with
radius less than $r_{0}$ is obtained. The results of the SLZ detection
algorithm in the image plane $I$ and in the homography projection to the
head’s plane $P$ are illustrated in Fig. 3.

Figure 3: Safe landing zones in the image plane (left) and its corresponding
projection in the head’s plane (right). The safe landing zones are found in
the head’s plane using the projected occupation map.

##  IV Multiple landing zone Tracking

By themselves, the safe landing zone proposals could be used to indicate where
a UAV could land. Nonetheless, due to the difficulty of the task, where not
only the camera moves, but also people in the scene move freely with unknown
dynamics within the plane $P$, in addition to the density maps not being
exempt to fails due to environmental variations, the SLZ candidates might
critically change from a frame in time $k$ to a frame in time $k+1$. This
could be catastrophic in the context of an emergency landing mission.
Therefore, as the second stage of our algorithm, we propose to use multiple
instance KF trackers to filter out outliers in the SLZ detection, smoothing
out the movement of the landing regions, preventing abrupt jumps and ensuring
temporal consistency along frames.

For a given instance $i$ of a SLZ encountered using the procedure described in
the previous section, let us consider a discrete time KF, with state vector at
time $k$ equal to

$\mathbf{x}_{k}^{(i)}=\left(x_{k}^{(i)},y_{k}^{(i)},r_{k}^{(i)},\dot{x}_{k}^{(i)},\dot{y}_{k}^{(i)},\dot{r}_{k}^{(i)}\right)^{T},$
(4)

where $x_{k}^{(i)}$ and $y_{k}^{(i)}$ represent the coordinates of the SLZ
circle’s center in the head plane $P$, $r_{k}^{(i)}$ is the circle radius,
while the velocity of the coordinates of the center and the rate of change of
the radius are given by $\dot{x}_{k}^{(i)}$, $\dot{y}_{k}^{(i)}$ and
$\dot{r}_{k}^{(i)}$, respectively.

For the dynamic model, we use a constant velocity model with acceleration
considered as part of the process noise $\mathbf{\omega}$ with normal
distribution [saho18_kalman_filter_movin_objec_track], that is,

$\mathbf{x}_{k+1}^{(i)}=\begin{bmatrix}\mathbf{I}_{3}&\Delta
t\mathbf{I}_{3}\\\
\mathbf{0}_{3}&\mathbf{I}_{3}\end{bmatrix}\mathbf{x}_{k}^{(i)}+\mathbf{\omega}\
\ \ |\ \ \ \mathbf{\omega}\sim\mathcal{N}(\mathbf{0},\,\mathbf{Q}),$ (5)

with $\mathbf{0}_{3}\in\mathbb{R}^{3\times 3}$ a matrix of zeros,
$\mathbf{I}_{3}\in\mathbb{R}^{3\times 3}$ the identity matrix, and $\Delta t$
the time increment. The process noise covariance matrix $\mathbf{Q}$ is
obtained as:

$\mathbf{Q}=\sigma_{a}\begin{bmatrix}\frac{\Delta
t^{4}}{4}\mathbf{I}_{3}&\frac{\Delta t^{3}}{2}\mathbf{I}_{3}\\\ \frac{\Delta
t^{3}}{2}\mathbf{I}_{3}&\Delta t^{2}\mathbf{I}_{3}\end{bmatrix},$ (6)

where $\sigma_{a}$ is the acceleration uncertainty and is selected
empirically.

The measurement model considers the observation vector $\mathbf{z}^{(i)}_{k}$
as the coordinates of the center of the SLZ and its radius, as obtained by the
algorithm described in Section III, with a measurement noise $\mathbf{\nu}$
with normal distribution, therefore, the measurement model is given by:

$\mathbf{z}^{(i)}_{k}=\begin{bmatrix}\mathbf{I}_{3}&\mathbf{0}_{3}\end{bmatrix}\mathbf{x}^{(i)}_{k}+\mathbf{\nu}\
\ \ |\ \ \ \mathbf{\nu}\sim\mathcal{N}(\mathbf{0},\,\mathbf{R}),$ (7)

with $\mathbf{R}$ standing for the measurement covariance matrix. Note that
this approach is based on a simplified model, but other alternative models and
tracking algorithms could be used instead at a computational cost.

To handle multiple target tracking and associate the main SLZ candidates with
their corresponding trackers, the Hungarian algorithm is employed, a
combinatorial optimization algorithm that solves the assignment problem in
polynomial time [bruff2005assignment]. We associate a SLZ proposal with a KF
by creating a cost matrix using the Intersection over Union (IoU) criteria,
between all combinations of the SLZ measurements and the KFs estimates. For
two circle areas $A_{i}$ and $A_{j}$, the IoU is computed as

$IoU=\frac{A_{i}\cap A_{j}}{A_{i}\cup A_{j}}.$ (8)

Then, the cost matrix is fed to the Hungarian algorithm, which associates the
“most similar” SLZ proposal for each KF tracker. If no match is found for a
KF, it is assumed that no consistent observation is available for that filter,
and a not-found counter is increased by one. When the counter is greater than
a threshold $\mu_{1}$, the respective KF is eliminated. If an unmatched SLZ
proposal is found for a number of $\mu_{2}$ consecutive iterations, and the
number of KFs is below $N_{p}$, a new KF is initialized.

At the beginning of a trial, a maximum of $N_{p}$ KFs are created by
initializing their states using the first SLZ proposals. In the following
frames, the propagation step is performed, and when available, the predicted
state is corrected through the corresponding observation.

|  |
---|---|---
|  |

a)

(a)

b)

(b)

c)

(c)

d)

(d)

e)

(e)

f)

(f)

Figure 4: Experiments results. The red zones are where the crowd is located,
the green circles are the current frame safe landing proposals and the blue
circles are the trackers. As shown in the images, the goal is to find suitable
spots for the drone to land in a variety of scenarios, from sparse crowds
(Fig. 4(a)) to different illumination conditions (Fig. 4(e)).

##  V Experimental Validation

The proposed strategy was validated in a wide variety of real scenarios
recorded from a drone in public squares, as well as with videos from the
Venice [Liu_2019_CVPR] and VisDrone [zhu2020vision] datasets. The former
scenarios were recorded by our lab members using a drone Parrot Bebop 2.

###  V-A Qualitative Analysis

As shown in Fig. 4, the scenes come from different points of view of distinct
public plazas taken from a UAV in flight, with people moving at will. These
recordings display a variety of characteristics that make them specially
challenging to test safe landing zone proposal algorithms. All of them contain
moving crowds with changes in illumination, heterogeneous surfaces on land,
different drone’s altitude and movement, and different backgrounds.

In the dataset experiments, we used three scenarios from the Venice dataset
[Liu_2019_CVPR] recorded in the San Marcos plaza. These videos were recorded
using a cellphone above the crowd at a resolution of $720$p at $25$ fps. Since
no camera pose was obtained, the homographies were estimated from the image
alone at every $60$ frames mark. Since the homography has an update rate of
$0.42$ Hz, we fix it in place for short periods of time where it does not
appear to be significant camera movement. All of the scenes present a crowd
moving on a flat surface, which make them suitable to test the SLZ proposal
algorithm despite being created for the crowd counting task.

Fig. 4 also shows the performance of our algorithm in different scenes. We see
in green the current frame SLZ candidates as computed by the algorithm
described in Section III and in blue the KF trackers estimate. In all figures,
it can be seen that the algorithm, for that frame, locates suitable SLZ
proposals and trackers in zones free of people. The challenge that each scene
offers varies depending on the characteristics described above. For example,
in Figure 4(e) we see one of the most challenging scenarios; where the camera
is near to the ground in a narrow passage, capturing people in a constant
flow, making it difficult to keep track of the landing zone proposals, while
presenting high illumination contrasts due to the sun. Despite of this, the
oldest tracker (signaled by a lower number) keeps track of umbrellas and a
small green area where the drone could land without hurting people. Also, in
Fig. 4(f), we display a scene from the Venice dataset, showing good
performance since all trackers cover a valid SLZ free of people.

Another interesting aspect can be observed in Figs. 4(a), 4(b), 4(c), where
sparser crowds are present in a wide scene. We can observe how the rich
textures on the floor may confuse the density map appearing as false positive
crowds. This is produced since the density map generators are explicitly
trained to overestimate the crowds to increase the algorithm’s security.
Nevertheless, the proposed strategy was able to find suitable SLZ in the
scenes. A more robust density map generator, such as DM-Count
[wang20_distr_match_crowd_count], could help to mitigate this phenomena at a
computational cost, compromising also the capacity to implement them in real-
time applications.

Finally, to better appreciate the performance of the strategy in dynamic
scenarios, we provide a video in the multimedia material of this paper,
containing the results of our algorithm in various scenarios, including SLZ
detection in real-time during landing missions using a human-piloted drone.

###  V-B Performance Analysis

Figure 5: Performance analysis on the VisDrone dataset. We evaluate the number of people entering the target SLZ (warning region), entering the risk area of $1m$ radius around the landing center (inner circle), the SLZ area and the computation time. Three different solutions are studied, our current proposal with the larger SLZ (blue circle), and with the older SLZ (green circle), as well as the solution from [gonzalez-trejo21_light_densit_map_archit_uavs] (red circle). Table I: Performance comparison for SLZ detection. work / Dataset | warning | danger | SLZ area | execution time [$s$] | IoU
---|---|---|---|---|---
| avg | $\sigma$ | avg | $\sigma$ | [$m^{2}$] | avg | max | min |
[gonzalez-trejo21_light_densit_map_archit_uavs] / Venice | $1.502$ | $6.239$ | $0.089$ | $0.435$ | $29.749$ | $1.0798$ | $2.495$ | $0.413$ | $0.358$
current (larger SLZ) / Venice | $0.185$ | $0.906$ | $0.029$ | $0.254$ | $\mathbf{191.645}$ | $\mathbf{0.105}$ | $\mathbf{0.208}$ | $\mathbf{0.092}$ | $\mathbf{0.611}$
current (older SLZ) / Venice | $\mathbf{0.077}$ | $\mathbf{0.451}$ | $\mathbf{0.005}$ | $\mathbf{0.077}$ | $124.470$ | $\mathbf{0.106}$ | $0.227$ | $\mathbf{0.090}$ | $0.537$
[gonzalez-trejo21_light_densit_map_archit_uavs] / VisDrone | $2.752$ | $4.658$ | $0.015$ | $0.177$ | $\mathbf{1796.85}$ | $1.051$ | $4.372$ | $0.193$ | $\mathbf{0.438}$
current (larger SLZ) / VisDrone | $2.646$ | $4.651$ | $\mathbf{0.004}$ | $0.072$ | $\mathbf{1716.06}$ | $0.143$ | $0.312$ | $\mathbf{0.109}$ | $\mathbf{0.432}$
current (older SLZ) / VisDrone | $\mathbf{2.123}$ | $\mathbf{4.085}$ | $\mathbf{0.004}$ | $\mathbf{0.066}$ | $1499.81$ | $\mathbf{0.127}$ | $\mathbf{0.255}$ | $\mathbf{0.102}$ | $0.428$

We further evaluate the performance of our SLZ detection algorithm and compare
it with our previous solution in [gonzalez-
trejo21_light_densit_map_archit_uavs], to our knowledge, the only other work
available aiming at solving the problem of SLZ detection in crowded scenarios.
To evaluate the proposed SLZS safety, we use two public crowds datasets,
Venice [Liu_2019_CVPR] and VisDrone [zhu2020vision]. The former containing
$167$ annotated frames from four videos, two minutes long, at resolution
$1280\times 720$, while the VisDrone dataset contains $2460$ annotated frames
at resolution $1920\times 1080$, from $82$ real scenarios taken from a drone.

We introduce a set of metrics and compare three different solutions, the one
presented at [gonzalez-trejo21_light_densit_map_archit_uavs], and two
different criteria to select “the best” SLZ in our current approach, taking
advantage of the availability of multiple SLZ candidates. The first criterion
consists in selecting the largest available SLZ, while the second corresponds
to the SLZ with the oldest tracker, aiming at finding more stable regions.

We are interested in measuring the level of safety of the SLZ that a drone
should target for landing. Hence, we define two levels of risk: warning,
corresponding to people entering the SLZ, and danger, corresponding to people
entering the safety minimum radius of $r_{0}=1m$ around the center of the SLZ.
Then, we count the average and standard deviation of the number of people
entering each risk area. To obtain the metric, we use the heads’ annotations
or coordinates of the people’s heads, and measure how many of them lay inside
a risk region. For a fair comparison, since our previous solution was only
found in the image plane $I$, we evaluate this metric on the image plane by
projecting the proposed SLZ of our current approach. The obtained results are
presented in Table I, where we can observe that the current strategy clearly
surpasses our previous solution in terms of people entering the danger and
warning regions. Moreover, we can observe that the oldest tracked SLZ improves
the results obtained with the largest one. This was expected since the oldest
tracker tends to be assigned to more stable regions in the image, where people
barely enter, not to mention the time consistency provided by the KF trackers.
These are two of the main advantages of our current solution, namely, being
able to select different criteria to find the most suitable SLZ. An example of
the results obtained in the VisDrone dataset is depicted in Fig. 5, where our
current approach is in the left, with the oldest SLZ represented with a green
circle, and the largest one with a blue circle, meanwhile, on the right image,
the solution from [gonzalez-trejo21_light_densit_map_archit_uavs] is presented
with a red circle.

We further evaluate the average size of the area of the target SLZ in $m^{2}$,
along each video sequence, provided that a larger SLZ implies a higher level
of security, specially in the case of people moving towards the SLZ center.
Because our previous solution does not propose the SLZ in the real world, we
project it and then measure its approximate area. VisDrone does not contain
the homography to correctly obtain the heads plane; however, the images are
taken from a camera looking straight down. Thus, we approximate the homography
as a simple scale transformation, where the scale is estimated from the size
of well-known objects such as cars. As expected, our current approach using
the larger SLZ performed the best considering both datasets.

Another interesting metric is the IoU with respect to one of the available
ground truth SLZ. To generate the ground truth SLZ, we use the head
annotations of the datasets and dilate them to a square area of $0.16m^{2}$,
that we call the personal space for each person. From this ground truth
occupation map, we obtain the 10 biggest SLZ using the distance transform
algorithm, as explained in Section III. Afterwards, for each frame, we obtain
the IoU for each ground truth SLZ and report the best of them, such that we
can get an idea of how close is the solution to the most similar ground truth
SLZ. As we can appreciate from Table I, our current proposal with the largest
SLZ was the best solution considering both datasets.

Finally, the computation time for the three different solutions is also
presented in Table I, where the three approaches were implemented in the same
computer described in the next subsection. In both solutions, the main
bottleneck is the computation of the density map by our Pruned BL CCNN. In
average, the model takes $0.044$s in the Venice dataset and $0.106$s in the
VisDrone dataset. Otherwise, an important improvement is accomplished by our
current proposal, in spite of the addition of the perspective transformation
and the KF tracking. This is mainly due to the expensive computation of the
Polylabel algorithm in our previous approach. These computation times are
suitable for real-time applications.

###  V-C Tele-operated Landing Task

Figure 6: Experimental results during a tele-operated landing task. The
proposed algorithm is implemented in real-time on a ground computer. The
strategy is able to propose suitable SLZ at different altitudes.

This section presents a proof of concept of a complete landing task. The
landing zones detection method and Kalman filtering were implemented to run
automatically in real-time in the ground computer utilizing the video streamed
from a Bebop 2 and its on-board sensor measurements, while the vehicle
guidance was performed by tele-operation. The drone was remotely controlled
using ROS from a ground station computer. The ground computer is equipped with
a Intel core i7 9750H paired with a RTX 2070 Mobile, 16 GB of Ram memory at
2400 MHz. It recovers data streamed from the UAV via WiFi communication, and
computes the SLZ detection algorithms in real-time at around $18$ Hz, having
as bottleneck the density map generation at $0.019$s. The video resolution was
set to $856\times 480$ at $15$ to $30$ fps, depending on the communication
quality. We used the IMU (Inertial Measurement Unit), GPS (Global Positioning
System), and altimeter filtered data available at $5$ Hz to obtain the drone’s
pose and use it to project the image to inertial frame coordinates. The camera
is mounted in the front of the drone with a wide-angle lens. Additionally, the
camera allows to select a region of interest with distortion correction by
software using virtual tilt and pan angles. The tilt angles range from
$80^{\circ}$, almost looking straight to the sky, to $-80^{\circ}$, almost
looking down straight to the floor.

During the trials, the drone’s camera was set with a tilt angle of
$-80^{\circ}$. Fig. 6 shows a sequence of images taken from the drone, at
different altitudes, as the vehicle descends to finally land. As it was done
before, the landing zones have a numeric label that serves as a unique
identifier; by means of this identifier, the evolution of each SLZ can be
tracked in Fig. 5. In the trial, the proposed algorithm is able to propose
valid SLZs during the whole landing mission, ranging from around $10m$ to
$2m$. Once the drone crosses the head plane, the algorithm is not longer
needed. A piece of the landing trails with a real drone is presented in the
provided video.

##  VI Conclusions and Future Work

In this paper, we present a novel algorithm for visual-based autonomous safe
landing of UAVs in populated environments. We believe that this kind of
algorithms is a key aspect for the successful application of UAVs in urban
areas, provided that they will help to prevent injuring people in case of
emergency landing.

The proposed algorithm, in conjunction with an adequate density map generator,
is capable of finding SLZ proposals in a variety of real uncontrolled and
challenging scenes, where the crowds move with unknown dynamics, the
illumination changes, and different crowd densities and rich image features
are present. Moreover, with the help of the KFs, the algorithm is able to keep
track of Safe Landing Zones even when the density map fails to correctly
detect parts of the crowds in consecutive frames.

Although the present work shows promising results, it is not exempt from
failure, and due to the hard security constraints involved and the wide
diversity of challenging situations that may occur during an emergency landing
mission, there is still room for improvement and further work is still
required to obtain a fully reliable solution. Future efforts will be devoted
to define a suitable criteria to choose the best SLZ available. Moreover, it
is required to impose a horizontality condition on the main planes of the
scene in order to discriminate unsuited regions for landing, such as walls or
steep terrains. Another interesting approach to further increase the people’s
security would be to add a human body detector as the UAV comes closer to the
ground. In addition, at present we focus all efforts in avoiding to hurt
people, but preventing the UAV from crashing into other obstacles would be
also of interest, this may be accomplished by recovering depth information
using stereo cameras or range sensors. Another interesting axis of research
would be to explore other kinds of filters, such as the Extended Kalman Filter
(EKF) or the Unscented Kaman Filter (UKF), coupled with more sophisticated
models better suited for the complex dynamics of the people on the scene.
Finally, the real-time capabilities could be further improved by having an
end-to-end deep learning solution trained using real or synthetic data.


}L^{\infty}(0,T_{0};H^{2}(\Omega_{\varepsilon})),\quad\partial_{t}\rho_{n}\rightharpoonup^{*}\partial_{t}\rho\text{
weakly-$*$ in }L^{\infty}(0,T_{0};H^{1}(\Omega_{\varepsilon})),\\\
\rho_{n}\rightarrow\rho\text{ strongly in
}L^{2}((0,T_{0})\times\Omega_{\varepsilon}),\quad
u_{n}\rightharpoonup^{*}u\text{ weakly-$*$ in
}L^{\infty}(0,T_{0};H^{2}(\Omega_{\varepsilon})),\\\ u_{n}\rightharpoonup
u\text{ weakly in
}L^{2}(0,T_{0};H^{3}(\Omega_{\varepsilon})),\quad\partial_{t}u_{n}\rightharpoonup\partial_{t}u\text{
weakly in }L^{2}(0,T_{0};H^{1}(\Omega_{\varepsilon})),\quad\mbox{and}\\\
u_{n}\rightarrow u\text{ strongly in
}L^{2}(0,T_{0};H^{1}(\Omega_{\varepsilon}))\cap
L^{\infty}(0,T_{0};L^{2}(\Omega_{\varepsilon})).\end{gathered}$

Thus, the limit pair $(\rho,u)$ belongs to
$\mathcal{X}_{2}(T_{0})\times\mathcal{Y}_{2}(T_{0})$, concluding the proof of
the existence part.

##### Step IV

We finally adapt the Cauchy estimates in (Step III) to obtain the uniqueness
of the constructed solutions.

##  Appendix B Proof of Proposition 4.1: approximation of the initial data

In this part, we provide the details of proof for Proposition 4.1. We begin
with some quantitative estimates on $R_{\varepsilon}$ defined as in (4.9) that
will be crucially used for the convergence estimates on approximations of
initial data:

(B.1)
$\varepsilon^{\frac{3}{4(\gamma-1)}}R_{\varepsilon}^{\widetilde{\kappa}}=\left(\varepsilon^{\frac{3}{4(\gamma-1)\widetilde{\kappa}}}R_{\varepsilon}\right)^{\widetilde{\kappa}}\leq
c\left(\varepsilon|\log\varepsilon|\right)^{\frac{3}{4(\gamma-1)}}\rightarrow
0\quad\mbox{as}\quad\varepsilon\to 0$

and

(B.2)
$\varepsilon^{2\beta}\left(\varepsilon\exp{(-R_{\varepsilon}^{2})}\right)^{\frac{2(\alpha-\gamma)}{{\gamma-\frac{1}{2}}}}\leq\varepsilon^{2}\exp{(|\log\varepsilon|^{2\theta})}\leq\varepsilon,$

where $\beta$ is appeared in (4.1b). Also, we have (5.9), and more precisely,
$\varepsilon^{\beta}R_{\varepsilon}^{1-\alpha}\to 0$ as
$\varepsilon\rightarrow 0$, which is an important hypothesis in Lemma 6.2.

#### Approximation of the initial density

Here we provide several convergences related to the initial density
$\rho^{\varepsilon}_{0}$ in the lemma below.

######  Lemma B.1.

Let $R_{\varepsilon}$ be given as in (4.9). Then, we have the following
estimates and convergences:

(B.3)
$\displaystyle\inf_{x\in\Omega_{\varepsilon}}\rho_{0}^{\varepsilon}(x)\geq
c\left(\varepsilon\exp{(-R_{\varepsilon}^{2})}\right)^{\frac{1}{{\gamma-\frac{1}{2}}}}>0$
(B.4) $\displaystyle\rho_{0}^{\varepsilon}\rightarrow\rho_{0}\text{ strongly
in }L^{1}\cap L^{\gamma}(\mathbb{R}),$ (B.5)
$\displaystyle|x|^{\kappa+2}\rho_{0}^{\varepsilon}\rightarrow|x|^{\kappa+2}\rho_{0}\text{
strongly in }L^{1}(\mathbb{R}).$

###### Proof.

By definition of $\rho_{0}^{\varepsilon}$, (B.3) is clear.

We observe that $\rho_{0}\in L^{p}(\mathbb{R})\subset L^{1}\cap
L^{\infty}(\mathbb{R})$ for all $p\in[1,\infty]$, which together with (4.4),
gives $\rho_{0}^{\gamma-\frac{1}{2}}\in H^{1}(\mathbb{R})$. Thus, from the
definition of $\widetilde{\rho}_{0}^{\varepsilon}$, we deduce

$\displaystyle\begin{aligned}
&(\widetilde{\rho}_{0}^{\varepsilon})^{\gamma-\frac{1}{2}}\rightarrow\rho_{0}^{\gamma-\frac{1}{2}}\text{
strongly in }L^{p}(\mathbb{R})\text{ for all }p\geq 1\text{ and
}\cr&\partial_{x}((\widetilde{\rho}_{0}^{\varepsilon})^{\gamma-\frac{1}{2}})\rightarrow\partial_{x}(\rho_{0}^{\gamma-\frac{1}{2}})\text{
strongly in }L^{2}(\mathbb{R}).\end{aligned}$

Now we consider two cases: (${\gamma-\frac{1}{2}}\leq 1$) and
(${\gamma-\frac{1}{2}}>1$).

##### (Case I: ${\gamma-\frac{1}{2}}\leq 1$)

Note that

(B.6) $\displaystyle|a^{p}-b^{p}|\leq 2^{p}|a-b||a^{p-1}+b^{p-1}|\text{ for
}a,b>0\text{ and }p\geq 1.$

This together with noticing $\frac{1}{\gamma-\frac{1}{2}}\geq 1$ gives

$\displaystyle|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|$
$\displaystyle\quad=\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)^{\frac{1}{{\gamma-\frac{1}{2}}}}-\rho_{0}\right|$
$\displaystyle\quad\leq
2^{\frac{1}{\gamma-\frac{1}{2}}}\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)-\rho_{0}^{{\gamma-\frac{1}{2}}}\right|\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)^{\frac{1}{{\gamma-\frac{1}{2}}}-1}+\rho_{0}^{\frac{1}{{\gamma-\frac{1}{2}}}-1}\right|.$

Thus, we get for any $p\in[1,\infty]$,

$\|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}\|_{L^{p}}\leq
c\|\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}-\rho_{0}^{{\gamma-\frac{1}{2}}}\|_{L^{p}}+c\varepsilon\to
0\qquad\text{as $\varepsilon\to 0$}.$

In particular, we obtain

$\widetilde{\rho}_{0}^{\varepsilon}\rightarrow\rho_{0}\text{ in }L^{1}\cap
L^{\gamma}(\mathbb{R}),$

and subsequently

(B.7) $Z_{\varepsilon}\to 1\quad\mbox{as}\quad\varepsilon\to 0.$

Then, by using the above convergences, we conclude that

$\rho_{0}^{\varepsilon}\rightarrow\rho_{0}\text{ in }L^{1}\cap
L^{\gamma}(\mathbb{R})\quad\mbox{as}\quad\varepsilon\to 0.$

Next, we show the convergence (B.5). We first observe that

(B.8)
$\displaystyle\left(|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\rho_{0}^{{\gamma-\frac{1}{2}}}\right)\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}\rightarrow|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\rho_{0}^{{\gamma-\frac{1}{2}}}\text{
strongly in
}L^{\frac{1}{{\gamma-\frac{1}{2}}}}(\mathbb{R})\quad\mbox{as}\quad\varepsilon\to
0$

due to $\displaystyle|x|^{\kappa+2}\rho_{0}\in L^{1}(\mathbb{R})$. We then
claim that

(B.9) $\lim_{\varepsilon\to
0}\displaystyle\left\|\left(|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\rho_{0}^{{\gamma-\frac{1}{2}}}\right)\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}-|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}\right)\right\|_{L^{\frac{1}{{\gamma-\frac{1}{2}}}}}=0.$

Note that if (B.9) holds, it follows from (B.8) that

$\displaystyle|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}\right)\rightarrow|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\rho_{0}^{{\gamma-\frac{1}{2}}}\mbox{
strongly in
}L^{\frac{1}{{\gamma-\frac{1}{2}}}}(\mathbb{R})\quad\mbox{as}\quad\varepsilon\to
0.$

Furthermore, by using (B.6) and the fact that

$\displaystyle\left\||x|^{\kappa+2}\left(\varepsilon
e^{-x^{2}}\right)^{\frac{1}{{\gamma-\frac{1}{2}}}}\right\|_{L^{1}}\leq
c\varepsilon^{\frac{1}{\gamma-\frac{1}{2}}},$

we find

$\displaystyle|x|^{\kappa+2}\widetilde{\rho}_{0}^{\varepsilon}\rightarrow|x|^{\kappa+2}\rho_{0}\text{
in }L^{1}(\mathbb{R}).$

Then again we use (B.7) and monotone convergence theorem to have (B.5).

###### Proof of Claim (B.9).

Note that

$||x-y|^{p}-|x|^{p}|\leq c\left(|y|^{p}+|y||x-y|^{p-1}\right)\quad\mbox{for
$p>1$}.$

Since $(\kappa+2)(\gamma-\frac{1}{2})>1$ and
$\|\eta_{\varepsilon^{\gamma-\frac{1}{2}}}\|_{L^{1}}=1$, we estimate

$\displaystyle\int_{\mathbb{R}}{\left|\left(|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\rho_{0}^{{\gamma-\frac{1}{2}}}\right)\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}-|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}\right)\right|^{\frac{1}{{\gamma-\frac{1}{2}}}}}\,{\rm{d}}x$
$\displaystyle\quad\leq\int_{\mathbb{R}}{\left(\int_{\mathbb{R}}{\left|\left(|x-y|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}-|x|^{{(\kappa+2)}{({\gamma-\frac{1}{2}})}}\right)\rho_{0}^{{\gamma-\frac{1}{2}}}(x-y)\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}(y)\;\right|}\,{\rm{d}}y\right)^{\frac{1}{{\gamma-\frac{1}{2}}}}}\,{\rm{d}}x$
$\displaystyle\quad\leq
c\int_{\mathbb{R}}\left(\int_{\mathbb{R}}\left(|y|^{(\kappa+2)(\gamma-\frac{1}{2})}+|y||x-y|^{(\kappa+2)(\gamma-\frac{1}{2})-1}\right)\rho^{\gamma-\frac{1}{2}}_{0}(x-y)\eta_{\varepsilon^{\gamma-\frac{1}{2}}}(y)\,\mathrm{d}y\right)^{\frac{1}{\gamma-\frac{1}{2}}}\,\mathrm{d}x$
$\displaystyle\quad\leq
c\iint_{\mathbb{R}\times\mathbb{R}}\left(|y|^{(\kappa+2)}+|y|^{\frac{1}{\gamma-\frac{1}{2}}}|x-y|^{(\kappa+2)-\frac{1}{\gamma-\frac{1}{2}}}\right)\rho_{0}(x-y)\eta_{\varepsilon^{\gamma-\frac{1}{2}}}(y)\,\mathrm{d}y\,\mathrm{d}x$
$\displaystyle\quad=:I+II,$

due to Jensen’s inequality. Here,

$\displaystyle|I|\leq
c\varepsilon^{(\kappa+2)(\gamma-\frac{1}{2})}\iint_{\mathbb{R}\times\mathbb{R}}\rho_{0}(x-y)\eta_{\varepsilon^{\gamma-\frac{1}{2}}}(y)\,\mathrm{d}y\,\mathrm{d}x=c\varepsilon^{(\kappa+2)(\gamma-\frac{1}{2})}.$

For $II$, we split it into two terms:

$\displaystyle|II|$ $\displaystyle\leq
c\varepsilon\iint_{\mathbb{R}\times\mathbb{R}}|x-y|^{(\kappa+2)-\frac{1}{\gamma-\frac{1}{2}}}({\bf
1}_{|x-y|\leq 1}+{\bf 1}_{|x-y|\geq
1})\rho_{0}(x-y)\eta_{\varepsilon^{\gamma-\frac{1}{2}}}(y)\,\mathrm{d}y\,\mathrm{d}x=:II_{1}+II_{2}.$

Note that $(\kappa+2)(\gamma-\frac{1}{2})-1>0$ and thus

$II_{1}\leq
c\varepsilon\iint_{\mathbb{R}\times\mathbb{R}}\rho_{0}(x-y)\eta_{\varepsilon^{\gamma-\frac{1}{2}}}(y)\,\mathrm{d}y\,\mathrm{d}x=c\varepsilon$

and

$II_{2}\leq
c\varepsilon\iint_{\mathbb{R}\times\mathbb{R}}|x-y|^{\kappa+1}\rho_{0}(x-y)\eta_{\varepsilon^{\gamma-\frac{1}{2}}}(y)\,\mathrm{d}y\,\mathrm{d}x=c\varepsilon\|\rho_{0}\|_{L^{1}_{\kappa+1}}.$

Hence, by combining all of the above estimates, we conclude (B.9). ∎

##### (Case II: ${\gamma-\frac{1}{2}}>1$)

Since

(B.10)
$\displaystyle|a^{\frac{1}{p}}-b^{{\frac{1}{p}}}|\leq|a-b|^{\frac{1}{p}}\text{
for }a,b>0\text{ and }p\geq 1,$

we obtain

(B.11)
$|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|=\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)^{\frac{1}{{\gamma-\frac{1}{2}}}}-\rho_{0}\right|\leq\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)-\rho_{0}^{{\gamma-\frac{1}{2}}}\right|^{{\frac{1}{{\gamma-\frac{1}{2}}}}}.$

This together with $\widetilde{\rho}^{\varepsilon}_{0},\,\rho_{0}\in
L^{\infty}(\mathbb{R})$ implies that for any $p\geq\gamma-\frac{1}{2}$

$\|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}\|_{L^{p}}\rightarrow
0\quad\text{as}\quad\varepsilon\rightarrow 0.$

In particular, we get

$\widetilde{\rho}_{0}^{\varepsilon}\rightarrow\rho_{0}\text{ in
}L^{\gamma}(\mathbb{R}).$

Next, we claim that

$\displaystyle\lim_{\varepsilon\to
0}\|(\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0})\mathbf{1}_{{\overline{\Omega}_{\varepsilon}}}\|_{L^{1}}=0.$

For this, we use (B.11) to estimate

$\displaystyle\int_{\Omega_{\varepsilon}}{|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|}\,{\rm{d}}x$
$\displaystyle\leq\int_{\Omega_{\varepsilon}}{\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)-\rho_{0}^{{\gamma-\frac{1}{2}}}\right|^{{\frac{1}{{\gamma-\frac{1}{2}}}}}}\,{\rm{d}}x$
$\displaystyle\leq\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}\int_{\Omega_{\varepsilon}}{1}\,{\rm{d}}x+\frac{c}{\delta^{2({\gamma-\frac{1}{2}})}}\int_{\Omega_{\varepsilon}}{\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)-\rho_{0}^{{\gamma-\frac{1}{2}}}\right|^{2}}\,{\rm{d}}x$
$\displaystyle\leq
2\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}R_{\varepsilon}+\frac{\varepsilon^{2}}{\delta^{2({\gamma-\frac{1}{2}})}}+\frac{c}{\delta^{2({\gamma-\frac{1}{2}})}}\|\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}-\rho_{0}^{{\gamma-\frac{1}{2}}}\|_{L^{2}}^{2},$

where $\delta>0$ will be determined later. On the other hand, we observe

$\displaystyle\|\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}-\rho_{0}^{{\gamma-\frac{1}{2}}}\|_{L^{2}}\leq
C\varepsilon^{{\gamma-\frac{1}{2}}}\|\rho_{0}^{{\gamma-\frac{1}{2}}}\|_{H^{1}},$

and thus

$\displaystyle\int_{\Omega_{\varepsilon}}{|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|}\,{\rm{d}}x$
$\displaystyle\leq
2\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}R_{\varepsilon}+\frac{\varepsilon^{2}}{\delta^{2({\gamma-\frac{1}{2}})}}+\frac{c}{\delta^{2({\gamma-\frac{1}{2}})}}\varepsilon^{2({\gamma-\frac{1}{2}})}\leq
2\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}R_{\varepsilon}+\frac{c\varepsilon^{2}}{\delta^{2({\gamma-\frac{1}{2}})}}$

due to $\varepsilon<1$. Then, we choose
$\delta=\varepsilon^{\frac{3}{4({\gamma-\frac{1}{2}})}}>0$ to get from (B.1)
that

$\displaystyle\int_{\Omega_{\varepsilon}}{|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|}\,{\rm{d}}x$
$\displaystyle\leq
2\varepsilon^{\frac{3}{4(\gamma-1)}}R_{\varepsilon}+c\varepsilon^{1/2}\leq
2\varepsilon^{\frac{3}{4(\gamma-1)}}R_{\varepsilon}^{\tilde{\kappa}}+c\varepsilon^{1/2}\to
0\quad\mbox{as}\quad\varepsilon\to 0.$

On the other hand, due to the monotone convergence theorem, we find

$\displaystyle\rho_{0}\mathbf{1}_{\overline{\Omega}_{\varepsilon}}\rightarrow\rho_{0}\text{
in }L^{1}(\mathbb{R}).$

Hence, the triangle inequality gives
$\widetilde{\rho}_{0}^{\varepsilon}\mathbf{1}_{\overline{\Omega}_{\varepsilon}}\rightarrow\rho_{0}$
in $L^{1}(\mathbb{R})$. In particular, this shows $Z_{\varepsilon}\to 1$ and
subsequently, we have $\rho^{\varepsilon}_{0}\to\rho_{0}$ in $L^{1}\cap
L^{\gamma}(\mathbb{R})$ as $\varepsilon\to 0$.

For the convergence (B.5), similarly as before, by monotone convergence
theorem, we first observe

$\displaystyle|x|^{2+\kappa}\rho_{0}\mathbf{1}_{{\overline{\Omega}_{\varepsilon}}}\rightarrow|x|^{2+\kappa}\rho_{0}\text{
in }L^{1}(\mathbb{R})\quad\mbox{as}\quad\varepsilon\to 0.$

Then, for any $\delta>0$, we again use (B.11) to estimate

$\displaystyle\int_{\Omega_{\varepsilon}}{|x|^{2+\kappa}|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|}\,{\rm{d}}x$
$\displaystyle\quad\leq\int_{\Omega_{\varepsilon}}{|x|^{2+\kappa}\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)-\rho_{0}^{{\gamma-\frac{1}{2}}}\right|^{{\frac{1}{{\gamma-\frac{1}{2}}}}}}\,{\rm{d}}x$
$\displaystyle\quad\leq\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}\int_{\Omega_{\varepsilon}}{|x|^{(2+\kappa)\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}}\,{\rm{d}}x+\frac{c}{\delta^{2({\gamma-\frac{1}{2}})}}\int_{\Omega_{\varepsilon}}{\left|\left(\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}+\varepsilon
e^{-x^{2}}\right)-\rho_{0}^{{\gamma-\frac{1}{2}}}\right|^{2}}\,{\rm{d}}x$
$\displaystyle\quad\leq
c\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}R_{\varepsilon}^{(2+\kappa)\frac{{\gamma-\frac{1}{2}}}{\gamma-1}+1}+\frac{\varepsilon^{2}}{\delta^{2({\gamma-\frac{1}{2}})}}+\frac{c}{\delta^{2({\gamma-\frac{1}{2}})}}\|\rho_{0}^{{\gamma-\frac{1}{2}}}\ast\eta_{\varepsilon^{{\gamma-\frac{1}{2}}}}-\rho_{0}^{{\gamma-\frac{1}{2}}}\|_{L^{2}}^{2}$
$\displaystyle\quad\leq
c\delta^{\frac{{\gamma-\frac{1}{2}}}{\gamma-1}}R_{\varepsilon}^{(2+\kappa)\frac{{\gamma-\frac{1}{2}}}{\gamma-1}+1}+\frac{c\varepsilon^{2}}{\delta^{2({\gamma-\frac{1}{2}})}}.$

By choosing $\delta=\varepsilon^{\frac{3}{4({\gamma-\frac{1}{2}})}}$, we have

$\displaystyle\int_{\Omega_{\varepsilon}}{|x|^{2+\kappa}|\widetilde{\rho}_{0}^{\varepsilon}-\rho_{0}|}\,{\rm{d}}x$
$\displaystyle\leq
2\varepsilon^{\frac{3}{4(\gamma-1)}}R_{\varepsilon}^{\widetilde{\kappa}}+c\varepsilon^{1/2}\to
0\quad\mbox{as}\quad\varepsilon\to 0,$

where $\widetilde{\kappa}=(2+\kappa)\frac{{\gamma-\frac{1}{2}}}{\gamma-1}+1$
due to (B.1). Then again by using $Z_{\varepsilon}\to 1$ as $\varepsilon\to
0$, we conclude the convergence (B.5). ∎

#### Approximation of initial velocity

In the following lemma, we show the convergence of kinetic energy and uniform-
in-$\varepsilon$ bound estimate on the velocity moments.

######  Lemma B.2.

The following convergence and bound estimate hold:

$\lim_{\varepsilon\to
0}\int_{\mathbb{R}}\rho^{\varepsilon}_{0}|u^{\varepsilon}_{0}|^{2}\,\mathrm{d}x=\int_{\mathbb{R}}\rho_{0}|u_{0}|^{2}\,\mathrm{d}x\quad\mbox{and}\quad\limsup_{\varepsilon\to
0}\int_{\mathbb{R}}\rho^{\varepsilon}_{0}|u^{\varepsilon}_{0}|^{2+\kappa}\,\mathrm{d}x\leq\int_{\mathbb{R}}\rho_{0}|u_{0}|^{2+\kappa}\,\mathrm{d}x.$

###### Proof.

By using the triangle inequality, Hölder’s inequality, and (B.10), we obtain

$\displaystyle\|w_{0}^{\varepsilon}-w_{0}\|_{L^{2}}$
$\displaystyle\leq\|{(\rho_{0}^{\varepsilon})}^{\frac{1}{2}-\frac{1}{2+\kappa}}\left(\tilde{w}_{0}^{\varepsilon}-\tilde{w}_{0}\right)\|_{L^{2}}+\|({(\rho_{0}^{\varepsilon})}^{\frac{1}{2}-\frac{1}{2+\kappa}}-{(\rho_{0})}^{\frac{1}{2}-\frac{1}{2+\kappa}})\tilde{w}_{0}\|_{L^{2}}$
$\displaystyle\leq\|{(\rho_{0}^{\varepsilon})}^{\frac{1}{p}}\|_{L^{p}}\|\tilde{w}_{0}^{\varepsilon}-\tilde{w}_{0}\|_{L^{2+k}}+\|{(\rho_{0}^{\varepsilon})}^{\frac{1}{p}}-{(\rho_{0})}^{\frac{1}{p}}\|_{L^{p}}\|\tilde{w}_{0}\|_{L^{2+\kappa}}$
$\displaystyle\leq\|\rho^{\varepsilon}_{0}\|_{L^{1}}^{\frac{1}{p}}\|\tilde{w}_{0}^{\varepsilon}-\tilde{w}_{0}\|_{L^{2+k}}+\|\rho^{\varepsilon}_{0}-\rho_{0}\|_{L^{1}}\|\tilde{w}_{0}\|_{L^{2+\kappa}}\
$

where $\frac{1}{p}=\frac{1}{2}-\frac{1}{2+\kappa}$. This together with Lemma
B.1 and (4.6) yields the convergence of initial kinetic energy.

For the uniform bound estimate on the velocity moments, by definition of
$u^{\varepsilon}_{0}$, we easily find

$\int_{\mathbb{R}}\rho^{\varepsilon}_{0}|u^{\varepsilon}_{0}|^{2+\kappa}\,\mathrm{d}x\leq\int_{\mathbb{R}}|\tilde{w}^{\varepsilon}_{0}|^{2+\kappa}\,\mathrm{d}x.$

Thus, we get from (4.6) that

$\limsup_{\varepsilon\to
0}\int_{\mathbb{R}}\rho^{\varepsilon}_{0}|u^{\varepsilon}_{0}|^{2+\kappa}\,\mathrm{d}x\leq\int_{\mathbb{R}}\rho_{0}|u_{0}|^{2+\kappa}\,\mathrm{d}x.$

This completes the proof. ∎

We are now in a position to give the details of proof for Proposition 4.1.

###### Proof of Proposition 4.1.

As a direct consequence of Lemma B.2, we get

$\limsup_{\varepsilon\to
0}\int_{\mathbb{R}}\rho^{\varepsilon}_{0}|u^{\varepsilon}_{0}|^{2+\kappa}\,\mathrm{d}x\leq\int_{\mathbb{R}}\rho_{0}|u_{0}|^{2+\kappa}\,\mathrm{d}x.$

For the convergence of $\mathcal{J}_{\tau,\lambda,\varepsilon}$, by definition
of that, it suffices to show that

(B.12)
$\int_{\mathbb{R}}\rho^{\varepsilon}_{0}|u^{\varepsilon}_{0}|^{2}\,\mathrm{d}x\to\int_{\mathbb{R}}\rho_{0}|u_{0}|^{2}\,\mathrm{d}x,\quad\mathcal{F}_{\varepsilon}(\rho^{\varepsilon}_{0})\to\mathcal{F}(\rho_{0}),$

and

(B.13)
$\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|u^{\varepsilon}_{0}+\partial_{x}\delta\mathcal{F}_{\varepsilon}(\rho_{0}^{\varepsilon})|^{2}}\,{\rm{d}}x\to\int_{\Omega_{\varepsilon}}{\rho_{0}|u_{0}+\partial_{x}\delta\mathcal{F}(\rho_{0})|^{2}}\,{\rm{d}}x$

as $\varepsilon\to 0$. Note that the convergences (B.12) can be directly
obtained from the results of Lemmas B.1 and B.2. For (B.13), we notice that

$\displaystyle\frac{1}{2}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|u^{\varepsilon}_{0}+\partial_{x}\delta\mathcal{F}_{\varepsilon}(\rho_{0}^{\varepsilon})|^{2}}\,{\rm{d}}x$
$\displaystyle\quad=\frac{1}{2}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|u^{\varepsilon}_{0}|^{2}}\,{\rm{d}}x+\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}u^{\varepsilon}_{0}\left(\partial_{x}\delta\mathcal{F}_{\varepsilon}(\rho_{0}^{\varepsilon})\right)}\,{\rm{d}}x+\frac{1}{2}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|\partial_{x}\delta\mathcal{F}_{\varepsilon}(\rho_{0}^{\varepsilon})|^{2}}\,{\rm{d}}x,$

where

$\displaystyle\frac{1}{2}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|\partial_{x}\delta\mathcal{F}_{\varepsilon}(\rho_{0}^{\varepsilon})|^{2}}\,{\rm{d}}x$
$\displaystyle=\frac{1}{2}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|\partial_{x}\varphi(\rho_{0}^{\varepsilon})|^{2}}\,{\rm{d}}x+\frac{\varepsilon^{\beta}\alpha}{\alpha-1}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}(\partial_{x}\varphi(\rho_{0}^{\varepsilon}))(\partial_{x}(\rho_{0}^{\varepsilon})^{\alpha-1})}\,{\rm{d}}x$
$\displaystyle\quad+\varepsilon^{2\beta}\frac{1}{2}\left(\frac{\alpha}{\alpha-1}\right)^{2}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|\partial_{x}(\rho_{0}^{\varepsilon})^{\alpha-1}|^{2}}\,{\rm{d}}x.$

Here we use (B.3) to estimate the third term on the right-hand side of the
above equality as

$\displaystyle\varepsilon^{2\beta}\int_{\Omega_{\varepsilon}}{\rho_{0}^{\varepsilon}|\partial_{x}(\rho_{0}^{\varepsilon})^{\alpha-1}|^{2}}\,{\rm{d}}x$
$\displaystyle=c_{\alpha}Z_{\varepsilon}^{1-2\gamma}\varepsilon^{2\beta}\int_{\Omega_{\varepsilon}}{(\rho_{0}^{\varepsilon})^{2(\alpha-\gamma)}|\partial_{x}(\widetilde{\rho}_{0}^{\varepsilon})^{\gamma-\frac{1}{2}}|^{2}}\,{\rm{d}}x$
$\displaystyle\leq
c_{\alpha}\varepsilon^{2\beta}\left(\varepsilon\exp{(-R_{\varepsilon}^{2})}\right)^{\frac{2(\alpha-\gamma)}{{\gamma-\frac{1}{2}}}}\to
0\quad\mbox{as}\quad\varepsilon\to 0,$

due to (B.2) and $Z_{\varepsilon}\to 1$ as $\varepsilon\to 0$. This together
with the convergence results of Lemmas B.1 and B.2 concludes the desired
result. ∎

## Acknowledgement

The research of YPC was supported by NRF grant no. 2022R1A2C1002820. The
research of NC and EZ was supported by the EPSRC Early Career Fellowship
EP/V000586/1. Also, the work of NC was partly supported by the “Excellence
Initiative Research University (IDUB)” program at the University of Warsaw. OT
acknowledges support from NWO grant OCENW.M.21.012.

## References

  * [1] L. Ambrosio, N. Fusco and D. Pallara.  Functions of bounded variation and free discontinuity problems. Oxford university press, 2000.
  * [2] M. Bhatnagar, H. Liu, and C. Tan.  Critical thresholds in the Euler-Poisson-alignment system. J. Differential Equations, 375, 82–119, 2023.
  * [3] D. Bresch and B. Desjardins.  On viscous shallow-water equations (Saint-Venant model) and the quasi-geostrophic limit. C. R. Math. Acad. Sci. Paris, 335, 1079–108, 2002.
  * [4] D. Bresch and B. Desjardins.  Some diffusive capillary models of Korteweg type. C. R. Math. Acad. Sci. Paris, Section Mécanique, 332 , 881–886, 2004.
  * [5] D. Bresch and B. Desjardins.  On the existence of global weak solutions to the Navier-Stokes equations for viscous compressible and heat conducting fluids. J. Math. Pures Appl., 87, 57–90, 2007.
  * [6] J. A. Carrillo and Y.-P. Choi.  Mean-field limits: from particle descriptions to macroscopic equations. Arch. Ration. Mech. Anal., 241, 1529–1573, 2021.
  * [7] J. A. Carrillo, Y.-P. Choi, and J. Jung.  Quantifying the hydrodynamic limit of Vlasov-type equations with alignment and nonlocal forces. Math. Models Methods Appl. Sci., 31, 327–408, 2021.
  * [8] J. A. Carrillo, Y.-P. Choi, and S. P. Perez.  A review on attractive-repulsive hydrodynamics for consensus in collective behavior. In _Active particles. Vol. 1._ , Model. Simul. Sci. Eng. Technol., 259–298. Birkhäuser/Springer, Cham, 2017.
  * [9] J. A. Carrillo, Y.-P. Choi, E. Tadmor, and C. Tan.  Critical thresholds in 1D Euler equations with nonlocal forces. _Math. Mod. Meth. Appl. Sci._ , 26, 185–206, 2016.
  * [10] J. A. Carrillo, Y.-P. Choi, and O. Tse.  Convergence to equilibrium in Wasserstein distance for damped Euler equations with interaction forces. Comm. Math. Phys., 365, 329–361, 2019.
  * [11] J. A. Carrillo, Y.-P. Choi and E. Zatorska.  On the pressureless damped Euler-Poisson equations with quadratic confinement: critical thresholds and large-time behaviour. _Math. Models Methods Appl. Sci._ , 26, 2311–2340, 2016.
  * [12] J. A. Carrillo and S. T. Galtung.  Equivalence of entropy solutions and gradient flows for pressureless 1D Euler systems. arXiv:2312.04932, 2023.
  * [13] J. A. Carrillo, A. Wróblewska-Kami$\acute{\mathrm{n}}$ska, and E. Zatorska. On long-time asymptotics for viscous hydrodynamic models of collective behavior with damping and nonlocal interactions. _Math. Models Methods Appl. Sci._ , 29, 31-63, 2019.
  * [14] J. A. Carrillo, A. Wróblewska-Kami$\acute{\mathrm{n}}$ska, and E. Zatorska.  Pressureless Euler with nonlocal interactions as a singular limit of degenerate Navier-Stokes system. J. Math. Anal. Appl., 492, 124400, 27, 2020.
  * [15] J. A. Carrillo and R. Shu.  Existence of radial global smooth solutions to the pressureless Euler-Poisson equations with quadratic confinement. Arch. Ration. Mech. Anal., 247, 2023.
  * [16] J. A. Carrillo, G.-Q. Chen, D. Yuan, E. Zatorska.  Global solutions of the one-dimensional compressible Euler equations with non-local interactions via the inviscid limit. _In preparation_ 2023\.
  * [17] G.-Q. Chen and M. Perepelitsa.  Vanishing viscosity limit of the Navier-Stokes equations to the Euler equations for compressible fluid flow. Comm. Pure Appl. Math., 63, 1469–1504, 2010.
  * [18] G.-Q. Chen and M. Perepelitsa.  Vanishing viscosity limit solutions of the compressible Euler equations with spherically symmetry and large initial data. Comm. Math. Phys., 338, 771–800, 2015.
  * [19] G.-Q. Chen, L. He, Y. Wang, D. Yuan.  Global solutions of the compressible Euler-Poisson equations with large initial data of spherical symmetry. Comm. Pure Appl. Math., https://doi.org/10.1002/cpa.22149, 2023.
  * [20] Y.-P. Choi, S.-Y. Ha, and Z. Li.  Emergent dynamics of the Cucker-Smale flocking model and its variants. In _Active particles. Vol. 1._ , Model. Simul. Sci. Eng. Technol., 299–331, Birkhäuser/Springer, Cham, 2017.
  * [21] Y.-P. Choi and B.-H. Hwang.  From BGK-alignment model to the pressured Euler-alignment system with singular communication weights. J. Differential Equations, 379, 363–412, 2024.
  * [22] Y.-P. Choi and J. Jung.  Global well-posedness for the Euler-alignment system with singular communication weights in multi-dimensions. Nonlinear Anal. Real World Appl., 76, 104028, 2024.
  * [23] Y.-P. Choi and J. Kim.  Rigorous derivation of the Euler-alignment model with singular communication weights from a kinetic Fokker-Planck-alignment model. Math. Models Methods Appl. Sci., 33, 31–65, 2023.
  * [24] P. Constantin, T. D. Drivas, H. Q. Nguyen, and F. Pasqualotto.  Compressible fluids and active potentials, Ann. Inst. H. Poincaré Anal. Non Linéaire, 37, 145–180, 2020.
  * [25] P. Constantin, T. D. Drivas, and R. Shvydkoy.  Entropy Hierarchies for equations of compressible fluids and self-organized dynamics. SIAM J. Math. Anal., 52, 3073–3092, 2020.
  * [26] F. Cucker and S. Smale.  On the mathematics of emergence. _Japan. J. Math._ , 2, 197–227, 2007.
  * [27] F. Cucker and S. Smale.  Emergent behavior in flocks. _IEEE Trans. Automat. Control_ , 52, 852–862, 2007.
  * [28] T. Do, A. Kiselev, L. Ryzhik, and C. Tan.  Global regularity for the fractional Euler alignment system. _Arch. Ration. Mech. Anal._ , 228, 1–37, 2018.
  * [29] M.-J. Kang and A. Figalli.  A rigorous derivation from the kinetic Cucker-Smale model to the pressureless Euler system with nonlocal alignment. _A_ nalysis & PDE, 12, 843–866, 2019.
  * [30] S.-Y. Ha and E. Tadmor.  From particle to kinetic and hydrodynamic descriptions of flocking. _Kinet. Relat. Models_ , 1, 415–435, 2008.
  * [31] B. Haspot.  Weak-strong uniqueness for compressible Navier-Stokes system with degenerate viscosity coefficient and vacuum in one dimension. Commun. Math. Sci., 15, 587–591, 2017.
  * [32] B. Haspot and E. Zatorska.  From the highly compressible Navier-Stokes equations to the porous medium equation – rate of convergence.  Discrete Contin. Dyn. Syst., 36, 3107–3123, 2016.
  * [33] Q. Jiu and Z. Xin.  The Cauchy problem for 1D compressible flows with density-dependent viscosity coefficients.  Kinet. Relat. Models, 1, 313–330, 2008.
  * [34] R. B. Kelman and T. J. Rivlin.  Classroom notes: conditions for integrand of an improper integral to be bounded or tend to zero.  Am. Math. Mon., 67, 1019–1022, 1960.
  * [35] A. Klar and O. Tse.  An entropy functional and explicit decay rates for a nonlinear partially dissipative hyperbolic system.  ZAMM, 95(5), 469–475, 2015.
  * [36] T. M. Leslie and C. Tan.  Sticky particle Cucker-Smale dynamics and the entropic selection principle for the 1D Euler-alignment system. _Comm. Partial Differential Equations_ , 48, 753–791, 2023.
  * [37] H.-L. Li, J. Li, and Z. Xin.  Vanishing of vacuum states and blow-up phenomena of the compressible Navier-Stokes equations.  Comm. Math. Phys., 281, 401–444, 2008.
  * [38] J. Li and Z. Xin. Global existence of weak solutions to the barotropic compressible Navier-Stokes flows with degenerate viscosities. _arXiv:1504.06826_ , 2015.
  * [39] A. Mellet and A. Vasseur.  On the barotropic compressible Navier-Stokes equations.  Comm. PDEs, 32, 431–452, 2007.
  * [40] P. Minakowski, P. B. Mucha, J. Peszek, and E. Zatorska.  Singular Cucker-Smale dynamics. In _Active Particles, Vol. 2_ , Springer, 2019.
  * [41] R. Shvydkoy.  Dynamics and analysis of alignment models of collective behavior. Nećas Center Series, Birkhäuser/Springer, Cham, 2021.
  * [42] R. Shvydkoy.  Environmental averaging. _a_ rXiv:2211.00117, 2022.
  * [43] R. Shvydkoy and E. Tadmor.  Eulerian dynamics with a commutator forcing II: flocking. _Disc. and Cont. Dyn. Sys._ , 37, 5503–5520, 2017.
  * [44] R. Shvydkoy and E. Tadmor.  Eulerian dynamics with a commutator forcing. _Transactions of Mathematics and Its Applications_ , 1, tnx001, 2017.
  * [45] R. Shvydkoy and E. Tadmor.  Eulerian dynamics with a commutator forcing III. fractional diffusion of order $0<\alpha<1$. _Physica D: Nonlinear Phenomena_ , Special Issue: Nonlinear Partial Differential Equations in Mathematical Fluid Dynamics. 376-377, 131–137, 2018.
  * [46] S. Serfaty.  Mean field limit for Coulomb-type flows. Duke Math. J., 169, 2887–2935, 2020.
  * [47] E. Tadmor and C. Tan.  Critical thresholds in flocking hydrodynamics with non-local alignment. _Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci._ , 372, 20130401, 22, 2014.
  * [48] A. Vasseur and C. Yu.  Existence of global weak solutions for 3D degenerate compressible Navier-Stokes equations. Invent. Math. 206, 935–974, 2016.


# Possibility and Impossibility of the Entropy Balance in Lattice Boltzmann
Collisions

Alexander N. Gorban and Dave Packwood  Department of Mathematics, University
of Leicester, United Kingdom

###### Abstract

We demonstrate that in the space of distributions operated on by lattice
Boltzmann methods that there exists a vicinity of the equilibrium where
collisions with entropy balance are possible and, at the same time, there
exist an area of nonequilibrium distributions where such collisions are
impossible. We calculate and graphically represent these areas for some simple
entropic equilibria using single relaxation time models. Therefore it is shown
that the definition of an entropic LBM is incomplete without a strategy to
deal with certain highly nonequilibrium states. Such strategies should be
explicitly stated as they may result in the production of additional entropy.

##  I Introduction

Lattice Boltzmann schemes are a type of discrete algorithm which can be used
to simulate fluid dynamics and more Succi ; Benzi . Although such a method can
be derived as a discretization of the fully continuous Boltzmann equation,
some thermodynamics properties may be lost in this process. The Entropic
lattice Boltzmann method (ELBM) was invented first in 1998 as a tool for the
construction of single relaxation time lattice Boltzman models which respects
a $H$-theorem Htheorem ; SucciRevModPhys . For this purpose, instead of the
mirror image with a local equilibrium as the reflection center, the entropic
involution was proposed, which preserves the entropy value. Later, it was
called the Karlin-Succi involution gorban06 .

Nevertheless, controlling the proper entropy balance remained until recently a
challenging problem for many lattice Boltzmann models Nonexist . Some
discussions of modern ELBM implementations and results were published recently
KarlinSucciComment .

The distribution functions at the centre of lattice Boltzmann methods are
often referred to and understood as particle densities. Of course for such an
interpretation to be meaningful the distribution function should be strictly
positive. Despite this some lattice Boltzmann implementations may, as a
numerical scheme, tolerate negative population values. An ELBM usually
involves an evaluation of a Boltzmann type entropy function, which does not
exist for negative populations, hence such an ELBM cannot ever tolerate a
negative population value. Due to this there are population values for which
an entropic involution cannot be performed. A complete definition of an ELBM
must include a strategy for what to do in such a situation. The choice of such
a strategy should be explicitly given in any definition of an ELBM as it may
have side-effects with modification of dissipation which should be understood
separately from the influence of the proper entropy balance.

In this paper we study the regions in the spaces of distributions
(populations) where collisions with entropy preservation are possible (near
the equilibrium) and where they are impossible (sufficiently far from the
equilibrium) and demonstrate that both such areas always exist apart some
trivial degenerated cases.

##  II Single Relaxation Time LB Schemes

For fluids, LB systems can be derived as a discretization of the Boltzmann
Equation

${\partial}_{t}f+\mathbf{v}\cdot{\partial}_{\mathbf{x}}f=Q(f)$ (1)

where $f\equiv f(\mathbf{x},\mathbf{v},t)$ is a one particle distribution
function over space, velocity space and time and $Q(f)$ represents the
interaction between particles, sometimes called a collision operation. A
particular example of the interaction $Q(f)$ is the Bhatnagar-Gross-Krook
equation

$Q(f)=-\frac{1}{\tau}(f-f^{\mathrm{eq}}).$ (2)

The BGK operation represents a relaxation towards the local equilibrium
$f^{\mathrm{eq}}$ with rate $1/\tau$. The distribution $f^{\mathrm{eq}}$ is
given by the Maxwell Boltzmann distribution,

$f^{\mathrm{eq}}=\frac{\rho}{(2\pi
T)^{D/2}}\exp\left(\frac{-(\mathbf{v}-\mathbf{u})^{2}}{2T}\right).$ (3)

The macroscopic quantities are available as integrals over velocity space of
the distribution function,

$\rho=\int
f\;\mathrm{d}\mathbf{v},\;\rho\mathbf{u}=\int\mathbf{v}f\;\mathrm{d}\mathbf{v},\;\rho\mathbf{u}^{2}+\rho
T=\int\mathbf{v}^{2}f\;\mathrm{d}\mathbf{v}.$

A discrete approximation to these integrals is the first ingredient to
discretize this system. The scalar field of the population function (over
space, vector space and time) becomes a sequence of vector fields (over space)
in time $f_{i}(\mathbf{x},n_{t}\;\epsilon),n_{t}\in\mathbb{Z}$, where the
elements of the vector each correspond with an element of the quadrature.
Explicitly the macroscopic moments are given by,

$\rho=\sum_{i=1}^{n}f_{i},\;\rho\mathbf{u}=\sum_{i=1}^{n}\mathbf{v}_{i}f_{i},\;\rho\mathbf{u}^{2}+\rho
T=\sum_{i=1}^{n}\mathbf{v}_{i}^{2}f_{i}.$

The complete discrete scheme is given by

$f_{i}(\mathbf{x}+\epsilon\mathbf{v}_{i},t+\epsilon)=f_{i}(\mathbf{x},t)+\omega(f^{\mathrm{eq}}_{i}(\mathbf{x},t)-f_{i}(\mathbf{x},t))$
(4)

where $\epsilon$ is the time step. For this system a discrete equilibrium must
be used. The choice of the velocity set
$\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{n}\\}$ and the discrete equilibrium
distribution $f^{\mathrm{eq}}_{i}$ should provide the best approximation of
the transport equations for the moments by the discrete scheme (4).

##  III ELBM

In the continuous case the Maxwellian distribution maximizes entropy, as
measured by the Boltzmann $H$ function, and therefore also has zero entropy
production. In the context of lattice Boltzmann methods a discrete form of the
$H$-theorem has been suggested as a way to introduce thermodynamic control to
the system Htheorem ; Boghosian .

A variation on the LBGK is the ELBGK ELBM . In this family of methods, the
equilibria are defined as the conditional entropy maximizers under given
values of macroscopic variables (entropic equilibria). The entropies have been
constructed in a lattice dependent fashion in LatticeEntropies . A slightly
different notation is used for the lattice Boltzmann algorithm,

$f_{i}(\mathbf{x}+\epsilon\mathbf{v}_{i},t+\epsilon)=f_{i}(\mathbf{x},t)+\alpha\beta(f^{\mathrm{eq}}_{i}(\mathbf{x},t)-f_{i}(\mathbf{x},t)).$
(5)

The single parameter $\omega$ is replaced by a composite parameter
$\alpha\beta$. In this case $\beta$ controls the viscosity and $\alpha$ is
varied to ensure a constant entropy condition according to the discrete
$H$-theorem. With knowledge of the entropy function $S$, $\alpha$ is found as
the non-trivial root of the equation

$S(\mathbf{f})=S(\mathbf{f}+\alpha(\mathbf{f}^{\ast}-\mathbf{f})).$ (6)

The trivial root $\alpha=0$ returns the entropy value of the original
populations. ELBGK then finds the non-trivial $\alpha$ such that (6) holds.
This version of the BGK collision one calls entropic BGK (or EBGK) collision.
A solution of (6) must be found at every time step and lattice site. The EBGK
collision obviously respects the Second Law (if $\beta\leq 1$), and simple
analysis of entropy dissipation gives the proper evaluation of viscosity.

In general the entropy function is based upon the lattice. For example, in the
case of the simple one dimensional lattice with velocities
$\mathbf{v}=(-c,0,c)$ and corresponding populations
$\mathbf{f}=(f_{-},f_{0},f_{+})$ an explicit Boltzmann style entropy function
is known LatticeEntropies :

$S(\mathbf{f})=-f_{-}\log(f_{-})-f_{0}\log(f_{0}/4)-f_{+}\log(f_{+}).$ (7)

##  IV Regions of Existence and Non-existence of Entropic Involution

Let us study the entropic involution in the distribution simplex $\Sigma$
given by $\sum f_{i}=const>0$, $f_{i}\geq 0$.

Figure 1: The simplex $\Sigma$ is given by the white background. (1)
Populations relax through the equilibrium given by the single point to an
equal entropy point, if possible. The boundary of this possibility is given.
(2) The regions $A$ (the entropic involution is possible) and $B$ (the
involution is impossible) as subsets of the simplex divided by this boundary
are presented.  Figure 2: The simplex $\Sigma$ is given by the white
background. (1) Populations relax through the their corresponding equilibrium
point along the line given by constant $u$ to an equal entropy point, if
possible. The boundary of this possibility is given. (2) The regions $A$ and
$B$ separated by this boundary are presented.

Let us prove that under very natural assumptions about some properties of the
entropy that the simplex of distributions can be split into two subsets $A$
and $B$: in the set $A$ the entropic involution exists, and for distributions
from the set $B$ equation (6) has no non-trivial solutions. Both sets $A$ and
$B$ have non-empty interior (apart of a trivial symmetric degenerated case).

Let the entropy $S$ be a strictly concave continuous function in the
distribution simplex $\Sigma$. We assume also that $S$ is twice
differentiable, the Hessian of $S$, $\partial^{2}S/\partial f_{i}\partial
f_{j}$, is negative definite in the interior of the simplex, $\Sigma_{+}$,
where $\sum f_{i}=const$, $f_{i}>0$ and the global maximizer of $S$, the
equilibrium, belongs to the interior of the simplex.

For example, the relative Boltzmann entropy, $S=-\sum
f_{i}(\ln(f_{i}/W_{i})-1)$, $W_{i}>0$, satisfies these conditions, because
$f\ln f\to 0$ when $f\to 0$ and $\partial^{2}S/\partial f_{i}\partial
f_{j}=-\delta_{ij}/f_{i}$, whereas the relative Burg entropy $S=\sum
W_{i}(\ln(f_{i}/W_{i}))$ does not satisfy these conditions because it does not
exist on the border of the simplex.

Macroscopic variables are linear functions of $\mathbf{f}$. The sets with
given values of the macroscopic variables in the simplex $\Sigma$ are
polyhedra, intersections of $\Sigma$ with linear manifolds with the given
values of moments. We assume that in any such a polyhedron the entropy
achieves its (conditionally) global maximum at an internal point. This
assumption holds for the Boltzmann relative entropy because of the logarithmic
singularity of the “chemical potentials” $\mu_{i}=\ln(f_{i}/W_{i})$ on the
border of positivity. These maximizers are equilibria. If $\mathbf{f}$ is
sufficiently close to a positive equilibrium then, due to the implicit
function theorem, the nontrivial solution to equation (6) exists and it gives
$\alpha=2+o(\mathbf{f}-\mathbf{f}^{*})$. The value $\alpha=2$ corresponds to
the mirror image, the small term $o(\mathbf{f}-\mathbf{f}^{*})$ gives the
corrections to the value $\alpha=2$. Therefore, in some vicinity of the
equilibrium the entropic involution exists.

To prove the existence of the area where entropic involution is impossible,
let us consider one polyhedron with given values of the macroscopic variables
and a positive equilibrium. The local minima of the entropy in this polyhedron
are situated at the vertices. At least one of them is a global minimum. Let
this vertex be $\mathbf{f}^{\mathbf{v}}$. Let us draw a straight line $l$
through points $\mathbf{f}^{\mathbf{v}}$ and $\mathbf{f}^{*}$. The
intersection $l\cap\Sigma$ is an interval and $S$ achieves its global minimum
on this interval at the point $\mathbf{f}^{\mathbf{v}}$. If the dimension of
the polyhedron is more than one then the opposite end of this interval is not
even a local minimum of $S$ in the polyhedron and the entropic involution does
not exists for $\mathbf{f}^{\mathbf{v}}$ and some vicinity around it.

A special degeneration is possible when the polyhedra are one-dimensional,
i.e. intervals, and the values of the entropy at both ends of each interval
coincide. For example, for two-dimensional distributions, $f_{+},f_{-}$, the
entropy $S==f_{+}\ln f_{+}-f_{-}\ln f_{-}$ and the macroscopic variable
$\rho=f_{+}+f_{-}$. Apart from such symmetric one-dimensional cases there
exists an area near the maximally non-equilibrium vertex
$\mathbf{f}^{\mathbf{v}}$ where the entropic involution cannot be defined.
Such an area may also exist near some other vertices, where local entropy
minima are reached.

For the Burg entropy, the entropic involution is always possible Boghosian
because it tends to $-\infty$ at the border of positivity. The same is true,
for the relative entropy of the form $S=-\beta^{-1}\sum
W_{i}((W_{i}/f_{i})^{\beta}-1)$ that tends to the Burg entropy when $\beta\to
0$ GoGoJudge2011 . This negative brunch of the relative Tsallis entropy is
less known. The standard Tsallis entropy Tsallis  is finite at the border of
positivity, hence, collisions with entropy preservation are not always
possible for it.

We now demonstrate the population function values where the involution cannot
be performed for some simple examples. We use the standard 1-D lattice
described in Section III with the discrete equilibrium given in Eq 7. We begin
with an LBM with only one conserved moment in collision, namely density. The
equilibrium is
$f_{-}^{*}=\frac{\rho}{6},\;\;f_{0}^{*}=\frac{2\rho}{3},\;\;f_{+}^{*}=\frac{\rho}{6}.$

In Fig. 1, the simplex $\Sigma$ of positive populations with a fixed density
$\rho=1$ is the triangle given by the intersection of three half-planes,
$f^{+}>0,\,f^{-}>0$, and $1-f^{+}-f^{-}>0$. Within that region we plot several
entropy level contours $S(\mathbf{f})=c$ and the unique equilibrium point. The
region is divided into the parts where the entropic involution is possible
(around the equilibrium) and where it is impossible.

A more common use of lattice Boltzmann involves a second fixed moment,
momentum. The entropic equilibria used by the ELBGK are available explicitly
as the maximum of the entropy function (7),

$f_{\mp}^{*}=\frac{\rho}{6}(\mp
3u-1+2\sqrt{1\\!+\\!3u^{2}}),\;f_{0}^{*}=\frac{2\rho}{3}(2-\sqrt{1\\!+\\!3u^{2}}).$

In this case the dimension of the equilibrium is one greater. In Fig. 2 all
relaxation occurs parallel to the lines of constant $u$. The region where
entropic involution is possible is again given.

In each experiment the region is discretized into many individual points. For
each point a value for $\alpha$ is attempted to be found. The method used is
simply to begin with a guess of $\alpha=1$ and then add increments of
$10^{-3}$ until a solution of Eq. 6 occurs, or the edge of the positivity
domain is reached. This method would be inappropriate to use in a usual ELBM,
due to the very large computational cost, but it is very robust and hence
useful for this experiment with many higly non-equilibrium distributions.
Another approach (with the same result) implies calculation of the entropic
involution for all the boundary points where it exists. In this method we draw
a straight line $l$ through a boundary point $\mathbf{f}$ and the equilibrium
and find the intersection $l\cap\Sigma$ which consists of all points on $l$
with non-negative coordinates. One end of this interval is $\mathbf{f}$,
another end is also a boundary point, $\mathbf{f}^{\prime}$. The entropic
involution for $\mathbf{f}$ exits if and only if $S(\mathbf{f})^{\prime}\leq
S(\mathbf{f})$. After we check this inequality, we can solve Eq. (6). The
images of these involutions form the border that separates sets $A$ and $B$
(see Figs).

##  V Conclusion

The entropic involution is not always possible to perform. We have
demonstrated that apart some special one-dimensional spaces of distributions
with additional symmetry there exist domains where collisions with the
preservation of entropy are not possible. We illustrated this statement by
some simple and well known examples of ELBGK systems for which we directly
calculated the areas where entropic collisions exist and where they do not
exist.

Such phenomena should be observable in all ELBM schemes with the classical
entropies: there exists a vicinity of the equilibrium where the entropic
involution is possible but for some areas of non-equilibrium distributions
there exists no non-trivial root of equation (6). A collision which preserves
entropy does not exist for this area. Therefore, for the regimes close to
equilibrium (the vicinities $A$ of equilibria, Figs 1,2), ELBM schemes
guaranty the precise balance of the entropy and for more nonequilibrium
regimes, when at some sites the distribution belonges to sets $B$, ELBM
schemes work as limiters Limiters . with additional dissipation. It is
necessary for any complete definition of an ELBM algorithm to prescribe what
to do when the involution is not possible. A reasonable choice would be to
over-relax the maximum amount possible while maintaining positive population
values. Such a technique is independently in use as a stabilizer for lattice
Boltzmann schemes, sometimes called the ‘positivity limiter’ BGJ ; Limiters ;
Li ; Servan ; Tosi . An effect of this operation is a local increase in
viscosity/entropy production. Hence, if an ELBM were to apply such a scheme it
would necessarily break the proper entropy balance. In this sense, ELBM
belongs to a large family of add-ons that regularise LBM by the management of
the addtional dissipation Add-ons .

## References

  * (1) S. Ansumali and I. V. Karlin, Phys. Rev. E, 62 (6), 7999–8003 (2000).
  * (2) R. Benzi, S. Succi, and M. Vergassola, Phys. Reports, 222, 145–197 (1992).
  * (3) B. M. Boghosian, J. Yepez, P. V. Coveney, A. Wagner, Proc. R. Soc. Lond. A, 457, 717–766 (2001).
  * (4) R. A. Brownlee, A. N. Gorban, and J. Levesley, Phys. Rev. E, 75, 036711 (2007).
  * (5) R. A. Brownlee, A. N. Gorban, and J. Levesley, Physica A, 387 (2-3), 385–406 (2008).
  * (6) R. A. Brownlee, J. Levesley, D. Packwood, A.N. Gorban, arXiv:1110.0270 [physics.comp-ph].
  * (7) A. N. Gorban, in _Model Reduction and Coarse-Graining Approaches for Multiscale Phenomena_ , (Springer, Berlin-Heidelberg, New York 2006), 117–176; arXiv:cond-mat/0602024 [cond-mat.stat-mech].
  * (8) A. N. Gorban, P. A. Gorban, G. Judge, arXiv:1003.1377 [physics.data-an].
  * (9) I. V. Karlin, A. Ferrante, and H. C. Öttinger, Europhys. Lett. 47, 182–188 (1999).
  * (10) I. V. Karlin, A. N. Gorban, S. Succi, and V. Boffi,  Phys. Rev. Lett., 81, 6–9 (1998).
  * (11) I. V. Karlin, S. Succi, arXiv:1107.3025 [cond-mat.stat-mech]
  * (12) Y. Li, R. Shock, R. Zhang, H. Chen. J. Fluid Mech., 519, 273–300 (2004)
  * (13) B. Servan-Camas, FT-C. Tsai. J Comput Physics 228 (1), 236–256 (2009).
  * (14) S. Succi, _The lattice Boltzmann equation for fluid dynamics and beyond_ (Oxford University Press, New York 2001).
  * (15) S. Succi, I. V. Karlin, and H. Chen, Rev. Mod. Phys. 74, 1203 (2002).
  * (16) F. Tosi, S. Ubertini, S. Succi, H. Chen, I.V. Karlin, Math Comput Simulation, 72, 227–231 (2006).
  * (17) C. Tsallis, J. Stat. Phys., 52, 479–487 (1988).
  * (18) W.-A. Yong and L.-S. Luo, Phys. Rev. E, 67, 051105 (2003).