Communication: Green-Kubo approach to the average swim speed in active Brownian systems Free

Assemblies of active, interacting Brownian particles (ABPs) are intrinsically nonequilibrium systems. In contrast to equilibrium, for which the statistical mechanics of Boltzmann and Gibbs enables the calculation of average properties, there is no analogous framework out-of-equilibrium. However, useful exact expressions exist, which enable average quantities to be calculated in the nonequilibrium system by integrating an appropriate time correlation function: the Green-Kubo formulae of linear response theory.^1–3 Transport coefficients, such as the diffusion coefficient or shear viscosity, are thus conveniently related to equilibrium autocorrelation functions. Given the utility of the approach, it is surprising that the application of Green-Kubo-type methods to active Brownian systems has received little attention.⁴ 

The primary aim of the present work is to extend Green-Kubo-type methods to treat ABPs. This approach has two appealing features. First, information about the active system can be obtained from equilibrium simulations. Second, the exact expressions derived provide a solid starting point for the development of approximation schemes and first-principles theory. The method we employ is a variation of the integration-through-transients approach, originally developed for treating interacting Brownian particles subject to external flow.^5–8 

A fundamental feature of ABPs is the persistent character of the particle trajectories. For strongly interacting many-particle systems, the interplay between persistent motion and interparticle interactions can generate a rich variety of collective phenomena, such as motility-induced phase separation (see Ref. 9 for a recent overview). A quantity which features prominently in many theories of ABPs^9–13 is the density-dependent average swim speed, which describes how the motion of each particle is obstructed by its neighbours. Given the ubiquity of the average swim speed in the literature on ABPs, we choose it as a relevant observable with which to illustrate our general Green-Kubo-type approach. We demonstrate that this quantity can be obtained from a history integral over the equilibrium autocorrelation of tagged-particle force fluctuations, which we investigate in detail using Brownian dynamics (BD) simulation.

We consider a three dimensional system of N active, interacting, spherical Brownian particles with coordinate r_i and orientation specified by an embedded unit vector p_i. A time-dependent self-propulsion of speed v₀(t) acts in the direction of orientation. Allowing for time-dependence of this quantity both clarifies the general structure of the theory and leaves open the possibility to model physical systems for which the amount of fuel available to the particles is not constant (see, e.g., Refs. 14 and 15). Omitting hydrodynamic interactions, the motion can be modelled by the Langevin equations

where γ is the friction coefficient and the force on particle i is generated from the total potential energy according to F_i = − ∇_iU_N. The stochastic vectors ξ_i(t) and η_i(t) are Gaussian distributed with zero mean and have time correlations 〈ξ_i(t)ξ_j(t′)〉 = 2D_t1δ_ijδ(t − t′) and 〈η_i(t)η_j(t′)〉 = 2D_r1δ_ijδ(t − t′). The translational and rotational diffusion coefficients, D_t and D_r, are treated in this work as independent model parameters.

It follows exactly from (1) that the joint probability distribution, P(r^N, p^N, t), evolves according to¹⁶ 

The time-evolution operator can be split into a sum of two terms, Ω_a(t) = Ω_eq + δΩ_a(t), where the equilibrium contribution is given by

with rotation operator R = p × ∇_p (see, e.g., Ref. 17) and the time-dependent, active part of the dynamics is described by the operator

To solve (2) we define a nonequilibrium part of the distribution function, δP(t) = P(t) − P_eq,⁷ where P_eq is the equilibrium distribution of position and orientation. Using Ω_eqP_eq = 0 yields the equation of motion

Treating the last term as an inhomogeneity and solving for δP(t), we obtain a formal solution for the nonequilibrium distribution

where e₊(⋅) is a positively ordered exponential function (see the appendix in Ref. 8) and we have used δΩ_a(t) P_eq = − βv₀(t) F^pP_eq, with “projected force” fluctuation

The projected force emerges as a central quantity within our approach and indicates to what extent the interparticle interaction forces act in the direction of orientation, either assisting or hindering the self-propulsion. We will show that this quantity is closely related to the average swim speed in the active system.

Introducing a test function, f, on the space of positions and orientations and integrating (6) by parts yield a formally exact expression for a nonequilibrium average

where e₋(⋅) denotes a negatively ordered exponential⁸ and 〈 ⋅ 〉_eq is an equilibrium average over positional and orientational degrees of freedom. The adjoint operator is given by $Ω a † ( t ) = Ω eq † −δ Ω a ( t )$⁠, where

generates the equilibrium dynamics. The integrand appearing in (8) involves the equilibrium correlation between the projected force at time t′ and the observable f, which evolves from t′ to t according to the full dynamics. The average is nonlinear in v₀(t) because of the activity dependence of the adjoint operator.

The response of the system to linear order in v₀(t) is obtained by replacing the full time-evolution operator $Ω a † ( t )$ in (8) by the time-independent equilibrium operator $Ω eq †$⁠. Further simplification occurs if the activity is constant in time, v₀(t) → v₀, leading to

which can be used to define a general active transport coefficient $α= lim v 0 → 0 ( 〈 f 〉 lin − 〈 f 〉 eq ) / v 0$⁠. Equation (10) is the desired Green-Kubo relation for calculating the linear response of ABPs to a time-independent activity.

As mentioned previously, a quantity of current interest is the average, density-dependent swim speed, v(ρ). This describes how the bare swim speed, v₀, is influenced by interparticle interactions and is an important quantity in many of the various theories addressing ABPs.^9–13 In particular, the tendency of the system to undergo motility-induced phase-separation is determined by the rate of decrease of v(ρ) with increasing density; a positive feedback mechanism can result when increasing the local density leads to a sufficiently strong reduction of the local average swim velocity. The average swim speed is defined as the nonequilibrium average

where v_i is the velocity of particle i. Using (1) to eliminate the velocity in favour of the forces and using the fact that the Brownian force ξ_i is uncorrelated with the orientation p_i, it follows that

For a time-independent v₀, we can employ (10) to calculate the average in (12) to linear order

where the integrand is the equilibrium autocorrelation of projected force fluctuations

Spatial and orientational degrees of freedom decouple in equilibrium, which enables the orientational integrals in (14) to be evaluated exactly. This yields

where F is the interaction force acting on an arbitrarily chosen (“tagged”) particle, $Ω eq,s † = ∑ i D t ∇ i + β F i ⋅ ∇ i$ is the spatial part of the time-evolution operator and 〈 ⋅ 〉_eq,s indicates an equilibrium average over spatial degrees of freedom. The initial value is given by H(0) = β²〈|F|²〉_eq/3. If we consider pairwise additive interaction potentials, then the Yvon theorem²⁰ leads to

where ρ is the number density, u(r) is the passive pair potential, and g_eq(r) is the corresponding equilibrium radial distribution function.

Equation (15) shows that the nontrivial physics underlying the linear response of the system to activity is contained in the tagged-particle force-autocorrelation function. This function was encountered many years ago by Klein and co-workers¹⁸ in a study of the velocity autocorrelation in overdamped Brownian systems. By manipulation of the operator (3), it was shown that

where Z_eq(t) is the velocity autocorrelation function, defined in terms of the tagged particle velocity, v(t), according to the familiar relation

The velocity autocorrelation function is a quantity of fundamental interest in describing the dynamics of interacting liquids and is closely related to other important quantities (e.g., the mean-squared displacement and self-diffusion coefficient). Substituting (17) into (15) yields

thus providing, via (13), a direct connection between Z_eq(t) and v(ρ). The latter can thus be determined to linear order in v₀ using a standard, equilibrium BD simulation. The δ function appearing in Eq. (19) is cancelled by an equal contribution from Z_eq(t)¹⁸ ensuring that H(t) remains finite at t = 0. Finally, we note that H(t) remains integrable in all spatial dimensions because of the exponential in (15). There is thus no principal difficulty in calculating v(ρ) in two dimensions, in contrast to the situation for transport coefficients, such as the self-diffusion coefficient, for which the relevant Green-Kubo time-integral diverges.⁵ 

In a recent study of the pressure in active systems, Solon et al.¹⁹ express the density-dependent average swim speed in the form v(ρ) = v₀ + I₂/ρ, where ρ is the bulk number density. The interaction potential is encoded in the quantity I₂ via its dependence on a static structural correlation between density and polarization, which are given, respectively, by the first and second harmonic moments of the orientation-resolved single particle density. This leads to the identification of $I 2 =− D t ρ v 0 ∫ 0 ∞ dt H ( t )$⁠. An advantage of the present Green-Kubo formulation over that of Solon et al. is that it enables identification of the relevant relaxation processes contributing to the decrease of v(ρ). Moreover, we anticipate that (13) will prove more convenient for the development of approximations.

In order to test the range of validity of the linear response result (13), we perform BD simulations on a three-dimensional system of N = 1000 particles interacting via the pair-potential βu(r) = 4ε((σ/r)¹² − (σ/r)⁶), where σ sets the length scale and we set ε = 1. The potential is truncated at its minimum r = 2^1/6σ to yield a softly repulsive interaction. The system size L is determined as L = (N/ρ)^1/3 in order to obtain the desired density. The integration time step is fixed to dt = 10⁻⁵τ_B where τ_B = d²/D_t is the time scale of translational diffusion. The equation for time evolution of orientation vector (Eq. (1)) is evaluated as an Ito integral. Measurements are made after a minimum time of 20τ_B to ensure equilibration. In order to measure time-correlations, the system is sampled every τ_p/100 s, where τ_p = 1/2D_R is the rotational diffusion time scale. The total run time is 300τ_B. We choose the ratio of diffusion coefficients as D_r/D_t = 20, although there is nothing special about this particular choice.

In Fig. 1(a) we show the correlator H(t) as a function of time for a number of different densities, the largest of which is close to the freezing transition for our model interaction potential. Aside from the strong increase of H(0) with increasing density (shown in the inset), the most striking aspect of the correlator is that the decay of H(t) is much faster than the time scale of rotational diffusion (note that time is scaled with τ_p in the figure). Indeed, very large values of the ratio D_r/D_t would be required for the exponential factor in (15) to significantly influence the decay of H(t). In the limit of large D_r, we obtain H(t) = H(0)exp(−2D_rt) and thus v(ρ)/v₀ = 1 − H(0) D_t/(2D_r). We conclude that, provided the value of D_r is not extremely large, the relevant relaxation process is the decorrelation of the tagged particle interaction force.

In the inset of Fig. 1(a), we show the initial value, H(0), as a function of the density. To check expression (16), we have confirmed that using g_eq(r) from our equilibrium simulations to evaluate the r.h.s. indeed reproduces the t → 0 limit of our dynamical H(t) data. Moreover, we have also employed an approximate liquid-state integral equation theory (Percus-Yevick theory)²⁰ to calculate g_eq(r) and evaluate H(0). Very good agreement of the predicted H(0) with simulation data is obtained.

In Fig. 1(b) we replot the data on a semi-logarithmic scale, with the initial value scaled out. This representation makes clear that H(t) is non-exponential and that the decay occurs more rapidly as the density is increased, in contrast to the structural relaxation of the system, which slows down with increasing density. The latter observation can be rationalized by considering that small positional changes can give rise to large changes in the force for closely packed particles residing in regions of strong interaction-force gradient. The fact that H(t) is non-exponential is not surprising, given that it can be expressed in terms of the velocity autocorrelation function, a quantity which famously exhibits power law asymptotic behaviour (“long-time tails”).^18,20 Klein et al. have shown analytically that for a dilute system of Brownian hard-spheres $〈F ( t ) ⋅F ( 0 ) 〉 eq,s ∼ t − 5 2$ for long times.

In Fig. 2 we show simulation data for the average swim speed as a function of density. The red diamonds show the linear response prediction obtained by using the data of Fig. 1 in the integral expression (13). This yields a result for v(ρ)/v₀ which is independent of v₀. The remaining curves show data obtained by direct evaluation of (11) using active BD simulations at three different values of v₀. As one might expect, deviations from linear response occur at lower density for larger values of v₀.

The above observation can be made more concrete by estimating a region in the (v₀, ρ) plane where linear response breaks down. In Fig. 3 we use our simulation data to map the locus of points for which the error in the linear response result, relative to the full active BD simulations, equals 5%. Although the chosen criterion is somewhat arbitrary, it at least gives a visual impression of the range of validity of linear response within the space of our control parameters. The locus of points shown in Fig. 3 is correlated with the onset of strong spatial inhomogeneities and phase separation. However, an analysis of active phase separation would go beyond the scope of the present work. The linear response formula (13) thus appears to be reliable for parameter values away from phase separation, but, beyond this, higher orders in v₀ will become important in determining v(ρ).

To summarize our main findings, we have derived a formally exact expression (8) for calculating averages in a system of interacting Brownian particles, subject to a time-dependent activity v₀(t). From this we obtain the linear-response expression (10) for a time-independent activity. Application of this result to calculate the average swim speed yields (13) and identifies the relevant time-correlation function, H(t), as given by (15). We find that linear response provides an accurate account of v(ρ) over a large parameter range, except for those regions of parameter space where phase separation occurs.

Although we have focused our attention on the linear-response regime, our exact results could in principle be used to develop nonlinear theories in the spirit of Refs. 6–8, which address Brownian particles under external flow. An important aspect that will be pursued in a future study is the identification of the range of validity of the linear-response regime, in particular, how it depends on v₀ and D_r. It would also be interesting to use (8) to investigate the transient dynamics arising from time-dependent activity, but we defer this line of enquiry until an experimentally relevant protocol can be identified. Aside from using an equilibrium integral equation theory to determine H(0) (inset to Fig. 1(a)), all of the data presented come from BD simulation. A clear next step is to investigate approximations to H(t) which enable predictions to be made from first-principles, without simulation input. Given relation (19), it seems likely that existing approximations to the velocity autocorrelation function (e.g., projection operator approaches) could be usefully exploited.