Green-Kubo approach to the average swim speed in active Brownian systems
GGreen-Kubo approach to the average swim speed in active Brownian systems
A. Sharma and J.M. Brader
Department of Physics, University of Fribourg, CH-1700 Fribourg, Switzerland
We develop an exact Green-Kubo formula relating nonequilibrium averages in systems of interact-ing active Brownian particles to equilibrium time-correlation functions. The method is applied tocalculate the density-dependent average swim speed, which is a key quantity entering coarse grainedtheories of active matter. The average swim speed is determined by integrating the equilibrium au-tocorrelation function of the interaction force acting on a tagged particle. Analytical results arevalidated using Brownian dynamics simulations.
PACS numbers: 82.70.Dd,64.75.Xc,05.40.-a
Assemblies of active, interacting Brownian particles(ABPs) are intrinsically nonequilibrium systems. In con-trast to equilibrium, for which the statistical mechanicsof Boltzmann and Gibbs enables the calculation of av-erage properties, there is no analogous framework out-of-equilibrium. However, useful exact expressions ex-ist, which enable average quantities to be calculated inthe nonequilibrium system by integrating an appropri-ate time correlation function; the Green-Kubo formulaeof linear response theory [1–3]. Transport coefficients,such as the diffusion coefficient or shear viscosity, arethus conveniently related to equilibrium autocorrelationfunctions. Given the utility of the approach it is surpris-ing that the application of Green-Kubo-type methods toactive Brownian systems has received little attention [4].The primary aim of the present work is to extendGreen-Kubo-type methods to treat ABPs. This approachhas two appealing features. Firstly, information aboutthe active system can be obtained from equilibrium simu-lations. Secondly, the exact expressions derived provide asolid starting point for the development of approximationschemes and first-principles theory. The method we em-ploy is a variation of the integration-through-transientsapproach, originally developed for treating interactingBrownian particles subject to external flow [5–8].A fundamental feature of ABPs is the persistent char-acter of the particle trajectories. For strongly interact-ing many-particle systems the interplay between persis-tent motion and interparticle interactions can generatea rich variety of collective phenomena, such as motility-induced phase separation (see [9] for a recent overview).A quantity which features prominantly in many theoriesof ABPs [9–13] is the density-dependent average swimspeed, which describes how the motion of each particleis obstructed by its neighbours. Given the ubiquity ofthe average swim speed in the literature on ABPs wechoose it as a relevant observable with which to illustrateour general Green-Kubo-type approach. We demonstratethat this quantity can be obtained from a history integralover the equilibrium autocorrelation of tagged-particleforce fluctuations, which we investigate in detail usingBrownian 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 vec-tor 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 generalstructure of the theory and leaves open the possibility tomodel physical systems for which the the amount of fuelavailable to the particles is not constant (see e.g. [14, 15]).Omitting hydrodynamic interactions the motion can bemodelled by the Langevin equations˙ r i = v ( t ) p i + γ − F i + ξ i , ˙ p i = η i × p i , (1)where γ is the friction coefficient and the force on par-ticle i is generated from the total potential energy ac-cording to F i = −∇ i U N . The stochastic vectors ξ i ( t )and η i ( t ) are Gaussian distributed with zero mean andhave time correlations (cid:104) ξ i ( t ) ξ j ( t (cid:48) ) (cid:105) = 2 D t δ ij δ ( t − t (cid:48) ) and (cid:104) η i ( t ) η j ( t (cid:48) ) (cid:105) = 2 D r δ ij δ ( t − t (cid:48) ). The translational androtational diffusion coefficients, D t and D r , are treatedin this work as independent model parameters.It follows exactly from (1) that the joint probabilitydistribution, P ( r N , p N , t ), evolves according to [16] ∂P ( t ) ∂t = Ω a ( t ) P ( t ) . (2)The time-evolution operator can be split into a sum oftwo terms, Ω a ( t ) = Ω eq + δ Ω a ( t ), where the equilibriumcontribution is given byΩ eq = N (cid:88) i =1 ∇ i · (cid:2) D t ( ∇ i − β F i ) (cid:3) + D r R i , (3)with rotation operator R = p ×∇ p (see, e.g. [17]) and thetime-dependent, active part of the dynamics is describedby the operator δ Ω a ( t ) = − N (cid:88) i =1 v ( t ) ∇ i · p i . (4)To solve (2) we define a nonequilibrium part of the dis-tribution function, δP ( t ) = P ( t ) − P eq [18], where P eq is a r X i v : . [ c ond - m a t . s o f t ] A ug the equilibrium distribution of position and orientation.Using Ω eq P eq = 0 yields the equation of motion ∂∂t δP ( t ) = Ω a ( t ) δP ( t ) + δ Ω a ( t ) P eq . (5)Treating the last term as an inhomogeneity and solvingfor δP ( t ) we obtain a formal solution for the nonequilib-rium distribution P ( t ) = P eq − (cid:90) t −∞ dt (cid:48) v ( t (cid:48) ) e (cid:82) tt (cid:48) ds Ω a ( s )+ βF p P eq , (6)where e + ( · ) is a positively ordered exponential function(see the appendix in [8]) and we have used δ Ω a ( t ) P eq = − βv ( t ) F p P eq , with ‘projected force’ fluctuation F p = (cid:88) i p i · F i . (7)The projected force emerges as a central quantity withinour approach and indicates to what extent the interpar-ticle interaction forces act in the direction of orientation,either assisting or hindering the self-propulsion. We willshow that this quantity is closely related to the averageswim speed in the active system.Introducing a test function, f , on the space of positionsand orientations and integrating (6) by parts yields aformally exact expression for a nonequilibrium average (cid:104) f (cid:105) ( t ) = (cid:104) f (cid:105) eq − (cid:90) t −∞ dt (cid:48) v ( t (cid:48) ) (cid:104) βF p e (cid:82) tt (cid:48) ds Ω † a ( s ) − f (cid:105) eq , (8)where e − ( · ) denotes a negatively ordered exponential [8]and (cid:104)·(cid:105) eq is an equilibrium average over positional andorientational degrees of freedom. The adjoint operator isgiven by Ω † a ( t ) = Ω † eq − δ Ω a ( t ), whereΩ † eq = (cid:88) i D t ( ∇ i + β F i ) · ∇ i + D r R i (9)generates the equilibrium dynamics. The integrand ap-pearing in (8) involves the equilibrium correlation be-tween the projected force at time t (cid:48) and the observable f ,which evolves from t (cid:48) to t according to the full dynamics.The average is nonlinear in v ( t ), because of the activitydependence 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 oper-ator Ω † eq . Further simplification occurs if the activity isconstant in time, v ( t ) → v , leading to (cid:104) f (cid:105) lin = (cid:104) f (cid:105) eq − v (cid:90) ∞ dt (cid:104) βF p e Ω † eq t f (cid:105) eq , (10)which can be used to define a general active transportcoefficient α = lim v → ( (cid:104) f (cid:105) lin −(cid:104) f (cid:105) eq ) /v . Equation (10)is the desired Green-Kubo relation for calculating thelinear response of ABPs to a time-independent activity. t/τ p H ( t ) (a) ρ=0.10ρ=0.19ρ=0.31ρ=0.40ρ=0.50ρ=0.57ρ=0.67ρ=0.76ρ=0.95 ρ H ( ) Equilibrium BDPercus Yevick -2 -1 t/τ p -3 -2 -1 H ( t ) / H ( ) (b) ρ=0.10ρ=0.19ρ=0.31ρ=0.40ρ=0.50ρ=0.57ρ=0.67ρ=0.76ρ=0.95 FIG. 1. (a) The correlator H ( t ) for a system of soft spheres.The arrow indicates the direction of increasing density. Inset:The initial value H (0) as a function of the density. Squares:BD simulation. Circles: using equation (16) and the Percus-Yevick g eq ( r ) as input. (b) The same data as in (a) on a logscale. The relaxation becomes faster with increasing density. As mentioned previously, a quantity of current inter-est is the average, density-dependent swim speed, v ( ρ ).This describes how the bare swim speed, v , is influencedby interparticle interactions and is an important quan-tity in many of the various theories addressing ABPs [9–13]. In particular, the tendency of the system to undergomotility-induced phase-separation is determined by therate of decrease of v ( ρ ) with increasing density; a posi-tive feedback mechanism can result when increasing thelocal density leads to a sufficiently strong reduction ofthe local average swim velocity. The average swim speedis defined as the nonequilibrium average v ( ρ ) = 1 N (cid:42)(cid:88) i v i · p i (cid:43) (11)where v i is the velocity of particle i . Using (1) to elim-inate the velocity in favour of the forces and using thefact that the Brownian force ξ i is uncorrelated with theorientation p i , it follows that v ( ρ ) = v + γ − N (cid:104) F p (cid:105) . (12)For a time-independent v we can employ (10) to calcu-late the average in (12) to linear order v ( ρ ) = v (cid:18) − D t (cid:90) ∞ dt H ( t ) (cid:19) , (13)where the integrand is the equilibrium autocorrelation ofprojected force fluctuations H ( t ) = 1 N (cid:104) βF p e Ω † eq t βF p (cid:105) eq . (14)Spatial and orientational degrees of freedom decouple inequilibrium, which enables the orientational integrals in(14) to be evaluated exactly. This yields H ( t ) = 13 e − D r t β (cid:68) F · e Ω † eq , s t F (cid:69) eq , s , (15)where F is the interaction force acting on an arbitrarilychosen (‘tagged’) particle, Ω † eq , s = (cid:80) i D t ( ∇ i + β F i ) · ∇ i is the spatial part of the time-evolution operator and (cid:104)·(cid:105) eq , s indicates an equilibrium average over spatial de-grees of freedom. The initial value is given by H (0) = β (cid:104)| F | (cid:105) eq /
3. If we consider pairwise additive interac-tion potentials, then the Yvon theorem [21] leads to H (0) = 13 ρ (cid:90) d r g eq ( r ) ∇ βu ( r ) , (16)where ρ is the number density, u ( r ) is the passive pair po-tential and g eq ( r ) is the corresponding equilibrium radialdistribution function.Equation (15) shows that the nontrivial physics un-derlying the linear response of the system to activityis contained in the tagged-particle force-autocorrelationfunction. This function was encountered many years agoby Klein and coworkers [19] in a study of the velocityautocorrelation in overdamped Brownian systems. Bymanipulation of the operator (3) it was shown that (cid:104) F ( t ) · F (0) (cid:105) eq , s = 3( βD t ) (cid:0) D t δ ( t ) − Z eq ( t ) (cid:1) , (17)where Z eq ( t ) is the velocity autocorrelation function, de-fined in terms of the tagged particle velocity, v ( t ), ac-cording to the familiar relation Z eq ( t ) = 13 (cid:104) v ( t ) · v (0) (cid:105) eq , s . (18)The velocity autocorrelation function is a quantity of fun-damental interest in describing the dynamics of inter-acting liquids and is closely related to other importantquantities (e.g. the mean-squared displacement and selfdiffusion coefficient). Substituting (17) into (15) yields H ( t ) = 1 D t e − D r t (cid:0) D t δ ( t ) − Z eq ( t ) (cid:1) , (19)thus providing, via (13), a direct connection between Z eq ( t ) and v ( ρ ). The latter can thus be determined tolinear order in v using a standard, equilibrium BD sim-ulation. Finally, we note that H ( t ) remains integrable inall spatial dimensions, because of the exponential in (15).There is thus no principal difficulty in calculating v ( ρ ) intwo dimensions, in contrast to the situation for transportcoefficents, such as the self-diffusion coefficient, for whichthe relevant Green-Kubo time-integral diverges [5].In a recent study of the pressure in active systemsSolon et al. [20] express the density-dependent aver-age swim speed in the form v ( ρ ) = v + I /ρ , where ρ is the bulk number density. The interaction potentialis encoded in the quantity I via its dependence on a static structural correlation between density and polar-ization, which are given, respectively, by the first andsecond harmonic moments of the orientation-resolvedsingle particle density. This leads to the identification ρ v ( ρ ) / v Active BD v =5Active BD v =20Active BD v =80Linear response FIG. 2. The scaled average swim speed as a function of den-sity. Lines with symbols: data from direct calculation of (11)using active BD simulations. Diamonds: the linear responseresult (independent of v ) calculated using the equilibriumtime correlation function data from Fig.1 as input to (13). I = − D t ρ v (cid:82) ∞ dt H ( t ). An advantage of the presentGreen-Kubo formulation over that of Solon et al. is thatit enables identification of the relevant relaxation pro-cesses contributing to the decrease of v ( ρ ). Moreover, weanticipate that (13) will prove more convenient for thedevelopment of approximations.In order to test the range of validity of the linear re-sponse result (13) we perform BD simulations on a three-dimensional system of N = 1000 particles interacting viathe pair-potential βu ( r ) = 4 ε (( σ/r ) − ( σ/r ) ), where σ sets the length scale and we set ε = 1. The poten-tial is truncated at its minimum, r = 2 / σ to yield asoftly repulsive interaction. The system size L is de-termined as L = ( N/ρ ) / in order to obtain the desireddensity. The integration time step is fixed to dt = 10 − τ B where τ B = d /D t is the time-scale of translational diffu-sion. Measurements are made after a minimum time of20 τ B to ensure equilibration. In order to measure time-correlations the system is sampled every τ p /
100 s, where τ p = 1 / D R is the rotational diffusion time scale. Thetotal run time is 300 τ B . We choose the ratio of diffu-sion coefficients as D r /D t = 20, although there is nothingspecial about this particular choice.In Fig. 1a we show the correlator H ( t ) as a functionof time for a number of different densities, the largestof which is close to the freezing transition for our modelinteraction potential. Aside from the strong increase of H (0) with increasing density (shown in the inset), themost striking aspect of the correlator is that the de-cay of H ( t ) is much faster than the timescale of rota-tional diffusion (note that time is scaled with τ p in thefigure). Indeed, very large values of the ratio D r /D t would be required for the exponential factor in (15) tosignificantly influence the decay of H ( t ). In the limit oflarge D r we obtain H ( t ) = H (0) exp( − D r t ) and thus
10 20 30 40 50 60 70 80 v ρ ∆vv(ρ) ≤5% FIG. 3. The region for which linear response (13) agreeswith the result of active-BD simulations to a relative errorless than 5%. ∆ v is the difference between linear responseand the active BD simulation result. The breakdown of lin-ear response is related to the onset of activity-induced phaseseparation. v ( ρ ) /v = 1 − H (0) D t / (2 D r ). We conclude that, pro-vided the value of D r is not extremely large, the relevantrelaxation process is the decorrelation of the tagged par-ticle interaction force.In the inset to Fig. 1a we show the initial value, H (0),as a function of the density. To check the expression (16)we have confirmed that using g eq ( r ) from our equilibriumsimulations to evaluate the r.h.s. indeed reproduces the t → H ( t ) data. Moreover, wehave also employed an approximate liquid-state integralequation theory (Percus-Yevick theory) [21] to calculate g eq ( r ) and evaluate H (0). Very good agreement of thepredicted H (0) with simulation data is obtained.In Fig. 1b we replot the data on a semi-logarithmicscale, with the initial value scaled out. This representa-tion makes clear that H ( t ) is non-exponential and thatthe 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 ob-servation can be rationalized by considering that smallpositional changes can give rise to large changes in theforce for closely packed particles residing in regions ofstrong interaction-force gradient. The fact that H ( t ) isnon-exponential is not surprising, given that it can be ex-pressed in terms of the velocity autocorrelation function,a quantity which famously exhibits power law asymptoticbehaviour (‘long-time tails’) [19, 21]. Klein et al. haveshown analytically that for a dilute system of Brownianhard-spheres (cid:104) F ( t ) · F (0) (cid:105) eq , s ∼ t − for long times.In Fig. 2 we show simulation data for the average swimspeed as a function of density. The red diamonds showthe linear response prediction obtained by using the dataof Fig. 1 in the integral expression (13). This yields aresult for v ( ρ ) /v which is independent of v . The re-maining curves show data obtained by direct evaluation of (11) using active BD simulations at three different val-ues of v . As one might expect, deviations from linearresponse occur at lower density for larger values of v .The above observation can be made more concrete byestimating a region in the ( v , φ ) plane where linear re-sponse breaks down. In Fig. 3 we use our simulationdata to map the locus of points for which the error inthe linear response result, relative to the full active BDsimulations, equals 5%. Although the chosen criterion issomewhat arbitrary, it at least gives a visual impressionof the range of validity of linear response within the spaceof our control parameters. The locus of points shown inFig. 3 is correlated with the onset of strong spatial inho-mogeneities and phase separation. However, an analysisof active phase separation would go beyond the scope ofthe present work. The linear response formula (13) thusappears to be reliable for parameter values away fromphase separation, but, beyond this, higher orders in v will become important in determining v ( ρ ).To summarize our main findings: we have derived aformally exact expression (8) for calculating averages ina system of interacting Brownian particles, subject to atime-dependent activity v ( t ). From this we obtain thelinear-response expression (10) for a time-independentactivity. Application of this result to calculate the av-erage swim speed yields (13) and identifies the relevanttime-correlation function, H ( t ), as given by (15). We findthat linear response provides an accurate account of v ( ρ )over a large parameter range, except for those regions ofparameter space where phase separation occurs.Although we have focused our attention on the linear-response regime, our exact results could in principle beused to develop nonlinear theories in the spirit of Refs. [6–8], which address Brownian particles under external flow.It would also be interesting to use (8) to investigate thetransient dynamics arising from time-dependent activity,but we defer this line of enquiry until an experimentallyrelevant protocol can be identified. Aside from using anequilibrium integral equation theory to determine H (0)(inset to Fig. 1a), all of the data presented comes fromBD simulation. A clear next step is to investigate approx-imations to H ( t ) which enable predictions to be madefrom first-principles, without simulation input. Given therelation (19) it seems likely that existing approximationsto the velocity autocorrelation function (e.g. projectionoperator approaches) could be usefully exploited. [1] M. S. Green, J.Chem.Phys. 398 (1954).[2] R. Kubo, J.Phys.Soc.Jpn. 570 (1957).[3] R. Kubo, M. Toda and N. Hashitsume, Statisti-cal Physics II: Nonequilibrium Statistical Mechanics (Springer, Berlin, 1991).[4] U. Seifert, Phys. Rev. Lett.
Statistical mechanics of nonequilibrium liquids (Cambridge University Press,2008).[6] M. Fuchs, M.E. Cates, Phys.Rev.Lett. S1681(2005).[8] J. M. Brader, M. E. Cates, and M. Fuchs, Phys.Rev.E
219 (2015).[10] M. E. Cates, D. Marenduzzo, I. Pagonabarraga and J.Tailleur, Proc. Natl. Acad. Sci. U.S.A.
439 (2013).[15] F. Schweitzer, W. Ebeling and B. Tilch, Phys. Rev. Lett. Handbook of stochastic methods (Springer,Berlin, 1985).[17] P. M. Morse and H. Feshbach
Methods of TheoreticalPhysics , page 33 (McGraw-Hill, New York, 1953).[18] M. Fuchs and M. E. Cates, J. Phys.: Condens. Matter, S1681 (2005).[19] S. Hanna, W. Hess and R. Klein, J. Phys. A: Math. Gen. L493 (1981).[20] A. P. Solon et al. , Phys. Rev. Lett.
Theory of simple liq-uids , 3 rdrd