The emergence of spatiotemporal structures in suspensions of active particles is an important topic in the current research in nonequilibrium statistical physics (Romanczuk et al., 2012; Marchetti et al., 2013; Costanzo et al., 2014; Grossmann et al., 2014; O’Keeffe et al., 2017; Granek et al., 2024). Interestingly, in such settings, order in terms of separation of different sorts of particles may emerge other than by repelling or attracting forces but via the mutual control of motility.

A particularly striking example is the demixing of two different strains of Escherichia coli bacteria with reciprocal motility interaction, studied experimentally and in a reaction-diffusion model by Curatolo et al. (2020). A central part of the mechanism of particle separation was a mutual amplification of the diffusion strength caused by the presence of the respective other bacteria. This could be mimicked in the paper by Schimansky-Geier et al. (2021), which presented, simulated, and studied analytically a minimal model in which two populations of overdamped Brownian particles mutually control their strength of diffusivity (the noise intensity of driving fluctuations increases with the concentration of the respective other sort of particles). The dependence of the diffusion coefficient has to be strong enough (qualitatively and quantitatively) to observe a demixing of the two population densities.

Here, we study a variant of the model suggested by Schimansky-Geier et al. (2021) and demonstrate a novel bistability between a well-mixed and a demixed state that results in a system with finite particles in stochastic transitions between an ordered and a disordered state. This result might be interesting and also observable experimentally when the reciprocal motility interaction is particularly tailored.

Our paper is organized as follows. In the next section, we introduce the studied model, focusing on the novel aspect, the sigmoidal diffusion coefficient. In Section 3, we derive the conditions under which a single demixed state, a single homogeneous (well-mixed) state, and bistability between both states can be observed in the model and present the system bifurcation diagram. In Section 4, we show and discuss the results of simulations in the bistable parameter regime, in which we measure the mean times that the system spends in the two bistable states and inspect how these times depend on the bifurcation parameter. We conclude in Section 5 with a discussion of our results.

We consider two populations of Brownian particles, A and B, each of size N , that diffuse on an interval [ − 1,1 ] with reflecting boundary conditions at x = ± 1 according to d x A , i d t = f p ̃ B x A , i , t   ξ A , i t , d x B , j d t = f p ̃ A x B , j , t   ξ B , j t . (1) Here i , j = 1 , … , N , ξ A , i ( t ) and ξ B , j ( t ) are white Gaussian noise sources with 〈 ξ F , i ( t ) ξ G , i ( t ′ ) 〉 = δ F , G δ i , j δ ( t − t ′ )   w i t h   F , G ∈ A , B , and the equations are interpreted in the Ito sense (Gardiner, 1985). The intensity by which these fluctuations enter at a point x depends on the density of the respective other particles in a neighborhood (sensing region) of the point x . Specifically, we use the so-called sensing radius, r s , so that in the limit N → ∞ the spatially averaged densities, p ̃ A , B ( x ) in Equation 1 are (Schimansky-Geier et al., 2021) p ̃ A , B x = 1 2 r s ∫ − r s r s d x ′ p A , B x + x ′ . (2) Averaging close to the boundaries is implemented in a reduced window over an inside range of the respective densities; see below for an explanation. The case of local interaction where p ̃ A , B ( x ) = p A , B ( x ) corresponds to r s = 0 .

Our model belongs to the broad class of active matter systems, ensembles of particles that continuously consume energy and therefore operate out of equilibrium (Romanczuk et al., 2012; Marchetti et al., 2013). In Equation 1, noise represents fluctuations in the internal drive of the particles. The time-dependent local density of the opposite species modulates each particle’s effective diffusivity (noise intensity). This reciprocal motility control represents energy influx at the microscopic scale, breaking detailed balance. Consequently, the observed steady states are maintained by a sustained flux of energy.

For the nonlinear function, f ( p ) , that determines the dependence of the noise intensity on the density of the other particle species, we choose a logistic function (Heinen et al., 2024), f p = 1 1 + exp − α p − p 0 , (3) a function that monotonically grows with its argument (hence, the fluctuations driving the A particles at x become stronger when there are a lot of B particles) but shows a saturation: the noise intensity cannot grow unbounded, as was the case for the quadratic function used in Schimansky-Geier et al. (2021). We set the maximum diffusion coefficient (or noise intensity) to one, noticing that a value different from one can be absorbed by simply rescaling time. The positive parameters α and p 0 in Equation 3 determine the slope and the inflection point of the sigmoid, respectively.

Implementing the average of the densities as estimated from simulations of a finite number of particles deserves some comments, see Schimansky-Geier et al. (2021) for details. In our numerical simulations, the sensing radius accounts for a moving average of binned probability density functions (PDF). For each bin and a given time instant, the PDF is estimated by the number of particles that have fallen into this bin. To estimate p ̃ A , B , these PDF values are averaged over s bins to the left and s bins to the right, i.e., over 2 s + 1 bins, except for bins close to the boundaries. With the bin size Δ x = 2 / ( m − 1 ) , this implies a sensing radius of r s = s ⋅ Δ x in Equation 2. The averaging window size is shrunk near the endpoints to include only existing bins.

In the limit of infinite particle numbers, N → ∞ , the dynamics of the system is described by coupled nonlinear and nonlocal Smoluchowski equations, which describe the evolution of the probability densities (Schimansky-Geier et al., 2021) ∂ t p A x , t = ∂ x 2 f p ̃ B x , t p A x , t , ∂ t p B x , t = ∂ x 2 f p ̃ A x , t p B x , t . (4) We emphasize that these nonlinear diffusion equations may permit more than one steady-state solution. We implement a numerical solution of these equations by a finite-difference scheme described in detail in Appendix B of Schimansky-Geier et al. (2021), testing different initial conditions to determine the number and stability of asymptotic solutions. The integration of the Smoluchowski equations was terminated once the probability distributions at two consecutive time steps differed by less than 1 0 − 8 at every spatial point x . The stationary value of the order parameter (time-independent) was then computed from the stationary probability densities.

As in Schimansky-Geier et al. (2021), the system admits a homogeneous steady-state solution p A , h ( x ) = p B , h ( x ) ≡ 1 / 2 (the interval length is 2). However, this homogeneous state may be unstable. To capture the emergence of inhomogeneity, we define the time-dependent order parameter I h t = ∫ − 1 1 p A x , t − 1 2 2 d x , (5) which becomes significantly larger than zero whenever the system is demixed. This quantity can be computed from both finite-population simulations and from numerical solutions of the coupled Smoluchowski equations.

As an additional measure of order in the system, we consider the instantaneous mean position of the A particles: ⟨ x A ⟩ = 1 N ∑ i = 1 N x A , i t , (6) which should remain close to zero when the A-particle distribution is homogeneous, but will deviate if A particles preferentially occupy one side of the interval (negative for the left side, positive for the right). Compared to I h , an advantage of this measure is that it explicitly distinguishes between the two possible demixed states. Finally, we note that these order parameters are symmetric with respect to an exchange of A and B particles, i.e., they remain invariant (with respect to their absolute value) under the substitutions p A → p B in Equation 5 and x A → x B in Equation 6.

For finite N , we employed the Euler–Maruyama scheme in the Itô interpretation to simulate the stochastic differential Equation 1 numerically. With the sensing radius, the dependence on the integration of the time step, Δ t , is weak: the larger r s (or s ), the weaker the dependence on Δ t (Schimansky-Geier et al., 2021). For s = 10 used in the following, the results for Δ t = 1 0 − 3 and 1 0 − 4 can hardly be distinguished. In contrast, omitting the averaging via sensing radius altogether and using directly the PDF itself in Equation 1, leads to an exceptionally strong dependence on the time step, such that extreme simulation times are required, especially to measure the jump statistics between metastable demixed and mixed states. In the following, we use the parameters as given in Table 1 if not explicitly stated otherwise.

The simulation time window was chosen to ensure a sufficient number of transitions between metastable states in the bistable parameter regime for reliable estimation of the mean residence times, the probability density, and the correlation function of the order parameter. For exponentially distributed residence times, the standard error of the sample mean τ ̂ of M switching events is SE ( τ ̂ ) = τ / M , where τ is the true mean residence time. By the central limit theorem, the relative margin of error of τ ̂ at confidence level 1 − α is ϵ = z 1 − α / 2 / M , with z 1 − α / 2 = Φ − 1 ( 1 − α / 2 ) the standard normal quantile, see e.g., Bendat and Piersol (2011). Thus, the required number of switching events is M ≥ ( z 1 − α / 2 / ϵ ) 2 . Accordingly, the expected simulation window is T max ≥ τ M . For ϵ = 0.1 and 1 − α = 0.95 this gives M ≥ 384 . We used M = 1 0 3 and τ = 1 0 4 , which results in T max = τ M = 1 0 7 , the value adopted in our simulations, accommodating for possible deviations from exponential statistics and occasional long residence times. To allow the system to approach a steady state, we first integrated for a relaxation period of T relax = 1 0 4 . The system was then integrated for an additional T max = 1 0 7 , and the time-dependent order parameter Equation 5 was approximated using the trapezoidal rule to compute the integral.

We start with the numerical analysis of Equation 4 and find the asymptotic state(s) of the system by solving the nonlinear equations for various initial conditions on a grid.

When we vary the parameters α and p 0 of the sigmoidal nonlinearity, we encounter parameter regimes that have been seen before with a quadratic nonlinearity in (Schimansky-Geier et al., 2021). Specifically, when the inflection point p 0 is low, and the function is not particularly steep (small α ), the homogeneous state in which both densities are well-mixed is the only stable state. When the inflection point is enlarged (but α is still small), we observe that the system demixes. There are regions where the A particles are in excess and the B particles are diminished in numbers, and there are other regions where it is the other way around; this is the demixed state. Now, beyond these two known regimes in which either a stable homogeneous or a stable demixed state exists, there is at large α and small p 0 also a regime of coexistence of both states. The lines of bifurcations are shown in Figure 1a; we also illustrate that the bistability does not hinge on a finite sensing radius (the red lines show the bifurcation for the case r s = 0 ). In Figure 1b, the new coexistence of demixed and homogeneous states is demonstrated for a specific parameter set: A and B particles might here be distributed in a homogeneous manner (black line) or by an increased density of A particles on the left and B particles on the right (because of the symmetry of the problem, there is the symmetric demixed solution with an excess of B particles on the left, etc.). We note that the densities, because of the finite averaging, do not add up to a constant, but there is an overall reduction of particles around x = 0 .

As a consequence of the different parameter regimes, we see different types of transitions when the inflection point p 0 is increased, depending on whether the logistic function is shallow ( α ≲ 7 ) or steep ( α ≳ 7 ) . For a shallow logistic function, the transition shown in Figure 1c is continuous: close to the bifurcation, the order parameter I h is still close to zero, and demixing starts as a slight difference in the densities from the homogeneous state. For a steeper nonlinearity as illustrated in Figure 1d, we see a first-order-like transition with a homogeneous state ( I h = 0 ) that remains stable when p 0 is increased (orange line), then jumps abruptly to a pronouncedly demixed state, and continues further from there. Upon reversal of the parameter variation (blue line), we see a hysteresis effect. Namely, the system stays in the bistable state where previously the abrupt jump occurred before it jumps to I h = 0 at a smaller value of p 0 than the one at which the jump in the forward direction was observed.

The observed bistability in the system with an infinite number of A and B particles (corresponding to the Smoluchowski equations) raises the interesting question of what will happen in a system with a finite number of particles that is prone to fluctuations. We explore the bifurcation in a system with many but finite particles next and then turn to smaller systems in which finite-size fluctuations may trigger transitions between bistable states.

We now turn to a large system with N = 1 0 4 of each A and B particle and integrate the stochastic differential equations Equation 1 with an α sufficiently large to encounter the above-discussed bistability. We use parameter continuation with different initial conditions: (i) the mixed state, in which x A and x B are uniformly distributed in ( − 1,1 ) , and (ii) the demixed state, where x A was distributed solely in ( − 1,0 ) , and x B in (0,1). From the two distributions of particles, we obtain time series I h ( t ) from which a time average I h ̄ as well as its probability density, P ( I h ) , can be estimated. We illustrate the mean value in the case of a parameter continuation p 0 for the two initial conditions (i) and (ii) in Figure 2a by solid and dashed lines and for the values of the sensing radius r s (red, blue, green). It becomes apparent that the exact transition point p 0 * as well as the width of the hysteresis depend on r s . With growing sensing radius, the critical inflection points move up, and the size of the hysteresis range shrinks. This is further explored in Figure 2b, where we show the time-averaged order parameter as a function of sensing radius for different values of p 0 . Also, here we observe an abrupt jump of I h ̄ with increasing r s ; it is quite plausible that too much averaging (a too large value of r s ) will destroy the demixing mechanism. Interestingly, the order parameter for moderate increases may also cause a modest growth in the order parameter (cf., for instance, the maximum of the green curve in Figure 2b).

How does the bifurcation plot depend on the number of used particles? This is exhibited in Figure 3 for our standard value of sensing radius, r s = 0.1 , and four different system sizes. The figure clearly shows hysteresis for large N ( 1 0 4 and 5 × 1 0 3 ), while for smaller N , the averaged order parameter displays a smoother dependence on p 0 . The latter dependence is expected if the system can jump between the bistable states within the averaging time window. We next explore the behavior of small particle numbers in more detail.

For small N , e.g., for N = 500 , the system exhibits switching between mixed and demixed states for parameter values within the range of bistability. Figure 4 exemplifies time traces of I h and ⟨ x A ⟩ (panel a) and corresponding probability distributions P ( I h ) (b) for three parameter values. For p 0 = 0.357 (red line), the switching results in a bimodal distribution P ( I h ) .

From Figure 4a we can also conclude that both I h and ⟨ x A ⟩ are viable measures of demixing. Furthermore, transitions between a demixed state with ⟨ x A ⟩ < 0 (A particles preferentially on the left) and a demixed state with ⟨ x A ⟩ > 0 (A particles preferentially on the right) or the other way around seem to have to pass through the demixed state first ( ⟨ x A ⟩ ≈ 0 ) ; this is definitely the case for smaller values of p 0 but becomes less well visible for p 0 = 0.36 for which the mixed state is least stable.

Figure 5 shows a complete picture of P ( I h ) vs. p 0 , demonstrating a region of switching dynamics. Superimposed are the locations of maxima and minima of P ( I h ) (illustrating the bimodality in the intermediate region) and the smoothly varying time-averaged order parameter, I ̄ h .

The stochastic transitions can be most appropriately characterized by the mean residence times in the mixed and demixed states. To this end, we extract many ( ≥ 1 0 3 ) realizations of residence times from the time course I h ( t ) . The residence times were estimated as time intervals between a threshold crossing events lasting for at least T min = 5 , to avoid fast threshold recrossing. The threshold was determined as the local minimum of the probability distribution of the order parameter (black circles in Figure 5). Figure Figure 6 shows the averages over the two types of intervals for N = 300 , 500, and 1000. The chosen range of p 0 values (with α = 8 ) falls entirely into the bistable regime shown in Figure 1a.

Clearly, a larger number of particles implies ‘more waiting’: for N = 1000 , the mean residence times are exponentially larger than for the smaller particle numbers. Furthermore, the mean residence time in the mixed (or homogeneous) state decreases upon increasing p 0 - this is quite plausible because we move the system away from the bifurcation boundary where this homogeneous state becomes the only stable state in the system and hence go the other way, we destabilize this state in the sense that we increase the system’s chance to leave this state. In contrast, the mean residence time of the demixed state increases - we stabilize this state by moving towards the bifurcation boundary at which the demixed state turns into the only stable state of the system.

There is a particular value of p 0 (which depends on N ) at which both mean residence times are equal. Here, both mixed and demixed states are roughly equally likely, and the histogram of I h may exhibit a pronounced bimodality. The particle number not only controls the finite-size noise in the system (in the sense of enlarging the mean residence times because fewer fluctuations can induce a transition) but also shifts the residence-time curves with respect to p 0 .

The slow back and forth between the different states of the system (especially at larger particle numbers and amid the bistable regime) may lead to long temporal correlations. We inspect their dependence on the system parameter by calculating the normalized correlation function of the order parameter R τ = ⟨ I h t I h t + τ ⟩ − ⟨ I h t ⟩ 2 ⟨ I h t 2 ⟩ − ⟨ I h t ⟩ 2 and extracting its typical decay time, the correlation time by the following estimate t cor = ∫ 0 T d τ | R τ | . (7) We note that the correlation time is easier to calculate than the residence times, as one does not need to define and use a threshold to determine the residence epochs. In Figure 7a, the correlation time Equation 7 displays a pronounced maximum at the sweet spot p 0 * , corresponding roughly to the equal residence of the mixed and demixed states. Making the system more asymmetric, i.e., biasing it towards the mixed or demixed state, will reduce the correlation time. Similar effects have been observed for the diffusion coefficient of particles with a bimodal velocity distribution (that is actually proportional to the correlation time of the velocity process) (Lindner and Nicola, 2008; Lindner and Sokolov, 2016) and for the count variability in bursting neurons (Kullmann et al., 2022). In many of these systems, it was observed that the system size controls the height of the maximum, and the maximal correlation time grows exponentially with N . This is also true for our system, cf. Figure 7b.

Furthermore, when the number of particles N is increased, the maximum becomes more peaked such that outside a critical range of asymmetry, the correlation time decreases with increasing N (or is almost constant as is the case for the two extreme examples shown in Figure 7b in green and blue) while inside this range t corr increases strongly with increasing N . The range itself has a moderate dependence on N , as has also been observed for some of the other systems mentioned (Lindner and Sokolov, 2016; Kullmann et al., 2022).

We have inspected a minimalistic model of two populations of Brownian particles that mutually enhance by their presence the diffusion coefficient of the respective other population. Building on the results from Schimansky-Geier et al. (2021), as a new model ingredients we used here a sigmoidal nonlinearity characterized by an inflection point p 0 and a slope parameter α , replacing a simple quadratic nonlinearity previously considered in Schimansky-Geier et al. (2021). We showed that, besides the already known regimes in which either only a homogeneous (well-mixed) state is stable or a demixed state is stable, this model exhibits a novel regime in which these two states (mixed and demixed) can coexist. This was first demonstrated using the nonlinear Smoluchowski equations corresponding to a system with infinite particles.

Exploring the case of finite particle numbers using Langevin equations revealed the emergence of repeated stochastic transitions between the homogeneous (disordered) mixed state and a demixed state in which one kind of particle is in excess in one region while the other one is in excess in another region. Hence, in a strong formulation, the system switches back and forth between order and chaos, which we quantified by an order parameter that measures the deviation from a uniform distribution. We note that another measure would be the entropy of the two distributions: we have verified that the entropy indeed displays stochastic transitions between two values corresponding to the more ordered, demixed state and the more disordered, homogeneously distributed state (not shown). Repeated back-and-forth transitions between order and disorder have been observed experimentally in several systems (Pufall et al., 2004; Paul et al., 2022); interestingly, in the context of our work, the simplicity of the model leads to such transitions between high- and low-entropy states.

We have furthermore explored the transitions using the residence times and the correlation time of the order parameter and, specifically, their dependence on the bifurcation parameters and the system size. Typical for a symmetric bistable setup, the correlation time is maximized and displays an exponential growth in the system size.

It would be fascinating to observe this type of transition between mixed and demixed states in experimental systems of active matter. The new feature of our nonlinear dependence of the diffusion coefficient on the probability density seems to be rather realistic: diffusion coefficients can certainly not grow unbounded with the number of opposite particles (at some point, the presence of other particles may even hinder and limit diffusion), hence a saturation of D ( p ) is quite plausible. Whether this would also be possible in more complicated models of mutual motility dependence and in the experiments presented in Curatolo et al. (2020), remains an exciting question for future research. In particular, the model can be extended by incorporating a quorum-sensing effect by introducing a self-dependent diffusivity, in which, in addition to cross-density dependence, the motility of a species is affected by its own density. Our preliminary results indicate that quorum-sensing self-dependent diffusivity can drive clustering, the so-called motility-induced phase separation of active Brownian particles (Cates and Tailleur, 2015). However, the bistable mixing/demixing reported here relies on reciprocal motility control between two species.

Our logistic mutual diffusion model captures self-organized demixing and noise-induced transitions between mixed and demixed states, behavior that resembles a thermodynamic first-order phase separation transition as in liquid-liquid phase separation (LLPS) (Banani et al., 2017; Zeigler et al., 2021). Further, cross-diffusion systems can be recast formally as Cahn–Hilliard–type equations (Berendsen et al., 2017) and therefore may share phenomenological features with LLPS models (Laghmach and Potoyan, 2020; Li and Hou, 2022). Conceptually, however, our model is an active media system and is not an LLPS model. In our model, the demixing is driven by reciprocal modulation of diffusion rather than by intermolecular interaction energies, as in LLPS.

The observed behavior in our model is a form of self-organization—more precisely, a transient self-organization that repeatedly emerges and decays. It may be useful to regard the spontaneous formation of high-concentration domains as a dynamical network whose nodes exchange particles and either adjust their size (approaching a single-interface metastable state) or disaggregate (transitioning to the disordered state). Methods from the theory of adaptive networks (Gross and Sayama, 2009; Berner et al., 2023) could be used to extend and analyze our framework in settings where state-dependent (context-dependent) interactions are expected. In our model, “state-dependent coupling” refers to the cross-diffusion gains depending on the locally averaged density of the partner species. Beyond microbial consortia, analogous activity-dependent couplings occur in several physiological contexts, for example, cancer–immune crosstalk in tissues (Hsieh et al., 2022; Marzban et al., 2024) and organ–organ interactions discussed in Network Physiology (Bashan et al., 2012; Qi et al., 2024; Lehnertz et al., 2020; Ivanov, 2021). We emphasize that these connections are qualitative, as quantitative mapping would require system-specific variables and timescales.

Our results may be viewed in the broader framework of synergetics, originally proposed by Hermann Haken, wherein emergent macroscopic patterns and states in nonlinear systems are cast via a small set of order parameters that “enslave” microscopic degrees of freedom (Haken, 1977; Haken 1983; Haken, 2012). In particular, our model’s demixing phenomenon, driven by reciprocal diffusivity feedback, exemplifies the hallmark of self-organization: above a threshold in the control parameter (here, the sigmoidal strength), a new collective state emerges and coexists with the homogeneous state, leading to noise-induced transitions between these metastable configurations.