Hydrodynamic Cucker-Smale model with time delay and obstacle avoidance

We study a hydrodynamic Cucker-Smale-type model incorporating both time delays and obstacle potentials. The model governs the evolution of velocity and density fields of the system, where delayed interactions drive alignment and obstacle potentials account for responses to obstacles or predators. We further extend the framework to two-species systems. To numerically solve the model, we design a high-order finite volume method based on a Lax–Friedrichs numerical flux with fifth-order weighted essentially non-oscillatory reconstruction and third-order Runge–Kutta time discretization, ensuring numerical stability and high-order accuracy. Numerical experiments confirm the stability and accuracy of the proposed scheme and illustrate how time delays and obstacle potentials, under specific communication kernels and initial conditions, affect the emergence of flocking or non-flocking behavior.

1 Introduction

Flocking refers to the spontaneous emergence of coordinated collective motion among self-propelled agents, driven by local alignment of velocities and spatial cohesion [1,2]. This phenomenon is widely observed in nature, including bird flocks, fish schools, and insect swarms [3–6]. It also inspires algorithmic design in robotics, autonomous vehicles, and swarm control systems [7,8].

Classical models capturing flocking dynamics often use a particle-based framework. Reynolds’ Boids model [9] introduced simple behavioral rules (alignment, cohesion, separation) to generate realistic group motion. The Vicsek model [8] further simplified this by using noisy velocity alignment, revealing phase transitions between ordered and disordered states. The Cucker–Smale (CS) model [10] formalized alignment interactions using nonlocal communication kernels, laying the foundation for rich mathematical theory, including kinetic and hydrodynamic continuum descriptions [11–15].

While many Cucker–Smale (CS) type models assume instantaneous interactions among agents, real-world agents typically respond to information from their surroundings with a certain processing or reaction delay $\tau >0$ [16]. Recent studies have increasingly focused on the influence of delays, with the majority of them conducted within the framework of particle-based models [16–20]. For instance, Du [21] analyzed the delayed CS system in a harmonic potential and showed that delay affects the conditions for flocking. Choi and Haskovec [18] established sufficient conditions for global regularity and flocking, where for non-integrable kernels the delay can be arbitrary, whereas for integrable kernels the results generally require small delays [18]. shows that sufficiently small delays are necessary to satisfy the flocking condition. While delayed CS models have been studied at the hydrodynamic level [22], most results remain purely analytical. In contrast, existing simulations primarily focus on the particle level [19,20], but numerical investigations at the continuum scale remain limited.

In addition to delayed alignment, repulsive potentials are key to realistic swarm models, maintaining spacing, avoiding obstacles, and simulating predator evasion [23–30]. Empirical studies show that fish schools and bird flocks react to static obstacles by deflecting, slowing down, and regrouping, while preserving cohesion and direction [28,30]. While previous studies have analyzed obstacle-induced behaviors such as group splitting and path deviation in swarm models [26,29], the impact of obstacles in hydrodynamic Cucker–Smale systems remains less understood. Aung et al. [31] showed that obstacles and spatial heterogeneity can enhance local interactions and global order. In this work, we introduce stationary obstacles as localized repulsive potentials and focus on their interaction with time delay in shaping aggregation and pattern formation.

Incorporating both time delay and obstacle-induced repulsion potentials into the hydrodynamic Cucker–Smale (CS) model significantly increases the mathematical and numerical challenges. Specifically, such coupled models may develop singularities such as finite-time density concentration, blow-up of velocity gradients, or loss of regularity in the solution. Although particle-based simulations can capture some detailed behaviors, they typically suffer from high computational costs and limited scalability in high-dimensional or large-scale systems [32]. Therefore, continuum or hydrodynamic modeling offers a more feasible and efficient framework for simulating large-scale collective behavior. Several works have proposed high-order numerical methods for hydrodynamic CS models in the absence of time delay and obstacle effects [24,33,34]. However, numerical methods that simultaneously address both time-delayed interactions and obstacle-induced repulsionpotentials remain largely unexplored.

This work aims to contribute to this gap by studying hydrodynamic Cucker–Smale models that combine time delay and obstacle effects through the numerical simulation. Compared with existing studies that often consider these factors separately, we provide a numerical investigation of their joint influence and furthermore explore the solution behavior of the two-species scenario. In this work, we design a high-order finite volume scheme to simulate hydrodynamic Cucker–Smale models incorporating both time delays and obstacle-induced repulsive potentials. The method combines Lax–Friedrichs flux, fifth-order weighted essentially non-oscillatory (WENO) reconstruction, and third-order Runge–Kutta time discretization. Nonlocal alignment terms are computed efficiently by using the fast Fourier transforms (FFT). The proposed framework provides a numerical framework to study delayed and obstacle-induced potential hydrodynamic flocking models. We apply the scheme to one-dimensional single- and two-species systems, as well as a two-dimensional single-species setting. The simulations highlight how delay and obstacle interactions affect flocking, avoidance, and the onset of singular behavior.

Numerical experiments demonstrate that time delay can suppress flocking and may induce finite-time singularities, such as divergence in velocity gradients or blow-up of density. Repulsive interactions promote obstacle avoidance but tend to disrupt alignment, often contributing to singular behavior. Moreover, the presence of fixed obstacles combined with delay accelerates singularity onset and spatially shifts aggregation closer to the obstacle. These findings reveal the intricate, nonlinear interactions between delay and repulsion in collective behavior. These results illustrate the nonlinear interaction of delay and obstacle effects in shaping collective behavior.

The remainder of the paper is organized as follows. Section 2 introduces the hydrodynamic Cucker–Smale model with time delays and repulsive obstacle potentials, including its extension to two-species systems. Section 3 describes our high-order finite volume scheme, detailing the WENO reconstruction, Runge–Kutta time discretization, and FFT acceleration. Section 4 presents numerical simulations illustrating the effects of delays and repulsive potentials on flocking, obstacle avoidance, and singularity formation. Section 5 concludes the paper.

2 The model

In this section, we introduce hydrodynamic Cucker–Smale models incorporating both time delays and repulsive obstacle potentials, formulated for both single-species and two-species systems.

2.1 Single-species model with delay and obstacle

The single-species hydrodynamic Cucker–Smale system with time delay and obstacle-induced repulsionpotential is given by

$\left\{\begin{array}{cc}\hfill & {\partial }_t\rho +\nabla \cdot \left(\rho \mathbf{u}\right)=0,\hfill \\ \hfill & {\partial }_t\left(\rho \mathbf{u}\right)+\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}\right)={\int }_{{\mathbb{R}}^d}\varphi \left(|\mathbf{x}-\mathbf{y}|\right){\rho }_{\tau }\left(\mathbf{y}\right)\rho \left(\mathbf{x}\right)\left({\mathbf{u}}_{\tau }\left(\mathbf{y}\right)-\mathbf{u}\left(\mathbf{x}\right)\right)\mathrm{d}\mathbf{y}+\rho \left(\mathbf{x}\right)\mathbf{R}\left(\mathbf{x}\right),\hfill \end{array}\right.(1)$

for $t\ge 0$ and $\mathbf{x}\in {\mathbb{R}}^d$, where $d$ denotes the spatial dimension. The functions $\rho (\mathbf{x},t)$ and $\mathbf{u}(\mathbf{x},t)$ represent the density and velocity fields, respectively. The delayed quantities are defined as ${\rho }_{\tau }(\mathbf{y})≔\rho (\mathbf{y},t-\tau )$ and ${\mathbf{u}}_{\tau }(\mathbf{y})≔\mathbf{u}(\mathbf{y},t-\tau )$, where $\tau \ge 0$ is a constant denoting the time delay in communication and information processing. The communication kernel $\varphi (|\mathbf{x}-\mathbf{y}|)$ is a smooth, radially symmetric function encoding the strength of velocity alignment between agents. The term $\mathbf{R}(\mathbf{x})$ models a repulsive force field generated by static obstacles in the environment.

In this work, we adopt the fat-tailed communication kernel

$\varphi \left(r\right)=\frac{\delta }{{\left(1+r^2\right)}^{\theta }},(2)$

where $\theta >0$ controls the decay rate of interactions with respect to distance, and the constant $\delta$ is chosen such that ${\int }_0^{\infty }\varphi (r)dr=1$, and we model the obstacle force $\mathbf{R}(\mathbf{x})$ as the negative gradient of a localized Gaussian-type repulsive potential:

$\mathbf{R}\left(\mathbf{x}\right)=-\nabla \mathrm{\Phi }\left(\mathbf{x}\right),\text{with}\mathrm{\Phi }\left(\mathbf{x}\right)={\eta }_0\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}_0{|}^2}{l_0^2}\right),(3)$

where ${\eta }_0$ and $l_0$ control the strength and spatial scale of the repulsion, respectively, and ${\mathbf{x}}_0$ denotes the center of the obstacle.

While these specific forms are chosen in this study, alternative kernels and potentials (see Remark 2.1) may be considered depending on the modeling context or application.

Remark 2.1. (Alternative communication kernels and obstacle potentials). For communication kernels, besides fat-tailed ones, compactly supported kernels such as

$\varphi \left(r\right)={\left(1+r\right)}^2{\left(1-r\right)}^2{\chi }_{\left[\mathrm{0,1}\right]}\left(r\right),$

the singular communication kernels like $\varphi (r)=r^{-\alpha }$, $\alpha >0$, and the exponentially decaying kernels have been widely used [1,34,35].

For obstacle potentials, Gaussian repulsion is one of several possible choices. A widely used alternative is the attractive–repulsive potential function

$\mathrm{\Phi }\left(x\right)=\frac{|x{|}^a}{a}-\frac{|x{|}^b}{b},$

proposed in [24,36], offers a flexible shape controlling repulsion strength and range, and is often used in the study of obstacle problems. Attractive–repulsive hydrodynamics for collective consensus. These alternative choices affect the dynamics and numerical methods, and should be selected according to the specific application context.

2.2 Two-species model with delay and obstacle

We extend the model to a two-species system, where each species experiences both intra-species and inter-species alignment, possibly with different time delays.

Let ${\rho }_{\alpha }(\mathbf{x},t)$ and ${\mathbf{u}}_{\alpha }(\mathbf{x},t)$ denote the density and velocity fields of species $\alpha$. The governing equations are:

$\left\{\begin{array}{cc}\hfill & {\partial }_t{\rho }_{\alpha }+\nabla \cdot \left({\rho }_{\alpha }{\mathbf{u}}_{\alpha }\right)=0,\hfill \\ \hfill & {\partial }_t\left({\rho }_{\alpha }{\mathbf{u}}_{\alpha }\right)+\nabla \cdot \left({\rho }_{\alpha }{\mathbf{u}}_{\alpha }\otimes {\mathbf{u}}_{\alpha }\right)={\mathbf{S}}_{\alpha }\left(\mathbf{x},t\right),\hfill \end{array}\right.(4)$

where the source term ${\mathbf{S}}_{\alpha }(\mathbf{x},t)$ captures the nonlocal interactions and external obstacle force, given by

${\mathbf{S}}_{\alpha }\left(\mathbf{x},t\right)≔{\displaystyle \sum _{\beta =1}^2}{\int }_{{\mathbb{R}}^d}{\varphi }_{\alpha \beta }\left(|\mathbf{x}-\mathbf{y}|\right){\rho }_{\alpha }\left(\mathbf{x}\right){\rho }_{\beta ,{\tau }_{\beta }}\left(\mathbf{y}\right)\left({\mathbf{u}}_{\beta ,{\tau }_{\beta }}\left(\mathbf{y}\right)-{\mathbf{u}}_{\alpha }\left(\mathbf{x}\right)\right)\mathrm{d}\mathbf{y}+{\rho }_{\alpha }\left(\mathbf{x}\right){\mathbf{R}}_{\alpha }\left(\mathbf{x}\right),$

where $\alpha ,\beta \in \{\mathrm{1,2}\}$, ${\rho }_{\alpha }(\mathbf{x})≔{\rho }_{\alpha }(\mathbf{x},t)$ and ${\mathbf{u}}_{\alpha }(\mathbf{x})≔{\mathbf{u}}_{\alpha }(\mathbf{x},t)$ denote the density and velocity of species $\alpha$, respectively. The delayed terms are given by ${\rho }_{\beta ,{\tau }_{\beta }}(\mathbf{y})≔{\rho }_{\beta }(\mathbf{y},t-{\tau }_{\beta })$ and ${\mathbf{u}}_{\beta ,{\tau }_{\beta }}(\mathbf{y})≔{\mathbf{u}}_{\beta }(\mathbf{y},t-{\tau }_{\beta })$, with ${\tau }_{\beta }>0$ the communication delay of species $\beta$. The kernel ${\varphi }_{\alpha \beta }(|\mathbf{x}-\mathbf{y}|)$ characterizes the influence of species $\beta$ on $\alpha$: ${\varphi }_{\alpha \alpha }$ and ${\varphi }_{\beta \beta }$ describe intra-species interactions, while ${\varphi }_{\alpha \beta }$ and ${\varphi }_{\beta \alpha }$ $(\alpha \ne \beta )$ correspond to inter-species alignment. Existing theoretical analyses have mostly focused on the symmetric interaction case $({\varphi }_{\alpha \beta }={\varphi }_{\beta \alpha })$. In our numerical simulations, we will adopt the symmetric interactions ${\varphi }_{\alpha \beta }(r)={\varphi }_{\beta \alpha }(r)$ in all examples except Example 4.8. The term ${\mathbf{R}}_{\alpha }(\mathbf{x})$ accounts for local repulsive effects from external obstacles. As in the single-species case, it is modeled by the gradient of a repulsive potential:

${\mathbf{R}}_{\alpha }\left(\mathbf{x}\right)=-\nabla {\mathrm{\Phi }}_{\alpha }\left(\mathbf{x}\right),{\mathrm{\Phi }}_{\alpha }\left(\mathbf{x}\right)={\eta }_{\alpha }\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}_{\alpha }{|}^2}{l_{\alpha }^2}\right),(5)$

where ${\mathbf{x}}_{\alpha }$ is the center of the obstacle, and ${\eta }_{\alpha },l_{\alpha }$ control the strength and range of repulsion.

Remark 2.2. (Conservation Properties). The model preserves the total mass of each species under suitable boundary conditions. Define the total mass of species $\alpha$ as

$\mathbf{M}\left(t\right)≔{\displaystyle \sum _{\alpha \in \mathcal{I}}}{\mathbf{M}}_{\alpha }\left(t\right),\text{with}{\mathbf{M}}_{\alpha }\left(t\right)≔{\int }_{{\mathbb{R}}^d}{\rho }_{\alpha }\left(\mathbf{x},t\right)d\mathbf{x},\alpha \in \mathcal{I},$

where $\mathcal{I}=\{1\}$ for a single-species system and $\mathcal{I}=\{\mathrm{1,2}\}$ for a two-species system. Then each ${\mathbf{M}}_{\alpha }(t)$ is conserved over time:

$\frac{d}{dt}{\mathbf{M}}_{\alpha }\left(t\right)=0,\forall t>0,\forall \alpha \in \mathcal{I}.$

Define the total momentum as

$\mathbf{P}\left(t\right)≔{\displaystyle \sum _{\alpha \in \mathcal{I}}}{\mathbf{P}}_{\alpha }\left(t\right),\text{with}{\mathbf{P}}_{\alpha }\left(t\right)≔{\int }_{{\mathbb{R}}^d}{\rho }_{\alpha }\left(\mathbf{x},t\right){\mathbf{u}}_{\alpha }\left(\mathbf{x},t\right)d\mathbf{x}.$

In the absence of communication delays and obstacles, the symmetry of the interaction kernels ${\varphi }_{\alpha \beta }$ guarantees conservation of total momentum:

$\frac{d}{dt}\mathbf{P}\left(t\right)=0,\forall t>0.$

3 Numerical scheme

In this section, we develop a numerical scheme for the hydrodynamic model with time delays and repulsive obstacle potentials. To accurately capture sharp solution features while minimizing numerical dissipation, we adopt a fifth-order WENO reconstruction combined with a local Lax–Friedrichs numerical flux. Time discretization is performed using a third-order Runge–Kutta method, which ensures stability and maintains high-order accuracy. To efficiently compute the nonlocal alignment term, we employ the FFT-based convolution, significantly reducing computational cost.

3.1 Finite volume method for one-dimensional case

Consider a one-dimensional periodic domain $\mathrm{\Omega }=[a,b]$, divided into $N$ uniform cells of width $\mathrm{\Delta }x=(b-a)/N$. Let $x_j=a+(j-1/2)\mathrm{\Delta }x$ denote the center of cell $j$ for $j=1,\dots ,N$, and define the corresponding cell interval by $I_j=\left[x_{j-1/2},x_{j+1/2}\right],\text{where\hspace{0.17em}}{x}_{j\pm 1/2}=x_j\pm \frac{\mathrm{\Delta }x}{2}.$

${\rho }_j\left(t\right)\approx \frac{1}{\mathrm{\Delta }x}{\int }_{I_j}\rho \left(x,t\right)dx,m_j\left(t\right)\approx \frac{1}{\mathrm{\Delta }x}{\int }_{I_j}\rho \left(x,t\right)u\left(x,t\right)dx,$

with the velocity defined by $u_j(t)=m_j(t)/{\rho }_j(t)$. We denote the cell-averaged variables by $U_j(t)=\left(\begin{array}{c}\hfill {\rho }_j(t)\hfill \\ \hfill m_j(t)\hfill \end{array}\right)$. The semi-discrete finite volume formulation reads:

$\frac{d{U}_j}{dt}=-\frac{1}{\mathrm{\Delta }x}\left({\mathbf{F}}_{j+1/2}-{\mathbf{F}}_{j-1/2}\right)+\left(\begin{array}{c}\hfill 0\hfill \\ \hfill S_j\hfill \end{array}\right),$

where ${\mathbf{F}}_{j\pm 1/2}$ are numerical fluxes, and $S_j$ is the source term incorporating nonlocal alignment and obstacle effects. For convenience, we define the spatial operator $\mathcal{L}(U_j)$ as

$\mathcal{L}\left(U_j\right)≔-\frac{1}{\mathrm{\Delta }x}\left({\mathbf{F}}_{j+\frac{1}{2}}-{\mathbf{F}}_{j-\frac{1}{2}}\right)+\left(\begin{array}{c}\hfill 0\hfill \\ \hfill S_j\hfill \end{array}\right).(6)$

3.1.1 Lax–Friedrichs flux with fifth-order WENO reconstruction

In our scheme, we consider the pressureless Euler-type flux function:

$\mathbf{f}\left(U\right)=\left(\begin{array}{c}\hfill \rho u\hfill \\ \hfill \rho u^2\hfill \end{array}\right),U=\left(\begin{array}{c}\hfill \rho \hfill \\ \hfill m\hfill \end{array}\right),u=\frac{m}{\rho }.$

The Lax–Friedrichs numerical flux at the cell interface $x_{j+\frac{1}{2}}$ is defined by

${\mathbf{F}}_{j+\frac{1}{2}}=\frac{1}{2}\left(\mathbf{f}\left(U_{j+\frac{1}{2}}^{-}\right)+\mathbf{f}\left(U_{j+\frac{1}{2}}^{+}\right)\right)-\frac{{\lambda }_{j+\frac{1}{2}}}{2}\left(U_{j+\frac{1}{2}}^{+}-U_{j+\frac{1}{2}}^{-}\right),(7)$

where $U_{j+\frac{1}{2}}^{-}$ and $U_{j+\frac{1}{2}}^{+}$ are the fifth-order WENO reconstructed values from the left and right, respectively, and ${\lambda }_{j+\frac{1}{2}}$ is the maximum absolute eigenvalue of the Jacobian at the interface, defined by

${\lambda }_{j+\frac{1}{2}}=\underset{U\in \left\{U_{j+\frac{1}{2}}^{-},U_{j+\frac{1}{2}}^{+}\right\}}{\max }\left|eig\left(\frac{\partial \mathbf{f}}{\partial U}\right)\right|=\underset{U\in \left\{U_{j+\frac{1}{2}}^{-},U_{j+\frac{1}{2}}^{+}\right\}}{\max }|u|,(8)$

where $\frac{\partial \mathbf{f}}{\partial U}$ denotes the Jacobian matrix of the flux function $\mathbf{f}(U)$ with respect to the conserved variables $U$, evaluated at the intercell state $U_{j+\frac{1}{2}}$, and $eig(\cdot )$ denotes the eigenvalue of a matrix.

Remark 3.1. The wave speed ${\lambda }_{j+\frac{1}{2}}$ in Equation 7 can be chosen locally, as in Equation 8, to reduce numerical diffusion and capture local flow features, or globally as $\lambda ={\max }_j|u_j|$, which is more diffusive but more robust near strong gradients or discontinuities. Unless stated otherwise, we use the local wave speed to balance accuracy and stability.

To construct a high-order accurate scheme, we apply the fifth-order WENO scheme component-wise to each conservative variable in $U$ using a five-point stencil $\{U_{j-2},U_{j-1},U_j,U_{j+1},U_{j+2}\}$, which forms three three-point substencils:

$\left\{U_{j-2},U_{j-1},U_j\right\},\left\{U_{j-1},U_j,U_{j+1}\right\},\left\{U_j,U_{j+1},U_{j+2}\right\}.$

Each substencil produces a third-order approximation $q_k=p_k(x_{j+1/2})$. The final reconstructed value at the interface is a weighted sum:

$U_{j+1/2}^{-}={\displaystyle \sum _{k=0}^2}{\omega }_k{q}_k,$

where the nonlinear weights ${\omega }_k$ are computed by

${\omega }_k=\frac{{\tilde{\omega }}_k}{{\sum }_{l=0}^2{\tilde{\omega }}_l},{\tilde{\omega }}_k=\frac{d_k}{{\left(\epsilon +{\beta }_k\right)}^2}.$

Here, $d_k$ are linear weights, $\epsilon$ is a small positive parameter, and ${\beta }_k$ are smoothness indicators as defined in [37]:

$\begin{array}{cc}\hfill {\beta }_0& =\frac{13}{12}{\left(U_{j-2}-2U_{j-1}+U_j\right)}^2+\frac{1}{4}{\left(U_{j-2}-4U_{j-1}+3U_j\right)}^2,\hfill \\ \hfill {\beta }_1& =\frac{13}{12}{\left(U_{j-1}-2U_j+U_{j+1}\right)}^2+\frac{1}{4}{\left(U_{j-1}-U_{j+1}\right)}^2,\hfill \\ \hfill {\beta }_2& =\frac{13}{12}{\left(U_j-2U_{j+1}+U_{j+2}\right)}^2+\frac{1}{4}{\left(3U_j-4U_{j+1}+U_{j+2}\right)}^2.\hfill \end{array}$

For the right-biased value $U_{j+1/2}^{+}$, the same procedure is applied to the reversed stencil $\{U_{j+3},U_{j+2},U_{j+1},U_j,U_{j-1}\}$. The WENO scheme achieves fifth-order accuracy in smooth regions while maintaining stability near discontinuities. The spatial order of accuracy has been confirmed numerically using smooth initial conditions; see Section 4.1.

3.1.2 Source term: nonlocal alignment and obstacle force

The discrete source term evaluated at the cell interval $I_j$ is given by

$S_j={\displaystyle \sum _k}\varphi \left(|x_j-x_k|\right){\rho }_j{\rho }_k^{\tau }\left(u_k^{\tau }-u_j\right)\mathrm{\Delta }x+{\rho }_j{R}_j,(9)$

where ${\rho }_k^{\tau }≔{\rho }_k(t-\tau )$, $u_k^{\tau }≔u_k(t-\tau )$, and $\mathrm{\Delta }x$ denotes the uniform cell width of the finite volume mesh. The function $\varphi (\cdot )$ is the communication kernel defined in Equation 2. To facilitate efficient computation of the nonlocal alignment term in Equation 9, we introduce the auxiliary quantities:

$A_j≔{\displaystyle \sum _k}\varphi \left(|x_j-x_k|\right){\rho }_k^{\tau }{u}_k^{\tau }\mathrm{\Delta }x,B_j≔{\displaystyle \sum _k}\varphi \left(|x_j-x_k|\right){\rho }_k^{\tau }\mathrm{\Delta }x.$

Using these, the discrete source term can be reformulated as

$S_j={\rho }_j\left(A_j-u_j{B}_j\right)+{\rho }_j{R}_j.(10)$

We first focus on the efficient computation of the alignment contribution ${\rho }_j(A_j-u_j{B}_j)$ in Equation 10. Under periodic boundary conditions, these convolutions can be efficiently evaluated using the discrete Fast Fourier Transform (FFT). Specifically, the convolution of two periodic sequences $f_j$ and $g_j$ is given by

${\left(f\ast g\right)}_j={\displaystyle \sum _k}{f}_k{g}_{j-k},\text{which}\text{can}\text{be}\text{computed}\text{via}f\ast g={\mathcal{F}}^{-1}\left(\mathcal{F}\left(f\right)\cdot \mathcal{F}\left(g\right)\right),$

where $\mathcal{F}$ and ${\mathcal{F}}^{-1}$ denote the discrete Fourier transform and its inverse, respectively, and the multiplication $\cdot$ is performed element-wise in Fourier space. This approach reduces the computational complexity of evaluating each convolution from $\mathcal{O}(N^2)$ to $\mathcal{O}(N\log N)$, which is especially beneficial for large-scale simulations.

The second component of the source term, $R_j$, accounts for obstacle-induced repulsion potential and is modeled as the discrete gradient of a smooth obstacle potential field:

$R_j≔-{\partial }_x\mathrm{\Phi }\left(x_j\right),$

where $\mathrm{\Phi }(x)$ is given by Equation 3, and the gradient ${\partial }_x\mathrm{\Phi }$ is discretized using central finite differences.

3.1.3 Time discretization

For the time discretization, we apply the third-order Runge–Kutta method of Shu and Osher (SSP-RK3) [38–40], which offers third-order temporal accuracy and is known for its numerical stability to the semi-discrete system To advance the semi-discrete system in time, we adopt the third-order strong stability-preserving Runge–Kutta (SSP-RK3) method [38], applied to the system.

$\frac{d{U}_j}{dt}=\mathcal{L}\left(U_j\right),$

where $\mathcal{L}(U_j)$ is the spatial operator defined in Equation 6.

The SSP-RK3 scheme advances $U_j$ from $t^n$ to $t^{n+1}=t^n+\mathrm{\Delta }t$ via the following three-stage procedure:

$\begin{array}{ll}\hfill U_j^{\left(1\right)}& =U_j^n+\mathrm{\Delta }t\mathcal{L}\left(U_j^n\right),\hfill \\ \hfill U_j^{\left(2\right)}& =\frac{3}{4}{U}_j^n+\frac{1}{4}\left(U_j^{\left(1\right)}+\mathrm{\Delta }t\mathcal{L}\left(U_j^{\left(1\right)}\right)\right),\hfill \\ \hfill U_j^{n+1}& =\frac{1}{3}{U}_j^n+\frac{2}{3}\left(U_j^{\left(2\right)}+\mathrm{\Delta }t\mathcal{L}\left(U_j^{\left(2\right)}\right)\right).\hfill \end{array}$

For stability, the time step is chosen to satisfy the CFL condition

$\mathrm{\Delta }t\le \frac{\mathrm{C}\mathrm{F}\mathrm{L}\cdot \mathrm{\Delta }x}{{\max }_j{\lambda }_{j+\frac{1}{2}}},(11)$

where $\mathrm{C}\mathrm{F}\mathrm{L}\in (\mathrm{0,1})$ is a prescribed constant (e.g., $\mathrm{C}\mathrm{F}\mathrm{L}=0.7$) and local wave speed ${\lambda }_{j+\frac{1}{2}}$ as defined in Equation 8. This high-order scheme ensures third-order accuracy in time and is validated in Section 4.1 via convergence tests with smooth initial data.

3.2 Finite volume scheme for two-dimensional case

We consider a two-dimensional periodic domain $\mathrm{\Omega }=[a_x,b_x]\times [a_y,b_y]$, which is discretized into $N_x\times N_y$ uniform cells with mesh sizes $\mathrm{\Delta }x=(b_x-a_x)/N_x$, $\mathrm{\Delta }y=(b_y-a_y)/N_y$. Denote the center of cell $(i,j)$ by ${\mathbf{x}}_{i,j}=(x_i,y_j)$, where $x_i=a_x+\left(i-1/2\right)\mathrm{\Delta }x,y_j=a_y+\left(j-1/2\right)\mathrm{\Delta }y,$ for $i=1,\dots ,N_x$, $j=1,\dots ,N_y$. The control volume $C_{i,j}$ is defined as the rectangular cell centered at ${\mathbf{x}}_{i,j}$:

$C_{i,j}=\left[x_i-\frac{\mathrm{\Delta }x}{2},x_i+\frac{\mathrm{\Delta }x}{2}\right]\times \left[y_j-\frac{\mathrm{\Delta }y}{2},y_j+\frac{\mathrm{\Delta }y}{2}\right].$

Let the cell averages of the conserved variables be defined as

${\mathbf{U}}_{i,j}\left(t\right)=\left(\begin{array}{c}\hfill {\rho }_{i,j}\left(t\right)\hfill \\ \hfill m_{i,j}^x\left(t\right)\hfill \\ \hfill m_{i,j}^y\left(t\right)\hfill \end{array}\right)\approx \frac{1}{\mathrm{\Delta }x\mathrm{\Delta }y}{\iint }_{C_{i,j}}\left(\begin{array}{c}\hfill \rho \left(x,y,t\right)\hfill \\ \hfill \rho \left(x,y,t\right)u\left(x,y,t\right)\hfill \\ \hfill \rho \left(x,y,t\right)v\left(x,y,t\right)\hfill \end{array}\right)dxdy.$

The semi-discrete finite volume scheme is given by

$\frac{d}{dt}{\mathbf{U}}_{i,j}=-\frac{1}{\mathrm{\Delta }x}\left({\mathbf{F}}_{i+\frac{1}{2},j}^x-{\mathbf{F}}_{i-\frac{1}{2},j}^x\right)-\frac{1}{\mathrm{\Delta }y}\left({\mathbf{F}}_{i,j+\frac{1}{2}}^y-{\mathbf{F}}_{i,j-\frac{1}{2}}^y\right)+\left(\begin{array}{c}\hfill 0\hfill \\ \hfill {\mathbf{S}}_{i,j}\hfill \end{array}\right),$

where ${\mathbf{F}}^x,{\mathbf{F}}^y$ are numerical fluxes in the $x$- and $y$-directions, and ${\mathbf{S}}_{i,j}$ is the discrete source term. For convenience, we define the two-dimensional spatial operator $\mathcal{D}({\mathbf{U}}_{ij})$ as

$\mathcal{D}\left({\mathbf{U}}_{ij}\right)≔-\frac{1}{\mathrm{\Delta }x}\left({\mathbf{F}}_{i+\frac{1}{2},j}^x-{\mathbf{F}}_{i-\frac{1}{2},j}^x\right)-\frac{1}{\mathrm{\Delta }y}\left({\mathbf{F}}_{i,j+\frac{1}{2}}^y-{\mathbf{F}}_{i,j-\frac{1}{2}}^y\right)+\left(\begin{array}{c}\hfill 0\hfill \\ \hfill {\mathbf{S}}_{i,j}\hfill \end{array}\right).(12)$

3.2.1 2D Lax–Friedrichs flux with fifth-order WENO reconstruction

In our scheme, we consider pressureless Euler-type fluxes with the conserved variable and velocity field:

$\mathbf{U}=\left(\begin{array}{c}\hfill \rho \hfill \\ \hfill m^x\hfill \\ \hfill m^y\hfill \end{array}\right),\mathbf{u}=\frac{1}{\rho }\left(\begin{array}{c}\hfill m^x\hfill \\ \hfill m^y\hfill \end{array}\right)=\left(\begin{array}{c}\hfill u\hfill \\ \hfill v\hfill \end{array}\right).$

The physical fluxes in the $x$- and $y$-directions are given by

${\mathbf{F}}^x\left(\mathbf{U}\right)=\left(\begin{array}{c}\hfill m^x\hfill \\ \hfill \frac{{\left(m^x\right)}^2}{\rho }\hfill \\ \hfill \frac{m^x{m}^y}{\rho }\hfill \end{array}\right),{\mathbf{F}}^y\left(\mathbf{U}\right)=\left(\begin{array}{c}\hfill m^y\hfill \\ \hfill \frac{m^x{m}^y}{\rho }\hfill \\ \hfill \frac{{\left(m^y\right)}^2}{\rho }\hfill \end{array}\right).$

The Lax–Friedrichs numerical fluxes at cell interfaces are:

$\begin{array}{ll}\hfill {\mathbf{F}}_{i+\frac{1}{2},j}^x& =\frac{1}{2}\left[{\mathbf{F}}^x\left({\mathbf{U}}_{i+\frac{1}{2},j}^{-}\right)+{\mathbf{F}}^x\left({\mathbf{U}}_{i+\frac{1}{2},j}^{+}\right)\right]-\frac{{\lambda }_{i+\frac{1}{2},j}^x}{2}\left({\mathbf{U}}_{i+\frac{1}{2},j}^{+}-{\mathbf{U}}_{i+\frac{1}{2},j}^{-}\right),\hfill \\ \hfill {\mathbf{F}}_{i,j+\frac{1}{2}}^y& =\frac{1}{2}\left[{\mathbf{F}}^y\left({\mathbf{U}}_{i,j+\frac{1}{2}}^{-}\right)+{\mathbf{F}}^y\left({\mathbf{U}}_{i,j+\frac{1}{2}}^{+}\right)\right]-\frac{{\lambda }_{i,j+\frac{1}{2}}^y}{2}\left({\mathbf{U}}_{i,j+\frac{1}{2}}^{+}-{\mathbf{U}}_{i,j+\frac{1}{2}}^{-}\right).\hfill \end{array}$

The reconstructed states ${\mathbf{U}}_{i\pm \frac{1}{2},j}^{\pm }$, ${\mathbf{U}}_{i,j\pm \frac{1}{2}}^{\pm }$ at each interface are obtained using fifth-order WENO reconstruction in the respective spatial direction. The local characteristic wave speeds are given by the largest eigenvalues of the Jacobian matrices $\frac{\partial {\mathbf{F}}^x}{\partial \mathbf{U}},\frac{\partial {\mathbf{F}}^y}{\partial \mathbf{U}}$:

${\lambda }_{i+\frac{1}{2},j}^x=\max \left\{|u_{i+\frac{1}{2},j}^{-}|,|u_{i+\frac{1}{2},j}^{+}|\right\},{\lambda }_{i,j+\frac{1}{2}}^y=\max \left\{|v_{i,j+\frac{1}{2}}^{-}|,|v_{i,j+\frac{1}{2}}^{+}|\right\}.(13)$

Remark 3.2. Alternatively, a global maximum wave speed

$\lambda =\underset{i,j}{\max }\left\{|u_{i,j}|,|v_{i,j}|\right\}$

which is more diffusive but more robust near strong gradients or discontinuities.

We now describe the treatment of the source term in the semi-discrete scheme. As in the one-dimensional case, the discrete source term at the control volume $C_{i,j}$ consists of a nonlocal alignment interaction and an obstacle-induced repulsion force:

${\mathbf{S}}_{i,j}={\displaystyle \sum _{k,l}}\varphi \left(|{\mathbf{x}}_{i,j}-{\mathbf{x}}_{k,l}|\right){\rho }_{i,j}{\rho }_{k,l}^{\tau }\left({\mathbf{u}}_{k,l}^{\tau }-{\mathbf{u}}_{i,j}\right)\mathrm{\Delta }x\mathrm{\Delta }y+{\rho }_{i,j}{\mathbf{R}}_{i,j},(14)$

where the delayed fields are defined by ${\rho }_{k,l}^{\tau }≔{\rho }_{k,l}(t-\tau ),{\mathbf{u}}_{k,l}^{\tau }≔{\mathbf{u}}_{k,l}(t-\tau ).$ To facilitate efficient computation of the nonlocal alignment term in Equation 14, we introduce the auxiliary quantities:

${\mathbf{A}}_{i,j}={\displaystyle \sum _{k,l}}\varphi \left(|{\mathbf{x}}_{i,j}-{\mathbf{x}}_{k,l}|\right){\rho }_{k,l}^{\tau }{\mathbf{u}}_{k,l}^{\tau }\mathrm{\Delta }x\mathrm{\Delta }y,B_{i,j}={\displaystyle \sum _{k,l}}\varphi \left(|{\mathbf{x}}_{i,j}-{\mathbf{x}}_{k,l}|\right){\rho }_{k,l}^{\tau }\mathrm{\Delta }x\mathrm{\Delta }y.$

Using these, the source term Equation 14 can be written more compactly as:

${\mathbf{S}}_{i,j}={\rho }_{i,j}\left({\mathbf{A}}_{i,j}-{\mathbf{u}}_{i,j}{B}_{i,j}\right)+{\rho }_{i,j}{\mathbf{R}}_{i,j}.(15)$

We first focus on the alignment contribution ${\rho }_{i,j}\left({\mathbf{A}}_{i,j}-{\mathbf{u}}_{i,j}{B}_{i,j}\right)$ in Equation 15. Under periodic boundary conditions, it can be efficiently evaluated using the Fast Fourier Transform (FFT). Next, we consider the obstacle-induced repulsion force ${\mathbf{R}}_{i,j}$ is modeled as the negative gradient of a smooth potential field $\mathrm{\Phi }(\mathbf{x})$, evaluated at the grid point ${\mathbf{x}}_{i,j}$. It is discretized using central finite differences:

${\mathbf{R}}_{i,j}=-\nabla \mathrm{\Phi }\left({\mathbf{x}}_{i,j}\right).$

Consistent with the one-dimensional case (see Section 3.1.3), we adopt the third-order strong stability-preserving Runge–Kutta method (SSP-RK3) [38] for temporal discretization at each cell $(i,j)$:

${\mathbf{U}}_{i,j}^{\left(1\right)}={\mathbf{U}}_{i,j}^n+\mathrm{\Delta }t\mathcal{D}\left({\mathbf{U}}_{i,j}^n\right),$

${\mathbf{U}}_{i,j}^{\left(2\right)}=\frac{3}{4}{\mathbf{U}}_{i,j}^n+\frac{1}{4}\left({\mathbf{U}}_{i,j}^{\left(1\right)}+\mathrm{\Delta }t\mathcal{D}\left({\mathbf{U}}_{i,j}^{\left(1\right)}\right)\right),$

${\mathbf{U}}_{i,j}^{n+1}=\frac{1}{3}{\mathbf{U}}_{i,j}^n+\frac{2}{3}\left({\mathbf{U}}_{i,j}^{\left(2\right)}+\mathrm{\Delta }t\mathcal{D}\left({\mathbf{U}}_{i,j}^{\left(2\right)}\right)\right),$

where $\mathcal{D}({\mathbf{U}}_{i,j})$ is the spatial operator defined in Equation 12. The time step $\mathrm{\Delta }t$ is constrained by the CFL condition:

$\mathrm{\Delta }t\le \frac{\mathrm{C}\mathrm{F}\mathrm{L}\cdot \min \left(\mathrm{\Delta }x,\mathrm{\Delta }y\right)}{\underset{i,j}{\max }\left\{{\lambda }_{i+\frac{1}{2},j}^x,{\lambda }_{i,j+\frac{1}{2}}^y\right\}},$

where the wave speed ${\lambda }_{i+\frac{1}{2},j}^x,{\lambda }_{i,j+\frac{1}{2}}^y$ defined in Equation 13 and $\mathrm{C}\mathrm{F}\mathrm{L}\in (\mathrm{0,1})$.

3.3 Numerical conservation properties

In this section, we analyze the conservation properties of the proposed finite volume scheme. For clarity, the analysis is restricted to the one-dimensional case with periodic boundary conditions.

3.3.1 Mass conservation

Theorem 3.3. (Discrete mass conservation). The proposed finite volume scheme with fifth-order WENO reconstruction, Lax–Friedrichs flux, and SSP-RK3 time discretization exactly preserves the total mass under periodic boundary conditions:

${\displaystyle \sum _i}{\rho }_i^{n+1}\mathrm{\Delta }x={\displaystyle \sum _i}{\rho }_i^n\mathrm{\Delta }x.$

Proof. For convenience, we denote the first (density) component of the update operator $\mathcal{L}(U_j)$ defined in Equation 6 by ${\mathcal{L}}^1(U_j)$, i.e.,

${\mathcal{L}}^1\left(U_j\right)≔-\frac{1}{\mathrm{\Delta }x}\left({\widehat{\mathbf{F}}}_{j+\frac{1}{2}}^1-{\widehat{\mathbf{F}}}_{j-\frac{1}{2}}^1\right),$

where ${\widehat{\mathbf{F}}}_{j\pm \frac{1}{2}}^1$ denotes the first component of the numerical flux ${\mathbf{F}}_{j\pm \frac{1}{2}}$ defined in Equation 7, corresponding to the mass flux.

Let ${\rho }_i^n$ denote the cell average of density at time $t^n$. The SSP-RK3 updates for density read:

${\rho }_i^{\left(1\right)}={\rho }_i^n+\mathrm{\Delta }t{\mathcal{L}}^{\left(1\right)}\left(U_i^n\right),(16)$

${\rho }_i^{\left(2\right)}=\frac{3}{4}{\rho }_i^n+\frac{1}{4}\left({\rho }_i^{\left(1\right)}+\mathrm{\Delta }t{\mathcal{L}}^{\left(1\right)}\left(U_i^{\left(1\right)}\right)\right),(17)$

${\rho }_i^{n+1}=\frac{1}{3}{\rho }_i^n+\frac{2}{3}\left({\rho }_i^{\left(2\right)}+\mathrm{\Delta }t{\mathcal{L}}^{\left(1\right)}\left(U_i^{\left(2\right)}\right)\right).(18)$

Define the total mass at each stage as

$M^{\left(*\right)}≔{\displaystyle \sum _i}{\rho }_i^{\left(*\right)}\mathrm{\Delta }x,\text{for\hspace{0.17em}}\ast \in \left\{n,\left(1\right),\left(2\right),n+1\right\}.$

By using the periodic boundary condition, we have

${\displaystyle \sum _i}\left({\widehat{\mathbf{F}}}_{i+1/2}^{\left(1\right)}-{\widehat{\mathbf{F}}}_{i-1/2}^{\left(1\right)}\right)=0,$

summing Equations 16–18 over all $i$, we find:

$M^{\left(1\right)}=M^n,M^{\left(2\right)}=M^n,M^{n+1}=M^n.$

Thus, the scheme is exactly mass conservative.

3.3.2 Momentum conservation

Theorem 3.4. (Momentum Conservation). Under periodic boundary conditions and in the absence of delays and external forces, the proposed finite volume scheme with fifth-order WENO reconstruction, Lax–Friedrichs flux, and SSP-RK3 time discretization exactly preserves the total discrete momentum at each time step:

${\displaystyle \sum _i}{\left(\rho u\right)}_i^{n+1}\mathrm{\Delta }x={\displaystyle \sum _i}{\left(\rho u\right)}_i^n\mathrm{\Delta }x.$

Proof. We denote the second (momentum) component of the update operator $\mathcal{L}(U_i)$ defined in Equation 6 by ${\mathcal{L}}^2(U_i)$, i.e.,

${\mathcal{L}}^2\left(U_i\right)≔-\frac{1}{\mathrm{\Delta }x}\left({\widehat{\mathbf{F}}}_{i+\frac{1}{2}}^2-{\widehat{\mathbf{F}}}_{i-\frac{1}{2}}^2\right)+S_i,$

where ${\widehat{\mathbf{F}}}_{i\pm \frac{1}{2}}^2$ denotes the second component of the numerical flux ${\mathbf{F}}_{i\pm \frac{1}{2}}$ defined in Equation 7, and $S_i$ is the discrete source term defined in Equation 9.

Let ${(\rho u)}_i^n$ be the cell average of momentum at time $t^n$. The SSP-RK3 scheme updates for momentum are given by:

${\left(\rho u\right)}_i^{\left(1\right)}={\left(\rho u\right)}_i^n+\mathrm{\Delta }t{\mathcal{L}}_i^2\left(U^n\right),(19)$

${\left(\rho u\right)}_i^{\left(2\right)}=\frac{3}{4}{\left(\rho u\right)}_i^n+\frac{1}{4}\left({\left(\rho u\right)}_i^{\left(1\right)}+\mathrm{\Delta }t{\mathcal{L}}_i^2\left(U^{\left(1\right)}\right)\right),(20)$

${\left(\rho u\right)}_i^{n+1}=\frac{1}{3}{\left(\rho u\right)}_i^n+\frac{2}{3}\left({\left(\rho u\right)}_i^{\left(2\right)}+\mathrm{\Delta }t{\mathcal{L}}_i^2\left(U^{\left(2\right)}\right)\right).(21)$

Define the total momentum at each stage as

$P^{\left(*\right)}≔{\displaystyle \sum _i}{\left(\rho u\right)}_i^{\left(*\right)}\mathrm{\Delta }x,\text{for\hspace{0.17em}}\ast \in \left\{n,\left(1\right),\left(2\right),n+1\right\}.$

Under periodic boundary conditions, it is easy to see that:

${\displaystyle \sum _i}\left({\widehat{\mathbf{F}}}_{i+\frac{1}{2}}^{\left(2\right)}-{\widehat{\mathbf{F}}}_{i-\frac{1}{2}}^{\left(2\right)}\right)=0.(22)$

Therefore, the total contribution from the update operator reduces to the source term:

${\displaystyle \sum _i}{\mathcal{L}}^{\left(2\right)}\left(U_i\right)\mathrm{\Delta }x={\displaystyle \sum _i}{S}_i\mathrm{\Delta }x.$

In the special case when $\tau =0$ and $R_i=0$, the source term $S_i$ defined in Equation 9 simplifies to

$S_i={\displaystyle \sum _j}{\varphi }_{ij}{\rho }_i{\rho }_j\left(u_j-u_i\right)\mathrm{\Delta }x,(23)$

where ${\varphi }_{ij}≔\varphi (|x_i-x_j|)$. Since ${\varphi }_{ij}={\varphi }_{ji}$ is symmetric in $i,j$, and $(u_j-u_i)=-(u_i-u_j)$ is antisymmetric, the double sum satisfies

${\displaystyle \sum _{i,j}}{\varphi }_{ij}{\rho }_i{\rho }_j\left(u_j-u_i\right)=-{\displaystyle \sum _{i,j}}{\varphi }_{ji}{\rho }_j{\rho }_i\left(u_i-u_j\right)=-{\displaystyle \sum _{i,j}}{\varphi }_{ij}{\rho }_i{\rho }_j\left(u_j-u_i\right).$

This implies

${\displaystyle \sum _{i,j}}{\varphi }_{ij}{\rho }_i{\rho }_j\left(u_j-u_i\right)=0.$

Therefore, from Equation 23 we have

${\displaystyle \sum _i}{S}_i\mathrm{\Delta }x={\displaystyle \sum _{i,j}}{\varphi }_{ij}{\rho }_i{\rho }_j\left(u_j-u_i\right){\left(\mathrm{\Delta }x\right)}^2=0.(24)$

Summing Equations 19–21 over all $i$, in view of Equation 22 and Equation 24, we see that the total momentum is conserved during the time update, i.e.,

$P^{\left(1\right)}=P^n,P^{\left(2\right)}=P^n,P^{n+1}=P^n.$

Thus we obtain the conclusion.

Remark 3.5. When delay $(\tau >0)$ is present$\tau >0$, the source term involves past states ${\rho }_j^{\tau }{u}_j^{\tau }$, breaking the symmetry required for momentum conservation. Similarly, nonzero obstacle forces $R_i\ne 0$ introduce additional asymmetries. As a result, the total momentum is no longer conserved. This loss of conservation is clearly observed in the numerical results of Example 4.3.

Remark 3.6. (Energy fluctuation). We introduce the energy fluctuation [32].

$\mathrm{\Delta }E\left(t\right)≔{\iint }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\frac{|\mathbf{u}(\mathbf{x})-\mathbf{u}(\mathbf{y}){|}^2}{2}\rho \left(\mathbf{x}\right)\rho \left(\mathbf{y}\right)d\mathbf{x}d\mathbf{y}=\mathbf{M}{\int }_{{\mathbb{R}}^d}|\mathbf{u}(\mathbf{x})-{\mathbf{u}}_{\mathit{\infty }}{|}^2\rho \left(t,\mathbf{x}\right)d\mathbf{x},$

where $\mathbf{M}={\int }_{{\mathbb{R}}^d}\rho (\mathbf{x})d\mathbf{x}$ is the total mass and ${\mathbf{u}}_{\mathit{\infty }}≔\frac{1}{\mathbf{M}}{\int }_{{\mathbb{R}}^d}\mathbf{u}(t,\mathbf{x})\rho (t,\mathbf{x})d\mathbf{x}$ denotes the mass–averaged velocity. In numerical experiments, we plot the energy fluctuation in Example 4.3 (see Figure 1) and observe that it tends to zero. This indicates that the numerical solution is stable. It should be emphasized that in this work we testify the stability of the numerical solution only from a experimental perspective, while a rigorous theoretical proof of stability will be left for future study.

Figure 1

Figure 1. Time evolution of total mass, momentum, and energy fluctuation for the three configurations considered in Example 4.3.

4 Numerical experiments

This section presents numerical simulations. Throughout all simulations, the interaction kernel is chosen as $\varphi (r)=\delta /{(1+r^2)}^{\theta },\theta =3$, where the constant $\delta$ is chosen such that ${\int }_0^{\infty }\varphi (r)dr=1$, and the CFL number is taken as 0.7, unless explicitly stated otherwise. Periodic boundary conditions are imposed in all cases. In Section 4.1, we first validate the spatial fifth-order and temporal third-order accuracy of the proposed scheme through convergence tests under two different configurations. In addition, we verify the conservation of mass and momentum through a dedicated numerical experiment, and assess the numerical stability of the scheme by monitoring the energy fluctuation. Then, in Section 4.2, we apply the numerical scheme to two models, leading to several interesting numerical observations. In total, seven numerical experiments are presented to investigate the effects of the delay parameter $\tau$ and the presence of obstacles on the emergence of collective behaviors such as flocking and non-flocking in both single-species and two-species systems.

4.1 Convergence and conservation tests

Example 4.1. (Spatial convergence).

To test the spatial accuracy of the fifth-order WENO reconstruction, we consider a smooth initial condition for the single-species system Equation 1, defined on the periodic domain $[-L/2,L/2]$ with $L=2$:

${\rho }_0\left(x\right)={\gamma }_1\left(2+\sin \left(\frac{\pi x}{L}\right)\right),u_0\left(x\right)\equiv 1,$

where the normalization constant ${\gamma }_1$ is chosen such that ${\int }_{-L/2}^{L/2}{\rho }_0(x)dx=1$. We vary the number of cells as $N_x=\mathrm{10,20,40,80}$ and set the time step according to $\mathrm{\Delta }t\sim \mathrm{\Delta }{x}^{5/3}$, where $\mathrm{\Delta }x=L/N_x$. We compute the cell averages and compare the numerical solution at time $T=2$ with the exact solution ${\rho }_{\text{ref}}$ and ${(\rho u)}_{\text{ref}}$.

As shown in Table 1, the $L^1$ errors for $\rho$ and $\rho u$ converge with an order approaching 5 as the grid is refined, where $\Vert \rho -{\rho }_{\text{ref}}{\Vert }_{L^1}={\sum }_i\mathrm{\Delta }x|{\rho }_i-{\rho }_{\text{ref},i}|$ and $\Vert \rho u-{(\rho u)}_{\text{ref}}{\Vert }_{L^1}={\sum }_i\mathrm{\Delta }x|{(\rho u)}_i-{(\rho u)}_{\text{ref},i}|$. This confirms that the scheme achieves fifth-order spatial accuracy, consistent with the WENO reconstruction employed in the spatial discretization.

Table 1

Table 1. Accuracy test for Example 4.1 at $t=2$.

Example 4.2. (Temporal convergence). To assess the temporal accuracy of the third-order Runge–Kutta method, we consider a smooth initial condition for the single-species system Equation 1, defined on the periodic domain $[-L/2,L/2]$ with $L=2$:

${\rho }_0\left(x\right)={\gamma }_2\left[\exp \left(-\frac{{\left(x+L/4\right)}^2}{{\sigma }^2}\right)+\exp \left(-\frac{{\left(x-L/4\right)}^2}{{\sigma }^2}\right)\right],u_0\left(x\right)\equiv 1,$

where $\sigma =0.2$, and the normalization constant ${\gamma }_2$ is chosen such that ${\int }_{-L/2}^{L/2}{\rho }_0(x)dx=1$. We fix the number of spatial cells at $N_x=2000$ and vary the time step $\mathrm{\Delta }t$ by adjusting the CFL number defined in Equation 11 accordingly.

As shown in Table 2, the $L^1$ errors converge with an order approaching 3 as the CFL number decreases, confirming the third-order temporal accuracy of the Runge–Kutta scheme. The observed convergence order aligns well with the theoretical order of the time discretization.

Table 2

Table 2. Accuracy test for Example 4.2 at $t=0.2$.

Example 4.3. (Verification of conservation properties). To verify the theoretical results in Theorems 3.3 and 3.4, we present numerical simulations for the single-species system Equation 1 under three configurations: (a) without delay or obstacle, (b) with delay but no obstacle, and (c) with an obstacle force but no delay.

The initial data are set as

$\rho \left(x,0\right)={\gamma }_3\cos \left(\frac{\pi x}{L}\right),u\left(x,0\right)=-c\cos \left(\frac{\pi x}{L}\right),$

where the normalization constant ${\gamma }_3$ is chosen so that ${\int }_{-L/2}^{L/2}{\rho }_0(x)dx=1$, and $c=0.1$ controls the velocity amplitude. The simulation is conducted on a periodic domain $x\in [-\mathrm{1,1}]$ (i.e., $L=2$), with spatial resolution $N_x=1800$.

In the third scenario, we introduce an obstacle modeled as the gradient of a repulsive potential centered at ${\mathbf{x}}_0=0$, as defined in Equation 3:

$\mathbf{R}\left(\mathbf{x}\right)=-\nabla \mathrm{\Phi }\left(\mathbf{x}\right),\text{with}\mathrm{\Phi }\left(\mathbf{x}\right)={\eta }_0\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}_0{|}^2}{l_0^2}\right),$

with repulsion parameters ${\eta }_0=0.01$ and $l_0=0.5$.

Figure 1 shows the time evolution of total mass, momentum, and energy fluctuation $\mathrm{\Delta }E$ (defined in Remark 3.6) for each configuration, with the top, middle, and bottom panels displaying mass $M(t)≔{\int }_{-L/2}^{L/2}\rho (x,t)dx$, momentum $P(t)≔{\int }_{-L/2}^{L/2}\rho (x,t)u(x,t)dx$, and $\mathrm{\Delta }E$, respectively. Figure 1 displays the time evolution of total mass and momentum for each configuration. In each subfigure, the top panel shows the total mass $M(t)≔{\int }_{-L/2}^{L/2}\rho (x,t)dx$, and the bottom panel shows the total momentum $P(t)≔{\int }_{-L/2}^{L/2}\rho (x,t)u(x,t)dx$. The results confirm the theoretical predictions: mass is conserved with high accuracy in all cases, in agreement with Theorem 3.3. Momentum is preserved only in the absence of delay and obstacle, as stated in Theorem 3.4. In the presence of delay or obstacle force, the loss of symmetry in the source term leads to gradual momentum deviation, as anticipated in the accompanying Remark 3.5. Furthermore, Figure 1 illustrates the evolution of $\mathrm{\Delta }E$. In panels (a) and (b), $\mathrm{\Delta }E$ approaches to zero by $T=10$, while in panel (c) it tends to zero around $T=60$. These results suggest that the scheme exhibits stable behavior under the considered scenarios.

4.2 Effects of delay and obstacle on flocking dynamics

4.2.1 Effects of delay and obstacle for single-species

First, we examine the effects of the time delay $\tau$ and obstacle force on the collective behavior and regularity of solutions in the single-species system Equation 1. The initial data is given by

$\rho \left(x,0\right)={\gamma }_3\cos \left(\frac{\pi x}{L}\right),u\left(x,0\right)=-c\sin \left(\frac{2\pi x}{L}\right),$

where the normalization constant ${\gamma }_3$ is chosen such that ${\int }_{-L/2}^{L/2}\rho (x,0)dx=1$, and the parameter $c>0$ controls the amplitude of the initial velocity. The simulation is performed on the periodic domain $[-\mathrm{1,1}]$, i.e., $L=2$, with the number of cell $N_x=1800$.

Two Examples 4.4 and 4.5 are presented below. In the first, no obstacle is introduced, and we focus on the effects of delay with varying initial velocity amplitudes. In the second, an obstacle is included to study the combined effects of delay and obstacle forces.

Example 4.4. (Effects of delay without obstacle). We consider two initial velocity amplitudes: $c=0.10$ and $c=0.17$. For each case, we vary the delay parameter $\tau \in \{\mathrm{0,0.5},3\}$ to examine its influence on the long-time behavior.

Figure 2 presents the density and velocity $(\rho ,u)$ at final time $T=15$ for $c=0.10$. In all delay settings, the solution remains globally smooth, and the density stays bounded. As the delay increases, the alignment process becomes slower and the density profile exhibits mild increase, but no signs of instability or singularity are observed.

In contrast, Figure 3 presents the results for $c=0.17$. When $\tau =0$, the solution remains smooth throughout the simulation. However, as $\tau$ increases, the density begins to concentrate more sharply. In particular, for $\tau =3$, we observe that the solution develops a near-singular profile: the density becomes highly concentrated around $t\approx 4.79$, and the velocity field exhibits a steep gradient. These observations suggest the emergence of near-singular behavior, potentially indicating the onset of instability or breakdown driven by the delay effect.

These numerical results show that the influence of delay $\tau$ on the single-species system strongly depends on the choice of initial conditions. For mild initial configurations, delay mainly slows down the flocking convergence without causing instability or loss of regularity. However, when the initial velocity is sufficiently large, the same delay can intensify velocity gradients and density concentrations, eventually leading to a breakdown of smoothness or finite-time singularity beyond a critical delay threshold. This highlights a nonlinear interplay between initial conditions and time delay in determining the long-time stability and regularity of the system.

Figure 2

Figure 2. Time evolution of density $\rho$ and velocity $u$ for $\tau =0$, 0.5, 3 with fixed $c=0.1$ (Example 4.4).

Figure 3

Figure 3. Time evolution of density $\rho$ and velocity $u$ for $\tau =0$, 0.5, 3 with fixed $c=0.17$ (Example 4.4).

Example 4.5. (Effects of delay with obstacle). In this example, we introduce an obstacle modeled as the gradient of a repulsive potential centered at ${\mathbf{x}}_0=0$, as defined in Equation 3:

$\mathbf{R}\left(\mathbf{x}\right)=-\nabla \mathrm{\Phi }\left(\mathbf{x}\right),\text{with}\mathrm{\Phi }\left(\mathbf{x}\right)={\eta }_0\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}_0{|}^2}{l_0^2}\right),$

with repulsion parameters ${\eta }_0=0.02$ and $l_0=0.2$. The initial velocity amplitude is set to $c=0.17$, and we examine the impact of varying the delay parameter $\tau \in \{\mathrm{0,0.5},3\}$.

Figure 4 shows the density and velocity fields at the moments when the system reaches its peak aggregation for different delays. In all three cases, the density is observed to concentrate sharply at two distinct points symmetrically located on either side of the obstacle. Correspondingly, the velocity field exhibits steep gradients precisely at these points.

Notably, the time at which this concentrated state emerges becomes earlier as the delay increases: $t=3.95$ for $\tau =0$, $t=3.68$ for $\tau =0.5$, and $t=3.00$ for $\tau =3$. Furthermore, the aggregation points move progressively closer to the obstacle as the delay increases, suggesting that larger delays enhance the effective attraction toward the obstacle and accelerate the localization process. Under the current setting of the communication kernel and repulsive potential, these findings indicate a synergistic interaction between delay and obstacle: the presence of delay not only triggers earlier aggregation but also intensifies density concentration near the obstacle.

Figure 4

Figure 4. Time evolution of density $\rho$ and velocity $u$ for $\tau =0$, 0.5, 3 with fixed $c=0.17$ (Example 4.5).

4.2.2 Effects of delay and obstacle for two-species

To investigate the effects of time delay and obstacles in a two-species setting, we simulate the system Equation 4 with the following different initial conditions for the two species:

$\begin{array}{cccc}\hfill {\rho }_1\left(x,0\right)& ={\lambda }_1\left[0.1+\exp \left(-x^2\right)\right],\hfill & \hfill u_1\left(x,0\right)& =-c_1\sin \left(\frac{2\pi x}{L}\right),\hfill \\ \hfill {\rho }_2\left(x,0\right)& ={\lambda }_2\cos \left(\frac{\pi x}{L}\right),\hfill & \hfill u_2\left(x,0\right)& =-c_2\sin \left(\frac{2\pi x}{L}\right),\hfill \end{array}$

with $L=2$ and constants $c_1=0.3,c_2=0.32$. The normalization constants ${\lambda }_1,{\lambda }_2$ are chosen such that ${\int }_{-L/2}^{L/2}{\rho }_i(x)dx=1,i=\mathrm{1,2}$. The simulation is conducted using a periodic domain discretized with the number of cell $N_x=1800$ and a CFL number of 0.7, running from $t=0$ to $t=15$.

For this setting, we provide three examples. Examples 4.6 and 4.7 are presented first. In the first Example 4.6, no obstacle is introduced, and we focus on the effect of delay on the collective dynamics of the two-species system. In the second Example 4.7, an obstacle is included to study the combined effects of delay and obstacle forces.

To further extend the study, the last Example 4.8 is designed similarly to Example 4.6, but with asymmetric inter-species kernels.

Example 4.6. (Effects of delay without obstacle). In this example, we vary the delay parameter $\tau \in \{0,0.5,3\}$ to examine its influence on the long-time behavior of the two-species system.

According to the numerical analysis in Example 4.4, both of the initial conditions $({\rho }_1,u_1)$ and $({\rho }_2,u_2)$, if evolved independently in a single-species setting, would lead to finite-time singularity, regardless of the presence of delay. However, our numerical results show that when $\tau =0$, the two-species system remains smooth and well-behaved for all time (see Figures 5a-d). This indicates that cross-species interaction can regularize the system.

As the delay $\tau$ increases, the solution exhibits progressively sharper spatial structures: velocity gradients steepen, and the density profile becomes increasingly concentrated over time (see Figures 5a–c for density; Figures 5d–f for velocity). In particular, at time $t=2.70$, Figures 5c-f show pronounced density localization near a single point and steep velocity gradients. These phenomena indicate a growing tendency toward instability and suggest a potential breakdown of collective coherence.

Our numerical experiments indicate that there may exist a critical delay threshold ${\tau }_c$, beyond which the system exhibits a loss of regularity and a breakdown of coordinated motion.

Figure 5

Figure 5. Time evolution of density $\rho$ and velocity $u$ for $\tau =0$, 0.5, 3 (Example 4.6).

Remark 4.1. A possible mechanism for this loss of regularity is that the use of delayed velocity information in the interaction term introduces misalignment between agents, potentially weakening the stabilizing of the alignment mechanism. Moreover, the time-lagged nonlocal feedback causes the system’s response to deviate from the current state, amplifying local gradients and driving the system toward instability.

Example 4.7. (Effects of delay with obstacle). In this example, we introduce an obstacle modeled as the gradient of a repulsive potential centered at ${\mathbf{x}}_{\alpha }=0$, as defined in Equation 5:

${\mathbf{R}}_{\alpha }\left(\mathbf{x}\right)=-\nabla {\mathrm{\Phi }}_{\alpha }\left(\mathbf{x}\right),{\mathrm{\Phi }}_{\alpha }\left(\mathbf{x}\right)={\eta }_{\alpha }\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}_{\alpha }{|}^2}{l_{\alpha }^2}\right),$

with repulsion parameters ${\eta }_{\alpha }=0.1$ and $l_{\alpha }=0.2$, $\alpha =\mathrm{1,2}$. The obstacle force is taken to be the same for both species, i.e., ${\mathbf{R}}_1={\mathbf{R}}_2$. Figure 6 illustrates the evolution of the two-species density $\rho (x,t)$ and velocity field $u(x,t)$ in the presence of a fixed obstacle force and varying delay parameters $\tau =0$, 0.5, and 3.

Across all cases, we observe the emergence of symmetric, highly concentrated density peaks on both sides of the obstacle. These are accompanied by steep gradients in the velocity field, suggesting the formation of localized structures with near-singular behavior, as in the single-species case (Example 4.5).

As the delay parameter $\tau$ increases, the formation of these concentrated structures occurs earlier in time: at $t=1.60$ for $\tau =0$, $t=1.43$ for $\tau =0.5$, and $t=1.34$ for $\tau =3$. Moreover, the aggregation regions become more sharply localized and shift closer to the obstacle center as $\tau$ increases. These trends indicate that the time delay not only enhances the spatial focusing induced by the obstacle but also accelerates the onset of singular-like configurations.

The numerical results reveal that, similar to the single-species scenario and under the current settings of the communication kernel and repulsive potential, the interaction between delay and the obstacle plays a significant role in promoting rapid localization and potential loss of regularity. Although the two-species system introduces additional inter-species interactions and complexity, the observed trends persist: increasing delay leads to earlier aggregation and sharper density concentration near the obstacle.

Figure 6

Figure 6. Time evolution of density $\rho$ and velocity $u$ for $\tau =0$, 0.5, 3 (Example 4.7).

Example 4.8. (Effects of delay with asymmetric inter-species kernels). In this example, we adopt the same setting as in Example 4.6, but modify the communication kernels by choosing ${\theta }_{11}={\theta }_{12}={\theta }_{22}=3$ and ${\theta }_{21}=1.2$. Hence, species 1 applies the same exponent $\theta =3$ for both self- and cross-interactions, while species 2 interacts with species 1 through the smaller exponent ${\theta }_{21}=1.2$. This asymmetry highlights the effect of asymmetric inter-species interactions on the collective dynamics.

For the different initial conditions $({\rho }_1(x,0)\ne {\rho }_2(x,0))$ chosen for the two species, we observe that when $\tau =0$ and 0.5, the solution behavior is similar to that in Example 4.6. However, when $\tau =3$, Examples 4.6 and 4.8 show that the aggregation levels of ${\rho }_1$ and ${\rho }_2$ differ significantly (see Figures 5c, 7c). We also find that asymmetric and symmetric kernels affect the solution behavior differently, and a more detailed study will be left for future work.

Figure 7

Figure 7. Time evolution of density $\rho$ and velocity $u$ for $\tau =0$, 0.5, 3 (Example 4.8).

4.2.3 Effects of delay and obstacle in 2D

Now, we consider the two-dimensional single-species system Equation 1, incorporating both communication delay and a repulsive obstacle potential. The simulation is conducted on a periodic square domain ${\left[-\frac{L}{2},\frac{L}{2}\right]}^2$ with $L=1.2$, uniformly discretized using $N_x=N_y=120$ grid points. The CFL number is set to 0.012. The initial conditions are specified as follows:

${\rho }_0\left(x,y\right)=\gamma \cos \left(\frac{\pi x}{L}\right)\cos \left(\frac{\pi y}{L}\right),$

$u_0\left(x,y\right)=-c\sin \left(\frac{2\pi x}{L}\right),v_0\left(x,y\right)=-c\sin \left(\frac{2\pi y}{L}\right).$

where the normalization constant $\gamma$ is selected such that ${\iint }_{{[-\frac{L}{2},\frac{L}{2}]}^2}{\rho }_0(x,y)dxdy=1$ and the initial velocity parameter $c=0.12$. We investigate the long-time evolution of the system for various delay values $\tau =\mathrm{0,0.5},\mathrm{1,5}$.

Example 4.9. (Effects of delay without obstacle). In this example, we examine the system without repulsive potential. This analysis demonstrates the influence of increasing delay on collective dynamics and pattern formation in the 2D setting.

For $c=0.12$, the two-dimensional single-species system remains smooth and appears to converge toward a steady state in the absence of delay $(\tau =0)$, as shown in Figure 8. Both the density and velocity fields gradually stabilize, with the velocity field exhibiting alignment and the density approaching a spatially symmetric configuration.

As the delay parameter $\tau$ increases, the system displays more pronounced dynamical behavior. As illustrated in Figures 9a–c, the density gradually becomes more concentrated near the center of the domain. This tendency appears to be enhanced with larger values of $\tau$, suggesting that communication delay enhances aggregation in the density field. Meanwhile, the velocity fields shown in Figures 9d–f exhibit growing misalignment and irregularities over time, in contrast to the more coordinated behavior observed when $\tau =0$.

These numerical observations indicate the possible existence of a critical delay threshold ${\tau }_c$, beyond which the global regularity of the solution may deteriorate and the system may begin to exhibit signs of instability. This behavior is consistent with similar trends reported in one-dimensional single- and two-species models, where increased communication delay can exert a destabilizing influence on the dynamics. The observed mismatch between the rapid central concentration of density and the slower alignment of the velocity fields suggests that increasing delay may impair the system’s ability to maintain coherent collective dynamics and stable spatial patterns.

Figure 8

Figure 8. Time evolution of density $\rho$ and velocity $u$ field $\mathbf{u}$ for $\tau =0$ (Example 4.9). In 2D, the velocity field $\mathbf{u}$ is represented by its speed $|\mathbf{u}|$ (color) and direction (arrows).

Figure 9

Figure 9. Final-time density $\rho$ and velocity $u$ field $\mathbf{u}$ for $\tau =0.5$, 1, 5 (Example 4.9).

Example 4.10. (Effects of delay with obstacle). We introduce an obstacle force $\mathbf{R}(\mathbf{x})$ modeled as the gradient of a repulsive potential centered at the obstacle location, as defined in Equation 3:

$\mathbf{R}\left(\mathbf{x}\right)=-\nabla \mathrm{\Phi }\left(\mathbf{x}\right),\text{with}\mathrm{\Phi }\left(\mathbf{x}\right)={\eta }_0\exp \left(-\frac{|\mathbf{x}-{\mathbf{x}}_0{|}^2}{l_0^2}\right).$

In our numerical example, we consider a single obstacle located at ${\mathbf{x}}_0=(\mathrm{0,0})$, with repulsion parameters ${\eta }_0=0.01$ and $l_0=0.1$.

Figure 10 presents the full temporal evolution of the one-species density field under zero delay $(\tau =0)$, while Figure 11 displays the final-time density distributions for increasing delays $\tau =0.5$, 1, and 5. In all cases, the obstacle potential remains fixed, and only the delay parameter varies.

The introduction of the obstacle induces a characteristic spatial segregation: the density splits into two symmetric high-density regions located on either side of the obstacle. This behavior is consistent with earlier observations in the one-dimensional single-species case with obstacle (Example 4.5) and the one-dimensional two-species case with obstacle (Example 4.7), indicating a persistent effect of obstacle-induced localization across dimensions and system complexity.

As the delay $\tau$ increases, two systematic trends are observed. First, the system develops sharply concentrated density states at progressively earlier times. Specifically, the onset of significant aggregation shifts forward from $t=8.68$ for $\tau =0$ to $t=7.02$ for $\tau =0.5$, $t=5.89$ for $\tau =1$, and $t=4.17$ for $\tau =5$, as shown in Figure 11. This confirms that communication delay accelerates the aggregation process. Second, the spatial location of these dense regions moves consistently closer to the obstacle as $\tau$ increases. As depicted in Figures 10c, 11a–c, the high-density zones not only become more sharply localized but also shift toward the center of the obstacle.

Taken together, under the current setting of the communication kernel and repulsive potentials, the findings from all three settings—1D single species with obstacle (Example 4.5), 1D two species with obstacle (Example 4.7), and 2D single species with obstacle—reveal a consistent mechanism: increasing communication delay (i) advances the onset of strong aggregation and (ii) enhances spatial localization near repulsive obstacles. This delay-obstacle interplay appears to persist across spatial dimensions and system configurations.

Theoretically, in the hydrodynamic setting, Choi and Haskovec [22] established sufficient conditions for global regularity and flocking under normalized communication weights: (a) For a fixed integral influence function, it is necessary to choose $\tau$ sufficiently small in order to satisfy the flocking condition; (b) For the fat-tailed kernel $1/{(1+r^2)}^{\theta }$ with $\theta \in (\mathrm{0,1}/2)$, any $\tau$ can satisfy the flocking condition. We extend the theoretical case to the scenario with obstacles and two species. In our work, we intend to investigate the case (a) of an integral kernel of the form $\varphi (r)=\delta /{(1+r^2)}^{\theta }$ with $\theta >1/2$. Our numerical experiments show that small delays do not affect system alignment, but large delays can lead to steep gradients and high-density peaks, which appear numerically as near-singular behavior. Moreover, in the presence of obstacles, delays accelerate the aggregation of the system, and as the delay increases, the aggregation center moves closer to the obstacle. Furthermore, we also examine the two-species case with delays and obstacles. These observations are based on numerical evidence under specified parameters (e.g., potential shape, strength, and initial conditions). Nonetheless, it should be noted that these results are obtained under specified parameters (e.g., potential shape, strength, and initial conditions) in each example. While the trends are consistent across scenarios, further investigation is needed to assess their generality under broader modeling assumptions and in more realistic biological or physical systems. Rigorous analysis, such as establishing precise thresholds for finite-time singularity or delay-induced blow-up, is left for future work.

Figure 10

Figure 10. Time evolution of density $\rho$ and velocity $u$ field $\mathbf{u}$ for $\tau =0$ (Example 4.10).

Figure 11

Figure 11. Final-time density $\rho$ and velocity $u$ field $\mathbf{u}$ for $\tau =0.5$, 1, 5 (Example 4.10).

Remark 4.2. This work focuses on fat-tailed kernels. We also test compactly supported kernels (such as $\varphi (r)=\delta {\mathbf{1}}_{\{|r|<d\}}$), and we find that as the delay increases, the density tends to concentrate locally and the velocity gradient becomes large, similar to Example 4.4 in the single-species case. However, the behavior of the solutions remains different under the fat-tailed and compactly supported kernels, we will investigate compactly supported kernels in more detail in future work.

5 Conclusion

This study examines the influence of time delay and obstacle on the collective dynamics of non-local kinetic models in one- and two-dimensional settings for both single- and two-species systems. The numerical results of six representative cases reveal that the long-term behavior of the system is highly sensitive to initial conditions, with different initial states that lead to different outcomes, such as global regularity, aggregation or finite-time singularity formation. Increasing the time delay generally reduces stability, promoting earlier formation of singularities. Moreover, the presence of static obstacles combined with delay accelerates singularity onset and spatially shifts aggregation closer to the obstacle. These findings highlight the intricate interplay between initial data, delay, species interactions, and environmental heterogeneity in shaping emergent patterns. While our numerical experiments focus on fat-tailed kernels and isotropic Gaussian obstacles, the framework is readily applicable to other types of interaction kernels and potential functions. Future work will extend the model to include singular kernels, dynamic or reactive obstacles such as moving predators, and asymmetric inter-species interactions and anisotropic obstacles to explore richer collective dynamics, so as to better capture realistic ecological scenarios. Analytical characterization of critical delay thresholds and singularity formation also remains an important avenue for further research.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

TZ: Investigation, Writing – review and editing, Formal Analysis, Data curation, Methodology, Writing – original draft. GZ: Writing – review and editing, Funding acquisition, Methodology, Supervision. YE: Writing – review and editing, Funding acquisition, Supervision, Methodology.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. YE was supported by JSPS KAKENHI Grant Number 25K07137. GZ was supported by NSFC General Project No.12171071, and the Natural Science Foundation of Sichuan Province: No. 2023NSFSC0055.

Acknowledgments

The authors would like to thank the reviewers for their valuable comments and suggestions, which greatly helped to improve the quality of this paper. We also express our sincere gratitude to Jinrui Zhou for her assistance with the numerical methods.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Minakowski P, Mucha PB, Peszek J, Zatorska E. Singular Cucker–Smale dynamics, active particles. In: Advances in Theory, Models, and Applications, Vol. 2. Cham: Springer (2019). p. 201–43.

Google Scholar

2. Toner J, Tu Y. Flocks, herds, and schools: a quantitative theory of flocking. Phys Rev E (1998) 58(4):4828–58. doi:10.1103/physreve.58.4828

CrossRef Full Text | Google Scholar

3. Becco C, Vandewalle N, Delcourt J, Poncin P. Experimental evidences of a structural and dynamical transition in fish school. Phys A (2006) 367:487–93. doi:10.1016/j.physa.2005.11.041

CrossRef Full Text | Google Scholar

4. Bialek W, Cavagna A, Giardina I, Mora T, Silvestri E, Viale M, et al. Statistical mechanics for natural flocks of birds. Proc Natl Acad Sci U S A (2012) 109(13):4786–91. doi:10.1073/pnas.1118633109

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Couzin ID, Franks NR. Self-organized Lane formation and optimized traffic flow in army ants. Proc R Soc London, Ser B (2003) 270(1511):139–46. doi:10.1098/rspb.2002.2210

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Partridge BL. The structure and function of fish schools. Sci Am (1982) 246(6):114–23. doi:10.1038/scientificamerican0682-114

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Brambilla M, Ferrante E, Birattari M, Dorigo M. Swarm robotics: a review from the swarm engineering perspective. Swarm Intell (2013) 7(1):1–41. doi:10.1007/s11721-012-0075-2

CrossRef Full Text | Google Scholar

8. Vicsek T, Czirók A, Ben-Jacob E, Cohen I, Shochet O. Novel type of phase transition in a system of self-driven particles. Phys Rev Lett (1995) 75(6):1226–9. doi:10.1103/physrevlett.75.1226

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Reynolds CW. Flocks, herds and schools: a distributed behavioral model. Proc 14th Annu Conf Comput Graph Interact Tech (1987) 25–34. doi:10.1145/37401.37406

CrossRef Full Text | Google Scholar

10. Cucker F, Smale S. Emergent behavior in flocks. IEEE Trans Autom Control (2007) 52(5):852–62. doi:10.1109/tac.2007.895842

CrossRef Full Text | Google Scholar

11. Carrillo JA, Fornasier M, Rosado J, Toscani G. Asymptotic flocking dynamics for the kinetic Cucker–Smale model. SIAM J Math Anal (2010) 42(1):218–36. doi:10.1137/090757290

CrossRef Full Text | Google Scholar

12. Cucker F, Smale S. On the mathematics of emergence. Jpn J Math (2007) 2:197–227. doi:10.1007/s11537-007-0647-x

CrossRef Full Text | Google Scholar

13. Ha SY, Liu JG. A simple proof of the Cucker-Smale flocking dynamics and mean-field limit. Commun Math Sci (2009) 7(2):297–325. doi:10.4310/cms.2009.v7.n2.a2

CrossRef Full Text | Google Scholar

14. Ha SY, Tadmor E. From particle to kinetic and hydrodynamic descriptions of flocking. Kinet Relat Models (2008) 1(3):415–35. doi:10.3934/krm.2008.1.415

CrossRef Full Text | Google Scholar

15. Tadmor E, Tan C. Critical thresholds in flocking hydrodynamics with non-local alignment. le Philos Trans R Soc A (2014) 372(2028):20130401. doi:10.1098/rsta.2013.0401

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Haskovec J, Markou I. Asymptotic flocking in the Cucker–Smale model with reaction-type delays in the non-oscillatory regime. Kinet Relat Models (2020) 13(4):795–813. doi:10.3934/krm.2020027

CrossRef Full Text | Google Scholar

17. Cartabia MR. Cucker-Smale model with time delay. Dyn Syst (2022) 42(5):2409–32. doi:10.3934/dcds.2021195

CrossRef Full Text | Google Scholar

18. Choi YP, Haskovec J. Cucker-Smale model with normalized communication weights and time delay. Kinet Relat Models (2017) 10(4):1011–33. doi:10.3934/krm.2017040

CrossRef Full Text | Google Scholar

19. Erban R, Haskovec J, Sun Y. A Cucker-Smale model with noise and delay. SIAM J Appl Math (2016) 76(4):1535–57. doi:10.1137/15M1030467

CrossRef Full Text | Google Scholar

20. Zhang Z, Yin X, Gao Z. Non-flocking and flocking for the Cucker-Smale model with distributed time delays. J Franklin Inst (2023) 360(12):8788–805. doi:10.1016/j.jfranklin.2022.03.028

CrossRef Full Text | Google Scholar

21. Du L, Han X, Wang Y. Collective behavior for the delayed Cucker-Smale system in a harmonic potential field. Proc Am Math Soc (2024) 152(01):423–34. doi:10.1090/proc/16471

CrossRef Full Text | Google Scholar

22. Choi YP, Haskovec J. Hydrodynamic Cucker-Smale model with normalized communication weights and time delay. SIAM J Math Anal (2019) 51(3):2660–85. doi:10.1137/17m1139151

CrossRef Full Text | Google Scholar

23. Canizo JA, Carrillo JA, Rosado J. A well-posedness theory in measures for some kinetic models of collective motion. Math Models Methods Appl Sci (2011) 21(3):515–39. doi:10.1142/s0218202511005131

CrossRef Full Text | Google Scholar

24. Carrillo JA, Choi YP, Perez SP. A review on attractive–repulsive hydrodynamics for consensus in collective behavior, active particles, volume 1: advances in theory, models, and applications. Cham: Springer (2017). p. 259–98.

Google Scholar

25. Chen Y, Kolokolnikov T. A minimal model of predator–swarm interactions. J R Soc Interf (2014) 11(94):20131208. doi:10.1098/rsif.2013.1208

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Croft S. Modelling collective motion and obstacle avoidance to assess avian collision risk with wind turbines. Heslington: University of York (2014). Ph.D. thesis.

Google Scholar

27. Johansson A, Helbing D, Shukla P. Specification of the social force pedestrian model by evolutionary adjustment to video tracking data. Adv Complex Syst (2007) 10(Suppl. 02):271–88. doi:10.1142/s0219525907001355

CrossRef Full Text | Google Scholar

28. Lin HT, Ros IG, Biewener AA. Through the eyes of a bird: modelling visually guided obstacle flight. J R Soc Interf (2014) 11(96):20140239. doi:10.1098/rsif.2014.0239

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Mecholsky NA, Ott E, Antonsen TM. Obstacle and predator avoidance in a model for flocking. Phys D (2010) 239(12):988–96. doi:10.1016/j.physd.2010.02.007

CrossRef Full Text | Google Scholar

30. Nguyen LTH, Tạ VT, Yagi A. Obstacle avoiding patterns and cohesiveness of fish school. J Theor Biol (2016) 406:116–23. doi:10.1016/j.jtbi.2016.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Aung E, McClure JE, Abaid N. Local interactions in active matter are reinforced by spatial structure. Phys Rev E (2024) 110(6):L062102. doi:10.1103/physreve.110.l062102

PubMed Abstract | CrossRef Full Text | Google Scholar

32. He S, Tadmor E. A game of alignment: collective behavior of multi-species. Ann Inst H Poincaré Anal Non Linaire (2021) 38(4):1031–53. doi:10.1016/j.anihpc.2020.10.003

CrossRef Full Text | Google Scholar

33. Carrillo JA, Castro MJ, Kalliadasis S, Perez SP. High-order well-balanced finite-volume schemes for hydrodynamic equations with nonlocal free energy. SIAM J Sci Comput (2021) 43(2):A828–58. doi:10.1137/20m1332645

CrossRef Full Text | Google Scholar

34. Mavridis CN, Tirumalai A, Baras JS. Learning swarm interaction dynamics from density evolution. IEEE Trans Control Netw Syst (2022) 10(1):214–25. doi:10.1109/tcns.2022.3198784

CrossRef Full Text | Google Scholar

35. Motsch S, Tadmor E. Heterophilious dynamics enhances consensus. SIAM Rev (2014) 56(4):577–621. doi:10.1137/120901866

CrossRef Full Text | Google Scholar

36. Carrillo JA, Delgadino MG, Mellet A. Regularity of local minimizers of the interaction energy via obstacle problems, commun. Math Phys (2016) 343:747–81.

Google Scholar

37. Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phys (1996) 126(1):202–28. doi:10.1006/jcph.1996.0130

CrossRef Full Text | Google Scholar

38. Gottlieb S, Shu CW. Total variation diminishing Runge-Kutta schemes. Math Comp (1998) 67(221):73–85. doi:10.1090/s0025-5718-98-00913-2

CrossRef Full Text | Google Scholar

39. Gottlieb S, Shu CW, Tadmor E. Strong stability-preserving high-order time discretization methods. SIAM Rev (2001) 43(1):89–112. doi:10.1137/s003614450036757x

CrossRef Full Text | Google Scholar

40. Shu CW, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys (1988) 77(2):439–71. doi:10.1016/0021-9991(88)90177-5

CrossRef Full Text | Google Scholar