Multidimensional stationary probability distribution for interacting active particles.

We derive the stationary probability distribution for a non-equilibrium system composed by an arbitrary number of degrees of freedom that are subject to Gaussian colored noise and a conservative potential. This is based on a multidimensional version of the Unified Colored Noise Approximation. By comparing theory with numerical simulations we demonstrate that the theoretical probability density quantitatively describes the accumulation of active particles around repulsive obstacles. In particular, for two particles with repulsive interactions, the probability of close contact decreases when one of the two particle is pinned. Moreover, in the case of isotropic confining potentials, the radial density profile shows a non trivial scaling with radius. Finally we show that the theory well approximates the "pressure" generated by the active particles allowing to derive an equation of state for a system of non-interacting colored noise-driven particles.

A generic system, that is in thermal equilibrium at a temperature T, will be found in the neighbourhood of a configuration of energy E with a probability density given by the the Boltzmann factor exp[− E/kT]12. As an example, two Brownian colloidal particles, interacting through a conservative attractive force, will show an increased probability density for the low energy bound state. A quite different behaviour is observed in active particle systems3. Generally speaking, active matter is composed by biological or synthetic objects that are capable of absorbing energy from the environment and convert it into different kinds of persistent motions. Even when stationary states are reached, the probability distributions can display large deviations from their equilibrium Boltzmann counterparts. Those deviations are not just a matter of quantity, but a radically different qualitative behaviour may be observed, like the widespread tendency to accumulate around repulsive objects. This “attraction for repulsion” is responsible for phenomena like, particle accumulation at solid walls4 or the formation of bound states between repulsive objects56. Our intuitive notion that particles like to stay where external forces attract them to is biased by our familiarity with equilibrium statistical mechanics and needs to be replaced by novel statistical mechanics concepts that are capable to describe the stationary probability distributions in systems of interacting active particles. In this context some schematic models have been proposed to model the dynamics of active particles as, for example, the “run and tumble” (RnT) model. The RnT dynamics is appropriate to describe the motion of bacteria such as E. coli78910 that swim along almost straight runs interrupted by random reorientations. In the case of colloids propelled by chemical reactions the “active Brownian” model describes the motion of particles pushed by a force of constant magnitude that gradually reorient by rotational Brownian motion11121314. However, despite the simplicity of the dynamics of these systems, it is hardly possible to find the analogue of the Boltzmann distribution. Indeed, the Boltzmann prescription assigns a precise weight to a given configuration of positions and momenta of particles at equilibrium. These particles are embedded in a space of dimensionality d, are subject to arbitrary external fields and mutually interact via whatsoever potential12. On the contrary, in the case of active particles the exact stationary probability distribution is known only in rare instances as, for example, in the 1-dimensional RnT model in an external force field9. The impossibility of writing explicitly the stationary probability density prevents one from applying the standard methods of statistical mechanics. A Gaussian colored-noise model (GCN) can be used to reproduce the dynamics of passive colloidal particles immersed in a bath of dense swimming bacteria1516. This model has been intensively studied in the past as the simplest model that could elucidate the basic physics of systems subject to time-correlated noise. Interestingly GCN was originally used to interpret the behaviour of very different physical systems such as noisy electronic circuits17 and dye-laser radiation18. The analytical study of GCN-driven systems resulted very challenging and led to the development of different approximation schemes aimed to reduce the complexity of the GCN-model to a tractable level1920. Among these approaches one emerges by having a number of advantages with respect to the others. This is the Unified Colored Noise Approximation (UCNA) developed by Hänggi and Jung20 that, under certain conditions, describes both the small and the large correlation-time regimes, both in the high and low-friction limit21. More importantly for the present work the UCNA scheme can be generalized to a phase space of arbitrary dimensionality22. In this work we report, for the first time, the explicit formula of the stationary probability (obtained within the UCNA) for a system that is subject to a generic conservative potential and that is composed by an arbitrary number of degrees of freedom. We name this “multidimensional unified colored noise approximated stationary probability” (MUCNASP). The MUCNASP plays basically the same role as the Boltzmann distribution for the approximated GCN-driven system. We show how the MUCNASP allows to predict several non-equilibrium properties of the active system in experimentally relevant cases where a simple external potential acts on a small number of degrees of freedom. We focus on the case of steep repulsive interactions and spherically symmetric external potentials. In all these situations we use numerical simulation to test the quality of the approximation and find that the GCN-driven and the RnT particles display a strikingly similar behavior. In particular we show how our approximated probability density captures very well the accumulation of the active particles around repulsive obstacles. Moreover the theory describes well the dependence on dimensionality of the probability density function when the active particles are confined by a circular repulsive wall. Understanding how the concept of pressure generalizes to active matter has become recently the subject of intense theoretical research 232425. In this context we show how our theoretical probability density allows us to derive the pressure that the active particles exert on the repulsive walls and leads to the derivation of equations of state for the non-interacting active particle system. Finally we discuss the most relevant limitations of the present theory and suggest new routes to follow in the theoretical study of active matter.

Results

We consider the following set of stochastic differential equations:where the position variables  = (x1,…,x) are determined by the deterministic velocity − ∇Φ generated by a conservative potential Φ(x) and by a set of stochastic processes  = (η1,…,η). We assume that these are N independent Gaussian processes with zero mean 〈η〉 = 0 and exponential time-correlation: 〈η(t)η(s)〉 = . Here D is the diffusion coefficient characterizing the amplitude of the noise and τ is its relaxation time. Note that here we absorb the mobility μ in the velocity field − ∇Φ = μ, where  =  − ∇U is the deterministic force generated by the potential energy function U(x) = Φ(x)/μ. By using the UCNA Eq. (1) reduces to the (Stratonovich) Langevin equation22:where is the identity matrix, is the Hessian associated with the potential Φ, and Γ is a set of independent white-noise sources having 〈Γ〉 = 0 and 〈Γ(t)Γ(s)〉 = 2δδ(t − s). We have found that, in the flow-free case, the steady state probability of finding the dynamical system of Eq. (2) in a specific configuration is always proportional to the weight:where |…| represents the norm of a vector and ||…|| indicates the absolute value of the determinant of a matrix. We have demonstrated the validity of Eq. (3) by deriving the corresponding Fokker-Planck equation from Eq. (2)26 and solving it in the zero-current case by using the Jacobi’s formula27 (see Supplemental Material).

One single degree of freedom

As Eq. (3) is used for specific choices of the potential it reveals several interesting non-equilibrium properties of the active system under study. We initially focus on a simple one-dimensional case and study GCN-driven particles when they are subject to a steep repulsive potential of the form Φ(x) = Ax−12 setting A = 1. Such a potential can be thought as a repulsive obstacle that perturbs the dynamics of the particles28. To verify the quality of the UCNA, we integrate numerically the stochastic equation of motion (1) in the presence of such a potential. To this aim we have implemented a code for Euler integration of Eq. (1), to be executed on GPU where the dynamics of many independent particles can be simulated in parallel29. We consider several different values of 0.1 ≤ τ ≤ 1 s and 0.1 ≤ D ≤ 100 μm2/s, ranges that cover the typical persistence times and diffusivities of colloids in bacterial baths, swimming bacteria such as E. coli3031 and of chemically self-propelled Janus-type particles13. The size of the simulation box is chosen to be L = 20μm (with periodic boundaries located at ± L/2). Fig. 1(a) shows that in equilibrium (τ = 0) the probability density decreases rapidly before the core of the repulsive potential is reached, whereas in the GCN the distribution peaks substantially in a region where the potential is very high before vanishing at the core. Note that the specific choice of the constant A = 1 defines the size of the repulsive “wall” created by the external potential. The thickness of this impenetrable region is about 2μm and it depends very little on the values of D and τ considered since the potential is steeply repulsive. In the GCN-driven system the exact probability distribution is unknown but it can be approximated by Eq. (3) that reduces, for a single degree of freedom, to the known form19:where the prime indicates differentiation with respect to x. The approximate probability, obtained by normalizing Eq. (4) , is plotted in Fig. 1(a) as a dashed line and it is found to reproduce well the numerical distribution at two well separated values of D and τ. Knowing the probability we can also compute all the average quantities of interest, such as the average value of the modulus of the velocity . The theoretical (approximated) 〈|Φ'|〉 is compared with the numerical value in Fig. 1(b) where one sees that the 〈|Φ'|〉 UCNA prediction is very close to the numerically result at all values of D and τ here investigated. The behaviour of the P(x) and of 〈|Φ'|〉 can be qualitatively understood by considering the strong peaking of the Ω(x) close to the repulsive barrier. When the potential is very steep, as in the case of Φ = x−12, the maxima of Ω(x) are found at x = x* where . It is clear that the probability peaks where the external potential balances the root mean-squared velocity of the particle induced by GCN . In this hard-wall limit is possible to approximate Ω(x) in the neighborhood of x* by a strongly peaked function k(x) whose integral is , as found by a saddle-point approximation, while far away from x* Eq. (4) reduces to unity (see Supplemental Material) and we can write: Ω(x) ≈ k(x − x*) + k(x + x*) + 1. Note that corresponds to the typical correlation length of the active motion. Integrating from − L/2 to L/2 we find the average velocity , where L* = L − 2x* is the overall length available to the particles and report this result as a dashed line in Fig. 1(b) where it captures the trend of 〈|Φ'|〉 obtained in simulations with the potential Φ = x−12. By considering the force exerted by the particles located only on the right-hand side of the potential (x > 0) we find the force f: , which corresponds to an equation of state being f the 1-dimensional pressure that the particles exert on the “wall” represented by the external potential. Note that if we consider N independent particles, in the limit τ → 0, and set D = μkT we arrive to the ideal gas law in 1d: 〈f〉 = NkT/L* and that this equilibrium value constitute an upper bound for the pressure of the active system (see dashed-dotted line in Fig. 1(b). Interestingly these results can be derived exactly for the RnT model, in particular the stationary probability distribution of the RnT model in presence of two hard walls is composed by two Dirac deltas plus a constant and the expression of 〈|Φ'|〉 has the same form of the one found in the UCNA (see Supplemental Material).

Figure 1

(a) Probability density function of the position of a single GCN-driven particle in presence of the external potential Φ = x−12 (shaded area). Full lines: simulations, dashed lines: theory, dashed-dotted line: Boltzmann distribution. The curve with the higher peak corresponds to D = 100μm2 /s, τ = 1 s and the one with the lower peak to D = 0.4μm2 /s, τ = 0.1 s, (zoomed in the inset) (b) Average value of |Φ'| as a function of D for three different values of τ = 0.1, 0.325, 1 s from top to bottom respectively. Points: simulations, full lines: theory, dashed lines: theory in the limit of a hard potential, dashed dotted line: white noise case with a hard potential. (c) Probability density function of the distance between two GCN-driven particles interacting via the potential Φ = Δx−12, same legend as Fig. (a). (d) Average value of |Φ'| as a function of D for three different values of τ = 0.1,0.325,1 s from top to bottom respectively, same legend as as Fig. (b).

Two interacting particles in one dimension

Up to this point we have derived results from the known formula (4) considering only one degree of freedom. We now use Eq. (2) to characterize the steady state properties of two interacting GCN-driven particles moving in 1d. We consider the positions x1 and x2 of two particles interacting via a pair potential that is a function of the distance Δx = |x1 − x2| between the particles: Φ(x1,x2) = Φ(Δx). We further assume that the particles are free to move in a 1d space of extension L with periodic boundaries located at ± L/2. In this case Eq. (2) can be used to compute the probability of finding the two particles separated by Δx:which is identical to Eq. (4) with τ replaced by 2τ. This is at variance with equilibrium statistical mechanics in which the probability of finding two particles at a given distance does not vary if one of the two particles is pinned at some fixed position323334. To be more specific let us consider a repulsive potential of the form Φ(x1, x2) = (x1 − x2)−12. In this case again the probability is well described by the MUCNASP as shown in Fig. 1(b). Similarly the deterministic velocity component experienced by one particle 〈|Φ'|〉 resulting from simulations is well approximated by the theory (see Fig. 1(c). Again, in the limit of an infinitely steep potential (see Supplemental Material), we find that the area of the peak of Ω(Δx) is approximated by and which is plotted in Fig. 1(c) as a dashed line. This can be physically interpreted as follows: when both particles are free to move they can be more often found in contact since they move coherently in the same direction, this happens without the particles pushing onto each other, yielding a lower value of the average interaction force. It is important to note that these theoretical results cannot be derived analytically in the RnT model since the coupled dynamics of more particles makes the problem far too complicated. Nevertheless we have found that the MUCNASP produces results for the average 〈|Φ'|〉 that are very similar those found numerically for the RnT model despite the stationary probability density has a very different form (see Supplemental Material). This suggests that the MUCNASP can be used as a convenient approximation also for calculating the averages in RnT dynamics at least in the case of steeply repulsive potentials. Note also that such a scenario is in agreement with the findings of Ref.5 where it was demonstrated, by combining experiments and simulations, that two colloids suspended in a bacterial bath tend to stay in contact because of the colored-noise forces induced by swimming bacteria.

Radially symmetric potentials

When the d-dimensional potential is spherically symmetric, i.e. Φ(x) = Φ(r), Eq. (3) simplifies towhere and Θ is the d-dimensional solid angle. Note that the Boltzmann distribution, obtained by setting τ = 0 in Eq. (6), depends on the dimensionality only via the trivial term r while in the GCN-case this dependence is more complicated. To understand this issue, we consider GCN-driven particles in d = 2 in the presence of a circular repulsive potential of radius R of the form Φ(r) = (r − R)−12 where and R = 5μm. Simulation results show that particles accumulate near the ring at r = R and the theoretical probability reproduces well this behaviour (see Fig. 2(a). From Eq. (6) we can compute the averages of interest as the radial component of the velocity field 〈|Φ'|〉 and compare it with simulation results in Fig. 2(d) showing a good agreement. In the limit of an infinitely steep potential we have that Ω(r) strongly peaks where and reduces to unity elsewhere. The area of the peak can be approximated by (see Supplemental Material), where R* = R − r* is the radial coordinate of the peak with r* ≈ 1 in the D-τ range considered. For the average radial velocity component we get which is plotted as dashed lines in Fig. 2(d) and follows nicely the trend displayed by the numerical data. This is practically a colored-noise version of the ideal gas law in a circular container. This is clear if we set τ = 0 to obtain 〈|Φ'|〉 = 2D/R* which is proportional to the average radial force 〈f〉 = 2D/(R*μ). Dividing this by the 2d “surface” 2πR* we get the ideal gas pressure p = N〈f〉/(2πR*) = NkT/(πR*2) for N independent particles and D = μkT. The unperturbed dynamics of the RnT model in 2d is well understood35, while the problem of a 2d symmetric potential is not tractable analytically. However we have found that the simulation results for the 〈|Φ'|〉 in RnT model are very similar to those produced by the MUCNASP (see Supplemental Material) suggesting that MUCNASP could actually describe the behaviour of a wider class of active particles. An experimental situation, that is close to the spherically symmetric case treated here, has been recently studied36. In these experiments it was observed a marked accumulation of swimming bacteria at the border of spherical liquid droplets. In this kind of experiment it would be interesting to check whether, in the dilute regime, the number of bacteria found in contact with the border scales with the droplet radius and the characteristic run length as predicted by the MUCNASP (Eq. (6)) in the spherical case.

Figure 2

(a) Probability density function of the radial distance of one GCN-driven particle in presence of the spherically symmetric external potential Φ = (r − R)−12 in 2d (shaded area). Full lines: simulations, dashed lines: theory, dashed-dotted line: Boltzmann distribution. The curve with the higher peak corresponds to D = 100μm2 /s, τ = 1 s, and the one with the lower peak to D = 0.4μm2 /s, τ = 0.1 s, (zoomed in the inset) (b) Average value of |Φ'| as a function of D for three different values of τ = 0.1, 0.325, 1 s from top to bottom respectively. Points: simulations, full lines: theory, dashed lines: theory in the limit of a hard potential, dashed dotted line: white noise case with a hard potential.

Discussion

By using the unified colored noise approximation, we have derived the explicit formula for the non-equilibrium stationary probability (MUCNASP, Eq. (3)) of a system composed by an arbitrary number of degrees of freedom subject to GCN. We have focused onto the case of steep repulsive potentials where the probability distribution of one single active particle tends to concentrate on the repulsive part of the potential oppositely to the case of a Brownian particle. Moreover we have verified that, as predicted by the MUCNASP, two active particles interacting repulsively behave differently with respect to equilibrium and are found in contact more often than if one particle is fixed. Finally the MUCNASP predicts that, when active particles are confined inside a repulsive ring-shaped barrier, the probability peaks on the boundary and the area of this peak increases with increasing radius and with increasing persistence length. Surprisingly the results obtained by the MUCNASP are very close to those obtained for RnT particles for which an analytical solution is not available. As discussed before192022, the UCNA is accurate in those portions of phase space where the all the eigenvalues of the Hessian matrix are positive. This restriction defines where Eq. (3) can be used as a valid approximation for the probability of GCN-driven system. Moreover it appears difficult to derive an explicitly formula for probability including also Brownian fluctuations. Nevertheless, to our knowledge, the MUCNASP is the only available explicit probability formula accounting for multiple active degrees of freedom and, provided that all eigenvalues are positive, becomes exact both in the limit of τ → 0 and τ → ∞. This makes the MUCNASP a valuable schematic model for tackling the many-body problem in active matter. For example it would be interesting to check whether, at the mean-field level, the MUCNASP can predict a motility-induced phase separation103738 or the colored-noise induced shift in the synchronization threshold of the noisy Kuramoto model39.

Additional Information

How to cite this article: Maggi, C. et al. Multidimensional Stationary Probability Distribution for Interacting Active Particles. Sci. Rep. 5, 10742; doi: 10.1038/srep10742 (2015).