We focus our study on an approximate regime of a lattice gas model. The key feature, as previously noted, is the exclusion principle, which prevents particles from overlapping or occupying the same cell, regardless of the specific rules governing the transition probabilities. This characteristic makes it particularly challenging to derive an approximate recurrence relation or equation of motion, as the occupation states of neighboring lattice cells are highly correlated. Nevertheless, we develop a recurrence relation for the cell occupation that remains consistent with the asynchronous updating scheme characteristic of lattice gas dynamics with nearest-neighbor interactions. Accordingly, we propose the following recurrence relation for a given species q: ρ t + 1 q r = 1 − ρ t r ∑ 〈 r ′ 〉 ρ t q r ′ P r r ′ → r q + ρ t q r P r r → r q + ∑ 〈 r ′ 〉 ρ t r ′ P r r → r ′ q , (3) where ρ t ( r ) ≡ ∑ k = 1 m ρ t ( k ) ( r ) is the total density of particles at a given cell r at time step t . The proposed equation decomposes the update process into two key contributions: (1) particle inflow into the cell and (2) particle persistence or outflow from the cell. In the lattice gas dynamics, the occupation state is binary, with ρ t ( q ) ( r ) = 1 if the cell is occupied and ρ t ( q ) ( r ) = 0 if it is vacant. However, in the mean-field (MF) approximation, ρ t ( r ) represents an average density field and does not directly reflect the actual number of particles n in the system. To bridge this gap and improve consistency between the MF and MC descriptions, we generalized the normalization condition of the MF recurrence relation. Specifically, we modified the global constraint so that the sum over all lattice sites scales as a power of n : ∑ r ρ t r = n β , where β is a tunable exponent. The introduction of β allows us to adjust the MF approximation to better reproduce the behavior observed in MC simulations. By treating β as a fitting parameter, we aim to compensate for the differences between the averaged MF treatment and the discrete, stochastic nature of the MC dynamics.

The first term on the right-hand side of Equation 3, the factor 1 − ρ t ( r ) , reflects the hard-body exclusion principle, which dictates that a particle can only occupy an empty cell. This factor multiplies the sum of possible sources of inflow of particles of species q from the neighboring cells of r , i.e., the notation ⟨ r ′ ⟩ . Each neighboring source of particles depends on the density of the species q at the neighboring site ρ t ( q ) ( r ′ ) and the transition probability P r r ′ → r ( q ) , reflecting its likelihood to hop from r ′ to r . The second term on the right-hand side of Equation 3 describes the contribution of particles already present in r . The first factor, ρ t ( q ) ( r ) , accounts for the actual occupation of it, whilst the factor inside the large brackets corresponds to a composition of two possibilities. The first possibility is the probability P r r → r ( q ) of the particle at r remaining in place, which contrasts with the second possibility that stands as the probability of the particle at r trying to move to an occupied neighboring cell r ′ .

With this framework, we explicitly separate the exclusion interaction, reflected in the terms 1 − ρ t ( r ) and ρ t ( r ) , from the transition probabilities, which can incorporate additional interaction mechanisms such as long-range fields or cooperative dynamics. As an example, Equation 3 can reproduce even deterministic models of particles, such as proposed in [21], which could allow us to describe the model proposed by Cividini et al., for instance, by a proper choice of parameters such as fixing m = 2 with one species moving eastbound and the other species moving northbound with transition probabilities equal to one towards their respective target directions.

We can rewrite Equation 3 as follows ρ t + 1 q r = ∑ ⟨ r ′ ⟩ ω t q r , r ′ − ω t q r ′ , r − ρ t r ω t q r , r ′ − ρ t r ′ ω t q r ′ , r + ρ t q , (4) where ω t ( q ) ( r , r ′ ) ≡ ρ t ( q ) ( r ′ ) P r r ′ → r ( q ) is the flow rate of particles of species q entering the cell at r from the cell r ′ . Alternatively, we could write Equation 4 as ρ t + 1 q r = Φ t q r − Ψ t q r + ρ t q , (5) where Φ t q r ≡ ∑ ⟨ r ′ ⟩ ω t q r , r ′ − ω t q r ′ , r (6) is the net flow rate of particles of species q into the cell at r and Ψ t q r ≡ ∑ ⟨ r ′ ⟩ ρ t r ω t q r , r ′ − ρ t r ′ ω t q r ′ , r (7) is the frustrated net flow rate of particles of species q into the cell at r . The frustrated net flow rate measures the particle interaction, i.e., measures the amount of q particles unable to access or leave the cell at r due to the its occupation and of its neighbouring cells Equations 5-7 provide an alternative representation that offers a more physical perspective on the mean-field dynamics.

When applying Equation 1 into the lattice gas recurrence relation given by Equation 3, we end up with the specific recurrence relation ρ t + 1 q r = 1 − ρ t r ∑ 〈 r ′ 〉 ρ t q r ′ p − α Δ r • û q r ′ + ρ t q r 1 − 4 p + ∑ 〈 r ′ 〉 ρ t r ′ p + α Δ r • û q r . (8)

After some algebra (see Supplementary Material A) we obtain the following partial differential equation ∂ ρ q ∂ t = c 1 ′ ∇ 2 ρ q − c 2 ′ ∇ • ρ q û q − c 3 ′ ∑ k ≠ q ρ k ∇ 2 ρ q − ρ q ∇ 2 ρ k + c 4 ′ ∑ k = 1 m ρ k ∇ • ρ q û q + ρ q ∇ ρ k • û q , (9) where c 1 ′ , c 3 ′ ∝ lim τ , ϵ → 0 p ϵ 2 τ and c 2 ′ , c 4 ′ ∝ lim τ , ϵ → 0 α ϵ τ .

In Equation 9, we notice that the constants c 3 ′ and c 4 ′ are species’ coupling factors and the mathematical path we took from Equation 3 makes them not to depend on each possible pair of species. However, a more general approach would consider these constants as dependent on each possible coupling pair of species, which we will not going to explore in this work.

To study how good a fit is the mean-field approximation is to the particle dynamics, we investigate how the spatial distributions of both methods differ from each other for different sets of the model parameters such as the total number of particles ( n ) , the total number of species ( m ) , and the bias level of movement ( α ) .

The simplest scenario we examined involves only a single species. As we defined previously by Equation 10, a single species tends to move according to the static floor field û = e x , with its initial distribution having μ x ( 1 ) = l / 4 and μ y ( 1 ) = l / 2 with σ x ( q ) = σ y ( q ) = l / 32 . In these simulations, we used β = 1 .

In Figure 3, we present plots of the marginal distribution of particles along x , where circles represent the Monte Carlo (MC) simulation results and solid lines correspond to predictions from the recurrence relation (mean-field approximation). We examine these distributions for various particle counts, ranging from n = 64 to n = 448 in increments of Δ n = 64 and shown at four distinct time steps: (a) t = 50 , (b) t = 250 , (c) t = 450 , and (d) t = 4950 . For this study, we used α = 0.15 , which corresponds to agents with mid range level of impetus as we already shown in [7, 8].

As a general remark, we notice a qualitatively good agreement between the two methods and that the distributions in more densely populated systems shows more positive skewness of the probability density function (PDF) in comparison to systems with fewer particles, revealing the exclusion effects. These exclusion effects are a result of agents leading the pack acting as obstacles for the particles in the bulk. As time passes, we notice that the distributions present increasing dispersion with the “front tail” maintaining synchronicity for different values of n . It is important to disclose that in the study that produced Figure 3, the recurrence relation exhibited numerical instability for n ≥ 512 , which the corresponding data is omitted here and throughout the manuscript for conciseness. To understand the source of the instabilities, note that Equation 8 includes the term 1 − ρ t ( r ) , which, in the context of discrete particle dynamics, only takes values 0 or 1, as previously explained. In this regime, the factor 1 − ρ t ( r ) behaves in a stable and well-defined manner. However, under the mean-field approximation, increasing either the total density N / V or the parameter α can lead to locally high average densities. As a result, 1 − ρ t ( r ) may become negative, causing ρ t + 1 ( r ) to also take on unphysical negative values. When this occurs, the normalization of the distribution becomes unstable, as the sum involves both large positive and negative terms—sometimes yielding normalized values with magnitudes as high as 1 0 20 . Thus, providing a quantitative indication that there exists a certain range of the parameters of the model beyond which the approximate method becomes unreliable and serves as justification for the assumption that a tuning parameter, β , might be necessary to also preserve the mean-field stability. The reader is advised to check the supplementary material to check the form of marginal distributions of particles right before numerical collapse.

Another key parameter we investigated is the bias movement level, α , which quantifies the particles’ tendency to move in their species’ preferred direction. For this analysis, we fixed n = 192 , corresponding to a generally low-density scenario as shown, while ensuring sufficient particle interactions to produce the positive skewness observed in the distribution.

In the lower panels of Figure 3 (plots e to h), we present the marginal distributions at time steps (e) t = 50 , (f) t = 250 , (g) t = 450 , and (h) t = 4950 for α = 0.050 , 0.100 , 0.150 , and 0.200. As expected, distributions with higher α values exhibit faster drift velocities in the + x direction. Additionally, increasing α amplifies the skewness, similar to the effect observed with higher n . Notably, at longer simulation times, larger α values lead to greater discrepancies between the two methods, as seen in Figure 3h.

Given the phenomena observed in Figure 3, we notice that systems with a greater total number of particles ( n ) tend to show less agreement, even if qualitatively their shapes follow the same pattern. However, this disparity is not observed when changing the level of movement bias ( α ) shown in Figure 3. With this in mind, we show, in Figure 4, how the exponent of normalization, β , influences the shape of the marginal distribution of x .

Figure 4a illustrates the influence of β on the distribution profile for a system with parameters m = 1 , n = 128 , and α = 0.249 . As discussed in the previous section, the normalization constraint in the recurrence relation can be interpreted through an ansatz in which the number of particles is raised to a power β , with β expected to be close to 1, if not exactly 1. As such, this parameter serves to tune the agreement between methods at the same time, it allows for systems with number of particles expected to show numerical collapse to evolve in a well-behaved manner.

To better measure the agreement between the two methods and the range of parameters where the mean-field shows numerical stability, we compared both methods using the time average of the squared difference of the spatial entropies of each method, i.e., Δ S 2 ̄ ≡ 1 t max ∑ t = 1 t max S t MC − S t MF 2 , (11) where S t = ∑ r ρ t ( r ) ln ρ t ( r ) . The loss function introduced in Equation 11 is one among several possible choices for quantifying the discrepancy between the MC and MF results. However, as we demonstrate below, it proved to be particularly suitable for optimization purposes. For comparison, we also tested an alternative loss function based on the squared differences of local densities, as illustrated in Supplementary Figure S2 and defined in its caption in the Supplementary Material. This alternative metric failed to produce a well-behaved optimization landscape, which made it ineffective for identifying optimal values of β . For this reason, we opted not to pursue it further.

In Figure 4b, we show Δ S 2 ̄ as a function of the exponent β for different values of n and for a fixed simulation time of t max = 1 0 4 time steps. We studied the range of β ∈ [ 0.5 : 1.5 ] . In this range, we observed numerical instability of the recurrence relation in all curves shown for values of β greater than certain values, which are reflected by the discontinuity of the curves. More specifically, greater n shows lower β ’s where the numerical collapse occurs. However, within the range of numerical stability, we observe that, for n ≤ 256 , methods show an optimal agreement at β ≡ β c > 1 , shown by the local minima which become closer to 1 as n increases. For n ≥ 256 , we notice that the numerical instability of the MF makes the lowest stable β assume the stance of optimal value, i.e., it becomes the value that minimizes Δ S 2 ̄ .

Figure 4c completes the analysis for m = 1 , by showing the linear dependence of the optimal value β c on the number of agents. The data is fitted with a linear function a x + b , yielding a = − 0.000588 ± 5.3 × 1 0 − 5 and b = 1.17 ± 0.017 . Such behavior will be different of upper values of m as it is previously presented in Figure 4d. However, we must first study the particle distribution for these cases under varying parameters, which will be performed in the next sections.

The case m = 2 is the simplest case within our proposed framework of initial conditions (see Figure 2) where two species of particles move in counterflow, which is a scenario vastly studied in the literature, once it can describe a two species systems with opposing responses to an external field such as oppositely charged coloids [9, 20] (external electrical field) or chromatographic columns [10, 26–29] (gravitational field), to cite a few. For this case, initial distributions assume μ x ( 1 ) = l / 4 and μ y ( 1 ) = l / 2 , μ x ( 2 ) = 3 l / 4 and μ y ( 2 ) = l / 2 as we defined in the previous section. Static floor fields are then û ( 1 ) = e x and û ( 2 ) = − e x as stated by Equation 10, which still carries a rather trivial symmetry in the y -direction (diffusive profile) making observations of the marginal distribution on the x coordinate rather more interesting because is where the confrontation occurs.

In Figure 5, we show the marginal distributions of x of the MC simulations with circles and lines for the MF approach. We studied two values of the total number of particles shown simultaneously: n = 128 shown in purple and n = 256 shown in green for the time steps (a) t = 50 , (b) t = 250 , (c) t = 450 , and (d) t = 4950 , where both species appear as different peaks with same color for each case. In this first case, we used α = 0.15 and β = 1.0 , and we observe that, similarly to the m = 1 case, a greater total number of particles makes the distributions more asymmetric due to the exclusion interaction of particles of the same species. We do not observe, however, a clear influence of the opposing species on the shape of the distribution. We can interpret this “lack” of interaction deformity occurring because of the exclusion effects of the nearest neighbor interaction, suggesting that a sort of shield is formed by the same-species particles in a low-density regime. This warrants further investigation in future work.

While both methods show reasonably good agreement in Figure 5, we find that the agreement can be further improved by optimizing β . Revisiting Figure 4d and extending our earlier analysis for m = 1 to m = 2 , we observe that the critical value β c minimizing Δ S 2 ̄ follows a power-law dependence, a x − b . For m = 2 , the fit yields a = 3.3 ± 0.11 and b = 0.193 ± 0.0065 , contrasting with the linear behavior found for m = 1 . This indicates that increasing the model complexity (from m = 1 to m = 2 ) modifies the parameter dependence. However, no significant changes are observed beyond m = 2 up to m = 4 , as seen in Figure 4D—a point we will address in detail later. Importantly, we note that lower-density cases do not provide the best agreement between the two methods, revealing non-trivial model-specific characteristics that lack complete theoretical understanding.

The last case we discuss in this work is the case with four species moving along the x − and y − directions. More precisely, their initial distributions have mean value: μ x ( 1 ) = l / 4 and μ y ( 1 ) = l / 2 , μ x ( 2 ) = 3 l / 4 and μ y ( 2 ) = l / 2 , μ x ( 3 ) = l / 2 and μ y ( 3 ) = l / 4 , and μ x ( 4 ) = l / 2 and μ y ( 4 ) = 3 l / 2 as we expected to be by our previous definition. Similarly, the static floor fields assume the form û ( 1 ) = e x , û ( 2 ) = e y , û ( 3 ) = − e x , and û ( 4 ) = − e y (see Equation 10). Similarly to the two prior cases, having an initial condition with radial symmetry about the lattice center, also makes the four species present a symmetric “cross section” of the distribution, i.e., the distribution on the direction perpendicular to û ( q ) shows no asymmetry. Thus, studying the marginal distribution of x of the species q = 1 or q = 3 would be no different from studying the marginal distribution of y of the species q = 2 and q = 4 given the system’s initial arrangement and static floor fields (see Figure 2).

With this in mind, Figure 6 displays the particle distribution ρ ( 1 ) ( x ) for different values of n using α = 0.249 at four time points: (a) t = 50 , (b) t = 250 , (c) t = 450 , and (d) t = 4950 . Our analysis reveals two key observations:

(1) First, the discrepancy between methods becomes more pronounced for larger n values in this regime of highly driven agents. Despite this disagreement, the MC simulations demonstrate the formation of transient condensates for n = 320 , as evidenced by the emerging peak at t = 250 .

Building on our analysis of the m = 4 case, Figure 6 presents the evolution of ρ ( 1 ) ( x ) distributions for various α values in a system with n = 256 particles. We also examine the same four time points: (e) t = 50 , (f) t = 250 , (g) t = 450 , and (h) t = 4950 . The results reveal that while the wave fronts remain synchronized between MC simulations and mean-field theory, the MC distributions exhibit significantly greater dispersion. Returning to Figure 4d, we apply the same minimization procedure for Δ S 2 ̄ used for the m = 1 and m = 2 cases. For m = 4, we obtain an excellent power-law fit identical in form to the m = 2 case, with fitted parameters a = 3.04 ± 0.23 and b = 0.176 ± 0.014.

It is important to note that, despite some differences in the evolution of particle profiles, a qualitatively similar behavior is observed in both cases when the systems are initialized with a normal distribution of particles in the environment. The interactions among particles—especially the interaction with counterflowing particles—tend to deform the initial distribution shape. Over time, however, after several rounds of confrontation and interaction, the particles begin to reorganize themselves, ultimately restoring a Gaussian-like behavior.

This transition from one Gaussian state to another is illustrated in Figure 7. In Figures 7a–e, we show the results from MC simulations, while Figure 7f corresponds to the mean field approximation. All plots refer to the case of m = 4 species, n = 256 particles, and interaction parameter α = 0.249 .

In Figure 7a, we observe the particle distribution at t = 50 . We observe that the initial Gaussian behavior is disrupted due to particle interactions. The red and blue curves represent normal and log-normal fits, respectively, with the latter providing a more accurate fit at this stage. The inset plot shows the same trend on a semi-log scale. By t = 150 , the distribution becomes strongly non-Gaussian, with a pronounced peak emerging, as shown in Figure 7b.

At t = 250 , the shape of the distribution becomes nearly triangular (Figure 7c), signaling the beginning of a return to Gaussianity. A good Gaussian fit is eventually recovered in Figure 7d, and a noisy but stable Gaussian distribution is observed at t = 400 , as shown in Figure 7e. This reflects the establishment of a steady state after extensive particles interactions.

Finally, in Figure 7f, we analyze the mean-field approximation. Unlike the MC simulations, which display the full evolution of the particle profile, we summarize the results using the coefficient of determination R 2 , which quantifies the quality of a normal distribution fit over time. These values are represented by circles. Initially, we observe a sharp decline in R 2 , indicating the breakdown of Gaussianity (illustrated in the inset). As t increases, the system goes through a transient (points inside blue rectangle) phase before recovering a Gaussian profile, again shown in the inset. This confirms that the same Gaussian-to-Gaussian behavior occurs in the mean-field scenario, in agreement with the MC simulations presented in Figures 7a–e.