Abstract
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 [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 byfor and , where denotes the spatial dimension. The functions and represent the density and velocity fields, respectively. The delayed quantities are defined as and , where is a constant denoting the time delay in communication and information processing. The communication kernel is a smooth, radially symmetric function encoding the strength of velocity alignment between agents. The term models a repulsive force field generated by static obstacles in the environment.
In this work, we adopt the fat-tailed communication kernelwhere controls the decay rate of interactions with respect to distance, and the constant is chosen such that , and we model the obstacle force as the negative gradient of a localized Gaussian-type repulsive potential:where and control the strength and spatial scale of the repulsion, respectively, and 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 asthe singular communication kernels like , , 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 functionproposed 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 and denote the density and velocity fields of species . The governing equations are:where the source term captures the nonlocal interactions and external obstacle force, given bywhere , and denote the density and velocity of species , respectively. The delayed terms are given by and , with the communication delay of species . The kernel characterizes the influence of species on : and describe intra-species interactions, while and correspond to inter-species alignment. Existing theoretical analyses have mostly focused on the symmetric interaction case . In our numerical simulations, we will adopt the symmetric interactions in all examples except Example 4.8. The term accounts for local repulsive effects from external obstacles. As in the single-species case, it is modeled by the gradient of a repulsive potential:where is the center of the obstacle, and 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 aswhere for a single-species system and for a two-species system. Then each is conserved over time:Define the total momentum asIn the absence of communication delays and obstacles, the symmetry of the interaction kernels guarantees conservation of total momentum:
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 , divided into uniform cells of width . Let denote the center of cell for , and define the corresponding cell interval by with the velocity defined by . We denote the cell-averaged variables by . The semi-discrete finite volume formulation reads:where are numerical fluxes, and is the source term incorporating nonlocal alignment and obstacle effects. For convenience, we define the spatial operator as
3.1.1 Lax–Friedrichs flux with fifth-order WENO reconstruction
In our scheme, we consider the pressureless Euler-type flux function:
The Lax–Friedrichs numerical flux at the cell interface is defined bywhere and are the fifth-order WENO reconstructed values from the left and right, respectively, and is the maximum absolute eigenvalue of the Jacobian at the interface, defined bywhere denotes the Jacobian matrix of the flux function with respect to the conserved variables , evaluated at the intercell state , and denotes the eigenvalue of a matrix.
Remark 3.1The wave speed in Equation 7 can be chosen locally, as in Equation 8, to reduce numerical diffusion and capture local flow features, or globally as , 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 using a five-point stencil , which forms three three-point substencils:Each substencil produces a third-order approximation . The final reconstructed value at the interface is a weighted sum:where the nonlinear weights are computed byHere, are linear weights, is a small positive parameter, and are smoothness indicators as defined in [37]:For the right-biased value , the same procedure is applied to the reversed stencil . 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 is given bywhere , , and denotes the uniform cell width of the finite volume mesh. The function 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:
Using these, the discrete source term can be reformulated as
We first focus on the efficient computation of the alignment contribution 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 and is given bywhere and denote the discrete Fourier transform and its inverse, respectively, and the multiplication is performed element-wise in Fourier space. This approach reduces the computational complexity of evaluating each convolution from to , which is especially beneficial for large-scale simulations.
The second component of the source term, , accounts for obstacle-induced repulsion potential and is modeled as the discrete gradient of a smooth obstacle potential field:where is given by Equation 3, and the gradient 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.where is the spatial operator defined in Equation 6.
The SSP-RK3 scheme advances from to via the following three-stage procedure:
For stability, the time step is chosen to satisfy the CFL conditionwhere is a prescribed constant (e.g., ) and local wave speed 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 , which is discretized into uniform cells with mesh sizes , . Denote the center of cell by , where for , . The control volume is defined as the rectangular cell centered at :
Let the cell averages of the conserved variables be defined as
The semi-discrete finite volume scheme is given bywhere are numerical fluxes in the - and -directions, and is the discrete source term. For convenience, we define the two-dimensional spatial operator as
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:
The physical fluxes in the - and -directions are given by
The Lax–Friedrichs numerical fluxes at cell interfaces are:
The reconstructed states , 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 :
Remark 3.2Alternatively, a global maximum wave speedwhich 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 consists of a nonlocal alignment interaction and an obstacle-induced repulsion force:where the delayed fields are defined by To facilitate efficient computation of the nonlocal alignment term in Equation 14, we introduce the auxiliary quantities:Using these, the source term Equation 14 can be written more compactly as:We first focus on the alignment contribution 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 is modeled as the negative gradient of a smooth potential field , evaluated at the grid point . It is discretized using central finite differences: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 :where is the spatial operator defined in Equation 12. The time step is constrained by the CFL condition:where the wave speed defined in Equation 13 and .
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
(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:
Proof. For convenience, we denote the first (density) component of the update operator defined in Equation 6 by , i.e.,where denotes the first component of the numerical flux defined in Equation 7, corresponding to the mass flux.
Let denote the cell average of density at time . The SSP-RK3 updates for density read:
Define the total mass at each stage as
By using the periodic boundary condition, we havesumming Equations 16–18 over all , we find:
Thus, the scheme is exactly mass conservative.
3.3.2 Momentum conservation
(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:
Proof. We denote the second (momentum) component of the update operator defined in Equation 6 by , i.e.,where denotes the second component of the numerical flux defined in Equation 7, and is the discrete source term defined in Equation 9.
Let be the cell average of momentum at time . The SSP-RK3 scheme updates for momentum are given by:
Define the total momentum at each stage as
Under periodic boundary conditions, it is easy to see that:
Therefore, the total contribution from the update operator reduces to the source term:
In the special case when and , the source term defined in Equation 9 simplifies towhere . Since is symmetric in , and is antisymmetric, the double sum satisfies
This implies
Therefore, from Equation 23 we have
Summing Equations 19–21 over all , in view of Equation 22 and Equation 24, we see that the total momentum is conserved during the time update, i.e.,
Thus we obtain the conclusion.
Remark 3.5When delay is present, the source term involves past states , breaking the symmetry required for momentum conservation. Similarly, nonzero obstacle forces 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].where is the total mass and 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
4 Numerical experiments
This section presents numerical simulations. Throughout all simulations, the interaction kernel is chosen as , where the constant is chosen such that , 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 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 with :where the normalization constant is chosen such that . We vary the number of cells as and set the time step according to , where . We compute the cell averages and compare the numerical solution at time with the exact solution and .As shown in Table 1, the errors for and converge with an order approaching 5 as the grid is refined, where and . This confirms that the scheme achieves fifth-order spatial accuracy, consistent with the WENO reconstruction employed in the spatial discretization.
TABLE 1
| Order | Order | |||
|---|---|---|---|---|
| 10 | 1.5205E-02 | – | 1.5205E-02 | – |
| 20 | 7.3828E-04 | 4.36 | 7.3828E-04 | 4.36 |
| 40 | 2.3166E-05 | 4.99 | 2.3166E-05 | 4.99 |
| 80 | 7.2296E-07 | 5.00 | 7.2299E-07 | 5.00 |
Accuracy test for Example 4.1 at .
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 with :where , and the normalization constant is chosen such that . We fix the number of spatial cells at and vary the time step by adjusting the CFL number defined in Equation 11 accordingly.As shown in Table 2, the 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
| CFL | Order | Order | ||
|---|---|---|---|---|
| 0.8 | 1.8514E-07 | – | 1.8514E-07 | – |
| 0.4 | 2.9040E-08 | 2.67 | 2.9040E-08 | 2.67 |
| 0.2 | 3.7531E-09 | 2.95 | 3.7531E-09 | 2.95 |
| 0.1 | 4.6885E-10 | 3.00 | 4.6885E-10 | 3.00 |
Accuracy test for Example 4.2 at .
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 aswhere the normalization constant is chosen so that , and controls the velocity amplitude. The simulation is conducted on a periodic domain (i.e., ), with spatial resolution .In the third scenario, we introduce an obstacle modeled as the gradient of a repulsive potential centered at , as defined in Equation 3:with repulsion parameters and .Figure 1 shows the time evolution of total mass, momentum, and energy fluctuation (defined in Remark 3.6) for each configuration, with the top, middle, and bottom panels displaying mass , momentum , and , 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 , and the bottom panel shows the total momentum . 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 . In panels (a) and (b), approaches to zero by , while in panel (c) it tends to zero around . 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 and obstacle force on the collective behavior and regularity of solutions in the single-species system Equation 1. The initial data is given bywhere the normalization constant is chosen such that , and the parameter controls the amplitude of the initial velocity. The simulation is performed on the periodic domain , i.e., , with the number of cell .
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: and . For each case, we vary the delay parameter to examine its influence on the long-time behavior.Figure 2 presents the density and velocity at final time for . 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 . When , the solution remains smooth throughout the simulation. However, as increases, the density begins to concentrate more sharply. In particular, for , we observe that the solution develops a near-singular profile: the density becomes highly concentrated around , 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 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 3
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 , as defined in Equation 3:with repulsion parameters and . The initial velocity amplitude is set to , and we examine the impact of varying the delay parameter .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: for , for , and for . 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
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:with and constants . The normalization constants are chosen such that . The simulation is conducted using a periodic domain discretized with the number of cell and a CFL number of 0.7, running from to .
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 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 and , 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 , 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 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 , 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 , beyond which the system exhibits a loss of regularity and a breakdown of coordinated motion.
FIGURE 5
Remark 4.1A 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 , as defined in Equation 5:with repulsion parameters and , . The obstacle force is taken to be the same for both species, i.e., . Figure 6 illustrates the evolution of the two-species density and velocity field in the presence of a fixed obstacle force and varying delay parameters , 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 increases, the formation of these concentrated structures occurs earlier in time: at for , for , and for . Moreover, the aggregation regions become more sharply localized and shift closer to the obstacle center as 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
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 and . Hence, species 1 applies the same exponent for both self- and cross-interactions, while species 2 interacts with species 1 through the smaller exponent . This asymmetry highlights the effect of asymmetric inter-species interactions on the collective dynamics.For the different initial conditions chosen for the two species, we observe that when and 0.5, the solution behavior is similar to that in Example 4.6. However, when , Examples 4.6 and 4.8 show that the aggregation levels of and 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
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 with , uniformly discretized using grid points. The CFL number is set to 0.012. The initial conditions are specified as follows:where the normalization constant is selected such that and the initial velocity parameter . We investigate the long-time evolution of the system for various delay values .
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 , the two-dimensional single-species system remains smooth and appears to converge toward a steady state in the absence of delay , 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 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 , 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 .These numerical observations indicate the possible existence of a critical delay threshold , 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 9
Example 4.10(Effects of delay with obstacle). We introduce an obstacle force modeled as the gradient of a repulsive potential centered at the obstacle location, as defined in Equation 3:In our numerical example, we consider a single obstacle located at , with repulsion parameters and .Figure 10 presents the full temporal evolution of the one-species density field under zero delay , while Figure 11 displays the final-time density distributions for increasing delays , 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 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 for to for , for , and for , 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 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 sufficiently small in order to satisfy the flocking condition; (b) For the fat-tailed kernel with , any 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 with . 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 11
Remark 4.2This work focuses on fat-tailed kernels. We also test compactly supported kernels (such as ), 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.
Statements
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.
MinakowskiPMuchaPBPeszekJZatorskaE. Singular Cucker–Smale dynamics, active particles. In: Advances in Theory, Models, and Applications, Vol. 2. Cham: Springer (2019). p. 201–43.
2.
TonerJTuY. Flocks, herds, and schools: a quantitative theory of flocking. Phys Rev E (1998) 58(4):4828–58. 10.1103/physreve.58.4828
3.
BeccoCVandewalleNDelcourtJPoncinP. Experimental evidences of a structural and dynamical transition in fish school. Phys A (2006) 367:487–93. 10.1016/j.physa.2005.11.041
4.
BialekWCavagnaAGiardinaIMoraTSilvestriEVialeMet alStatistical mechanics for natural flocks of birds. Proc Natl Acad Sci U S A (2012) 109(13):4786–91. 10.1073/pnas.1118633109
5.
CouzinIDFranksNR. Self-organized Lane formation and optimized traffic flow in army ants. Proc R Soc London, Ser B (2003) 270(1511):139–46. 10.1098/rspb.2002.2210
6.
PartridgeBL. The structure and function of fish schools. Sci Am (1982) 246(6):114–23. 10.1038/scientificamerican0682-114
7.
BrambillaMFerranteEBirattariMDorigoM. Swarm robotics: a review from the swarm engineering perspective. Swarm Intell (2013) 7(1):1–41. 10.1007/s11721-012-0075-2
8.
VicsekTCzirókABen-JacobECohenIShochetO. Novel type of phase transition in a system of self-driven particles. Phys Rev Lett (1995) 75(6):1226–9. 10.1103/physrevlett.75.1226
9.
ReynoldsCW. Flocks, herds and schools: a distributed behavioral model. Proc 14th Annu Conf Comput Graph Interact Tech (1987) 25–34. 10.1145/37401.37406
10.
CuckerFSmaleS. Emergent behavior in flocks. IEEE Trans Autom Control (2007) 52(5):852–62. 10.1109/tac.2007.895842
11.
CarrilloJAFornasierMRosadoJToscaniG. Asymptotic flocking dynamics for the kinetic Cucker–Smale model. SIAM J Math Anal (2010) 42(1):218–36. 10.1137/090757290
12.
CuckerFSmaleS. On the mathematics of emergence. Jpn J Math (2007) 2:197–227. 10.1007/s11537-007-0647-x
13.
HaSYLiuJG. A simple proof of the Cucker-Smale flocking dynamics and mean-field limit. Commun Math Sci (2009) 7(2):297–325. 10.4310/cms.2009.v7.n2.a2
14.
HaSYTadmorE. From particle to kinetic and hydrodynamic descriptions of flocking. Kinet Relat Models (2008) 1(3):415–35. 10.3934/krm.2008.1.415
15.
TadmorETanC. Critical thresholds in flocking hydrodynamics with non-local alignment. le Philos Trans R Soc A (2014) 372(2028):20130401. 10.1098/rsta.2013.0401
16.
HaskovecJMarkouI. Asymptotic flocking in the Cucker–Smale model with reaction-type delays in the non-oscillatory regime. Kinet Relat Models (2020) 13(4):795–813. 10.3934/krm.2020027
17.
CartabiaMR. Cucker-Smale model with time delay. Dyn Syst (2022) 42(5):2409–32. 10.3934/dcds.2021195
18.
ChoiYPHaskovecJ. Cucker-Smale model with normalized communication weights and time delay. Kinet Relat Models (2017) 10(4):1011–33. 10.3934/krm.2017040
19.
ErbanRHaskovecJSunY. A Cucker-Smale model with noise and delay. SIAM J Appl Math (2016) 76(4):1535–57. 10.1137/15M1030467
20.
ZhangZYinXGaoZ. Non-flocking and flocking for the Cucker-Smale model with distributed time delays. J Franklin Inst (2023) 360(12):8788–805. 10.1016/j.jfranklin.2022.03.028
21.
DuLHanXWangY. Collective behavior for the delayed Cucker-Smale system in a harmonic potential field. Proc Am Math Soc (2024) 152(01):423–34. 10.1090/proc/16471
22.
ChoiYPHaskovecJ. Hydrodynamic Cucker-Smale model with normalized communication weights and time delay. SIAM J Math Anal (2019) 51(3):2660–85. 10.1137/17m1139151
23.
CanizoJACarrilloJARosadoJ. A well-posedness theory in measures for some kinetic models of collective motion. Math Models Methods Appl Sci (2011) 21(3):515–39. 10.1142/s0218202511005131
24.
CarrilloJAChoiYPPerezSP. 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.
25.
ChenYKolokolnikovT. A minimal model of predator–swarm interactions. J R Soc Interf (2014) 11(94):20131208. 10.1098/rsif.2013.1208
26.
CroftS. Modelling collective motion and obstacle avoidance to assess avian collision risk with wind turbines. Heslington: University of York (2014). Ph.D. thesis.
27.
JohanssonAHelbingDShuklaP. Specification of the social force pedestrian model by evolutionary adjustment to video tracking data. Adv Complex Syst (2007) 10(Suppl. 02):271–88. 10.1142/s0219525907001355
28.
LinHTRosIGBiewenerAA. Through the eyes of a bird: modelling visually guided obstacle flight. J R Soc Interf (2014) 11(96):20140239. 10.1098/rsif.2014.0239
29.
MecholskyNAOttEAntonsenTMJr. Obstacle and predator avoidance in a model for flocking. Phys D (2010) 239(12):988–96. 10.1016/j.physd.2010.02.007
30.
NguyenLTHTạVTYagiA. Obstacle avoiding patterns and cohesiveness of fish school. J Theor Biol (2016) 406:116–23. 10.1016/j.jtbi.2016.07.017
31.
AungEMcClureJEAbaidN. Local interactions in active matter are reinforced by spatial structure. Phys Rev E (2024) 110(6):L062102. 10.1103/physreve.110.l062102
32.
HeSTadmorE. A game of alignment: collective behavior of multi-species. Ann Inst H Poincaré Anal Non Linaire (2021) 38(4):1031–53. 10.1016/j.anihpc.2020.10.003
33.
CarrilloJACastroMJKalliadasisSPerezSP. High-order well-balanced finite-volume schemes for hydrodynamic equations with nonlocal free energy. SIAM J Sci Comput (2021) 43(2):A828–58. 10.1137/20m1332645
34.
MavridisCNTirumalaiABarasJS. Learning swarm interaction dynamics from density evolution. IEEE Trans Control Netw Syst (2022) 10(1):214–25. 10.1109/tcns.2022.3198784
35.
MotschSTadmorE. Heterophilious dynamics enhances consensus. SIAM Rev (2014) 56(4):577–621. 10.1137/120901866
36.
CarrilloJADelgadinoMGMelletA. Regularity of local minimizers of the interaction energy via obstacle problems, commun. Math Phys (2016) 343:747–81.
37.
JiangGSShuCW. Efficient implementation of weighted ENO schemes. J Comput Phys (1996) 126(1):202–28. 10.1006/jcph.1996.0130
38.
GottliebSShuCW. Total variation diminishing Runge-Kutta schemes. Math Comp (1998) 67(221):73–85. 10.1090/s0025-5718-98-00913-2
39.
GottliebSShuCWTadmorE. Strong stability-preserving high-order time discretization methods. SIAM Rev (2001) 43(1):89–112. 10.1137/s003614450036757x
40.
ShuCWOsherS. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys (1988) 77(2):439–71. 10.1016/0021-9991(88)90177-5
Summary
Keywords
Cucker–Smale model, time delays, hyperbolic systems, collective behavior, flocking, obstacle avoidance
Citation
Zheng T, Zhou G and Enatsu Y (2025) Hydrodynamic Cucker-Smale model with time delay and obstacle avoidance. Front. Phys. 13:1657927. doi: 10.3389/fphy.2025.1657927
Received
02 July 2025
Accepted
22 September 2025
Published
17 October 2025
Volume
13 - 2025
Edited by
Ryosuke Yano, Tokio Marine dR Co., Ltd., Japan
Reviewed by
Ndolane Sene, Cheikh Anta Diop University, Senegal
Poornachandra Sekhar Burada, Indian Institute of Technology Kharagpur, India
Updates
Copyright
© 2025 Zheng, Zhou and Enatsu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guanyu Zhou, wind_geno@live.com
Disclaimer
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.