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

HAL
open science

Learning Extremely High Density Crowds as Active
Matters

Feixiang He, Jiangbei Yue, Jialin Zhu, Armin Seyfried, Dan Casas, Julien

Pettré, He Wang

► To cite this version:

Feixiang He, Jiangbei Yue, Jialin Zhu, Armin Seyfried, Dan Casas, et al. Learning Extremely High Density Crowds as Active Matters. CVPR 2025 - IEEE/CVF Conference on Computer Vision and Pattern Recognition, Jun 2025, Nashville, United States. pp.540-550, 10.1109/CVPR52734.2025.00059. hal-05244322

HAL Id: hal-05244322

https://inria.hal.science/hal-05244322v1

Submitted on 8 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>
Learning Extremely High Density Crowds as Active Matters

Feixiang He², Jiangbei Yue³, Jialin Zhu², Armin Seyfried⁴, Dan Casas⁵,* Julien Pettré⁶, He Wang ¹,²†

¹AI Centre, University College London, UK      ²University College London, UK

³University of Leeds, UK      ⁴Forschungszentrum Jülich, Germany

⁵Universidad Rey Juan Carlos, Spain      ⁶INRIA Rennes, France

Abstract

Video-based high-density crowd analysis and prediction has been a long-standing topic in computer vision. It is notori- ously difficult due to, but not limited to, the lack of high- quality data and complex crowd dynamics. Consequently, it has been relatively under studied. In this paper, we propose a new approach that aims to learn from in-the-wild videos, often with low quality where it is difficult to track individu- als or count heads. The key novelty is a new physics prior to model crowd dynamics. We model high-density crowds as active matter, a continuum with active particles subject to stochastic forces, named ‘crowd material’. Our physics model is combined with neural networks, resulting in a neu- ral stochastic differential equation system which can mimic the complex crowd dynamics. Due to the lack of similar research, we adapt a range of existing methods which are close to ours for comparison. Through exhaustive evalua- tion, we show our model outperforms existing methods in analyzing and forecasting extremely high-density crowds. Furthermore, since our model is a continuous-time physics model, it can be used for simulation and analysis, provid- ing strong interpretability. This is categorically different from most deep learning methods, which are discrete-time models and black-boxes.

# 1. Introduction

Video-based crowd forecasting/analysis has been a long-
standing problem in computer vision [30], and is espe-
cially crucial in public safety and event organization [19].
In this area, one of the most difficult scenarios is ex-
tremely high-density crowds. We are particularly inter-
ested in crowds where the density is at least 5p/m² (people
per square meter). At this density level, forecast/analysis
becomes notoriously challenging as the dynamics can be-
come violently chaotic [14], resulting in life-threatening

crushes [47], which is under-studied in computer vision.

Early research has primarily focused on relatively low
density, with empirical modeling [5, 21, 40, 52, 55] and
data-driven methods [3, 16, 31, 32, 37, 67]. While the for-
mer is explicit and explainable but inaccurate in prediction,
the latter is accurate in prediction but lacks of explainabil-
ity [56, 65, 66]. A recent trend is hybrid models, combin-
ing black-box neural networks (NNs) with explicit models,
leveraging the advantages of both methodologies. These
methods combine NNs with differential equations [71], or
with physics conditioned on social interactions [9, 68], or
with physics and symbolic regression [39]. However, these
approaches are designed mainly for low-density crowds,
where the microscopic motion is simple and trajectory data
can be reliably obtained, which is unlikely in high-density
crowds from in-the-wild videos. Furthermore, merely fo-
cused on local interactions, they lack the ability to capture
large scale/global flow patterns [57]. A very recent attempt
introduced a hydrodynamics-informed neural network for
simulating dense crowds [76] utilizing optical flows. How-
ever, it fails to account for the self-propelled characteristics
of crowds and is not suitable for agent-based simulation.

Extremely high-density crowds present unique chal-
lenges for prediction and analysis. The foremost is the
scarcity of high-quality data. For high-density crowds,
often the only available data is Closed-Circuit Television
(CCTV) videos, because cameras are the *only* sensors al-
lowed in many public spaces due to their non-invasiveness.
As a result, it is imperative to be able to work directly with
these videos. However, in public spaces, the cameras are
often far from the crowd and the data is noisy therefore un-
likely to provide accurate individual motion, if at all, even
with crowd counting/localization methods [1, 33, 35, 36].
This makes agent-based analysis [18, 58–60, 71, 72],
which relies on trajectory estimation, impractical. Further-
more, complex dynamics arise when crowd density exceeds
5p/m². Vastly different from low-density crowds, high-
density crowds behave like an active matter [41], e.g. a con-
tinuum with particles moving non-deterministically caused
by autonomous behaviors and self-propulsion even when

<footnote>*work done prior to joining Amazon</footnote>

<footnote>†Corresponding author, he_wang@ucl.ac.uk</footnote>
being physically constrained by the surroundings. Consequently, the local interactions become significantly more pronounced, sometimes causing spatiotemporal fluctuations similar to waves [2, 22] leading to life-threatening states in crowds [44]. Accurately modeling such dynamics remains an open problem. In the past, there has been parallel drawn between high-density crowds and continuum [40]. Such parallels lead to explicit models but they are not learnable, i.e. hard to be adapted to specific crowds.

To address the aforementioned challenges, we propose a new learnable active matter model for high-density crowds. Since crowds show distinctive dynamics compared with common homogeneous materials e.g. water, we design a new material model, named “crowd material”, which has a special stress-strain relation to reflect the empirical observations on individuals in high-density crowds. Additionally, to mimic active matters, this new material model is augmented by integrating the collective dynamics of self-propelled agents, based on the Toner-Tu (TT) equation [51], enabling us to model and learn random individual motions. Next, our model is described by a new neural stochastic differential equation system. To directly learn from raw videos, we treat the optical flows from videos as noisy observations of the underlying velocity field of the continuum. Meanwhile, we also model individuals as active 2D particles. Note that treating people as 2D particles in a continuum is a common practice in macroscopic high-density crowd research [28]. This dual continuum-particle perspective naturally leads us to a new Material Point Method (MPM) [27], where the Eulerian grid discretizes the space, and the Lagrangian particles represent the individuals. We refer to this new MPM framework as CrowdMPM. Finally, we utilize NNs to learn the active forces and key material parameters.

Our model is designed to handle optical flows estimated from real-world videos. In addition, due to the scarcity of high-density crowd data and our ambition to learn from in the wild data, we use data from lab environments as well as public spaces. These datasets are all high-density scenarios with diverse dynamics such as evacuation, music festivals, etc. Next, due to the lack of similar research, we adapt both physics-based and data-driven methods for comparison. The experiments show that our model outperforms the alternatives in capturing, analyzing and predicting high-density crowd dynamics. Our contributions include:

* a new video-based analysis and prediction framework for extremely high-density crowds.

* a new differentiable physics model for high-density crowd dynamics.

* a new crowd material model capable of capturing crowds as active matters.

* a new Eulerian-Lagrangian scheme, CrowdMPM.

## 2. Related Work

### 2.1. Behavior Understanding from Videos

Traditional approaches for crowd behavior analysis are individual-based [8, 17, 58-60] and global-based [6, 24, 70]. The former treats a crowd as a collection of individuals and analyzes behaviors by segmenting or detecting individuals; conversely, the latter treats the crowd as a single entity without segmentation. In dense scenes, it is very difficult to employ individual-based approaches as tracking individuals is extremely challenging. Benefited from the data-fitting ability of deep learning (DL), DL-based approaches are proposed to understand human behaviors in crowded scenarios, such as behavior classification [66], anomaly detection [45] or predicting pedestrian trajectories [3, 16]. However, they are not designed for predicting high-density crowd dynamics directly from optical flow and cannot be used for simulation with changed environments and densities.

### 2.2. Differentiable Physics

Solving differentiable equations (DE) with deep learning has recently spiked strong interests. The research can be grouped into different categories according to how deep learning is involved. NNs can be used to generate finite element meshes [73, 74]. Neural differential equation models [43, 71, 72] attempt to parameterize part of the DEs by neural networks. Differentiable physics [15, 34, 64] focus on making the whole simulation process differentiable. Physics-informed neural networks [42, 46, 76] aim to bypass the DE solve using NNs for prediction. Highly inspired by the research above, we propose a novel differentiable physics model for forecasting high-density crowds.

### 2.3. Video Prediction

Recently, video prediction has become a prominent area of research, focusing on forecasting future frames based on past observations, which aligns closely with our objective of predicting future optical flow frames. Most existing video prediction methods like recurrent neural networks [62, 63], convolutional neural networks [12, 75], or attention mechanisms [49] are tailored for RGB inputs rather than optical flow. Shifting the input from RGB frames to optical flow enables us to focus directly on motion dynamics, a critical factor in crowd movement prediction. However, these methods typically require a long input sequence and rely on deep learning architectures with large numbers of parameters, which limits interpretability. In contrast, our model requires only a single optical flow frame as input and maintains interpretability, providing a more efficient and explainable approach to predicting crowd dynamics.
<figure><img src="image_4.png" /><figcaption>Figure 1. Overview. From left to right: optical flow estimation, velocity field generation, initial particle sampling, crowd simulation and loss calculation.</figcaption></figure>

## 3. Methodology

Our whole framework is shown in Fig. 1. Our model takes optical flows as input from videos, simulate the crowd movements, then learn the model parameters based on the difference between the predicted and the observed crowd motions. For simplicity, we only give the key equations below and refer the readers to the supplementary materials (SMs) for details.

### 3.1. Crowd as an Active Matter

High-density crowds bear strong similarity with continuum [52]. We start by modeling the continuum with mass and momentum conservation:

$$ \frac{D\rho}{Dt} + \rho\nabla \cdot \mathbf{v} = 0 \quad \text{and} \quad \rho\frac{D\mathbf{v}}{Dt} = \nabla \cdot \sigma^{\text{cm}} + \rho\mathbf{b} + \mathbf{f}^{\text{act}} \quad (1) $$

where $\rho$ is density and $\frac{D}{Dt}$ is the material derivative [50]. $\mathbf{v}$ is the velocity field. $\rho\mathbf{b}$ denotes the body force (e.g. gravity). The main difference between Eq. (1) and a standard formulation is the two terms, the new crowd material stress $\sigma^{\text{cm}}$, and the new stochastic active force $\mathbf{f}^{\text{act}}$. The new learnable $\sigma^{\text{cm}}$ is key to our “crowd material” property explained in Sec. 3.3, and the new learnable active force $\mathbf{f}^{\text{act}}$ captures the random motions caused by active self-propulsion explained in Sec. 3.4. We keep $\rho\mathbf{b}$ for forces like global attraction, if needed. Overall since we learn the parameters in $\sigma^{\text{cm}}$ and $\mathbf{f}^{\text{act}}$ by neural networks, Eq. (1) is a neural stochastic differential equation system.

### 3.2. CrowdMPM

To solve Eq. (1), Eulerian methods and Lagrangian methods are the two main approaches. Eulerian methods approximate quantities (e.g. velocity) on a spatial grid over a continuous domain, while Lagrangian methods discretize the continuum into mass particles and trace the quantities on them. However, the challenge is we only have Eulerian data but need to model Lagrangian behaviors. The only data available is video, where we can estimate the optical flows. The optical flows can be seen as a noisy observation of the velocity field of the underlying continuum, i.e. Eulerian data. Meanwhile, we need to model active forces on individuals, i.e. Lagrangian behaviors. Eulerian methods are

limited in modeling the individual behaviors explained later which is crucial for active forces, and Lagrangian methods do not work well with our Eulerian data as they can only compute quantities in areas when there are particles which cannot guarantee to cover the whole space.

Therefore, we choose material point methods (MPM) [25]. MPM is a hybrid Eulerian-Lagrangian method, which uses both spatial grids and mass particles. However, in our CrowdMPM, the particles are different from those in standard MPM which only serve as collocation points for integration and not corresponding to any entities. Our particles are individuals. This perspective is seemingly trivial but crucial in modeling the crowd material introduced later in Sec. 3.3 and Sec. 3.4.

Typically, MPM takes three steps to solve an equation: (1) Particle-to-grid transfer (P2G) (2) Grid Operation (GO) and (3) Grid-to-particle transfer (G2P). First, P2G transfers mass, velocity, and pressure from particles to the grid nodes, then GO solves the momentum equation under boundary conditions on the grid nodes, and finally, the updated velocity and deformation are transferred back to particles from the grid nodes. Here, we directly give the final discretized update scheme and leave the details in the SM.

**Notation** A standard MPM discretization [25] is shown in Fig. 2. Formally, we use subscript $i$ to index variables on a 2D grid with cell size $\Delta x$, and $p$ to index variables on particles, and superscript indicating time steps. Node $i$ and particle $p$ have attributes including mass $m_i$, $m_p \in \mathbb{R}$, position $\mathbf{x}_i$, $\mathbf{x}_p \in \mathbb{R}^2$, and velocity $\mathbf{v}_i$, $\mathbf{v}_p \in \mathbb{R}^2$. Each particle has an initial (undeformed) volume $V_p^0$ at $t=0$. $C_p$ is the affine velocity gradient [26] and $\mathbf{P}_p = J_p \sigma^{\text{cm}} \mathbf{F}_p^{-T}$ is the PK1 (the first Piola-Kirchhoff) stress, where $\mathbf{F}_p$ is the deformation gradient and $J_p$ is the determinant of $\mathbf{F}_p$ [27]. $\phi$ is a quadratic B-spline function [25]. We use $\Delta t$ to denote the time interval between time steps.

At time $n$, we compute the P2G step via:

$$ \begin{aligned} m_i^n &= \sum_p w_{ip}^n m_p \\ m_i^n \mathbf{v}_i^n &= \sum_p w_{ip}^n [m_p \mathbf{v}_p^n + m_p C_p^n (\mathbf{x}_i - \mathbf{x}_p^n)], \end{aligned} \quad (2) $$
<figure><img src="image_5.png" /><figcaption>Figure 2. (a)-(d) Standard MPM. (e): particle-to-particle traction modeling.</figcaption></figure>

where $w_{ip}^n = \phi(\mathbf{x}_i - \mathbf{x}_p^n)$. Then, we update velocities in the GO step via:

$$ \mathbf{v}_i^{n+1} = \mathbf{v}_i^n + \Delta t \mathbf{f}_i^n / m_i, \quad (3) $$

where $\mathbf{f}_i^n = \mathbf{f}_i^{n,st} + \mathbf{f}_i^{n,act} + \mathbf{f}_i^{n,bd}$ is specified by our new stress and active forces explained later. Boundary conditions are also considered in GO to refine $\mathbf{v}_i^{n+1}$:

$$ \mathbf{v}_i^{n+1} = \mathbf{B}\mathbf{C}(\mathbf{v}_i^{n+1}) = \mathbf{v}_i^{n+1} - \gamma \mathbf{n} \langle \mathbf{n}, \mathbf{v}_i^{n+1} \rangle, \quad (4) $$

where $\mathbf{n}$ is the surface normal of the boundary, $\gamma$ is a scalar. **BC** is a (dampened) no-slip boundary [10].

Lastly, we update the particle velocity $\mathbf{v}_p^{n+1}$, velocity gradient $\mathbf{C}_p^{n+1}$, deformation gradient $\mathbf{F}_p^{n+1}$ and position $\mathbf{x}_p^{n+1}$ by gathering information from their neighboring grid nodes, i.e. the G2P step:

$$ \mathbf{v}_p^{n+1} = \sum_i w_{ip}^n \mathbf{v}_i^{n+1}, \quad \mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \mathbf{v}_p^{n+1}, \quad (5) $$

$$ \mathbf{C}_p^{n+1} = \frac{4}{\Delta x^2} \sum_i w_{ip}^n \mathbf{v}_i^{n+1} (\mathbf{x}_i - \mathbf{x}_p^n)^T, \quad (6) $$

$$ \mathbf{F}_p^{n+1} = (\mathbf{I} + \Delta t \mathbf{C}_p^{n+1}) \mathbf{F}_p^n. \quad (7) $$

Overall, solving Eq. (2)-Eq. (7) solves Eq. (1) in every step.

### 3.3. New Learnable Stress $\sigma^{cm}$

Empirically, the 'crowd material' is significantly different from homogeneous materials like water. Firstly, crowds can easily scatter but not be compressed severely [40], due to people cannot be superpositioned, i.e. easily stretched but not as easily compressed, which we refer to as elastic asymmetry. Next, during compression, there is a multi-phasic resistance change based on the inter-personal distance [20]. The resistance is small when the distance is larger than a comfortable threshold, but exponentially increased within the threshold until it reaches incompressibility (e.g. people in full contact). We refer to this property as exponential resistance. Finally, while resisting compression, people can be more tolerant to relative motions such as passing each other at a close distance. This suggests that 'crowd material' has less resistance in shearing and rotations, to which we refer as compression dominance.

The aforementioned observations dictate that we need to design a new material model that can capture these features and is also learnable. To this end, we propose a new

learnable stress inspired by the above observations. First, unlike particles in standard MPM as collocation points and not corresponding to any entities, we need to model the strain-stress between particles directly, not via the grid nodes as in standard MPM. So we treat each individual as a Lagrangian particle. This leads to a new stress $\sigma^{cm}$ below for the force $\mathbf{f}_i^{st}$ in Eq. (3).

To model the elastic asymmetry, we choose a Cauchy stress for weakly compressible fluid [48] $\sigma_p = \epsilon(1 - \frac{1}{J_p})\mathbf{I}$, where $\epsilon$ is the Young's modulus, therefore, the force $\mathbf{f}_i^{st}$ in Eq. (3) can be denoted as:

$$ \mathbf{f}_i^{st} = \sum_p \omega_{ip} \mathbf{G}_p (\mathbf{x}_i - \mathbf{x}_p), \quad (8) $$

$$ \mathbf{G}_p = -\frac{4}{\Delta x^2} V_p^0 \sigma_p J_p = -\frac{4}{\Delta x^2} \epsilon V_p^0 (J_p - 1) \mathbf{I}, \quad (9) $$

where we omit the time superscripts $n$ for simplicity. Eq. (8) follows the standard MPM, which ensures the overall small compressibility.

Next, to model the exponential resistance, we augment each particle with a volume, shown in Fig. 1 (Initial Particles), which consists of an incompressible volume (green disc) and a comfortable zone (white disc). We assume that when a particle $p'$ is within the comfortable zone of particle $p$, there is a repulsive force based on their distance. We use a comfort distance $d_c = 2(r_b - r_a)$ where $r_a$ and $r_b$ are the radii of the incompressible volume and the comfortable zone respectively. Then we model the repulsive force as:

$$ \mathbf{f}_r^{pp'} = \begin{cases} -k \log(d_{pp'}) \mathbf{e}_{p'p} + \mathbf{f}_t^{pp'} & \text{if } 0 < d_{pp'} < 1 \\ \mathbf{f}_t^{pp'} & \text{otherwise} \end{cases} \quad (10) $$

where $d_{pp'} = \frac{\|\mathbf{x}_p - \mathbf{x}_{p'}\| - 2r_a}{d_c}$ is the ratio between the current distance and the comfortable distance and $\mathbf{e}_{p'p} = \frac{\mathbf{x}_p - \mathbf{x}_{p'}}{|\mathbf{x}_p - \mathbf{x}_{p'}|}$. $k$ is a positive learnable parameter. Therefore, $-k \log(d_{pp'}) \mathbf{e}_{p'p}$ models the exponential resistance.

To model the compression dominance, we separate the inter-particle compression from other stress components related to shear/rotation. We introduce a term $\mathbf{f}_t^{pp'}$ in Eq. (10). $\mathbf{f}_t^{pp'}$ is the traction at $p$ (due to $\mathbf{G}_p$) projected along the direction from $p'$ to $p$, which can be calculated as $\mathbf{f}_t^{pp'} = \langle \mathbf{G}_p(\mathbf{x}_i - \mathbf{x}_p), \mathbf{e}_{p'p} \rangle \mathbf{e}_{p'p}$, shown in Fig. 2 (e).
Then we incorporate Eq. (9) and Eq. (10) into Eq. (8) and modify $\mathbf{f}_i^{st}$ for grid update in Eq. (3) to:

$$ \mathbf{f}_i^{st} = \sum_p w_{ip} \{ \mathbf{G}_p(\mathbf{x}_i - \mathbf{x}_p) + \sum_{p' \in N_p} (\mathbf{f}_r^{pp'} - \mathbf{f}_t^{pp'}) \}, \quad (11) $$

where $N_p$ is the set of the neighboring particles of $p$. In Eq. (11), we see that when $\mathbf{f}_r^{pp'} = \mathbf{f}_t^{pp'}$, Eq. (11) is simply Eq. (8), i.e. the standard MPM; otherwise, the right hand side becomes $\sum_p w_{ip}\{\mathbf{G}_p(\mathbf{x}_i - \mathbf{x}_p) + \sum_{p' \in N_p} -k\log(d_{pp'})e_{p'p}\}$. Essentially, we model an exponential residual force which only impacts compression, not shearing or rotation, between particles.

Overall, our model learns parameters $k$ in Eq. (10) and Young's modulus $\epsilon$ in Eq. (8) through neural networks, equivalent to $\sigma^{cm}$ in Eq. (1) being learnable. Considering that these parameters vary for each particle at every time step, we estimate them via:

$$ k_p = NN_k(\mathbf{x}_p, \mathbf{v}_p, N_p), \quad \epsilon_p = NN_\epsilon(\mathbf{x}_p, \mathbf{v}_p, N_p). \quad (12) $$

where $N_p$ represents the neighborhood of particle $p$. $NN_k$ and $NN_\epsilon$ include multiple Continuous Convolution (CConv) and fully-connected (FC) layers. The detailed architectures are in the SM.

## 3.4. New Learnable Active Forces $\mathbf{f}^{act}$

Besides the new stress, we know that crowds are also influenced by active forces from individuals [4], e.g. when one person stiffens the body to restore balance, or move with the neighbors. Overall, these local actions manifest as random active forces from individuals to the whole system, which can be described by the Toner-Tu (TT) equation [51]:

$$ \frac{\partial \mathbf{v}}{\partial t} + \lambda(\mathbf{v} \cdot \nabla)\mathbf{v} = \alpha \mathbf{v} - \beta |\mathbf{v}|^2 \mathbf{v} - \nabla P + D_L \nabla (\nabla \cdot \mathbf{v}) + D_1 \nabla^2 \mathbf{v} + D_2 (\mathbf{v} \cdot \nabla)^2 \mathbf{v} + \tilde{\mathbf{f}}, \quad (13) $$

where $\beta$, $D_1$, $D_2$, $D_L$ are all positive, and $\alpha < 0$ in disordered phases and $\alpha > 0$ in ordered state. Specifically, $\alpha \mathbf{v}$ and $-\beta |\mathbf{v}|^2 \mathbf{v}$ represent self-driven motions and velocity saturation, essential for modeling the collective behavior in high-density active systems. The diffusion terms $D_L \nabla(\nabla \cdot \mathbf{v})$, $D_1 \nabla^2 \mathbf{v}$ and $D_2 (\mathbf{v} \cdot \nabla)^2 \mathbf{v}$ handle the spread of velocity perturbations, while $\tilde{\mathbf{f}}$ is a Gaussian random force.

However, we do not directly use the TT equation for the active force $\mathbf{f}^{act}$ in Eq. (1), because we need to learn the free parameters ($\alpha$, $\beta$, $D_1$, $D_2$, $D_L$). Therefore, we first remove $\nabla P$ as it is caused by the Cauchy stress and is learned in Sec. 3.3. Then we separate the right hand side of TT into two parts, $\alpha \mathbf{v}$ and the rest. This is because $\alpha \mathbf{v}$ dominates motion alignment caused by active forces. We learn $\alpha$ via $\alpha = NN_\alpha(\mathbf{v})$, where the neural network $NN_\alpha$ is detailed in the SM.

To learn the rest, we observe that although $\tilde{\mathbf{f}}$ is a Gaussian random force, the overall distribution of the rest, i.e.

<figure><img src="image_6.png" /><figcaption>Figure 3. Four crowd scenes with diverse dynamics. The yellow arrows indicate the main movement.</figcaption></figure>

$-\beta |\mathbf{v}|^2 \mathbf{v} + D_L \nabla(\nabla \cdot \mathbf{v}) + D_1 \nabla^2 \mathbf{v} + D_2 (\mathbf{v} \cdot \nabla)^2 \mathbf{v} + \tilde{\mathbf{f}}$ is not Gaussian. Therefore, we assume its distribution is a Gaussian in a latent space, captured by a latent variable $z \sim Normal$, so that:

$$ \mathbf{f}^{act} = NN_{act}(|\mathbf{v}|^2 \mathbf{v}, \nabla(\nabla \cdot \mathbf{v}), \nabla^2 \mathbf{v}, (\mathbf{v} \cdot \nabla)^2 \mathbf{v}), z) \quad (14) $$

We formulate $NN_{act}$ as the decoder of a Conditional Variational Autoencoder which is detailed in the SM.

## 3.5. Body Forces

Other than the previous forces, there can be global forces that resembles body forces. For instance, the crowd can be attracted to a goal, or circulate around an area, which can be formulated as a global force that is evenly applied to every unit in the continuum. This is modeled by the $f_p = \rho b$ in Eq. (1). For example, for a goal attraction force [21], we use $\mathbf{f}_p^{gl} = \frac{1}{\tau}(\mathbf{v}_d \tilde{e}_{pg} - \mathbf{v}_p)$ where $\frac{1}{\tau} = \frac{m_p}{\Delta t}$ and $e_{pg} = \frac{\mathbf{x}_g - \mathbf{x}_p}{||\mathbf{x}_g - \mathbf{x}_p||}$ where $\mathbf{x}_g$ is the goal and $\mathbf{v}_d$ is the preferred speed. For circular motions, we can define a centripetal force: $\mathbf{f}_p^{cen} = m_p \frac{\mathbf{v}_p^2}{r}$ where $r$ is the circle radius.

## 3.6. Training

Overall, the learnable parameters are the Young's modulus $\epsilon$ and $k$ in Eq. (12), the $NN_\alpha$, and the $NN_{act}$ in Eq. (14). Owing to the full differentiability of our model, we learn the model via Adam. More details are in the SM.

## 4. Experiments

**Datasets** We employ Evacuation drill [13], Hajj [11], Hellfest [23] and Marathon [5] (Fig. 3). Evacuation drill include multiple videos, each containing ~100 participants who were instructed to leave through an exit as quickly as they could. This dataset are divided into three subsets based on competitiveness: Low Competitiveness (Drill1), Moderate Competitiveness (Drill2) and High Competitiveness (Drill3). Hajj and Hellfest are two sets of videos from YouTube. Pilgrims in Hajj walk around the Kaaba in a counter-clockwise direction. In Hellfest, two groups of people intentionally run into one another at a high speed. The density of Evacuation, Hajj and Hellfest easily reaches $6-8p/m^2$. Marathon is selected from UCF Dataset[5], where people run in a U-shaped road. We exclude synthetic datasets as existing simulators fail to generate realistic crowd dynamics in extremely high-density scenarios.
ios [40]. Consequently, models that can learn well on synthetic data cannot guarantee to work well on real crowds.

Similar to other differentiable physics models [15], our model does not need large amounts of training data, because it essentially learns a parameterized PDE. In Drill$_1$, Drill$_2$ and Drill$_3$, we randomly select 4 videos for training, 2 videos for validation and 2 videos for testing. For Hajj and Marathon, we divide the video frames into training, validation and testing splits using a ratio of 60%, 20% and 20% respectively. For the Hellfest dataset, we use a split ratio of 40% for training, 30% for validation, and 30% for testing. Initially, we compute pixel-wise optical flow between consecutive frames using denseflow [61]. Then we discretize the image space into grids and treat every optical flow vector as a velocity vector of a particle and use the P2G transfer (introduced in 3.2) to estimate the velocity field.

**Metrics** Since forecasting high-density crowds from in-the-wild videos is a new topic, there is no established evaluation protocol and metrics to our best knowledge. Therefore, we employ both quantitative and qualitative comparison. On the raw data, we employ the Mean Square Error (MSE) on the optical flow (*Err*<sub>flow</sub>). The predicted velocity field is converted into optical flow via G2P transfer (Sec. 3.2). Furthermore, to evaluate the field locally, we use MSE on the velocity field (*Err*<sub>vel</sub>).

## 4.1. Prediction

**Comparison baselines** Since there is no similar existing method, we adapt several closest methods as baselines. They include physics models (Baseline<sub>1</sub> [48]), physics-informed neural network (HINN [76]), deep learning methods for scene parsing (PredFlow [29]), sequence prediction (PredRNNv2 [63]), and video prediction (SimVP [12], TAU [49]). These baselines span a diverse range of architectures and tasks. We give details of how we adapt them for comparison in the SM.

Given an input up to frame $t$, we test the prediction on the frame $t+n$. If $n$ is too large, the prediction task is ill-defined as there is unobserved external impact on the crowds; if it is too small, there is too little meaningful dynamics. So after considering the total length of the data and the crowd dynamics, we choose a reasonably large $n$, corresponding to 8.2s in Drill$_{1-3}$, 2s, 1.6s, and 2.4s in Hajj, Hellfest and Marathon respectively.

We show the prediction results in Tab. 1. Our results are the average and the best result based on 10 trials. Overall, our method achieves the highest performance across both metrics on all datasets except Hajj where our method is slightly worse than SimVP. In Drill$_1$ and Marathon, we obtain the best Err$_{vel}$ and the second best Err$_{flow}$. For further investigation, we compare predictions across varying horizons in Fig. 4. Overall, our method achieves better predictions

especially in long-horizon prediction. There are two exceptions. One exception is the optical flow in Marathon. A further investigation suggests that the main reason is the main moving crowds in this scene only occupies a portion of the whole space. However, the optical flow can exist in other places due to noises. A pure data-driven approach tend to fit the optical flows in these areas better. Therefore, our method predicts the best on the velocity field but not in optical flows, as these noises are filtered in the P2G step. The other exception is Hajj. This is because Hajj contains a very peaceful crowd slowly moving around one point. Therefore, this is actually an easier prediction task, where nearly all methods perform reasonably well.

**Continuous vs discrete time modeling** One distinctive feature of our method is that it is a continuous-time model, which does not require data to be observed on evenly discretized timeline. This is different from the other deep learning methods which are discrete-time models requiring observations to be evenly spaced by a fixed $\Delta t$. We show more results in the SM.

**Input requirements and parameter efficiency** Our approach demonstrates a significant advantage in both input requirements and parameter efficiency. PredFlow requires four consecutive RGB frames and has 42.534M parameters, while PredRNNv2, SimVP, and TAU demand ten consecutive frames, with 24.216M, 18.604M, and 44.657M parameters, respectively. Although Baseline$_1$ operates on a single optical flow frame, it needs laborious hand-tuning. Our model achieves mostly the best long-horizon prediction with only a single optical flow frame as input, and merely 3.06M parameters. Vastly different from other models, our input contains little dynamics. The substantial reduction in parameter count, coupled with the minimal input requirements, highlights our model's efficiency, especially when data and computational resources are limited.

## 4.2. Prediction under Unseen Behavior

For further generalization tests, we conduct cross-behavior testing, i.e. trained on one behavior and tested on another. Specifically, we use the Drill as they are conducted in the same environment, with only the behavior changed. People were asked to be peaceful, neural and aggressive in *Drill*$_1$, *Drill*$_2$ and *Drill*$_3$ respectively. These scenes share some common behaviors and cross-behavior prediction can test whether a model can learn such behaviors.

Here, we conduct two groups of experiments by training the model on *Drill*$_1$ and test the model on *Drill*$_2$ and *Drill*$_3$. Quantitative results are reported in Fig. 5. Obviously, as behavior shifts from peaceful to neutral and then to aggressive, model performance declines consistently across both velocity field and optical flow metrics. This suggests
<table><caption>Table 1. Prediction comparison. The red, orange and yellow colors represent the top three results. Our method achieves up to 18.69% and 12.11% improvement in Err<sub>vel</sub> and Err<sub>flow</sub> respectively on Drill<sub>2</sub>.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Drill<sub>1</sub></th><th colspan="2">Drill<sub>2</sub></th><th colspan="2">Drill<sub>3</sub></th><th colspan="2">Hajj</th><th colspan="2">Hellfest</th><th colspan="2">Marathon</th></tr><tr><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th></tr></thead><tbody><tr><td>PredFlow</td><td>1.0669</td><td>1.6944</td><td>10.260</td><td>13.944</td><td>2.4555</td><td>3.7547</td><td>2246.6</td><td>3581.2</td><td>8.8040</td><td>14.040</td><td>122.86</td><td>621.70</td></tr><tr><td>Baseline<sub>1</sub></td><td>0.7555</td><td>1.1455</td><td>1.3319</td><td>1.8826</td><td>2.1150</td><td>3.0985</td><td>0.9354</td><td>1.2511</td><td>3.5151</td><td>5.2745</td><td>2.8778</td><td>11.820</td></tr><tr><td>PredRNNv2</td><td>1.3239</td><td>1.7329</td><td>1.3170</td><td>1.9586</td><td>1.6613</td><td>2.8734</td><td>0.6805</td><td>1.0022</td><td>3.3374</td><td>5.4267</td><td>2.0510</td><td>10.983</td></tr><tr><td>SimVP</td><td>2.2364</td><td>2.6801</td><td>5.6415</td><td>6.8968</td><td>2.9760</td><td>4.2312</td><td>0.6212</td><td>0.9160</td><td>4.9703</td><td>7.3458</td><td>1.6636</td><td>10.633</td></tr><tr><td>TAU</td><td>0.5330</td><td>0.9466</td><td>2.8687</td><td>3.6093</td><td>2.6546</td><td>3.9093</td><td>0.6791</td><td>0.9778</td><td>5.3888</td><td>7.8822</td><td>1.9015</td><td>10.966</td></tr><tr><td>HINN</td><td>0.5618</td><td>0.9224</td><td>1.1187</td><td>1.6603</td><td>2.6590</td><td>3.6835</td><td>1.1600</td><td>1.5402</td><td>7.1427</td><td>9.5330</td><td>4.3488</td><td>14.210</td></tr><tr><td>Ours(mean)</td><td>0.5284</td><td>0.9293</td><td>1.0721</td><td>1.6558</td><td>1.6461</td><td>2.7354</td><td>0.6591</td><td>0.9708</td><td>3.0457</td><td>4.8823</td><td>1.4927</td><td>10.794</td></tr><tr><td>Ours(best)</td><td>0.5259</td><td>0.9263</td><td>1.0708</td><td>1.6547</td><td>1.6435</td><td>2.7326</td><td>0.6502</td><td>0.9609</td><td>3.0422</td><td>4.8790</td><td>1.4581</td><td>10.759</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 4. Comparison on velocity field and optical flow. The right and left vertical axes are for PredFlow and the rest methods respectively, due to the scale of PredFlow results is much larger. Vertical axes are logscaled.</figcaption></figure>

that the model's capacity to generalize across behaviors decreases as the level of crowd competitiveness rises, indicating that aggressive crowd dynamics present more complex interactions that are harder for the model to predict.

## 4.3. Qualitative Evaluation

One distinctive feature of our method is, once learned, it is possible to use it as a simulator for analysis, which is not possible for all the pure data-driven baseline methods. This is especially important for high-density crowd management as it often requires to predict in 'what-if' situations which can be dangerous in real world and there is no data for reference. However, it is still an open question in

terms of how to quantitatively evaluate the simulation of high-density crowds especially if there is no data reference. Therefore, following the practice [53, 69], we conduct qualitative evaluation.

### 4.3.1. Crowd Material

We first design a simulation scenario where people are pushed into an enclosed area, to show how our method captures the empirical observations about 'crowd material' mentioned in Sec. 3.2, in comparison with a standard MPM model. It is a half-open area where people keep coming from the left and push the whole crowd to the right, shown in Fig. 6. Fig. 6 Left shows the stress map. Normally, to
<figure><img src="image_8.png" /><figcaption>Figure 5. Generalization. Trained on *Drill*₁, and tested on *Drill*₂ and *Drill*₃. Vertical axes are logscaled.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6. Crowd moving to the right with an initial velocity towards the right. Left: the stress colormap at t = 0.5s. Right: the crowd state at t = 1.8s. There are three experiments in the left and right group respectively. They are: (1) small stiffness in standard MPM $\epsilon$=1, (2) large stiffness in standard MPM $\epsilon$=100, (3) Ours.</figcaption></figure>

simulate such crowds with a physics model, one needs to hand-tune the key parameters, which in standard MPM is the Young's modulus $\epsilon$. However, no matter whether the stiffness is low or high in standard MPM (Fig. 6 Left 1 and 2) where the Young's modulus $\epsilon$ = 1 and $\epsilon$ = 100, it does not model the exponential resistance as our model (Fig. 6 Left 3). The overall stress increases steeply near the right end where the inter-person distance becomes smaller. This indicates even hand-tuning MPM will not mimic the crowd material as it does not exist in its parameter space.

Next, Fig. 6 Right shows the compression dominance. Fig. 6 Right 1 is too soft to resist compression, while Fig. 6 Right 2 resists compression, shear and rotation all together. Comparatively, our model (Fig. 6 Right 3) resists compression between particles but still allows shearing/rotation as highlighted in the red box. Again, this shows that the crowd material does not exist in the parameter space of standard MPM. Furthermore, the red box shows a ‘wave’ formation due to compression induced by the left side. Simulating wave formation is only recently addressed by [57], but [57] is not learnable, unlike our method.

### 4.3.2. Simulation in a New Environment

We also show simulation where the behaviors are first learned in one environment but simulated with changes to the environment. This is a common scenario in designing a public space for public safety.

We show several examples in Fig. 7, where the first row

<figure><img src="image_10.png" /><figcaption>Figure 7. Evacuation simulations in the original scene (top), with a wider exit (middle) and a pillar (bottom). The color indicates the magnitude of stress.</figcaption></figure>

shows the simulation on the original *Drill*₂ scene, the second and third rows are two different scenarios, one with a wider exit (second row) and one with a pillar in front of the exit (third row). Intuitively, a wider exit reduces the stress at the exit. This is shown in the second row, showing our method simulates the correct behaviors. More interestingly, an obstacle in front of the exit, if placed well, can also reduce the stress (third row). This is observed in some circumstance e.g. a certain exit width and crowd density [7]. We are able to replicate similar effects. These simulations demonstrate that the behavior learned by our method can be effectively used for simulating alternative scenarios.
## 5. Discussion, Conclusion and Future Work

In this paper, we introduced a neural stochastic differential equation system to learn extreme high-density crowds as active matter. Our model is specifically designed for high-density crowds and is not suitable for low-density crowds. In crowd modeling, it is common practice that different densities require us to choose vastly different physical models. To our knowledge, no unified model exists for all crowd densities. Our method fills the gap of learning physical models for extremely high-density crowds, which is first of its kind, to the best of our knowledge. Additionally, we focus on scenarios involving large crowd movement, as slow or static crowds are typically less critical from a safety perspective.

While our model shows strong capabilities in capturing high-density crowd dynamics, it has certain limitations. A primary limitation is that its current single-view video input design restricts cross-view generalization, which can be improved by incorporating multiview perspectives when estimating the velocity field. Another challenge is that it does not explicitly model complex behaviors like queuing, grouping, lane formation, nor account for other space users like vehicles. These limitations can be mitigated by introducing additional inter-particle forces (e.g., grouping forces) or object detection modules. Finally, predicting entirely unseen scenes and crowds remains an open challenge. In the future, we will further generalize the model to incorporate the aforementioned scenarios, aiming for a wider range of real-world applications.

## Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No 899739 CrowdDNA.

## References

[1] Shahira Abousamra, Minh Hoai, Dimitris Samaras, and Chao Chen. Localization in the crowd with topological constraints. In *Proceedings of the AAAI Conference on Artificial Intelligence*, pages 872–881, 2021. 1

[2] Juliane Adrian, Armin Seyfried, and Anna Sieben. Crowds in front of bottlenecks at entrances from the perspective of physics and social psychology. *Interface*, 17(165):20190871 –, 2020. 2

[3] Alexandre Alahi, Kratarth Goel, Vignesh Ramanathan, Alexandre Robicquet, Li Fei-Fei, and Silvio Savarese. Social lstm: Human trajectory prediction in crowded spaces. In *Proceedings of the IEEE conference on computer vision and pattern recognition*, pages 961–971, 2016. 1, 2

[4] Saad Ali. *Taming crowded visual scenes*. University of Central Florida, 2008. 5

[5] Saad Ali and Mubarak Shah. A lagrangian particle dynamics approach for crowd flow segmentation and stability analysis.

In 2007 IEEE Conference on Computer Vision and Pattern Recognition, pages 1–6. IEEE, 2007. 1, 5

[6] Saad Ali and Mubarak Shah. Floor fields for tracking in high density crowd scenes. In *European conference on computer vision*, pages 1–14. Springer, 2008. 2

[7] Tom Axblad and Alvaro Ortega Gonzalez. Evacuation with obstacles in real-time using crowd simulation, 2018. 8

[8] Yizheng Cai, Nando de Freitas, and James J Little. Robust visual tracking for multiple targets. In *European conference on computer vision*, pages 107–118. Springer, 2006. 2

[9] Hongyi Chen, Jingtao Ding, Yong Li, Yue Wang, and Xiao-Ping Zhang. Social physics informed diffusion model for crowd simulation. In *Proceedings of the AAAI Conference on Artificial Intelligence*, pages 474–482, 2024. 1

[10] Michael A Day. The no-slip condition of fluid dynamics. *Erkenntnis*, 33(3):285–296, 1990. 4

[11] Dreamerr86. Azan & muslim pilgrims "hajj" prayer in mecca kaaba. https://www.youtube.com/watch?v=KGukAoiGhZU&t=5s, 2013. 5

[12] Zhangyang Gao, Cheng Tan, Lirong Wu, and Stan Z Li. Simvp: Simpler yet better video prediction. In *Proceedings of the IEEE/CVF conference on computer vision and pattern recognition*, pages 3170–3180, 2022. 2, 6, 13

[13] Angel Garcimartín, Daniel R Parisi, Jose M Pastor, César Martín-Gómez, and Iker Zuriguel. Flow of pedestrians through narrow doors with different competitiveness. *Journal of Statistical Mechanics: Theory and Experiment*, 2016 (4):043402, 2016. 5

[14] Abhinav Golas, Rahul Narain, and Ming C Lin. Continuum modeling of crowd turbulence. *Physical review E*, 90(4):042816, 2014. 1

[15] Deshan Gong, Zhanxing Zhu, Andrew J Bulpitt, and He Wang. Fine-grained differentiable physics: a yarn-level model for fabrics. *arXiv preprint arXiv:2202.00504*, 2022. 2, 6

[16] Agrim Gupta, Justin Johnson, Li Fei-Fei, Silvio Savarese, and Alexandre Alahi. Social gan: Socially acceptable trajectories with generative adversarial networks. In *Proceedings of the IEEE conference on computer vision and pattern recognition*, pages 2255–2264, 2018. 1, 2

[17] Feixiang He, Yuanhang Xia, Xio Zhao, and He Wang. Informative scene decomposition for crowd analysis, comparison and simulation guidance. *ACM Transaction on Graphics (TOG)*, 4(39), 2020. 2

[18] Feixiang He, Yuanhang Xiang, Xi Zhao, and He Wang. Informative scene decomposition for crowd analysis, comparison and simulation guidance. *ACM Transactions on Graphics (TOG)*, 39(4):50–1, 2020. 1

[19] Gao-qi He, Yu Yang, Zhi-hua Chen, Chun-hua Gu, and Zhi-geng Pan. A review of behavior mechanisms and crowd evacuation animation in emergency exercises. *Journal of Zhejiang University SCIENCE C*, 14(7):477–485, 2013. 1

[20] Laure Heigeas, Annie Luciani, Joelle Thollot, and Nicolas Castagné. A physically-based particle model of emergent crowd behaviors. *arXiv preprint arXiv:1005.4405*, 2010. 4

[21] Dirk Helbing and Peter Molnar. Social force model for pedestrian dynamics. *Physical review E*, 51(5):4282, 1995. 1, 5
[22] Dirk Helbing, Anders Johansson, and Habib Zein Al-Abideen. Dynamics of crowd disasters: An empirical study. *Phys. Rev. E*, 75:046109, 2007. 2

[23] HoToox. Wall of death (extreme) - with full force 2014. https://www.youtube.com/watch?v=ySPlanMCmM4, 2014. 5

[24] Min Hu, Saad Ali, and Mubarak Shah. Learning motion patterns in crowded scenes using motion flow field. In *2008 19th International Conference on Pattern Recognition*, pages 1–5. IEEE, 2008. 2

[25] Yuanming Hu, Yu Fang, Ziheng Ge, Ziyin Qu, Yixin Zhu, Andre Pradhana, and Chenfanfu Jiang. A moving least squares material point method with displacement discontinuity and two-way rigid body coupling. *ACM Transactions on Graphics* (TOG), 37(4):1–14, 2018. 3, 18

[26] Chenfanfu Jiang, Craig Schroeder, Andrew Selle, Joseph Teran, and Alexey Stomakhin. The affine particle-in-cell method. *ACM Transactions on Graphics* (TOG), 34(4):1–10, 2015. 3, 18

[27] Chenfanfu Jiang, Craig Schroeder, Joseph Teran, Alexey Stomakhin, and Andrew Selle. The material point method for simulating continuum materials. In *Acm siggraph 2016 courses*, pages 1–52. 2016. 2, 3, 17, 18, 19, 20

[28] Hao Jiang, Wenbin Xu, Tianlu Mao, Chunpeng Li, Shihong Xia, and Zhaoqi Wang. Continuum crowd simulation in complex environments. *Computers & Graphics*, 34(5):537–544, 2010. 2

[29] Xiaojie Jin, Huaxin Xiao, Xiaohui Shen, Jimei Yang, Zhe Lin, Yunpeng Chen, Zequn Jie, Jiashi Feng, and Shuicheng Yan. Predicting scene parsing and motion dynamics in the future. *Advances in neural information processing systems*, 30, 2017. 6, 13

[30] Julio Cezar Silveira Jacques Junior, Soraia Raupp Musse, and Claudio Rosito Jung. Crowd analysis using computer vision techniques. *IEEE Signal Processing Magazine*, 27(5): 66–77, 2010. 1

[31] Ioannis Karamouzas, Nick Sohre, Ran Hu, and Stephen J Guy. Crowd space: a predictive crowd analysis technique. *ACM Transactions on Graphics (TOG)*, 37(6):1–14, 2018. 1

[32] Kang Hoon Lee, Myung Geol Choi, Qyoun Hong, and Jehee Lee. Group behavior from video: a data-driven approach to crowd simulation. In *Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation*, pages 109–118, 2007. 1

[33] Dingkang Liang, Wei Xu, and Xiang Bai. An end-to-end transformer model for crowd localization. In *Computer Vision-ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part 1*, pages 38–54. Springer, 2022. 1

[34] Junbang Liang, Ming Lin, and Vladlen Koltun. Differentiable cloth simulation for inverse problems. *Advances in Neural Information Processing Systems*, 32, 2019. 2

[35] Chengxin Liu, Hao Lu, Zhiguo Cao, and Tongliang Liu. Point-query quadtree for crowd counting, localization, and more. In *Proceedings of the IEEE/CVF International Conference on Computer Vision*, pages 1676–1685, 2023. 1

[36] Xinyan Liu, Guorong Li, Yuankai Qi, Zhenjun Han, Anton van den Hengel, Nicu Sebe, Ming-Hsuan Yang, and Qingming Huang. Consistency-aware anchor pyramid network for crowd localization. *IEEE transactions on pattern analysis and machine intelligence*, 2024. 1

[37] Axel López, François Chaumette, Eric Marchand, and Julien Pettré. Character navigation in dynamic environments based on optical flow. In *Computer Graphics Forum*, pages 181–192. Wiley Online Library, 2019. 1

[38] A NSAV Marana, Sergio A Velastin, L da F Costa, and RA Lotufo. Automatic estimation of crowd density using texture. *Safety Science*, 28(3):165–175, 1998. 13

[39] Zhaobin Mo, Yongjie Fu, and Xuan Di. Pi-neugode: Physics-informed graph neural ordinary differential equations for spatiotemporal trajectory prediction. In *Proceedings of the 23rd International Conference on Autonomous Agents and Multiagent Systems*, pages 1418–1426, 2024. 1

[40] Rahul Narain, Abhinav Golas, Sean Curtis, and Ming C Lin. Aggregate dynamics for dense crowd simulation. In *ACM SIGGRAPH Asia 2009 papers*, pages 1–8. 2009. 1, 2, 4, 6

[41] Erçağ Pinçe, Sabareesh KP Velu, Agnese Callegari, Parviz Elahi, Sylvain Gigan, Giovanni Volpe, and Giorgio Volpe. Disorder-mediated crowd control in an active matter system. *Nature communications*, 7(1):10907, 2016. 1

[42] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686–707, 2019. 2

[43] Siyuan Shen, Yang Yin, Tianjia Shao, He Wang, Chenfanfu Jiang, Lei Lan, and Kun Zhou. High-order differentiable autoencoder for nonlinear model reduction. *arXiv preprint arXiv:2102.11026*, 2021. 2

[44] Anna Sieben and Armin Seyfried. Inside a life-threatening crowd: Analysis of the love parade disaster from the perspective of eyewitnesses. *Safety Science*, 166:106229, 2023. 2

[45] Kuldeep Singh, Shantanu Rajora, Dinesh Kumar Vishwakarma, Gaurav Tripathi, Sandeep Kumar, and Gurjit Singh Walia. Crowd anomaly detection using aggregation of ensembles of fine-tuned convnets. *Neurocomputing*, 371: 188–198, 2020. 2

[46] Yanjie Song, He Wang, He Yang, Maria Luisa Taccari, and Xiaohui Chen. Loss-attentional physics-informed neural networks. *Journal of Computational Physics*, page 112781, 2024. 2

[47] Prof. Dr. G. Keith Still. Crowd safety and crowd risk analysis. 1

[48] Andre Pradhana Tampubolon, Theodore Gast, Gergely Klár, Chuyuan Fu, Joseph Teran, Chenfanfu Jiang, and Ken Museth. Multi-species simulation of porous sand and water mixtures. *ACM Transactions on Graphics (TOG)*, 36(4): 1–11, 2017. 4, 6, 13, 16

[49] Cheng Tan, Zhangyang Gao, Lirong Wu, Yongjie Xu, Jun Xia, Siyuan Li, and Stan Z Li. Temporal attention unit: Towards efficient spatiotemporal predictive learning. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, pages 18770–18782, 2023. 2, 6, 13
[50] Roger Temam. Navier-Stokes equations: theory and numerical analysis. American Mathematical Soc., 2001. 3, 17

[51] John Toner and Yuhai Tu. Long-range order in a two-dimensional dynamical xy model: how birds fly together. Physical review letters, 75(23):4326, 1995. 2, 5

[52] Adrien Treuille, Seth Cooper, and Zoran Popović. Continuum crowds. ACM Transactions on Graphics (TOG), 25(3): 1160–1168, 2006. 1, 3, 16

[53] Branislav Ulicny and Daniel Thalmann. Towards interactive real-time crowd behavior simulation. In Computer Graphics Forum, pages 767–775. Wiley Online Library, 2002. 7

[54] Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian fluid simulation with continuous convolutions. In International Conference on Learning Representations, 2019. 12

[55] Jur Van den Berg, Ming Lin, and Dinesh Manocha. Reciprocal velocity obstacles for real-time multi-agent navigation. In 2008 IEEE international conference on robotics and automation, pages 1928–1935. IEEE, 2008. 1

[56] Wouter Van Toll and Julien Pettré. Algorithms for microscopic crowd simulation: Advancements in the 2010s. In Computer Graphics Forum, pages 731–754. Wiley Online Library, 2021. 1

[57] Wouter van Toll, Thomas Chatagnon, Cédric Braga, Barbara Solenthaler, and Julien Pettré. Sph crowds: Agent-based crowd simulation up to extreme densities using fluid dynamics. Computers & Graphics, 98:306–321, 2021. 1, 8

[58] He Wang and Carol O'Sullivan. Globally continuous and non-markovian crowd activity analysis from videos. In European conference on computer vision, pages 527–544. Springer, 2016. 1, 2

[59] He Wang, Jan Ondřej, and Carol O'Sullivan. Path patterns: Analyzing and comparing real and simulated crowds. In ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games 2016, pages 49–57, 2016.

[60] He Wang, Jan Ondřej, and Carol O'Sullivan. Trending paths: A new semantic-level metric for comparing simulated and real crowd data. IEEE Transactions on Visualization and Computer Graphics, PP(99):1–1, 2016. 1, 2

[61] Shiguang* Wang, Zhizhong* Li, Yue Zhao, Yuanjun Xiong, Limin Wang, and Dahua Lin. denseflow. https://github.com/open-mmlab/denseflow, 2020. 6

[62] Yunbo Wang, Zhifeng Gao, Mingsheng Long, Jianmin Wang, and S Yu Philip. Predrnn++: Towards a resolution of the deep-in-time dilemma in spatiotemporal predictive learning. In International conference on machine learning, pages 5123–5132. PMLR, 2018. 2

[63] Yunbo Wang, Haixu Wu, Jianjin Zhang, Zhifeng Gao, Jianmin Wang, S Yu Philip, and Mingsheng Long. Predrnn: A recurrent neural network for spatiotemporal predictive learning. IEEE Transactions on Pattern Analysis and Machine Intelligence, 45(2):2208–2225, 2022. 2, 6, 13

[64] Keenon Werling, Dalton Omens, Jeongseok Lee, Ioannis Exarchos, and C Karen Liu. Fast and feature-complete differentiable physics for articulated rigid bodies with contact. arXiv preprint arXiv:2103.16021, 2021. 2

[65] David Wolinski, S J. Guy, A-H Olivier, Ming Lin, Dinesh Manocha, and Julien Pettré. Parameter estimation and comparative evaluation of crowd simulations. In Computer Graphics Forum, pages 303–312. Wiley Online Library, 2014. 1

[66] Shuang Wu, Hua Yang, Shibao Zheng, Hang Su, Yawen Fan, and Ming-Hsuan Yang. Crowd behavior analysis via curl and divergence of motion trajectories. International Journal of Computer Vision, 123(3):499–519, 2017. 1, 2, 14, 15, 16

[67] Beihao Xia, Conghao Wong, Qinmu Peng, Wei Yuan, and Xinge You. Cscnet: Contextual semantic consistency network for trajectory prediction in crowded spaces. Pattern Recognition, 126:108552, 2022. 1

[68] Wei Xiang, YIN Haoteng, He Wang, and Xiaogang Jin. Socialvae: Predicting pedestrian trajectory via interaction conditioned latents. In Proceedings of the AAAI Conference on Artificial Intelligence, pages 6216–6224, 2024. 1

[69] Shanwen Yang, Tianrui Li, Xun Gong, Bo Peng, and Jie Hu. A review on crowd simulation and modeling. Graphical Models, 111:101081, 2020. 7

[70] Yang Yang, Jingen Liu, and Mubarak Shah. Video scene understanding using multi-scale analysis. In 2009 IEEE 12th International Conference on Computer Vision, pages 1669–1676. IEEE, 2009. 2

[71] J. Yue, D. Manocha, and H. Wang. Human trajectory prediction via neural social physics. In Proceedings of the European Conference on Computer Vision (ECCV), 2022. 1, 2

[72] Jiangbei Yue, Dinesh Manocha, and He Wang. Human trajectory forecasting with explainable behavioral uncertainty, 2023. 1, 2

[73] Zheyan Zhang, Yongxing Wang, Peter K. Jimack, and He Wang. Meshingnet: A new mesh generation method based on deep learning. In Computational Science – ICCS 2020, pages 186–198, Cham, 2020. Springer International Publishing. 2

[74] Zheyan Zhang, Peter K. Jimack, and He Wang. MeshingNet3D: Efficient generation of adapted tetrahedral meshes for computational mechanics. Advances in Engineering Software, 157-158, 2021. 2

[75] Yiqi Zhong, Luming Liang, Ilya Zharkov, and Ulrich Neumann. Mmvp: Motion-matrix-based video prediction. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 4273–4283, 2023. 2

[76] Yanshan Zhou, Pingrui Lai, Jiaqi Yu, Yingjie Xiong, and Hua Yang. Hydrodynamics-informed neural network for simulating dense crowd motion patterns. In Proceedings of the 32nd ACM International Conference on Multimedia, pages 4553–4561, 2024. 1, 2, 6
# Learning Extremely High Density Crowds as Active Matters
## Supplemental Material

Feixiang He², Jiangbei Yue³, Jialin Zhu², Armin Seyfried⁴, Dan Casas⁵, Julien Pettré⁶, He Wang ¹,²

¹AI Centre, University College London, UK ²University College London, UK

³University of Leeds, UK ⁴Forschungszentrum Jülich, Germany

⁵Universidad Rey Juan Carlos, Spain ⁶INRIA Rennes, France

# 1. Additional Technical Details

## 1.1. Loss function

During training, We employ the Mean Square Error loss between the predicted velocity field {$\hat{v}^t$}$_{t=1}^T$ and the ground truth {$v^t$}$_{t=0}^T$ to optimize the model:

$$ \min \frac{1}{T} \int_{t=1}^{T} \|v^t - \hat{v}^t\|_2^2 \approx \frac{1}{T} \sum_{t=1}^{T} \|v^t - \hat{v}^t\|_2^2 \quad (15) $$

where we can only approximate the first integral as the summation over observed frames. But the nature of the approximation does not dictate how the frames are distributed on the timeline. They can be evenly or unevenly distributed.

## 1.2. Error Metrics

To evaluate our proposed model, we employ two metrics: Err$_{vel}$ and Err$_{flow}$. As presented in Sec. 4 of the main paper, Err$_{vel}$ is used to assess the velocity field locally and is calculated as:

$$ Err_{vel} = \frac{1}{T} \sum_{t=1}^{T} \|v^t - \hat{v}^t\|_2^2 \quad (16) $$

Differently, Err$_{flow}$ operates directly on the raw optical flow data and is defined as:

$$ Err_{flow} = \frac{1}{T} \sum_{t=1}^{T} \|o^t - \hat{o}^t\|_2^2 \quad (17) $$

Here, {$\hat{o}^t$}$_{t=1}^T$ and {$o^t$}$_{t=0}^T$ represent the predicted optical flow and the ground truth respectively. Note that our model predicts the velocity field instead of directly predicting optical flow. Therefore, the predicted velocity field is converted into optical flow via G2P transfer, as detailed in sec.3.2 of the main paper.

## 1.3. Neural Networks for Crowd Material

As illustrated in Sec 3.6 in the main paper, there are two types of learnable parameters: Young's modulus $\epsilon$ and $k$. These parameters are estimated using neural networks, specifically $NN_\epsilon$ for $\epsilon$ and $NN_k$ for $k$. Both networks share the same architecture, as shown in Fig. 8. The design is mainly to consider the motions of particles and its neighboring particles. Since we model a continuum where particles

<figure><img src="image_11.png" /><figcaption>Figure 8. Parameter estimation network. Operations annotated by * are followed by a ReLU activation.</figcaption></figure>

can be anywhere and not necessarily evenly spaced, one of the component in the network is the Continuous Convolution or CCov, inspired by the convolution layer designed originally for a Lagrangian method [54]. The continuous convolution captures the relative motions of the particles in the neighborhood, which dictates the local deformation of the field. In addition, we assume that the material parameter is different for every particle (pedestrian) in space and time. Therefore, we also use a fully-connected component (FC) to allow the individual particle motion to play a bigger role compared with its neighbors whose motions are already considered in the continuous convolution of the neighborhood. After several repeated layers of CCov and FC components, the network finally outputs the target parameter.

## 1.4. Conditional Variational Autoencoder for Random Active Forces

As outlined in the main paper, active force learning is divided into two parts. The first part, motion alignment (represented by $\alpha v$), is learned through a straightforward 6-layer convolutional neural network $NN_\alpha$, where the kernel size, stride, and padding are set to 3, 1, and 1, respectively. Each layer, except the final one, is followed by a Tanh activation function. The channel configuration for the layers is 32, 64, 128, 64, 32, and 1, respectively. The input of $NN_\alpha$ is the velocity field at the current step, and its output is the corresponding $\alpha$.

Secondly, to model the rest (expressed as $-\beta|v|^2v + D_L\nabla(\nabla \cdot v) + D_1\nabla^2v + D_2(v \cdot \nabla)^2v + \hat{f}$), a conditional variational autoencoder is utilized, as illustrated in Fig. 9. The CVAE takes two inputs: $x$, representing the ground truth of the remaining active forces, and $y$, serving as the condition. The condition is composed of four terms $|v|^2v, \nabla(\nabla \cdot v), \nabla^2v$ and $(v \cdot \nabla)^2v$, and processed through an embedding layer. Note that the decoder consists of multiple components ($D_0, D_1, D_2, \dots$) which generate intermediate
<figure><img src="image_12.png" /><figcaption>Figure 9. Architecture of CVAE.</figcaption></figure>

outputs $(\hat{x}_0, \hat{x}_1, \hat{x}_2, \dots)$. These intermediate results are com-
bined using a weighted sum with weights $(\omega_0, \omega_1, \omega_2, \dots)$
to produce the final output $\hat{x}$. This architecture allows the
CVAE to effectively model the complex relationships be-
tween the input and the conditional data, making it suitable
for learning the intricate dynamics.

1.5. Optical flows to/from grid

Optical flows can be transferred to and from a grid repre-
sentation to for training, testing, analysis and comparison.
When transferring optical flow data to a grid, we make use
of the particle-to-grid operation in the MPM (P2G). Each
pixel in the optical flow field is treated as an individual par-
ticle, with its velocity mapped to the corresponding grid cell
using a B-spline function. Conversely, when transforming
the velocity field from a grid back to its corresponding opti-
cal flows, we use the grid-to-particle operation in the MPM
(G2P). G2P computes the velocity at grid nodes back to
each pixel, which is treated as the optical flow at that pixel.

This P2G then G2P mapping between optical flows
and the velocity field, aside from facilitate learning, has a
smoothening effect due to the B-spline based interpolation,
effectively removing random noises from the noisy optical
flows. In addition, this mapping also preserve the spatial
and temporal dynamics. This dual representation enables
the comparisons of the particle-level details and the grid-
level velocity.

1.6. Hyperparameters

In all datasets, the density is fixed at 1. This is because it is
not possible estimate each individual’s mass, and if it were,
there is no commonly agreed way of defining the density in
high-density crowds [38]. By defining the density to be 1,
the particle mass depends on the particle radius. The parti-
cle radii are set to 20 pixels on *Drill*₁₋₃, 5 pixels on *Hajj*
and *Hellfest*, 10 pixels on *Marathon*. Note this setting
will not affect the learning of the dynamics as the dynamics
are learned mainly through the material parameters and the
active force.

For training, we employ Adam optimizer with an initial
learning rate of 1×10⁻⁴. To dynamically adjust the learning
rate during training, we apply a LambdaLR scheduler with

a decay function defined as $lr = 1 \times 10^{-4} \times 0.9^{epoch/50}$,
epoch is the current iteration. The batch size is set to 4.

1.7. Adaptation of baseline methods

To demonstrate the superiority of out method, we adapt sev-
eral closest methods as baselines. **Baseline**<sub>*I*</sub> is a physics-
based method that models crowds using a weakly incom-
pressible fluid framework [48]. PredFlow [29] is a deep
learning method designed for predicting future scene pars-
ing and motion dynamics. Additionally, we include several
sequence prediction models, PredRNNv2 [63], SimVP [12],
and TAU [49], which are also pure deep learning ap-
proaches and intended for video prediction tasks. PredFlow
takes four consecutive RGB frames as input and outputs the
future optical flows in an autoregressive fashion, **Baseline**<sub>*I*</sub>,
however, simulates velocity fields but we need to exhaus-
tively tune the parameters. For the video prediction meth-
ods, we adapt them to take optical flow as input and fore-
cast future optical flow sequences. Specifically, during the
training, we take 1.2s optical flow as input and forecast the
next 1.2s optical flow. These baselines span a diverse range
of architectures, from physics-based models to cutting-edge
deep learning frameworks, ensuring comprehensive cover-
age across methodologies.

1.8. Simulation

Our simulation runs ~ 40 FPS for 75 individuals. Like
any simulator, our method experiences a decrease in per-
formance when scaling to a large number of agents. Also,
albeit now shown in the paper, our method allows dynamic
adding or removing objects (dynamic or static), environ-
mental changes or events (e.g. explosion, evacuation). The
first two can be easily handled by boundary conditions,
while the latter can be realized by adding particle forces to
repel people from certain places, e.g. explosion, or attract
them, e.g. evacuation.

2. Additional Experimental Results

## 2.1. Continuous-time Prediction

As aforementioned in the main paper, one distinctive feature
of our method is that it is a continuous-time model, which
does not require data to be observed on evenly discretized
timeline. This is because the core of our method is a PDE
which can be solved for any horizon in future. This feature
provides good flexibility in terms of data collection. For
CCTV camera videos, sometimes there can be objects e.g.
tree leaves flying across the view or certain sudden change
of the lighting condition. Under these conditions, the opti-
cal flows will be extremely unreliable. This requires certain
data pre-processing step, e.g. abandoning frames that are of
too low quality. However, in this situation, the other base-
line methods will not be able to be directly trained on the
<table><caption>Table 2. Evaluation on temporal discontinuities on *Drill2*.</caption><thead><tr><th rowspan="2">Percent</th><th colspan="7">Err<sub>vel</sub></th><th colspan="7">Err<sub>flow</sub></th></tr><tr><th>1.2s</th><th>2.4s</th><th>3.6s</th><th>4.8s</th><th>6.0s</th><th>7.2s</th><th>8.4s</th><th>1.2s</th><th>2.4s</th><th>3.6s</th><th>4.8s</th><th>6.0s</th><th>7.2s</th><th>8.4s</th></tr></thead><tbody><tr><td>0%</td><td>1.1099</td><td>1.1748</td><td>1.1600</td><td>1.1332</td><td>1.1128</td><td>1.0954</td><td>1.0721</td><td>1.8016</td><td>1.8463</td><td>1.8072</td><td>1.7574</td><td>1.7185</td><td>1.6883</td><td>1.6558</td></tr><tr><td>50%</td><td>1.1078</td><td>1.1632</td><td>1.1419</td><td>1.1155</td><td>1.0966</td><td>1.0850</td><td>1.0640</td><td>1.7958</td><td>1.8303</td><td>1.7843</td><td>1.7349</td><td>1.6974</td><td>1.6729</td><td>1.6435</td></tr><tr><td>70%</td><td>1.1121</td><td>1.1637</td><td>1.1381</td><td>1.1069</td><td>1.0853</td><td>1.0705</td><td>1.0729</td><td>1.8086</td><td>1.8386</td><td>1.7891</td><td>1.7339</td><td>1.6942</td><td>1.6653</td><td>1.6378</td></tr></tbody></table>

<table><caption>Table 3. Evaluation with noises (1-3) on *Drill2*. 1: a zero-mean Gaussian noise with a Std (standard deviation); 2: a 2D uniform noise with Prob being the probability of values within [-0.7, 0.7] × [-0.8, 0.8]; 3: a mixture noise as the weighted sum of the Gaussian (G) and the uniform (U) noise.</caption><thead><tr><th>Std</th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Prob</th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>G/U</th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th></tr></thead><tbody><tr><td>0.1</td><td>1.0437</td><td>1.6281</td><td>0.1</td><td>1.0544</td><td>1.6380</td><td>0.3/0.7</td><td>1.0709</td><td>1.6476</td></tr><tr><td>1</td><td>1.0615</td><td>1.6438</td><td>0.7</td><td>1.0491</td><td>1.6347</td><td>0.5/0.5</td><td>1.0808</td><td>1.6523</td></tr><tr><td>10</td><td>1.0639</td><td>1.6484</td><td>1</td><td>1.0759</td><td>1.6526</td><td>0.7/0.3</td><td>1.0669</td><td>1.6451</td></tr></tbody></table>

pre-processed data, while our method can still be trained on data unevenly distributed on the timeline. We refer to data unevenly distributed in time as temporal discontinuities.

To verify the robustness of our model under temporal discontinuities, we conduct three groups of experiments by randomly masking 0%, 50% and 70% of the frames respectively. Accordingly, we also can easily adapt our loss function to:

$$ \min \frac{1}{T} \int_{t=1}^{T} \|\mathbf{v}^t - \hat{\mathbf{v}}^t\|_2^2 \approx \frac{1}{T} \sum_{t \in \Phi} \|\mathbf{v}^t - \hat{\mathbf{v}}^t\|_2^2 \quad (18) $$

where $\Phi$ is the set of the observed frames and $T = |\Phi|$.

The corresponding experimental results are presented in Tab. 2. For both velocity prediction errors (*Err*<sub>vel</sub>) and flow prediction errors (*Err*<sub>flow</sub>), the performance remains stable across varying masking levels, with only minor variations observed. Notably, even under the 70% masking scenario, the errors are comparable to those without masking. These findings highlight the model's capacity to maintain accuracy and adaptability despite significant temporal discontinuities, making it well-suited for real-world applications where data observations are sparse or unevenly distributed in time.

## 2.2. Robustness to Input Noise

CCTV camera data is by far the major data on high-density crowd particularly because of the uninvasiveness of the data collection process. However, these videos often suffer from low quality because they are captured from distant cameras and tend to be blurred, making optical flow estimation unreliable. Owing to the physics model in Crowd MPM, our model incorporates a self-correction strategy (P2G operation) that effectively handles spatially noisy.

Since the noise in optical flows is unknown, we conducted an experiment on *Drill2* to show that our method

<figure><img src="image_13.png" /><figcaption>Figure 10. a–e Five simple motion vector fields corresponding to c and d (where c and d denote the value of curl and divergence respectively)[66].(a) no rotation or dispersal/gathering, (b) the rotation is clockwise and there is no divergence, (c) the rotation is counterclockwise without any dispersal/gathering, (d) the motion is pure dispersal without any curl component, (e) the motion is pure gathering without any curl component, (f) the motion is relatively complex since it has rotation and convergence simultaneously.</figcaption></figure>

is robust to several noises with various shapes and magnitudes, e.g. Gaussian, uniform, and their mixtures, in Tab. 3. All results are similar, demonstrating the robustness of our model. To understand this, when noisy data are given, the P2G projects it onto the grid (i.e. equivalent to one step smoothing) then the mass and momentum are preserved the whole time, essentially filtering out high-frequency noises. During training, even when the ground-truth is noisy, Crowd MPM computes its closest physically sensible prediction, by conserving the mass and momentum. Admittedly, if the optical flow is completely obscure, e.g. a white noise, Crowd MPM would not learn the underlying dynamics. In this case, Crowd MPM can still be hand-tuned to visually mimic the data.

## 2.3. Visual Analysis Based on Learned Crowd Materials

In addition to prediction and simulation, our method is also a good tool for visual analysis of high-density crowds, especially in terms of its global flows. Previously, operators such as curl (vorticity) and divergence have been proposed as a good visualization tool for analysis [66]. However, [66] proposes to compute vorticity and divergence on optical flows while our method compute them on the velocity fields. Curl is defined as $c = \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}$, representing the tendency of the flow to exhibit rotational motion or vorticity. Divergence, on the other hand, is defined as $d = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}$, representing the rate at which a quantity (such as mass or velocity) expands or contracts at a point. Note that v is a motion vector field, which can either be optical flow or ve-
<figure><img src="image_14.png" /><figcaption>Figure 11. Qualitative comparison on curl and divergence map. iii and iv are generated from [66], while v and vi are obtained by our method. Zoom in for a better view.</figcaption></figure>
locity field. Several examples of global flow patterns are
provided in Fig. 10, illustrating how these flows can be an-
alyzed and interpreted using curl and divergence.

In our scenario, curl reflects steering behaviors, showing
how individuals adjust their movements, particularly at the
periphery or near obstacles or other people. Divergence in-
dicate areas where the crowd density is decreasing (positive
divergence) or increasing (negative divergence), such as the
region near an exit in Drill where people move towards. We
show a qualitative comparison on Drill in Fig. 11 at 3s, 9s
and 15s.

[66] is designed for crowd behavior classification, such
as lane formation, clockwise arch and so on, so it aims to
recognize motion patterns based on normalized velocities
(per pixel), shown in iii and iv in Fig. 11 (zoom in for a
better view). Therefore, the curl and divergence computed
based on normalized velocities can only indicate the general
patterns of motions, not the actual motions. Visually, the
global flow patterns are not obvious.

In contrast, our model can act as a more detailed ana-
lyzer based on the velocity and interpret how the behavior
evolves. For example, the areas pointed out by arrows in
Fig. 11a vi, Fig. 11b vi and Fig. 11c vi have relatively high
positive divergence. Physically, this means there are masses
flowing out of these areas. In time, these areas move toward
the exit. This captures that the fact that people come into
the scene and walk towards the exit. After a short while,
there is no people entering the scene. Correspondingly, the
high divergence areas first appear near the top entrance then
gradually move towards the exit. The divergence in front
of the exit is always negative, i.e. masses flowing in, due to
that for this area the exit rate is lower than the entering rate
(i.e. people aggregating in front of the exit but not many can
get through), hence the blue regions.

In terms of the curl, it reflects the local rotations of
masses in a region. This is reflected by the red and blue
areas near the exit at *t* = 3s in Fig. 11a v. This is when
the first batch of people arrive at the exit. At this time, peo-
ple can still go through the exit relatively easily, especially
the people right in front of the exit walking relatively faster
than the people on both sides of the exit. This caused a
relative local rotation of the mass flow, *i.e.* masses in the
middle flowing out faster than the ones on both sides, and
masses on both sides getting into the middle flow. Both the
clockwise (red region) and counter-clockwise (blue region)
vorticity reflect this. This curl becomes less prominent at
*t* = 9s and *t* = 12, when many people are crowded in front
of the exit and start to block the way, so that the general
movements become slow.

Fig. 11 demonstrates the importance of being able to
extract and stably simulate the velocity fields for analy-
sis. Capturing the underlying velocity field provides clear
global flow trends, then direct analysis on the optical flows

which might be polluted by noises [66]. Our method pro-
vides a tool which not only extracts and simulates the veloc-
ity field, but also learns the dynamics from specific crowds.

2.4. Ablation Study

We analyze the impact of different components, by evalu-
ating four variants of our model, including different combi-
nations of learnable parameters:

* Baseline_I, a weakly incompressible fluid model [48], with a global, non-learnable parameter $\epsilon$ whose value is obtained by a grid search.

* Baseline<sub>II</sub> is the same as Baseline<sub>I</sub> but with $\epsilon$ learnable for each particle.

* Baseline<sub>III</sub> has a new strain-stress tensor and new learnable parameter *k*, based on Baseline<sub>II</sub>.

* Baseline<sub>IV</sub> is Baseline<sub>III</sub> with $NN_\alpha$ (introduced in Sec.3.4 in main paper).

* *Ours* is the full model.

Table 4 clearly shows a continuously improvement across multiple metrics, when more components are added. It has been argued that high-density crowds behave like fluids [52]. This is confirmed by the hand-tuned *Baseline*I, whose performance is reasonable. *Baseline*_II shows the necessity of making the material learnable. With learnable $\epsilon$, the material becomes more heterogeneous and fits the data better. Furthermore, an addition of a learnable compression component (*Baseline*<sub>III</sub>) can capture the dynamics better. We notice that the difference between *Baseline*<sub>*I*−*III*</sub> seems small. After further investigation, this is likely to be caused by two factors. First, *Baseline*I is a already good baseline as argued by previous research. Hajj crowds form stable flows circling around the center of the scene. Therefore, it is possible to hand-tune *Baseline*I to mimic the general dynamics, although this relies on heavy human labor. Second, the metrics are based on the optical flow or the velocity field estimated from optical flow, where the scale of values is small.

The results are significantly after active forces are intro-
duced, even with only the motion alignment component in
*Baseline*<sub>IV</sub> with no stochasticity. This further proves such
crowds can be interpreted as active matters and learning
such active forces greatly enhance the model's ability to
predict dynamics. Further adding the stochastic part again
significantly improves the results (*Ours*).

3. Crowd Material Point Method

There are two mainstream methods to solve PDEs. One is
the Lagrangian method which treats the material as a set
of particles and tracks them to determine their properties
such as mass, positions, and velocities. The other one is the
Eulerian method which discretizes the space into grids con-
taining nodes with various properties and employs the grids
to represent materials. Material Point Method (MPM) is a
<table><caption>Table 4. Ablation study on *Drill*₁ and *Hajj*. Results involving active force are based on 10 trials. The training and testing duration are 60 frames on *Drill*₁ and 30 frames on *Hajj*.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Drill<sub>1</sub></th><th colspan="2">Hajj</th></tr><tr><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th><th>Err<sub>vel</sub></th><th>Err<sub>flow</sub></th></tr></thead><tbody><tr><td><b>Baseline<sub>I</sub></b></td><td>0.60740</td><td>0.98405</td><td>1.45878</td><td>1.68066</td></tr><tr><td><b>Baseline<sub>II</sub></b></td><td>0.60727</td><td>0.98396</td><td>1.45709</td><td>1.67930</td></tr><tr><td><b>Baseline<sub>III</sub></b></td><td>0.60483</td><td>0.98160</td><td>1.45679</td><td>1.67906</td></tr><tr><td><b>Baseline<sub>IV</sub></b></td><td>0.48987</td><td>0.88263</td><td>0.76161</td><td>1.07112</td></tr><tr><td><b>Ours(mean)</b></td><td>0.48276</td><td>0.87165</td><td>0.69135</td><td>0.99511</td></tr><tr><td><b>Ours(best)</b></td><td>0.48161</td><td>0.87033</td><td>0.68281</td><td>0.98604</td></tr></tbody></table>

hybrid Eulerian-Lagrangian method combining the advantages of Eulerian and Lagrangian methods. As our model involves Eulerian data and Lagrangian behaviors, we enhance the standard MPM to propose Crowd MPM to simulate extremely high-density crowds.

Our derivations are mainly following [27]. The main changes are brought by our newly introduced learnable stress and active forces. Similar to the standard MPM, our Crowd MPM also incorporates both the Lagrangian view and the Eulerian view. Under the Lagrangian view, we focus on the particle *p* with mass *m*<sub>*p*</sub>, position **x**<sub>*p*</sub>, and velocity **v**<sub>*p*</sub>. The motion of material is defined by a deformation map $\Phi(\cdot, t) : \Omega^0 \to \Omega^t$, where $\Omega^0, \Omega^t \subset \mathbb{R}^2$. $\Omega^t$ denotes the set of positions of particles at time *t*. As a result, $\Omega^0$ is the set of initial positions of particles. Particularly, we use **X**<sub>*p*</sub> and **x**<sub>*p*</sub> to denote the initial position and the current position for any particle in the simulated material, respectively. We can note that **X**<sub>*p*</sub> = **x**<sub>*p*</sub><sup>0</sup>. The deformation map $\Phi$ further enables us to determine the velocity (we ignore the subscript *p* for notation simplicity in this section):

$$ \mathbf{V}(\mathbf{X}, t) = \frac{\partial \Phi}{\partial t}(\mathbf{X}, t). \qquad (19) $$

The Jacobian of the deformation map $\Phi$ is significantly important and will be used often later:

$$ \mathbf{F}(\mathbf{X}, t) = \frac{\partial \Phi}{\partial \mathbf{X}}(\mathbf{X}, t). \qquad (20) $$

The most common way to employ **F** is to use its determinant *J*. Finally, we give the governing equations based on the conservation of mass and conservation of momentum:

$$ R(\mathbf{X}, t)J(\mathbf{X}, t) = R(\mathbf{X}, 0), \qquad (21) $$

$$ R(\mathbf{X}, 0) \frac{\partial \mathbf{V}}{\partial t} = \nabla^{\mathbf{X}} \cdot \mathbf{P}^{\mathbf{cm}} + R(\mathbf{X}, 0)(\mathbf{B}^{\mathbf{bd}} + \mathbf{B}^{\mathbf{act}}), \quad (22) $$

where *R* is the Lagrangian mass density, **P**<sup>cm</sup> is the crowd material stress under the Lagrangian view, and **B**<sup>bd</sup> and **B**<sup>act</sup> denote the accelerations deriving from the body force

and the active force on **X**, respectively. To be more specific, the Lagrangian mass density is defined via:

$$ R(\mathbf{X}, t) = \rho(\Phi(\mathbf{X}, t), t) = \rho(\mathbf{x}, t), \qquad (23) $$

$$ \rho(\mathbf{x}, t) = \lim_{\epsilon \to +0} \frac{\mathrm{mass}(B_{\epsilon}^{t})}{\int_{B_{\epsilon}^{t}} d\mathbf{x}}, \qquad (24) $$

where $B_{\epsilon}^{t} \subset \Omega^{t}$ is a ball with radius $\epsilon$ and center $x \in \Omega^{t}$. We note that the Lagrangian view builds governing equations on material space $\Omega^{0}$. We offer the proof for Eq. (21) and Eq. (22) in Sec. 4.

Under the Eulerian view, we pay attention to the world space $\Omega^{t}$. In addition, grids consisting of a series of nodes with index **i** are built to estimate material motion. Each note **i** has physics properties including mass $m_{\epsilon}$, position $\mathbf{x}_{\epsilon}$, and velocity $\mathbf{v}_{\epsilon}$, etc. We also give the governing equations based on the conservation of mass and conservation of momentum:

$$ \frac{D}{Dt}\rho(\mathbf{x}, t) + \rho(\mathbf{x}, t)\nabla^{\mathbf{x}} \cdot \mathbf{v}(\mathbf{x}, t) = 0 \qquad (25) $$

$$ \rho(\mathbf{x}, t) \frac{D\mathbf{v}}{Dt} = \nabla^{\mathbf{x}} \cdot \sigma^{\mathbf{cm}} + \rho(\mathbf{x}, t)(\mathbf{b}^{\mathbf{bd}} + \mathbf{b}^{\mathbf{act}}) \qquad (26) $$

where $\frac{D}{Dt}$ is the material derivative [50], $\mathbf{v}(\mathbf{x}, t) = \mathbf{V}(\Phi^{-1}(\mathbf{x}, t), t)$, $\sigma^{\mathbf{cm}}$ is the crowd material stress under the Eulerian view, and $\mathbf{b}^{\mathbf{bd}}$ and $\mathbf{b}^{\mathbf{act}}$ denote the accelerations deriving from the body force and the active force on **x**, respectively. We also leave the proof for Eq. (25) and Eq. (26) in Sec. 4.

Finally, we give the weak form, which is crucial for the discretization, of the conservation of momentum in the Eulerian view:

$$ \int_{\Omega^t} q_i \rho a_i d\mathbf{x} = \int_{\partial \Omega^t} q_i t_i ds(\mathbf{x}) - \int_{\Omega^t} q_{i,k} \sigma_{ik}^{cm} d\mathbf{x} + \int_{\Omega^t} q_i \rho b_i^{bd} d\mathbf{x} + \int_{\Omega^t} q_i \rho b_i^{act} d\mathbf{x}. \qquad (27) $$

where $q_i$ is the ith component of an arbitrary function $\mathbf{q}(\cdot, t) : \Omega^t \to \mathbb{R}^d$, $q_{i,k} = \frac{\partial q_i}{\partial x_k}$, $a_i$ is the ith component of $\mathbf{a} = \frac{D\mathbf{v}}{Dt}$, and $t_i$ is the ith component of the boundary force per unit reference area $\mathbf{t}(\mathbf{x}, t)$. Similarly, $\sigma_{ik}^{cm}$, $b_i^{bd}$ and $b_i^{act}$ are the components of $\sigma^{\mathbf{cm}}$, $\mathbf{b}^{\mathbf{bd}}$, and $\mathbf{b}^{\mathbf{act}}$, respectively. $ds$ denotes a tiny area. For efficient expression, the summation is implied on the repeated index. The proof of the weak form can be found in Sec. 4.

Subsequently, we can clarify three steps in our Crowd MPM individually: (1) Particle-to-grid transfer in Sec. 3.1 (2) Grid Operation in Sec. 3.2, and (3) Grid-to-particle transfer in Sec. 3.3.

### 3.1. Particle-to-grid Transfer

This step aims to transfer particle mass and momentum to the grid at each timestamp *n*. Each particle contributes to
the transfer process based on weights relying on particles and grids. Specifically, the transfer process in defined by:

$$m_i^n = \int_{\mathbf{x} \in \Omega^n} W_i^n(\mathbf{x}) \rho(\mathbf{x}) d\mathbf{x}, \quad (28)$$

$$m_i^n \mathbf{v}_i^n = \int_{\mathbf{x} \in \Omega^n} W_i^{n'}(\mathbf{x}) \rho(\mathbf{x}) \mathbf{v}(\mathbf{x}) d\mathbf{x}, \quad (29)$$

where $W_i^n$ and $W_i^{n'}$ are weight functions, $\Omega^n$ is the set of all particles at timestamp $n$, $\rho$ denotes density.

We follow the Affine Particle-In-Cell (APIC) method [26] to discretize Eq. (28) and Eq. (29) due to its excellent numerical properties:

$$m_i^n = \sum_p w_{ip}^n m_p, \quad (30)$$

$$m_i^n \mathbf{v}_i^n = \sum_p w_{ip}^n [m_p \mathbf{v}_p^n + m_p \mathbf{C}_p^n (\mathbf{x}_i - \mathbf{x}_p^n)], \quad (31)$$

where $w_{ip}^n = \phi(\mathbf{x}_i - \mathbf{x}_p^n)$ ($\phi$ is a quadratic B-spline function [25]), $\mathbf{C}_p^n$ is the affine velocity gradient [26]. $m_p$ and $\mathbf{x}_i$ don't have superscript $n$ because $m_p$ doesn't change with time and the grids are rebuilt at each MPM iteration.

## 3.2. Grid Operation

We update velocities on the grid by:

$$\mathbf{v}_i^{n+1} = \mathbf{v}_i^n + \Delta t \frac{\mathbf{f}_i^n}{m_i}, \quad (32)$$

Then boundary conditions are employed to refine the updated velocities:

$$\mathbf{v}_i^{n+1} = \mathbf{B}\mathbf{C}(\mathbf{v}_i^{n+1}) = \mathbf{v}_i^{n+1} - \gamma \mathbf{n} \langle \mathbf{n}, \mathbf{v}_i^{n+1} \rangle. \quad (33)$$

We explain how to obtain $\mathbf{v}_i^n$ and $m_i$ in the last section. Next, I will introduce how to estimate $\mathbf{f}_i^n$.

Like Sec. 3.1, we transfer forces on particles to the grid at each timestamp $n$. We give the transfer process for each component of forces. To avoid confusion, we use Greek letters like $\alpha$ and $\beta$ to denote the component index.

$$\mathbf{f}_{i\alpha}^n = \int_{\mathbf{x} \in \Omega^n} W_i^n(\mathbf{x}) \rho(\mathbf{x}, t) \mathbf{a}_\alpha(\mathbf{x}, t) d\mathbf{x} \quad (34)$$

where other components of $\mathbf{f}_i^n$ has a similar transfer process. According to the weak form of the conservation of momentum in the Eulerian view Eq. (27), we have:

$$\mathbf{f}_{i\alpha} = \int_{\Omega^t} W_i \rho a_\alpha d\mathbf{x} = \int_{\partial\Omega^t} W_i t_\alpha ds(\mathbf{x}) - \int_{\Omega^t} W_{i,\gamma} \sigma_{\alpha\gamma}^{cm} d\mathbf{x} \\ + \int_{\Omega^t} W_i \rho b_\alpha^{bd} d\mathbf{x} + \int_{\Omega^t} W_i \rho b_\alpha^{act} d\mathbf{x}. \quad (35)$$

where the superscript is ignored for notation simplicity. Here, we set $q_\alpha = W_i$ and $q_\gamma = 0$ because $\mathbf{q}$ is arbitrary. Following previous work [27], we omit the $\int_{\partial\Omega^t} W_i t_\alpha ds(\mathbf{x})$ in Eq. (35) for simplicity given its minor effect. Given our estimation of $\sigma_p^{cm} = \sigma(\mathbf{x}_p, t)$ at every Lagrangian particle $\mathbf{x}_p$, we have:

$$\begin{aligned} f_{i\alpha}^{st} &= - \int_{\Omega^t} W_{i,\gamma} \sigma_{\alpha\gamma}^{cm} d\mathbf{x} \approx - \sum_p \sigma_{p\alpha\gamma}^{cm} w_{ip,\gamma} V_p \\ &= - \sum_p \sigma_{p\alpha}^{cm} \nabla w_{ip} V_p, \end{aligned} \quad (36)$$

where $\sigma_{p\alpha}^{cm}$ is the $\alpha$th row of $\sigma_p^{cm}$, $V_p$ is the volume of $B_{\Delta x}^t$ at the Lagrangian particle $\mathbf{x}_p$. We follow [26] to estimate $V_p$. Because we use the quadratic B-spline function to represent $w_{ip}$, we can obtain its gradient by:

$$\nabla w_{ip} = \frac{4}{\Delta x^2} w_{ip} (\mathbf{x}_i - \mathbf{x}_p) \quad (37)$$

Combining Eq. (36) and Eq. (37) results in:

$$f_{i\alpha}^{st} \approx - \sum_p \frac{4}{\Delta x^2} w_{ip} \sigma_{p\alpha}^{cm} (\mathbf{x}_i - \mathbf{x}_p) V_p \quad (38)$$

By defining:

$$\frac{4}{\Delta x^2} \sigma_{p\alpha}^{cm} (\mathbf{x}_i - \mathbf{x}_p) V_p = \mathbf{G}_p (\mathbf{x}_i - \mathbf{x}_p) + \sum_{p' \in N_p} (\mathbf{f}_r^{pp'} - \mathbf{f}_t^{pp'}), \quad (39)$$

we can get the Eq.11 in the main paper. For the body force, we take the centripetal force as an example. Other body forces such as the goal attraction force have a similar derivation. Therefore, $\mathbf{b}^{bd}(\mathbf{x}, t) = \frac{\mathbf{v}(\mathbf{x}, t)^2}{r}$, where $r$ is the length of the radius of the moving circle. Then, we have:

$$\begin{aligned} f_{i\alpha}^{bd} &= \int_{\Omega^t} W_i \rho b_\alpha^{bd} d\mathbf{x} = \int_{\Omega^t} W_i \rho \frac{\mathbf{v}_\alpha^2}{r} d\mathbf{x} \\ &\approx \sum_p w_{ip} m_p \frac{\mathbf{v}_{p\alpha}^2}{r}. \end{aligned} \quad (40)$$

As for the active force, we have $\mathbf{b}^{act} = a\mathbf{v} + \mathbf{c}$, where $\mathbf{c}$ is the stochastic part of the active force, i.e. $\mathbf{c} = -b|\mathbf{v}|^2\mathbf{v} + D_L\nabla(\nabla \cdot \mathbf{v}) + D_1\nabla^2\mathbf{v} + D_2(\mathbf{v} \cdot \nabla)^2\mathbf{v} + \tilde{\mathbf{f}}$. $a$ and $\mathbf{c}$ are predicted by neural networks. Then, we have:

$$\begin{aligned} f_{i\alpha}^{act} &= \int_{\Omega^t} W_i \rho b_\alpha^{act} d\mathbf{x} = \int_{\Omega^t} W_i \rho (a\mathbf{v}_\alpha + c_\alpha) d\mathbf{x} \\ &= \sum_p w_{ip} m_p (a\mathbf{v}_{p\alpha} + c_{p\alpha}). \end{aligned} \quad (41)$$
Combining Eq. (38), Eq. (40) and Eq. (41), we have:

$$
\begin{equation}
\begin{split}
\mathbf{f}_{i\alpha} = f_{i\alpha}^{st} + f_{i\alpha}^{bd} + f_{i\alpha}^{act} &\approx \sum_{p} \frac{4}{\Delta x^2} w_{ip} \sigma_{p\alpha}^{cm} (\mathbf{x}_i - \mathbf{x}_p) V_p \\
&\quad + w_{ip} m_p \frac{\mathbf{v}_{p\alpha}^2}{r} + w_{ip} m_p (a \mathbf{v}_{p\alpha} + c_{p\alpha})
\end{split}
\tag{42}
\end{equation}
$$

**3.3. Grid-to particle Transfer**

Finally, we transfer the updated velocities $\mathbf{v}_i^{n+1}$ back to particles and update positions of these particles:

$$
\mathbf{v}_p^{n+1} = \sum_i w_{ip}^n \mathbf{v}_i^{n+1}, \quad \mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \mathbf{v}_p^{n+1}, \quad (43)
$$

We update the velocity gradient $C_p^{n+1}$ for next step:

$$
\mathbf{C}_p^{n+1} = \frac{4}{\Delta x^2} \sum_i w_{ip}^n \mathbf{v}_i^{n+1} (\mathbf{x}_i - \mathbf{x}_p^n)^T . \quad (44)
$$

Additionally, we also update the deformation gradient $\mathbf{F}_p^{n+1}$ for estimating $V_p$:

$$
\mathbf{F}_p^{n+1} = (\mathbf{I} + \Delta t \mathbf{C}_p^{n+1}) \mathbf{F}_p^n. \quad (45)
$$

**4. Proof Details**

**Push Forward and Pull Back.** We introduce the push forward and pull back operations based on the deformation map $\Phi$ for efficient formula derivation and understanding. The deformation map $\Phi: \Omega^0 \to \Omega^t$ typically is assumed to be bijective, which means that any function defined on one set like $\Omega^0$ can naturally be regarded as another function defined on another set like $\Omega^t$ by changing variables. Push forward and pull back are two ways of variable substitution here. To be formal, given the time $t$ and a function $G(\cdot, t) : \Omega^0 \to \mathbb{R}$, we define the push forward $g(\cdot, t) : \Omega^t \to \mathbb{R}$ as $g(\mathbf{x}, t) = G(\Phi^{-1}(\mathbf{x}, t), t)$ with $\mathbf{x} \in \Omega^t$. Similarly, given the time $t$ and a function $g(\cdot, t) : \Omega^t \to \mathbb{R}$, the pull back $G(\cdot, t) : \Omega^0 \to \mathbb{R}$ is defined as $G(\mathbf{X}, t) = g(\Phi(\mathbf{X}, t), t)$ with $\mathbf{X} \in \Omega^0$. Particularly, the pull back of the pull forward of a function such as $G(\cdot, t)$ is itself. We commonly think of a function of $\mathbf{X}/\mathbf{x}$ as Lagrangian/Eulerian. Therefore, the push forward (pull back) of a function which is Lagrangian (Eulerian) is Eulerian (Lagrangian).

**Conservation of Mass under the Lagrangian View.**
For an arbitrary **x** ∈ Ω<sup>t</sup>, we consider the ball B<sub>ε</sub><sup>t</sup> ⊂ Ω<sup>t</sup> with radius ε and center **x**. The mass in B<sub>ε</sub><sup>t</sup> shouldn't vary over time when ε is small enough. Therefore, we have mass(B<sub>ε</sub><sup>t</sup>) = mass(B<sub>ε</sub><sup>0</sup>). According to the definition of ρ, R, and J, we have:

$$
\text{mass}(B_{\epsilon}^{t}) = \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) d\mathbf{x} = \int_{B_{\epsilon}^{0}} R(\mathbf{x}, t) J d\mathbf{x}, \quad (46)
$$

$$
\mathrm{mass}(B_{\epsilon}^{0}) = \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) d\mathbf{X}, \quad (47)
$$

for all $B_\epsilon^t \subset \Omega^t$ ($t \ge 0$), where the second equality in Eq. (46) is obtained by the formula of changing variables. From Eq. (46) and Eq. (47), we can derive:

$$
\int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) J d\mathbf{X} = \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) d\mathbf{X}. \quad (48)
$$

Given that $B_{\epsilon}^{t}$ is arbitrary, we can obtain:

$$
R(\mathbf{X}, t)J(\mathbf{X}, t) = R(\mathbf{X}, 0), \quad (49)
$$

for any $\mathbf{X} \in \Omega^0$.

**Conservation of Mass under the Eulerian View.** We start from Eq. (49), i.e. Eq. (21), to derive the conservation of mass equation under the Eulerian view. First, we have:

$$
\frac{\partial}{\partial t} (R(\mathbf{X}, t) J(\mathbf{X}, t)) = \frac{\partial}{\partial t} R(\mathbf{X}, 0) = 0, \quad (50)
$$

according to Eq. (49). Then, it is easy to get:

$$
\frac{\partial}{\partial t}(RJ) = \frac{\partial R}{\partial t}J + R\frac{\partial J}{\partial t} = 0, \quad (51)
$$

where we neglect $(\mathbf{X}, t)$ for notation simplicity. Then, we use the result from [27]:

$$
\frac{\partial J}{\partial t} = J \frac{\partial v_1}{\partial x_1} + J \frac{\partial v_2}{\partial x_2}, \quad (52)
$$

where $\mathbf{x} = (x_1, x_2) \in \Omega^t$ and $\mathbf{v}(\mathbf{x}, t) = (v_1, v_2)$. Combining Eq. (51) and Eq. (52) results in:

$$
\frac{\partial R}{\partial t} J + R J \left( \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2} \right) = 0. \quad (53)
$$

We push forward on both sides of Eq. (53) to obtain:

$$
\frac{D}{Dt} \rho(\mathbf{x}, t) + \rho(\mathbf{x}, t) \nabla^{\mathbf{x}} \cdot \mathbf{v}(\mathbf{x}, t) = 0. \quad (54)
$$

**Conservation of Momentum under the Eulerian**
**View.** Given an arbitrary $B_{\epsilon}^{t} \subset \Omega^{t}$, the momentum change on $B_{\epsilon}^{t}$ can be expressed as:

$$
\begin{equation}
\begin{split}
\frac{d}{dt} \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) \mathbf{v}(\mathbf{x}, t) d\mathbf{x} &= \int_{\partial B_{\epsilon}^{t}} \sigma^{\mathrm{cm}} \mathbf{n} ds(\mathbf{x}) + \\
&\quad \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) (\mathbf{b}^{\mathrm{bd}} + \mathbf{b}^{\mathrm{act}}) d\mathbf{x},
\end{split}
\tag{55}
\end{equation}
$$

where $\mathbf{n}(\mathbf{x})$ is the unit outward normal of $\partial B_{\epsilon}^t$ at $\mathbf{x}$ and $ds$ denotes a tiny area. We further have:

$$
\begin{align}
& \frac{d}{dt} \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) \mathbf{v}(\mathbf{x}, t) d\mathbf{x} = \frac{d}{dt} \int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) \mathbf{V}(\mathbf{X}, t) J d\mathbf{X} \nonumber \\
& = \int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) J(\mathbf{X}, t) \mathbf{A}(\mathbf{X}, t) d\mathbf{X}, \tag{56}
\end{align}
$$
where $\mathbf{A}(\mathbf{X}, t) = \frac{\partial^2 \Phi}{\partial t^2}(\mathbf{X}, t) = \frac{\partial \mathbf{V}}{\partial t}(\mathbf{X}, t)$, the first equality comes from the formula of changing variables, and the second equality obtained by swapping the order of the taking derivative symbol and the integral symbol. Combining Eq. (55) and Eq. (56) results in:

$$ \int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) J(\mathbf{X}, t) \mathbf{A}(\mathbf{X}, t) d\mathbf{X} = \int_{\partial B_{\epsilon}^{t}} \sigma^{\text{cm}} \mathbf{n} ds(\mathbf{x}) + \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) (\mathbf{b}^{\text{bd}} + \mathbf{b}^{\text{act}}) d\mathbf{x} \quad (57) $$

We push forward the left side of Eq. (57) to get:

$$ \int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) J(\mathbf{X}, t) \mathbf{A}(\mathbf{X}, t) d\mathbf{X} = \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) \mathbf{a}(\mathbf{x}, t) d\mathbf{x}, \quad (58) $$

where $\mathbf{a}(\mathbf{x}, t) = \mathbf{A}(\Phi^{-1}(\mathbf{x}, t), t)$. Further, we also have $\mathbf{A}(\Phi^{-1}(\mathbf{x}, t), t) = \frac{D\mathbf{v}}{Dt}$, which has proof in [27]. Then, Eq. (57) becomes:

$$ \begin{aligned} \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) \mathbf{a}(\mathbf{x}, t) d\mathbf{x} &= \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) \frac{D\mathbf{v}}{Dt} d\mathbf{x} = \\ & \int_{\partial B_{\epsilon}^{t}} \sigma^{\text{cm}} \mathbf{n} ds(\mathbf{x}) + \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) (\mathbf{b}^{\text{bd}} + \mathbf{b}^{\text{act}}) d\mathbf{x} = \\ & \int_{B_{\epsilon}^{t}} \nabla^{\mathbf{x}} \cdot \sigma^{\text{cm}} d\mathbf{x} + \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) (\mathbf{b}^{\text{bd}} + \mathbf{b}^{\text{act}}) d\mathbf{x} \end{aligned} \quad (59) $$

Since $B_{\epsilon}^{t}$ is arbitrary, we have:

$$ \rho(\mathbf{x}, t) \frac{D\mathbf{v}}{Dt} = \nabla^{\mathbf{x}} \cdot \sigma^{\text{cm}} + \rho(\mathbf{x}, t) (\mathbf{b}^{\text{bd}} + \mathbf{b}^{\text{act}}). \quad (60) $$

**Conservation of Momentum under the Lagrangian View.** We can also derive the conservation of momentum under the Lagrangian view from Eq. (57). Instead of pushing forward, we pull back the right side of Eq. (57). The pull back of the first integral on the right side of Eq. (57) is:

$$ \begin{aligned} \int_{\partial B_{\epsilon}^{t}} \sigma^{\text{cm}}(\mathbf{x}, t) \mathbf{n} ds(\mathbf{x}) &= \\ \int_{\partial B_{\epsilon}^{0}} J(\mathbf{X}, t) \sigma^{\text{cm}}(\Phi(\mathbf{X}, t), t) \mathbf{F}^{-T}(\mathbf{X}, t) \mathbf{N}(\mathbf{X}) ds(\mathbf{X}), \end{aligned} \quad (61) $$

where $\mathbf{N}(\mathbf{X})$ is the unit outward normal of $\partial B_{\epsilon}^{0}$ at $\mathbf{X}$. We introduce our the Crowd first Piola Kirchoff stress $\mathbf{P}^{\text{cm}} = J\sigma^{\text{cm}}\mathbf{F}^{-T}$ and get:

$$ \begin{aligned} \int_{\partial B_{\epsilon}^{t}} \sigma^{\text{cm}}(\mathbf{x}, t) \mathbf{n} ds(\mathbf{x}) &= \int_{\partial B_{\epsilon}^{0}} \mathbf{P}^{\text{cm}}(\mathbf{X}, t) \mathbf{N} ds(\mathbf{X}) \\ &= \int_{B_{\epsilon}^{0}} \nabla^{\mathbf{x}} \cdot \mathbf{P}^{\text{cm}}(\mathbf{X}, t) d\mathbf{X}. \end{aligned} \quad (62) $$

Subsequently, we pull back the second integral on the right side of Eq. (57):

$$ \begin{aligned} \int_{B_{\epsilon}^{t}} \rho(\mathbf{x}, t) (\mathbf{b}^{\text{bd}} + \mathbf{b}^{\text{act}}) d\mathbf{x} &= \\ \int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) J(\mathbf{X}, t) (\mathbf{B}^{\text{bd}} + \mathbf{B}^{\text{act}}) d\mathbf{X} &= \\ \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) (\mathbf{B}^{\text{bd}} + \mathbf{B}^{\text{act}}) d\mathbf{X}, \end{aligned} \quad (63) $$

where we use Eq. (49) to get the second equality. We note that the left side of Eq. (57) can be written as using Eq. (49):

$$ \begin{aligned} \int_{B_{\epsilon}^{0}} R(\mathbf{X}, t) J(\mathbf{X}, t) \mathbf{A}(\mathbf{X}, t) d\mathbf{X} &= \\ \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) \mathbf{A}(\mathbf{X}, t) d\mathbf{X} &= \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) \frac{\partial \mathbf{V}}{\partial t}(\mathbf{X}, t) d\mathbf{X}. \end{aligned} \quad (64) $$

Combining Eq. (57), Eq. (62), Eq. (63), and Eq. (64) results in:

$$ \begin{aligned} \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) \frac{\partial \mathbf{V}}{\partial t}(\mathbf{X}, t) d\mathbf{X} &= \\ \int_{B_{\epsilon}^{0}} \nabla^{\mathbf{x}} \cdot \mathbf{P}^{\text{cm}}(\mathbf{X}, t) d\mathbf{X} + \int_{B_{\epsilon}^{0}} R(\mathbf{X}, 0) (\mathbf{B}^{\text{bd}} + \mathbf{B}^{\text{act}}) d\mathbf{X}. \end{aligned} \quad (65) $$

Because $B_{\epsilon}^{0}$ is arbitrary, we can derive:

$$ R(\mathbf{X}, 0) \frac{\partial \mathbf{V}}{\partial t} = \nabla^{\mathbf{x}} \cdot \mathbf{P}^{\text{cm}} + R(\mathbf{X}, 0) (\mathbf{B}^{\text{bd}} + \mathbf{B}^{\text{act}}). \quad (66) $$

**Weak Form of the Conservation of Momentum in Both Views.** The weak form of an equation means that we can derive its weak form from the equation but can't derive the equation from its weak form. The weak forms of the conservation of momentum in both views are extremely important to derive our discretized Crowd MPM. We start with the Lagrangian view. We ignore the body force and the active force for simplicity:

$$ R(\mathbf{X}, 0) \mathbf{A}(\mathbf{X}, t) = \nabla^{\mathbf{x}} \cdot \mathbf{P}^{\text{cm}}. \quad (67) $$

Then, we have:

$$ R_0 A_i = \sum_j \frac{\partial P_{ij}^{cm}}{\partial X_j} = \sum_j P_{ij,j}^{cm}, \quad (68) $$

where $R_0 = R(\mathbf{X}, 0)$, $A_i$ is the ith component of $\mathbf{A}(\mathbf{X}, t)$, $P_{ij}^{cm}$ is the element at the ith row and jth column in $\mathbf{P}^{\text{cm}}$. For efficient expression, the summation is implied on the repeated index. Therefore, we have:

$$ R_0 A_i = P_{ij,j}^{cm}. \quad (69) $$
Then, we derive the weak form of Eq. (67) by computing the dot product between an arbitrary function $Q(\cdot, t) : \Omega^0 \to \mathbb{R}^d$ and two sides of Eq. (67), respectively, and integrate over $\Omega^0$:

$$ 
\begin{aligned}
\int_{\Omega^0} Q_i(\mathbf{X}, t) R(\mathbf{X}, 0) A_i(\mathbf{X}, t) d\mathbf{X} &= \\
\int_{\Omega^0} Q_i(\mathbf{X}, t) P_{ij,j}^{cm}(\mathbf{X}, t) d\mathbf{X} &= \\
\int_{\Omega^0} ((Q_i(\mathbf{X}, t) P_{ij}^{cm}(\mathbf{X}, t)), j - Q_{i,j}(\mathbf{X}, t) P_{ij}^{cm}(\mathbf{X}, t)) d\mathbf{X} &= \\
\int_{\partial\Omega^0} Q_i P_{ij}^{cm} N_j ds(\mathbf{X}) - \int_{\Omega^0} Q_{i,j} P_{ij}^{cm} d\mathbf{X}.
\end{aligned} 
\quad (70) $$

The $P_{ij}^{cm} N_j$ would be determined by a boundary condition. Assuming that $T(\mathbf{X}, t)$ is the boundary force per unit reference area, we have $T_i = P_{ij}^{cm} N_j$. As a result, we have that for an arbitrary $Q(\cdot, t) : \Omega^0 \to \mathbb{R}^d$:

$$ \int_{\Omega^0} Q_i R_0 A_i d\mathbf{X} = \int_{\partial\Omega^0} Q_i T_i ds(\mathbf{X}) - \int_{\Omega^0} Q_{i,j} P_{ij}^{cm} dX. \quad (71) $$

Considering the full Eq. (66), i.e. Eq. (22), we have the final weak form of the conservation of momentum in the Lagrangian view:

$$ 
\begin{aligned}
\int_{\Omega^0} Q_i R_0 A_i d\mathbf{X} = \int_{\partial\Omega^0} Q_i T_i ds(\mathbf{X}) - \int_{\Omega^0} Q_{i,j} P_{ij}^{cm} d\mathbf{X} \\
+ \int_{\Omega^0} Q_i R_0 B_i^{bd} d\mathbf{X} + \int_{\Omega^0} Q_i R_0 B_i^{act} d\mathbf{X}.
\end{aligned} 
\quad (72) $$

Then we push forward Eq. (72) to get the weak form of the conservation of momentum in the Eulerian view. The push forward of the left side of Eq. (72) is derived by introducing $\mathbf{q}$ that is the push forward of $\mathbf{Q}$:

$$ \int_{\Omega^0} Q_i R_0 A_i d\mathbf{X} = \int_{\Omega^t} q_i(\mathbf{x}, t) \rho(\mathbf{x}, t) a_i(\mathbf{x}, t) d\mathbf{x}. \quad (73) $$

Similarly, let $\mathbf{t}$ be the push forward of $\mathbf{T}$, we have:

$$ \int_{\partial\Omega^0} Q_i T_i ds(\mathbf{X}) = \int_{\partial\Omega^t} q_i(\mathbf{x}, t) t_i(\mathbf{x}, t) ds(\mathbf{x}). \quad (74) $$

Before giving the push forward of the second integral on the right side of Eq. (72), we analyze $Q_{i,j}$:

$$ Q_{i,j} = \frac{\partial Q_i}{\partial X_j} = \frac{\partial q_i}{\partial x_k} \frac{\partial x_k}{\partial X_j} = q_{i,k} F_{kj}, \quad (75) $$

where the summation is implied on the repeated index $k$, $F_{kj}$ is the element in $\mathbf{F}$ defined by Eq. (20), and $\mathbf{x} = \Phi(\mathbf{X}, t)$. Recall that $\mathbf{P}^{cm} = J\sigma^{cm}\mathbf{F}^{-T}$, then we have

$\sigma^{cm} = \frac{1}{J} \mathbf{P}^{cm} \mathbf{F}^T$. Therefore, we the relationship between $\sigma^{cm}$ and $\mathbf{P}^{cm}$:

$$ \sigma_{ik}^{cm} = \frac{1}{J} P_{ij}^{cm} F_{kj}. \quad (76) $$

Combining Eq. (75) and Eq. (76) let us derive the push forward of the second integral on the right side of Eq. (72):

$$ \int_{\Omega^0} Q_{i,j} P_{ij}^{cm} d\mathbf{X} = \int_{\Omega^t} q_{i,k}(\mathbf{x}, t) \sigma_{ik}^{cm}(\mathbf{x}, t) d\mathbf{x}. \quad (77) $$

As $\mathbf{b}^{bd}$ and $\mathbf{b}^{act}$ are the natural push forward of $\mathbf{B}^{bd}$ and $\mathbf{B}^{act}$, it is easy to derive the push forward of the remaining integrals in Eq. (72):

$$ 
\begin{aligned}
\int_{\Omega^0} Q_i R_0 B_i^{bd} d\mathbf{X} &= \int_{\Omega^t} q_i(\mathbf{x}, t) \rho(\mathbf{x}, t) b_i^{bd}(\mathbf{x}, t) d\mathbf{x} \\
\int_{\Omega^0} Q_i R_0 B_i^{act} d\mathbf{X} &= \int_{\Omega^t} q_i(\mathbf{x}, t) \rho(\mathbf{x}, t) b_i^{act}(\mathbf{x}, t) d\mathbf{x}.
\end{aligned} 
\quad (78) $$

Combining Eq. (73), Eq. (74), Eq. (77), and Eq. (78) results in the final weak form of the conservation of momentum in the Eulerian view:

$$ 
\begin{aligned}
\int_{\Omega^t} q_i \rho a_i d\mathbf{x} = \int_{\partial\Omega^t} q_i t_i ds(\mathbf{x}) - \int_{\Omega^t} q_{i,k} \sigma_{ik}^{cm} d\mathbf{x} \\
+ \int_{\Omega^t} q_i \rho b_i^{bd} d\mathbf{x} + \int_{\Omega^t} q_i \rho b_i^{act} d\mathbf{x}.
\end{aligned} 
\quad (79) $$